notebook.community

Edit and run

Initialization





In [1]:



    


import numpy as np
import matplotlib
from matplotlib import pyplot as plt
import ensemble as we
import ssad
import util
from imp import reload



    











In [2]:



    


reload(we)
reload(ssad)
reload(util)



    


















    
Out[2]:








<module 'util' from '/home/max/Documents/Research/wessa-s2014/python/util.py'>

Uniform Paving





In [2]:



    


paving = we.UniformPaving([-10, -3], [15, 9], [5, 4])



    











In [13]:



    


paving.get_bin_num([1, 1]) # 9



    


















    
Out[13]:








9

















In [14]:



    


paving.get_bin_num([1, 4]) # 10



    


















    
Out[14]:








10

















In [15]:



    


paving.get_bin_num([1, 3]) # 10



    


















    
Out[15]:








10

















In [16]:



    


paving.get_bin_num([0, 0]) # 9



    


















    
Out[16]:








9

















In [17]:



    


paving.get_bin_num([4, -4]) # 8



    


















    
Out[17]:








8

















In [18]:



    


paving.get_bin_num([4, -4000]) # 8



    


















    
Out[18]:








8

















In [19]:



    


paving.get_bin_num([5, -4]) # 12



    


















    
Out[19]:








12

















In [20]:



    


paving.get_bin_num([-10, -4]) # 0



    


















    
Out[20]:








0

















In [21]:



    


paving.get_bin_num([-10000, -4]) # 0



    


















    
Out[21]:








0

















In [22]:



    


paving.get_bin_num([15, 9]) # 19



    


















    
Out[22]:








19

















In [23]:



    


paving.get_bin_num([15, 6]) # 19



    


















    
Out[23]:








19

















In [24]:



    


paving.get_bin_num([14, 6]) # 19



    


















    
Out[24]:








19

















In [25]:



    


paving.get_bin_num([0, 9]) # 11



    


















    
Out[25]:








11

















In [3]:



    


coord_set = np.array([[1, 1],
                      [1, 4],
                      [1, 3],
                      [0, 0],
                      [4, -4],
                      [4, -4000],
                      [5, -4],
                      [-10, -4],
                      [-10000, -4],
                      [15, 9],
                      [15, 6],
                      [14, 6],
                      [0, 9]])



    











In [4]:



    


paving.get_bin_num(coord_set)



    


















    
Out[4]:








array([ 9, 10, 10,  9,  8,  8, 12,  0,  0, 19, 19, 19, 11])

Some Edge Cases





In [45]:



    


paving = we.UniformPaving([0], [60], (60))



    











In [46]:



    


paving.get_bin_num(33)



    


















    
Out[46]:








33

















In [47]:



    


# This is a nice, quick, visual way to find any irregularities in the binning
xs = np.array(range(61))
xs = xs[:,np.newaxis]
plt.plot(xs, paving.get_bin_num(xs))
plt.show()



    











In [13]:



    


nbins = 30
binrange = (0, 60)
bin_xs, bin_width = np.linspace(*binrange, num=nbins, endpoint=False, retstep=True)
paving = we.UniformPaving(*binrange, bin_counts=nbins)
xs = np.array(range(61))
plt.plot(xs, paving.get_bin_num(xs[:,np.newaxis]))
plt.show()

Self-Implemented Binomial Coefficient





In [28]:



    


ssad._combinations(5, 6)



    


















    
Out[28]:








0

















In [29]:



    


ssad._combinations(5, 5)



    


















    
Out[29]:








1

















In [30]:



    


ssad._combinations(1000, 1000)



    


















    
Out[30]:








1

















In [31]:



    


ssad._combinations(1000, 6)



    


















    
Out[31]:








1368173298991500

















In [32]:



    


ssad._combinations(1000, 994)



    


















    
Out[32]:








1368173298991500

















In [33]:



    


ssad._combinations(30, 0)



    


















    
Out[33]:








1

















In [37]:



    


ssad._combinations(30, -2)



    


















    
Out[37]:








0

















In [38]:



    


ssad._combinations(0, 0)



    


















    
Out[38]:








1

Trajectory Cloning (with weights)





In [15]:



    


species = ['X']
rxns = [ssad.Reaction([1], [-1], 1),  # Production
        ssad.Reaction([0], [+1], 25)] # Degradation
trj1 = we.WeightedTrajectory([0], rxns, 1.0)



    











In [16]:



    


trj1.run_dynamics(20)
states = np.asarray(trj1.hist_states)
times = np.asarray(trj1.hist_times)



    











