(4) Coagulation-Fragmentation Model - Derivation, Code Generation, and Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Coagulation-Fragmentation:\n",
    "## Full derivation with Moment Closure and Code Generation\n",
    "$$\n",
    "\\def\\n{\\mathbf{n}}\n",
    "\\def\\x{\\mathbf{x}}\n",
    "\\def\\N{\\mathbb{\\mathbb{N}}}\n",
    "\\def\\X{\\mathbb{X}}\n",
    "\\def\\NX{\\mathbb{\\N_0^\\X}}\n",
    "\\def\\C{\\mathcal{C}}\n",
    "\\def\\Jc{\\mathcal{J}_c}\n",
    "\\def\\DM{\\Delta M_{c,j}}\n",
    "\\newcommand\\diff{\\mathop{}\\!\\mathrm{d}}\n",
    "\\def\\Xc{\\mathbf{X}_c}\n",
    "\\def\\Yc{\\mathbf{Y}_c}\n",
    "\\newcommand{\\muset}[1]{\\dot{\\{}#1\\dot{\\}}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we illustrate the full treatment of another system, from model definition to code generation and simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove matplotlib deprecation warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=UserWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize sympy printing (for latex output)\n",
    "from sympy import init_printing, Symbol\n",
    "init_printing()\n",
    "\n",
    "# import functions and classes for compartment models\n",
    "from compartor import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\begin{align} \\emptyset&\\overset{h_{I}}{\\longrightarrow}\\left[y\\right] && h_{I} = \\pi_{Poiss}(y; \\lambda) k_{I}\\\\\\left[x\\right]&\\overset{h_{E}}{\\longrightarrow}\\emptyset && h_{E} = k_{E} n{\\left(x \\right)}\\\\\\left[x\\right] + \\left[y\\right]&\\overset{h_{C}}{\\longrightarrow}\\left[x + y\\right] && h_{C} = \\frac{k_{C} \\left(n{\\left(y \\right)} - \\delta_{x y}\\right) n{\\left(x \\right)}}{\\delta_{x y} + 1}\\\\\\left[x\\right]&\\overset{h_{F}}{\\longrightarrow}\\left[x - y\\right] + \\left[y\\right] && h_{F} = \\pi_{F(y|x)} k_{F} n{\\left(x \\right)} {x}_{0} \\end{align}$"
      ],
      "text/plain": [
       "                     EmptySet ---> [y],  h_I = \\pi_{Poiss}(y; \\lambda)*k_I    \n",
       "                               [x] ---> EmptySet,  h_E = k_E*n(x)             \n",
       "[x] + [y] ---> [x + y],  h_C = k_C*(n(y) - KroneckerDelta(x, y))*n(x)/(Kroneck\n",
       "                    [x] ---> [x - y] + [y],  h_F = \\pi_F(y|x)*k_F*n(x)*x[0]   \n",
       "\n",
       "                  \n",
       "                  \n",
       "erDelta(x, y) + 1)\n",
       "                  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = Content('x')\n",
    "y = Content('y')\n",
    "\n",
    "# define intake distribution\n",
    "pi_I = OutcomeDistribution.Poisson(Symbol(\"\\pi_{Poiss}(y; \\lambda)\"), y[0], Symbol(\"\\lambda\"))\n",
    "# define fragmentation distribution\n",
    "pi_F = OutcomeDistribution.Uniform(Symbol(\"\\pi_F(y|x)\"), y[0], 0, x[0])\n",
    "\n",
    "Intake        = TransitionClass( {}        -to> [y],         'k_I',         pi=pi_I, name='I')\n",
    "Exit          = TransitionClass( [x]       -to> {},          'k_E',                  name='E')\n",
    "Coagulation   = TransitionClass( [x] + [y] -to> [x+y],       'k_C',                  name='C')\n",
    "Fragmentation = TransitionClass( [x]       -to> [y] + [x-y], 'k_F', g=x[0], pi=pi_F, name='F')\n",
    "\n",
    "transitions = [Intake, Exit, Coagulation, Fragmentation]\n",
    "display_transition_classes(transitions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deriving the moment equations\n",
    "\n",
    "Next we specify the moments we want to derive, for example, $N$, $N^2$, $M^1$, and $(M^1)^2$.\n",
    "\n",
    "To denote a particular moment, use the `Moment` class. It takes as constructor arguments the integer exponents in each species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ M^{0}, \\  {\\left(M^{0}\\right)}^{\\mathtt{\\text{2}}}, \\  M^{1}, \\  {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right]$"
      ],
      "text/plain": [
       "⎡                 2                      2   ⎤\n",
       "⎣Moment(0), Moment (0), Moment(1), Moment (1)⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "moments = [\n",
    "    Moment(0),\n",
    "    Moment(0)**2,\n",
    "    Moment(1),\n",
    "    Moment(1)**2\n",
    "]\n",
    "display(moments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With the automated procedure\n",
    "\n",
    "The fastest way to generate a closed system of moment equations for the desired set of moment is to call the function `automated_moment_equations()` as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "Computed moment equations for desired moments $\\left< N\\right> $, $\\left< N^{2}\\right> $, $\\left< M^{1}\\right> $, and $\\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right> $."
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "Equations were iteratively added for $\\left< M^{2}\\right> $ and $\\left< N M^{1}\\right> $."
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "Gamma closures were substituted for $\\left< M^{3}\\right> $, $\\left< N^{2} M^{1}\\right> $, and $\\left< N^{3}\\right> $."
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N\\right>  = \\frac{k_{C} \\left< N\\right> }{2} - \\frac{k_{C} \\left< N^{2}\\right> }{2} - k_{E} \\left< N\\right>  + k_{F} \\left< M^{1}\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                                                         ⎛ 2⎞                 \n",
       "d                    k_C⋅Expectation(N)   k_C⋅Expectation⎝N ⎠                 \n",
       "──(Expectation(N)) = ────────────────── - ─────────────────── - k_E⋅Expectatio\n",
       "dt                           2                     2                          \n",
       "\n",
       "                                       \n",
       "                                       \n",
       "n(N) + k_F⋅Expectation(Moment(1)) + k_I\n",
       "                                       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N^{2}\\right>  = - k_{C} \\left(- \\left< N\\right>  \\left< N^{2}\\right>  + \\frac{2 {\\left< N^{2}\\right> }^{\\mathtt{\\text{2}}}}{\\left< N\\right> }\\right) - \\frac{k_{C} \\left< N\\right> }{2} + \\frac{3 k_{C} \\left< N^{2}\\right> }{2} + k_{E} \\left< N\\right>  - 2 k_{E} \\left< N^{2}\\right>  + 2 k_{F} \\left< N M^{1}\\right>  + k_{F} \\left< M^{1}\\right>  + 2 k_{I} \\left< N\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                            ⎛                                                2\n",
       "d ⎛           ⎛ 2⎞⎞         ⎜                            ⎛ 2⎞   2⋅Expectation \n",
       "──⎝Expectation⎝N ⎠⎠ = - k_C⋅⎜- Expectation(N)⋅Expectation⎝N ⎠ + ──────────────\n",
       "dt                          ⎝                                     Expectation(\n",
       "\n",
       "⎛ 2⎞⎞                                         ⎛ 2⎞                            \n",
       "⎝N ⎠⎟   k_C⋅Expectation(N)   3⋅k_C⋅Expectation⎝N ⎠                            \n",
       "────⎟ - ────────────────── + ───────────────────── + k_E⋅Expectation(N) - 2⋅k_\n",
       "N)  ⎠           2                      2                                      \n",
       "\n",
       "                                                                              \n",
       "             ⎛ 2⎞                                                             \n",
       "E⋅Expectation⎝N ⎠ + 2⋅k_F⋅Expectation(N⋅Moment(1)) + k_F⋅Expectation(Moment(1)\n",
       "                                                                              \n",
       "\n",
       "                              \n",
       "                              \n",
       ") + 2⋅k_I⋅Expectation(N) + k_I\n",
       "                              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< M^{1}\\right>  = \\lambda k_{I} - k_{E} \\left< M^{1}\\right> $"
      ],
      "text/plain": [
       "d                                                                    \n",
       "──(Expectation(Moment(1))) = \\lambda⋅k_I - k_E⋅Expectation(Moment(1))\n",
       "dt                                                                   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  = \\lambda^{2} k_{I} + 2 \\lambda k_{I} \\left< M^{1}\\right>  + \\lambda k_{I} - 2 k_{E} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  + k_{E} \\left< M^{2}\\right> $"
      ],
      "text/plain": [
       "d ⎛           ⎛      2   ⎞⎞          2                                        \n",
       "──⎝Expectation⎝Moment (1)⎠⎠ = \\lambda ⋅k_I + 2⋅\\lambda⋅k_I⋅Expectation(Moment(\n",
       "dt                                                                            \n",
       "\n",
       "                                     ⎛      2   ⎞                             \n",
       "1)) + \\lambda⋅k_I - 2⋅k_E⋅Expectation⎝Moment (1)⎠ + k_E⋅Expectation(Moment(2))\n",
       "                                                                              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< M^{2}\\right>  = \\lambda^{2} k_{I} + \\lambda k_{I} + k_{C} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  - k_{C} \\left< M^{2}\\right>  - k_{E} \\left< M^{2}\\right>  - \\frac{k_{F} \\left(\\frac{2 {\\left< M^{2}\\right> }^{\\mathtt{\\text{2}}}}{\\left< M^{1}\\right> } - \\frac{\\left< M^{1}\\right>  \\left< M^{2}\\right> }{\\left< N\\right> }\\right)}{3} + \\frac{k_{F} \\left< M^{2}\\right> }{3}$"
      ],
      "text/plain": [
       "                                                                              \n",
       "                                                                              \n",
       "                                                                              \n",
       "d                                   2                                    ⎛    \n",
       "──(Expectation(Moment(2))) = \\lambda ⋅k_I + \\lambda⋅k_I + k_C⋅Expectation⎝Mome\n",
       "dt                                                                            \n",
       "\n",
       "                                                                        ⎛     \n",
       "                                                                        ⎜2⋅Exp\n",
       "                                                                    k_F⋅⎜─────\n",
       "  2   ⎞                                                                 ⎝  Exp\n",
       "nt (1)⎠ - k_C⋅Expectation(Moment(2)) - k_E⋅Expectation(Moment(2)) - ──────────\n",
       "                                                                              \n",
       "\n",
       "        2                                                           ⎞         \n",
       "ectation (Moment(2))   Expectation(Moment(1))⋅Expectation(Moment(2))⎟         \n",
       "──────────────────── - ─────────────────────────────────────────────⎟         \n",
       "ectation(Moment(1))                    Expectation(N)               ⎠   k_F⋅Ex\n",
       "───────────────────────────────────────────────────────────────────── + ──────\n",
       "                             3                                                \n",
       "\n",
       "                    \n",
       "                    \n",
       "                    \n",
       "pectation(Moment(2))\n",
       "────────────────────\n",
       "      3             "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N M^{1}\\right>  = \\lambda k_{I} \\left< N\\right>  + \\lambda k_{I} - \\frac{k_{C} \\left(- \\left< N^{2}\\right>  \\left< M^{1}\\right>  + \\frac{2 \\left< N^{2}\\right>  \\left< N M^{1}\\right> }{\\left< N\\right> }\\right)}{2} + \\frac{k_{C} \\left< N M^{1}\\right> }{2} - 2 k_{E} \\left< N M^{1}\\right>  + k_{E} \\left< M^{1}\\right>  + k_{F} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  + k_{I} \\left< M^{1}\\right> $"
      ],
      "text/plain": [
       "                                                                              \n",
       "                                                                              \n",
       "                                                                          k_C⋅\n",
       "d                                                                             \n",
       "──(Expectation(N⋅Moment(1))) = \\lambda⋅k_I⋅Expectation(N) + \\lambda⋅k_I - ────\n",
       "dt                                                                            \n",
       "\n",
       "⎛                                                        ⎛ 2⎞                 \n",
       "⎜             ⎛ 2⎞                          2⋅Expectation⎝N ⎠⋅Expectation(N⋅Mo\n",
       "⎜- Expectation⎝N ⎠⋅Expectation(Moment(1)) + ──────────────────────────────────\n",
       "⎝                                                         Expectation(N)      \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                         2                                    \n",
       "\n",
       "        ⎞                                                                     \n",
       "ment(1))⎟                                                                     \n",
       "────────⎟                                                                     \n",
       "        ⎠   k_C⋅Expectation(N⋅Moment(1))                                      \n",
       "───────── + ──────────────────────────── - 2⋅k_E⋅Expectation(N⋅Moment(1)) + k_\n",
       "                         2                                                    \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "                                                                              \n",
       "                                          ⎛      2   ⎞                        \n",
       "E⋅Expectation(Moment(1)) + k_F⋅Expectation⎝Moment (1)⎠ + k_I⋅Expectation(Momen\n",
       "                                                                              \n",
       "\n",
       "     \n",
       "     \n",
       "     \n",
       "     \n",
       "t(1))\n",
       "     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "equations = automated_moment_equations(1, transitions, [Moment(0),Moment(0)**2,Moment(1),Moment(1)**2])\n",
    "display_moment_equations(equations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step-by-step\n",
    "\n",
    "Alternatively, one can opt for a step-by-step derivation to manually control the construction of the ODE system and the application of moment closure schemes.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call the function `compute_moment_equations()`, giving the reaction network and the desired moment expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "equations = compute_moment_equations(transitions, moments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To display the derived equations, use the `display_moment_equations()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N\\right>  = \\frac{k_{C} \\left< N\\right> }{2} - \\frac{k_{C} \\left< N^{2}\\right> }{2} - k_{E} \\left< N\\right>  + k_{F} \\left< M^{1}\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                                                         ⎛ 2⎞                 \n",
       "d                    k_C⋅Expectation(N)   k_C⋅Expectation⎝N ⎠                 \n",
       "──(Expectation(N)) = ────────────────── - ─────────────────── - k_E⋅Expectatio\n",
       "dt                           2                     2                          \n",
       "\n",
       "                                       \n",
       "                                       \n",
       "n(N) + k_F⋅Expectation(Moment(1)) + k_I\n",
       "                                       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N^{2}\\right>  = - \\frac{k_{C} \\left< N\\right> }{2} + \\frac{3 k_{C} \\left< N^{2}\\right> }{2} - k_{C} \\left< N^{3}\\right>  + k_{E} \\left< N\\right>  - 2 k_{E} \\left< N^{2}\\right>  + 2 k_{F} \\left< N M^{1}\\right>  + k_{F} \\left< M^{1}\\right>  + 2 k_{I} \\left< N\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                                                              ⎛ 2⎞            \n",
       "d ⎛           ⎛ 2⎞⎞     k_C⋅Expectation(N)   3⋅k_C⋅Expectation⎝N ⎠            \n",
       "──⎝Expectation⎝N ⎠⎠ = - ────────────────── + ───────────────────── - k_C⋅Expec\n",
       "dt                              2                      2                      \n",
       "\n",
       "                                                                              \n",
       "      ⎛ 3⎞                                         ⎛ 2⎞                       \n",
       "tation⎝N ⎠ + k_E⋅Expectation(N) - 2⋅k_E⋅Expectation⎝N ⎠ + 2⋅k_F⋅Expectation(N⋅\n",
       "                                                                              \n",
       "\n",
       "                                                                    \n",
       "                                                                    \n",
       "Moment(1)) + k_F⋅Expectation(Moment(1)) + 2⋅k_I⋅Expectation(N) + k_I\n",
       "                                                                    "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< M^{1}\\right>  = \\lambda k_{I} - k_{E} \\left< M^{1}\\right> $"
      ],
      "text/plain": [
       "d                                                                    \n",
       "──(Expectation(Moment(1))) = \\lambda⋅k_I - k_E⋅Expectation(Moment(1))\n",
       "dt                                                                   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  = \\lambda^{2} k_{I} + 2 \\lambda k_{I} \\left< M^{1}\\right>  + \\lambda k_{I} - 2 k_{E} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  + k_{E} \\left< M^{2}\\right> $"
      ],
      "text/plain": [
       "d ⎛           ⎛      2   ⎞⎞          2                                        \n",
       "──⎝Expectation⎝Moment (1)⎠⎠ = \\lambda ⋅k_I + 2⋅\\lambda⋅k_I⋅Expectation(Moment(\n",
       "dt                                                                            \n",
       "\n",
       "                                     ⎛      2   ⎞                             \n",
       "1)) + \\lambda⋅k_I - 2⋅k_E⋅Expectation⎝Moment (1)⎠ + k_E⋅Expectation(Moment(2))\n",
       "                                                                              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "**The system is not closed!** Moment equations are missing for:"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{M^{0} M^{1}, {\\left(M^{0}\\right)}^{\\mathtt{\\text{3}}}, M^{2}\\right\\}$"
      ],
      "text/plain": [
       "⎧                           3              ⎫\n",
       "⎨Moment(0)⋅Moment(1), Moment (0), Moment(2)⎬\n",
       "⎩                                          ⎭"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_moment_equations(equations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `display_moment_equations()` function also notifies us that the system is not closed, and which moments are missing. (This output can be suppressed by passing the optional second argument `display_moment_equations(equations, print_missing=False)`).\n",
    "\n",
    "These \"missing\" moments are moment expressions whose expectations occur on the right-hand side of equations but which are not governed by the system. They can also be computed using the `get_missing_moments()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{M^{0} M^{1}, {\\left(M^{0}\\right)}^{\\mathtt{\\text{3}}}, M^{2}\\right\\}$"
      ],
      "text/plain": [
       "⎧                           3              ⎫\n",
       "⎨Moment(0)⋅Moment(1), Moment (0), Moment(2)⎬\n",
       "⎩                                          ⎭"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "missing = get_missing_moments(equations)\n",
    "missing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We add the moments below order 3 to the list of moments to derive and iterate the derivation. (We apply gamma closure for 3rd order moments below.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ M^{0}, \\  {\\left(M^{0}\\right)}^{\\mathtt{\\text{2}}}, \\  M^{1}, \\  {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}, \\  M^{0} M^{1}, \\  M^{2}\\right]$"
      ],
      "text/plain": [
       "⎡                 2                      2                                   ⎤\n",
       "⎣Moment(0), Moment (0), Moment(1), Moment (1), Moment(0)⋅Moment(1), Moment(2)⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "moments = [\n",
    "    Moment(0),\n",
    "    Moment(0)**2,\n",
    "    Moment(1),\n",
    "    Moment(1)**2,\n",
    "    Moment(0)*Moment(1),\n",
    "    Moment(2)\n",
    "]\n",
    "display(moments)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N\\right>  = \\frac{k_{C} \\left< N\\right> }{2} - \\frac{k_{C} \\left< N^{2}\\right> }{2} - k_{E} \\left< N\\right>  + k_{F} \\left< M^{1}\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                                                         ⎛ 2⎞                 \n",
       "d                    k_C⋅Expectation(N)   k_C⋅Expectation⎝N ⎠                 \n",
       "──(Expectation(N)) = ────────────────── - ─────────────────── - k_E⋅Expectatio\n",
       "dt                           2                     2                          \n",
       "\n",
       "                                       \n",
       "                                       \n",
       "n(N) + k_F⋅Expectation(Moment(1)) + k_I\n",
       "                                       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N^{2}\\right>  = - \\frac{k_{C} \\left< N\\right> }{2} + \\frac{3 k_{C} \\left< N^{2}\\right> }{2} - k_{C} \\left< N^{3}\\right>  + k_{E} \\left< N\\right>  - 2 k_{E} \\left< N^{2}\\right>  + 2 k_{F} \\left< N M^{1}\\right>  + k_{F} \\left< M^{1}\\right>  + 2 k_{I} \\left< N\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                                                              ⎛ 2⎞            \n",
       "d ⎛           ⎛ 2⎞⎞     k_C⋅Expectation(N)   3⋅k_C⋅Expectation⎝N ⎠            \n",
       "──⎝Expectation⎝N ⎠⎠ = - ────────────────── + ───────────────────── - k_C⋅Expec\n",
       "dt                              2                      2                      \n",
       "\n",
       "                                                                              \n",
       "      ⎛ 3⎞                                         ⎛ 2⎞                       \n",
       "tation⎝N ⎠ + k_E⋅Expectation(N) - 2⋅k_E⋅Expectation⎝N ⎠ + 2⋅k_F⋅Expectation(N⋅\n",
       "                                                                              \n",
       "\n",
       "                                                                    \n",
       "                                                                    \n",
       "Moment(1)) + k_F⋅Expectation(Moment(1)) + 2⋅k_I⋅Expectation(N) + k_I\n",
       "                                                                    "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< M^{1}\\right>  = \\lambda k_{I} - k_{E} \\left< M^{1}\\right> $"
      ],
      "text/plain": [
       "d                                                                    \n",
       "──(Expectation(Moment(1))) = \\lambda⋅k_I - k_E⋅Expectation(Moment(1))\n",
       "dt                                                                   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  = \\lambda^{2} k_{I} + 2 \\lambda k_{I} \\left< M^{1}\\right>  + \\lambda k_{I} - 2 k_{E} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  + k_{E} \\left< M^{2}\\right> $"
      ],
      "text/plain": [
       "d ⎛           ⎛      2   ⎞⎞          2                                        \n",
       "──⎝Expectation⎝Moment (1)⎠⎠ = \\lambda ⋅k_I + 2⋅\\lambda⋅k_I⋅Expectation(Moment(\n",
       "dt                                                                            \n",
       "\n",
       "                                     ⎛      2   ⎞                             \n",
       "1)) + \\lambda⋅k_I - 2⋅k_E⋅Expectation⎝Moment (1)⎠ + k_E⋅Expectation(Moment(2))\n",
       "                                                                              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N M^{1}\\right>  = \\lambda k_{I} \\left< N\\right>  + \\lambda k_{I} + \\frac{k_{C} \\left< N M^{1}\\right> }{2} - \\frac{k_{C} \\left< N^{2} M^{1}\\right> }{2} - 2 k_{E} \\left< N M^{1}\\right>  + k_{E} \\left< M^{1}\\right>  + k_{F} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  + k_{I} \\left< M^{1}\\right> $"
      ],
      "text/plain": [
       "                                                                              \n",
       "d                                                                         k_C⋅\n",
       "──(Expectation(N⋅Moment(1))) = \\lambda⋅k_I⋅Expectation(N) + \\lambda⋅k_I + ────\n",
       "dt                                                                            \n",
       "\n",
       "                                          ⎛ 2          ⎞                      \n",
       "Expectation(N⋅Moment(1))   k_C⋅Expectation⎝N ⋅Moment(1)⎠                      \n",
       "──────────────────────── - ───────────────────────────── - 2⋅k_E⋅Expectation(N\n",
       "         2                               2                                    \n",
       "\n",
       "                                                                              \n",
       "                                                          ⎛      2   ⎞        \n",
       "⋅Moment(1)) + k_E⋅Expectation(Moment(1)) + k_F⋅Expectation⎝Moment (1)⎠ + k_I⋅E\n",
       "                                                                              \n",
       "\n",
       "                     \n",
       "                     \n",
       "xpectation(Moment(1))\n",
       "                     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< M^{2}\\right>  = \\lambda^{2} k_{I} + \\lambda k_{I} + k_{C} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  - k_{C} \\left< M^{2}\\right>  - k_{E} \\left< M^{2}\\right>  + \\frac{k_{F} \\left< M^{2}\\right> }{3} - \\frac{k_{F} \\left< M^{3}\\right> }{3}$"
      ],
      "text/plain": [
       "d                                   2                                    ⎛    \n",
       "──(Expectation(Moment(2))) = \\lambda ⋅k_I + \\lambda⋅k_I + k_C⋅Expectation⎝Mome\n",
       "dt                                                                            \n",
       "\n",
       "  2   ⎞                                                             k_F⋅Expect\n",
       "nt (1)⎠ - k_C⋅Expectation(Moment(2)) - k_E⋅Expectation(Moment(2)) + ──────────\n",
       "                                                                              \n",
       "\n",
       "ation(Moment(2))   k_F⋅Expectation(Moment(3))\n",
       "──────────────── - ──────────────────────────\n",
       "  3                            3             "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "**The system is not closed!** Moment equations are missing for:"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{{\\left(M^{0}\\right)}^{\\mathtt{\\text{2}}} M^{1}, {\\left(M^{0}\\right)}^{\\mathtt{\\text{3}}}, M^{3}\\right\\}$"
      ],
      "text/plain": [
       "⎧      2                     3              ⎫\n",
       "⎨Moment (0)⋅Moment(1), Moment (0), Moment(3)⎬\n",
       "⎩                                           ⎭"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "equations = compute_moment_equations(transitions, moments)\n",
    "display_moment_equations(equations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gamma closure\n",
    "Next we compute the Gamma closures for the missing moments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left< M^{3}\\right>  = \\frac{2 {\\left< M^{2}\\right> }^{\\mathtt{\\text{2}}}}{\\left< M^{1}\\right> } - \\frac{\\left< M^{1}\\right>  \\left< M^{2}\\right> }{\\left< N\\right> }$"
      ],
      "text/plain": [
       "                                      2                                       \n",
       "                         2⋅Expectation (Moment(2))   Expectation(Moment(1))⋅Ex\n",
       "Expectation(Moment(3)) = ───────────────────────── - ─────────────────────────\n",
       "                           Expectation(Moment(1))                    Expectati\n",
       "\n",
       "                    \n",
       "pectation(Moment(2))\n",
       "────────────────────\n",
       "on(N)               "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left< N^{2} M^{1}\\right>  = - \\left< N^{2}\\right>  \\left< M^{1}\\right>  + \\frac{2 \\left< N^{2}\\right>  \\left< N M^{1}\\right> }{\\left< N\\right> }$"
      ],
      "text/plain": [
       "                                                                              \n",
       "           ⎛ 2          ⎞                ⎛ 2⎞                          2⋅Expec\n",
       "Expectation⎝N ⋅Moment(1)⎠ = - Expectation⎝N ⎠⋅Expectation(Moment(1)) + ───────\n",
       "                                                                              \n",
       "\n",
       "      ⎛ 2⎞                         \n",
       "tation⎝N ⎠⋅Expectation(N⋅Moment(1))\n",
       "───────────────────────────────────\n",
       "       Expectation(N)              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left< N^{3}\\right>  = - \\left< N\\right>  \\left< N^{2}\\right>  + \\frac{2 {\\left< N^{2}\\right> }^{\\mathtt{\\text{2}}}}{\\left< N\\right> }$"
      ],
      "text/plain": [
       "                                                                  2⎛ 2⎞\n",
       "           ⎛ 3⎞                               ⎛ 2⎞   2⋅Expectation ⎝N ⎠\n",
       "Expectation⎝N ⎠ = - Expectation(N)⋅Expectation⎝N ⎠ + ──────────────────\n",
       "                                                       Expectation(N)  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "missing = get_missing_moments(equations)\n",
    "closures = gamma_closures(missing)\n",
    "display_closures(closures)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Currently, only Gamma closure for 3rd order moments is implemented.\n",
    "If the closure for an expression can not be computed, the function raises an error.\n",
    "Alternatively, a list can be supplied as second argument `gamma_closures(expressions, could_not_close)`, and the expressions that could not be closed will be added there (instead or raising an error).\n",
    "\n",
    "In this example, all of the `missing` expression could be closed, so..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABAAAAAVCAYAAABPPm7SAAAAd0lEQVR4nO2UwQmAMAxFn8URvDuKHcVNjJvVLTw6gOde60UxFtNLr/0QCP+HRykknYhQI6d6DyRVu8qGLEtP0P9ANyAAp/IisN79DIwlQAAk86LyvAY4KtUADdAANmChfA8mPazX+eDdebDvwUc5QIyXRiur/oMLtccckPBLdKYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[ \\right]$"
      ],
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "could_not_close=[]\n",
    "closures = gamma_closures(missing, could_not_close)\n",
    "could_not_close"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We substitute the closures into the moment equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N\\right>  = \\frac{k_{C} \\left< N\\right> }{2} - \\frac{k_{C} \\left< N^{2}\\right> }{2} - k_{E} \\left< N\\right>  + k_{F} \\left< M^{1}\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                                                         ⎛ 2⎞                 \n",
       "d                    k_C⋅Expectation(N)   k_C⋅Expectation⎝N ⎠                 \n",
       "──(Expectation(N)) = ────────────────── - ─────────────────── - k_E⋅Expectatio\n",
       "dt                           2                     2                          \n",
       "\n",
       "                                       \n",
       "                                       \n",
       "n(N) + k_F⋅Expectation(Moment(1)) + k_I\n",
       "                                       "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N^{2}\\right>  = - k_{C} \\left(- \\left< N\\right>  \\left< N^{2}\\right>  + \\frac{2 {\\left< N^{2}\\right> }^{\\mathtt{\\text{2}}}}{\\left< N\\right> }\\right) - \\frac{k_{C} \\left< N\\right> }{2} + \\frac{3 k_{C} \\left< N^{2}\\right> }{2} + k_{E} \\left< N\\right>  - 2 k_{E} \\left< N^{2}\\right>  + 2 k_{F} \\left< N M^{1}\\right>  + k_{F} \\left< M^{1}\\right>  + 2 k_{I} \\left< N\\right>  + k_{I}$"
      ],
      "text/plain": [
       "                            ⎛                                                2\n",
       "d ⎛           ⎛ 2⎞⎞         ⎜                            ⎛ 2⎞   2⋅Expectation \n",
       "──⎝Expectation⎝N ⎠⎠ = - k_C⋅⎜- Expectation(N)⋅Expectation⎝N ⎠ + ──────────────\n",
       "dt                          ⎝                                     Expectation(\n",
       "\n",
       "⎛ 2⎞⎞                                         ⎛ 2⎞                            \n",
       "⎝N ⎠⎟   k_C⋅Expectation(N)   3⋅k_C⋅Expectation⎝N ⎠                            \n",
       "────⎟ - ────────────────── + ───────────────────── + k_E⋅Expectation(N) - 2⋅k_\n",
       "N)  ⎠           2                      2                                      \n",
       "\n",
       "                                                                              \n",
       "             ⎛ 2⎞                                                             \n",
       "E⋅Expectation⎝N ⎠ + 2⋅k_F⋅Expectation(N⋅Moment(1)) + k_F⋅Expectation(Moment(1)\n",
       "                                                                              \n",
       "\n",
       "                              \n",
       "                              \n",
       ") + 2⋅k_I⋅Expectation(N) + k_I\n",
       "                              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< M^{1}\\right>  = \\lambda k_{I} - k_{E} \\left< M^{1}\\right> $"
      ],
      "text/plain": [
       "d                                                                    \n",
       "──(Expectation(Moment(1))) = \\lambda⋅k_I - k_E⋅Expectation(Moment(1))\n",
       "dt                                                                   "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  = \\lambda^{2} k_{I} + 2 \\lambda k_{I} \\left< M^{1}\\right>  + \\lambda k_{I} - 2 k_{E} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  + k_{E} \\left< M^{2}\\right> $"
      ],
      "text/plain": [
       "d ⎛           ⎛      2   ⎞⎞          2                                        \n",
       "──⎝Expectation⎝Moment (1)⎠⎠ = \\lambda ⋅k_I + 2⋅\\lambda⋅k_I⋅Expectation(Moment(\n",
       "dt                                                                            \n",
       "\n",
       "                                     ⎛      2   ⎞                             \n",
       "1)) + \\lambda⋅k_I - 2⋅k_E⋅Expectation⎝Moment (1)⎠ + k_E⋅Expectation(Moment(2))\n",
       "                                                                              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< N M^{1}\\right>  = \\lambda k_{I} \\left< N\\right>  + \\lambda k_{I} - \\frac{k_{C} \\left(- \\left< N^{2}\\right>  \\left< M^{1}\\right>  + \\frac{2 \\left< N^{2}\\right>  \\left< N M^{1}\\right> }{\\left< N\\right> }\\right)}{2} + \\frac{k_{C} \\left< N M^{1}\\right> }{2} - 2 k_{E} \\left< N M^{1}\\right>  + k_{E} \\left< M^{1}\\right>  + k_{F} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  + k_{I} \\left< M^{1}\\right> $"
      ],
      "text/plain": [
       "                                                                              \n",
       "                                                                              \n",
       "                                                                          k_C⋅\n",
       "d                                                                             \n",
       "──(Expectation(N⋅Moment(1))) = \\lambda⋅k_I⋅Expectation(N) + \\lambda⋅k_I - ────\n",
       "dt                                                                            \n",
       "\n",
       "⎛                                                        ⎛ 2⎞                 \n",
       "⎜             ⎛ 2⎞                          2⋅Expectation⎝N ⎠⋅Expectation(N⋅Mo\n",
       "⎜- Expectation⎝N ⎠⋅Expectation(Moment(1)) + ──────────────────────────────────\n",
       "⎝                                                         Expectation(N)      \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                         2                                    \n",
       "\n",
       "        ⎞                                                                     \n",
       "ment(1))⎟                                                                     \n",
       "────────⎟                                                                     \n",
       "        ⎠   k_C⋅Expectation(N⋅Moment(1))                                      \n",
       "───────── + ──────────────────────────── - 2⋅k_E⋅Expectation(N⋅Moment(1)) + k_\n",
       "                         2                                                    \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "                                                                              \n",
       "                                          ⎛      2   ⎞                        \n",
       "E⋅Expectation(Moment(1)) + k_F⋅Expectation⎝Moment (1)⎠ + k_I⋅Expectation(Momen\n",
       "                                                                              \n",
       "\n",
       "     \n",
       "     \n",
       "     \n",
       "     \n",
       "t(1))\n",
       "     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d}{d t} \\left< M^{2}\\right>  = \\lambda^{2} k_{I} + \\lambda k_{I} + k_{C} \\left< {\\left(M^{1}\\right)}^{\\mathtt{\\text{2}}}\\right>  - k_{C} \\left< M^{2}\\right>  - k_{E} \\left< M^{2}\\right>  - \\frac{k_{F} \\left(\\frac{2 {\\left< M^{2}\\right> }^{\\mathtt{\\text{2}}}}{\\left< M^{1}\\right> } - \\frac{\\left< M^{1}\\right>  \\left< M^{2}\\right> }{\\left< N\\right> }\\right)}{3} + \\frac{k_{F} \\left< M^{2}\\right> }{3}$"
      ],
      "text/plain": [
       "                                                                              \n",
       "                                                                              \n",
       "                                                                              \n",
       "d                                   2                                    ⎛    \n",
       "──(Expectation(Moment(2))) = \\lambda ⋅k_I + \\lambda⋅k_I + k_C⋅Expectation⎝Mome\n",
       "dt                                                                            \n",
       "\n",
       "                                                                        ⎛     \n",
       "                                                                        ⎜2⋅Exp\n",
       "                                                                    k_F⋅⎜─────\n",
       "  2   ⎞                                                                 ⎝  Exp\n",
       "nt (1)⎠ - k_C⋅Expectation(Moment(2)) - k_E⋅Expectation(Moment(2)) - ──────────\n",
       "                                                                              \n",
       "\n",
       "        2                                                           ⎞         \n",
       "ectation (Moment(2))   Expectation(Moment(1))⋅Expectation(Moment(2))⎟         \n",
       "──────────────────── - ─────────────────────────────────────────────⎟         \n",
       "ectation(Moment(1))                    Expectation(N)               ⎠   k_F⋅Ex\n",
       "───────────────────────────────────────────────────────────────────── + ──────\n",
       "                             3                                                \n",
       "\n",
       "                    \n",
       "                    \n",
       "                    \n",
       "pectation(Moment(2))\n",
       "────────────────────\n",
       "      3             "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "equations = substitute_closures(equations, closures)\n",
    "display_moment_equations(equations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code generation\n",
    "From the closed equations, we can generate code for simulating the system.\n",
    "Currently, Python or Julia code can be generated.\n",
    "Both derive from a common `AbstractCodeGenerator` class, and it should be relatively straightforward to tweak the code generation for other languages/solvers.\n",
    "\n",
    "Here, we use the Python generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "generator = GeneratePython()\n",
    "code = generator.generate(equations, function_name=\"coagulation_fragmentation\")\n",
    "exec(code) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This creates code for two functions, `coagulation_fragmentation_ODEs` and `coagulation_fragmentation_initial`.\n",
    "The arguments to the functions are listed in the doc strings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function coagulation_fragmentation_initial in module __main__:\n",
      "\n",
      "coagulation_fragmentation_initial(initial)\n",
      "    Create inital expected moments vector\n",
      "    \n",
      "    :param initial: tuple of initial values for expectations of (Moment(0,), Moment(1,), Moment(2,))\n",
      "    :return: initial expected moments vector\n",
      "\n",
      "Help on function coagulation_fragmentation_ODEs in module __main__:\n",
      "\n",
      "coagulation_fragmentation_ODEs(M, dM, parameters)\n",
      "    Evaluate derivatives of expected moments\n",
      "    \n",
      "    Indices in M and dM vectors are \n",
      "      M[0] = Moment(0,)\n",
      "      M[1] = Moment(0)**2\n",
      "      M[2] = Moment(1,)\n",
      "      M[3] = Moment(1)**2\n",
      "      M[4] = Moment(0)*Moment(1)\n",
      "      M[5] = Moment(2,)\n",
      "    \n",
      "    :param M: expected moments\n",
      "    :param dM: result, the derivative dM/dt is stored here\n",
      "    :param parameters: tuple of values for constants (\\lambda, k_C, k_E, k_F, k_I)\n",
      "    :return: dM\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(coagulation_fragmentation_initial)\n",
    "help(coagulation_fragmentation_ODEs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving and plotting results\n",
    "We use the generated functions with `scipy.solve_ivp` to solve the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[100.          98.09964031  96.31871186 ... 100.04489169 100.04479285\n",
      " 100.04470542]\n"
     ]
    }
   ],
   "source": [
    "from scipy.integrate import solve_ivp\n",
    "import numpy as np\n",
    "\n",
    "# vector of timepoints where to solve ODE problem\n",
    "timepoints = np.linspace(0, 100, 1001)\n",
    "\n",
    "# initial condition\n",
    "# 100 compartments with 10 molecules each\n",
    "M0 = coagulation_fragmentation_initial((100, 1000, 10000))\n",
    "\n",
    "# rate constants and other parameters\n",
    "kC = 0.005\n",
    "kF = 0.005\n",
    "kI = 10.0\n",
    "kE = 0.1\n",
    "Lambda = 50\n",
    "parameters = (Lambda, kC, kE, kF, kI)\n",
    "\n",
    "# solve\n",
    "sol = solve_ivp(\n",
    "    lambda t, M: coagulation_fragmentation_ODEs(M, np.zeros(6), parameters),\n",
    "    (np.min(timepoints), np.max(timepoints)),\n",
    "    M0,\n",
    "    method='BDF',\n",
    "    t_eval=timepoints)\n",
    "\n",
    "# print solution of N\n",
    "print(sol.y[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we use `matplotlib` to plot the solution for mean and standard deviation of $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIcAAAAVCAYAAABsSf1CAAAD10lEQVR4nO3aWWhdVRQG4C91rFFSiajg0Ei00QdxQqtIG4K1Cr4oCn2xICIoipWiCA7greLw4FCcFXEgPuhDUREUBS1EKw6IPgjRViVVsRKnRKtWtNaHdU5zenJOuFdzc++h94fLunfts/+998o6e6+1drpqtZoOOijCvBL9MxhH9xzOpYPW4BTswGX5hiLnOBUrcRd+y+jHEpKiz3czDH44nsS3+DPhWYsDG1lBE7jqxUV4AG/hF7HeZ2d4fky17PQhXsRt2D/bsGfBw7cLIzxS0DaZDJLH1pJJ9uMdHIyX8ClOwzU4F2fix5K+zeRqBDfjBLHGb3BsHX2qZqc78R5W4Y5UmXeORViGJ/BHwcATqNU5SXg4meQq8faluBerhSNe0QKuRrBaOMXnGMT6OvpMqJad3hdOdLk4Mf5h+rFyKbrwfJ0TmQn9WC62tIdybbeII2ul+uKa2eRqFOuxSRwLzUC72Ok5HImzU0XeOZZhO94tGXwfXIwbxTY1hD1Knh1K5OsST8zgV2zAfji9pH+zuOYCVbTThkQWOkc3TsSoXQPRLA7FsNia1uJN8VYNFjw7kMiNJVybErmopL1ZXHOBKtrpg0QuTRVZ5zhMePeWEuKncJZYeDeOx2Pow6siaMuiJ5GTJXypfkFJe7O4mo2q2mkS28TRgl0D0t5E/lxCvCb3+xMR2GzFtSIAu6Ck7+6EKtvpJxyS/sjuHGl2sm+DhI8mcmlOn3ppj2Kk+ok6xphNrlahCnaaL5OlZp1jPJG9GsP3icxHwJ8lsuysPCaRZedjs7hahXa30zxx3IxnFSm2iAUMaAxp5PtlTp/WA5abnhUdIIoxvyvPjJrF1Sq0u50GRBnj41SRJdiBERyEo3Mdj1OcG/fhweR7vqT8hUip+nBVrm1NwjdsembUL6qQe80C19NiXZcUzL0ZqKqdmHLenUW+fIV0HS7EOaIimGKFCKZGsFnkzP04T8Qor+DuggGvFKXc+0UEP4rFIh/fiJsK+ryBhThKFHP+D1fq/H8XtNWL85MPkYHAGcLx4Adcl3yvqp2I3Wa7KLmDrtyV/d74OhlscUY/KCLuk0ylaBNiCxpOPmUVxCNwq6jr94rj6wXhyUWZ0ZjiRf8Xro/EH2dhSXs9qInqYhk2izeV6tqpR1wKvmbqRZjmHHCDuHw5WRi3qlggLpjuwfWtnUrb42qx0yzB26my6Mr+PnwlPK/KWIK/xIVTB+WYLzaEdTKOQfGV/TZxOTMktsWyUnq742WN12x2R/ThcVMx1E4UOQcRUI00bz4dtBFGlfx7Qdm/CXbQgX8BR3pji2T3U9AAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left( 50.0, \\  150.0\\right)$"
      ],
      "text/plain": [
       "(50.0, 150.0)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEZCAYAAABiu9n+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuqElEQVR4nO3deZhddZ3v+/d31zxXKjVkrhAIkBCmGAS6FbBBxRwUD40MtgiI0vb1dNvqedT29Dna97a3j/d2t1cfbBtaBpFJpDsytIA0DQKKQBIIIYQhZKwkVamk5nFP3/vHWlXsKmolNe9K5fPKs5+95vVbe6XWZ6/fWnv9zN0REREZSSzbBRARkZlLISEiIpEUEiIiEkkhISIikRQSIiISSSEhIiKRFBJyVDOz68zsuWyXQ6aG9m/2KSSOMWa208x6zawr43VTFsvztJl9Plvrn4nM7DtmdlcW1+9mdkK21i8zS262CyBZ8XF3/49sF0Ley8z0NzlFzCzX3ZPZLsfRRmcSMsjMfmxm/5rR/z0ze9ICF5hZg5l9y8wOhmckf5IxbYGZ/b2Z7TazJjP7ZzMryhh/qZm9YmYdZvaOmV1sZt8FPgjclHlGY2Ynm9kTZtZiZm+a2RUZy5lrZg+Fy3kROP4I2/QBM/udmbWZ2R4zuy4cXmFmd5pZs5ntMrO/NrNYOO46M/utmX0/nG+7mf1BOHyPmR0ws2sz1nFHuL1PmFmnmf3GzOozxv8gnK/DzDaY2Qczxn3HzB4ws7vMrAP4IvAt4MrwM9kUTve0mf1tuC1dZvZw+FncHS73JTNbmrHcw32Gd5jZj8zs38PyvmBmx4fjngkn2xSu58oRPtPrzOy5cH+3mtkOM/tYxvidZnbRsG28K+xeGp6pXB9+Jq1m9kUzO8vMXg0/7+FntmZmN5lZu5m9YWYXZoyoMLNbzWy/me0NP6OcEfbjIeA7h/u/IhHcXa9j6AXsBC6KGFcMvAVcR3DwPggsCsddACSBfwQKgPOBbuCkcPz3gYeAKqAMeBj4u3Dc+4F24MMEX0wWAieH454GPp9RhhJgD3A9wZnumWE5Vobj7wPuD6dbBewFnovYnnqgE7gayAPmAmeE4+4EHgzLujTc7hvCcdeF23o9kAP8LbAb+FG47R8Jl1saTn9H2H9eOP4HmWUCPhOuOxf4GtAIFIbjvgMkgE+Gn01ROOyuYdvyNLCNIBQrgNfDMl8ULvdO4PZRfoZ3AIfC/ZIL3A3cl7EuB044zP+h68IyfyH8fP4M2AfYSP/HMrcn/Kwd+GegMPws+4BfArUE/zcOAOcP2xdfCffhlQT/l6rC8euAm8NtrgVeBP502Lx/Hm5nUbb//o7GV9YLoNc07/DgD7gLaMt4fSFj/NlAC7ALuDpj+AXhH1xJxrD7gf8JGEFgHJ8x7lxgR9h9M/D9iPI8zdCQuBJ4dtg0NwPfDg9ICcKACcf930SHxF8B60YYngPEBw6a4bA/BZ4Ou68D3s4Yd2p4YKvLGHaIdwPnDoYeZEuBFLA4olytwOlh93eAZ4aNHzyoDvuc/kdG/z8Aj2b0fxx45UifYUZ5f5Ixbi3wRkb/aEJiW0Z/cTjPvIz/Y0cKiYXDPssrM/r/FfjLjHUNBlA47EXgGqAO6Cfj4E/wheCpjHl3Z/tv7mh/qf7z2PRJj7gm4e4vmNl2gm9l9w8b3eru3Rn9u4AFQA3BgWKDmQ2MM4KDMcBi4FejLFs9cLaZtWUMywV+Fq4nl+BbcmYZoiwG3hlheDXBt9LMeXcRfIsd0JTR3Qvg7sOHlWb0D5bJ3bvMrIXgs9ljZv8duCHsd6A8LMN75j2C4euPKs/hPsMBjRndPcO2ZTQG53f3nnC/j2UZo90WgL0eHvVDA//v6gn24/6M/3cxhn6eo/1sJYJCQoYwsy8RVJnsA74O/F3G6DlmVpIRFEuA1wiqMnqBU9x97wiL3UP0tYPhjyHeA/zG3T88QtlyCM5mFgNvZJQhyh6CKpXhDhKckdQTVNsMLGekso/W4oxylhJUu+0Lrz98HbgQ2OLuaTNrJQjRAcM/g4k+mjnyM5wm3QRfGgbMm+DyFpqZZQTFEoKqzT0EZxLVHn1BWo+5niBduJZBZnYiQf37ZwhO579uZmcMm+xvzCw/PPhdAvzC3dPAvwDfN7PacFkLzeyj4Ty3Ateb2YVmFgvHnRyOawKWZSz/EeBEM7vGzPLC11lmtsLdU8C/Ad8xs2IzWwlcS7S7gYvM7Aozyw0v9J4RLud+4LtmVhZeZP4qMJHbTtdacJE8H/i/gN+7+x6Cax5JoBnINbP/RXAmcThNwFILL6SPQ+RnOMr5h++TsXoFuCpc7xrg8gksC4Kz2r8Il/cpYAXwK3ffD/wa+AczKw//bx1vZudPcH2SQSFxbHrYhv5OYp0Ft17eBXzP3Te5+9sEd9n8zMwKwvkaCerT9xEcgL/o7gPf6L9BcGH19xbcpfMfwEkA7v4iwUXU7xNcdPwNwbd4CC7yXh7e5fJDd+8kuJh5VbieRuB7BGc3AP+NoCqikaBu/faojXT33QT17V8juM7yCnB6OPrPCb7xbgeeA+4Bbhv1J/he9xBcN2kB3kcQtACPA48RXGTeRXCR9khVIL8I3w+Z2caxFmQUn+GRfAf4aXin0RVHmngE/5PgzLEV+BuCz2YiXgCWE5wBfhe43N0PheM+C+QTnBG2Ag8A8ye4PslgQ6v6REZmZhcQXHxclOWizDhmdgfQ4O5/ne2yiEw2nUmIiEgkhYSIiERSdZOIiETSmYSIiESadb+TqK6u9qVLl2a7GCIiR40NGzYcdPeakcbNupBYunQp69evz3YxRESOGmYW+eQCVTeJiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIpGkNCTO7zcwOmNlrI4z7mpm5mVWH/WZmPzSzbWb2qpmtns6yiojI9J9J3AFcPHygmS0GPgLszhj8MWB5+LoR+PE0lE9ERDJMa0i4+zNAywijvg98HfCMYZcCd3rg90Clmc2fhmKKiEgo69ckzOxSYK+7bxo2aiGwJ6O/IRw20jJuNLP1Zra+ubl5ikoqInLsyWpImFkx8C3gf01kOe5+i7uvcfc1NTU1k1M4EREhN8vrPx44DthkZgCLgI1m9n5gL7A4Y9pF4TAREZkmWT2TcPfN7l7r7kvdfSlBldJqd28EHgI+G97ldA7Q7u77s1leEZFjzXTfAnsv8Dxwkpk1mNkNh5n8V8B2YBvwL8D/MQ1FFBGRDNNa3eTuVx9h/NKMbge+NNVlEhGRaFm/u0lERGYuhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRJrWkDCz28zsgJm9ljHs/zWzN8zsVTNbZ2aVGeP+ysy2mdmbZvbR6SyriIhM/5nEHcDFw4Y9Aaxy99OAt4C/AjCzlcBVwCnhPP9kZjnTV1QREZnWkHD3Z4CWYcN+7e7JsPf3wKKw+1LgPnfvd/cdwDbg/dNWWBERmXHXJD4HPBp2LwT2ZIxrCIeJiMg0mTEhYWb/A0gCd49j3hvNbL2ZrW9ubp78womIHKNmREiY2XXAJcCfuLuHg/cCizMmWxQOew93v8Xd17j7mpqamiktq4jIsSTrIWFmFwNfBz7h7j0Zox4CrjKzAjM7DlgOvJiNMoqIHKtyp3NlZnYvcAFQbWYNwLcJ7mYqAJ4wM4Dfu/sX3X2Lmd0PvE5QDfUld09NZ3lFRI519m7tzuywZs0aX79+fbaLISJy1DCzDe6+ZqRxWa9uEhGRmUshISIikRQSIiISSSEhIiKRFBIiIhJJISEiIpEUEiIiEkkhISIikRQSIiISSSEhIiKRFBIiIhJJISEiIpEUEiIiEkkhISIikRQSIiISSSEhIiKRFBIiIhJpWpsvFTmWpD1NKp0i5anB7rSnB/uP9HL3d7vxwWHDuweHpdOkCecL33EGuwfnI2iN0tPh/AP9Gcsb0p8xDBicHn+3O3MZA/2Z0w3OO9L4EYw4LmLygXKOZjkjbscI46LGH7bMI8w/0vSRrYGOtZFQG9qbn5PPFadcQWFu4RgXdHgKCTnmJNNJ4qk4yXSSRCpBIp0gmU5GvvqSfcRT8eCVjNOf6n+3Px2nP9lPIhUsI3NZ7g4G5oaZDR4wzGzwgGDhX/qQg4kx5ACcyXHMDWzoAdcI1pG5TOzdbss4ogxM9+7q7LDDhxvtdFHTR043ScuZyDpGu76xLGc0y5sMjZ2NfOKkTygkRAD6k/30p/rpS/YN6Y6n4nTFu+iOd9OT6AleyR56E73BK9lLKp3CzN5zMB3yTTHjIB2LxYgRIyeWQ8xixCxGjmV0x3LIy8mjILcAwwaHT/VBQSRTe1/7lCxXISEzSiKVoCveRWe8k+54N53xTlp6W2jrbaOtr42OeAdd8S7SngaGfisfqBbJieWQG8slN5ZLjr3bXVZQxpyiOcRMl+JERkshIdOuP9lPW18b7f3ttPa2cqD7AE3dTTR3N9Od6B5aBeOQn5tPfk7wKsotorygXAd6kWmikJAp4+609bVxsOcgTV1N7Gzfyb7OfbT3tQ85AyjILaAwt5CygjLmFs/NbqFFZAiFhEyavmQfjV2N7O3cyzst77CzbSfxZHywzr8kv4TivGIqKyuzXVQRGSWFhIxbX7KPfZ372NG2g60HttLY3QgEgVCaX8rcornk5eRluZQiMhEKCRk1d6epu4mdbTvZ3LSZ3e27cZyYxagsrGRx+WLd0SMyyygk5LCS6SR72vew9eBWNjVuojvRDUBFQQULyxfqArLILKeQkPdIppPsbt/N5gOb2dS4if5kP3k5eVQVVenCssgxRiEhgxq7GtnUuImX9r5Eb7KXgtwCXVcQOcYpJI5xfck+3jj4Bs/tfo79nfvJjeVSU1JDbU5ttosmIjOAQuIY1dLbwkt7X+L5hudJpBLMKZpDfWV9toslIjOMQuIYs69zH8/uepZXm14lJ5ZDXUmdqpNEJNKYQsLMyoCPAacCvwae88jn3spMsr9zP0/ueJLXD7xOYV4hiysW684kETmiI4aEmS0APgH8V+B04ElgM/A3wMlm9jjwS+DX7t47dUWV8WjtbeXJHU+ycd9GivKKqK+s128ZRGTURnMm8RzwEPD/AL9x92Q4/H+bWSXwX4A/AT4JXD8FZZRxiKfi/Hb3b3lyx5PkxnJZUrlEZw4iMmZHDAl3X3aYcW3A3eHriMzsNuAS4IC7rwqHVQE/B5YCO4Er3L3Vgq+7PwDWAj3Ade6+cTTrOdbtaN3BA1sfoK23jQVlC3TNQUTG7YhfLc3sH8zsGjNbZWY5E1zfHcDFw4Z9E3jS3ZcTVGV9Mxz+MWB5+LoR+PEE1z3r9Sf7eeStR7h5w83gUF9Zr4AQkQkZTXVTMfAlgovVZmZbgFeAl8P3Te7ePZqVufszZrZ02OBLgQvC7p8CTwPfCIffGV4Y/72ZVZrZfHffP5p1jYW785l1n2FlzUq++YffJCc20Sycfvs79/PzLT/nYM9B6ivqj8ptOJoNtEed8tRgu9bD35Pp5MjDPBW0Tx22Pz2wrOHdg+1ekx5x+JDujOmilvGeNqqP0J/ZXOqQhp6ONG3QccRph693xHazI6Y97L4Zc+PRIyxjEu7POWI5x7iOzGuLAy0ifu0Pvjaush3OaKqb/iwsUAxYQXDx+kzgMuD/BCqAiXxdrcs48DcCdWH3QmBPxnQN4bD3hISZ3UhwtsGSJUvGXAAz49G3H2V763a++L4vHlWPnnB3Nu7fyLo31lGaX8qSirFv/7Eg7WniqaA96r5k3+BroNnTgWZQB4b1J4N2rBPpoA3sgbawh7wPdIf9k3EwmkwD7V7HLDZ4EBlotjWz+dbh3UPmD57zPuFpM6ePmnai6z3s5zCONqmHLeSIN3xMeB2Mvh3s6bypdNS3wLp7GthiZjuBBLAAyAFem6zCuLub2Zi33t1vAW4BWLNmzbg+vaqiKjr6O+iMdx41IRFPxfn3t/6dF/a+wMKyhRTkFmS7SNMimU7S0d9Be187HfEOeuI9dCe6g1fYtvVAd3ci6O9N9I76IF6QU0BBbgH5OfnkxfLIjeWSn5NPbiyXwoJC8mP55ObkkhfLC145wTR5sTxyYjnkWE7wntkdNqOa2T9kmrDN7KiD+pG6o4JAjh172vcceaJxGFVImFktQfXPfwU+ALwIPAh8y913TbAMTQPVSGY2HzgQDt8LLM6YblE4bEpUF1ezv2s/Hf0dU7WKSdXW18Y9m+9hX8c+llYunRV3LqU9TXtfOy29LbT0ttDe3057Xzvt/e1BKIT9A0+iHc4wivOKKckvoSSvhLKCMuaVzqMkr4SivKLBFvAKcwqD99zCwWEFOQWD/bPhsxSZLKP5ncRvgZXAE8BdwKfDu5omy0PAtcD/Dt8fzBj+38zsPuBsoH0qrkcMmFs8l+2t22ntbZ2qVUyaPe17uHPTnaQ9zZLKo6t6KZFKDLZp3dTdRFNXE4d6DtHSFwRD2tNDps+N5VJRUEF5QTm1xbUsr1pOeUF5MKwweC/JK6Ekv4TC3EId4EUm2WjOJM4luFZQAJwE9JrZy+6+c6wrM7N7CS5SV5tZA/BtgnC438xuAHYBV4ST/4rg9tdtBLfATulvMGqKa+hN9nKw5+BUrmbCXjvwGve9dh8VBRVUFFZkuziR+pP97O3cS0NHA/s69w0GQktvy5Bqn4qCCqqLq1lWuYw189dQVVxFVWEVVUVVVBZWUpxXrGoTkSwaTUhUAWeErzMJfjS3wsy6gVeBV9z9y6NZmbtfHTHqwhGmdYK7qqZFTXENvYlemrubp2uVY/bC3hdY9/o65pfNpyivKNvFGdSX7GNH6w62t22noaOBho4GmrubB8OgIKeAutI6ls1ZxrmLzqWutI66kjrqSusozC3MculF5HBGc3dTG8FtqU8PDDOzfGAVsJrgbqejXk1JDSlPsa9zX7aL8h7uzlM7n+Lxdx5nccVi8nPys1qelt4WtrVsY1vLNra3BsEwEAg1xTUsKl/E2QvPZlH5IhaVL2Ju0VydDYgcpUZzTeIfCZ7N9OzAw/zcPQ5sBDaaWZmZXQlUuftR+4O32pKg/YTmnmbiqXjWD8QD0p7mV2/9iud2P0d9ZT25sel/cG9fso83D73J1uatvN78Ok3dTUBwhrBszjLWLl/LCVUncFzlcTPqDEdEJm40R5ydBNcOVoQP83uQ4AF/HyGoejqN4JfSt01NEadHdXE1AN2JbrriXVQVVWW5RMGtng++8SDr961n6ZzpvYOptbeVV5pe4eX9L/N2y9ukPU1eLI8T557IB+s/yElzT2Jh2UL9aE9klhtNddMPgR+GD/NbC1wNfJfgUeHfI3joX2oqCzkdBkKiJ9FDZ39n1kMimU7ywOsPsKlpE/WV9dMSEO197by490U27N/AjrYdAMwvnc+Hl32YlTUrOX7O8XrMh8gxZiw/pmsD7glfs85ASHT1d9EZ78xqWQYC4tXGV1lauXRK6/MTqQSvNr3K8w3Ps6V5S3BbbcUSLj3pUs6cdybzy+ZP2bpFZOZTy3ShgZDoTfVyqOdQ1sqR9jQPvvEgm5o2TWlAtPa28vTOp3l297N0J7qpLKzkI8d/hHMXncu80nlTsk4ROfooJEIVBRXkxfLoS/TR2N2YlTK4O49te4yX9r00ZQGxo3UHT+54kg37N+DunDHvDM6rP4+Tq0/WD9FE5D0UEiEzo7Kwkt5kL/s7p+yH3Yf1uz2/45ldz1BfMfnXIN5peYeH33qYrQe3UpRbxIXHXcgFSy8YPIMSERmJQiJDVVEV3YluDvUcIpVOTeudO1ubt/Lwmw+zuGLxpK53R+sOHnzzQbYe3EpZfhl/vOKPOa/+PP2ITURGRSGRoaa4hrcOvYXjtPe3T9sdTvs693HP5nuYVzZv0u4eaultYd3Wdby470XK8su4fMXlnFd/3jHzpFgRmRwKiQx1pXVs2L8Bw2jra5uWkOiKd3HXq3dRml9KcV7xhJcXT8V5fNvjPP7O4wCsXb6Wjx7/UZ05iMi4KCQy1JXU0ZvsJZlO0tLTwrI5kc17T4q0p1m3dR3d8W4Wli+c8PK2Nm/lrs13cbDnIGsWrOGyky87atrGEJGZSSGRYeA3Acl0kobOBtawZkrX9+LeF9lyYAvHzTluQsvpjnfzi9d/wfMNz1NbUstXz/kqJ1WfNEmlFJFjmUIiw/zSICQS6QR7O6asfSMAGrsaeeStR1hUsWhCt7q+degtbnv5Ntr72/nYCR9j7fK1M+a5UyJy9FNIZBio8ulP9dPU3UQynZySB+olUgl+seUXFOcVj/uAnkqneOTtR3j07UepKa7hG3/4DZZWLp3cgspRzd2HtN2R2S7ywPCRhkXNP3z64fMcrt3lqKZjD1fGqHJGLe9w2xJlpHazo5Y51mWPZf2D/Yf5wni4NrTNjJ5Ez6SVK5NCIkN9RT0AHX0dVBdVc7Dn4JT8+vg3u35DY1cj9ZX145q/tbeVWzbewvbW7fzB4j/gylOu1IXpCUp7mlQ6RcpT73lPe3pwfJo07k7a00P+aDP/uA/3x4zB4DFmoDtj2EgHIMPee7DG3513YJrwgOb4kDauB35zY2bEiIEF0w9Mk9k9fPqBf4PzhMNGHG9Dhw/ZhmEHvxhBW9zvKUtEed6zzoF1ZU6b8VkM//wGhjnDQundDx53f+9B2jI7D38AP+x+z1jfiAGXUQ4YGp6ZQXm40F9csZii3Ml/CrNCIsPiisXELEZLbwtLK5dyoPvApIdEQ0cD/7njP1lUvmhc829r2cY/r/9n4qk4nz/z85y18KxJLd/Ryt1JppMk0gkSqQSJdIJkOkkynSTlqaHfGIf9jaY9Tcxig+1eD2/7uiCngPzcfPJzgldeLI/cWC4xiwUHMDNiFiPHcoL3WM7guMFpsCHTDhzYBg7IgwfmzAP3SAfEYd0jHTwzpxOZKIVEhqLcIsoLymnpbaEwt5AdrTs4re60SVt+PBXnF1t+QUVBxbiqsZ7Z9Qz3vnYv1cXVfPXcr7KgbMGklW2mcncS6QTxVPzdVzI++A0Uwm9n7hTlFVGaX0pVcRVl+WWU5ZdRml9KSX5JcKDPyScvJ4+8WB55OXmDB/38nPystNMhcjTQX0aGnFgOlYWVHOo9RFlBGdtbt0/q8p/e8TSHeg+xpGLJmOZzd3755i95bNtjrKpdxQ1n3jApv6mYCeKpOL2JXvqSfcRT8cGqEnj3NLskr4SKwgrml85nTtEcqoqqgoN/Xgkl+SUU5RZRlFekZ0+JTAGFxDBzi+byTus7FOUWsbt9N72J3klpbW1P+x6e2vnUmAMilU7xs1d/xvMNz3PekvO4atVVR1VDP+5Of6qfvmQffck++pP9g1UiaU9TkldCbWktdSV1zC2eS3lB+eDBvzivmOK8Yh38RbJIITFMXWkdG/dvDC5MmtHU3TThu4aS6STr3lhHZWHlmA7w/cl+bt5wM1uat/CJEz/B2uVrZ2xdc9rT9CZ66U5005vsHbyQ5zjlBeXUlgRBUFtSS0VhBRUFFZQXlOsxISIznEJimAWlC3Cctr42YhZjR+uOCYfEs7uepbGrcUzL6U/2c9NLN/H2obe55rRr+MCSD0yoDJPJ3YMW/OKd9Cf7B7/p15bUsqp2FQvKFlBVVEVlYWXwCHa1Zidy1FJIDDNQHXSo9xALyhaw+cBmPnTch8a9vD3te3hi+xMsLBv9Yzf6kn3c9OJNbGvZxufO/BzvX/j+ca9/MqQ9TVe8i46+juAWUJza4lreN/991FfWU1Ncw9ziubr4KzIL6a96mOVVywFo7mlmedVydnfsprO/k7KCsjEvqzvezT2b72FO4ZxRf5uOp+Lc9OJNvNP6DjesvoGzFmTnFtd4Kk5rbyv9yX7cnEVlizhrwVksqVjCvNJ5k3KdRkRmPoXEMCfMPYGYxTjQdWDwnvqGjgZW1KwY03JS6RQPvP4APYmeUT+8L5VOccuGW9jWsi0rAZFIJWjuaSaRSlCYW8iptaeysnYli8oXzZq7qURkbBQSw5TllzGncA4Hug8AUJhbyNaDW8ccEr/b8zveOPjGqB/e5+7ctfkuNh/YzKdP/fS0BUTa07T0ttAd76Ygt4CzFpzFqtpVLCpfdFTdRSUiU0MhMUxRXhFVRVU0dTcBQWt1m5s2c8mJl4z6OUu72nbx6LZHx/Sr6nVvrON3e37HJSdewvn154+r7GMRT8U50H2AVDrFSdUncc6icziu8jhdZBaRIRQSwxTlBiGxs20naU+Tl5NHf6qfHa07RvX47ebuZu7cdCdVRVWjPuA+vfNpHn/ncc6vP59Lll8y0U04rO54Nwd7DlKQW8CHln6IM+adwZyiOVO6ThE5eikkhhk4k0ikE4Ot05UXlPPs7mePGBItvS3c9vJt5MRyKC8oH9X6tjZv5edbfs5pdadx1aqrpux3EJ39nRzqOcScojlctuIyVtWu0m8UROSIFBLDDJxJAOzv3E9VURVzCuewvXU7jV2NkQ/8a+lt4daNt5JMJ6kpqRnVupq6mrhl4y3MK53HDWfeMCW/LO6Kd9Hc3Ux1cTWfPvXTnFxzsm5VFZFR0/MOhsmJ5Qw+MryhowEInrBZkFPAM7ueGXGe5u5mfrLxJ8RT8VEHRE+ihx+99CNiFuNLZ31p0h/13ZfsY1fbLhKpBJ8+9dP85Tl/yaq6VQoIERkTHTFGML9sPpUFlTR0NgwOqyutY8O+DayoXsGpdacODn/70Nvcs/ke8nLyRh0QaU/zLxv/hYM9B/nKOV+hurh60sqeTCfZ37mf3JxcLjnpEt43/31qqU5Exk0hMYLKwkrmlc4bPJOA4Hn/C8sXcu/me+mMd1JZWMmmxk1satpEbUktpfmlo17+Q28+xOvNr3PNadewfO7ySSmzu9PU3UR/sp8P1n+QDy75ICX5JZOybBE5dikkRjCncA41JTW81fIWiVRi8C6lwtxC5pXN41dv/wqAgpwCllYuHdO1hE1Nm3h026N8YPEHJu15TN3xbpq6m1hRvYK1y9eO+oxGRORIZkxImNlXgM8TtBu2GbgemA/cB8wFNgDXuHt8qstSVVRFTXENaU/T0NEw5AdxhbmFY37c94AD3Qe4/eXbWVKxhKtWXTXhcqY9zf7O/eTn5HPNadewsmbljH1KrIgcnWbEhWszWwj8BbDG3VcBOcBVwPeA77v7CUArcMN0lKe8oHzwh3DvtL4zKcuMp+LcvP5mYhbjT9/3pxP+0VpPoocdbTs4tfZUvnzOlzml9hQFhIhMuhkREqFcoMjMcoFiYD/wR8AD4fifAp+cjoKU5pdSXlA+2ADRRLk7d796N3s79/K5Mz834QvVB7oO0N7XzmdP+yyXn3L5mK6HiIiMxYyobnL3vWb298BuoBf4NUH1Upu7J8PJGoDRP297AgYOusfPOZ43D72Ju0/oW/ozu57h93t/z8dP/DiraleNezlpT7O7fTdLypdw5aor9UtpEZlyM+JMwszmAJcCxwELgBLg4jHMf6OZrTez9c3NzRMuT2l+KWnSLJ+7nPb+9sHnOI3HtpZt3LflPlbVrmLt8rXjXk48FWdH2w7OWXQON6y+QQEhItNiRoQEcBGww92b3T0B/Bvwh0BlWP0EsAjYO9LM7n6Lu69x9zU1NRO/s6cgt4D8nPzBx3BsPrB5XMtp7W3l5g03U11cPaFfVPcketjbsZfLVlzGJ076hB7CJyLTZqaExG7gHDMrtqBe50LgdeAp4PJwmmuBB6erQNVF1ZTmlTK/dD5bDmwZ8/yJVIKbN9xMf7KfP1vzZ+Nuj6Er3kVzTzPXnn4tZy88e1zLEBEZrxkREu7+AsEF6o0Et7/GgFuAbwBfNbNtBLfB3jpdZaourqYv2cepdafy5qE36ezvHPW87s69r93LjrYdXH/G9SwoWzCuMnT0d9DW18YXVn+Bk2tOHtcyREQmYkaEBIC7f9vdT3b3Ve5+jbv3u/t2d3+/u5/g7p9y9/7pKk91SRAS5yw8h7SneWnfS6Oe96mdT/HbPb9l7fK1nDn/zHGtv6O/g87+Tr6w+gssrVw6rmWIiEzUjAmJmaamuIakJ1lYvpAlFUt4bvdzuPsR59u4fyP3b7mfM+rO4OMnfnxc6+6Kd9HR18HnV3+exRWLx7UMEZHJoJCIUF5QTiz8eD609EPs7dzLq02vHnaerc1bufXlWzluznHcsHp8F6p7Ej209LZw/ZnXKyBEJOsUEhEqCytxgjOHsxeeTXVxNeveWEcilRhx+lebXuWml26irqSOL531pXE9ebUv2ceB7gNcc9o1o24bW0RkKikkIpQXlBOzGGlPkxPL4epVV7O/az/3bbmPtKcHp0ulUzy27TH+6aV/YkHZAr567lfH9QvoRCrBvs59XHnKlaNqJlVEZDrMiF9cz0Qxi1FXUkdPoofS/FJW1a7i4uMv5rF3HqOtt41zFp1DV7yLZ3c/y97Ovbxv/vv47OmfHVfjQal0it0du7n0xEs5fd7pU7A1IiLjo5A4jMXli9nUtGnwzOCTJ3+SsoIyHnnrEV5rfg2AeaXzuHH1jayev3pcj+5wd3a37+b8+vM5d/G5k1p+EZGJUkgcxqLyRby478XBfjPjomUXcX79+TR1N1GYW8jcorkTeq7Tno49nD7vdD56/Ef1FFcRmXEUEocxt3guxnsP3Hk5eYOPEp+IhvYGjp9zPJetuIycWM6ElyciMtl04fowqourB+9wmmz7OvexoGwBV596tdqgFpEZSyFxGCX5JVQUVNCb6J3U5TZ3N1NeUM5nTv/MuC50i4hMF4XEEZw490Q6+jsmbXktvS3kxHK4/ozr1ViQiMx4CokjOHHuifSl+iZlWe197cRTcT535ufUHoSIHBUUEkewsHwhk3FZoiveRWe8k8+d+TlqS2onvkARkWmgkDiC8oJy5pXOoyveNe5l9CR6ONRziOvOuG5S7ooSEZkuColRWD1/NS29LeOaN56K09jVyGdP/yzL5iyb5JKJiEwthcQoLJ+7fMjzmkYrlU6xp2MPf7zij/U8JhE5KikkRqG2pJb6ivoxnU24O7vad3HRcRdx1sKzprB0IiJTRyExShctu4j2vvZRTTsQEKvnr+bCZRdOcclERKaOQmKUls1ZxglzT6C5u/mw0w08sO+0utO4bMVl42p4SERkptARbJTMjLUnrKUn0UNfcuTfTQycQaysWcnlKy8nN6ZHY4nI0U0hMQbzy+bzqZWfYl/nPuKp+JBxyXSSnW07ObX2VK445QoFhIjMCjqSjdEZ888g6Ul++cYviRGjOL+YnkQPqXSKi46/iAvqL9ATXUVk1lBIjMOaBWtYWrmUVxpfoamriQVlCzi17lSqi6uzXTQRkUmlkBin6uJqLlp2UbaLISIypXRNQkREIikkREQkkkJCREQiKSRERCSSQkJERCIpJEREJJJCQkREIikkREQkkkJCREQiKSRERCSSQkJERCLNmJAws0oze8DM3jCzrWZ2rplVmdkTZvZ2+D4n2+UUETmWzJiQAH4APObuJwOnA1uBbwJPuvty4MmwX0REpsmMCAkzqwDOA24FcPe4u7cBlwI/DSf7KfDJbJRPRORYNSNCAjgOaAZuN7OXzewnZlYC1Ln7/nCaRqBupJnN7EYzW29m65ubD98GtYiIjN5MCYlcYDXwY3c/E+hmWNWSuzvgI83s7re4+xp3X1NTUzPlhRUROVbMlJBoABrc/YWw/wGC0Ggys/kA4fuBLJVPROSYNCNCwt0bgT1mdlI46ELgdeAh4Npw2LXAg1konojIMWsmNV/658DdZpYPbAeuJwix+83sBmAXcEUWyycicsyZMSHh7q8Aa0YYdeE0F0VEREIzorpJRERmJoWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiEQyd892GSaVmTUTtGI3HtXAwUksztFA2zz7HWvbC9rmsap395qRRsy6kJgIM1vv7iO1jjdraZtnv2Nte0HbPJlU3SQiIpEUEiIiEkkhMdQt2S5AFmibZ79jbXtB2zxpdE1CREQi6UxCREQiKSRERCSSQgIws4vN7E0z22Zm38x2eaaCmS02s6fM7HUz22JmXw6HV5nZE2b2dvg+J9tlnWxmlmNmL5vZI2H/cWb2Qri/f25m+dku42Qys0oze8DM3jCzrWZ27mzfz2b2lfD/9Wtmdq+ZFc62/Wxmt5nZATN7LWPYiPvVAj8Mt/1VM1s93vUe8yFhZjnAj4CPASuBq81sZXZLNSWSwNfcfSVwDvClcDu/CTzp7suBJ8P+2ebLwNaM/u8B33f3E4BW4IaslGrq/AB4zN1PBk4n2PZZu5/NbCHwF8Aad18F5ABXMfv28x3AxcOGRe3XjwHLw9eNwI/Hu9JjPiSA9wPb3H27u8eB+4BLs1ymSefu+919Y9jdSXDgWEiwrT8NJ/sp8MmsFHCKmNki4L8APwn7Dfgj4IFwklm1zWZWAZwH3Arg7nF3b2OW72cgFygys1ygGNjPLNvP7v4M0DJscNR+vRS40wO/ByrNbP541quQCA6UezL6G8Jhs5aZLQXOBF4A6tx9fziqEajLVrmmyP8HfB1Ih/1zgTZ3T4b9s21/Hwc0A7eHVWw/MbMSZvF+dve9wN8DuwnCoR3YwOzezwOi9uukHdcUEscYMysF/hX4S3fvyBznwf3Qs+aeaDO7BDjg7huyXZZplAusBn7s7mcC3QyrWpqF+3kOwTfn44AFQAnvrZaZ9aZqvyokYC+wOKN/UThs1jGzPIKAuNvd/y0c3DRwGhq+H8hW+abAHwKfMLOdBNWIf0RQX18ZVkvA7NvfDUCDu78Q9j9AEBqzeT9fBOxw92Z3TwD/RrDvZ/N+HhC1XyftuKaQgJeA5eGdEPkEF7weynKZJl1YF38rsNXd/zFj1EPAtWH3tcCD0122qeLuf+Xui9x9KcF+/U93/xPgKeDycLLZts2NwB4zOykcdCHwOrN4PxNUM51jZsXh//OBbZ61+zlD1H59CPhseJfTOUB7RrXUmOgX14CZrSWou84BbnP372a3RJPPzD4APAts5t36+W8RXJe4H1hC8Ij1K9x9+MWxo56ZXQD8d3e/xMyWEZxZVAEvA59x9/4sFm9SmdkZBBfq84HtwPUEXwhn7X42s78BriS4i+9l4PMEdfCzZj+b2b3ABQSPBG8Cvg38khH2axiWNxFUu/UA17v7+nGtVyEhIiJRVN0kIiKRFBIiIhJJISEiIpEUEiIiEkkhISIikRQSIpPAzJaYWVf4wEiRWUMhITJOZrbTzC4CcPfd7l7q7qlsl0tkMikkREQkkkJCZBzM7GcEv3J9OKxm+rqZ+cCzgszsaTP7WzP7XTj+YTOba2Z3m1mHmb0UPo13YHknh43GtFjQANYVWdo0kSEUEiLj4O7XEDwz6OPuXkrwaIThrgKuIXg8xPHA88DtBI+J2ErwWAXCR3k/AdwD1Ibz/dMsbfxKjjIKCZGpc7u7v+Pu7cCjwDvu/h9hGwe/IGjTA+ASYKe73+7uSXd/meBpvZ/KTrFF3pV75ElEZJyaMrp7R+gvDbvrgbPNrC1jfC7wsyktncgoKCRExm+yno65B/iNu394kpYnMmlU3SQyfk3AsklYziPAiWZ2jZnlha+zzGzFJCxbZEIUEiLj93fAX4fVRJcfYdpI7t4JfITggvU+graKvwcUTEIZRSZE7UmIiEgknUmIiEgkhYSIiERSSIiISCSFhIiIRFJIiIhIJIWEiIhEUkiIiEgkhYSIiET6/wHO5lRTgYvwiQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "y = sol.y\n",
    "N = y[0]\n",
    "stdN = np.sqrt(y[1] - N**2)\n",
    "plt.plot(timepoints, y[0], color=\"g\")\n",
    "plt.ylabel(r'$\\left< N \\right>$', size=12)\n",
    "plt.xlabel('time', size=12)\n",
    "plt.title(\"Expected compartment number\", size=12)\n",
    "plt.fill_between(timepoints, N-stdN, N+stdN, alpha=0.4, color=\"g\")\n",
    "plt.ylim((50,150))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}