binary_polynomial_model.hpp

Go to the documentation of this file.

//    Copyright 2022 Jij Inc.
 
//    Licensed under the Apache License, Version 2.0 (the "License");
//    you may not use this file except in compliance with the License.
//    You may obtain a copy of the License at
 
//        http://www.apache.org/licenses/LICENSE-2.0
 
//    Unless required by applicable law or agreed to in writing, software
//    distributed under the License is distributed on an "AS IS" BASIS,
//    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//    See the License for the specific language governing permissions and
//    limitations under the License.
 
#pragma once
 
#include <algorithm>
#include <bitset>
#include <cstdint>
#include <iostream>
#include <set>
#include <string>
#include <tuple>
#include <typeinfo>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
#include <nlohmann/json.hpp>
 
#include "cimod/hash.hpp"
#include "cimod/utilities.hpp"
#include "cimod/vartypes.hpp"
 
 
namespace cimod {
 
  template<typename IndexType, typename FloatType>
  using Polynomial = std::unordered_map<std::vector<IndexType>, FloatType, vector_hash>;
 
  template<typename IndexType>
  using PolynomialKeyList = std::vector<std::vector<IndexType>>;
 
  template<typename FloatType>
  using PolynomialValueList = std::vector<FloatType>;
 
  template<typename IndexType>
  using Sample = std::unordered_map<IndexType, int32_t>;
 
  template<typename IndexType, typename FloatType>
  class BinaryPolynomialModel {
 
  public:
    BinaryPolynomialModel( const Polynomial<IndexType, FloatType> &poly_map, const Vartype vartype ) : vartype_( vartype ) {
      if ( vartype_ == Vartype::NONE ) {
        throw std::runtime_error( "Unknown vartype detected" );
      }
      AddInteractionsFrom( poly_map );
      UpdateVariablesToIntegers();
    }
 
    BinaryPolynomialModel(
        PolynomialKeyList<IndexType> &key_list,
        const PolynomialValueList<FloatType> &value_list,
        const Vartype vartype ) :
        vartype_( vartype ) {
      if ( vartype_ == Vartype::NONE ) {
        throw std::runtime_error( "Unknown vartype detected" );
      }
      AddInteractionsFrom( key_list, value_list );
      UpdateVariablesToIntegers();
    }
 
    BinaryPolynomialModel(
        const PolynomialKeyList<IndexType> &key_list,
        const PolynomialValueList<FloatType> &value_list,
        const Vartype vartype ) :
        vartype_( vartype ) {
      if ( vartype_ == Vartype::NONE ) {
        throw std::runtime_error( "Unknown vartype detected" );
      }
      AddInteractionsFrom( key_list, value_list );
      UpdateVariablesToIntegers();
    }
 
    BinaryPolynomialModel(
        const std::vector<IndexType> &variables,
        const PolynomialKeyList<std::size_t> &poly_key_distance_list,
        const PolynomialValueList<FloatType> &poly_value_list,
        const Vartype vartype ) :
        vartype_( vartype ) {
      if ( vartype_ == Vartype::NONE ) {
        throw std::runtime_error( "Unknown vartype detected" );
      }
 
      if ( poly_key_distance_list.size() != poly_value_list.size() ) {
        throw std::runtime_error( "The sizes of key_list and value_list must match each other" );
      }
 
      variables_ = std::unordered_set<IndexType>( variables.begin(), variables.end() );
 
      if ( variables_.size() != variables.size() ) {
        throw std::runtime_error( "Unknown error. It seems that the input variables contain the same variables" );
      }
 
      std::size_t num_interactions = poly_key_distance_list.size();
      poly_key_list_.resize( num_interactions );
      poly_value_list_.resize( num_interactions );
 
#pragma omp parallel for
      for ( int64_t i = 0; i < ( int64_t )num_interactions; ++i ) {
        std::vector<IndexType> temp;
        for ( const auto &it : poly_key_distance_list[ i ] ) {
          temp.push_back( variables[ it ] );
        }
        std::sort( temp.begin(), temp.end() );
        poly_key_list_[ i ] = temp;
        poly_value_list_[ i ] = poly_value_list[ i ];
      }
 
      for ( std::size_t i = 0; i < num_interactions; ++i ) {
        poly_key_inv_[ poly_key_list_[ i ] ] = i;
        for ( const auto &it : poly_key_list_[ i ] ) {
          each_variable_num_[ it ]++;
        }
      }
 
      UpdateVariablesToIntegers();
    }
 
