{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 4. How to Run a Simulation\n", "\n", "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`." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from ecell4_base.core import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.1. How to Setup a Simulator\n", "\n", "Except for the initialization (so-called constructor function) with arguments specific to the algorithm, all `Simulator`s have the same APIs." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from ecell4_base import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before constructing a `Simulator`, parepare a `Model` and a `World` corresponding to the type of `Simulator`." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from ecell4 import species_attributes, reaction_rules, get_model\n", "\n", "with species_attributes():\n", " A | B | C | {'D': 1, 'radius': 0.005}\n", "\n", "with reaction_rules():\n", " A + B == C | (0.01, 0.3)\n", "\n", "m = get_model()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "w1 = gillespie.World()\n", "w2 = ode.World()\n", "w3 = spatiocyte.World()\n", "w4 = bd.World()\n", "w5 = meso.World()\n", "w6 = egfrd.World()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Simulator` requires both `Model` and `World` in this order at the construction." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "sim1 = gillespie.Simulator(w1, m)\n", "sim2 = ode.Simulator(w2, m)\n", "sim3 = spatiocyte.Simulator(w3, m)\n", "sim4 = bd.Simulator(w4, m)\n", "sim5 = meso.Simulator(w5, m)\n", "sim6 = egfrd.Simulator(w6, m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you bind the `Model` to the `World`, you need only the `World` to create a `Simulator`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "w1.bind_to(m)\n", "w2.bind_to(m)\n", "w3.bind_to(m)\n", "w4.bind_to(m)\n", "w5.bind_to(m)\n", "w6.bind_to(m)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "sim1 = gillespie.Simulator(w1)\n", "sim2 = ode.Simulator(w2)\n", "sim3 = spatiocyte.Simulator(w3)\n", "sim4 = bd.Simulator(w4)\n", "sim5 = meso.Simulator(w5)\n", "sim6 = egfrd.Simulator(w6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, the `Model` and `World` bound to a `Simulator` can be drawn from `Simulator` in the way below:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", " \n", " \n", " \n", " \n", " \n" ] } ], "source": [ "print(sim1.model(), sim1.world())\n", "print(sim2.model(), sim2.world())\n", "print(sim3.model(), sim3.world())\n", "print(sim4.model(), sim4.world())\n", "print(sim5.model(), sim5.world())\n", "print(sim6.model(), sim6.world())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After updating the `World` by yourself, you must initialize the internal state of a `Simulator` before running simulation." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "w1.add_molecules(Species('C'), 60)\n", "w2.add_molecules(Species('C'), 60)\n", "w3.add_molecules(Species('C'), 60)\n", "w4.add_molecules(Species('C'), 60)\n", "w5.add_molecules(Species('C'), 60)\n", "w6.add_molecules(Species('C'), 60)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "sim1.initialize()\n", "sim2.initialize()\n", "sim3.initialize()\n", "sim4.initialize()\n", "sim5.initialize()\n", "sim6.initialize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For algorithms with a fixed step interval, the `Simulator` also requires `dt`." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "sim2.set_dt(1e-6) # ode.Simulator. This is optional\n", "sim4.set_dt(1e-6) # bd.Simulator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.2. Running Simulations\n", "\n", "For running simulations, `Simulator` provides two APIs, `step` and `run`.\n", "\n", "`step()` advances a simulation for the time that the `Simulator` expects, `next_time()`." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0 0.03782982587603311 0.03782982587603311\n", "0.0 1e-06 1e-06\n", "0.0 6.666666666666667e-05 6.666666666666667e-05\n", "0.0 1e-06 1e-06\n", "0.0 0.0025188519074043213 0.0025188519074043213\n", "0.0 0.0 0.0\n" ] } ], "source": [ "print(sim1.t(), sim1.next_time(), sim1.dt())\n", "print(sim2.t(), sim2.next_time(), sim2.dt()) # => (0.0, 1e-6, 1e-6)\n", "print(sim3.t(), sim3.next_time(), sim3.dt())\n", "print(sim4.t(), sim4.next_time(), sim4.dt()) # => (0.0, 1e-6, 1e-6)\n", "print(sim5.t(), sim5.next_time(), sim5.dt())\n", "print(sim6.t(), sim6.next_time(), sim6.dt()) # => (0.0, 0.0, 0.0)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "sim1.step()\n", "sim2.step()\n", "sim3.step()\n", "sim4.step()\n", "sim5.step()\n", "sim6.step()" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.03782982587603311\n", "1e-06\n", "6.666666666666667e-05\n", "1e-06\n", "0.0025188519074043213\n", "0.0\n" ] } ], "source": [ "print(sim1.t())\n", "print(sim2.t()) # => 1e-6\n", "print(sim3.t())\n", "print(sim4.t()) # => 1e-6\n", "print(sim5.t())\n", "print(sim6.t()) # => 0.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`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." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[(, )]\n", "[]\n", "[]\n", "[]\n", "[]\n" ] } ], "source": [ "print(sim1.last_reactions())\n", "# print(sim2.last_reactions())\n", "print(sim3.last_reactions())\n", "print(sim4.last_reactions())\n", "print(sim5.last_reactions())\n", "print(sim6.last_reactions())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`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`." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True 0.1531251186177207\n", "True 2e-06\n", "True 0.00013333333333333334\n", "True 2e-06\n", "True 0.003502192717675309\n", "True 0.0\n" ] } ], "source": [ "print(sim1.step(1.0), sim1.t())\n", "print(sim2.step(1.0), sim2.t())\n", "print(sim3.step(1.0), sim3.t())\n", "print(sim4.step(1.0), sim4.t())\n", "print(sim5.step(1.0), sim5.t())\n", "print(sim6.step(1.0), sim6.t())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To run a simulation just until the time, `upto`, call `step(upto)` while it returns `True`." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "while sim1.step(1.0): pass\n", "while sim2.step(0.001): pass\n", "while sim3.step(0.001): pass\n", "while sim4.step(0.001): pass\n", "while sim5.step(1.0): pass\n", "while sim6.step(0.001): pass" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0\n", "0.001\n", "0.001\n", "0.001\n", "1.0\n", "0.001\n" ] } ], "source": [ "print(sim1.t()) # => 1.0\n", "print(sim2.t()) # => 0.001\n", "print(sim3.t()) # => 0.001\n", "print(sim4.t()) # => 0.001\n", "print(sim5.t()) # => 1.0\n", "print(sim6.t()) # => 0.001" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is just what `run` does. `run(tau)` advances a simulation upto `t()+tau`." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "sim1.run(1.0)\n", "sim2.run(0.001)\n", "sim3.run(0.001)\n", "sim4.run(0.001)\n", "sim5.run(1.0)\n", "sim6.run(0.001)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.0\n", "0.002\n", "0.002\n", "0.002\n", "2.0\n", "0.002\n" ] } ], "source": [ "print(sim1.t()) # => 2.0\n", "print(sim2.t()) # => 0.002\n", "print(sim3.t()) # => 0.002\n", "print(sim4.t()) # => 0.002\n", "print(sim5.t()) # => 2.0\n", "print(sim6.t()) # => 0.02" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`num_steps` returns the number of steps during the simulation." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "28\n", "2001\n", "28\n", "2001\n", "835\n", "738\n" ] } ], "source": [ "print(sim1.num_steps())\n", "print(sim2.num_steps())\n", "print(sim3.num_steps())\n", "print(sim4.num_steps())\n", "print(sim5.num_steps())\n", "print(sim6.num_steps())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.3. Capsulizing Algorithm into a Factory Class\n", "\n", "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.\n", "\n", "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.\n", "\n", "`Factory` just provides two functions, `world` and `simulator`." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "def singlerun(f, m):\n", " w = f.world(Real3(1, 1, 1))\n", " w.bind_to(m)\n", " w.add_molecules(Species('C'), 60)\n", " \n", " sim = f.simulator(w)\n", " sim.run(0.01)\n", " print(sim.t(), w.num_molecules(Species('C')))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`singlerun` above is free from the algorithm. Thus, by just switching `Factory`, you can easily compare the results." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.01 60\n", "0.01 59\n", "0.01 60\n", "0.01 60\n", "0.01 60\n", "0.01 60\n" ] } ], "source": [ "singlerun(gillespie.Factory(), m)\n", "singlerun(ode.Factory(), m)\n", "singlerun(spatiocyte.Factory(), m)\n", "singlerun(bd.Factory(bd_dt_factor=1), m)\n", "singlerun(meso.Factory(), m)\n", "singlerun(egfrd.Factory(), m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When you need to provide several parameters to initialize `World` or `Simulator`, `run_simulation` also accepts `Factory` instead of `solver`." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.01, 0.0, 0.0, 60.0]\n", "[0.01, 0.17972919304002502, 0.17972919304002502, 59.820270806960046]\n", "[0.01, 1.0, 1.0, 59.0]\n", "[0.01, 0.0, 0.0, 60.0]\n", "[0.01, 1.0, 1.0, 59.0]\n", "[0.01, 1.0, 1.0, 59.0]\n" ] } ], "source": [ "from ecell4.util import run_simulation\n", "print(run_simulation(0.01, model=m, y0={'C': 60}, solver=gillespie.Factory()).as_array()[-1])\n", "print(run_simulation(0.01, model=m, y0={'C': 60}, solver=ode.Factory()).as_array()[-1])\n", "print(run_simulation(0.01, model=m, y0={'C': 60}, solver=spatiocyte.Factory()).as_array()[-1])\n", "print(run_simulation(0.01, model=m, y0={'C': 60}, solver=bd.Factory(bd_dt_factor=1)).as_array()[-1])\n", "print(run_simulation(0.01, model=m, y0={'C': 60}, solver=meso.Factory()).as_array()[-1])\n", "print(run_simulation(0.01, model=m, y0={'C': 60}, solver=egfrd.Factory()).as_array()[-1])" ] } ], "metadata": { "celltoolbar": "Edit Metadata", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.1" } }, "nbformat": 4, "nbformat_minor": 1 }