continuous_time_ising.hpp

Go to the documentation of this file.

//    Copyright 2023 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 <cassert>
#include <utility>
#include <vector>
 
#include "openjij/graph/all.hpp"
#include "openjij/system/system.hpp"
#include "openjij/utility/eigen.hpp"
 
namespace openjij {
namespace system {
 
template <typename GraphType> struct ContinuousTimeIsing;
 
template <typename FloatType>
struct ContinuousTimeIsing<graph::Sparse<FloatType>> {
  using GraphType = graph::Sparse<FloatType>;
  using system_type = transverse_field_system;
  using TimeType = FloatType;
  using CutPoint = std::pair<TimeType, graph::Spin>;
 
  using SparseMatrixXx = Eigen::SparseMatrix<FloatType, Eigen::RowMajor>;
 
  using SpinConfiguration = std::vector<std::vector<CutPoint>>;
 
  ContinuousTimeIsing(const SpinConfiguration &init_spin_config,
                      const GraphType &init_interaction, const double gamma)
      : spin_config(init_spin_config), num_spins(init_spin_config.size() + 1),
        interaction(
            utility::gen_matrix_from_graph<Eigen::RowMajor>(init_interaction)),
        gamma(gamma) {
 
    assert(init_spin_config.size() == init_interaction.get_num_spins());
 
    // insert numbers to diagnonal elements in the interaction
    SparseMatrixXx diag(init_interaction.get_num_spins() + 1,
                        init_interaction.get_num_spins() + 1);
 
    for (typename SparseMatrixXx::InnerIterator it(
             interaction, init_interaction.get_num_spins());
         it; ++it) {
      std::size_t j = it.index();
      FloatType v = it.value();
      if (j != init_interaction.get_num_spins()) {
        diag.insert(j, j) = v;
      } else {
        diag.insert(j, j) = -1;
      }
    }
 
    interaction += diag;
 
    spin_config.push_back(std::vector<CutPoint>{CutPoint(0.0, 1)});
    // initialize auxiliary spin with 1 along entire timeline
  }
 
  ContinuousTimeIsing(const graph::Spins &init_spins,
                      const GraphType &init_interaction, const double gamma)
      : ContinuousTimeIsing(convert_to_spin_config(init_spins),
                            init_interaction, gamma) {
  } // constructor delegation
 
  /* initialization helper functions */
private:
  static SpinConfiguration convert_to_spin_config(const graph::Spins &spins) {
    SpinConfiguration spin_config;
    spin_config.reserve(spins.size());
 
    /* set spin configuration for each site, which has only one cut at time zero
     */
    for (auto spin : spins) {
      spin_config.push_back(std::vector<CutPoint>{
          CutPoint(TimeType(), spin)}); // TimeType() is zero value of the type
    }
 
    return spin_config;
  }
 
  /* member functions*/
public:
  void reset_spins(const SpinConfiguration &init_spin_config) {
    assert(init_spin_config.size() == this->num_spins - 1);
 
    this->spin_config = init_spin_config;
    this->spin_config.push_back(std::vector<CutPoint>{CutPoint(TimeType(), 1)});
    // add auxiliary timeline
  }
 
  void reset_spins(const graph::Spins &classical_spins) {
    assert(classical_spins.size() == this->num_spins - 1);
 
    for (size_t i = 0; i < this->num_spins - 1; i++) {
      this->spin_config[i] = std::vector<CutPoint>{
          CutPoint(TimeType(),
                   classical_spins[i]) // TimeType() is zero value of the type
      };
    }
    this->spin_config[this->num_spins - 1] =
        std::vector<CutPoint>{CutPoint(TimeType(), 1)};
  }
 
  size_t get_temporal_spin_index(graph::Index site_index,
                                 TimeType time_point) const {
    static const auto first_lt = [](CutPoint x, CutPoint y) {
      return x.first < y.first;
    };
    // function to compare two time points (lt; less than)
 
    const auto &timeline = this->spin_config[site_index];
    const auto dummy_cut =
        CutPoint(time_point, 0); // dummy variable for binary search
    auto found_itr =
        std::upper_bound(timeline.begin(), timeline.end(), dummy_cut, first_lt);
 
    if (found_itr ==
        timeline.begin()) { // if the time_point lies before any time points
      found_itr = timeline.end() - 1; // periodic boundary condition
    } else {
      found_itr--;
    }
 
    return std::distance(timeline.begin(), found_itr);
  }
 
  /*
   * @brief return spin configuration at given temporal slice, not containing
   * auxiliary spin
   *
   * @param slice_time
   */
  graph::Spins get_slice_at(TimeType slice_time) const {
    graph::Spins slice;
 
    for (graph::Index i = 0; i < this->spin_config.size() - 1; i++) {
      auto temporal_index = get_temporal_spin_index(i, slice_time);
      slice.push_back(this->spin_config[i][temporal_index].second);
    }
 
    return slice;
  }
 
  /*
   * @brief return auxiliary spin state at given temporal slice
   *
   * @param slice_time
   */
  graph::Spin get_auxiliary_spin(TimeType slice_time) const {
    auto last_index = this->spin_config.size() - 1;
    auto temporal_index = get_temporal_spin_index(last_index, slice_time);
    return this->spin_config[last_index][temporal_index].second;
  }
 
  /* Member variables */
 
  SpinConfiguration spin_config;
 
  const std::size_t num_spins;
 
