Update (July 9 2025)

 

Update 7/9/25: By The End Of Summer 2025

Main Updates: Looking In The Near Future

(1)Detailed Ising model and spin glass system replications

(2) Noisy Coin

(3)Original Research, New Simulation technique on spin glass systems with (2)

I want to get to (3) before my birthday on the 13th although that's unlikely :)

More Updates:

I've got a social media account dedicated to this page: https://bsky.app/profile/gaurav-kamath.bsky.social

I've abandoned by personal codex approach due to it's lack of replication and small scale experimentation 

needed to approach new ideas. It has however helped me generate new ideas that I wish to build in the

near future.

I also apologize for any indentations problems that may arise in my released code as copying it from

 Jupyter Notebook into blogspot was done a little poorly.

Replicating Sandpile Research: Nature Fig 1a) From Greedy Control Of Cascading Failures In Interdependent Networks

Mu Controlled Sandpile Model 

#I switched to python from C to make the sandpile model to make my life easier

#replicating Nature fig 1a: https://www.nature.com/articles/s41598-021-82843-8/figures/1

import numpy as np

import matplotlib.pyplot as plt

import random

def make_area(x, y):

    return [[0]*x for _ in range(y)]

def find_at_capacity_sites(height, XMAX, YMAX):

    return [(x,y) for x in range(XMAX) for y in range(YMAX) if height[y][x] >= 3] 

def find_not_at_capacity_sites(height, XMAX, YMAX):

    return [(x,y) for x in range(XMAX) for y in range(YMAX) if height[y][x] < 3] 

def flip(mu):

    return random.random() < mu

def add_grain(area, mu, XMAX, YMAX):

    at_capacity = find_at_capacity_sites(area, XMAX, YMAX)

    not_at_capacity = find_not_at_capacity_sites(area, XMAX, YMAX)

    if not at_capacity:

        rx, ry = random.choice(not_at_capacity)

    else:

        if flip(mu):

            rx, ry = random.choice(at_capacity)

        else:

            rx, ry = random.choice(not_at_capacity)

    area[ry][rx] += 1

    return (rx, ry)

def topple(area, XMAX, YMAX): 

    avalanche_size = 0

    unstable_sites = []

    # Find initial unstable sites

    for y in range(YMAX):

        for x in range(XMAX):

            if area[y][x] >= 4:

                unstable_sites.append((x, y))

    

    while unstable_sites:

        x, y = unstable_sites.pop()

        if area[y][x] < 4:

            continue

        avalanche_size += 1

        area[y][x] -= 4

        # Update neighbors

        for dx, dy in [(1,0), (-1,0), (0,1), (0,-1)]:

            nx, ny = x + dx, y + dy

            if 0 <= nx < XMAX and 0 <= ny < YMAX:

                area[ny][nx] += 1

                if area[ny][nx] >= 4:  # Check if neighbor becomes unstable

                    unstable_sites.append((nx, ny))

    return avalanche_size

def run_sandpile(mu, grains_to_add, XMAX, YMAX):

    area = make_area(XMAX, YMAX)

    avalanche_sizes = []

    for _ in range(grains_to_add):

        add_grain(area, mu, XMAX, YMAX)

        aval_size = topple(area, XMAX, YMAX)

        if aval_size > 0:

            avalanche_sizes.append(aval_size)

    return avalanche_sizes

def plot_loglog_cascade_size(avalanche_data, mu_values):

    plt.figure(figsize=(8,6))

    

    # Use markers/colors from paper

    markers = {0.20: 'o', 0.37: 's', 0.90: '^'}  # circle, square, triangle

    colors = {0.20: 'blue', 0.37: 'green', 0.90: 'red'}

    

    for mu, avalanches in zip(mu_values, avalanche_data):

        if len(avalanches) == 0:

            continue

            

        max_s = max(avalanches)

        bins = np.logspace(0, np.log10(max_s), 30)

        hist, bin_edges = np.histogram(avalanches, bins=bins, density=True)

        bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

        

        # Plot only non-zero points

        non_zero = hist > 0

        plt.loglog(bin_centers[non_zero], hist[non_zero], 

                   marker=markers[mu], linestyle='-', lw=1.5, ms=6,

                   color=colors[mu], label=f'μ = {mu}')

    

    plt.xlabel("Cascade size, s", fontsize=14)

    plt.ylabel("Probability, P(s)", fontsize=14)

    plt.title("Avalanche Size Distribution", fontsize=16)

    plt.legend(fontsize=12)

    plt.grid(True, which="both", ls="--", alpha=0.7)

    plt.xlim(1, 1e5)  # Match paper's axis limits

    plt.ylim(1e-6, 1)

    plt.tight_layout()

    plt.show()

