我具有以下功能，可以随时间模拟生化系统。代码背后的想法是对时间步进行采样，然后将其添加到当前时间。
我有一个全局定义的开始时间：tao = 0.0
和以下功能：
def gillespie_tau_leaping(propensity_calc, popul_num, popul_num_all, tao_all, rxn_vector, tao, epsi):
t = simulation_timer()
t.start()
while tao < tmax:
evaluate_propensity = propensity_calc(LHS, popul_num, stoch_rate)
a0 = (sum(propensity)) # a0 is a numpy.float64
if a0 == 0.0:
break
if popul_num.any() < 0: # This isnt working
print("Number of molecules {} is too small for reaction to fire".format(popul_num))
break
# code equation 26 here!
for x in range(len(popul_num)):
# equation 22:
expctd_net_change = a0*state_change_array[x]
# equation 24:
part_diff = derivative(evaluate_propensity, popul_num[x])
# find the partial derivative of propensity with respect to popul_num (number of discrete molecules)
# need to find a a way to SELECT delta_t before I can code it!
# equation 26:
t_step = epsi*a0 / sum(expctd_net_change*part_diff)
delta_t = optimize.fmin(t_step, 0.00012)
print (delta_t)
lam = (evaluate_propensity*delta_t)
rxn_vector = np.random.poisson(lam) # probability of a reaction firing in the given time period
if tao + delta_t > tmax:
break
tao = tao + delta_t
leap_counter = tao / delta_t
if tao >= 2/a0: # if tao is big enough
for j in range(len(rxn_vector)):
state_change_lambda = np.squeeze(np.asarray(state_change_array[j])*rxn_vector[j])
#new_propensity = evaluate_propensity
propensity_check = evaluate_propensity.copy()
propensity_check[0] += state_change_lambda[0]
propensity_check[1:] += state_change_lambda
for n in range(len(propensity_check)):
if propensity_check[n] - evaluate_propensity[n] >= epsi*a0:
print("The value of delta_t {} choosen is too large".format(delta_t))
break
else:
popul_num = popul_num + state_change_lambda
popul_num_all.append(popul_num)
tao_all.append(tao)
else:
t = np.random.exponential(1/a0)
rxn_probability = evaluate_propensity / a0
num_rxn = np.arange(rxn_probability.size)
if tao + t > tmax:
tao = tmax
break
j = stats.rv_discrete(values=(num_rxn, rxn_probability)).rvs()
tao = tao + t
popul_num = popul_num + np.squeeze(np.asarray(state_change_array[j]))
popul_num_all.append(popul_num)
tao_all.append(tao)
print("tao:\n", tao)
print("Molecule numbers:\n", popul_num)
print("Number of leaps taken:\n", leap_counter)
t.stop()
return popul_num_all.append(popul_num), tao_all.append(tao)
打印（gillespie_tau_leaping（倾向计算，popul_num，popul_num_all，tao_all，rxn_vector，tao，epsi））
仅出现错误：UnboundLocalError: local variable 'leap_counter' referenced before assignment尝试打印或返回模拟过程中执行的时间步数时，保存在变量jump_counter中，我不知道为什么。我已经尝试过更改缩进量，但这是行不通的。怎么回事，如何解决？
编辑：以下代码指定全局环境中的delta_t的值。我已经尝试过在上面进行更新，因此我在每次迭代时都要计算delta_t的值。但是后来我遇到了跳数计数器的问题。
def gillespie_tau_leaping(propensity_calc, popul_num, popul_num_all, tao_all, rxn_vector, tao, delta_t, epsi):
t = simulation_timer()
t.start()
while tao < tmax:
propensity = propensity_calc(LHS, popul_num, stoch_rate)
print("propensity:\n", type(propensity))
a0 = (sum(propensity))
if a0 == 0.0:
break
# if reaction cannot fire corresponding element in rxn_vector should be zero --> tau leaping method
if popul_num.any() < 0:
print("Number of molecules {} is too small for reaction to fire".format(popul_num))
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
print("tao:\n", tao)
print("Molecule numbers:\n", popul_num)
# divide tao by delta_t to calculate number of leaps
leap_counter = tao / delta_t
if tao >= 2/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
new_propensity = propensity_calc(LHS, popul_num, stoch_rate)
for n in range(len(new_propensity)):
propensity_check = propensity + state_change_lambda
if propensity_check[n] - new_propensity[n] >= epsi*a0:
print("The value of delta_t {} choosen is too large".format(delta_t))
break
else:
popul_num = popul_num + state_change_lambda
popul_num_all.append(popul_num)
tao_all.append(tao)
else:
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]))
t.stop()
return popul_num_all.append(popul_num), tao_all.append(tao), leap_counter
此代码完全按照我想要的方式工作，我无法理解的是第一个功能突然间不起作用
欢呼声
[leap_counter仅在while循环内可见（因此可访问），请尝试在循环外创建（最好使用合理的默认值）。