一度1. E-Cell4を用いたシミュレーションの概要を読んだ後は、World と Simulator を使うのは難しくありません。
volume と {'C' : 60} は World に相当し、ソルバーは以下の Simulator に相当します。
In [1]:
%matplotlib inline
from ecell4 import *
with reaction_rules():
A + B == C | (0.01, 0.3)
run_simulation(10.0, {'C': 60}, volume=1.0)
ここでは run_simulation がその内部で行っていることを分解して示します。
run_simulation はデフォルトでODEシミュレータを使用するので、ode.World を段階的に作成してみましょう。
In [2]:
from ecell4_base.core import *
from ecell4_base import *
In [3]:
w = ode.World(Real3(1, 1, 1))
Real3 は座標ベクトルです。
この例では、 ode.World コンストラクタの最初の引数は立方体を表しています。
run_simulation 引数のように、ode.World 引数にはvolumeを使用できないことに注意してください。
今度はシミュレーション用の立方体の空間を作成した後、分子をその立方体中に投げ込みましょう。
In [4]:
w = ode.World(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
print(w.t(), w.num_molecules(Species('C'))) # must return (0.0, 60)
分子を追加するには add_molecules、分子を削除するにはremove_molecules、分子の数を知るには num_molecules を使います。
各メソッドの第一引数は、あなたが知りたい Species です。
t メソッドで現在時刻を取得できます。
しかし、ODEソルバーの分子数は実数であり、これらの _molecules 関数では整数に対してのみ動作します。
ODEで実数を扱うときは、 set_value と get_value を使います。
In [5]:
pos = Real3(1, 2, 3)
print(pos) # must print like <ecell4.core.Real3 object at 0x7f44e118b9c0>
print(tuple(pos)) # must print (1.0, 2.0, 3.0)
Real3 オブジェクトの内容を直接printすることはできません。
Real3 をPythonタプルまたはリストに一度変換する必要があります。
In [6]:
pos1 = Real3(1, 1, 1)
x, y, z = pos[0], pos[1], pos[2]
pos2 = pos1 + pos1
pos3 = pos1 * 3
pos4 = pos1 / 5
print(length(pos1)) # must print 1.73205080757
print(dot_product(pos1, pos3)) # must print 9.0
dot_product のような基本的な関数も使えます。
もちろん、Real3 をnumpyの配列に変換することもできます。
In [7]:
import numpy
a = numpy.asarray(tuple(Real3(1, 2, 3)))
print(a) # must print [ 1. 2. 3.]
Integer3 は整数のトリプレットになります。
In [8]:
g = Integer3(1, 2, 3)
print(tuple(g))
もちろん、単純な算術演算を Integer3 に適用することもできます。
In [9]:
print(tuple(Integer3(1, 2, 3) + Integer3(4, 5, 6))) # => (5, 7, 9)
print(tuple(Integer3(4, 5, 6) - Integer3(1, 2, 3))) # => (3, 3, 3)
print(tuple(Integer3(1, 2, 3) * 2)) # => (2, 4, 6)
print(dot_product(Integer3(1, 2, 3), Integer3(4, 5, 6))) # => 32
print(length(Integer3(1, 2, 3))) # => 3.74165738677
In [10]:
with reaction_rules():
A + B > C | 0.01 # equivalent to create_binding_reaction_rule
C > A + B | 0.3 # equivalent to create_unbinding_reaction_rule
m = get_model()
sim = ode.Simulator(w, m)
sim.run(10.0)
run メソッドを呼び出すと、シミュレーションが実行されます。
この例では、シミュレーションは10秒間実行されます。
また World の状態は下記のようにしてチェックできます。
In [11]:
print(w.t(), w.num_molecules(Species('C'))) # must return (10.0, 30)
Species C の数が60から30に減少していることがわかります。
World はある時点での状態を表しているため、World でシミュレーションの遷移を見ることはできません。
時系列の結果を得るには、 Observer を使います。
In [12]:
w = ode.World(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
sim = ode.Simulator(w, m)
obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)
print(obs.data()) # must return [[0.0, 0.0, 60.0], ..., [10.0, 29.994446899691276, 30.005553100308752]]
E-Cell4にはいくつかのタイプの Observer があります。
FixedIntervalNumberObserver は、時系列の結果を得るための最もシンプルな Observer です。
その名前が示唆するように、この Observer は各時間ステップの分子の数を記録します。
第1引数は時間ステップ、第2引数は分子種です。
その結果は data メソッドで確認できますが、これにはショートカットがあります。
In [13]:
show(obs)
上記は時系列の結果を簡単にプロットします。
ここでは run_simulation 関数の内部について説明しました。
Simulator を作成した後に World を変更した場合は Simulator にそれを示す必要があります。
またそのあとで sim.initialize() を呼び出すことを忘れないでください。
In [14]:
from ecell4 import *
with reaction_rules():
A + B == C | (0.01, 0.3)
m = get_model()
# ode.World -> gillespie.World
w = gillespie.GillespieWorld(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)
# ode.Simulator -> gillespie.Simulator
sim = gillespie.GillespieSimulator(w, m)
obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)
viz.plot_number_observer(obs)
World と Simulator は Model 自体を決して変更しないので、一つの Model に対して 複数の Simulator を差し替えることができます。