Simulated Annealing: Ising Model Project

Copycopy code to clipboard
class Ising2D:
    """2D Ising model with periodic boundary conditions."""
    
    def __init__(self, L, J=1.0, random=True):
        self.L = L
        self.J = J
        self.N = L * L
        
        if random:
            self.spins = np.random.choice([-1, 1], size=(L, L)).astype(np.int8)
        else:
            self.spins = np.ones((L, L), dtype=np.int8)
    
    def get_neighbors(self, i, j):
        """Get 4 neighbors with PBC."""
        return [
            self.spins[i, (j+1) % self.L],  # Right
            self.spins[i, (j-1) % self.L],  # Left
            self.spins[(i+1) % self.L, j],  # Down
            self.spins[(i-1) % self.L, j]   # Up
        ]
    
    def delta_energy(self, i, j):
        """Fast energy change: O(1)."""
        neighbors = self.get_neighbors(i, j)
        sum_neighbors = sum(neighbors)
        return 2 * self.J * self.spins[i, j] * sum_neighbors
    
    def flip_spin(self, i, j):
        """Flip spin at (i, j)."""
        self.spins[i, j] *= -1
    
    def get_energy(self):
        """Compute total energy."""
        E = 0.0
        for i in range(self.L):
            for j in range(self.L):
                neighbors = self.get_neighbors(i, j)
                E -= self.J * self.spins[i, j] * sum(neighbors)
        return E / 2
    
    def get_energy_per_spin(self):
        return self.get_energy() / self.N
    
    def get_magnetization_per_spin(self):
        return abs(np.sum(self.spins)) / self.N


def metropolis_sweep_2d(ising, T, seed=None):
    """One Metropolis sweep (2D)."""
    if seed is not None:
        np.random.seed(seed)
    
    accepted = 0
    for _ in range(ising.N):
        i = np.random.randint(0, ising.L)
        j = np.random.randint(0, ising.L)
        delta_E = ising.delta_energy(i, j)
        if delta_E < 0 or np.random.random() < np.exp(-delta_E / T):
            ising.flip_spin(i, j)
            accepted += 1
    return accepted / ising.N


def run_fixed_T_2d(ising, T, n_sweeps, n_equil=100, seed=None):
    """Run Metropolis at fixed temperature T (2D)."""
    if seed is not None:
        np.random.seed(seed)
    
    # Equilibration
    for _ in range(n_equil):
        metropolis_sweep_2d(ising, T)
    
    # Production
    energies = []
    magnetizations = []
    for _ in range(n_sweeps):
        metropolis_sweep_2d(ising, T)
        energies.append(ising.get_energy_per_spin())
        magnetizations.append(ising.get_magnetization_per_spin())
    
    return {
        'energies': np.array(energies),
        'magnetizations': np.array(magnetizations)
    }

Copycopy code to clipboard
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from tqdm import tqdm

class Ising3D:
    """3D Ising model with periodic boundary conditions."""
    
    def __init__(self, L, J=1.0, random=True):
        """Initialize 3D Ising lattice.
        
        Args:
            L: Linear size (L×L×L lattice)
            J: Coupling constant
            random: If True, random spins; if False, all +1 (ground state)
        """
        self.L = L
        self.J = J
        self.N = L * L * L  # Total number of spins
        
        # Initialize spins: 3D array of ±1
        if random:
            self.spins = np.random.choice([-1, 1], size=(L, L, L)).astype(np.int8)
        else:
            self.spins = np.ones((L, L, L), dtype=np.int8)
    
    def get_neighbors(self, i, j, k):
        """Get 6 neighbors with periodic boundary conditions (PBC).
        
        Uses modulo operator: (i+1) % L automatically wraps around.
        """
        neighbors = []
        # 6 neighbors in 3D: ±x, ±y, ±z
        for di, dj, dk in [(-1,0,0), (1,0,0), (0,-1,0), (0,1,0), (0,0,-1), (0,0,1)]:
            ni = (i + di) % self.L  # PBC: wraps around
            nj = (j + dj) % self.L
            nk = (k + dk) % self.L
            neighbors.append((ni, nj, nk))
        return neighbors
    
    def delta_energy(self, i, j, k):
        """Fast energy change calculation: O(1) instead of O(N).
        
        Only neighbors contribute to energy change when flipping one spin.
        Formula: ΔE = 2J s_i Σ_j s_j (sum over neighbors)
        """
        neighbors = self.get_neighbors(i, j, k)
        sum_neighbors = sum(self.spins[ni, nj, nk] for ni, nj, nk in neighbors)
        return 2 * self.J * self.spins[i, j, k] * sum_neighbors
    
    def flip_spin(self, i, j, k):
        """Flip spin at position (i, j, k)."""
        self.spins[i, j, k] *= -1
    
    def get_energy(self):
        """Compute total energy (slow, for verification)."""
        E = 0.0
        for i in range(self.L):
            for j in range(self.L):
                for k in range(self.L):
                    neighbors = self.get_neighbors(i, j, k)
                    for ni, nj, nk in neighbors:
                        E -= self.J * self.spins[i, j, k] * self.spins[ni, nj, nk]
        return E / 2  # Divide by 2 to avoid double counting
    
    def get_energy_per_spin(self):
        """Get energy per spin."""
        return self.get_energy() / self.N
    
    def get_magnetization_per_spin(self):
        """Get magnetization per spin."""
        return np.sum(self.spins) / self.N
    
    def visualize_slice(self, axis=2, slice_idx=None, ax=None):
        """Visualize a 2D slice of the 3D configuration."""
        if slice_idx is None:
            slice_idx = self.L // 2
        
        if ax is None:
            fig, ax = plt.subplots(1, 1, figsize=(8, 8))
        
        # Extract slice
        if axis == 0:
            slice_data = self.spins[slice_idx, :, :]
        elif axis == 1:
            slice_data = self.spins[:, slice_idx, :]
        else:
            slice_data = self.spins[:, :, slice_idx]
        
        cmap = ListedColormap(['blue', 'red'])
        im = ax.imshow(slice_data, cmap=cmap, vmin=-1, vmax=1)
        ax.set_title(f'3D Ising Slice (axis={axis}, slice={slice_idx})', fontsize=11)
        return ax


