3
\$\begingroup\$

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()
```
\$\endgroup\$
0

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.