A3 - 古典イジング模型¶
この章ではOpenJijのcore interface (core python interface)の使い方を説明し、簡単な計算のデモンストレーションを行います。
core interfaceは前回までのチュートリアルよりも下部のレイヤーのAPIです。対象読者としては前回までのOpenJijチュートリアルを一通り終えて、イジングモデルやモンテカルロ法などの用語を知っている方を想定します。次のような目的を持った読者を想定しています。
最適化問題だけでなくサンプリングや研究用途などより専門的な用途にOpenJijを用いたい
アニーリングスケジュールの設定や使用するアルゴリズム等を直接触りたい
[1]:
import cxxjij.graph as G
#問題サイズを100とします。
N = 100
graph = G.Dense(N)
#sparseの場合は以下を使用します。
#graph = G.Sparse(N)
\(J_{ij}, h_i\)を設定します。今回は平均0、標準偏差1のGauss分布から生成される値を設定します。
[ ]:
!pip install numpy #乱数生成にnumpyを使います。
[2]:
import numpy as np
mu, sigma = 0, 1
for i in range(N):
for j in range(N):
#Jijの値が大きくなりすぎてしまうので、1/Nで規格化を行なっています。
graph[i,j] = 0 if i == j else np.random.normal()/N
for i in range(N):
graph[i] = np.random.normal()/N
縦磁場に関しては、graph[i]でも、graph[i,i]でもどちらでもアクセスできます。また、イジングモデルの定義上、\(J_{ij}\)と\(J_{ji}\)は自動で同じ値となります。試しに以下のように出力を行なってみましょう。
[3]:
graph[20] = 0.5
print(graph[20,20])
print(graph[20])
graph[12,34] = -0.6
print(graph[12,34])
print(graph[34,12])
0.5
0.5
-0.6
-0.6
システムの設定 - System -¶
続いて計算を行うためのシステムを定義します。ここでは古典イジング模型か横磁場イジング模型か、また別の模型にするか等を選べます。
まずは古典イジング模型のシステムを作成してみます。system.make_classical_isingで作成できます。
[4]:
import cxxjij.system as S
mysystem = S.make_classical_ising(graph.gen_spin(), graph)
ここで、1つ目の引数にはランダムに生成したスピン、2つめにはGraphそのものを代入します。これにより初期のスピン配位がgraph.gen_spin()となる古典イジング模型のシステムの作成ができます。
システムに直接アクセスして、値を読むことも可能です。
[5]:
print(mysystem.spin)
[-1. 1. -1. -1. 1. -1. -1. 1. -1. 1. 1. 1. -1. -1. 1. -1. 1. -1.
-1. 1. -1. 1. -1. -1. 1. 1. 1. -1. 1. -1. -1. -1. 1. -1. 1. -1.
1. -1. 1. -1. 1. -1. 1. 1. -1. -1. -1. 1. -1. 1. 1. -1. 1. 1.
-1. -1. -1. 1. 1. 1. -1. 1. 1. 1. -1. -1. -1. -1. -1. 1. 1. -1.
-1. 1. 1. 1. 1. -1. 1. 1. 1. 1. 1. -1. -1. 1. -1. -1. -1. -1.
1. 1. -1. 1. -1. -1. -1. -1. -1. -1. 1.]
古典イジング模型以外にも様々なシステムが用意されており、これを用途別に使うことが出来ます。用いるSystemによって初期化の方法は多少異なります。これは後ほど少しずつ紹介していきます。
アルゴリズムの実行 -Updater, Algorithm-¶
Systemを定義した後はUpdaterを選択し、Algorithmを実行していきます。
Updater¶
Systemに対して使用できるUpdaterは決められており、古典イジング模型に対するUpdaterは主に
SingleSpinFlip (メトロポリス・ヘイスティング法によるスピン1つずつのアップデート)
SwendsenWang (SwendsenWang法によるクラスターアップデート)
の2つが用意されています。 Algorithmを実行するには、スケジュールリストが必要になります。まずスケジュールリストを作成するところから始めましょう。
Algorithm¶
スケジュールリスト¶
スケジュールリストは(パラメータ, モンテカルロステップ数)のリストで与えられるものです。パラメータ部分に入力する値はSystemによって異なります。例えば、古典イジング模型ならばパラメータとして温度の逆数である逆温度\(\beta\)を設定します。例として以下のように設定してみましょう。
[6]:
schedule_list = [(0.01, 10),(10, 80),(0.1, 30)]
utilityにあるmake_classical_schedule_listを使うとより便利です。[7]:
import cxxjij.utility as U
schedule_list = U.make_classical_schedule_list(0.1, 50, 20, 10)
print(schedule_list)
[((beta: 0.100000) mcs: 20), ((beta: 0.199474) mcs: 20), ((beta: 0.397897) mcs: 20), ((beta: 0.793701) mcs: 20), ((beta: 1.583223) mcs: 20), ((beta: 3.158114) mcs: 20), ((beta: 6.299605) mcs: 20), ((beta: 12.566053) mcs: 20), ((beta: 25.065966) mcs: 20), ((beta: 50.000000) mcs: 20)]
上の例では\(\beta=0.1\)から\(\beta=50\)まで、各温度で20モンテカルロステップ計算しながら10段階で温度を変えていく設定例です。計200モンテカルロステップの計算を行います。
Algorithmの実行¶
続いて、Algorithmを実行します。Algorithm_[Updater]_runのように書くことで、指定したUpdaterで計算を行います。次の例ではSingleSpinFlipを実行しています。
[8]:
import cxxjij.algorithm as A
A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list)
一瞬で処理が終わりましたが、この間に計200モンテカルロステップの計算が行われています。 > A.Algorithm_SingleSpinFlip_run(mysystem, seed, schedule_list)とすることで、seedを固定したまま計算を行うことができます。結果に再現性をもたせたい際に使うことができます。
[9]:
energies = []
def callback_log_energy(system, beta):
#graphは以前にGraphモジュールにて定義したオブジェクトです
energies.append(graph.calc_energy(system.spin))
このcallbackを用いて同じAlgorithmを実行します。
[10]:
#スケジュールをもっと長く取ります (計20000モンテカルロステップ)
schedule_list = U.make_classical_schedule_list(0.1, 50, 200, 100)
A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list, callback_log_energy)
記録したシステムのエネルギーを、横軸をモンテカルロステップ、縦軸をエネルギーでプロットすると次のようになります。
[ ]:
!pip install matplotlib
[13]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(range(len(energies)), energies)
plt.xlabel('Monte Carlo step')
plt.ylabel('energy')
plt.show()
アニーリングが進むに連れ徐々にエネルギーが低くなっているのが分かります。Algorithmの動作中にシステムの様子を知りたい時に有用です。
結果の取得 -Result-¶
result.get_solutionで計算結果であるスピン列を取得できます。古典イジング模型の場合は直接mysystem.spinを参照することで、スピン列を取得も可能です。しかし、result.get_solutionはそれ以外のシステムについてもスピン列を取得できる便利なメソッドです。
[14]:
import cxxjij.result as R
print(R.get_solution(mysystem))
[-1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, -1]
このスピン列がアニーリングによって得られた答えです。ハミルトニアンの基底状態 (に近い状態)であることが期待されます。
C++ core interface¶
多少の違いはありますが、C++ core interfaceでも上記とほぼ同じことが可能です。以下の点に留意しましょう。
seed値を入れる引数には、乱数生成器 (C++11 random)を代入する必要があります。
Graphクラスで、\(J_{ij}, h_i\)へのアクセス方法が多少異なります。
今までの内容をC++ core interfaceで記述すると次のようになります。
#include <graph/all.hpp>
#include <system/all.hpp>
#include <updater/all.hpp>
#include <algorithm/all.hpp>
#include <result/all.hpp>
#include <utility/schedule_list.hpp>
#include <utility/random.hpp>
#include <random>
#include <iostream>
using namespace openjij;
int main(void){
//generate dense graph with size N=100
constexpr std::size_t N = 100;
auto dense = graph::Dense<double>(N);
//generate random engine
auto rand_engine = std::mt19937(0x1234);
//of course you can specify another random generator that is compatible with C++ random generator, say utility::Xorshift,
//auto rand_engine = utility::Xorshift(0x1234);
//Gaussian distribution
auto gauss = std::normal_distribution<>{0, 1};
//set interactions
for(std::size_t i=0; i<N; i++){
for(std::size_t j=0; j<N; j++){
dense.J(i, j) = (i == j) ? 0 : gauss(rand_engine)/N;
}
}
//set local fields
for(std::size_t i=0; i<N; i++){
dense.h(i) = gauss(rand_engine);
}
//create classical Ising system
auto system = system::make_classical_ising(dense.gen_spin(rand_engine), dense);
//generate schedule list
//from beta=0.1 to beta=50, 20 Monte Carlo step for each tempearature, 10 samples
auto schedule_list = utility::make_classical_schedule_list(0.1, 50, 20, 10);
//or in the case of classical system, you can specify arbitrary schedules by using a std::vector of utility::Schedule<system::classical_system>
//std::vector<utility::Schedule<system::classical_system>> schedule_list;
//utility::Schedule<system::classical_system> schedule;
//schedule.updater_parameter = {0.01};
//schedule.one_mc_step = 10; //number of Monte Carlo step per temperature
//schedule_list.push_back(schedule);
//
//schedule.updater_parameter = {10};
//schedule.one_mc_step = 80; //number of Monte Carlo step per temperature
//schedule_list.push_back(schedule);
//
//schedule.updater_parameter = {0.1};
//schedule.one_mc_step = 30; //number of Monte Carlo step per temperature
//schedule_list.push_back(schedule); //schedule_list -> [(0.01, 10), (10, 80), (0.1, 30)]
//do annealing (updater: SingleSpinFlip)
algorithm::Algorithm<updater::SingleSpinFlip>::run(system, rand_engine, schedule_list);
//show spins
std::cout << "The result spins are [";
for(auto&& elem : result::get_solution(system)){
std::cout << elem << " ";
}
std::cout << "]" << std::endl;
}