    Polynomial<IndexType, FloatType> GetPolynomial() const {
      Polynomial<IndexType, FloatType> poly_map;
      for ( std::size_t i = 0; i < poly_key_list_.size(); ++i ) {
        poly_map[ poly_key_list_[ i ] ] = poly_value_list_[ i ];
      }
      return poly_map;
    }
 
    FloatType GetPolynomial( std::vector<IndexType> &key ) const {
      FormatPolynomialKey( &key, vartype_ );
      if ( poly_key_inv_.count( key ) != 0 ) {
        return poly_value_list_[ poly_key_inv_.at( key ) ];
      } else {
        return 0;
      }
    }
 
    FloatType GetPolynomial( const std::vector<IndexType> &key ) const {
      std::vector<IndexType> copied_key = key;
      return GetPolynomial( copied_key );
    }
 
    const std::unordered_map<IndexType, int64_t> &GetVariablesToIntegers() {
      if ( relabel_flag_for_variables_to_integers_ ) {
        UpdateVariablesToIntegers();
      }
      return variables_to_integers_;
    }
 
    std::unordered_map<IndexType, int64_t> GetVariablesToIntegers() const {
      if ( relabel_flag_for_variables_to_integers_ ) {
        return GenerateVariablesToIntegers();
      } else {
        return variables_to_integers_;
      }
    }
 
    int64_t GetVariablesToIntegers( const IndexType &index ) {
      if ( relabel_flag_for_variables_to_integers_ ) {
        UpdateVariablesToIntegers();
      }
      if ( variables_to_integers_.count( index ) == 0 ) {
        return -1;
      } else {
        return variables_to_integers_.at( index );
      }
    }
 
    int64_t GetVariablesToIntegers( const IndexType &index ) const {
      if ( variables_.count( index ) == 0 ) {
        return -1;
      } else {
        std::vector<IndexType> sorted_variables = GetSortedVariables();
        return std::distance(
            sorted_variables.begin(), std::lower_bound( sorted_variables.begin(), sorted_variables.end(), index ) );
      }
    }
 
    const PolynomialKeyList<IndexType> &GetKeyList() const {
      return poly_key_list_;
    }
 
    const PolynomialValueList<FloatType> &GetValueList() const {
      return poly_value_list_;
    }
 
    const std::unordered_map<std::vector<IndexType>, std::size_t, vector_hash> &GetKeysInv() const {
      return poly_key_inv_;
    }
 
    const std::unordered_set<IndexType> &GetVariables() const {
      return variables_;
    }
 
    const std::vector<IndexType> &GetSortedVariables() {
      if ( relabel_flag_for_variables_to_integers_ ) {
        UpdateVariablesToIntegers();
      }
      return sorted_variables_;
    }
 
    std::vector<IndexType> GetSortedVariables() const {
      if ( relabel_flag_for_variables_to_integers_ ) {
        return GenerateSortedVariables();
      } else {
        return sorted_variables_;
      }
    }
 
    std::size_t GetDegree() const {
      std::size_t degree = 0;
      for ( const auto &it : poly_key_list_ ) {
        if ( degree < it.size() ) {
          degree = it.size();
        }
      }
      return degree;
    }
 
    FloatType GetOffset() const {
      return GetPolynomial( std::vector<IndexType>{} );
    }
 
    Vartype GetVartype() const {
      return vartype_;
    }
 
    std::size_t GetNumInteractions() const {
      return poly_key_list_.size();
    }
 
    std::size_t GetNumVariables() const {
      return variables_.size();
    }
 
    BinaryPolynomialModel Empty( const Vartype vartype ) const {
      return BinaryPolynomialModel( {}, vartype );
    }
 
    void Clear() {
      each_variable_num_.clear();
      variables_to_integers_.clear();
      PolynomialKeyList<IndexType>().swap( poly_key_list_ );
      PolynomialValueList<FloatType>().swap( poly_value_list_ );
      std::unordered_set<IndexType>().swap( variables_ );
      poly_key_inv_.clear();
      relabel_flag_for_variables_to_integers_ = true;
    }
 
