#!/usr/bin/env python
# coding: utf-8

# # 8. More about 1. Brief Tour of E-Cell4 Simulations
# 
# Once you read through [1. Brief Tour of E-Cell4 Simulations](tutorial1.ipynb), it is NOT difficult to use `World` and `Simulator`.
# `volume` and `{'C': 60}` is equivalent of the `World` and solver is the `Simulator` below.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
from ecell4.prelude import *

with reaction_rules():
    A + B == C | (0.01, 0.3)

run_simulation(10.0, {'C': 60}, volume=1.0)


# Here we give you a breakdown for `run_simulation`.
# `run_simulation` use ODE simulator by default, so we create `ode.World` step by step.

# ## 8.1. Creating ODEWorld
# 
# You can create `World` like this.

# In[2]:


from ecell4_base.core import *
from ecell4_base import *


# In[3]:


w = ode.World(Real3(1, 1, 1))


# `Real3` is a coordinate vector.
# In this example, the first argument for `ode.World` constructor is a cube.
# Note that you can NOT use volume for `ode.World` argument, like `run_simulation` argument.
# 
# Now you created a cube box for simulation, next let's throw molecules into the cube.

# 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)


# Use `add_molecules` to add molecules, `remove_molecules` to remove molecules, `num_molecules` to know the number of molecules.
# First argument for each method is the `Species` you want to know.
# You can get current time by `t` method.
# However the number of molecules in ODE solver is real number, in these `_molecules` functions work only for integer number.
# When you handle real numbers in ODE, use `set_value` and `get_value`.

# ## 8.2. How to Use Real3
# 
# Before the detail of `Simulator`, we explaing more about `Real3`.

# 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)


# You can not print the contents in `Real3` object directly.
# You need to convert `Real3` to Python tuple or list once.

# 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


# You can use basic function like `dot_product`.
# Of course, you can convert `Real3` to numpy array too.

# In[7]:


import numpy
a = numpy.asarray(tuple(Real3(1, 2, 3)))
print(a)  # must print [ 1.  2.  3.]


# `Integer3` represents a triplet of integers.

# In[8]:


g = Integer3(1, 2, 3)
print(tuple(g))


# Of course, you can also apply simple arithmetics to `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


# ## 8.3. Creating and Running ODE Simulator
# 
# You can create a `Simulator` with `Model` and `World` like 

# 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)


# then call `run` method, the simulation will run.
# In this example the simulation runs for 10 seconds.
# 
# You can check the state of the `World` like this.

# In[11]:


print(w.t(), w.num_molecules(Species('C')))  # must return (10.0, 30)


# You can see that the number of the `Species` `C` decreases from 60 to 30.
# 
# `World` describes the state at a timepoint, so you can NOT tack the transition during the simulation with the `World`.
# To obtain the time-series result, use `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]]


# There are several types of `Observer`s for E-Cell4.
# `FixedIntervalNumberObserver` is the simplest `Observer` to obtain the time-series result.
# As its name suggests, this `Observer` records the number of molecules for each time-step.
# The 1st argument is the time-step, the 2nd argument is the molecule types.
# You can check the result with `data` method, but there is a shortcut for this.

# In[13]:


show(obs)


# This plots the time-series result easily.
# 
# We explained the internal of `run_simulation` function.
# When you change the `World` after creating the `Simulator`, you need to indicate it to `Simulator`.
# So do NOT forget to call `sim.initialize()` after that.

# ## 8.4. Switching the Solver
# 
# It is NOT difficult to switch the solver to stochastic method, as we showed `run_simulation`.

# In[14]:


from ecell4 import *

with reaction_rules():
    A + B == C | (0.01, 0.3)

m = get_model()

# ode.World -> gillespie.World
w = gillespie.World(Real3(1, 1, 1))
w.add_molecules(Species('C'), 60)

# ode.Simulator -> gillespie.Simulator
sim = gillespie.Simulator(w, m)
obs = FixedIntervalNumberObserver(0.1, ('A', 'C'))
sim.run(10.0, obs)

show(obs, step=True)


# `World` and `Simulator` never change the `Model` itself, so you can switch several `Simulator`s for one `Model`.