<jupyter_start><jupyter_text>Experiments  PreliminaryIn the following, let $G=(V, E)$ be a graph, $P$ denotes the transition probability matrix for each Markov chain, $\pi$ denotes the stationary distribution with respect to $P$, and $\sigma, \tau, \eta\in\{-1, +1\}^V$.- Hamiltonian on the spin configuration space $\{\pm 1\}^{V}$$$    H(\sigma)    = -\sum_{\{x, y\}\in E}J_{x, y} \sigma_x \sigma_y - \sum_{x\in V} h_x \sigma_x    = -\frac{1}{2}\sum_{x, y\in V}J_{x, y} \sigma_x \sigma_y - \sum_{x\in V} h_x \sigma_x$$- Hamiltonian on $(\{\pm 1\}^{V^\mathrm{L}}, \{\pm 1\}^{V^\mathrm{R}})$ of the bipartite graph $(G^\mathrm{L}, G^\mathrm{R})$$$    \tilde H(\tau, \sigma) = -\frac{1}{2}\sum_{x, y\in V}J_{x, y} \tau_x \sigma_y - \frac{1}{2}\sum_{x\in V} h_x \left(\tau_x + \sigma_x\right) + \frac{1}{2}\sum_{x\in V}q_x \left( 1 - \tau_x \sigma_x\right)$$- Cavity field$$    \tilde h_x(\sigma) = \sum_{y\in V}J_{x, y} \sigma_y + h_x$$  Importing libraris<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
#import sys
#sys.path.append('../python')
#import simulator
import simulatorWithCpp as simulator
import math