In [17]:



    


plt.plot(times, states[:,0])
plt.show()



    











In [18]:



    


trj_clones = trj1.clone(3)



    











In [19]:



    


trj_clones



    


















    
Out[19]:








[<ensemble.WeightedTrajectory at 0x7fa074c90b38>,
 <ensemble.WeightedTrajectory at 0x7fa061defc88>,
 <ensemble.WeightedTrajectory at 0x7fa061e17668>]

















In [20]:



    


print([clone.weight for clone in trj_clones])



    


















    






[0.25, 0.25, 0.25]

















In [21]:



    


trj1.weight



    


















    
Out[21]:








0.25

















In [22]:



    


type(trj_clones[0])



    


















    
Out[22]:








ensemble.WeightedTrajectory

















In [24]:



    


states1 = np.asarray(trj_clones[1].hist_states)
times1 = np.asarray(trj_clones[1].hist_times)
plt.plot(times1, states1[:,0])
plt.show()

Make sure the trajectories aren't aliased to each other.





In [25]:



    


trj1.run_dynamics(20)



    











In [30]:



    


states0 = np.asarray(trj1.hist_states)
times0 = np.asarray(trj1.hist_times)
plt.plot(times0, states0[:,0])
plt.show()



    











In [29]:



    


states1 = np.asarray(trj_clones[2].hist_states)
times1 = np.asarray(trj_clones[2].hist_times)
plt.plot(times1, states1[:,0])
plt.show()

History Pruning

Set up a system with delayed reactions.





In [141]:



    


species = ['X']
rxns = [ssad.Reaction([0], [+1], 100),
        ssad.Reaction([1], [-1], 3),
        ssad.Reaction([1], [-1], 5, delay=20.0)]
trj = ssad.Trajectory([0], rxns)



    











In [142]:



    


trj.run_dynamics(60)



    











In [143]:



    


states = np.asarray(trj.hist_states)
times = np.asarray(trj.hist_times)
plt.plot(times, states[:,0])
plt.show()



    











In [144]:



    


trj.time



    


















    
Out[144]:








60.0

















In [145]:



    


trj.prune_history()



    











In [146]:



    


states = np.asarray(trj.hist_states)
times = np.asarray(trj.hist_times)
plt.plot(times, states[:,0])
plt.show()



    











In [147]:



    


trj.last_rxn_time



    


















    
Out[147]:








59.99968628625786

















In [148]:



    


min(trj.hist_times)



    


















    
Out[148]:








39.99344721303941

















In [150]:



    


trj.last_rxn_time - 20



    


















    
Out[150]:








39.99968628625786

Ensemble Bin Population Control

Growing





In [171]:



    


species = ['X']
rxns = [ssad.Reaction([1], [-1], 1),  # Production
        ssad.Reaction([0], [+1], 25)] # Degradation
trj1 = we.WeightedTrajectory([0], rxns, 1.0)
paving = we.UniformPaving([0], [50], 2)
ens = we.Ensemble(10, paving, (3, 5), [trj1])



    











In [172]:



    


print(ens.bins[0][0])



    


















    






Trajectory at time 0.0 . Current state: [0]. Weight 1.0.

















In [173]:



    


ens.run_step()



    











In [174]:



    


ens.bins



    


















    
Out[174]:








defaultdict(<function Ensemble._recompute_bins.<locals>.<lambda> at 0x7f98c3ad6440>, {0: [<ensemble.WeightedTrajectory object at 0x7f98c3a85a50>, <ensemble.WeightedTrajectory object at 0x7f98c3b46150>], 1: [<ensemble.WeightedTrajectory object at 0x7f98c3b2d890>]})

















In [175]:



    


for traj in ens:
    print(traj)



    


















    






Trajectory at time 10.0 . Current state: [21]. Next reaction scheduled for time 10.006347944735348. Weight 0.25.
Trajectory at time 10.0 . Current state: [24]. Next reaction scheduled for time 10.019015875139718. Weight 0.25.
Trajectory at time 10.0 . Current state: [32]. Next reaction scheduled for time 10.001066234502385. Weight 0.5.

















In [176]:



    


ens._resample()



    











In [177]:



    


for traj in ens:
    print(traj)



    


















    






