У меня есть большая функция, которая моделирует систему реакций во времени. Я добавил leap_couner, чтобы подсчитать количество шагов во времени.
def gillespie_tau_leaping(propensity_calc, popul_num, LHS, stoch_rate, popul_num_all, tao_all, rxn_vector, delta_t, tao, epsi):
t = simulation_timer()
t.start()
while tao < tmax:
propensity = propensity_calc(LHS, popul_num, stoch_rate)
a0 = (sum(propensity))
if a0 == 0.0:
print("Propensity sum is zero end execution")
break
if popul_num.any() < 0:
print("Molecules numbers below zero end execution")
break
lam = (propensity_calc(LHS, popul_num, stoch_rate)*delta_t)
rxn_vector = np.random.poisson(lam)
if tao + delta_t > tmax:
break
tao += delta_t
leap_counter = 0
if delta_t >= 1 / a0:
for j in range(len(rxn_vector)):
state_change_lambda = np.squeeze(np.asarray(state_change_array[j])*rxn_vector[j])
popul_num = popul_num + state_change_lambda
leap_counter += 1 # Not working!
popul_num_all.append(popul_num)
tao_all.append(tao)
leap_counter += 1
else: # else execute the ssa because it's faster
next_t = np.random.exponential(1/a0)
rxn_probability = propensity / a0
num_rxn = np.arange(rxn_probability.size)
if tao + next_t > tmax:
tao = tmax
break
j = stats.rv_discrete(values=(num_rxn, rxn_probability)).rvs()
tao = tao + next_t
popul_num = popul_num + np.squeeze(np.asarray(state_change_array[j]))
popul_num_all.append(popul_num)
tao_all.append(tao)
print("Molecule numbers:\n", popul_num)
print("Time of final simulation:\n", tao)
print("leap counter:\n", leap_counter)
t.stop()
return popul_num_all.append(popul_num), tao_all.append(tao), popul_num
Это полная функция, здесь указан c бит со счетчиком:
if tao + delta_t > tmax:
break
tao += delta_t
leap_counter = 0
if delta_t >= 1 / a0:
for j in range(len(rxn_vector)):
state_change_lambda = np.squeeze(np.asarray(state_change_array[j])*rxn_vector[j])
popul_num = popul_num + state_change_lambda
leap_counter += 1 # Not working!
popul_num_all.append(popul_num)
tao_all.append(tao)
tao - текущее время, delta_t - значение для увеличения by и tmax - максимальное время моделирования. Если tao + delta_t меньше tmax, текущее время увеличивается на delta_t.
Затем я использую delta_t, чтобы определить, какой метод прыжка я использую. Если я использую метод первого скачка if delta_t >= 1/a0, то я хочу подсчитывать каждый раз приращение (на delta_t) и возвращать количество приращений leap_counter
Только у меня проблемы со счетчиком на в тот момент, когда он просто возвращает ноль, и я не уверен, как изменить его, чтобы он работал
Ура