%matplotlib inline
np.set_printoptions(threshold=16, edgeitems=8)<jupyter_output><empty_output><jupyter_text>Constant numbers (parameters)We fix the parameters:- the number of nodes $= 256$,- the maximum monte carlo steps $= 1000$,- the number of trials of an annealing $= 200$,- the initial spin configuration given uniformly at random by either PCG64 in the Python library or MT19937_64 in the C++ library with the random seed $1024$.<jupyter_code>MaxSteps = int(4.e3)
MaxTrials = int(2.e2)
NumNodes = 256
SeedForConfiguration = 1024<jupyter_output><empty_output><jupyter_text>Subroutines to sample and draw a graphNote that `Energy` property means the original Hamiltonian for the Metropolis method or the Glauber dynamics, whereas it means the Hamiltonian on the bipartite graph for the SCA, the MA or the MMA.<jupyter_code>isingModel = simulator.IsingModel({node: 0.e0 for node in range(NumNodes)}, {})
isingModel.SetSeed(SeedForConfiguration)
isingModel.GiveSpins(simulator.ConfigurationsType.Uniform)
InitialConfiguration = isingModel.Spins
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 TryExperimentFor(isingModel, initialTemperature, title='Ising model'):
    finalEnergies = np.empty(0, dtype=np.float)
    samples = np.empty((0, 3), dtype=np.float)
    for i in range(MaxTrials):
        isingModel.Spins = InitialConfiguration
        isingModel.SetSeed()
        if i == 0:
            isingModel.Temperature = initialTemperature
            isingModel.Write()
            print()
        for n in range(MaxSteps + 1):
            isingModel.Temperature = initialTemperature * 0.99 ** n
            #isingModel.Temperature = 100 * np.exp(-0.005 * n)
            isingModel.Update()
            if i == MaxTrials - 1:
                samples = np.append(samples, np.array([n, CalcEnergy(isingModel), isingModel.Temperature], dtype=np.float).reshape((1, 3)), axis=0)
        finalEnergies = np.append(finalEnergies, CalcEnergy(isingModel))
        
    print('Mean: {}'.format(np.mean(finalEnergies)))
    print('Standard deviation: {}'.format(np.std(finalEnergies)))
    print('Mode: {}'.format(stats.mode(finalEnergies)))
    print('Minimum: {}'.format(np.min(finalEnergies)))

    fig = plt.figure(figsize=(7, 3), dpi=200)
    ax = fig.add_subplot(121, xlabel='MC step', ylabel='Energy')
    ax.grid()
    ax.plot(samples[:, 0], samples[:, 1])
    ax = fig.add_subplot(122, xlabel='Energy', ylabel='Frequency')
    ax.grid(which='both')
    ax.hist(finalEnergies, bins=30)
    plt.subplots_adjust(wspace=0.3)
    plt.title(title)
    #plt.savefig(title + '.pdf', bbox_inches='tight', pad_inches=0)
    plt.show()<jupyter_output><empty_output><jupyter_text>Algorithms  Metropolis method- Stationary distribution$$    \pi(\sigma) = \frac{\exp\bigl(-\beta H(\sigma)\bigr)}{\sum_{\sigma'}\exp\bigl(-\beta H(\sigma')\bigr)}.$$- Transition probability$$    P(\tau \mid \sigma) = \begin{cases}        \displaystyle \frac{1}{\lvert V\rvert} \min\left\{1, \exp\bigl(-\beta \Delta E\bigr)\right\} & [\tau = \sigma^x],\\        \displaystyle 1 - \sum_{x\in V} P_{\sigma, \sigma^x} & [\tau = \sigma],\\        0 & [\text{otherwise}],    \end{cases}$$where $\Delta E \equiv H(\sigma^x) - H(\sigma) = 2 \tilde h_x(\sigma) \sigma_x$.- Then, the detailed balance condition$$    \pi(\sigma) P(\tau \mid \sigma) = \pi(\tau) P(\sigma \mid \tau)$$is satisfied.  Glauber dynamics- Stationary distribution$$    \pi(\sigma) = \frac{\exp\bigl(-\beta H(\sigma)\bigr)}{\sum_{\sigma'}\exp\bigl(-\beta H(\sigma')\bigr)}.$$- Transition probability$$    P(\tau \mid \sigma) = \begin{cases}        \displaystyle \frac{1}{\lvert V\rvert}\frac{\exp\bigl(-\beta \tilde h_x(\sigma) \sigma_x\bigr)}{2 \cosh\bigl(\beta \tilde h_x(\sigma)\bigr)} & [\tau = \sigma^x],\\        \displaystyle 1 - \sum_{x\in V} P_{\sigma, \sigma^x} & [\tau = \sigma],\\        0 & [\text{otherwise}].    \end{cases}$$- Then, the detailed balance condition$$    \pi(\sigma) P(\tau \mid \sigma) = \pi(\tau) P(\sigma \mid \tau)$$is satisfied.  Stochastic cellular automata- Stationary distribution on the bipartite graph $(G^\mathrm{L}, G^\mathrm{R})$$$    \pi(\tau, \sigma) = \frac{\exp\bigl(-\beta\tilde H(\tau, \sigma)\bigr)}{\sum_{\tau', \sigma'}\exp\bigl(-\beta\tilde H(\tau', \sigma')\bigr)}.$$- Transition probability (Glauber dynamics on the bipartite graph $(G^\mathrm{L}, G^\mathrm{R})$)\begin{align*}    P(\tau, \sigma \mid \eta, \sigma) &= \frac{\exp\bigl(-\beta\tilde H(\tau, \sigma)\bigr)}{\sum_{\tau'}\exp\bigl(-\beta\tilde H(\tau', \sigma)\bigr)}\\    &= \prod_{x\in V} \frac{\exp\bigl(-\beta (\tilde h_x(\sigma) + q_x\sigma_x) \tau_x / 2\bigr)}{2\cosh\bigl(\beta (\tilde h_x(\sigma) + q_x\sigma_x) / 2\bigr)}\\    &= \prod_{x\in V} \frac{1}{1 + \exp\bigl(-\beta (\tilde h_x(\sigma) + q_x\sigma_x) \tau_x\bigr)}.\end{align*}- Then, the detailed balance condition$$    \pi(\eta, \sigma) P(\tau, \sigma \mid \eta, \sigma) = \pi(\tau, \sigma) P(\eta, \sigma \mid \tau, \sigma)$$is satisfied.- Let $L_x\sim\mathrm{Logistic}(0, 1)$ for each $x\in V$ and $T = 1 / \beta$.  We regard $\sigma$ as the current spin configuration and $\tau$ as the next spin configuration in the MCMC updating.  Then, the algorithm is represented by$$    \DeclareMathOperator{\sgn}{sgn}    \tau_x = \sgn\left(\tilde h_x(\sigma) + q_x\sigma_x - TL_x\right).$$  Momentum annealing(In this section, we do not control temperature, so that this MCMC algorithm is not "annealing".  However, since the MCMC do not have any name, we call it the "momentum annealing" for convinience)- Stationary distribution on the bipartite graph $(G^\mathrm{L}, G^\mathrm{R})$$$    \pi(\tau, \sigma) = \frac{\exp\bigl(-\beta\tilde H(\tau, \sigma)\bigr)}{\sum_{\tau', \sigma'}\exp\bigl(-\beta\tilde H(\tau', \sigma')\bigr)}.$$- When we update the spins on $V^\mathrm{L}$, the energy difference is expressed by\begin{align*}    \Delta E &\equiv \tilde H(\tau, \sigma) - \tilde H(\eta, \sigma)\\    &= \frac{1}{2} \sum_{x\in V} \left(\tilde h_x(\sigma) + q_x \sigma_x\right) \left(\eta_x - \tau_x\right)\\    &= \sum_{\substack{x\in V\\ (\tau_x = -\eta_x)}} \left(\tilde h_x(\sigma) + q_x \sigma_x\right) \eta_x\end{align*}for $\eta, \tau\in\{-1, +1\}^{V^\mathrm{L}}$ and $\sigma\in\{-1, +1\}^{V^\mathrm{R}}$.  Let $\Delta E_x = \left(\tilde h_x(\sigma) + q_x \sigma_x\right) \eta_x$.- Transition probability (Metropolis method on the bipartite graph $(G^\mathrm{L}, G^\mathrm{R})$)\begin{align*}    P(\tau, \sigma \mid \eta, \sigma)    &= \prod_{x\in V} \left(\min\left\{1, \mathrm{e}^{-\beta \Delta E_x}\right\} \delta_{\tau_x, -\eta_x} + \left( 1 - \min\left\{1, \mathrm{e}^{-\beta \Delta E_x}\right\}\right) \delta_{\tau_x, \eta_x}\right)\\    &= \left(\prod_{\substack{x\in V\\ (\tau_x = -\eta_x)}} \min\left\{1, \mathrm{e}^{-\beta \Delta E_x}\right\}\right)        \left(\prod_{\substack{y\in V\\ (\tau_y = \eta_y)}} \left( 1 - \min\left\{1, \mathrm{e}^{-\beta \Delta E_y}\right\}\right)\right).\end{align*}- Then, the detailed balance condition$$    \pi(\eta, \sigma) P(\tau, \sigma \mid \eta, \sigma) = \pi(\tau, \sigma) P(\eta, \sigma \mid \tau, \sigma)$$is satisfied.- Let $\varGamma_x\sim\mathrm{Exponential}(1)$ for each $x\in V$ and $T = 1 / \beta$.  We regard $\sigma$ as the current spin configuration, $\tau$ as the next spin configuration and $\eta$ as the previous spin configuration in the MCMC updating.  Then, the algorithm is represented by$$    \DeclareMathOperator{\sgn}{sgn}    \tau_x = \sgn\left(\tilde h_x(\sigma) + q_x\sigma_x - T\varGamma_x\eta_x\right).$$  Annealing scheduleFor $n=\mathbb{N}\cup\{0\}$,$$    T_n = T_0 \cdot 0.99^n.    %T_n = 100 \cdot \mathrm{e}^{-0.005 n}.$$  Square latticeOur experiment runs under the conditions:- no external magnetic field,- antiferromagnet,- the free boundary condition.<jupyter_code>def GenerateSquareLatticeEdges(numNodes):
    columns = math.ceil(math.sqrt(numNodes))
    result = {}
    for i in range(numNodes - 1):
        if (i + 1) % columns > 0:
            result[(i, i + 1)] = -1
        if (i + columns) < numNodes:
            result[(i, i + columns)] = -1
    return result