if __name__ == '__main__':

    XMAX, YMAX = 100, 100

    grains_to_add = 50000

    mu_values = [0.20, 0.37, 0.90]

    avalanche_data = []

    for mu in mu_values:

        print(f"Running sandpile for μ = {mu}...")

        avalanches = run_sandpile(mu, grains_to_add, XMAX, YMAX)

        avalanche_data.append(avalanches)

    

    plot_loglog_cascade_size(avalanche_data, mu_values)

Frenkel Kontorova Model On A 2-D Lattice

 # Frenkel-Kontorova Model

                                   

#Please keep in mind to run the model multiple times. Sometimes the restoration force or time steps

#aren't enough to restore the lattice back into the lattice back.

import random

import matplotlib.pyplot as plt

XMAX = 15

YMAX = 15

delta = 0.01  

dt = 0.1

tmax = 1500

def actualx(x): return x % XMAX

def actualy(y): return y % YMAX

def metric(x, x2): return min(abs(x - x2), XMAX - abs(x - x2))

class Particle:

    def __init__(self, x, y, vx, vy, q, m):

        self.x = x

        self.y = y

        self.vx = vx

        self.vy = vy

        self.q = q

        self.m = m

    def potential(self, x, y):

        dx = metric(x, self.x)

        dy = metric(y, self.y)

        r = (dx**2 + dy**2)**0.5

        return self.q / r if r != 0 else 0

class Defect:

    def __init__(self, x, y, strength):

        self.x = x

        self.y = y

        self.strength = strength

    def force_on(self, px, py):

        dx = metric(px, self.x)

        dy = metric(py, self.y)

        r2 = dx**2 + dy**2

        if r2 == 0: return 0, 0

        return (-0.5*self.strength * dx / r2, -0.5*self.strength * dy / r2)

def chargelistgenerator(m, q, xmax, ymax, vxmax, vymax, numberofparticles):

    return [Particle(random.uniform(0, xmax), random.uniform(0, ymax),

                     random.uniform(-vxmax, vxmax), random.uniform(-vymax, vymax), q, m)

            for _ in range(numberofparticles)]

def generate_defects(num_defects, strength):

    return [Defect(random.uniform(0, XMAX), random.uniform(0, YMAX), strength)

            for _ in range(num_defects)]

def gridPotentialAtPoint(x, y, chargelist):

    return sum(charge.potential(x, y) for charge in chargelist)

def gridForceAtPointX(x, y, q, chargelist, Edrive):

    return q*Edrive - (q * (gridPotentialAtPoint(actualx(x + delta), actualy(y), chargelist) -

                            gridPotentialAtPoint(actualx(x - delta), actualy(y), chargelist)) / (2 * delta))

def gridForceAtPointY(x, y, q, chargelist):

    return - (q * (gridPotentialAtPoint(actualx(x), actualy(y + delta), chargelist) -

                   gridPotentialAtPoint(actualx(x), actualy(y - delta), chargelist)) / (2 * delta))

def particleVelocityUpdate(particle, chargelist, defects, Edrive):

    fx = gridForceAtPointX(particle.x, particle.y, particle.q, chargelist, Edrive)

    fy = gridForceAtPointY(particle.x, particle.y, particle.q, chargelist)

#account for force from defects 

    for defect in defects:

        dfx, dfy = defect.force_on(particle.x, particle.y)

        fx += dfx

        fy += dfy

    particle.vx += fx / particle.m

    particle.vy += fy / particle.m

    particle.vx *= 0.999

    particle.vy *= 0.999

def particlePositionUpdate(particle, deltat):

    particle.x = actualx(particle.x + particle.vx * deltat)

    particle.y = actualy(particle.y + particle.vy * deltat)

def run_simulation(Edrive):

    particles = chargelistgenerator(m=1, q=0.5, xmax=XMAX, ymax=YMAX,

                                     vxmax=0.05, vymax=0.05, numberofparticles=20)

    defects = generate_defects(num_defects=11, strength=0.01)

    t = 0

    while t < tmax:

        for p in particles:

            particleVelocityUpdate(p, particles, defects, Edrive)

        for p in particles:

            particlePositionUpdate(p, dt)

        t += dt

    return particles, defects

particles_before, defects = run_simulation(Edrive=0)

particles_after, _ = run_simulation(Edrive=450)

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

