<jupyter_start><jupyter_text>Pythonによるデモ本稿ではPythonによるシミュレーターの使い方を示す。実際にはPythonのみで書かれた `simulator.py` ではなく、C++拡張ライブラリを用いる。前者は後者に比べて3倍以上遅いので、実装を確認する程度にすることを勧める。  準備  ダウンロードリポジトリのページの "Code" --> "Download ZIP" から圧縮ファイルをダウンロードし、それを解凍する。あるいは、`git clone https://github.com/Wandao123/ising_model.git` コマンドを実行する。![downloading](downloading.png)  Vanilla PythonPythonの処理系は[Windows Store](https://www.microsoft.com/ja-jp/p/python-38/9mssztt1n39l?activetab=pivot:overviewtab)や[公式ページ](https://www.python.org/)からダウンロード・インストールできる。ここでは仮想環境を作成し、そこへpipを用いて必要なライブラリをインストールする。PowerShellあるいはコマンドプロンプトを起動した上で、次のコマンドを実行する。```PowerShellPS> cd ダウンロードしたフォルダPS> cd pythonPS> python -m venv envPS> env/Scripts/activatePS> pip install numpy matplotlib jupyterlab```  Anaconda一般的な用途ではなく、データ分析や機械学習などに特化してPythonを用いる場合は[Anaconda](https://www.anaconda.com/)を利用する方法もある。Anacondaを利用する場合、必要なライブラリが既にインストールされている筈なので、前節の手順はほぼ不要である。しかしながら、独自ライブラリのインストールの関係上、仮想環境の作成は勧める。Anaconda PowerShell Promptを起動して、次を実行する。```PowerShellPS> conda create -n ising_modelPS> conda activate ising_modelPS> conda install pipPS> conda install jupyterlab   Jupyter Labを使う場合。```  C++拡張ライブラリのインストールリポジトリのページの "Releases"（「ダウンロード」節の画像を参照）からリリースのページへ移動する。そこで最新版のwheelパッケージをダウンロードする。![downloading-whl](downloading-whl.png)ダウンロードが完了したら、pipを用いてそのwheelパッケージをインストールする（次の手順では `env/Scripts/activate` あるいは `conda activate ising_model` で既に仮想環境が有効化されているとする）：```PowerShellPS> cd wheelパッケージをダウンロードしたフォルダPS> pip install simulatorWithCpp-*-*-*-*.whl```ただし、\* の部分には適当なバージョン名やアーキテクチャ名が入る。  デモンストレーション必要なライブラリの読み込みを行う。Pythonのみで書かれたライブラリとC++拡張ライブラリとを切り替えるには、使いたい方のコメントアウトを解除し、他方をコメントアウトすればよい。それぞれについて、下記のコードを全く変えることがなく実行できることに注意。<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import sys
sys.path.append('../python')
import simulator                      # Pythonのライブラリを使う場合。
#import simulatorWithCpp as simulator  # C++拡張ライブラリを使う場合。

# 各種設定。
%matplotlib inline
np.set_printoptions(threshold=16, edgeitems=8)<jupyter_output><empty_output><jupyter_text>頂点と辺の指定方法  3スピン系による例題グラフ $G=(V, E)$ について、それぞれの集合が$$    V = \{\,a, b, c\,\}, \quad E = \{\,\{a, b\}, \{b, c\}, \{a, c\}\,\}$$で与えられたとする。頂点 $a$ に $1$ の、頂点 $b$ に $-2$ の、頂点 $c$ に $2$ の、辺 $\{a, b\}$ に $-4$ の、辺 $\{b, c\}$ に $2$ の、辺 $\{a, c\}$ に $1$ のバイアスが印加されていたとする。このとき、Hamiltonian関数は$$    H(\sigma) = 4 \sigma_a \sigma_b - 2 \sigma_b \sigma_c - \sigma_a \sigma_c - \sigma_a + 2 \sigma_b - 2 \sigma_c$$となる。これを表すコードは次のようになる。ただし、`Write` メソッドはIsing模型の現在の状態を表示する。初期配置はデフォルトで全てアップ・スピンであるため、`Current spin configuration` で配列の要素が全て $1$ になっている。任意の配置を指定、あるいは現在の配置を取得するには `Spins` プロパティを使う。それに辞書型変数を代入することで `isingModel` インスタンスのフィールドが変更される。<jupyter_code>isingModel = simulator.IsingModel({'a': 1.e0, 'b': -2.e0, 'c': 2.e0}, {('a', 'b'): -4.e0, ('b', 'c'): 2.e0, ('a', 'c'): 1.e0})
isingModel.Write()
print()
isingModel.Spins = {'a': +1, 'b': -1, 'c': +1}
print(isingModel.Spins)<jupyter_output><empty_output><jupyter_text>全てのスピン配置に対するエネルギーを書き下すと次の表を得る。|Spin configuration $\sigma$     |$\sigma_a$|$\sigma_b$|$\sigma_b$|Energy $H(\sigma)$||--------------------------------|----------|----------|----------|------------------||$\uparrow\uparrow\uparrow$      |+1        |+1        |+1        |0                 ||$\uparrow\uparrow\downarrow$    |+1        |+1        |-1        |10                ||$\uparrow\downarrow\uparrow$    |+1        |-1        |+1        |-8                ||$\uparrow\downarrow\downarrow$  |+1        |-1        |-1        |-6                ||$\downarrow\uparrow\uparrow$    |-1        |+1        |+1        |-4                ||$\downarrow\uparrow\downarrow$  |-1        |+1        |-1        |2                 ||$\downarrow\downarrow\uparrow$  |-1        |-1        |+1        |4                 ||$\downarrow\downarrow\downarrow$|-1        |-1        |-1        |2                 |これは次のプログラムによっても確かめられる。`Energy` プロパティ（読み込み専用）によって、Hamiltonian関数の値（エネルギー）が取得できることに注意。<jupyter_code>import itertools

for configuration in reversed(list(itertools.product([-1, +1], repeat=3))):
    isingModel.Spins = dict(zip(('a', 'b', 'c'), configuration))
    print('spin configuration={0:>12s}, energy={1:3.0f}'.format(str(configuration), isingModel.Energy))<jupyter_output><empty_output><jupyter_text>仕様一般に、Hamiltonian関数が$$    H(\sigma) = -\sum_{\{x, y\}\in E}J_{x, y} \sigma_x \sigma_y - \sum_{x\in V} h_x \sigma_x$$と表されているとする。このシミュレーターでは、IsingModelクラスに外部磁場の強さ $\{h_x\}_{x\in V}$ とスピン-スピン結合係数 $\{J_{x, y}\}_{\{x,y\}\in E}$ とをPythonの辞書型変数の形式で渡す。外部磁場の強さは文字列あるいは整数をキーとする辞書である。また、スピン-スピン結合係数は文字列あるいは整数のタプルをキーとする辞書である。このタプルの要素は2つのみであり、かつ `(a, b)` に対して `a < b` でなければならない（然らざる場合は無視される）。何も指定しないときにはそこでの外部磁場の強さやスピン-スピン結合係数が $0$ であるものと解釈される。特に、空の辞書 `{}` を渡した場合は全ての頂点あるいは辺で $0$ になる。__例__&nbsp;(Erdős-Rényiランダムグラフ上の反強磁性Ising模型)&nbsp;辺の生成をプログラムに任せることで、ランダムグラフやスピングラスが生成できる。実行する度にスピン-スピン結合係数の行列が変化することを確認してみよ。<jupyter_code>maxNodes = 8
probability = 0.5e0
rng = np.random.default_rng()
quadratic = {(i, j): -1 if rng.random() <= probability else 0 for i in range(maxNodes) for j in range(i + 1, maxNodes)}
isingModel = simulator.IsingModel({}, quadratic)
isingModel.Write()<jupyter_output><empty_output><jupyter_text>__例__&nbsp;(自由境界正方格子上の強磁性Ising模型)&nbsp;正方格子の生成には次のようなコードを用いる。<jupyter_code>import math

maxNodes = 16
columns = math.ceil(math.sqrt(maxNodes))
quadratic = {}
for i in range(maxNodes - 1):
    if (i + 1) % columns > 0:
        quadratic[(i, i + 1)] = 1
    if (i + columns) < maxNodes:
        quadratic[(i, i + columns)] = 1
isingModel = simulator.IsingModel({}, quadratic)
isingModel.Write()<jupyter_output><empty_output><jupyter_text>Glauber力学によるアニーリング`IsingModel` クラスは `Temperature` プロパティを持つ。その値に応じて、`Update` メソッドによるスピンの更新確率が変化する。ここでは、更新アルゴリズムにGlauber力学を指定した上で、各モンテカルロ・ステップ $n$ で温度を下げてゆく。アニーリング・スケジュールには$$    T_n = \frac{T_0}{3 \log (n + 1) + 1}, \quad T_0 = 19$$を用いる。また、初期配置を一様ランダムにとる。再現性の確保のため、random seedを明示的に指定する。なお、更新には毎回違うものに到達したいため、改めてrandom seedを選び直していることにも注意。※C++拡張ライブラリを用いているときはbad allocation errorが発生することがある。その場合は上記の「必要なライブラリの読み込み」のセルを選択した上で "Run" --> "Run Selected Cell and All Below" をクリックする。<jupyter_code># 初期化。
isingModel = simulator.IsingModel({0: 1.e0, 1: -2.e0, 2: 2.e0}, {(0, 1): -4.e0, (1, 2): 2.e0, (0, 2): 1.e0})
isingModel.Algorithm = simulator.Algorithms.Glauber  # 更新アルゴリズムを指定。
T0 = 19.e0
isingModel.Temperature = T0  # 温度の設定。
rng = np.random.default_rng(32)
initialConfiguration = {i: rng.choice([-1, +1]) for i in range(3)}  # 初期配置を設定。
isingModel.Spins = initialConfiguration
isingModel.SetSeed()
isingModel.Write()
print()

# サンプリング。
#samples = np.empty((0, 3), dtype=np.float)  # Numpyを使う場合。若干複雑になるので、最後にPythonのリストを変換する方法も併記している。
samples = []
for n in range(2000 * 2):
    isingModel.Temperature = T0 / (3 * np.log(1 + n) + 1.e0)
    isingModel.Update()
    #samples = np.append(samples, np.array([n, isingModel.Energy, isingModel.Temperature], dtype=np.float).reshape((1, 3)), axis=0)
    samples.append([n, isingModel.Energy, isingModel.Temperature])

# データの表示。
output = np.array(samples, dtype=np.float)
print(output)
print()
print('The final configuration: {}'.format(isingModel.Spins))<jupyter_output><empty_output><jupyter_text>最後に出力されたデータの各列はそれぞれ、ステップ数、エネルギー、温度を表す。横軸をステップ数、縦軸をエネルギーとすると、次のグラフを得る。<jupyter_code>x = output[:, 0]  # 0列目を抽出。
y = output[:, 1]

fig = plt.figure(figsize=(4, 4), dpi=200)
ax = fig.add_subplot(111)
ax.plot(x, y)
plt.show()<jupyter_output><empty_output><jupyter_text>また、次のグラフは終状態のエネルギーについてのヒストグラムである。全てで $100$ 回シミュレーションを実行している。これにより、エネルギーが $-8$ の状態の実現確率が最も大きいことが判る。<jupyter_code>samples = []
for counter in range(100):
    isingModel.Spins = initialConfiguration
    for n in range(2000 * 2):
        isingModel.Temperature = T0 / (3 * np.log(1 + n) + 1.e0)
        isingModel.Update()
    samples.append(isingModel.Energy)

fig = plt.figure(figsize=(4, 4), dpi=200)
ax = fig.add_subplot(111)
ax.hist(samples)
plt.show()<jupyter_output><empty_output><jupyter_text>Metropolis法によるアニーリング以下、他のアルゴリズムを適用したときのグラフを描いてゆく。共通の処理を関数に纏めていることに注意。<jupyter_code># エネルギーを計算する関数。それぞれのアルゴリズム毎にHamiltonianを場合分けする。
def CalcEnergy(isingModel):
    if isingModel.Algorithm == simulator.Algorithms.Glauber or isingModel.Algorithm == simulator.Algorithms.Metropolis:
        return isingModel.Energy
    elif isingModel.Algorithm == simulator.Algorithms.SCA or isingModel.Algorithm == simulator.Algorithms.MA or isingModel.Algorithm == simulator.Algorithms.MMA:
        return isingModel.EnergyOnBipartiteGraph
    elif isingModel.Algorithm == simulator.Algorithms.fcSCA:
        return isingModel.EnergyOnBipartiteGraph
    else:
        raise ValueError('Illeagal choises')

def SampleAndDraw(isingModel, initialTemperature):
    isingModel.SetSeed()

    # モンテカルロ・ステップ毎のエネルギーの変化。
    isingModel.Spins = initialConfiguration
    energiesPerStep = np.empty((0, 2), dtype=np.float)
    for n in range(2000 * 2):
        isingModel.Temperature = initialTemperature / (3 * np.log(1 + n) + 1.e0)
        isingModel.Update()
        energiesPerStep = np.append(energiesPerStep, np.array([n, CalcEnergy(isingModel)], dtype=np.float).reshape((1, 2)), axis=0)

    # 終状態でのエネルギー。    
    finalEnergies = np.empty(0, dtype=np.float)
    for counter in range(100):
        isingModel.Spins = initialConfiguration
        for n in range(2000 * 2):
            isingModel.Temperature = initialTemperature / (3 * np.log(1 + n) + 1.e0)
            isingModel.Update()
        finalEnergies = np.append(finalEnergies, CalcEnergy(isingModel))

    # グラフの描画。
    fig = plt.figure(figsize=(4, 8), dpi=200)
    ax = fig.add_subplot(211)
    ax.plot(energiesPerStep[:, 0], energiesPerStep[:, 1])
    ax = fig.add_subplot(212)
    ax.hist(finalEnergies)
    plt.show()

isingModel.Algorithm = simulator.Algorithms.Metropolis
SampleAndDraw(isingModel, T0)<jupyter_output><empty_output><jupyter_text>SCAによるアニーリング同様のことをSCAで行うと次のようになる。初期温度を変えていることに注意。<jupyter_code>isingModel.Algorithm = simulator.Algorithms.SCA
isingModel.PinningParameter = isingModel.CalcLargestEigenvalue() / 2
SampleAndDraw(isingModel, T0 + 3 * isingModel.PinningParameter)
isingModel.Algorithm = simulator.Algorithms.fcSCA
isingModel.FlipTrialRate = 0.7
isingModel.PinningParameter = 0
SampleAndDraw(isingModel, T0 + 3 * isingModel.PinningParameter)<jupyter_output><empty_output><jupyter_text>MAによるアニーリング<jupyter_code>isingModel.Algorithm = simulator.Algorithms.MA
SampleAndDraw(isingModel, T0 + 3 * isingModel.PinningParameter)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1519992254816536, 'fraction_numerical': 0.023235550392099912, 'mean_word_length': 4.826283135927806, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 60, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '446545', 'license': 'MIT', 'n_tokens_mistral': 5500, 'n_tokens_neox': 5143, 'n_words': 608, 'path': 'docs/demo.ipynb', 'repo_name': 'Wandao123/ising_model', 'size': '10329'}