Trajectory at time 10.0 . Current state: [21]. Next reaction scheduled for time 10.006347944735348. Weight 0.125.
Trajectory at time 10.0 . Current state: [24]. Next reaction scheduled for time 10.019015875139718. Weight 0.25.
Trajectory at time 10.0 . Current state: [21]. Next reaction scheduled for time 10.006347944735348. Weight 0.125.
Trajectory at time 10.0 . Current state: [32]. Next reaction scheduled for time 10.001066234502385. Weight 0.125.
Trajectory at time 10.0 . Current state: [32]. Next reaction scheduled for time 10.001066234502385. Weight 0.25.
Trajectory at time 10.0 . Current state: [32]. Next reaction scheduled for time 10.001066234502385. Weight 0.125.

















In [178]:



    


sum(traj.weight for traj in ens)



    


















    
Out[178]:








1.0

Shrinking





In [179]:



    


species = ['X']
rxns = [ssad.Reaction([1], [-1], 1),  # Production
        ssad.Reaction([0], [+1], 25)] # Degradation
init_trjs = []
for idx in range(14):
    init_trjs.append(we.WeightedTrajectory([0], rxns, 1/14))
paving = we.UniformPaving([0], [50], 2)
ens = we.Ensemble(10, paving, (3, 5), init_trjs)



    











In [180]:



    


ens.bins[0]



    


















    
Out[180]:








[<ensemble.WeightedTrajectory at 0x7f98c39d3f50>,
 <ensemble.WeightedTrajectory at 0x7f98c39d3f90>,
 <ensemble.WeightedTrajectory at 0x7f98c39d3fd0>,
 <ensemble.WeightedTrajectory at 0x7f98c396a050>,
 <ensemble.WeightedTrajectory at 0x7f98c396a090>,
 <ensemble.WeightedTrajectory at 0x7f98c396a0d0>,
 <ensemble.WeightedTrajectory at 0x7f98c396a110>,
 <ensemble.WeightedTrajectory at 0x7f98c396a150>,
 <ensemble.WeightedTrajectory at 0x7f98c396a190>,
 <ensemble.WeightedTrajectory at 0x7f98c396a1d0>,
 <ensemble.WeightedTrajectory at 0x7f98c396a210>,
 <ensemble.WeightedTrajectory at 0x7f98c396a250>,
 <ensemble.WeightedTrajectory at 0x7f98c396a290>,
 <ensemble.WeightedTrajectory at 0x7f98c396a2d0>]

















In [181]:



    


ens._run_dynamics_all()



    











In [182]:



    


ens._recompute_bins()



    











In [183]:



    


for bin_no, trjs in ens.bins.items():
    for traj in trjs:
        print(' '.join(("Bin", str(bin_no) + ".", str(traj))))



    


















    






Bin 0. Trajectory at time 10.0 . Current state: [18]. Next reaction scheduled for time 10.0103690017322. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [24]. Next reaction scheduled for time 10.050623596973466. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [21]. Next reaction scheduled for time 10.017628845032627. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [21]. Next reaction scheduled for time 10.006253704543637. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [24]. Next reaction scheduled for time 10.000897320430624. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [24]. Next reaction scheduled for time 10.0052217458396. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [17]. Next reaction scheduled for time 10.072085704197564. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [26]. Next reaction scheduled for time 10.003276832413649. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [26]. Next reaction scheduled for time 10.01471000348363. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [28]. Next reaction scheduled for time 10.027017278519278. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [29]. Next reaction scheduled for time 10.002456325144493. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [25]. Next reaction scheduled for time 10.05292670292451. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [27]. Next reaction scheduled for time 10.004609185578333. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [27]. Next reaction scheduled for time 10.020622096796496. Weight 0.07142857142857142.

















In [184]:



    


ens._resample()



    











In [185]:



    


for bin_no, trjs in ens.bins.items():
    for traj in trjs:
        print(' '.join(("Bin", str(bin_no) + ".", str(traj))))



    


















    






Bin 0. Trajectory at time 10.0 . Current state: [24]. Next reaction scheduled for time 10.050623596973466. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [21]. Next reaction scheduled for time 10.006253704543637. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [21]. Next reaction scheduled for time 10.017628845032627. Weight 0.14285714285714285.
Bin 0. Trajectory at time 10.0 . Current state: [24]. Next reaction scheduled for time 10.000897320430624. Weight 0.07142857142857142.
Bin 0. Trajectory at time 10.0 . Current state: [17]. Next reaction scheduled for time 10.072085704197564. Weight 0.14285714285714285.
Bin 1. Trajectory at time 10.0 . Current state: [26]. Next reaction scheduled for time 10.01471000348363. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [29]. Next reaction scheduled for time 10.002456325144493. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [28]. Next reaction scheduled for time 10.027017278519278. Weight 0.14285714285714285.
Bin 1. Trajectory at time 10.0 . Current state: [25]. Next reaction scheduled for time 10.05292670292451. Weight 0.07142857142857142.
Bin 1. Trajectory at time 10.0 . Current state: [27]. Next reaction scheduled for time 10.020622096796496. Weight 0.14285714285714285.

