def metropolis_sweep_3d(ising, T, seed=None):
    """Perform one Metropolis sweep at temperature T (3D version).
    
    One sweep = N proposed spin flips (one per spin, in random order).
    """
    if seed is not None:
        np.random.seed(seed)
    
    accepted = 0
    
    # Propose N flips (one sweep)
    for _ in range(ising.N):
        # Random spin
        i = np.random.randint(0, ising.L)
        j = np.random.randint(0, ising.L)
        k = np.random.randint(0, ising.L)
        
        # Fast energy change calculation
        delta_E = ising.delta_energy(i, j, k)
        
        # Metropolis acceptance
        if delta_E < 0 or np.random.random() < np.exp(-delta_E / T):
            ising.flip_spin(i, j, k)
            accepted += 1
    
    return accepted / ising.N

def run_fixed_T_3d(ising, T, n_sweeps, n_equil=100, seed=None):
    """Run Metropolis at fixed temperature T (3D)."""
    if seed is not None:
        np.random.seed(seed)
    
    # Equilibration
    for _ in range(n_equil):
        metropolis_sweep_3d(ising, T)
    
    # Production
    energies = []
    magnetizations = []
    for _ in range(n_sweeps):
        metropolis_sweep_3d(ising, T)
        energies.append(ising.get_energy_per_spin())
        magnetizations.append(abs(ising.get_magnetization_per_spin()))
    
    return {
        'energies': np.array(energies),
        'magnetizations': np.array(magnetizations)
    }

Copycopy code to clipboard
class GeometricSchedule:
    """Geometric cooling schedule: T_{k+1} = alpha * T_k"""
    
    def __init__(self, T0, T_min, alpha=0.95):
        self.T0 = T0
        self.T_min = T_min
        self.alpha = alpha
        self.T = T0
        self.k = 0
    
    def get_temperature(self):
        return self.T
    
    def cool(self):
        self.T = max(self.T * self.alpha, self.T_min)
        self.k += 1
    
    def is_done(self):
        return self.T <= self.T_min


def simulated_annealing_3d(ising, schedule, sweeps_per_T=10, seed=None):
    """Run Simulated Annealing on 3D Ising model.
    
    Args:
        ising: Ising3D instance
        schedule: Cooling schedule (GeometricSchedule)
        sweeps_per_T: Number of sweeps at each temperature
        seed: Random seed
    
    Returns:
        Dictionary with history (energies, temperatures, acceptance_rates)
    """
    if seed is not None:
        np.random.seed(seed)
    
    energies = []
    temperatures = []
    acceptance_rates = []
    
    while not schedule.is_done():
        T = schedule.get_temperature()
        
        # Run sweeps at this temperature
        for _ in range(sweeps_per_T):
            acc_rate = metropolis_sweep_3d(ising, T)
            energies.append(ising.get_energy_per_spin())
            temperatures.append(T)
            acceptance_rates.append(acc_rate)
        
        schedule.cool()
    
    return {
        'energies': np.array(energies),
        'temperatures': np.array(temperatures),
        'acceptance_rates': np.array(acceptance_rates)
    }