    void RemoveInteraction( std::vector<IndexType> &key ) {
      FormatPolynomialKey( &key, vartype_ );
      if ( poly_key_inv_.count( key ) == 0 ) {
        return;
      }
 
      for ( const auto &index : key ) {
        if ( each_variable_num_[ index ] >= 2 ) {
          each_variable_num_[ index ]--;
        } else if ( each_variable_num_[ index ] == 1 ) {
          each_variable_num_.erase( index );
          variables_.erase( index );
          relabel_flag_for_variables_to_integers_ = true;
        }
      }
 
      std::size_t inv = poly_key_inv_[ key ];
 
      std::swap( poly_key_inv_[ key ], poly_key_inv_[ poly_key_list_.back() ] );
      poly_key_inv_.erase( key );
 
      std::swap( poly_key_list_[ inv ], poly_key_list_.back() );
      poly_key_list_.pop_back();
 
      std::swap( poly_value_list_[ inv ], poly_value_list_.back() );
      poly_value_list_.pop_back();
    }
 
    void RemoveInteraction( const std::vector<IndexType> &key ) {
      std::vector<IndexType> copied_key = key;
      RemoveInteraction( copied_key );
    }
 
    void RemoveInteractionsFrom( PolynomialKeyList<IndexType> &key_list ) {
      for ( auto &&key : key_list ) {
        RemoveInteraction( key );
      }
    }
 
    void RemoveInteractionsFrom( const PolynomialKeyList<IndexType> &key_list ) {
      for ( const auto &key : key_list ) {
        RemoveInteraction( key );
      }
    }
 
    void RemoveOffset() {
      RemoveInteraction( std::vector<IndexType>{} );
    }
 
    void RemoveVariable( const IndexType &index ) {
      for ( auto &&key : poly_key_list_ ) {
        if ( std::binary_search( key.begin(), key.end(), index ) ) {
          RemoveInteraction( key );
        }
      }
    }
 
    void RemoveVariablesFrom( const std::vector<IndexType> &key ) {
      for ( const auto &index : key ) {
        RemoveVariable( index );
      }
    }
 
    void AddInteraction( std::vector<IndexType> &key, const FloatType &value, const Vartype vartype = Vartype::NONE ) {
      if ( std::abs( value ) <= 0.0 ) {
        return;
      }
 
      if ( vartype_ == vartype || vartype == Vartype::NONE ) {
        FormatPolynomialKey( &key, vartype_ );
        SetKeyAndValue( key, value );
      } else {
        const std::size_t original_key_size = key.size();
        const std::size_t changed_key_list_size = IntegerPower( 2, original_key_size );
 
        if ( vartype_ == Vartype::SPIN && vartype == Vartype::BINARY ) {
          FormatPolynomialKey( &key, vartype );
          for ( std::size_t i = 0; i < changed_key_list_size; ++i ) {
            const auto changed_key = GenerateChangedKey( key, i );
            int sign = ( ( original_key_size - changed_key.size() ) % 2 == 0 ) ? 1.0 : -1.0;
            SetKeyAndValue( changed_key, value * IntegerPower( 2, changed_key.size() ) * sign );
          }
        } else if ( vartype_ == Vartype::BINARY && vartype == Vartype::SPIN ) {
          FormatPolynomialKey( &key, vartype );
          FloatType changed_value = value * ( 1.0 / changed_key_list_size );
          for ( std::size_t i = 0; i < changed_key_list_size; ++i ) {
            SetKeyAndValue( GenerateChangedKey( key, i ), changed_value );
          }
        } else {
          throw std::runtime_error( "Unknown vartype error" );
        }
      }
    }
 
    void AddInteraction( const std::vector<IndexType> &key, const FloatType &value, const Vartype vartype = Vartype::NONE ) {
      std::vector<IndexType> copied_key = key;
      AddInteraction( copied_key, value, vartype );
    }
 
    void AddInteractionsFrom( const Polynomial<IndexType, FloatType> &poly_map, const Vartype vartype = Vartype::NONE ) {
      for ( const auto &it : poly_map ) {
        AddInteraction( it.first, it.second, vartype );
      }
    }
 
    void AddInteractionsFrom(
        PolynomialKeyList<IndexType> &key_list,
        const PolynomialValueList<FloatType> &value_list,
        const Vartype vartype = Vartype::NONE ) {
      if ( key_list.size() != value_list.size() ) {
        throw std::runtime_error( "The sizes of key_list and value_list must match each other" );
      }
      for ( std::size_t i = 0; i < key_list.size(); ++i ) {
        AddInteraction( key_list[ i ], value_list[ i ], vartype );
      }
    }
 
