セクション2と3では、モデルを構築し、初期状態を設定する方法について説明しました。
それでは、シミュレーションを実行してみましょう。
「Spatiocyte.SpatiocyteSimulator」、「egfrd.EGFRDSimulator」、「bd.BDSimulator」、「meso.MesoscopicSimulator」、「gillespie.GillespieSimulator」、および「ode.ODESimulator」の6つの「シミュレータ」クラスが存在します。
それぞれの SimulatorクラスはWorldの対応する型だけを受け入れますが、それらのすべてが同じ Modelを受け入れます。
In [1]:
from ecell4_base.core import *
In [2]:
from ecell4_base import *
Simulatorを構築する前に、Simulatorのタイプに対応するModelとWorldを用意してください。
In [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()
In [4]:
w1 = gillespie.World()
w2 = ode.World()
w3 = spatiocyte.World()
w4 = bd.World()
w5 = meso.World()
w6 = egfrd.World()
シミュレータは、構築時に下記の順序でモデルとワールドの両方を必要とします。
In [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)
ModelをWorldにバインドした場合は、シミュレータを作成するために要るのはWorldだけです。
In [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)
In [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)
もちろん、SimulatorにバインドされたModelとWorldは、Simulatorから以下のように引き出すことができます。
In [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())
Worldを自分で更新した後は、シミュレーションを実行する前にSimulatorの内部状態を初期化する必要があります。
In [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)
In [10]:
sim1.initialize()
sim2.initialize()
sim3.initialize()
sim4.initialize()
sim5.initialize()
sim6.initialize()
一定のステップ間隔を有するアルゴリズムの場合は上記に加えdtも必要です。
In [11]:
sim2.set_dt(1e-6) # ode.Simulator. This is optional
sim4.set_dt(1e-6) # bd.Simulator
In [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)
In [13]:
sim1.step()
sim2.step()
sim3.step()
sim4.step()
sim5.step()
sim6.step()
In [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
last_reactions()は、最後のステップで発生するReactionRuleとReactionInfoのペアのリストを返します。
各アルゴリズムにはReactionInfoの独自の実装があります。
詳細については、help(module.ReactionInfo)を参照してください。
In [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())
step(upto)は、next_timeがuptoより小さい場合はnext_timeのシミュレーションを進め、それ以外の場合はuptoまで実行します。
またstep(upto)は時間が上限(upto)に達していないかどうかを返します。
In [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())
upto の時間ちょうどまでシミュレーションをは走らせるには, それがTrueを返す間 step(upto) を呼びます。
In [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
In [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
これはちょうど run の動作と同じものです。
run(tau)はt()+tauの時間までシミュレーションを進めます。
In [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)
In [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
num_stepsは、シミュレーション中のステップ数を返します。
In [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())
Modelの移植性とWorldと Simulatorの一貫したAPIのおかげで、アルゴリズムに共通するスクリプトを書くのは非常に簡単です。
しかし、アルゴリズムを切り替えるときには、コード内のクラスの名前を1つずつ書き直す必要があります。
この問題を回避するために、E-Cell4はアルゴリズムごとに Factoryクラスも提供しています。 FactoryはWorldと Simulatorを構造化に必要な引数でカプセル化します。 Factoryクラスを使うことで、あなたのスクリプトはアルゴリズムの変更に対してポータブルかつロバストになります。
Factoryはworldと simulatorの2つの関数を提供します。
In [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はアルゴリズムから独立しています。
したがって、単にFactoryを切り替えるだけで、結果を簡単に比較することができます。
In [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)
WorldやSimulatorを初期化するためにいくつかのパラメータを指定する必要があるとき、 run_simulationはsolverの代わりに Factoryを受け取ります。
In [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])