如何计算Python模拟中出现的次数？

Question

我正在Python中运行TASEP模拟，其中给定大小的晶格上的晶格站点可以是空的，也可以是被占用的（0或1）。你知道吗

这个模拟给了我一个给定模拟时间的晶格配置图（不管一个状态是否被占用），但没有给出被占用状态的数量。你知道吗

我无法让Python计算占用状态的数量，因为该图来自模拟，而不是列表。你知道吗

TASEP代码：

import random, pylab, math
random.seed()
L=100 # Number of lattice sites
alpha=.2 # Rate of entry
beta=.4 # Rate of exit

Ntime=200000 # Simulation steps
state=[0 for k in range(L+1)]
for iter in range(Ntime):
   k=random.randint(0,L)
   if k==0:
      if random.random() < alpha: state[1]=1
   elif k==L:
         if random.random() < beta: state[L]=0
   elif state[k]==1 and state[k+1]==0: 
      state[k+1]=1
      state[k]=0
   if iter%2000 == 0: 
      yaxis=[]
      for i in range(L):
          if state[i+1]==1: yaxis.append(i)
      xaxis=[iter for k in range(len(yaxis))]
      pylab.plot(xaxis,yaxis,'r.')
pylab.xlabel('Number of steps')
pylab.ylabel('System configuration')
pylab.show()

Here is a plot from the simulation

Answer 1

好的，所以我基本上修复了你的代码，因为，没有冒犯，但它之前有点混乱（见注释）。你知道吗

import random
import matplotlib.pyplot as plt

fig, axs = plt.subplots(nrows=2, sharex=True)

L = 100  # Number of lattice sites
alpha = 0.2  # Rate of entry
beta = 0.4  # Rate of exit

n_time = 200000  # Simulation steps
record_each = 2000
state = [0]*(L + 1)
record = []  # store a record of the total number of states occupied

for itr in range(n_time):

    rand_int = random.randint(0, L)

    if rand_int == 0:
        if random.random() < alpha:
            state[1] = 1
    elif rand_int == L:
        if random.random() < beta:
            state[L] = 0
    elif state[rand_int] == 1 and state[rand_int + 1] == 0:
        state[rand_int + 1] = 1
        state[rand_int] = 0

    if itr % record_each == 0:
        yaxis = [i for i in range(L) if state[i + 1] == 1]
        axs[1].plot([itr]*len(yaxis), yaxis, 'r.')
        record.append(sum(state))  # add the total number of states to the record

axs[0].plot(range(0, n_time, record_each), record)  # plot the record
axs[1].set_xlabel('Number of steps')
axs[1].set_ylabel('System configuration')
axs[0].set_ylabel('Number of states occupied')
plt.show()

这个输出

Answer 2

FHTMitchell提出的解决方案是正确的，但效率低下。sum操作要求在每次迭代中执行O（L）个工作，使整个程序成为O（L*n\u时间）。你知道吗

请注意：

• 开始时occupied_state_count为0
• occupied_state_count应该只在零状态转换为一状态时递增
• 只有当一个状态被切换到零时，occupied_state_count才应该递减
• 如果您已经处于目标状态，则不需要更改，短路可以避免对random()的不必要调用
• 最后，当两个状态在相反的方向上切换时（你最后的elif），没有必要改变occupied_state_count。你知道吗

应用以上所有方法可以得到以下O（nаu time）实现，这要快得多：

import random
import matplotlib.pyplot as plt

fig, axs = plt.subplots(nrows=2, sharex=True)

L = 100  # Number of lattice sites
alpha = 0.2  # Rate of entry
beta = 0.4  # Rate of exit

n_time = 200000  # Simulation steps
record_each = 2000
state = [0]*(L + 1)
occupied_state_count = 0
record = []  # store a record of the total number of states occupied

for itr in range(n_time):

    rand_int = random.randint(0, L)

    if rand_int == 0:
        if state[1] == 0 and random.random() < alpha:
            state[1] = 1
            occupied_state_count += 1
    elif rand_int == L:
        if state[L] == 1 and random.random() < beta:
            state[L] = 0
            occupied_state_count -= 1
    elif state[rand_int] == 1 and state[rand_int + 1] == 0:
        state[rand_int + 1] = 1
        state[rand_int] = 0

    if itr % record_each == 0:
        yaxis = [i for i in range(L) if state[i + 1] == 1]
        axs[1].plot([itr]*len(yaxis), yaxis, 'r.')
        record.append(occupied_state_count)  # add the total number of states to the record

axs[0].plot(range(0, n_time, record_each), record)  # plot the record
axs[1].set_xlabel('Number of steps')
axs[1].set_ylabel('System configuration')
axs[0].set_ylabel('Number of states occupied')
plt.show()