如何计算 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

好的，所以我基本上修复了您的代码，因为没有冒犯，但它之前有点乱（参见 cmets）。

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_time)。

请注意：

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

应用上述所有方法会产生以下 O(n_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()