    void AddInteractionsFrom(
        const PolynomialKeyList<IndexType> &key_list,
        const PolynomialValueList<FloatType> &value_list,
        const Vartype vartype = Vartype::NONE ) {
      if ( key_list.size() != value_list.size() ) {
        throw std::runtime_error( "The sizes of key_list and value_list must match each other" );
      }
      for ( std::size_t i = 0; i < key_list.size(); ++i ) {
        std::vector<IndexType> copied_key = key_list[ i ];
        AddInteraction( copied_key, value_list[ i ], vartype );
      }
    }
 
    void AddOffset( FloatType offset ) {
      AddInteraction( std::vector<IndexType>{}, offset );
    }
 
    FloatType Energy( const Sample<IndexType> &sample, bool omp_flag = true ) const {
      if ( sample.size() != GetNumVariables() ) {
        throw std::runtime_error( "The size of sample must be equal to num_variables" );
      }
 
      if ( GetNumInteractions() == 0 ) {
        return 0.0;
      }
 
      std::size_t num_interactions = GetNumInteractions();
      FloatType val = 0.0;
 
      if ( omp_flag ) {
#pragma omp parallel for reduction( + : val )
        for ( int64_t i = 0; i < ( int64_t )num_interactions; ++i ) {
          int32_t spin_multiple = 1;
          for ( const auto &index : poly_key_list_[ i ] ) {
            spin_multiple *= sample.at( index );
            if ( spin_multiple == 0.0 ) {
              break;
            }
          }
          val += spin_multiple * poly_value_list_[ i ];
        }
      } else {
        for ( std::size_t i = 0; i < num_interactions; ++i ) {
          int32_t spin_multiple = 1;
          for ( const auto &index : poly_key_list_[ i ] ) {
            spin_multiple *= sample.at( index );
            if ( spin_multiple == 0.0 ) {
              break;
            }
          }
          val += spin_multiple * poly_value_list_[ i ];
        }
      }
      return val;
    }
 
    FloatType Energy( const std::vector<int32_t> &sample_vec, bool omp_flag = true ) {
      if ( sample_vec.size() != GetNumVariables() ) {
        throw std::runtime_error( "The size of sample must be equal to num_variables" );
      }
 
      if ( GetNumInteractions() == 0 ) {
        return 0.0;
      }
 
      if ( relabel_flag_for_variables_to_integers_ ) {
        UpdateVariablesToIntegers();
      }
 
      std::size_t num_interactions = GetNumInteractions();
      FloatType val = 0.0;
 
      if ( omp_flag ) {
#pragma omp parallel for reduction( + : val )
        for ( int64_t i = 0; i < ( int64_t )num_interactions; ++i ) {
          int32_t spin_multiple = 1;
          for ( const auto &index : poly_key_list_[ i ] ) {
            spin_multiple *= sample_vec[ variables_to_integers_.at( index ) ];
            if ( spin_multiple == 0.0 ) {
              break;
            }
          }
          val += spin_multiple * poly_value_list_[ i ];
        }
      } else {
        for ( std::size_t i = 0; i < num_interactions; ++i ) {
          int32_t spin_multiple = 1;
          for ( const auto &index : poly_key_list_[ i ] ) {
            spin_multiple *= sample_vec[ variables_to_integers_.at( index ) ];
            if ( spin_multiple == 0.0 ) {
              break;
            }
          }
          val += spin_multiple * poly_value_list_[ i ];
        }
      }
      return val;
    }
 
    PolynomialValueList<FloatType> Energies( const std::vector<Sample<IndexType>> &samples ) const {
      PolynomialValueList<FloatType> val_list( samples.size() );
#pragma omp parallel for
      for ( int64_t i = 0; i < ( int64_t )samples.size(); ++i ) {
        val_list[ i ] = Energy( samples[ i ], false );
      }
      return val_list;
    }
 
    PolynomialValueList<FloatType> Energies( const std::vector<std::vector<int32_t>> &samples_vec ) {
      PolynomialValueList<FloatType> val_list( samples_vec.size() );
#pragma omp parallel for
      for ( int64_t i = 0; i < ( int64_t )samples_vec.size(); ++i ) {
        val_list[ i ] = Energy( samples_vec[ i ], false );
      }
      return val_list;
    }
 
