Reputation: 395
I've got a big function that simulates a systems of reactions over time. I've added in a leap_couner
to count the number of steps through time.
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
Thats the full function and heres the specific bit with the counter:
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
is the current time, delta_t
is the value to increase by and tmax
is the maximum simulation time. If tao + delta_t
is less than tmax
then the current time is incremented by delta_t
.
I then use delta_t
to determine which leap method I use. If I use the first leap method if delta_t >= 1/a0
then I want to count each time increment (by delta_t
) and return the number of increments leap_counter
Only I'm having trouble with the counter at the moment it just returns zero and I'm not sure how to change it to make it work
Cheers
Upvotes: 1
Views: 327
Reputation: 7981
You need to move the line
leap_counter = 0
Outside of the while loop. Currently it's getting set to 0
each iteration
Upvotes: 1