for i, (particles, title) in enumerate(zip([particles_before, particles_after],

                                           ['Before Edrive', 'After Edrive'])):

    x_vals = [p.x for p in particles]

    y_vals = [p.y for p in particles]

    dx_vals = [d.x for d in defects]

    dy_vals = [d.y for d in defects]

    axs[i].scatter(x_vals, y_vals, c='blue', label='Particles')

    axs[i].scatter(dx_vals, dy_vals, c='red', marker='x', s=100, label='Defects')

    axs[i].set_xlim(0, XMAX)

    axs[i].set_ylim(0, YMAX)

    axs[i].set_title(title)

    axs[i].set_aspect('equal')

    axs[i].grid(True)

plt.suptitle("Frenkel-Kontorova: Before and After Applying Driving Force", fontsize=16)

plt.show()

Lorenz Driven Operationally Deterministic Ising Model

#Lorenz Driven Operationally Deterministic Ising Model

#Please cite me for credit if you are planning to use my code

#Operation:

# Only the initial state is randomized. all else (core functions) operates deterministically. 

# Here we deviate from the traditional Ising Model which are probabilistic, stochastic in operation: 

#X(t),Y(t),Z(t) are coupled through a system of continuous non-linear differential equations.

# raw Y(t) -> used to select sites.

# normalized X(t) -> used to define spin flip conditions. 

# raw Z(t) + forcing -> equation where we add our forcing component (forcing = 0 if you want pattern formation to occur)

# The free parameters of the Lorenz equation(r,s,b) are material specific and temperature independent properties of the material 

# which I aptly call "spin-sensivity parameters" of a material on which the experiment is being conducted. 

# Model Motivation & Model Outcome:

#The goal here is to bridge the gap between theory and practice of solid state physics by abandoning random.random() 

#for the internal operations of the Ising model. 

#The Ising Model is an ontologically deterministic spin system, but we are limited in

#knowing the initial conditions due to finite precision from our measurement instruments. Different potential patterns that emerge vary. 

#depending on perturbing our inputs which we simulate due to said limitations in the knowledge of the initial condition of

 #the system.These are extrinsic limitations. In addition the model has built structural limitations that it acknowledges:

 #computers have finite resolution and time steps despite the Lorenz model being innately continuous rather than discrete

 #.So our model reflects this intrinsic limitation and should provide fruitful study for how extrinsic and intrinsic

 #limitations interact with one another in limiting predictability of the exact patterns formed in experiment.

#the coupled equations in real life. The coupled nature of the equations in addition to this intrinsic limitation with delta #time and fidelity could also give the model more utility to the modeler over simplistic substitution of random.random with #one dimensional discrete chaotic map where time steps are discrete  f: x(t) -> x(t+1)  whereas my model can do outputs for t+0.001.

# The consequence being that different pattern formation appears "stochastic" due to limitations in knowledge of the initial conditions.

#At the same time we have a "know how" that some kind of pattern of spin 'up' and 'down' will eventually emerge due to #dropping temperature. 

#when you add a forcing onto Z you are altering the geometry of the attractor and thus altering the effect of the spin-#sensitivity parameters on the spin of the material on it. You consider it akin to chemical or physical alteration of a #material. This means we may be able to study the Ising with time dependent forcing with reversible alterations of the  material 

#Potential Numerical Advantages:

# pre-packaged inputs would offer and speed and predictability in computation over random.random() operations.

import numpy as np

import matplotlib.pyplot as plt

from decimal import Decimal, getcontext, ROUND_HALF_UP

def lorenz_step(state, dt=0.01, s=15, r=23, b=2):

    x, y, z = state

    x_dot = s * (y - x)

    y_dot = r * x - y - x * z

    z_dot = x * y - b * z

    new_state = np.array([x + x_dot * dt, y + y_dot * dt, z + z_dot * dt])

    return new_state

class LorenzRNG:

    def __init__(self, initial_state=(0.0, 1.0, 1.05), dt=0.01, buffer_size=1000, skip_transient=500):

        self.state = np.array(initial_state)

        self.dt = dt

        self.x_buffer = []

        self.transient_counter = 0

        self.skip_transient = skip_transient

        self.buffer_size = buffer_size

    def next(self):

        self.state = lorenz_step(self.state, self.dt)

        x, y, _ = self.state

        if self.transient_counter < self.skip_transient:

            self.transient_counter += 1

            return 0.5, y  

        self.x_buffer.append(x)

        if len(self.x_buffer) > self.buffer_size:

            self.x_buffer.pop(0)

            

        x_min = min(self.x_buffer)

        x_max = max(self.x_buffer)

        if x_max != x_min:

            x_norm = (x - x_min) / (x_max - x_min)

        else:

            x_norm = 0.5 

        x_norm = np.clip(x_norm, 0, 1)

        return x_norm, y