    void Scale(
        const FloatType scalar,
        const PolynomialKeyList<IndexType> &ignored_interactions = {},
        const bool ignored_offset = false ) {
 
      std::size_t num_interactions = GetNumInteractions();
 
      for ( std::size_t i = 0; i < num_interactions; ++i ) {
        poly_value_list_[ i ] *= scalar;
      }
 
      FloatType scalar_inv = 1.0 / scalar;
      for ( const auto &key : ignored_interactions ) {
        if ( poly_key_inv_.count( key ) != 0 ) {
          poly_value_list_[ poly_key_inv_[ key ] ] *= scalar_inv;
        }
      }
 
      if ( ignored_offset == true && poly_key_inv_.count( std::vector<IndexType>{} ) != 0
           && std::count( ignored_interactions.begin(), ignored_interactions.end(), std::vector<IndexType>{} ) == 0 ) {
        poly_value_list_[ poly_key_inv_[ std::vector<IndexType>{} ] ] *= scalar_inv;
      }
    }
 
    void normalize(
        const std::pair<FloatType, FloatType> &range = { 1.0, 1.0 },
        const PolynomialKeyList<IndexType> &ignored_interactions = {},
        const bool ignored_offset = false ) {
 
      if ( GetNumInteractions() == 0 ) {
        return;
      }
 
      FloatType max_poly_value = poly_value_list_[ 0 ];
      FloatType min_poly_value = poly_value_list_[ 0 ];
 
      for ( const auto &poly_value : poly_value_list_ ) {
        if ( max_poly_value < poly_value ) {
          max_poly_value = poly_value;
        }
        if ( min_poly_value > poly_value ) {
          min_poly_value = poly_value;
        }
      }
 
      FloatType inv_scale = std::max( min_poly_value / range.first, max_poly_value / range.second );
 
      if ( inv_scale != 0.0 ) {
        Scale( 1.0 / inv_scale, ignored_interactions, ignored_offset );
      }
    }
 
    BinaryPolynomialModel ChangeVartype( const Vartype vartype, const bool inplace ) {
 
      if ( vartype == Vartype::SPIN ) {
        if ( inplace ) {
          *this = ToSpin();
          return *this;
        } else {
          return ToSpin();
        }
      } else if ( vartype == Vartype::BINARY ) {
        if ( inplace ) {
          *this = ToBinary();
          return *this;
        } else {
          return ToBinary();
        }
      } else {
        throw std::runtime_error( "Unknown vartype error" );
      }
    }
 
    void ChangeVartype( const Vartype vartype ) {
      if ( vartype == Vartype::SPIN ) {
        *this = ToSpin();
      } else if ( vartype == Vartype::BINARY ) {
        *this = ToBinary();
      } else {
        throw std::runtime_error( "Unknown vartype error" );
      }
    }
 
    bool HasVariable( const IndexType &index ) {
      if ( variables_.count( index ) != 0 ) {
        return true;
      } else {
        return false;
      }
    }
 
    Polynomial<IndexType, FloatType> ToHubo() const {
      if ( vartype_ == Vartype::BINARY ) {
        return GetPolynomial();
      }
      Polynomial<IndexType, FloatType> poly_map;
      std::size_t num_interactions = GetNumInteractions();
      for ( std::size_t i = 0; i < num_interactions; ++i ) {
        const std::vector<IndexType> &original_key = poly_key_list_[ i ];
        const FloatType original_value = poly_value_list_[ i ];
        const std::size_t original_key_size = original_key.size();
        const std::size_t changed_key_list_size = IntegerPower( 2, original_key_size );
 
        for ( std::size_t j = 0; j < changed_key_list_size; ++j ) {
          const auto changed_key = GenerateChangedKey( original_key, j );
          int sign = ( ( original_key_size - changed_key.size() ) % 2 == 0 ) ? 1.0 : -1.0;
          FloatType changed_value = original_value * IntegerPower( 2, changed_key.size() ) * sign;
          poly_map[ changed_key ] += changed_value;
          if ( poly_map[ changed_key ] == 0.0 ) {
            poly_map.erase( changed_key );
          }
        }
      }
      return poly_map;
    }
 
