# Monte Carlo Simulation of a Markov Process

A 1D binary string has the following dynamics $$1100 \to 1010,\quad 1101\to 1011$$ each with rate $$\p<1\$$. That is if we see $$\110\$$ in the string, then we change to $$\101\$$ with probability $$\p\$$. I am interested in the count of $$\00,01\$$ and $$\11,\$$ and $$\1100\$$ in the binary string. For this, I have written the following code which I want to optimize (time).

This model is known as the Restricted Asymmetric Simple Exclusion Process (or RASEP), a modified version of TASEP.

import matplotlib.pyplot as plt
import numpy as np
from random import random

sweeps = 1000
n = 1000
p = 0.5
rho = 0.5
ensemble_size = 100

"""preparing 1D lattice with n number of 1s and L-n number of 0s.
Convinient to work with rho = n/L."""
def initialize(rho = 0.25, L = 100):
config = np.array([1]*int(rho*L)+[0]*(L-int(rho*L)))
np.random.shuffle(config)
return config

#N0 = number of 00, N1 = number of 01 and 10, N3 = number of 11, and N12 = number of 1100
def count(state):
n = len(state)
n0 = n1 = n3 = n12 = 0
for j in range(n):
neighbour_sum = int(f"{state[j%n]}{state[(j+1)%n]}{state[(j+2)%n]}{state[(j+3)%n]}",2)
if neighbour_sum % 4 == 0:
n0 += 1
if neighbour_sum == 12:
n12 += 1
elif neighbour_sum % 4 == 3:
n3 += 1
n1 = n - n0 - n3
return n0, n1, n3, n12

"""Dynamics"""
def dynamics(state, sweeps, p):
n = len(state)
n0, n1, n3, n12 = count(state)
N0, N1, N3, N12  = np.zeros(sweeps), np.zeros(sweeps), np.zeros(sweeps), np.zeros(sweeps)
# random state indices
J = np.random.randint(0, n, size=(sweeps, n))
#loop
for t in range(sweeps):
for tt in range(n):
# random indices
j = J[t, tt]
neighbour_sum = int(f"{state[j%n]}{state[(j+1)%n]}{state[(j+2)%n]}{state[(j+3)%n]}",2)
if neighbour_sum == 12 and random() < p:
state[(j+1)%n], state[(j+2)%n] = 0, 1
n0 -= 1
n1 += 2
n3 -= 1
n12 -= 1
elif neighbour_sum == 13 and random() < p:
state[(j+1)%n], state[(j+2)%n] = 0, 1
if state[(j+4)%n] == 0 and state[(j+5)%n] == 0:
n12 += 1
N0[t], N1[t], N3[t], N12[t] = n0, n1, n3, n12
return N0/n, N1/n, N3/n, N12/n

def ensembleAvg(ensemble_size):
ensemble = np.zeros((4, ensemble_size, sweeps))
for _ in range(ensemble_size):
state = initialize(rho, n)
ensemble[0,_], ensemble[1,_], ensemble[2,_], ensemble[3,_] = dynamics(state, sweeps, p)
return np.average(ensemble,axis = 1)

avg_ensemble = ensembleAvg(ensemble_size)
N0, N1, N3, N12 = avg_ensemble[0], avg_ensemble[1], avg_ensemble[2], avg_ensemble[3]

plt.title(f"size:{n}", color='b');

plt.plot(np.arange(sweeps), N0, color='blue', label = "00")
plt.plot(np.arange(sweeps), N1, color='darkgreen', label = "01 & 01")
plt.plot(np.arange(sweeps), N3, color='purple', label = "11")
plt.plot(np.arange(sweeps), N12, color='red', label = "1100")

plt.legend()
plt.xlabel("Sweeps", fontsize=20);
plt.ylabel("Number density", fontsize=20);
plt.axis('tight');

plt.savefig(f"{n}.pdf")
plt.show()
$$$$
`