# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Flip-contrained SCA
#
# For SCA dynamics, we plot a sample of the modified Hamiltonian in "Step-Energy graph".  For flip-constrained SCA dynamics, we plot a sample of the (original) Hamiltonian in "Step-Energy graph".

# +
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)
# -

# ## Annealing

# +
maxSteps = int(3.e3)
MaxTrials = int(1.e2)
NumNodes = 256
SeedForConfiguration = 1024

def MakeUniformConfiguration(numNodes, seed=None):
    rng = np.random.default_rng(seed)
    return {i: rng.choice([-1, +1]) for i in range(numNodes)}

InitialConfiguration = MakeUniformConfiguration(NumNodes, SeedForConfiguration)

def TryExperimentFor(isingModel, initialTemperature):
    minimumEnergiesData = np.empty(0, dtype=np.float)
    samples = np.empty((maxSteps + 1, 4), dtype=np.float)
    #isingModel.Write()
    for i in range(MaxTrials):
        isingModel.Spins = InitialConfiguration
        isingModel.SetSeed()
        for n in range(maxSteps + 1):
            isingModel.Temperature = initialTemperature * 0.99 ** n
            #isingModel.Temperature = 100 * np.exp(-0.005 * n)
            isingModel.Update()
            samples[n, 0] = n
            samples[n, 1] = isingModel.Temperature
            samples[n, 2] = isingModel.Energy
            samples[n, 3] = isingModel.EnergyOnBipartiteGraph
        minimumEnergiesData = np.append(minimumEnergiesData, samples[:, 2].min())

    print('Mean: {}'.format(np.mean(minimumEnergiesData)))
    print('Standard deviation: {}'.format(np.std(minimumEnergiesData)))
    print('Mode: {}'.format(stats.mode(minimumEnergiesData)))
    print('Minimum: {}'.format(np.min(minimumEnergiesData)))

    fig = plt.figure(figsize=(7, 3), dpi=200)
    ax = fig.add_subplot(121, xlabel='MC steps', ylabel='Energy')
    ax.grid()
    ax.plot(samples[:, 0], samples[:, 2], label='Original Hamiltonian')
    ax.plot(samples[:, 0], samples[:, 3], label='Double Hamiltonian')
    ax.legend()
    ax = fig.add_subplot(122, xlabel='Energy', ylabel='Frequency')
    ax.grid(which='both')
    ax.hist(minimumEnergiesData, bins=30)
    fig.suptitle(isingModel.Algorithm.name)
    plt.subplots_adjust(wspace=0.3)
    plt.show()


# -

# ### Square lattice

# +
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

quadratic = GenerateSquareLatticeEdges(NumNodes)
isingModel = simulator.IsingModel({}, quadratic)
T0 = 2.e0 * np.sum([np.abs(J) for J in quadratic.values()])
# -

isingModel.Algorithm = simulator.Algorithms.SCA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter)

isingModel.Algorithm = simulator.Algorithms.fcSCA
isingModel.PinningParameter = 0.125e0 * isingModel.CalcLargestEigenvalue()
isingModel.FlipTrialRate = 0.8e0
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter)


# ### Complete graph

# +
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()])
# -

isingModel.Algorithm = simulator.Algorithms.SCA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter)

isingModel.Algorithm = simulator.Algorithms.fcSCA
isingModel.PinningParameter = 0.e0
isingModel.FlipTrialRate = 0.3e0
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter)

# ### Erdős-Rényi random graph

# +
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()])
# -

isingModel.Algorithm = simulator.Algorithms.SCA
isingModel.PinningParameter = 0.5e0 * isingModel.CalcLargestEigenvalue()
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter)

isingModel.Algorithm = simulator.Algorithms.fcSCA
isingModel.PinningParameter = 0.e0
isingModel.FlipTrialRate = 0.2e0
TryExperimentFor(isingModel, T0 + NumNodes * isingModel.PinningParameter)


# ## Stationary distribution

# ## Dynamics for the propotion to spins

