How can I get a python counter to work in a function with multiple if statements?

Question

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

Answer 1

You need to move the line

leap_counter = 0

Outside of the while loop. Currently it's getting set to 0 each iteration