  SparseMatrixXx interaction;
 
  const FloatType gamma;
};
 
template <typename FloatType>
struct ContinuousTimeIsing<graph::CSRSparse<FloatType>> {
  using GraphType = graph::CSRSparse<FloatType>;
  using system_type = transverse_field_system;
  using TimeType = FloatType;
  using CutPoint = std::pair<TimeType, graph::Spin>;
 
  using SparseMatrixXx = Eigen::SparseMatrix<FloatType, Eigen::RowMajor>;
 
  using SpinConfiguration = std::vector<std::vector<CutPoint>>;
 
  ContinuousTimeIsing(const SpinConfiguration &init_spin_config,
                      const GraphType &init_interaction, const double gamma)
      : spin_config(init_spin_config), num_spins(init_spin_config.size() + 1),
        interaction(init_interaction.get_interactions()),
        gamma(gamma) {
 
    assert(init_spin_config.size() == init_interaction.get_num_spins());
 
    // insert numbers to diagnonal elements in the interaction
    SparseMatrixXx diag(init_interaction.get_num_spins() + 1,
                        init_interaction.get_num_spins() + 1);
 
    for (typename SparseMatrixXx::InnerIterator it(
             interaction, init_interaction.get_num_spins());
         it; ++it) {
      std::size_t j = it.index();
      FloatType v = it.value();
      if (j != init_interaction.get_num_spins()) {
        diag.insert(j, j) = v;
      } else {
        diag.insert(j, j) = -1;
      }
    }
 
    interaction += diag;
 
    spin_config.push_back(std::vector<CutPoint>{CutPoint(0.0, 1)});
    // initialize auxiliary spin with 1 along entire timeline
  }
 
  ContinuousTimeIsing(const graph::Spins &init_spins,
                      const GraphType &init_interaction, const double gamma)
      : ContinuousTimeIsing(convert_to_spin_config(init_spins),
                            init_interaction, gamma) {
  } // constructor delegation
 
  /* initialization helper functions */
private:
  static SpinConfiguration convert_to_spin_config(const graph::Spins &spins) {
    SpinConfiguration spin_config;
    spin_config.reserve(spins.size());
 
    /* set spin configuration for each site, which has only one cut at time zero
     */
    for (auto spin : spins) {
      spin_config.push_back(std::vector<CutPoint>{
          CutPoint(TimeType(), spin)}); // TimeType() is zero value of the type
    }
 
    return spin_config;
  }
 
  /* member functions*/
public:
  void reset_spins(const SpinConfiguration &init_spin_config) {
    assert(init_spin_config.size() == this->num_spins - 1);
 
    this->spin_config = init_spin_config;
    this->spin_config.push_back(std::vector<CutPoint>{CutPoint(TimeType(), 1)});
    // add auxiliary timeline
  }
 
  void reset_spins(const graph::Spins &classical_spins) {
    assert(classical_spins.size() == this->num_spins - 1);
 
    for (size_t i = 0; i < this->num_spins - 1; i++) {
      this->spin_config[i] = std::vector<CutPoint>{
          CutPoint(TimeType(),
                   classical_spins[i]) // TimeType() is zero value of the type
      };
    }
    this->spin_config[this->num_spins - 1] =
        std::vector<CutPoint>{CutPoint(TimeType(), 1)};
  }
 
  size_t get_temporal_spin_index(graph::Index site_index,
                                 TimeType time_point) const {
    static const auto first_lt = [](CutPoint x, CutPoint y) {
      return x.first < y.first;
    };
    // function to compare two time points (lt; less than)
 
    const auto &timeline = this->spin_config[site_index];
    const auto dummy_cut =
        CutPoint(time_point, 0); // dummy variable for binary search
    auto found_itr =
        std::upper_bound(timeline.begin(), timeline.end(), dummy_cut, first_lt);
 
    if (found_itr ==
        timeline.begin()) { // if the time_point lies before any time points
      found_itr = timeline.end() - 1; // periodic boundary condition
    } else {
      found_itr--;
    }
 
    return std::distance(timeline.begin(), found_itr);
  }
 
  /*
   * @brief return spin configuration at given temporal slice, not containing
   * auxiliary spin
   *
   * @param slice_time
   */
  graph::Spins get_slice_at(TimeType slice_time) const {
    graph::Spins slice;
 
    for (graph::Index i = 0; i < this->spin_config.size() - 1; i++) {
      auto temporal_index = get_temporal_spin_index(i, slice_time);
      slice.push_back(this->spin_config[i][temporal_index].second);
    }
 
    return slice;
  }
 
  /*
   * @brief return auxiliary spin state at given temporal slice
   *
   * @param slice_time
   */
  graph::Spin get_auxiliary_spin(TimeType slice_time) const {
    auto last_index = this->spin_config.size() - 1;
    auto temporal_index = get_temporal_spin_index(last_index, slice_time);
    return this->spin_config[last_index][temporal_index].second;
  }
 
  /* Member variables */
 
  SpinConfiguration spin_config;
 
  const std::size_t num_spins;
 
  SparseMatrixXx interaction;
 
  const FloatType gamma;
};
 
template <typename GraphType>
ContinuousTimeIsing<GraphType>
make_continuous_time_ising(const graph::Spins &init_spins,
                           const GraphType &init_interaction, double gamma) {
  return ContinuousTimeIsing<GraphType>(init_spins, init_interaction, gamma);
}
} // namespace system
} // namespace openjij