# Consider the antiferromagnet $J_{x, y} = -1$ with $h_x = 0$ and $q_x = 0$ on a complete graph.
# Taking the low temerature limit,
# $$
#     \lim_{\beta\uparrow\infty} P_\epsilon(\sigma, \tau)
#     = \prod_{x\in D_{\sigma, \tau}} \left(\epsilon \left(\frac{1}{2} \mathbf{1}_{\{\sum_{y\sim x} \sigma_x \sigma_y = 0\}} + \mathbf{1}_{\{\sum_{y\sim x} \sigma_x \sigma_y > 0\}}\right)\right)
#         \prod_{x\in D_{\sigma, \tau}^\complement} \left(1 - \epsilon \left(\frac{1}{2} \mathbf{1}_{\{\sum_{y\sim x} \sigma_x \sigma_y = 0\}} + \mathbf{1}_{\{\sum_{y\sim x} \sigma_x \sigma_y > 0\}}\right)\right).
# $$
# Let $a_t$ be the proportion of up spins to all spins.
# When $\lvert V\rvert$ is even, he dynamics of $\{a_t\}_{t=0}^{\infty}$ is given as
# $$
#     a_{t+1} = \begin{cases}
#         \epsilon + \left(1 - \epsilon\right) a_t & \left[0 \leq a_t < \frac{1}{2}\right],\\
#         a_t & \left[a_t = \frac{1}{2}\right],\\
#         \left(1 - \epsilon\right) a_t & \left[\frac{1}{2} < a_t \leq 1\right].
#     \end{cases}
# $$
# When $\lvert V\rvert$ is odd, he dynamics of $\{a_t\}_{t=0}^{\infty}$ is given as
# $$
#     a_{t+1} = \begin{cases}
#         \left(1 - \frac{\epsilon}{2}\right) a_t & \left[\frac{1}{2} \leq a_t \leq \frac{1}{2} \left(1 + \frac{1}{\lvert V\rvert}\right)\right],\\
#         \left(1 - \epsilon\right) a_t & \left[\frac{1}{2} \left(1 + \frac{1}{\lvert V\rvert}\right) < a_t \leq 1\right],\\
#         \frac{\epsilon}{2} + \left(1 - \frac{\epsilon}{2}\right) a_t & \left[\frac{1}{2} \left(1 - \frac{1}{\lvert V\rvert}\right) \leq a_t < \frac{1}{2}\right],\\
#         \epsilon + \left(1 - \epsilon\right) a_t & \left[0 \leq a_t < \frac{1}{2} \left(1 - \frac{1}{\lvert V\rvert}\right)\right].
#     \end{cases}
# $$

def spinPropotion(previous, flipTrialRate, cardVertices):
    if cardVertices % 2 == 0:
        if 0 <= previous < 0.5e0:
            return flipTrialRate + (1 - flipTrialRate) * previous
        elif previous == 0.5e0:
            return previous
        elif 0.5e0 < previous <= 1:
            return (1 - flipTrialRate) * previous
        else:
            raise ValueError('The value must be in the interval [0, 1]')
    else:
        if 0 <= previous < 0.5e0 * (1 - 1 / cardVertices):
            return flipTrialRate + (1 - flipTrialRate) * previous
        elif 0.5e0 * (1 - 1 / cardVertices) <= previous < 0.5e0:
            return 0.5e0 * flipTrialRate + (1 - 0.5e0 * flipTrialRate) * previous
        elif 0.5e0 <= previous <= 0.5e0 * (1 + 1 / cardVertices):
            return (1 - 0.5e0 * flipTrialRate) * previous
        elif 0.5e0 * (1 + 1 / cardVertices) < previous <= 1:
            return (1 - flipTrialRate) * previous
        else:
            raise ValueError('The value must be in the interval [0, 1]')


# The horizontal axis means initial values $a_0$.

# +
# Fixed points.
from itertools import product

NumDivision = 100
numDivFTR = 5

fig = plt.figure(figsize=(numDivFTR * 3 + 2, 7), dpi=100)
plt.subplots_adjust(wspace=0.3, hspace=0.3)
index = 0
for cardVertices, flipTrialRate in product((10, 11), np.linspace(0, 1, numDivFTR)):
    # Time evolution
    uSpinPropotion = np.frompyfunc(lambda x: spinPropotion(x, flipTrialRate, cardVertices), 1, 1)
    x = np.linspace(0, 1, NumDivision)
    y = uSpinPropotion(x)
    for _ in range(NumDivision):
        try:
            y = uSpinPropotion(y)
        except ValueError:
            print(y)
            break

    # Plotting
    index += 1
    ax = fig.add_subplot(2, numDivFTR, index)
    ax.grid()
    ax.plot(x, x, color='gray', linestyle='dashed')
    ax.plot(x, uSpinPropotion(x), label='a(1)')
    ax.plot(x, y, label='a(t)')
    nextY = uSpinPropotion(y)
    ax.plot(x, nextY, label='a(t+1)')
    ax.set_title('ε={:8.4f}'.format(flipTrialRate) + ', |V|={:d}'.format(cardVertices))
    ax.legend()
