Hw10 solutions
The important differences between Gillespie's first reaction and direct methods:
In the direct method, we picked one random number (uniformly on (0,1)) to generate the time of the next reaction (with exponential distribution parametrized by a = sum(a_i) ):
# how long until the next reaction r1 = random.random() tau = (1.0/(a+eps))*math.log(1.0/r1)
Then we picked a second random number uniformly on (0,a) and identified which reaction's "bin" it fell into (where the width of each reactions's bin is proportional to its current propensity).
r2a = random.random()*a a_sum = 0 for i in a_i: if r2a < (a_i[i]+a_sum): mu = i break a_sum += a_i[i]
Where a_sum is the "right edge" of the current bin you're considering.
In the first reaction method, one picks a random number (uniform on 0,1) and uses it to generate a a time (exponentially distributed with parameter a_i) for each reaction, and selects the reaction with the smallest time to fire at that time. Tau's for all reactions are recalculated at every step, never saved.
mintau = t_max # which reaction will happen first? # caluculate each reaction's time based on a different random number for rxn in a_i.keys(): ai = a_i[rxn] ri = random.random() taui = (1.0/(ai+eps)) * math.log(1.0/ri) if taui < mintau: # "sort as you go" mu = rxn # the putative first rxn tau = taui # the putative first time mintau = tau # reset the min
This method of sorting seemed fastest to us, but maybe storing the tau's in a dictionary and sorting the values would have been faster, like this:
tau_i = {} # which reaction will happen first? # caluculate each reaction's time based on a different random number for rxn in a_i.keys(): ai = a_i[rxn] ri = random.random() taui = (1.0/(ai+eps)) * math.log(1.0/ri) tau_i[taui] = rxn tau = min(tau_i.keys()) mu = tau_i[tau]