Initialization



In [1]:

    
%pylab









    



Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib



In [2]:

    
mpl.rcParams.update({'font.family': 'serif',
                     'font.size': 12,
                     'font.sans-serif': [],})



In [3]:

    
import sys
from collections import defaultdict
from imp import reload

import numpy.fft as fft
from mpl_toolkits.mplot3d import Axes3D
from scipy import optimize
#import prettyplotlib

import ssad
import ensemble as we



In [4]:

    
reload(ssad)
reload(we)
from ssad import Reaction, Trajectory



In [5]:

    
# From http://nbviewer.ipython.org/github/ipython/ipython/blob/master/examples/notebooks/Progress%20Bars.ipynb
class ProgressBar:
    def __init__(self, iterations):
        self.iterations = iterations
        self.prog_bar = '[]'
        self.fill_char = '='
        self.width = 50
        self.__update_amount(0)

    def animate(self, iter):
        print('\r', self, end='')
        sys.stdout.flush()
        self.update_iteration(iter + 1)

    def update_iteration(self, elapsed_iter):
        self.__update_amount((elapsed_iter / float(self.iterations)) * 100.0)
        self.prog_bar += '  %d of %s complete' % (elapsed_iter, self.iterations)

    def __update_amount(self, new_amount):
        percent_done = int(round((new_amount / 100.0) * 100.0))
        all_full = self.width - 2
        num_hashes = int(round((percent_done / 100.0) * all_full))
        self.prog_bar = '[' + self.fill_char * num_hashes + ' ' * (all_full - num_hashes) + ']'
        pct_place = (len(self.prog_bar) // 2) - len(str(percent_done))
        pct_string = '%d%%' % percent_done
        self.prog_bar = self.prog_bar[0:pct_place] + \
            (pct_string + self.prog_bar[pct_place + len(pct_string):])

    def __str__(self):
        return str(self.prog_bar)

Model Definition



In [93]:

    
species = ['X']



In [12]:

    
A = 100.0
B = 4.5
C = 5.0
tau = 1.0
rxns = [Reaction([0], [+1], A),
        Reaction([1], [-1], B),
        Reaction([1], [-1], C, delay=tau)]



In [26]:

    
# Calculate the expected equilibrium value
A / (B + C)









    Out[26]:





12.5

Compute the Hopf bifurcation point in the Langevin limit. Tentatively bracket to the first minimum of the cosine function,



In [13]:

    
ubound = np.sqrt((np.pi / tau)**2 + B**2)
bif_root = lambda C: np.cos(tau * np.sqrt(C**2 - B**2)) + B / C
C_bif = optimize.brentq(bif_root, B, ubound)



In [14]:

    
C_bif









    Out[14]:





5.204711579609008



In [21]:

    
C_range = np.linspace(B, ubound, 1000)
plot(C_range, bif_root(C_range))









    Out[21]:





[<matplotlib.lines.Line2D at 0x7f9d79e0bcf8>]



In [14]:

    
plot(C_range, -B / C_range)
plot(C_range, np.cos(tau * np.sqrt(C_range**2 - B**2)))









    Out[14]:





[<matplotlib.lines.Line2D at 0x7feb830f0710>]



In [95]:

    
# The starting state can be chosen either to induce oscillations or to remain at equilibrium.
trj = Trajectory([0], rxns)



In [96]:

    
trj.run_dynamics(120)



In [29]:

    
trj.run_dynamics(duration=None, stop_time=820.0)



In [97]:

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









    Out[97]:





[<matplotlib.lines.Line2D at 0x7f070888a588>]

With a weighted ensemble



In [98]:

    
init_trjs = []
init_trjs_zero = []
ntrajs = 40
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 = random.randint(*binrange)
    init_time = random.random_sample() * rxns[2].delay
    init_trjs.append(we.WeightedTrajectory([init_state], rxns, 1.0 / ntrajs, init_time=init_time))
    init_trjs_zero.append(we.WeightedTrajectory([0], rxns, 1.0 / ntrajs, init_time=init_time))



In [99]:

    
print([trj.state for trj in init_trjs])









    



[array([28]), array([20]), array([18]), array([1]), array([24]), array([4]), array([38]), array([54]), array([4]), array([24]), array([30]), array([15]), array([26]), array([33]), array([55]), array([39]), array([42]), array([16]), array([54]), array([43]), array([47]), array([26]), array([29]), array([47]), array([4]), array([31]), array([2]), array([2]), array([36]), array([28]), array([35]), array([45]), array([18]), array([0]), array([2]), array([7]), array([37]), array([10]), array([7]), array([27])]

Quick sanity check...



In [33]:

    
we_test = we.Ensemble(10, paving, (10, 10), init_trjs)
niter_test = 5
pbar = ProgressBar(niter_test)
pbar.animate(0)
for stidx in range(niter_test):
    we_test.run_step()
    pbar.animate(stidx + 1)









    



 [==========            20%                       ]  1 of 5 complete





    



---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-33-0390fb01b441> in <module>()
      4 pbar.animate(0)
      5 for stidx in range(niter_test):
----> 6     we_test.run_step()
      7     pbar.animate(stidx + 1)

/home/max/Documents/Research/wessa-s2014/python/ensemble.py in run_step(self, resample, prune)
    239         if resample:
    240             self._resample()
--> 241         self._run_dynamics_all()
    242         self._recompute_bins()
    243         #self._record_state()

/home/max/Documents/Research/wessa-s2014/python/ensemble.py in _run_dynamics_all(self)
    381         stop_time = self.time + self.step_time
    382         for traj in self:
--> 383             traj.run_dynamics(duration=None, stop_time=stop_time)
    384         self.time = stop_time
    385 

/home/max/Documents/Research/wessa-s2014/python/ssad.py in run_dynamics(self, duration, max_steps, stop_time)
    220                 resume = False
    221             else:
--> 222                 next_rxn, wait_time = self._sample_next_reaction()
    223                 next_rxn_time = self.time + wait_time
    224 

/home/max/Documents/Research/wessa-s2014/python/ssad.py in _sample_next_reaction(self)
    281                     self.time - rxn.delay, use_init_state = True)
    282             propensities[ridx] = rxn.calc_propensity(state)