In [186]:



    


sum(trj.weight for trj in ens)



    


















    
Out[186]:








0.9999999999999998

Bin Population Counter





In [60]:



    


A = 100.0
B = 3.0
C = 1.0
tau = 20.0
rxns = [ssad.Reaction([0], [+1], A),
        ssad.Reaction([1], [-1], B),
        ssad.Reaction([1], [-1], C, delay=tau)]



    











In [61]:



    


init_trjs = []
ntrajs = 100
nbins = 30
binrange = (0, 60)
bin_xs, bin_width = np.linspace(*binrange, num=nbins, endpoint=False, retstep=True)
paving = we.UniformPaving(*binrange, bin_counts=nbins)
for idx in range(ntrajs):
    init_state = np.random.randint(*binrange)
    init_time = np.random.random_sample() * rxns[2].delay
    init_trjs.append(we.WeightedTrajectory([init_state], rxns, 1.0 / ntrajs, init_time=init_time))



    











In [62]:



    


ens = we.Ensemble(0.1, paving, (10, 10), init_trjs)



    











In [63]:



    


print(ens.get_bin_counts())
ens.run_step()
print(ens.get_bin_counts())
ens._resample()
print(ens.get_bin_counts())



    


















    






[ 5.  3.  4.  4.  2.  6.  1.  4.  4.  2.  1.  1.  7.  5.  3.  3.  3.  5.
  3.  3.  2.  2.  3.  2.  3.  7.  6.  2.  2.  2.]
[ 10.  10.  10.  10.  10.   8.  11.  10.  11.  10.  10.  10.  10.  10.  10.
  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.]
[ 10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.
  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.  10.]

















In [64]:



    


plt.bar(bin_xs, ens.get_pdist())
plt.show()



    











In [65]:



    


ens.time



    


















    
Out[65]:








0.1

















In [66]:



    


ens.run_time(20.0)



    


















    
Out[66]:








20.100000000000016

















In [67]:



    


plt.bar(bin_xs, ens.get_pdist())
plt.show()



    











In [70]:



    


ens.get_bin_counts()



    


















    
Out[70]:








array([  0.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,
        10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.,
        10.,  10.,  10.,  10.,  10.,  10.,  10.,  10.])

















In [69]:



    


ens._resample()

History Sampling





In [3]:



    


rxns = [ssad.Reaction([0], [+1], 100),
        ssad.Reaction([1], [-1], 3),
        ssad.Reaction([1], [-1], 5, delay=20.0)]
trj = ssad.Trajectory([0], rxns)



    











In [4]:



    


trj.run_dynamics(1)



    











In [12]:



    


states = np.asarray(trj.hist_states)
times = np.asarray(trj.hist_times)
plt.step(times, states[:,0], where='post')
plt.show()



    











In [6]:



    


trj.hist_states[1]



    


















    
Out[6]:








array([1])

















In [7]:



    


trj.hist_times[1]



    


















    
Out[7]:








0.002830865171941392

















In [8]:



    


trj.hist_times[2]



    


















    
Out[8]:








0.004359918437150452

















In [10]:



    


print(trj.get_hist_index(0.001))
print(trj.sample_state(0.001))



    


















    






0
[0]

















In [13]:



    


print(trj.get_hist_index(0.003))
print(trj.sample_state(0.003))



    


















    






1
[1]

















In [14]:



    


print(trj.get_hist_index(0.0044))
print(trj.sample_state(0.0044))



    


















    






2
[2]

Joint Probability Integral

Construct an artificial trajectory with a history whose joint distribution is known





In [3]:



    


paving = we.UniformPaving([0], [7], 7)
rxns = [ssad.Reaction([], [+1], 1.0)]
trj = ssad.Trajectory([1], rxns)
trj._execute_rxn(rxns[0], 4)
trj._execute_rxn(rxns[0], 6)
trj._execute_rxn(rxns[0], 8)
trj._execute_rxn(rxns[0], 10)
trj._execute_rxn(rxns[0], 12)
trj.time = 13.0



    











In [4]:



    


plt.step(trj.hist_times, trj.hist_states, where='post')
plt.show()



    











In [5]:



    


j_pdist = util.delay_joint_pdist(trj, 5.0, paving, paving)
print(j_pdist)



    


















    