XMAX, YMAX = 50, 50

N = XMAX * YMAX

def delta_energy(x, y, grid):

    return grid[x, y] * (

        grid[(x + 1) % XMAX, y] + grid[(x - 1) % XMAX, y] +

        grid[x, (y + 1) % YMAX] + grid[x, (y - 1) % YMAX]

    )

def switch_spin(x, y, grid):

    grid[x, y] *= -1

def grid_with_finder(rows, columns, fidelity, target_num):

    getcontext().prec = 28

    getcontext().rounding = ROUND_HALF_UP

    fidelity_decimal = Decimal(str(fidelity))

    target_decimal = Decimal(str(target_num))

    grid_max = rows * columns * fidelity_decimal

    if grid_max != 0:

        adjusted_target = (target_decimal % grid_max) if target_decimal >= 0 else (grid_max - (abs(target_decimal) % grid_max))

    else:

        adjusted_target = Decimal('0')

    grid = []

    current_value = Decimal('0')

    

    for y in range(rows):

        row = []

        for x in range(columns):

            wrapped_value = current_value % grid_max if grid_max != 0 else current_value

            quantized_value = wrapped_value.quantize(fidelity_decimal)

            row.append(quantized_value)

            current_value += fidelity_decimal

        

        if y % 2 == 1:

            row.reverse()

        grid.insert(0, row)

    closest = None

    min_diff = float('inf')

    closest_x, closest_y = -1, -1

    for y in range(rows):

        for x in range(columns):

            diff = abs(grid[y][x] - adjusted_target)

            if diff < min_diff:

                min_diff = diff

                closest = grid[y][x]

                closest_x, closest_y = x, y

                

    return grid, closest, (closest_x, rows - 1 - closest_y), adjusted_target

def metropolis_step(grid, temperature, lorenz_rng, scaling_factor=0.153):

    x_norm, raw_y = lorenz_rng.next()

    _, _, (x, y), _ = grid_with_finder(XMAX, YMAX, scaling_factor, raw_y)

    dE = delta_energy(x, y, grid)

    if dE <= 0:

        switch_spin(x, y, grid)

    else:

        if x_norm < np.exp(-2 * dE / temperature):

            switch_spin(x, y, grid)

def temperature_schedule(step, total_steps, T_start, T_end):

    return T_end + (T_start - T_end) * np.exp(-5 * step / total_steps)

def run_simulation_with_snapshots(T_start, T_end, snapshot_steps, lorenz_initial_state, scaling_factor=0.153, initial_grid=None):

    lorenz_rng = LorenzRNG(initial_state=lorenz_initial_state, dt=0.01)

    if initial_grid is None:

        grid = np.random.choice([-1, 1], size=(XMAX, YMAX))

    else:

        grid = np.copy(initial_grid)

    snapshots = {}

    total_steps = max(snapshot_steps)

    for step in range(1, total_steps + 1):

        temperature = temperature_schedule(step, total_steps, T_start, T_end)

        for _ in range(XMAX * YMAX):

            metropolis_step(grid, temperature, lorenz_rng, scaling_factor=scaling_factor)

        if step in snapshot_steps:

            snapshots[step] = np.copy(grid)

    return snapshots

np.random.seed(42)

initial_grid = np.random.choice([-1, 1], size=(XMAX, YMAX))

fixed_initial_condition = (0.15, 2.32, 6.04)

scaling_factors = [0.0001]  

snapshot_steps = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

T_start = 100

T_end = 0.1

fig, axes = plt.subplots(len(scaling_factors), len(snapshot_steps), figsize=(16, 3 * len(scaling_factors)))

for i, scale in enumerate(scaling_factors):

    snapshots = run_simulation_with_snapshots(

        T_start, T_end, snapshot_steps, fixed_initial_condition,

        scaling_factor=scale,

        initial_grid=initial_grid

    )

    for j, step in enumerate(snapshot_steps):

        ax = axes[i, j] if len(scaling_factors) > 1 else axes[j]

        ax.imshow(snapshots[step], cmap='gray', interpolation='nearest')

        ax.set_title(f"Step {step}")

        ax.axis('off')

plt.suptitle(f"Lorenz Driven Operationally Deterministic Ising Model. Scale: {scale}", fontsize=18)

plt.tight_layout(rect=[0, 0, 1, 0.96])

plt.show()