4. How to Run a Simulation¶
In sections 2 and 3, we explained the way to build a model and to setup the intial state. Now, it is the time to run a simulation. Corresponding to World
classes, six Simulator
classes are there: spatiocyte.SpatiocyteSimulator
, egfrd.EGFRDSimulator
, bd.BDSimulator
, meso.MesoscopicSimulator
, gillespie.GillespieSimulator
, and ode.ODESimulator
. Each Simulator
class only accepts the corresponding type of World
, but all of them allow the same Model
.
[1]:
from ecell4_base.core import *
4.1. How to Setup a Simulator¶
Except for the initialization (so-called constructor function) with arguments specific to the algorithm, all Simulator
s have the same APIs.
[2]:
from ecell4_base import *
Before constructing a Simulator
, parepare a Model
and a World
corresponding to the type of Simulator
.
[3]:
from ecell4 import species_attributes, reaction_rules, get_model
with species_attributes():
A | B | C | {'D': 1, 'radius': 0.005}
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
[4]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
Simulator
requires both Model
and World
in this order at the construction.
[5]:
sim1 = gillespie.Simulator(w1, m)
sim2 = ode.Simulator(w2, m)
sim3 = spatiocyte.Simulator(w3, m)
sim4 = bd.Simulator(w4, m)
sim5 = meso.Simulator(w5, m)
sim6 = egfrd.Simulator(w6, m)
If you bind the Model
to the World
, you need only the World
to create a Simulator
.
[6]:
w1.bind_to(m)
w2.bind_to(m)
w3.bind_to(m)
w4.bind_to(m)
w5.bind_to(m)
w6.bind_to(m)
[7]:
sim1 = gillespie.Simulator(w1)
sim2 = ode.Simulator(w2)
sim3 = spatiocyte.Simulator(w3)
sim4 = bd.Simulator(w4)
sim5 = meso.Simulator(w5)
sim6 = egfrd.Simulator(w6)
Of course, the Model
and World
bound to a Simulator
can be drawn from Simulator
in the way below:
[8]:
print(sim1.model(), sim1.world())
print(sim2.model(), sim2.world())
print(sim3.model(), sim3.world())
print(sim4.model(), sim4.world())
print(sim5.model(), sim5.world())
print(sim6.model(), sim6.world())
<ecell4_base.core.Model object at 0x145cf5fee7b0> <ecell4_base.gillespie.GillespieWorld object at 0x145cf5ff1720>
<ecell4_base.core.Model object at 0x145cf5feeab0> <ecell4_base.ode.ODEWorld object at 0x145cf5ff16f0>
<ecell4_base.core.Model object at 0x145cf5feeab0> <ecell4_base.spatiocyte.SpatiocyteWorld object at 0x145cf5ff1720>
<ecell4_base.core.Model object at 0x145cf5feea30> <ecell4_base.bd.BDWorld object at 0x145cf5ff16f0>
<ecell4_base.core.Model object at 0x145cf5feeab0> <ecell4_base.meso.MesoscopicWorld object at 0x145cf5ff1738>
<ecell4_base.core.Model object at 0x145cf5feea30> <ecell4_base.egfrd.EGFRDWorld object at 0x145cf5ff16f0>
After updating the World
by yourself, you must initialize the internal state of a Simulator
before running simulation.
[9]:
w1.add_molecules(Species('C'), 60)
w2.add_molecules(Species('C'), 60)
w3.add_molecules(Species('C'), 60)
w4.add_molecules(Species('C'), 60)
w5.add_molecules(Species('C'), 60)
w6.add_molecules(Species('C'), 60)
[10]:
sim1.initialize()
sim2.initialize()
sim3.initialize()
sim4.initialize()
sim5.initialize()
sim6.initialize()
For algorithms with a fixed step interval, the Simulator
also requires dt
.
[11]:
sim2.set_dt(1e-6) # ode.Simulator. This is optional
sim4.set_dt(1e-6) # bd.Simulator
4.2. Running Simulation¶
For running simulation, Simulator
provides two APIs, step
and run
.
step()
advances a simulation for the time that the Simulator
expects, next_time()
.
[12]:
print(sim1.t(), sim1.next_time(), sim1.dt())
print(sim2.t(), sim2.next_time(), sim2.dt()) # => (0.0, 1e-6, 1e-6)
print(sim3.t(), sim3.next_time(), sim3.dt())
print(sim4.t(), sim4.next_time(), sim4.dt()) # => (0.0, 1e-6, 1e-6)
print(sim5.t(), sim5.next_time(), sim5.dt())
print(sim6.t(), sim6.next_time(), sim6.dt()) # => (0.0, 0.0, 0.0)
0.0 0.06981637440659177 0.06981637440659177
0.0 1e-06 1e-06
0.0 6.666666666666667e-05 6.666666666666667e-05
0.0 1e-06 1e-06
0.0 0.0024540697478422522 0.0024540697478422522
0.0 0.0 0.0
[13]:
sim1.step()
sim2.step()
sim3.step()
sim4.step()
sim5.step()
sim6.step()
[14]:
print(sim1.t())
print(sim2.t()) # => 1e-6
print(sim3.t())
print(sim4.t()) # => 1e-6
print(sim5.t())
print(sim6.t()) # => 0.0
0.06981637440659177
1e-06
6.666666666666667e-05
1e-06
0.0024540697478422522
0.0
last_reactions()
returns a list of pairs of ReactionRule
and ReactionInfo
which occurs at the last step. Each algorithm have its own implementation of ReactionInfo
. See help(module.ReactionInfo)
for details.
[15]:
print(sim1.last_reactions())
# print(sim2.last_reactions())
print(sim3.last_reactions())
print(sim4.last_reactions())
print(sim5.last_reactions())
print(sim6.last_reactions())
[(<ecell4_base.core.ReactionRule object at 0x145cf5ff1780>, <ecell4_base.gillespie.ReactionInfo object at 0x145cf5ff1798>)]
[]
[]
[]
[]
step(upto)
advances a simulation for next_time
if next_time
is less than upto
, or for upto
otherwise. step(upto)
returns whether the time does NOT reach the limit, upto
.
[16]:
print(sim1.step(1.0), sim1.t())
print(sim2.step(1.0), sim2.t())
print(sim3.step(1.0), sim3.t())
print(sim4.step(1.0), sim4.t())
print(sim5.step(1.0), sim5.t())
print(sim6.step(1.0), sim6.t())
True 0.2285050559395263
True 2e-06
True 0.00013333333333333334
True 2e-06
True 0.004146250248941816
True 0.0
To run a simulation just until the time, upto
, call step(upto)
while it returns True
.
[17]:
while sim1.step(1.0): pass
while sim2.step(0.001): pass
while sim3.step(0.001): pass
while sim4.step(0.001): pass
while sim5.step(1.0): pass
while sim6.step(0.001): pass
[18]:
print(sim1.t()) # => 1.0
print(sim2.t()) # => 0.001
print(sim3.t()) # => 0.001
print(sim4.t()) # => 0.001
print(sim5.t()) # => 1.0
print(sim6.t()) # => 0.001
1.0
0.001
0.001
0.001
1.0
0.001
This is just what run
does. run(tau)
advances a simulation upto t()+tau
.
[19]:
sim1.run(1.0)
sim2.run(0.001)
sim3.run(0.001)
sim4.run(0.001)
sim5.run(1.0)
sim6.run(0.001)
[20]:
print(sim1.t()) # => 2.0
print(sim2.t()) # => 0.002
print(sim3.t()) # => 0.002
print(sim4.t()) # => 0.002
print(sim5.t()) # => 2.0
print(sim6.t()) # => 0.02
2.0
0.002
0.002
0.002
2.0
0.002
num_steps
returns the number of steps during the simulation.
[21]:
print(sim1.num_steps())
print(sim2.num_steps())
print(sim3.num_steps())
print(sim4.num_steps())
print(sim5.num_steps())
print(sim6.num_steps())
29
2001
28
2001
921
599
4.3. Capsulizing Algorithm into a Factory Class¶
Owing to the portability of a Model
and consistent APIs of World
s and Simulator
s, it is very easy to write a script common in algorithms. However, when switching the algorithm, still we have to rewrite the name of classes in the code, one by one.
To avoid the trouble, E-Cell4 also provides a Factory
class for each algorithm. Factory
encapsulates World
and Simulator
with their arguments needed for the construction. By using Factory
class, your script could be portable and robust agaist changes in the algorithm.
Factory
just provides two functions, world
and simulator
.
[22]:
def singlerun(f, m):
w = f.world(Real3(1, 1, 1))
w.bind_to(m)
w.add_molecules(Species('C'), 60)
sim = f.simulator(w)
sim.run(0.01)
print(sim.t(), w.num_molecules(Species('C')))
singlerun
above is free from the algorithm. Thus, by just switching Factory
, you can easily compare the results.
[23]:
singlerun(gillespie.Factory(), m)
singlerun(ode.Factory(), m)
singlerun(spatiocyte.Factory(), m)
singlerun(bd.Factory(bd_dt_factor=1), m)
singlerun(meso.Factory(), m)
singlerun(egfrd.Factory(), m)
0.01 60
0.01 59
0.01 60
0.01 60
0.01 59
0.01 60
When you need to provide several parameters to initialize World
or Simulator
, run_simulation
also accepts Factory
instead of solver
.
[24]:
from ecell4.util import run_simulation
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', solver=gillespie.Factory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', solver=ode.Factory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', solver=spatiocyte.Factory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', solver=bd.Factory(bd_dt_factor=1))[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', solver=meso.Factory())[-1])
print(run_simulation(0.01, model=m, y0={'C': 60}, return_type='array', solver=egfrd.Factory())[-1])
[0.01, 0.0, 0.0, 60.0]
[0.01, 0.17972919304002502, 0.17972919304002502, 59.820270806960046]
[0.01, 1.0, 1.0, 59.0]
[0.01, 0.0, 0.0, 60.0]
[0.01, 1.0, 1.0, 59.0]
[0.01, 0.0, 0.0, 60.0]