    Polynomial<IndexType, FloatType> ToHising() const {
      if ( vartype_ == Vartype::SPIN ) {
        return GetPolynomial();
      }
      Polynomial<IndexType, FloatType> poly_map;
      const std::size_t num_interactions = GetNumInteractions();
      for ( std::size_t i = 0; i < num_interactions; ++i ) {
        const std::vector<IndexType> &original_key = poly_key_list_[ i ];
        const FloatType original_value = poly_value_list_[ i ];
        const std::size_t original_key_size = original_key.size();
        const std::size_t changed_key_list_size = IntegerPower( 2, original_key_size );
        const FloatType changed_value = original_value * ( 1.0 / changed_key_list_size );
 
        for ( std::size_t j = 0; j < changed_key_list_size; ++j ) {
          const auto changed_key = GenerateChangedKey( original_key, j );
          poly_map[ changed_key ] += changed_value;
          if ( poly_map[ changed_key ] == 0.0 ) {
            poly_map.erase( changed_key );
          }
        }
      }
      return poly_map;
    }
 
    nlohmann::json ToSerializable() const {
      nlohmann::json output;
      if ( vartype_ == Vartype::BINARY ) {
        output[ "vartype" ] = "BINARY";
      } else if ( vartype_ == Vartype::SPIN ) {
        output[ "vartype" ] = "SPIN";
      } else {
        throw std::runtime_error( "Variable type must be SPIN or BINARY." );
      }
 
      std::size_t num_interactions = GetNumInteractions();
      PolynomialKeyList<std::size_t> poly_key_distance_list( num_interactions );
      std::vector<IndexType> sorted_variables = GetSortedVariables();
 
#pragma omp parallel for
      for ( int64_t i = 0; i < ( int64_t )num_interactions; ++i ) {
        std::vector<std::size_t> temp;
        for ( const auto &it : poly_key_list_[ i ] ) {
          auto it_index = std::lower_bound( sorted_variables.begin(), sorted_variables.end(), it );
          std::size_t index_distance = std::distance( sorted_variables.begin(), it_index );
          temp.push_back( index_distance );
        }
        poly_key_distance_list[ i ] = temp;
      }
 
      output[ "variables" ] = sorted_variables;
      output[ "poly_key_distance_list" ] = poly_key_distance_list;
      output[ "poly_value_list" ] = poly_value_list_;
      output[ "type" ] = "BinaryPolynomialModel";
 
      return output;
    }
 
    template<typename IndexType_serial = IndexType, typename FloatType_serial = FloatType>
    static BinaryPolynomialModel<IndexType_serial, FloatType_serial> FromSerializable( const nlohmann::json &input ) {
      if ( input.at( "type" ) != "BinaryPolynomialModel" ) {
        throw std::runtime_error( "Type must be \"BinaryPolynomialModel\".\n" );
      }
      Vartype vartype;
      if ( input.at( "vartype" ) == "SPIN" ) {
        vartype = Vartype::SPIN;
      } else if ( input.at( "vartype" ) == "BINARY" ) {
        vartype = Vartype::BINARY;
      } else {
        throw std::runtime_error( "Variable type must be SPIN or BINARY." );
      }
      return BinaryPolynomialModel<IndexType_serial, FloatType_serial>(
          input[ "variables" ], input[ "poly_key_distance_list" ], input[ "poly_value_list" ], vartype );
    }
 
    static BinaryPolynomialModel FromHubo( const Polynomial<IndexType, FloatType> &poly_map ) {
      return BinaryPolynomialModel<IndexType, FloatType>( poly_map, Vartype::BINARY );
    }
 
    static BinaryPolynomialModel
    FromHubo( const PolynomialKeyList<IndexType> &key_list, const PolynomialValueList<FloatType> &value_list ) {
      return BinaryPolynomialModel<IndexType, FloatType>( key_list, value_list, Vartype::BINARY );
    }
 
    static BinaryPolynomialModel
    FromHubo( PolynomialKeyList<IndexType> &key_list, const PolynomialValueList<FloatType> &value_list ) {
      return BinaryPolynomialModel<IndexType, FloatType>( key_list, value_list, Vartype::BINARY );
    }
 
    static BinaryPolynomialModel FromHising( const Polynomial<IndexType, FloatType> &poly_map ) {
      return BinaryPolynomialModel<IndexType, FloatType>( poly_map, Vartype::SPIN );
    }
 