plt.show()

# +
# %%time
# Bifurcation diagram.
from functools import reduce

NumDivision = 400
maxSteps = 200

def calcFixedPoints(flipTrialRate, cardVertices):
    def inner(initialValue):
        result = []
        propotion = initialValue
        for i in range(maxSteps // 2):
            propotion = spinPropotion(propotion, flipTrialRate, cardVertices)
        for i in range(maxSteps // 2, maxSteps):
            propotion = spinPropotion(propotion, flipTrialRate, cardVertices)
            result.append(propotion)
        return np.array(result, dtype=np.float)
    
    fixedPoints = np.empty(0, dtype=np.float)
    initialValues = np.linspace(0, 1, NumDivision)
    fixedPoints = reduce(np.union1d, np.frompyfunc(inner, 1, 1)(initialValues))
    return fixedPoints

fig = plt.figure(figsize=(7, 3), dpi=200)
plt.subplots_adjust(wspace=0.3)
for index, cardVertices in enumerate((10, 11), 1):
    # Computation
    x = np.linspace(0, 1, NumDivision)
    y = np.frompyfunc(lambda r: calcFixedPoints(r, cardVertices), 1, 1)(x)  # This function returns a jagged array.
    
    # Plotting
    ax = fig.add_subplot(1, 2, index)
    ax.grid()
    ax.set_xlabel('Flip trial rate')
    ax.set_ylabel('Fixed point')
    ax.set_title('|V|={:d}'.format(cardVertices))
    for i in range(x.size):
        ax.plot(np.full_like(y[i], x[i]), y[i], ',m')
plt.show()


# -

# ## Antiferromagnet on a complete graph for the low temperature limit

# The following graphs are the system energies per monte carlo step for $\lvert V\rvert = 100$ and each $\epsilon\in [0, 1]$.

def DrawGraphsForAntiferrOnK(numNodes, maxSteps, numDivFTR):
    fig = plt.figure(figsize=(numDivFTR // 2 * 3 + 2, 13), dpi=100)
    plt.subplots_adjust(wspace=0.3, hspace=0.3)

    # Initialization
    isingModel = simulator.IsingModel({}, GenerateCompleteGraphEdges(numNodes))
    isingModel.Algorithm = simulator.Algorithms.fcSCA
    initialTemperature = 2.e0 * np.sum([np.abs(J) for J in quadratic.values()])
    initialConfiguration = MakeUniformConfiguration(numNodes)
    print('|V|={:d}'.format(numNodes))
    print('a(0)={:g}'.format(sum([spin for spin in initialConfiguration.values() if spin == +1]) / numNodes))

    for index, flipTrialRate in enumerate(np.linspace(0, 1, numDivFTR), 1):
        # Time evolution
        samples = np.empty((maxSteps + 1, 3), dtype=np.float)
        isingModel.Spins = initialConfiguration
        isingModel.FlipTrialRate = flipTrialRate
        isingModel.SetSeed()
        for n in range(maxSteps + 1):
            isingModel.Temperature = initialTemperature * 0.9 ** n
            isingModel.Update()
            samples[n, 0] = n
            samples[n, 1] = isingModel.Energy
            samples[n, 2] = sum([spin for spin in isingModel.Spins.values() if spin == +1]) / numNodes

        # Plotting
        ax = fig.add_subplot(4, numDivFTR // 2, index)
        ax.grid()
        ax.plot(samples[:, 0], samples[:, 1], label='Original Hamiltonian')
        ax.set_title('ε={:8.4f}'.format(flipTrialRate))
        ax.legend()
        ax = fig.add_subplot(4, numDivFTR // 2, numDivFTR + index)
        ax.grid()
        ax.plot(samples[:, 0], samples[:, 2], ',r', label='Up spin propotion')
        ax.set_title('ε={:8.4f}'.format(flipTrialRate))
        ax.legend()
    plt.show()


DrawGraphsForAntiferrOnK(100, 10000, 10)

# The following graphs are the system energies per monte carlo step for $\lvert V\rvert = 101$ and each $\epsilon\in [0, 1]$.

DrawGraphsForAntiferrOnK(101, 10000, 10)



StarcoderdataJupyterScriptsDedupFiltered