Task_2.py

import numpy as np
import matplotlib.pyplot as plt 
import mcmc
import Ising_class
import itertools

def task_2_1(Ns=np.array([10,25,100]), Ts=np.array([0.5,2.5,5.0])):
    print("Visualizing configurations at different temperatures...")

    fig, axes = plt.subplots(len(Ns), len(Ts), figsize=(9, 9))
    fig.suptitle("Spin Configurations for Different N and T", fontsize=14)

    for i, N in enumerate(Ns):
        for j, T in enumerate(Ts):
            B = 1 / T
            
            Ising = Ising_class.Ising(N, seed=33)

            s, E_seq, acc_burn, acc_meas = mcmc.mcmc(Ising,
                    samples = 1, wait=20, burn_in=50000, beta=B,
                    seed = None, debug_delta_cost = False) # set to True to enable the check

            ax = axes[i, j]
            ax.imshow(Ising.s, cmap='binary')
            ax.set_title(f"N={N}, T={Ts[j]:.2f}")
            ax.axis('off')
            

    plt.tight_layout()
    plt.show()


def check_stability(energies, tolerance=0.05):

    if len(energies) < 2:
        return False
        
    k = len(energies) // 2
    m1 = np.mean(energies[:k])
    m2 = np.mean(energies[k:])
    
    denom = (m1 + m2) / 2.0
    if abs(denom) < 1e-9: 
        denom = 1.0
        
    rel_error = np.abs((m2 - m1) / denom)
    
    return rel_error <= tolerance

def compute_burn_in_map(Ns=[10, 25, 100], Ts = np.arange(0.1, 5.1, 0.2)): 
    """
    Finds burn-in by increasing wait time until post-burn-in samples are stable.
    This is how i investigated the behavior at Tc for the report
    """
    burn_in_map = {}
    print("Computing burn-in...")

    for N, T in itertools.product(Ns, Ts):
        beta_val = 1.0 / T
        
        burn_in_attempt = 2000   
        factor = 1.5             
        cap = 4000000      #limit
        found = False
        test_samples = 20      
        wait_steps = 50          
        
        while burn_in_attempt <= cap:
            model = Ising_class.Ising(N=N, seed=0)
            
            _, energies, _ = mcmc.mcmc(model, burn_in=int(burn_in_attempt), 
                                       samples=test_samples, wait=wait_steps, beta=beta_val,seed=0 )
            
            if check_stability(energies, tolerance=0.05):
                found = True
                break
            
            burn_in_attempt = int(burn_in_attempt * factor)
        
        final_burn_in = burn_in_attempt if found else cap
        burn_in_map[(N, T)] = final_burn_in
        
        status = "Converged" if found else "Capped"
        print(f"  > (N={N}, T={T}): {final_burn_in} steps ({status})")

    return burn_in_map
    
burn_ins = {
    (10, 0.5): 4500,
    (10, 2.5): 2000,
    (10, 5.0): 2000,
    
    (25, 0.5): 115323,
    (25, 2.5): 76882,
    (25, 5.0): 6750,
    
    (100, 0.5): 4000000, 
    (100, 2.5): 583821, 
    (100, 5.0): 115300, 
    }

def burn_in_viz(Ns=np.array([10,25,100]), Ts=np.array([0.5,2.5,5.0])):
    """
    Where the function seem constant is where we can take burn_in value
    """
    fig, axes = plt.subplots(len(Ns), len(Ts), figsize=(14, 14))

    for i in range(len(Ns)):
        for j in range(len(Ts)):
                N, B = Ns[i], 1 / Ts[j]
                Ising = Ising_class.Ising(N, seed=None)

                s, E_seq, acc_burn = mcmc.mcmc(Ising,
                    samples = 1, wait=20, burn_in=500000, beta=B,
                    seed = None, debug_delta_cost = False)
                # plt.plot(E_seq, acc_burn)
                # plt.title(f"E_seq vs acc_burn at N={N}, Beta = {B}")
                ax = axes[i, j]
            
                ax.plot(E_seq, linewidth=1)

                ax.set_title(
                f"N={N}, T={1/B:.2f}\nacc_burn={acc_burn:.3f}"
                )

    plt.tight_layout()
    plt.show() 
    