--> 283         prop_csum = np.cumsum(propensities)
    284         total_prop = prop_csum[-1]
    285         wait_time = exponential(1.0 / total_prop)

/usr/lib/python3.4/site-packages/numpy/core/fromnumeric.py in cumsum(a, axis, dtype, out)
   1981     except AttributeError:
   1982         return _wrapit(a, 'cumsum', axis, dtype, out)
-> 1983     return cumsum(axis, dtype, out)
   1984 
   1985 

KeyboardInterrupt:



In [93]:

    
plot_trjs = list(we_test)
sample_step = int(len(plot_trjs) / 7)
plot_trjs = plot_trjs[::sample_step]
fig = figure()
ax = fig.gca()
for trj in plot_trjs:
    trj_times = np.array(trj.hist_times)
    trj_states = np.array(trj.hist_states)
    prettyplotlib.plot(trj_times, trj_states[:,0], axes=ax)



In [100]:

    
step_time = 2.0
run_time = 60.0
niter = int(run_time / step_time)
we_res = we.Ensemble(step_time, paving, (4, 4), init_trjs_zero)
we_nores = we.Ensemble(step_time, paving, (4, 4), init_trjs_zero)



In [101]:

    
pbar = ProgressBar(niter)
pdist_res = np.empty((nbins, niter))



In [102]:

    
pbar.animate(0)
for stidx in range(niter):
    we_res.run_step(prune=False)
    pdist_res[:,stidx] = we_res.get_pdist()
    pbar.animate(stidx + 1)









    



 [=====================100%=======================]  30 of 30 complete

Let's plot some sample trajectories, to make sure the phases are suffiently randomized:



In [107]:

    
plot_trjs = list(we_res)
sample_step = int(len(plot_trjs) / 5)
plot_trjs = plot_trjs[::sample_step]
fig = figure()
ax = fig.gca()
for trj in plot_trjs:
    trj_times = np.array(trj.hist_times)
    trj_states = np.array(trj.hist_states)
    plot(trj_times, trj_states[:,0], axes=ax)



