Как считать вхождения в симуляции Python?

Question

Я выполняю симуляцию TASEP в Python, где узлы решетки на решетке заданного размера могут быть либо пустыми, либо занятыми (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()

Вот график из симуляции

Answer 1

Решение, предложенное 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()

Answer 2

Хорошо, поэтому я в основном исправил ваш код, потому что, без обид, но раньше это было что-то вроде беспорядка (см. Комментарии).

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()

Это выводит