[[ 0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.  0.  0.]
 [ 0.  2.  0.  0.  0.  0.  0.]
 [ 0.  1.  1.  0.  0.  0.  0.]
 [ 0.  0.  1.  1.  0.  0.  0.]
 [ 0.  0.  0.  1.  0.  0.  0.]]

















In [6]:



    


j_pdist[6, 3]



    


















    
Out[6]:








1.0

















In [7]:



    


np.sum(j_pdist)



    


















    
Out[7]:








8.0

Restarting capability:





In [88]:



    


species = ['X']
tau = 20.0
rxns = [ssad.Reaction([0], [+1], 100),
        ssad.Reaction([1], [-1], 3),
        ssad.Reaction([1], [-1], 5, delay=tau)]
trj = ssad.Trajectory([0], rxns)



    











In [89]:



    


trj.run_dynamics(60)



    











In [90]:



    


paving = we.UniformPaving([0], [60], 60)
joint_pdist = util.delay_joint_pdist(trj, tau, paving, paving)



    











In [91]:



    


np.sum(joint_pdist)



    


















    
Out[91]:








40.000000000000078

















In [92]:



    


last_upd_time = trj.time



    











In [93]:



    


trj.run_dynamics(40)



    











In [94]:



    


joint_pdist_all = util.delay_joint_pdist(trj, tau, paving, paving)
joint_pdist_new = util.delay_joint_pdist(trj, tau, paving, paving, from_time=last_upd_time, existing=joint_pdist)
np.any(joint_pdist_all - joint_pdist_new > 1.0E-10)



    


















    
Out[94]:








False

















In [97]:



    


np.sum(joint_pdist)



    


















    
Out[97]:








80.000000000000142

















In [98]:



    


last_upd_time = trj.time
trj.prune_history()
trj.run_dynamics(40)
joint_pdist_final = util.delay_joint_pdist(trj, tau, paving, paving, from_time=last_upd_time, existing=joint_pdist)



    











In [99]:



    


np.sum(joint_pdist)



    


















    
Out[99]:








120.00000000000011

















In [100]:



    


plt.pcolormesh(joint_pdist, cmap='hot')
plt.show()

Everything works as expected.

Joint Probability Tally

Note: Replaced by time integral, see above





In [31]:



    


import bisect



    











In [32]:



    


species = ['X']
tau = 20.0
rxns = [ssad.Reaction([0], [+1], 100),
        ssad.Reaction([1], [-1], 3),
        ssad.Reaction([1], [-1], 5, delay=tau)]
trj = ssad.Trajectory([0], rxns)



    











In [33]:



    


trj.run_dynamics(60)



    











In [34]:



    


states = np.asarray(trj.hist_states)
times = np.asarray(trj.hist_times)
plt.plot(times, states[:,0])
plt.show()



    











In [35]:



    


start_idx = bisect.bisect_right(trj.hist_times, tau)



    











In [36]:



    


start_idx



    


















    
Out[36]:








4030

















In [37]:



    


len(trj.hist_times)



    


















    
Out[37]:








13418

















In [38]:



    


paving = we.UniformPaving([0], [60], 60)
joint_pdist = util.delay_joint_pdist(trj, tau, paving, paving)



    











In [40]:



    


plt.matshow(joint_pdist)
plt.show()

Excellent, the basics seem to work! Now on to checkpoint-style sampling.





In [64]:



    


trj = ssad.Trajectory([0], rxns)



    











In [65]:



    


trj.run_dynamics(100)



    











In [66]:



    


len(trj.hist_times)



    


















    
Out[66]:








22486

















In [67]:



    


start_idx = bisect.bisect_right(trj.hist_times, tau)



    











In [68]:



    


pdist = util.delay_joint_pdist(trj, tau, paving, paving)



    











In [81]:



    


last_upd_time = trj.time



    











In [55]:



    


trj.prune_history()



    











In [82]:



    


trj.run_dynamics(100)



    











In [83]:



    


util.delay_joint_pdist(trj, tau, paving, paving, existing=pdist, from_time=last_upd_time)



    


















    
Out[83]:








array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])

















In [84]:



    


np.sum(pdist)



    


















    
Out[84]:








85383.0

















In [85]:



    


len(trj.hist_times[start_idx:])



    


















    
Out[85]:








85383

















In [87]:



    


plt.pcolormesh(pdist, cmap='hot')
plt.show()

Good, works just as expected.





In [ ]:

Content source: max-veit/wessa-delay

Similar notebooks:

notebook.community | gallery | about