Tutorial 9 (Observer)

This is a tutorial for E-Cell4. Here, we explain Observer classes for logging various types of data during simulation.



In [1]:

    
from ecell4.core import *
from ecell4 import gillespie

m = NetworkModel()
m.add_reaction_rule(create_unimolecular_reaction_rule(Species("A"), Species("B"), 1.0))
m.add_reaction_rule(create_unimolecular_reaction_rule(Species("B"), Species("A"), 0.5))

w = gillespie.GillespieWorld(Real3(1, 1, 1))
w.add_molecules(Species("A"), 60)

s = gillespie.GillespieSimulator(m, w)
s.initialize()

Prepare observers. Not necessarily needed to be created here.



In [2]:

    
obs1 = FixedIntervalNumberObserver(1.0, ["A"])
obs2 = NumberObserver(["A", "B"])
obs3 = TimingNumberObserver([0, 1e-3, 1e-2, 1e-1, 1e+0, 1e+1], ["A"])
obs4 = FixedIntervalHDF5Observer(5.0, "test%02d.h5")

Run simulation with observers. Warning: HDF5 is not supported as a default. World::save raises NotSupported error without the HDF5 feature.



In [3]:

    
s.run(16.0, [obs1, obs2, obs3])
# s.run(16.0, [obs1, obs2, obs3, obs4])

NumberObservers return data by calling the member function NumberObserver::data(). NumberObserver::targets() returns a list of species, which you specified as an argument of the constructor.



In [4]:

    
print(obs1.data())
print([sp.serial() for sp in obs1.targets()])









    



[[0.0, 60.0], [1.0, 31.0], [2.0, 22.0], [3.0, 21.0], [4.0, 24.0], [5.0, 22.0], [6.0, 22.0], [7.0, 24.0], [8.0, 20.0], [9.0, 20.0], [10.0, 23.0], [11.0, 23.0], [12.0, 15.0], [13.0, 19.0], [14.0, 22.0], [15.0, 18.0], [16.0, 17.0]]
[u'A']

The normal NumberObserver stores the numbers of given Species just after every steps. It doesn't work properly with ODESimulator because ODESimulator has no internal dt. FixedIntervalNumberObserver logs them with the fixed step interval. TimingNumberObserver allow you to give the times of logging as an argument of the constructor.



In [5]:

    
print(obs2.data()[: 5])
print(obs3.data())









    



[[0.0, 60.0, 0.0], [0.03863568120520103, 59.0, 1.0], [0.03964146514302787, 58.0, 2.0], [0.05802737572001471, 57.0, 3.0], [0.07918729219539744, 56.0, 4.0]]
[[0.0, 60.0], [0.001, 60.0], [0.01, 60.0], [0.1, 56.0], [1.0, 31.0], [10.0, 23.0]]

FixedIntervalHDF5Observer saves a World in the HDF5 format. FixedIntervalHDF5Observer::filename() returns the name of a file scheduled to be saved next. At most one format string like %02d is allowed to use a step count in the file name. When you do not use the format string, it overwrites the latest data to the file.



In [6]:

    
print(obs4.filename())









    



test00.h5

viz.plot_number_observer is the easiest way to plot a result in the NumberObserver, FixedIntervalNumberObserver and TimingNumberObserver. For the detailed usage of viz.plot_number_observer, see help(viz.plot_number_observer).



In [7]:

    
%matplotlib inline
from ecell4 import viz
viz.plot_number_observer(obs2)

Of course, you can analyze and plot the results by yourself.



In [8]:

    
%matplotlib inline
import numpy
import matplotlib.pylab as plt

labels = [sp.serial() for sp in obs2.targets()]
data = numpy.asarray(obs2.data()).T

plt.plot(data[0], data[2] / data[1], label="{1} / {0}".format(labels[0], labels[1]))
plt.plot([data[0][0], data[0][-1]], [2.0, 2.0], 'k--')
plt.xlabel("Time")
plt.ylabel("The ratio of numbers")
plt.legend(loc="best", shadow=True)
plt.show()

Warning: HDF5 is not supported as a default. Check HDF5 files.



In [9]:

    
# import glob
# print(glob.glob("test*.h5"))
# w.load('test03.h5')
# print(w.t(), w.num_molecules(Species("A")), w.num_molecules(Species("B")))

For logging positions of particles, FixedIntervalHDF5Observer, FixedIntervalCSVObserver and FixedIntervalTrajectoryObserver are available. Here is an example with the lattice module.



In [10]:

    
from ecell4 import lattice

m = NetworkModel()
m.add_species_attribute(Species("A", "0.005", "1"))

w = lattice.LatticeWorld(Real3(1, 1, 1), 0.005)
w.bind_to(m)
w.add_molecules(Species("A"), 3)

s = lattice.LatticeSimulator(w)
s.initialize()

The usage of FixedIntervalCSVObserver is almost same with that of FixedIntervalHDF5Observer. It saves positions (x, y, z) of particles with the radius (r) and serial number of Species (sid). FixedIntervalTrajectoryObserver keeps tracking of particles, which you specified with ParticleIDs. Once some ParticleID is lost for the reaction, it just stop to trace the particle any more.



In [11]:

    
obs5 = FixedIntervalCSVObserver(0.3, "test%02d.csv")
obs6 = FixedIntervalTrajectoryObserver(1e-3, [pid for pid, p in w.list_particles()])

The way to run a simulation with these Observers is the same with the above for NumberObservers.



In [12]:

    
s.run(1.0, [obs5, obs6])

Chek CSV files.



In [13]:

    
import glob
print(glob.glob("test*.csv"))
print(open("test00.csv").read())









    



['test00.csv', 'test03.csv', 'test02.csv', 'test01.csv']
x,y,z,r,sid
0.49806291436591293,0.20207259421636903,0.029999999999999999,0.0050000000000000001,0
0.65319726474218087,0.21650635094610965,0.115,0.0050000000000000001,0
0.66136223055145815,0.6957070743734991,0.89500000000000002,0.0050000000000000001,0

Generally, World assumes a periodic boundary for each plane. To avoid the big jump of a particle at the edge because of the boundary condition, FixedIntervalTrajectoryObserver tries to keep the shift of positions. Thus, the positions stored in the Observer are not necessarily limited in the cuboid given for the World. Of course, you can disable this option. See help(FixedIntervalTrajectoryObserver).



In [14]:

    
print([tuple(pos) for pos in obs6.data()[0][: 5]])
print([tuple(pos) for pos in obs6.data()[0][-5:]])









    



[(0.49806291436591293, 0.20207259421636903, 0.03), (0.4409081537009721, 0.1991858428704209, -0.025000000000000022), (0.4164132562731403, 0.11547005383792515, 0.09999999999999998), (0.4082482904638631, 0.11258330249197702, 0.135), (0.4082482904638631, 0.12124355652982141, 0.18)]
[(3.3102687007525358, -0.5236860279185587, 1.775), (3.3265986323710903, -0.4977052658050256, 1.77), (3.26944387170615, -0.5641005467618325, 1.7149999999999999), (3.3021037349432585, -0.6247223250267433, 1.67), (3.3265986323710903, -0.6622500925240689, 1.605)]

A user can easily plot the trajectory of particles by using viz.plot_trajectory. For the detailed usage, see the Visualizer section and help(viz.plot_trajectory).



In [15]:

    
# viz.plot_trajectory(obs6)



In [ ]: