# 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]

Here's the code implementing the whole first reaction method.
http://web.mit.edu/~soniat/Public/firstrxn.py