#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Author: 冯振华
# StudentNumber: 202221140015
# Verdion: V4.0
# Copyright © 2023 feng <feng@arch>
#
# Distributed under terms of the MIT license.

import numpy as np
import matplotlib.pyplot as plt

# Define the parameters
N = 20 # number of spins in each direction
J = -1 # coupling constant
beta = 1 # inverse temperature
nsteps = 1000 # number of Monte Carlo steps

# Initialize the spins randomly
spins = np.random.choice([-1, 1], size=(N, N))

# Define the energy function
def energy(spins):
    return -J * (np.sum(spins[:-1,:] * spins[1:,:]) + np.sum(spins[:,:-1] * spins[:,1:])) \
           - J * (np.sum(spins[0,:] * spins[-1,:]) + np.sum(spins[:,0] * spins[:,-1]))

# Define the Metropolis algorithm
def metropolis(spins, beta):
    for i in range(N):
        for j in range(N):
            # Choose a random spin to flip
            x, y = np.random.randint(0, N), np.random.randint(0, N)
            # Calculate the energy difference
            delta_E = 2 * J * spins[x, y] * \
                      (spins[(x+1)%N, y] + spins[(x-1)%N, y] + spins[x, (y+1)%N] + spins[x, (y-1)%N])
            # Flip the spin with probability according to the Boltzmann factor
            if np.random.uniform(0, 1) < np.exp(-beta * delta_E):
                spins[x, y] = -spins[x, y]
    return spins

# Run the Monte Carlo simulation
energies = []
for i in range(nsteps):
    spins = metropolis(spins, beta)
    energies.append(energy(spins))

# Calculate the average energy
E = np.mean(energies)

# Print the result
print(f"Average energy: {E:.4f}")

# Plot the energy as a function of Monte Carlo steps
plt.plot(energies)
plt.xlabel("Monte Carlo step")
plt.ylabel("Energy")
plt.show()

# Plot the spin configuration
plt.imshow(spins, cmap='binary')
plt.show()

#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Copyright © 2023 feng <feng@archlinux>
#
# Distributed under terms of the MIT license.
"""
作者：冯振华
学号：202221140015
版本：V4.1
日期：2023年05月20日
项目：Monte Carlo methods for Overcoming critical slowing down
"""
import numpy as np
import matplotlib.pyplot as plt

def initialize_spins(N):
    """
    Initializes the spin lattice with random spins of -1 or +1
    """
    return np.random.choice([-1, 1], size=(N, N))

def bond_probability(beta):
    """
    Calculates probability of forming a bond between two neighboring spins
    """
    return 1 - np.exp(-2*beta)

def update_spins(spins, beta):
    """
    a single update of the spins using the Wolff algorithm
    """
    N = spins.shape[0]
    visited = np.zeros((N, N), dtype=bool)
    seed = (np.random.randint(N), np.random.randint(N))
    cluster_spins = [seed]
    cluster_perimeter = [(seed[0], (seed[1]+1)%N), (seed[0], (seed[1]-1)%N),
                         ((seed[0]+1)%N, seed[1]), ((seed[0]-1)%N, seed[1])]
    while len(cluster_perimeter) > 0:
        i, j = cluster_perimeter.pop()
        if not visited[i,j] and np.random.random() < bond_probability(beta):
            cluster_spins.append((i, j))
            visited[i, j] = True
            cluster_perimeter += [(i, (j+1)%N), (i, (j-1)%N),
                                  ((i+1)%N, j), ((i-1)%N, j)]
    spins[tuple(zip(*cluster_spins))] *= -1
    return spins

def calculate_energy(spins,J):
    """
    Calculates the energy of the spin configuration
    """
    N = spins.shape[0]
    energy = 0
    for i in range(N):
        for j in range(N):
            spin = spins[i, j]
            neighbors = spins[(i+1)%N, j] + spins[i, (j+1)%N] + spins[(i-1)%N, j] + spins[i, (j-1)%N]
            energy += 2*J*spin * neighbors    
    return energy 

def run_simulation(N=20, J=-1,beta=1.0, num_steps=2000, num_measurements=100):
    """
    Performs a full simulation of the Wolff algorithm on the spin lattice
    """
    spins = initialize_spins(N)
    energies = []
    mags = []
    for step in range(num_steps):
        spins = update_spins(spins, beta)
        if step % (num_steps // num_measurements) == 0:
            energy = calculate_energy(spins,J)
            mag = np.sum(spins)
            energies.append(energy)
            mags.append(mag)
    return energies, mags 

energies, mags = run_simulation(N=20,J=-1,beta=1.0, num_steps=100000, num_measurements=100)
#
# Calculate the average energy
E = np.mean(energies)
# Print the result
print(f"Average energy: {E:.4f}")
# Plot the energy as a function of Monte Carlo steps
plt.plot(energies)
plt.xlabel("Monte Carlo step")
plt.ylabel("Energy")
plt.show()

# Plot the spin configuration
plt.plot(mags)
plt.xlabel("Monte Carlo step")
plt.ylabel("Spin")
plt.show()
"""
    To run the simulation, simply the `run_simulation` function with the desired parameters, such as `energies, mags = run_simulation(N=30,=1.0, num_steps=100000, num_measurements=100)`. This will simulate 100,000 Wolff updates on a 30x30 lattice at a temperature of 1. and take measurements of the energy and magnetization every 1,000 steps. The output is two arrays, `energies` and `mags`, that record the energy and magnetization at each measurement step. These arrays can be used to averages and variances to analyze the system properties.
"""