Hw10 solutions - Revision history

Ldavid at 18:19, 8 December 2006

Ldavid — Fri, 08 Dec 2006 18:19:46 GMT

←Older revision		Revision as of 18:19, 8 December 2006
Line 23:		Line 23:

	mintau = t_max		mintau = t_max
		+	eps = math.exp(-200)
	# which reaction will happen first?		# which reaction will happen first?
	# caluculate each reaction's time based on a different random number		# caluculate each reaction's time based on a different random number

Ldavid at 18:18, 8 December 2006

Ldavid — Fri, 08 Dec 2006 18:18:43 GMT

←Older revision		Revision as of 18:18, 8 December 2006
Line 3:		Line 3:
	The important differences between Gillespie's '''first reaction''' and '''direct''' methods:		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 <tt> a = sum(a_i) </tt>):	+	In the '''direct method,''' we consider the time and identity of the next reaction independently. First, we picked one random number (uniformly on (0,1)) to generate the ''time'' of the next reaction (with exponential distribution parametrized by <tt> a = sum(a_i) </tt>):

	# how long until the next reaction		# how long until the next reaction
Line 20:		Line 20:
	Where <tt>a_sum</tt> is the "right edge" of the current bin you're considering.		Where <tt>a_sum</tt> 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.	+	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. Note that the behavior of the '''first''' and '''direct''' methods is equivalent.

	mintau = t_max		mintau = t_max
Line 34:		Line 34:
	mintau = tau # reset the min		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~~:	+	This method of sorting seemed fastest to us, but here's another way of tackling the problem:

	tau_i = {}		tau_i = {}

Ldavid at 18:13, 8 December 2006

Ldavid — Fri, 08 Dec 2006 18:13:38 GMT

←Older revision		Revision as of 18:13, 8 December 2006
Line 9:		Line 9:
	tau = (1.0/(a+eps))*math.log(1.0/r1)		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).	+	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 equal to its current propensity).

	r2a = random.random()*a		r2a = random.random()*a

SoniaT: Hw11 solutions moved to Hw10 solutions: typo- oops!

SoniaT — Fri, 08 Dec 2006 17:49:18 GMT

Hw11 solutions moved to Hw10 solutions: typo- oops!

←Older revision

Revision as of 17:49, 8 December 2006

SoniaT at 17:47, 8 December 2006

SoniaT — Fri, 08 Dec 2006 17:47:34 GMT

←Older revision		Revision as of 17:47, 8 December 2006
Line 1:		Line 1:
		+	[http://web.mit.edu/~soniat/Public/firstrxn.py Solution in code form]
		+
	The important differences between Gillespie's '''first reaction''' and '''direct''' methods:		The important differences between Gillespie's '''first reaction''' and '''direct''' methods:

Line 42:		Line 44:
	taui = (1.0/(ai+eps)) * math.log(1.0/ri)		taui = (1.0/(ai+eps)) * math.log(1.0/ri)
	tau_i[taui] = rxn		tau_i[taui] = rxn
-
	tau = min(tau_i.keys())		tau = min(tau_i.keys())
	mu = tau_i[tau]		mu = tau_i[tau]
-
-
-	~~Here's the code implementing the whole first reaction method.~~
-	~~http://web.mit.edu/~soniat/Public/firstrxn.py~~

SoniaT at 17:45, 8 December 2006

SoniaT — Fri, 08 Dec 2006 17:45:22 GMT

New page

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 <tt> a = sum(a_i) </tt>):

# 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 <tt>a_sum</tt> 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