# Test for the function GenerateSquareLatticeEdges.
for i in range(16):
    if i != 0 and i % 4 == 0:
        print()
    print('{:3d}'.format(i), end='')
print('\n')
quadratic = GenerateSquareLatticeEdges(16)
print(quadratic)
print('\n')
isingModel = simulator.IsingModel({}, quadratic)
isingModel.Write()
# Actual simulation.
quadratic = GenerateSquareLatticeEdges(NumNodes)
isingModel = simulator.IsingModel({}, quadratic)
T0 = 2.e0 * np.sum([np.abs(J) for J in quadratic.values()])<jupyter_output><empty_output><jupyter_text>Metropolis method<jupyter_code>isingModel.Algorithm = simulator.Algorithms.Metropolis
TryExperimentFor(isingModel, T0, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Glauber dynamics<jupyter_code>isingModel.Algorithm = simulator.Algorithms.Glauber
TryExperimentFor(isingModel, T0, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Stochastic cellular automata<jupyter_code>isingModel.Algorithm = simulator.Algorithms.SCA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Momentum annealing<jupyter_code>isingModel.Algorithm = simulator.Algorithms.MA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Complete graphOur experiment runs under the conditions:- no external magnetic field,- antiferromagnet.<jupyter_code>def GenerateCompleteGraphEdges(numNodes):
    return {(i, j): -1 for i in range(numNodes) for j in range(i + 1, numNodes)}

quadratic = GenerateCompleteGraphEdges(NumNodes)
isingModel = simulator.IsingModel({}, quadratic)
T0 = 2.e0 * np.sum([np.abs(J) for J in quadratic.values()])<jupyter_output><empty_output><jupyter_text>Metropolis method<jupyter_code>isingModel.Algorithm = simulator.Algorithms.Metropolis
TryExperimentFor(isingModel, T0, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Glauber dynamics<jupyter_code>isingModel.Algorithm = simulator.Algorithms.Glauber
TryExperimentFor(isingModel, T0, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Stochastic cellular automata<jupyter_code>isingModel.Algorithm = simulator.Algorithms.SCA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Momentum annealing<jupyter_code>isingModel.Algorithm = simulator.Algorithms.MA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Erdős-Rényi random graphOur experiment runs under the conditions:- no external magnetic field,- antiferromagnet,- the occupation probability $=0.5$,- Bernoulli random numbers given by MT19937 with the random seed $2048$.<jupyter_code>OccupationProbability = 0.5e0
SeedForRandomGraph = 2048

def GenerateErdosRenyiEdges(numNodes, probability):
    rng = np.random.Generator(np.random.MT19937(SeedForRandomGraph))
    return {(i, j): -1 if rng.random() <= probability else 0 for i in range(numNodes) for j in range(i + 1, numNodes)}

quadratic = GenerateErdosRenyiEdges(NumNodes, OccupationProbability)
isingModel = simulator.IsingModel({}, quadratic)
T0 = 2.e0 * np.sum([np.abs(J) for J in quadratic.values()])<jupyter_output><empty_output><jupyter_text>Metropolis method<jupyter_code>isingModel.Algorithm = simulator.Algorithms.Metropolis
TryExperimentFor(isingModel, T0, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Glauber dynamics<jupyter_code>isingModel.Algorithm = simulator.Algorithms.Glauber
TryExperimentFor(isingModel, T0, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Stochastic cellular automata<jupyter_code>isingModel.Algorithm = simulator.Algorithms.SCA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output><jupyter_text>Momentum annealing<jupyter_code>isingModel.Algorithm = simulator.Algorithms.MA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter, title=simulator.AlgorithmToStr(isingModel.Algorithm))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1602529235790046, 'fraction_numerical': 0.014549904813706827, 'mean_word_length': 4.6270084162203515, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 74, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '466101', 'license': 'MIT', 'n_tokens_mistral': 5247, 'n_tokens_neox': 4931, 'n_words': 1103, 'path': 'docs/experiments.ipynb', 'repo_name': 'Wandao123/ising_model', 'size': '14708'}