In [56]:

    
def plot_trj(trj):
    times = np.array(trj.hist_times)
    states = np.array(trj.hist_states)
    plot(times, states[:,0])



In [58]:

    
plot_trj(plot_trjs[2])



In [59]:

    
pdist_res = we_res.get_pdist()
bar(bin_xs, pdist_res)









    Out[59]:





<Container object of 30 artists>



In [60]:

    
we_res.get_bin_counts()









    Out[60]:





array([  5.,   3.,   0.,   2.,   0.,   0.,   0.,   0.,   0.,   0.,   3.,
         2.,   4.,   3.,   8.,   6.,   7.,  13.,   4.,   2.,   6.,   0.,
         0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.])



In [104]:

    
pbar = ProgressBar(niter)
pdist_nores = np.empty((nbins, niter))



In [105]:

    
pbar.animate(0)
for stidx in range(niter):
    we_nores.run_step(prune=False, resample=False)
    pbar.animate(stidx + 1)









    



 [=====================100%=======================]  30 of 30 complete



In [112]:

    
we_nores._recompute_bins()
we_nores._resample()



In [113]:

    
for stidx in range(10):
    we_nores.run_step(prune=False, resample=False)

As before, make sure the phases are suffiently randomized:



In [114]:

    
plot_trjs = list(we_nores)
plot_trjs = plot_trjs[::8]
fig = figure()
ax = fig.gca()
for trj in plot_trjs:
    trj_times = np.array(trj.hist_times)
    trj_states = np.array(trj.hist_states)
    plot(trj_times, trj_states[:,0], axes=ax)



In [70]:

    
plot(trj_times, trj_states[:,0])









    Out[70]:





[<matplotlib.lines.Line2D at 0x7f07080d6c50>]



In [115]:

    
plot_trj(plot_trjs[-1])



In [116]:

    
plot_trjs[-1].sample_state(50.0)









    Out[116]:





array([27])



In [90]:

    
plot_trjs[-1].init_time









    Out[90]:





60.0



In [85]:

    
test_trj = plot_trjs[-2]
test_trj.time - tau









    Out[85]:





60.0



In [86]:

    
test_trj.sample_state(test_trj.time - tau)









    Out[86]:





array([34])



In [80]:

    
test_trj.run_dynamics(10)



In [88]:

    
plot_trj(test_trj)

Now with ensembles of ensembles



In [16]:

    
step_time = 2.0
niter = 30
run_time = step_time * niter
nens = 30
wes_res = []
wes_nores = []
for eidx in range(nens):
    wes_res.append(we.Ensemble(step_time, paving, (2, 2), init_trjs))
    wes_nores.append(we.Ensemble(step_time, paving, (2, 2), init_trjs))



In [17]:

    
prog_ctr = 0
pbar = ProgressBar(nens * niter * 2)



In [18]:

    
pbar.animate(0)
for stidx in range(niter):
    for ens in wes_res:
        ens.run_step(resample=True)
        prog_ctr += 1
        pbar.animate(prog_ctr)
    for ens in wes_nores:
        ens.run_step(resample=False)
        prog_ctr += 1
        pbar.animate(prog_ctr)









    



 [======================57%==                     ]  1024 of 1800 complete





    



---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-18-dd055e5a0bba> in <module>()
      2 for stidx in range(niter):
      3     for ens in wes_res:
----> 4         ens.run_step(resample=True)
      5         prog_ctr += 1
      6         pbar.animate(prog_ctr)

/home/max/Documents/Research/wessa-s2014/python/ensemble.py in run_step(self, resample, prune)
    239         if resample:
    240             self._resample()
--> 241         self._run_dynamics_all()
    242         self._recompute_bins()
    243         #self._record_state()

/home/max/Documents/Research/wessa-s2014/python/ensemble.py in _run_dynamics_all(self)
    345         stop_time = self.time + self.step_time
    346         for traj in self:
--> 347             traj.run_dynamics(duration=None, stop_time=stop_time)
    348         self.time = stop_time
    349 