# def plot_burn_in(Ts=np.arange(0.1, 5.1, 0.2)):
#     burn_in = [76882,76882,115323,115323,76882,172984,76882,115323,76882,
#              51255,115323,51255,76882,51255,34170,10125,22780,10125,6750,
#              6750,6750,10125,15187,10125,6750] #results from previous function
#     burn_in = np.array(burn_in)
#     plt.plot(Ts, burn_in)
    
def task_2_2(Ns=np.array([10,25,100]), Ts=np.array([0.5,2.5,5.0])):
    
    results = {N: {'T': [], 'acc': [], 'energy': []} for N in Ns}

    print("Running simulations with burn-in values previously found...")
    print(f"{'N':<5} {'T':<5} {'Burn-in':<10} {'Acc Rate':<10} {'Energy/Spin':<10}")

    for N in Ns:
        for T in Ts:

            burn_in_steps = burn_ins.get((N, T), 2000)
            beta_val = 1.0 / T
            model = Ising_class.Ising(N=N, )
            
            _, energies, accepted = mcmc.mcmc(model, burn_in=burn_in_steps, 
                                              samples=1000, wait=10, beta=beta_val)
            
            mean_acc_rate = np.mean(accepted)
            mean_energy = np.mean(energies)
            energy_per_spin = mean_energy / (N * N)

            results[N]['T'].append(T)
            results[N]['acc'].append(mean_acc_rate)
            results[N]['energy'].append(energy_per_spin)
            
            print(f"{N:<5} {T:<5.1f} {burn_in_steps:<10} {mean_acc_rate:<10.4f} {energy_per_spin:<10.4f}")

    #Acc rate 
    plt.figure(figsize=(8, 6))
    for N in Ns:
        plt.plot(results[N]['T'], results[N]['acc'], marker='o', label=f'N={N}')
    
    plt.title("Acceptance Rate vs Temperature")
    plt.xlabel("Temperature (T)")
    plt.ylabel("Acceptance Rate")
    plt.ylim(-0.05, 1.05)
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.show()

    #Cost plot
    plt.figure(figsize=(8, 6))
    for N in Ns:
        plt.plot(results[N]['T'], results[N]['energy'], marker='s', label=f'N={N}')
    
    plt.title("Mean Cost (Energy per Spin) vs Temperature")
    plt.xlabel("Temperature (T)")
    plt.ylabel("Energy / Spin")
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.show()

    return results
    


# task_2_1()
# compute_burn_in_map()
# burn_in_viz()
# task_2_2()



#OUTPUT burn_in_map()

#   > (N=10, T=0.1): 4500 steps (Converged)
#   > (N=10, T=0.30000000000000004): 4500 steps (Converged)
#   > (N=10, T=0.5000000000000001): 4500 steps (Converged)
#   > (N=10, T=0.7000000000000001): 4500 steps (Converged)
#   > (N=10, T=0.9000000000000001): 4500 steps (Converged)
#   > (N=10, T=1.1000000000000003): 4500 steps (Converged)
#   > (N=10, T=1.3000000000000003): 4500 steps (Converged)
#   > (N=10, T=1.5000000000000004): 4500 steps (Converged)
#   > (N=10, T=1.7000000000000004): 4500 steps (Converged)
#   > (N=10, T=1.9000000000000004): 6750 steps (Converged)
#   > (N=10, T=2.1000000000000005): 2000 steps (Converged)
#   > (N=10, T=2.3000000000000007): 4500 steps (Converged)
#   > (N=10, T=2.5000000000000004): 2000 steps (Converged)
#   > (N=10, T=2.7000000000000006): 2000 steps (Converged)
#   > (N=10, T=2.900000000000001): 4500 steps (Converged)
#   > (N=10, T=3.1000000000000005): 2000 steps (Converged)
#   > (N=10, T=3.3000000000000007): 3000 steps (Converged)
#   > (N=10, T=3.500000000000001): 2000 steps (Converged)
#   > (N=10, T=3.7000000000000006): 2000 steps (Converged)
#   > (N=10, T=3.900000000000001): 2000 steps (Converged)
#   > (N=10, T=4.1000000000000005): 2000 steps (Converged)
#   > (N=10, T=4.300000000000001): 4500 steps (Converged)
#   > (N=10, T=4.500000000000001): 15187 steps (Converged)
#   > (N=10, T=4.7): 10125 steps (Converged)
#   > (N=10, T=4.9): 4500 steps (Converged)
#   > (N=25, T=0.1): 76882 steps (Converged)
#   > (N=25, T=0.30000000000000004): 76882 steps (Converged)
#   > (N=25, T=0.5000000000000001): 115323 steps (Converged)
#   > (N=25, T=0.7000000000000001): 115323 steps (Converged)
#   > (N=25, T=0.9000000000000001): 76882 steps (Converged)
#   > (N=25, T=1.1000000000000003): 172984 steps (Converged)
#   > (N=25, T=1.3000000000000003): 76882 steps (Converged)
#   > (N=25, T=1.5000000000000004): 115323 steps (Converged)
#   > (N=25, T=1.7000000000000004): 76882 steps (Converged)
#   > (N=25, T=1.9000000000000004): 51255 steps (Converged)
#   > (N=25, T=2.1000000000000005): 115323 steps (Converged)
#   > (N=25, T=2.3000000000000007): 51255 steps (Converged)
#   > (N=25, T=2.5000000000000004): 76882 steps (Converged)
#   > (N=25, T=2.7000000000000006): 51255 steps (Converged)
#   > (N=25, T=2.900000000000001): 34170 steps (Converged)
#   > (N=25, T=3.1000000000000005): 10125 steps (Converged)
#   > (N=25, T=3.3000000000000007): 22780 steps (Converged)
#   > (N=25, T=3.500000000000001): 10125 steps (Converged)
#   > (N=25, T=3.7000000000000006): 6750 steps (Converged)
#   > (N=25, T=3.900000000000001): 6750 steps (Converged)
#   > (N=25, T=4.1000000000000005): 6750 steps (Converged)
#   > (N=25, T=4.300000000000001): 10125 steps (Converged)
#   > (N=25, T=4.500000000000001): 15187 steps (Converged)
#   > (N=25, T=4.7): 10125 steps (Converged)
#   > (N=25, T=4.9): 6750 steps (Converged)
#   > (N=100, T=0.1): 1000000 steps (Capped)
#   > (N=100, T=0.30000000000000004): 1000000 steps (Capped)
#   > (N=100, T=0.5000000000000001): 1000000 steps (Capped)
#   > (N=100, T=0.7000000000000001): 1000000 steps (Capped)
#   > (N=100, T=0.9000000000000001): 1000000 steps (Capped)
#   > (N=100, T=1.1000000000000003): 1000000 steps (Capped)
#   > (N=100, T=1.3000000000000003): 1000000 steps (Capped)
#   > (N=100, T=1.5000000000000004): 1000000 steps (Capped)
#   > (N=100, T=1.7000000000000004): 1000000 steps (Capped)
#   > (N=100, T=1.9000000000000004): 1000000 steps (Capped)
#   > (N=100, T=2.1000000000000005): 1000000 steps (Capped)
#   > (N=100, T=2.3000000000000007): 1000000 steps (Capped)
#   > (N=100, T=2.5000000000000004): 583821 steps (Converged)
#   > (N=100, T=2.7000000000000006): 389214 steps (Converged)
#   > (N=100, T=2.900000000000001): 259476 steps (Converged)
#   > (N=100, T=3.1000000000000005): 172984 steps (Converged)
#   > (N=100, T=3.3000000000000007): 259476 steps (Converged)
#   > (N=100, T=3.500000000000001): 259476 steps (Converged)
#   > (N=100, T=3.7000000000000006): 172984 steps (Converged)
#   > (N=100, T=3.900000000000001): 172984 steps (Converged)
#   > (N=100, T=4.1000000000000005): 172984 steps (Converged)
#   > (N=100, T=4.300000000000001): 172984 steps (Converged)
#   > (N=100, T=4.500000000000001): 259476 steps (Converged)
#   > (N=100, T=4.7): 115323 steps (Converged)
#   > (N=100, T=4.9): 115323 steps (Converged)