sparse.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 <algorithm>
#include <cassert>
#include <cstddef>
#include <type_traits>
#include <unordered_map>
#include <utility>
 
#include "openjij/graph/graph.hpp"
#include "openjij/graph/json/parse.hpp"
#include "openjij/utility/pairhash.hpp"
 
namespace openjij {
namespace graph {
 
template <typename FloatType> class Sparse : public Graph {
  static_assert(std::is_floating_point<FloatType>::value,
                "FloatType must be floating-point type.");
 
public:
  using Interactions =
      std::unordered_map<std::pair<Index, Index>, FloatType, utility::PairHash>;
 
  using value_type = FloatType;
 
private:
  Interactions _J;
 
  std::size_t _num_edges;
 
  std::vector<Nodes> _list_adj_nodes;
 
  inline bool set_adj_node(Index from, Index to) {
    assert(from < this->get_num_spins());
    assert(to < this->get_num_spins());
 
    // get adjacent nodes of node "from"
    Nodes &nodes = _list_adj_nodes[from];
    // check if the node "to" exists in the nodes
    if (std::find(nodes.begin(), nodes.end(), to) == nodes.end()) {
      // nodes size must be smaller than num_edges
 
#ifndef NDEBUG
      assert(nodes.size() < _num_edges);
#else
      // return false
      if (nodes.size() >= _num_edges)
        return false;
#endif
      // add node
      nodes.push_back(to);
      // add node from "to" to "from"
      set_adj_node(to, from);
    }
    return true;
  }
 
public:
  Sparse(std::size_t num_spins, std::size_t num_edges)
      : Graph(num_spins), _num_edges(std::min(num_spins, num_edges)),
        _list_adj_nodes(num_spins) {}
 
  explicit Sparse(std::size_t num_spins) : Sparse(num_spins, num_spins) {}
 
  Sparse(const json &j, std::size_t num_edges)
      : Sparse(static_cast<std::size_t>(j["num_variables"]), num_edges) {
 
    // define SparseMatrix and iterator
    using SparseMatrix = Eigen::SparseMatrix<FloatType, Eigen::RowMajor>;
    using SpIter = typename SparseMatrix::InnerIterator;
 
    // define bqm with ising variables
    auto bqm = json_parse<FloatType, cimod::Sparse>(j);
    // interactions
    // for(auto&& elem : bqm.get_quadratic()){
    //    const auto& key = elem.first;
    //    const auto& val = elem.second;
    //    J(key.first, key.second) += val;
    //}
    // local field
    // for(auto&& elem : bqm.get_linear()){
    //    const auto& key = elem.first;
    //    const auto& val = elem.second;
    //    h(key) += val;
    //}
 
    // insert elements
    SparseMatrix quadmat = bqm.interaction_matrix();
    size_t num_variables = quadmat.rows() - 1;
    for (int k = 0; k < quadmat.outerSize(); k++) {
      for (SpIter it(quadmat, k); it; ++it) {
        size_t r = it.row();
        size_t c = it.col();
        FloatType val = it.value();
 
        if (r == num_variables && c == num_variables)
          continue;
 
        if (c == num_variables) {
          // local field
          h(r) += val;
        } else {
          // quadratic
          J(r, c) += val;
        }
      }
    }
  }
 
  Sparse(const json &j) : Sparse(j, j["num_variables"]) {}
 
  Sparse(const Sparse<FloatType> &) = default;
 
  Sparse(Sparse<FloatType> &&) = default;
 
  const Nodes &adj_nodes(Index ind) const { return _list_adj_nodes[ind]; }
 
  std::size_t get_num_edges() const { return _num_edges; }
 
  FloatType calc_energy(const Spins &spins) const {
    return this->energy(spins);
  }
 
  FloatType calc_energy(const Eigen::Matrix<FloatType, Eigen::Dynamic, 1,
                                            Eigen::ColMajor> &spins) const {
    return this->energy(spins);
  }
 
  FloatType energy(const Spins &spins) const {
    if (!(spins.size() == this->get_num_spins())) {
      std::out_of_range("Out of range in energy in Sparse graph.");
    }
 
    FloatType ret = 0;
    for (std::size_t ind = 0; ind < this->get_num_spins(); ind++) {
      for (auto &adj_ind : _list_adj_nodes[ind]) {
        if (ind != adj_ind)
          ret += (1. / 2) * this->J(ind, adj_ind) * spins[ind] * spins[adj_ind];
        else
          ret += this->h(ind) * spins[ind];
      }
    }
 
    return ret;
  }
 
  FloatType energy(const Eigen::Matrix<FloatType, Eigen::Dynamic, 1,
                                       Eigen::ColMajor> &spins) const {
    graph::Spins temp_spins(get_num_spins());
    for (size_t i = 0; i < temp_spins.size(); i++) {
      temp_spins[i] = spins(i);
    }
    return energy(temp_spins);
  }
 
  FloatType &J(Index i, Index j) {
    assert(i < get_num_spins());
    assert(j < get_num_spins());
 
    // i <= j
    // add node if it does not exist
    set_adj_node(i, j);
    return _J[std::make_pair(std::min(i, j), std::max(i, j))];
  }
 
  const FloatType &J(Index i, Index j) const {
    assert(i < get_num_spins());
    assert(j < get_num_spins());
    return _J.at(std::make_pair(std::min(i, j), std::max(i, j)));
  }
 
  FloatType &h(Index i) {
    assert(i < get_num_spins());
    set_adj_node(i, i);
    return _J[std::make_pair(i, i)];
  }
 
  const FloatType &h(Index i) const {
    assert(i < get_num_spins());
    return _J.at(std::make_pair(i, i));
  }
};
} // namespace graph
} // namespace openjij