/home/max/Documents/Research/wessa-s2014/python/ssad.py in run_dynamics(self, duration, max_steps, stop_time)
    224 
    225             # Execute the reaction, if possible
--> 226             if self._can_run_rxn(next_rxn):
    227                 if next_rxn_time > stop_time:
    228                     self._save_run_state(next_rxn, next_rxn_time)

/home/max/Documents/Research/wessa-s2014/python/ssad.py in _can_run_rxn(self, rxn)
    259 
    260         """
--> 261         return np.all(self.state + rxn.state_vec >= 0)
    262 
    263     def _sample_next_reaction(self):

/usr/lib/python3.4/site-packages/numpy/core/fromnumeric.py in all(a, axis, out, keepdims)
   1908     arr = asanyarray(a)
   1909 
-> 1910     try:
   1911         return arr.all(axis=axis, out=out, keepdims=keepdims)
   1912     except TypeError:

KeyboardInterrupt:



In [37]:

    
ens_pdist_res = np.vstack(ens.get_pdist() for ens in wes_res)
ens_pdist_nores = np.vstack(ens.get_pdist() for ens in wes_nores)



In [21]:

    
sum(len(ens) for ens in wes_res)









    Out[21]:





888



In [24]:

    
sum(len(trj.hist_states) for trj in wes_res[0])









    Out[24]:





181510



In [38]:

    
np.savez('../results/del_deg_ens_of_ens.npz', pdist_res=ens_pdist_res, pdist_nores=ens_pdist_nores)



In [39]:

    
ens_pdist_res.shape









    Out[39]:





(20, 20)



In [40]:

    
ens_pdist_res_avg = np.average(ens_pdist_res, axis=0)
ens_pdist_nores_avg = np.average(ens_pdist_nores, axis=0)
ens_pdist_res_sdm = np.std(ens_pdist_res, axis=0) / sqrt(nens)
ens_pdist_nores_sdm = np.std(ens_pdist_nores, axis=0) / sqrt(nens)



In [43]:

    
fig = figure()
ax = fig.gca()
ax.bar(bin_xs, ens_pdist_nores_avg, yerr=ens_pdist_nores_sdm, alpha=0.6, color='red', ecolor='red', label='Without Resampling')
ax.bar(bin_xs, ens_pdist_res_avg, yerr=ens_pdist_res_sdm, alpha=0.6, color='blue', ecolor='blue', label='With Resampling')
ax.legend(loc='best')
ax.set_ylim(bottom=0.0)
ax.set_xlabel('Number of X')
ax.set_ylabel('Probability Density (per bin)')









    Out[43]:





<matplotlib.text.Text at 0x7f8e39536310>



In [44]:

    
fig.savefig('../results/delayed-deg-pdist.pgf')

Now, once again, the $\chi$ plot:



In [45]:

    
chis = (ens_pdist_res_avg - ens_pdist_nores_avg) / np.sqrt(ens_pdist_nores_sdm**2 + ens_pdist_res_sdm**2)
chis[np.abs(chis) > 1E9] = np.nan









    



-c:1: RuntimeWarning: invalid value encountered in true_divide



In [52]:

    
fig = figure()
ax = fig.gca()
ax.bar(bin_xs, chis)
ax.set_xlabel('Number of X')
ax.set_ylabel(r'Normalized Difference (Wtd - UnWtd)')
ax.axhline(0, color='black')
ax.grid()



In [53]:

    
fig.savefig('../results/delayed-deg-chi.pgf')

Power Spectrum



In [15]:

    
n_samples = 2**11
max_freq_num = int(n_samples / 2)



In [16]:

    
s_times = np.linspace(0, 200, n_samples)
s_states = trj.sample_state_seq(s_times)



In [17]:

    
amp_spec = fft.rfft(s_states, axis=1)
freqs = fft.fftfreq(s_times.size, d=(s_times[-1] - s_times[0])/n_samples)



In [18]:

    
plot(freqs[:max_freq_num], np.abs(amp_spec[0,:-1])**2)









    Out[18]:





[<matplotlib.lines.Line2D at 0x7f197b1b28d0>]



In [ ]: