{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. How to Setup the Initial Condition\n",
    "\n",
    "Here, we explain the basics of `World` classes. In E-Cell4, six types of World classes are supported now: `spatiocyte.SpatiocyteWorld`, `egfrd.EGFRDWorld`, `bd.BDWorld`, `meso.MesoscopicWorld`, `gillespie.GillespieWorld`, and `ode.ODEWorld`.\n",
    "\n",
    "In the most of softwares, the initial condition is supposed to be a part of `Model`. But, in E-Cell4, the initial condition must be set up as `World` separately from `Model`. `World` stores an information about the state, such as a current time, the number of molecules, coordinate of molecules, structures, and random number generator, at a time point. Meanwhile, `Model` contains the type of interactions between molecules and the common properties of molecules. `Model` is reusable among algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ecell4_base.core import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1. Common APIs of World\n",
    "\n",
    "Even though `World` describes the spatial representation specific to the corresponding algorithm, it has compatible APIs. In this section, we introduce the common interfaces of the six `World` classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ecell4_base import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`World` classes accept different sets of arguments in the constructor, which determine the parameters specific to the algorithm. However, at least, all `World` classes can be instantiated only with their size, named `edge_lengths`. The type of `edge_lengths` is `Real3`, which represents a triplet of `Real`s. In E-Cell4, all 3-dimensional positions are treated as a `Real3`. See also [8. More about 1. Brief Tour of E-Cell4 Simulations](tutorial8.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_lengths = Real3(1, 2, 3)\n",
    "w1 = gillespie.World(edge_lengths)\n",
    "w2 = ode.World(edge_lengths)\n",
    "w3 = spatiocyte.World(edge_lengths)\n",
    "w4 = bd.World(edge_lengths)\n",
    "w5 = meso.World(edge_lengths)\n",
    "w6 = egfrd.World(edge_lengths)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`World` has getter methods for the size and volume."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1.0, 2.0, 3.0) 6.0\n",
      "(1.0, 2.0, 3.0) 6.0\n",
      "(1.0124557603503803, 2.0091789367798976, 3.0) 6.102614364352381\n",
      "(1.0, 2.0, 3.0) 6.0\n",
      "(1.0, 2.0, 3.0) 6.0\n",
      "(1.0, 2.0, 3.0) 6.0\n"
     ]
    }
   ],
   "source": [
    "print(tuple(w1.edge_lengths()), w1.volume())\n",
    "print(tuple(w2.edge_lengths()), w2.volume())\n",
    "print(tuple(w3.edge_lengths()), w3.volume())\n",
    "print(tuple(w4.edge_lengths()), w4.volume())\n",
    "print(tuple(w5.edge_lengths()), w5.volume())\n",
    "print(tuple(w6.edge_lengths()), w6.volume())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`spatiocyte.World` (`w3`) would have a bit larger volume to fit regular hexagonal close-packed (HCP) lattice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's add molecules into the World. Here, you must give `Species` attributed with **radius** and **D** to tell the shape of molecules. In the example below **0.0025** corresponds to **radius** and **1** to **D**. Positions of the molecules are randomly determined in the `World` if needed. **10** in **add_molecules** function represents the number of molecules to be added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp1 = Species(\"A\", 0.0025, 1)\n",
    "w1.add_molecules(sp1, 10)\n",
    "w2.add_molecules(sp1, 10)\n",
    "w3.add_molecules(sp1, 10)\n",
    "w4.add_molecules(sp1, 10)\n",
    "w5.add_molecules(sp1, 10)\n",
    "w6.add_molecules(sp1, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a model is bound to the world, you do not need to rewrite the **radius** and **D** once set in `Species` (unless you want to change it)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = NetworkModel()\n",
    "m.add_species_attribute(Species(\"A\", 0.0025, 1))\n",
    "m.add_species_attribute(Species(\"B\", 0.0025, 1))\n",
    "\n",
    "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)\n",
    "\n",
    "w1.add_molecules(Species(\"B\"), 20)\n",
    "w2.add_molecules(Species(\"B\"), 20)\n",
    "w3.add_molecules(Species(\"B\"), 20)\n",
    "w4.add_molecules(Species(\"B\"), 20)\n",
    "w5.add_molecules(Species(\"B\"), 20)\n",
    "w6.add_molecules(Species(\"B\"), 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, `remove_molecules` and `num_molecules_exact` are also available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "w1.remove_molecules(Species(\"B\"), 5)\n",
    "w2.remove_molecules(Species(\"B\"), 5)\n",
    "w3.remove_molecules(Species(\"B\"), 5)\n",
    "w4.remove_molecules(Species(\"B\"), 5)\n",
    "w5.remove_molecules(Species(\"B\"), 5)\n",
    "w6.remove_molecules(Species(\"B\"), 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 10 10 10 10 10\n",
      "15 15 15 15 15 15\n"
     ]
    }
   ],
   "source": [
    "print(w1.num_molecules_exact(Species(\"A\")), w2.num_molecules_exact(Species(\"A\")), w3.num_molecules_exact(Species(\"A\")), w4.num_molecules_exact(Species(\"A\")), w5.num_molecules_exact(Species(\"A\")), w6.num_molecules_exact(Species(\"A\")))\n",
    "print(w1.num_molecules_exact(Species(\"B\")), w2.num_molecules_exact(Species(\"B\")), w3.num_molecules_exact(Species(\"B\")), w4.num_molecules_exact(Species(\"B\")), w5.num_molecules_exact(Species(\"B\")), w6.num_molecules_exact(Species(\"B\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike `num_molecules_exact`, `num_molecules` returns the numbers that match a given `Species` in rule-based fashion. When all `Species` in the `World` has a single `UnitSpecies` with no sites, `num_molecules` is same with `num_molecules_exact`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 10 10 10 10 10\n",
      "15 15 15 15 15 15\n"
     ]
    }
   ],
   "source": [
    "print(w1.num_molecules(Species(\"A\")), w2.num_molecules(Species(\"A\")), w3.num_molecules(Species(\"A\")), w4.num_molecules(Species(\"A\")), w5.num_molecules(Species(\"A\")), w6.num_molecules(Species(\"A\")))\n",
    "print(w1.num_molecules(Species(\"B\")), w2.num_molecules(Species(\"B\")), w3.num_molecules(Species(\"B\")), w4.num_molecules(Species(\"B\")), w5.num_molecules(Species(\"B\")), w6.num_molecules(Species(\"B\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`World` holds its simulation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 0.0 0.0 0.0 0.0 0.0\n",
      "1.0 1.0 1.0 1.0 1.0 1.0\n"
     ]
    }
   ],
   "source": [
    "print(w1.t(), w2.t(), w3.t(), w4.t(), w5.t(), w6.t())\n",
    "w1.set_t(1.0)\n",
    "w2.set_t(1.0)\n",
    "w3.set_t(1.0)\n",
    "w4.set_t(1.0)\n",
    "w5.set_t(1.0)\n",
    "w6.set_t(1.0)\n",
    "print(w1.t(), w2.t(), w3.t(), w4.t(), w5.t(), w6.t())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, you can `save` and `load` the state of a `World` into/from a HDF5 file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "w1.save(\"gillespie.h5\")\n",
    "w2.save(\"ode.h5\")\n",
    "w3.save(\"spatiocyte.h5\")\n",
    "w4.save(\"bd.h5\")\n",
    "w5.save(\"meso.h5\")\n",
    "w6.save(\"egfrd.h5\")\n",
    "del w1, w2, w3, w4, w5, w6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 (1.0, 1.0, 1.0) 1.0 0 0\n",
      "0.0 (1.0, 1.0, 1.0) 1.0 0 0\n",
      "0.0 (1.0124557603503803, 1.0045894683899488, 1.0) 1.0171023940587303 0 0\n",
      "0.0 (1.0, 1.0, 1.0) 1.0 0 0\n",
      "0.0 (1.0, 1.0, 1.0) 1.0 0 0\n",
      "0.0 (1.0, 1.0, 1.0) 1.0 0 0\n"
     ]
    }
   ],
   "source": [
    "w1 = gillespie.World()\n",
    "w2 = ode.World()\n",
    "w3 = spatiocyte.World()\n",
    "w4 = bd.World()\n",
    "w5 = meso.World()\n",
    "w6 = egfrd.World()\n",
    "print(w1.t(), tuple(w1.edge_lengths()), w1.volume(), w1.num_molecules(Species(\"A\")), w1.num_molecules(Species(\"B\")))\n",
    "print(w2.t(), tuple(w2.edge_lengths()), w2.volume(), w2.num_molecules(Species(\"A\")), w2.num_molecules(Species(\"B\")))\n",
    "print(w3.t(), tuple(w3.edge_lengths()), w3.volume(), w3.num_molecules(Species(\"A\")), w3.num_molecules(Species(\"B\")))\n",
    "print(w4.t(), tuple(w4.edge_lengths()), w4.volume(), w4.num_molecules(Species(\"A\")), w4.num_molecules(Species(\"B\")))\n",
    "print(w5.t(), tuple(w5.edge_lengths()), w5.volume(), w5.num_molecules(Species(\"A\")), w5.num_molecules(Species(\"B\")))\n",
    "print(w6.t(), tuple(w6.edge_lengths()), w6.volume(), w6.num_molecules(Species(\"A\")), w6.num_molecules(Species(\"B\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 (1.0, 2.0, 3.0) 6.0 10 15\n",
      "1.0 (1.0, 2.0, 3.0) 6.0 10 15\n",
      "1.0 (1.0124557603503803, 2.0091789367798976, 3.0) 6.102614364352381 10 15\n",
      "1.0 (1.0, 2.0, 3.0) 6.0 10 15\n",
      "1.0 (1.0, 2.0, 3.0) 6.0 10 15\n",
      "1.0 (1.0, 2.0, 3.0) 6.0 10 15\n"
     ]
    }
   ],
   "source": [
    "w1.load(\"gillespie.h5\")\n",
    "w2.load(\"ode.h5\")\n",
    "w3.load(\"spatiocyte.h5\")\n",
    "w4.load(\"bd.h5\")\n",
    "w5.load(\"meso.h5\")\n",
    "w6.load(\"egfrd.h5\")\n",
    "print(w1.t(), tuple(w1.edge_lengths()), w1.volume(), w1.num_molecules(Species(\"A\")), w1.num_molecules(Species(\"B\")))\n",
    "print(w2.t(), tuple(w2.edge_lengths()), w2.volume(), w2.num_molecules(Species(\"A\")), w2.num_molecules(Species(\"B\")))\n",
    "print(w3.t(), tuple(w3.edge_lengths()), w3.volume(), w3.num_molecules(Species(\"A\")), w3.num_molecules(Species(\"B\")))\n",
    "print(w4.t(), tuple(w4.edge_lengths()), w4.volume(), w4.num_molecules(Species(\"A\")), w4.num_molecules(Species(\"B\")))\n",
    "print(w5.t(), tuple(w5.edge_lengths()), w5.volume(), w5.num_molecules(Species(\"A\")), w5.num_molecules(Species(\"B\")))\n",
    "print(w6.t(), tuple(w6.edge_lengths()), w6.volume(), w6.num_molecules(Species(\"A\")), w6.num_molecules(Species(\"B\")))\n",
    "del w1, w2, w3, w4, w5, w6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the `World` classes also accept a HDF5 file path as an unique argument of the constructor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "print(gillespie.World(\"gillespie.h5\").t())\n",
    "print(ode.World(\"ode.h5\").t())\n",
    "print(spatiocyte.World(\"spatiocyte.h5\").t())\n",
    "print(bd.World(\"bd.h5\").t())\n",
    "print(meso.World(\"meso.h5\").t())\n",
    "print(egfrd.World(\"egfrd.h5\").t())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2. How to Get Molecule Positions\n",
    "\n",
    "`World` also has the common functions to access the coordinates of the molecules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": [
    "First, you can place a molecule at the certain position with `new_particle`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp1 = Species(\"A\", 0.0025, 1)\n",
    "pos = Real3(0.5, 0.5, 0.5)\n",
    "(pid1, p1), suc1 = w1.new_particle(sp1, pos)\n",
    "(pid2, p2), suc2 = w2.new_particle(sp1, pos)\n",
    "pid3 = w3.new_particle(sp1, pos)\n",
    "(pid4, p4), suc4 = w4.new_particle(sp1, pos)\n",
    "(pid5, p5), suc5 = w5.new_particle(sp1, pos)\n",
    "(pid6, p6), suc6 = w6.new_particle(sp1, pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`new_particle` returns a particle created and whether it's succeeded or not. The resolution in representation of molecules differs. For example, `GillespieWorld` has almost no information about the coordinate of molecules. Thus, it simply ignores the given position, and just counts up the number of molecules here.\n",
    "\n",
    "`ParticleID` is a pair of `Integer`s named `lot` and `serial`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "print(pid6.lot(), pid6.serial())\n",
    "print(pid6 == ParticleID((0, 1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particle simulators, i.e. `spatiocyte`, `bd` and `egfrd`, provide an interface to access a particle by its id. `has_particle` returns if a particles exists or not for the given `ParticleID`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# print(w1.has_particle(pid1))\n",
    "# print(w2.has_particle(pid2))\n",
    "print(w3.has_particle(pid3))  # => True\n",
    "print(w4.has_particle(pid4))  # => True\n",
    "# print(w5.has_particle(pid5))\n",
    "print(w6.has_particle(pid6))  # => True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After checking the existence, you can get the partcle by `get_particle` as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pid1, p1 = w1.get_particle(pid1)\n",
    "# pid2, p2 = w2.get_particle(pid2)\n",
    "pid3, p3 = w3.get_particle(pid3)\n",
    "pid4, p4 = w4.get_particle(pid4)\n",
    "# pid5, p5 = w5.get_particle(pid5)\n",
    "pid6, p6 = w6.get_particle(pid6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "`Particle` consists of `species`, `position`, `radius` and `D`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A (0.5062278801751902, 0.5080682368868706, 0.5) 0.0025 1.0\n",
      "A (0.5, 0.5, 0.5) 0.0025 1.0\n",
      "A (0.5, 0.5, 0.5) 0.0025 1.0\n"
     ]
    }
   ],
   "source": [
    "# print(p1.species().serial(), tuple(p1.position()), p1.radius(), p1.D())\n",
    "# print(p2.species().serial(), tuple(p2.position()), p2.radius(), p2.D())\n",
    "print(p3.species().serial(), tuple(p3.position()), p3.radius(), p3.D())\n",
    "print(p4.species().serial(), tuple(p4.position()), p4.radius(), p4.D())\n",
    "# print(p5.species().serial(), tuple(p5.position()), p5.radius(), p5.D())\n",
    "print(p6.species().serial(), tuple(p6.position()), p6.radius(), p6.D())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of `spatiocyte`, a particle position is automatically round to the center of the voxel nearest to the given position.\n",
    "\n",
    "You can even move the position of the particle. `update_particle` replace the particle specified with the given `ParticleID` with the given `Particle` and return `False`. If no corresponding particle is found, create new particle and return `True`. If you give a `Particle` with the different type of `Species`, the `Species` of the `Particle` will be also changed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "False\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "newp = Particle(sp1, Real3(0.3, 0.3, 0.3), 0.0025, 1)\n",
    "# print(w1.update_particle(pid1, newp))\n",
    "# print(w2.update_particle(pid2, newp))\n",
    "print(w3.update_particle(pid3, newp))\n",
    "print(w4.update_particle(pid4, newp))\n",
    "# print(w5.update_particle(pid5, newp))\n",
    "print(w6.update_particle(pid6, newp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`list_particles` and `list_particles_exact` return a list of pairs of `ParticleID` and `Particle` in the `World`. `World` automatically makes up for the gap with random numbers. For example, `GillespieWorld` returns a list of positions randomly distributed in the `World` size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(<ecell4_base.core.ParticleID object at 0x15330c38d1f0>, <ecell4_base.core.Particle object at 0x15330c38d228>)]\n",
      "[(<ecell4_base.core.ParticleID object at 0x15330c38d1b8>, <ecell4_base.core.Particle object at 0x15330c38d260>)]\n",
      "[(<ecell4_base.core.ParticleID object at 0x15330c38d1f0>, <ecell4_base.core.Particle object at 0x15330c38d308>)]\n",
      "[(<ecell4_base.core.ParticleID object at 0x15330c38d1b8>, <ecell4_base.core.Particle object at 0x15330c38d228>)]\n",
      "[(<ecell4_base.core.ParticleID object at 0x15330c38d1f0>, <ecell4_base.core.Particle object at 0x15330c38d260>)]\n"
     ]
    }
   ],
   "source": [
    "print(w1.list_particles_exact(sp1))\n",
    "# print(w2.list_particles_exact(sp1))  # ODEWorld has no member named list_particles\n",
    "print(w3.list_particles_exact(sp1))\n",
    "print(w4.list_particles_exact(sp1))\n",
    "print(w5.list_particles_exact(sp1))\n",
    "print(w6.list_particles_exact(sp1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can remove a specific particle with `remove_particle`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "False\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "# w1.remove_particle(pid1)\n",
    "# w2.remove_particle(pid2)\n",
    "w3.remove_particle(pid3)\n",
    "w4.remove_particle(pid4)\n",
    "# w5.remove_particle(pid5)\n",
    "w6.remove_particle(pid6)\n",
    "\n",
    "# print(w1.has_particle(pid1))\n",
    "# print(w2.has_particle(pid2))\n",
    "print(w3.has_particle(pid3))  # => False\n",
    "print(w4.has_particle(pid4))  # => False\n",
    "# print(w5.has_particle(pid5))\n",
    "print(w6.has_particle(pid6))  # => False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3. Lattice-based Coordinate\n",
    "\n",
    "In addition to the common interface, each `World` can have their own interfaces. As an example, we explain methods to handle lattice-based coordinate here. `SpatiocyteWorld` is based on a space discretized to hexiagonal close packing lattices, `LatticeSpace`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = spatiocyte.World(Real3(1, 2, 3), voxel_radius=0.01)\n",
    "w.bind_to(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The size of a single lattice, called `Voxel`, can be obtained by `voxel_radius()`. `SpatiocyteWorld` has methods to get the numbers of rows, columns, and layers. These sizes are automatically calculated based on the given `edge_lengths` at the construction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01\n",
      "(64, 152, 118)\n",
      "1147904\n"
     ]
    }
   ],
   "source": [
    "print(w.voxel_radius())  # => 0.01\n",
    "print(tuple(w.shape()))  # => (64, 152, 118)\n",
    "# print(w.col_size(), w.row_size(), w.layer_size())  # => (64, 152, 118)\n",
    "print(w.size())  # => 1147904 = 64 * 152 * 118"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A position in the lattice-based space is treated as an `Integer3`, column, row and layer, called a global coordinate. Thus, `SpatiocyteWorld` provides the function to convert the `Real3` into a lattice-based coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# p1 = Real3(0.5, 0.5, 0.5)\n",
    "# g1 = w.position2global(p1)\n",
    "# p2 = w.global2position(g1)\n",
    "# print(tuple(g1))  # => (31, 25, 29)\n",
    "# print(tuple(p2))  # => (0.5062278801751902, 0.5080682368868706, 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In `SpatiocyteWorld`, the global coordinate is translated to a single integer. It is just called a coordinate. You can also treat the coordinate as in the same way with a global coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# p1 = Real3(0.5, 0.5, 0.5)\n",
    "# c1 = w.position2coordinate(p1)\n",
    "# p2 = w.coordinate2position(c1)\n",
    "# g1 = w.coord2global(c1)\n",
    "# print(c1)  # => 278033\n",
    "# print(tuple(p2))  # => (0.5062278801751902, 0.5080682368868706, 0.5)\n",
    "# print(tuple(g1))  # => (31, 25, 29)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these coordinates, you can handle a `Voxel`, which represents a `Particle` object. Instead of `new_particle`, `new_voxel` provides the way to create a new `Voxel` with a coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# c1 = w.position2coordinate(Real3(0.5, 0.5, 0.5))\n",
    "# ((pid, v), is_succeeded) = w.new_voxel(Species(\"A\"), c1)\n",
    "# print(pid, v, is_succeeded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `Voxel` consists of `species`, `coordinate`, `radius` and `D`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(v.species().serial(), v.coordinate(), v.radius(), v.D())  # => (u'A', 278033, 0.0025, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, you can get a voxel and list voxels with `get_voxel` and `list_voxels_exact` similar to `get_particle` and `list_particles_exact`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(w.num_voxels_exact(Species(\"A\")))\n",
    "# print(w.list_voxels_exact(Species(\"A\")))\n",
    "# print(w.get_voxel(pid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can move and update the voxel with `update_voxel` corresponding to `update_particle`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# c2 = w.position2coordinate(Real3(0.5, 0.5, 1.0))\n",
    "# w.update_voxel(pid, Voxel(v.species(), c2, v.radius(), v.D()))\n",
    "# pid, newv = w.get_voxel(pid)\n",
    "# print(c2)  # => 278058\n",
    "# print(newv.species().serial(), newv.coordinate(), newv.radius(), newv.D())  # => (u'A', 278058, 0.0025, 1.0)\n",
    "# print(w.num_voxels_exact(Species(\"A\")))  # => 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, `remove_voxel` remove a voxel as `remove_particle` does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(w.has_voxel(pid))  # => True\n",
    "# w.remove_voxel(pid)\n",
    "# print(w.has_voxel(pid))  # => False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 Structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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": [
    "By using a `Shape` object, you can confine initial positions of molecules to a part of `World`. In the case below, 60 molecules are positioned inside the given `Sphere`. Diffusion of the molecules placed here is **NOT** restricted in the `Shape`. This `Shape` is only for the initialization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp1 = Species(\"A\", 0.0025, 1)\n",
    "sphere = Sphere(Real3(0.5, 0.5, 0.5), 0.3)\n",
    "w1.add_molecules(sp1, 60, sphere)\n",
    "w2.add_molecules(sp1, 60, sphere)\n",
    "w3.add_molecules(sp1, 60, sphere)\n",
    "w4.add_molecules(sp1, 60, sphere)\n",
    "w5.add_molecules(sp1, 60, sphere)\n",
    "w6.add_molecules(sp1, 60, sphere)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A property of `Species`, `'location'`, is available to restrict diffusion of molecules. `'location'` is not fully supported yet, but only supported in `spatiocyte` and `meso`. `add_structure` defines a new structure given as a pair of `Species` and `Shape`.\n",
    "\n",
    "NOTICE: To use `add_structure` with `spatiocyte`, you should define a model to describe the attributes of your `Species` and bind it to an instance of `spatiocyte.World`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The below codes defines a model and bind it to w3(spatiocyte world).\n",
    "# Here, the model contains a species 'B' for the following context.\n",
    "from ecell4 import species_attributes, get_model\n",
    "with species_attributes():\n",
    "    M | {'dimension': 2}\n",
    "    B | {'radius': 0.0025, 'D': 0.1, 'location': 'M'}\n",
    "model = get_model()\n",
    "w3.bind_to(model)\n",
    "\n",
    "membrane = SphericalSurface(Real3(0.5, 0.5, 0.5), 0.4)  # This is equivalent to call `Sphere(Real3(0.5, 0.5, 0.5), 0.4).surface()`\n",
    "w3.add_structure(Species(\"M\"), membrane)\n",
    "w5.add_structure(Species(\"M\"), membrane)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining a structure, you can bind molecules to the structure as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp2 = Species(\"B\", 0.0025, 0.1, \"M\")  # `'location'` is the fourth argument\n",
    "w3.add_molecules(sp2, 60)\n",
    "w5.add_molecules(sp2, 60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The molecules bound to a `Species` named `B` diffuse on a structure named `M`, which has a shape of `SphericalSurface` (a hollow sphere). In `spatiocyte`, a structure is represented as a set of particles with `Species` `M` occupying a voxel. It means that molecules not belonging to the structure is not able to overlap the voxel and it causes a collision. On the other hand, in `meso`, a structure means a list of subvolumes. Thus, a structure doesn't avoid an incursion of other particles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.5. Random Number Generator\n",
    "\n",
    "A random number generator is also a part of `World`. All `World`s except `ODEWorld` store a random number generator, and updates it when the simulation needs a random value. On E-Cell4, only one class `GSLRandomNumberGenerator` is implemented as a random number generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6, 1, 2, 6, 2, 3, 6, 5, 4, 5, 5, 4, 2, 5, 4, 2, 3, 3, 2, 2]\n"
     ]
    }
   ],
   "source": [
    "rng1 = GSLRandomNumberGenerator()\n",
    "print([rng1.uniform_int(1, 6) for _ in range(20)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With no argument, the random number generator is always initialized with a seed, `0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6, 1, 2, 6, 2, 3, 6, 5, 4, 5, 5, 4, 2, 5, 4, 2, 3, 3, 2, 2]\n"
     ]
    }
   ],
   "source": [
    "rng2 = GSLRandomNumberGenerator()\n",
    "print([rng2.uniform_int(1, 6) for _ in range(20)])  # => same as above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can initialize the seed with an integer as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6, 5, 2, 4, 1, 1, 3, 5, 2, 6, 4, 1, 2, 5, 2, 5, 1, 2, 2, 6]\n"
     ]
    }
   ],
   "source": [
    "rng2 = GSLRandomNumberGenerator()\n",
    "rng2.seed(15)\n",
    "print([rng2.uniform_int(1, 6) for _ in range(20)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you call the `seed` function with no input, the seed is drawn from the current time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3, 3, 1, 2, 2, 3, 4, 6, 3, 6, 4, 6, 5, 5, 3, 4, 1, 1, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "rng2 = GSLRandomNumberGenerator()\n",
    "rng2.seed()\n",
    "print([rng2.uniform_int(1, 6) for _ in range(20)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`GSLRandomNumberGenerator` provides several ways to get a random number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.03033520421013236\n",
      "33\n",
      "0.8935555455208181\n"
     ]
    }
   ],
   "source": [
    "print(rng1.uniform(0.0, 1.0))\n",
    "print(rng1.uniform_int(0, 100))\n",
    "print(rng1.gaussian(1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`World` accepts a random number generator at the construction. As a default, `GSLRandomNumberGenerator()` is used. Thus, when you don't give a generator, behavior of the simulation is always same (determinisitc)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = GSLRandomNumberGenerator()\n",
    "rng.seed()\n",
    "w1 = gillespie.World(Real3(1, 1, 1), rng=rng)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access the `GSLRandomNumberGenerator` in a `World` through `rng` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4002890670672059\n"
     ]
    }
   ],
   "source": [
    "print(w1.rng().uniform(0.0, 1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`rng()` returns a shared pointer to the `GSLRandomNumberGenerator`. Thus, in the example above, `rng` and `w1.rng()` point exactly the same thing."
   ]
  }
 ],
 "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
}