    static BinaryPolynomialModel
    FromHising( const PolynomialKeyList<IndexType> &key_list, const PolynomialValueList<FloatType> &value_list ) {
      return BinaryPolynomialModel<IndexType, FloatType>( key_list, value_list, Vartype::SPIN );
    }
 
    static BinaryPolynomialModel
    FromHising( PolynomialKeyList<IndexType> &key_list, const PolynomialValueList<FloatType> &value_list ) {
      return BinaryPolynomialModel<IndexType, FloatType>( key_list, value_list, Vartype::SPIN );
    }
 
  protected:
    std::unordered_set<IndexType> variables_;
 
    std::unordered_map<IndexType, std::size_t> each_variable_num_;
 
    std::unordered_map<IndexType, int64_t> variables_to_integers_;
 
    std::vector<IndexType> sorted_variables_;
 
    bool relabel_flag_for_variables_to_integers_ = true;
 
    PolynomialKeyList<IndexType> poly_key_list_;
 
    PolynomialValueList<FloatType> poly_value_list_;
 
    std::unordered_map<std::vector<IndexType>, std::size_t, vector_hash> poly_key_inv_;
 
    Vartype vartype_ = Vartype::NONE;
 
    void SetKeyAndValue( const std::vector<IndexType> &key, const FloatType &value ) {
      // key is assumed to be sorted
      if ( poly_key_inv_.count( key ) == 0 ) {
        poly_key_inv_[ key ] = poly_value_list_.size();
        poly_key_list_.push_back( key );
        poly_value_list_.push_back( value );
      } else {
        if ( poly_value_list_[ poly_key_inv_[ key ] ] + value == 0.0 ) {
          RemoveInteraction( key );
          return;
        }
        poly_value_list_[ poly_key_inv_[ key ] ] += value;
      }
      for ( const auto &index : key ) {
        each_variable_num_[ index ]++;
        variables_.emplace( index );
        relabel_flag_for_variables_to_integers_ = true;
      }
    }
 
    std::size_t IntegerPower( std::size_t base, std::size_t exponent ) const {
      std::size_t val = 1;
      for ( std::size_t i = 0; i < exponent; ++i ) {
        val *= base;
      }
      return val;
    }
 
    std::vector<IndexType>
    GenerateChangedKey( const std::vector<IndexType> &original_key, const std::size_t num_of_key ) const {
      if ( original_key.size() >= UINT16_MAX ) {
        throw std::runtime_error( "Too large degree of the interaction" );
      }
      const std::size_t original_key_size = original_key.size();
      std::bitset<UINT16_MAX> bs( num_of_key );
      std::vector<IndexType> changed_key;
      for ( std::size_t i = 0; i < original_key_size; ++i ) {
        if ( bs[ i ] ) {
          changed_key.push_back( original_key[ i ] );
        }
      }
      return changed_key;
    }
 
    BinaryPolynomialModel ToSpin() const {
      if ( vartype_ == Vartype::SPIN ) {
        return *this;
      }
      return BinaryPolynomialModel( ToHising(), Vartype::SPIN );
    }
 
    BinaryPolynomialModel ToBinary() const {
      if ( vartype_ == Vartype::BINARY ) {
        return *this;
      }
      return BinaryPolynomialModel( ToHubo(), Vartype::BINARY );
    }
 
    void UpdateVariablesToIntegers() {
      sorted_variables_ = GenerateSortedVariables();
      variables_to_integers_.clear();
      for ( std::size_t i = 0; i < sorted_variables_.size(); ++i ) {
        variables_to_integers_[ sorted_variables_[ i ] ] = i;
      }
      relabel_flag_for_variables_to_integers_ = false;
    }
 
    std::unordered_map<IndexType, int64_t> GenerateVariablesToIntegers() const {
      std::vector<IndexType> sorted_variables = GenerateSortedVariables();
      std::unordered_map<IndexType, int64_t> variables_to_integers;
      for ( std::size_t i = 0; i < sorted_variables.size(); ++i ) {
        variables_to_integers[ sorted_variables[ i ] ] = i;
      }
      return variables_to_integers;
    }
 
    std::vector<IndexType> GenerateSortedVariables() const {
      std::vector<IndexType> sorted_variables( variables_.begin(), variables_.end() );
      std::sort( sorted_variables.begin(), sorted_variables.end() );
      return sorted_variables;
    }
  };
 
} // namespace cimod