OPT_METROPOLISを４次項まで対応 (#502)

# 変更

`sample_huio`関数において、`updater=OPT_METROPOLIS`とすると、ある変数に対して最もエネルギーが下がるような最適遷移が行われます。
現在のmainの実装ではこの最適遷移が行われるのはz^2のような２乗項だけでしたが、これを４乗項まで行うように変更しました。

---------

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

Author: K-Suzuki-Jij

include/openjij/graph/integer_polynomial_model.hpp

@@ -83,29 +83,19 @@ class IntegerPolynomialModel {
       std::sort(interactions.begin(), interactions.end());
     }
 
-    // Create only_multilinear_index_set and under_quadratic_index_set
+    // Create max degree of each variable
+    this->max_variable_degree_.resize(this->num_variables_, 0);
     for (std::int64_t index = 0; index < this->num_variables_; ++index) {
-      bool is_multilinear = true;
+      std::int64_t max_degree = 0;
       for (const auto &[_, degree] : this->index_to_interactions_[index]) {
-        if (degree != 1) {
-          is_multilinear = false;
-          break;
+        if (degree > max_degree) {
+          max_degree = degree;
         }
       }
-      if (is_multilinear) {
-        this->only_multilinear_index_set_.insert(index);
-      }
-
-      bool is_under_quadratic = true;
-      for (const auto &[_, degree] : this->index_to_interactions_[index]) {
-        if (degree > 2) {
-          is_under_quadratic = false;
-          break;
-        }
-      }
-      if (is_under_quadratic) {
-        this->under_quadratic_index_set_.insert(index);
+      if (max_degree == 0) {
+        throw std::runtime_error("Variable with no interactions found.");
       }
+      this->max_variable_degree_[index] = max_degree;
     }
   }
 
@@ -182,11 +172,11 @@ class IntegerPolynomialModel {
   GetIndexToInteractions() const {
     return this->index_to_interactions_;
   }
-  const std::unordered_set<std::int64_t> &GetOnlyMultilinearIndexSet() const {
-    return this->only_multilinear_index_set_;
+  std::int64_t GetEachVariableDegreeAt(const std::int64_t index) const {
+    return this->max_variable_degree_[index];
   }
-  const std::unordered_set<std::int64_t> &GetUnderQuadraticIndexSet() const {
-    return this->under_quadratic_index_set_;
+  const std::vector<std::int64_t> &GetEachVariableDegree() const {
+    return this->max_variable_degree_;
   }
 
 private:
@@ -199,8 +189,7 @@ class IntegerPolynomialModel {
       key_value_list_;
   std::vector<std::vector<std::pair<std::size_t, std::int64_t>>>
       index_to_interactions_;
-  std::unordered_set<std::int64_t> only_multilinear_index_set_;
-  std::unordered_set<std::int64_t> under_quadratic_index_set_;
+  std::vector<std::int64_t> max_variable_degree_;
 };
 
 } // namespace graph

include/openjij/system/integer_polynomial_sa_system.hpp

@@ -15,6 +15,7 @@
 #pragma once
 
 #include "../utility/variable.hpp"
+#include "../utility/min_polynomial.hpp"
 #include "./sa_system.hpp"
 #include <cstdint>
 #include <random>
@@ -186,64 +187,69 @@ class IntegerSASystem<graph::IntegerPolynomialModel, RandType> {
   }
 
   std::pair<int, double> GetMinEnergyDifference(std::int64_t index) {
-    if (this->UnderQuadraticCoeff(index)) {
-      const auto &x = this->state_[index];
-      const std::int64_t dxl = x.lower_bound - x.value;
-      const std::int64_t dxu = x.upper_bound - x.value;
+    const auto &x = this->state_[index];
+    const std::int64_t dxl = x.lower_bound - x.value;
+    const std::int64_t dxu = x.upper_bound - x.value;
 
+    if (this->UnderQuadraticCoeff(index)) {
       auto it_a = this->base_energy_difference_[index].find(2);
       const double a = (it_a != this->base_energy_difference_[index].end())
                            ? it_a->second
                            : 0.0;
-
       auto it_b = this->base_energy_difference_[index].find(1);
-      const double b_base = (it_b != this->base_energy_difference_[index].end())
+      const double b = (it_b != this->base_energy_difference_[index].end())
                                 ? it_b->second
                                 : 0.0;
-      const double b = b_base + 2 * x.value * a;
-
-      if (a > 0) {
-        const double center = -b / (2 * a);
-        if (dxu <= center) {
-          return {x.upper_bound,
-                  this->GetEnergyDifference(index, x.upper_bound)};
-        } else if (dxl < center && center < dxu) {
-          const std::int64_t dx_left =
-              static_cast<std::int64_t>(std::floor(center));
-          const std::int64_t dx_right =
-              static_cast<std::int64_t>(std::ceil(center));
-          if (center - dx_left <= dx_right - center) {
-            return {x.value + dx_left,
-                    this->GetEnergyDifference(index, x.value + dx_left)};
-          } else {
-            return {x.value + dx_right,
-                    this->GetEnergyDifference(index, x.value + dx_right)};
-          }
-        } else if (dxl >= center) {
-          return {x.lower_bound,
-                  this->GetEnergyDifference(index, x.lower_bound)};
-        } else {
-          throw std::runtime_error("Invalid state in GetMinEnergyDifference");
-        }
-      } else if (a == 0) {
-        if (b > 0) {
-          return {x.lower_bound,
-                  this->GetEnergyDifference(index, x.lower_bound)};
-        } else if (b < 0) {
-          return {x.upper_bound,
-                  this->GetEnergyDifference(index, x.upper_bound)};
-        } else {
-          return {x.GenerateRandomValue(this->random_number_engine), 0.0};
-        }
-      } else { // a < 0
-        const double dE_lower = this->GetEnergyDifference(index, x.lower_bound);
-        const double dE_upper = this->GetEnergyDifference(index, x.upper_bound);
-        if (dE_lower <= dE_upper) {
-          return {x.lower_bound, dE_lower};
-        } else {
-          return {x.upper_bound, dE_upper};
-        }
-      }
+
+      const double aa = a;
+      const double bb = b + 2 * x.value * a;
+
+      return utility::FindMinimumIntegerQuadratic(aa, bb, dxl, dxu, x.value, this->random_number_engine);
+    }
+    else if (this->IsCubicCoeff(index)) {
+      auto it_a = this->base_energy_difference_[index].find(3);
+      const double a = (it_a != this->base_energy_difference_[index].end())
+                           ? it_a->second
+                           : 0.0;
+      auto it_b = this->base_energy_difference_[index].find(2);
+      const double b = (it_b != this->base_energy_difference_[index].end())
+                                ? it_b->second
+                                : 0.0;
+      auto it_c = this->base_energy_difference_[index].find(1);
+      const double c = (it_c != this->base_energy_difference_[index].end())
+                                ? it_c->second
+                                : 0.0;
+
+      const double aa = a;
+      const double bb = 3 * a * x.value + b;
+      const double cc = 3 * a * x.value * x.value + 2 * b * x.value + c;
+
+      return utility::FindMinimumIntegerCubic(aa, bb, cc, dxl, dxu, x.value, this->random_number_engine);
+    } else if (this->IsQuarticCoeff(index)) {
+      auto it_a = this->base_energy_difference_[index].find(4);
+      const double a = (it_a != this->base_energy_difference_[index].end())
+                           ? it_a->second
+                           : 0.0;
+      auto it_b = this->base_energy_difference_[index].find(3);
+      const double b = (it_b != this->base_energy_difference_[index].end())
+                                ? it_b->second
+                                : 0.0;
+      auto it_c = this->base_energy_difference_[index].find(2);
+      const double c = (it_c != this->base_energy_difference_[index].end())
+                                ? it_c->second
+                                : 0.0;
+      auto it_d = this->base_energy_difference_[index].find(1);
+      const double d = (it_d != this->base_energy_difference_[index].end())
+                                ? it_d->second
+                                : 0.0;
+
+      const double aa = a;
+      const double bb = 4 * a * x.value + b;
+      const double cc = 6 * a * x.value * x.value + 3 * b * x.value + c;
+      const double dd = 4 * a * x.value * x.value * x.value +
+                        3 * b * x.value * x.value + 2 * c * x.value + d;
+                        
+      return utility::FindMinimumIntegerQuartic(aa, bb, cc, dd, dxl, dxu, x.value, this->random_number_engine);
     } else {
       double min_dE = std::numeric_limits<double>::infinity();
       std::int64_t min_value = -1;
@@ -281,12 +287,24 @@ class IntegerSASystem<graph::IntegerPolynomialModel, RandType> {
 
   double GetEnergy() const { return this->energy_; }
 
-  bool OnlyMultiLinearCoeff(std::int64_t index) const {
-    return this->model.GetOnlyMultilinearIndexSet().count(index) > 0;
+  bool IsLinearCoeff(std::int64_t index) const {
+    return this->model.GetEachVariableDegreeAt(index) == 1;
   }
 
   bool UnderQuadraticCoeff(std::int64_t index) const {
-    return this->model.GetUnderQuadraticIndexSet().count(index) > 0;
+    return (this->model.GetEachVariableDegreeAt(index) == 2) || (this->model.GetEachVariableDegreeAt(index) == 1);
+  }
+
+  bool IsCubicCoeff(std::int64_t index) const {
+    return this->model.GetEachVariableDegreeAt(index) == 3;
+  }
+
+  bool IsQuarticCoeff(std::int64_t index) const {
+    return this->model.GetEachVariableDegreeAt(index) == 4;
+  }
+
+  bool CanOptMove(std::int64_t index) const {
+    return this->UnderQuadraticCoeff(index) || this->IsCubicCoeff(index) || this->IsQuarticCoeff(index);
   }
 
 public:

include/openjij/system/integer_quadratic_sa_system.hpp

@@ -173,11 +173,11 @@ class IntegerSASystem<graph::IntegerQuadraticModel, RandType> {
 
   double GetEnergy() const { return this->energy_; }
 
-  bool OnlyMultiLinearCoeff(std::int64_t index) const {
+  bool IsLinearCoeff(std::int64_t index) const {
     return this->model.GetOnlyBilinearIndexSet().count(index) > 0;
   }
 
-  bool UnderQuadraticCoeff(std::int64_t index) const { return true; }
+  bool CanOptMove(std::int64_t index) const { return true; }
 
   double GetLinearCoeff(std::int64_t index) const {
     return this->linear_coeff_[index];

include/openjij/updater/single_integer_move.hpp

@@ -39,8 +39,8 @@ struct OptMetropolisUpdater {
   std::int64_t GenerateNewValue(SystemType &sa_system, const std::int64_t index,
                                 const double T, const double progress) {
     // Metropolis Optimal Transition if possible
-    // This is used for systems with only quadratic coefficientsa
-    if (sa_system.UnderQuadraticCoeff(index) && dist(sa_system.random_number_engine) < progress) {
+    // This is used for systems with up to 4th power coefficients
+    if (sa_system.CanOptMove(index) && dist(sa_system.random_number_engine) < progress) {
       const auto [min_val, min_dE] = sa_system.GetMinEnergyDifference(index);
       if (min_dE <= 0.0 ||
           dist(sa_system.random_number_engine) < std::exp(-min_dE / T)) {
@@ -67,7 +67,7 @@ struct HeatBathUpdater {
   template <typename SystemType>
   std::int64_t GenerateNewValue(SystemType &sa_system, const std::int64_t index,
                                 const double T, const double _progress) {
-    if (sa_system.OnlyMultiLinearCoeff(index)) {
+    if (sa_system.IsLinearCoeff(index)) {
       return ForBilinear(sa_system, index, T, _progress);
     } else {
       return ForAll(sa_system, index, T, _progress);