Comparison of Python and Cython ODE integration

compare_cython_with_solve_ivp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare Python and Cython implementation of ODE integrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from ode_cython import solve_ivp as solve_ivp_cython, Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare Runge-Kutta methods on the restricted three body problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = 0.0122771\n",
    "nu = 1 - mu\n",
    "\n",
    "t0 = 0\n",
    "tf = 17.0652165601579625588917206249\n",
    "y0 = [0.994, 0, 0, -2.00158510637908252240537862224]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define function computing the system's rhs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple Python function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, y):\n",
    "    D1 = ((y[0] + mu)**2 + y[1]**2) ** 1.5\n",
    "    D2 = ((y[0] - nu)**2 + y[1]**2) ** 1.5\n",
    "\n",
    "    return np.array([\n",
    "        y[2],\n",
    "        y[3],\n",
    "        y[0] + 2 * y[3] - nu * (y[0] + mu) / D1 - mu * (y[0] - nu) / D2,\n",
    "        y[1] - 2 * y[2] - nu * y[1] / D1 - mu * y[1] / D2\n",
    "    ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python class inherited from `Function` extenstion type provided by `cython_ode`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class MyFunction(Function):\n",
    "    def evaluate(self, t, y):\n",
    "        return fun(t, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another extention type inherited from `Function`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "from ode_cython.function cimport Function\n",
    "from libc.math cimport sqrt\n",
    "import numpy as np\n",
    "\n",
    "cdef class FunCython(Function):\n",
    "    cdef double [::1] f\n",
    "    cdef double mu, nu\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.f = np.empty(4)\n",
    "        self.mu = 0.0122771\n",
    "        self.nu = 1 - self.mu\n",
    "        \n",
    "    \n",
    "    cpdef double[::1] evaluate(self, double t, double [:] y):        \n",
    "        cdef double x = sqrt((y[0] + self.mu)*(y[0] + self.mu) + y[1]*y[1])\n",
    "        cdef double D1 = x * x * x\n",
    "        x = sqrt((y[0] - self.nu)*(y[0] - self.nu) + y[1]*y[1])\n",
    "        cdef double D2 = x * x * x\n",
    "\n",
    "        self.f[0] = y[2]\n",
    "        self.f[1] = y[3]\n",
    "        self.f[2] = y[0] + 2 * y[3] - self.nu * (y[0] + self.mu) / D1 - self.mu * (y[0] - self.nu) / D2\n",
    "        self.f[3] = y[1] - 2 * y[2] - self.nu * y[1] / D1 - self.mu * y[1] / D2\n",
    "        \n",
    "        return self.f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fun_cython = FunCython()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fun_python = MyFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, [t0, tf], y0, method='RK45')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ts, ys = solve_ivp_cython(fun_python, [t0, tf], y0, method='RK45')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(res.y, ys.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x109382588>]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAFkCAYAAACThxm6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xdc1XX7x/HX1z0wR4oezRE4MleCWgaV7aUNqJts3lnZ\nurNsrzuzZfPWxt24K22bdkfzbtv6QZkFWVlaJq70BJoDQcXB5/fHJTESFTyH7znwfj4e56GcxQUK\n530+4/p4zjlEREREQqWe3wWIiIhI7aJwISIiIiGlcCEiIiIhpXAhIiIiIaVwISIiIiGlcCEiIiIh\npXAhIiIiIaVwISIiIiGlcCEiIiIhpXAhIiIiIRXWcOF53kGe573ped4yz/OKPc87YSf3P2Tb/cpe\ntnqeFxvOOkVERCR0wj1y0RyYDVwC7OohJg7oAXTYdgk45/LCU56IiIiEWoNwPrlz7j3gPQDP87wq\nPHSFcy4/PFWJiIhIOEXimgsPmO153nLP8z7wPO9AvwsSERGRXRfWkYtqCAIXAt8AjYELgE89zxvi\nnJu9vQd4nrcncDSwCNhYQ3WKiIjUBk2AbsD7zrk/QvWkERUunHO/AL+UuWqm53nxwFjgnEoedjTw\nYrhrExERqcXOAF4K1ZNFVLioxCwgaQe3LwJ44YUX6N27d40UJDB27FgmTpzodxl1ir7nNU/f85qn\n73nNmjt3LmeeeSZsey0NlWgIF/th0yWV2QjQu3dvEhISaqYioWXLlvp+1zB9z2uevuc1T99z34R0\nWUFYw4Xnec2B7tgiTYA4z/MGAKucc0s9z5sAdHTOnbPt/pcDC4EfsXmgC4BDgSPDWaeIiIiETrhH\nLgYBn2C9KxzwwLbrnwVGYX0sOpe5f6Nt9+kIrAe+Bw53zn0e5jpFREQkRMLd5+IzdrDd1Tl3boWP\n7wPuC2dNIiIiEl6R2OdCosDIkSP9LqHO0fe85ul7XvP0Pa8dPOd2tSt3ZPI8LwHIysrK0iIgERGR\nKsjOziYxMREg0TmXHarn1ciFiIiIhJTChYiIiISUwoWIiIiElMKFiIiIhJTChYiIiISUwoWIiIiE\nlMKFiIiIhJTChYiIiISUwoWIiIiElMKFiIiIhJTChYiIiISUwoWIiIiElMKFiIiIhJTChYiIiISU\nwoWIiIiElMKFiIiIhJTChYiIiISUwoWIiIiElMKFiIiIhJTChYiIiISUwoWIiIiElMKFiIiIhJTC\nhYiIiISUwoWIiIiElMKFiIiIhJTChYiIiISUwoWIiIiElMKFiIiIhJTChYjUmNyCXJInJxP/UDzJ\nk5PJK8zzuyQRCQOFCxGpEb//Doc+kUrm0kxyVueQuTSTlGkpfpclImGgcCEiIZeXB+++C7ffDied\nBHvtBYEAzF0SLHe/r34KcsMNMHs2OOdTsSIScg38LkBEotvKlZCVBd98Y5esLFi61G5r3RoSE+Gs\ns2DQIJiQFyArL+fPx+7ZKMB/noC774aePSEtDU47Dfbd16cvRkRCQuFCRHbZqlV/DRKLF9ttLVta\ngBg50v5MTIS99wbPK338QYXppExLIVgQJBATID0tndbjYcYMePlleOghG+3o29eCRloa9Ojhz9cq\nItXnuSgfi/Q8LwHIysrKIiEhwe9yRGqN1astPJSEiawsWLjQbttjDwsPiYmlQSI+vnyQqI6iInj/\nfZg2Dd58EwoKYOBACxl/+5uFFREJnezsbBITEwESnXPZoXpejVyICGvWQHZ2+SCxYIHdFhNj4eHk\nky1IDBpkQaJeGFZsNW4MJ5xglw0b4J13LGiMHw/XXw9DhpQGjb32Cv3nF5HQULgQqWPy8/8aJObP\nt9uaN4eEBBgxojRI9OgRniCxM02bQmqqXQoK4O23berkhhvgqqsgOdmCximnQIcONV+fiFRO4UKk\nFlu3Dr79tvw6iV9+sduaNbMph+OOK53e6NkT6tf3t+btiYmxhZ6nnQZr18Ibb9iIxtixcPnlcMgh\nFjRSU6FtW7+rFRGtuRCpJQoKbEtn2SDx88+2xbNpU9hvv9L1EYMGwT77RGaQqIpVq+C11yxofPyx\nXXf44RY0Tj7ZdquISOXCteZC4UIkCq1fb0GiZFrjm29g3jwoLoYmTWDAgPJBondvaFDLxylXrIBX\nX7Wpk88/t6/36KMtaJxwgi1CFZHyonJBp+d5BwHXAIlAADjJOffmTh4zDHgA6AMsAe50zj0bzjpF\nItmGDfDdd+WDxE8/WZBo1MiCxCGH2DqExETrEdGwod9V17x27eCii+yyfDn89782onHWWbZQ9Ljj\nLGgMH25rS0QkfML9XqY5MBt4Gkjf2Z09z+sGvA08CpwOHAE85Xnecufch+ErUyQybNwI339fPkj8\n+CNs3WqBYcAASEqydQaJidCnjwUMKa9jRxgzxi5LlsArr1jQOO00W2syfLgFjWOPtSkjEQmtsIYL\n59x7wHsAnrdLO+AvBnKcc9du+/hnz/OSgbGAwoXUKkVF8MMP5RtSzZkDW7bYkH7//rD//nDppRYk\n+va1d+BSNV262KjOVVdBTg5Mn25TJ6mptlD0xBMtdBx1lIKaSKhE2izsAcBHFa57H5joQy0iIZFb\nkEvK9FQW/xGkyeYAQ5ek89PXsfzwA2zebEGib18LEBdeaGsk+vVTkAiHuDjrl3H99bbYddo0u7z4\nIrRqZYtA09LgsMPq5tSSSKhEWrjoAORWuC4X2MPzvMbOuSIfahLZLanTU/liaea2j3IINkzhtP0y\nOO88CxL9+9siTKlZvXrBLbfYZc6c0qAxZQrsuScce2ou3/ZMZX29IB1bWKvy2OaxfpctEhUiLVxU\n29ixY2nZsmW560aOHMnIkSN9qkjE/BIsfxJohx5Bnh7jUzGyXX372uW222wXzrRpMCk/laJ8C4UL\n1+SQMi2FjFEZPlcqUn1Tp05l6tSp5a5bu3ZtWD5XpIWL34H2Fa5rD+TvbNRi4sSJ2ooqEWfJEli1\nOACdSk8CDcQEfKxIdsTzrLHYwIHwyoNBctaU3hYsCFb+QJEosL033GW2ooaUD019d+hL4PAK1x21\n7XqRqFJcDOecA+0/TWf/QBJxreNI6pxEetpON05JBAi0KB8CFQpFdl24+1w0B7oDJTtF4jzPGwCs\ncs4t9TxvAtDROXfOttsfBy71PO8eYDIWNE4BjgtnnSLhMGkSfPopfPxxLIcequH0aJOels5Bj6Tw\ny/IgQ/sFFApFqiDc0yKDgE8At+3ywLbrnwVGYQs4O5fc2Tm3yPO847HdIWOA34DznHMVd5CIRLQ5\nc+yArbFj4dBD/a5GqiO2eSzX7JnBBddDxlZ/Dm8TiVbh7nPxGTuYenHOnbud6z7HOnqKRKWiIjjz\nTOjeHe66y+9qZHds2GBbghUsRKom0hZ0ikS9W2+19txffaUtptFu40Z18BSpDoULkRDKyIB77rER\ni4ED/a5GdteGDQoXItWhwT6REMnPt0OyDjwQrrnG72okFDZs0OiTSHVo5EIkRMaOhZUrYcYMqF/f\n72okFDRyIVI9ChciIfD66zB5Mjz1lJ1fIbWD1lyIVI+mRUR2U24uXHABnHACjBrldzUSShq5EKke\nhQuR3eAcnH++tY1+8kn7U2oPrbkQqR5Ni4jshqefhrffhjfegFgdmFnraFpEpHo0ciFSTQsWwBVX\n2MjFCSf4XY2Eg6ZFRKpH4UKkGrZssW2n7dvDv/7ldzUSDrkFuXzVJ5l3esSTPDmZvMI8v0sSiRoK\nFyLVcO+91oHzueegRQu/q5FwGP58KutaZ1LQMIfMpZmkTEvxuySRqKFwIVJF2dkwbhxcdx0kJfld\njYTD++9D9vxgueuCBcFK7i0iFSlciFTBhg12KFnfvnaGiNQuxcUwfjwceyzs4QXK3RaICVTyKBGp\nSLtFRKrghhsgJweysqBRI7+rkVBaudKC4wcfWMC44Ip0TnklhWBBkEBMgPS0dL9LFIkaChciu2jG\nDHjwQZg4Efr08bsaCaVZs+DUU6Gw0KZEjjwSIJaMURl+lyYSlTQtIrILVq+Gv/8dDjsMxozxuxoJ\nFefgsccgORkCAVtPY8FCRHaHwoXILvjHP2DdOpgyBerpp6ZWKCyEs8+GSy6BCy+Ezz+HLl38rkqk\ndtC0iMhOvPwyvPQSvPCCXnxqi59/htRUWLjQ/m1HjvS7IpHaRe/BRCqRW5DLkMeTOeOreNpek8wR\nJ6qJUm3w6qsweLA1Qps1S8FCJBwULkQqkTo9la9zMylulcPK5pmkTlcTpWi2eTNcdRWccoptNf36\nay3MFQkXTYuIVKJi06SFK9VEKVotXw5paTBzJkyaZItydYKtSPho5EKkEhWbJi3/OcAdd8DWrT4V\nJNXy6acwcKD1J/nsM7j8cgULkXBTuBCpRHpaOkmdk4hrHceBeyVxZad0xo2z7ahLl/pdneyMc3DP\nPXD44dZR9dtv4cAD/a5KpG7QtIhIJWKb/7WJ0olHWBfH/v3hySdt/l4iz5o11pfkjTfgxhvhttug\nfn2/qxKpOzRyIVIFBx8M330HRxxhHR3PP9/6JUjkmD0bBg2yKZA334Q771SwEKlpChciVdS6NUyf\nDk8/DVOnQkKCdXYU/02ZAkOHwh572PkvI0b4XZFI3aRwIVINngejRlmoiImBAw6A+++3UzWl5m3c\nCBdcYP8mZ54JX3wBcXF+VyVSdylciOyGXr3gyy/hiivgmmvgmGMgqB2rNSonxxZqvvACTJ5sa2Ga\nNPG7KpG6TeFCZDc1agT33gsffghz5thiz7fe8ruquuHttyExEdautZB37rl+VyQioHAhEjJHHGGL\nPYcOhRNOsMPONmzwu6raJ7cgl6Snk2l9azwj0pPZ/7A8srJgv/38rkxESihciIRQu3a2/fHf/7YF\nn4MHww8/+F1V7XL8c6l88Vsma7wc6JrJuuNTaNXK76pEpCyFC5EQ8zw7xvubb+x49sGD4ZFHrKmT\nVN+mTXDXXZD1S/lFLb8XaJGLSKRRuBAJkz597NTN0aPhsstsqmTFCr+rik6ffWbTHuPGQac9yrdl\nr9imXUT8p3AhEkZNmsBDD9nCw5kzoc/+ufT9VzLxD8WTPDmZvEId474jK1ZYp81hw6y/SHY2ZF9X\n2pY9qXMS6WnpfpcpIhWo/bdIDTj+eFt70fu+VH5clwlAzuocBt+fwgN9MthnH+jRAxo39rnQCFFc\nbA2xrr3WppOefNJ6WNSrB/DXtuwiElkULkRqSIcO0KZLkDVrSq/7bW2QU0+1v9evb42f9tkHeve2\nS8nfW7b0p2Y/zJkDF18MGRlwzjlw3322UFZEoofChUgNCrQIkLMm58+Ph/YN8PoKmDsX5s2zP+fO\ntfbiixaVeVxg+6GjY8fac3x4YSHcfjs88ADEx8Mnn9h0iIhEH4ULkRqUnpZOyrQUggVBAjEB0tPS\nadscDjrILmWtXw8//1w+dHz+OTz1lO2cAGjRYvuhIz4eGkTRT/f//geXXgq5uXDrrXD11ZoiEolm\nUfTrRyT6be8Y98o0awYDB9qlrC1bYOHC8qFj3jzrr7F2rd2nYUPo3v2voaNXLzsLJVL89htcfjmk\np8NRR8GMGRaMRCS6KVyIRJkGDWzxZ48e5U/9dA5+//2voeOZZ2DZstL7de7819CxZ9dcLvo4tdyI\nSmzz2LB9DVu2wMMPwy23WNiZNs2OsK8tUzwidV3Yw4XneZcCVwMdgO+Ay5xzX1dy30OATypc7YCA\nc0579kR2wPNsbUYgAIceWv62/HybYikbOt57z5p7bd0KnJsKXUt3sfS9PYUzN2XQrh20bWsLKste\nWrWqfhCYNQsuvNBapV96KdxxR91asCpSF4Q1XHielwY8AIwGZgFjgfc9z+vpnFtZycMc0BNY9+cV\nChYiu2WPPaxT6ODB5a/ftAkWLIBDXwuSu7n0+nUuyLvvWp+JVav+2l20QQMLHdsLHiWXsrftuScU\nFMCNN8Ljj9tUz6xZMGhQ+L92Eal54R65GAs84Zx7DsDzvIuA44FRwL07eNwK51x+mGsTqfMaNbJp\nke5fBshdWrqLJbFXgIy59vctWyxgrFhhl5UrS/9e9jJ3bultW7bs+PN26gT/+Q+89lrloUTHpotE\nr7CFC8/zGgKJwF0l1znnnOd5HwFDd/RQYLbneU2AOcCtzrkvwlWniGx/F0uJBg0gNtYuu8I5W1i6\nYgV88YV12CyRlmaBZsUK67ZZEky2d3psTMyOR0MqXmJibKomtyCX1Ok1t35ERP4qnCMXbYH6QG6F\n63OBXpU8JghcCHwDNAYuAD71PG+Ic252uAoVqeuqsotlZzwPmjaFl16CCRNg771tbcdxx1X+mMLC\nykdESi4//2yNtVasKN0VU1bjxhY+1qSkUrhn6fqRo55OYeboDI2EiNSgiNot4pz7BfilzFUzPc+L\nx6ZXzvGnKhGpihkz7FTYhQvhmmvgpptsW+2ONG9ul65dd+1zbNpUPoyU/fuDrvwpqd8tCNKihR0k\nN3AgJCTYnwMGWJ8QEQm9cIaLlcBWoH2F69sDv1fheWYBSTu709ixY2lZYcn5yJEjGTlyZBU+lYhU\nV24uXHUVvPiiNQR77TXYd9/wfK5Gjaw7aceOf73t48kBMsusH+kfF+Cih+Dbb20q5qWXLJx4nm3n\nTUgoDRwDB9riU5HaaOrUqUydOrXcdWu3NwwYAp6ruAw8lE/ueTOBr5xzl2/72AOWAA855+7bxef4\nAMh3zp1Sye0JQFZWVhYJCQkhqlxEdlVxsR0sdv31dj7K/ffbmSB+9azIK8z7y/qRsmsuNm+Gn36y\noFESOGbPtqkZgC5dygeOhATb3qseHFIbZWdnk5iYCJDonMsO1fOGe1rkX8AznudlUboVtRnwDIDn\neROAjs65c7Z9fDmwEPgRaIKtuTgUODLMdYpINXz3HVx0kR0nf955cM89/r/z39n6kYYNbUpkwAA4\n91y7butW+PXX0rCRnQ2TJtkuGbDFrBUDx957K3CIVCas4cI5N93zvLbAbdh0yGzgaOfcim136QB0\nLvOQRlhfjI7AeuB74HDn3OfhrFNEqqagwM4AmTTJunz+3/9BcrLfVVVf/frWGr1XLzjtNLvOOViy\npHzgeOYZuGvb/reWLcuv4UhIsMfXr+/blyESMcI6LVITNC0SHUq2By5fF6RjC20PjGavvw6XXQZ/\n/AHjxsHYsbYGoq7IzS0fOL79FnK2LfFo2tRGRErCRkKCLSTVIWwSqaJ1WkQE5+CwJ1L5qcC2By5c\nk0PivSm8enwGgwZBvXo+Fyi7ZPFiGDMG3nwTjj/etpd26+Z3VTWvfXs45hi7lFizxtZtlASOzz6D\nJ56w9SgNGkDfvuUDx4ABtjtGpLZSuJCw2bABnn/ehs7nHhmENqW3Lc8Psv/+0KGDvVCNGAFHHKFf\nuJFo82b7N7z1VmjdGl59FU4+WesNymrVCoYNs0uJ9evh++/Lj3C88IJ9Pz3PplDKTqsMHAht2lT2\nGUSii8KFhFwwCI8+Co89ZgviTjwR6u8dYM7a0u2BB/QNcNen8NZbdnn6aRs6PvxwCxrDh8Nee/n2\nJcg2X3xhCzZ//NGORh8/Xr0hdlWzZnDAAXYpsWlT6U6VksDxxhsWRMBGgiqu4wgEfClfZLdozYWE\nzOzZMHEiTJ1qQWHUKBtGj4/f+fbAX34pDRoZGbZ6f+BACxojRtgvWU2f1JxVq2xr6ZNP2mFnTzxh\n/x4Selu3wvz55QNHdrZNtYCN7lUMHN26aeRIQiNcay4ULmS3FBfD//4H//oXfPqp9QgYM8a2JbZq\nVb3nXLXKjgN/6y14911r9dyxo41mjBhhoxtNm4b0y5BtnLOh+6uusnfZEybA6NHaAVHTnLM1LmXD\nRnY2/L6t/WCrVn8NHC075vK3V3WmilSNwkUlFC78UVho2/IefNDedR1wgO0aSEmxBWyhsnmzbXMs\nGdVYsMCCxRFHlE6faNh495Ts5FmyKsjaZQHyn0xn5Amx/Otf9q5ZIkcwaGGjbOBYtMhuq3d+MsV7\nZf5536TOSSE7L0ZqL4WLSihc1KzffoOHH7bjstetg9RUCxVl55XDxTmYN680aHzxhY2cDBpUOn2y\n334aLq6q5MnJZC4tfVHq0yKJOVfqRSlarFplU5KnZsSzypWua2pUGMe/ey5g5EgtlJbKhStcaBZb\ndsnXX8Ppp1tXwscft2mPBQtg2rSaCRZgoaF3b7j2WhvNyM2F556zmu6/34aGu3a1Q7PefRc2bqyZ\nuqKZczBncfmDvjY0CFZyb4lEbdrAYYdB773KD+E1Lw4werRNKV52GcyZ41OBUicpXEiltm61bYfJ\nyTBkCHz1FTzwgI1e3H//rp9gGS5t28JZZ8H06XYq5ocf2hbJd9+1473btrWPJ0+2ICLlFRfbi87a\npeVflAIxmmeKRulp6SR1TiKudRxJnZOYNz6dnBz4xz/glVegXz87UO7FF6GoyO9qpbbTtIj8RX6+\nvSA/+KDN5x58sE19jBgRHQv7nLOtkyXTJzNn2vVDhsChw3N5r2Uq+cVBAnW4U2hxMVx4oW0Bvu/R\nPF5rVPlOHol+mzbZltfHH4ePP7bgfe65tli3e3e/qxM/ac1FJRQuQmfRInjoIXjqKWuAddppcMUV\nYP/voldeHrzzjgWN19to0duWLbZN+MUXYcoUOPtsvyuSmvTzz7a1+JlnYPVqOPJIuPhie/MQysXY\nEh205kLCwjlbGHnKKdaP4tlnbRh10SLrrhntwQLsRMu//92meLr2Lb+eYNnaurW+YPNmOOMMeOkl\nuyhY1D29etnW8WXL7Oe9oMB2eXXtamfFLF3qd4VSGyhc1FGbN8PLL9tizKQkW+z173/bL5a77oJO\nnfyuMDw6tii/nmDNbwG2bPGpmBpWVASnngqvvWZz8Glpflckfmra1MLlF1/YbpMTTrDQ0a0bnHSS\n9ZopLva7SolWChd1zOrVcO+9EBcHI0daK+e337aWxBddZC2La7Oyi976tEhi3VPpjBljIzi12YYN\ntrj1vffsVNOTT/a7IokkAwZYu/7ly611/6JFcOyxth7j7rttalGkKjTDVkfMn28LNJ95pnRo/Ior\noH9/vyurWbHNY8utsXhqD7jgAhsqvvxyHwsLo8JCO9/liy8sSB5xhN8VSaRq0cIW+o4ebbvDHn/c\nzpO55RbraXPRRbbAW71kZGc0clGLOWctuU880V48p0+Hq6+GJUtsN0hdCxbbc/75cM01cOWV9sJb\n2+Tn2zvQr76yUQsFC9kVnmdTps88Y2sz7r3XuoIOGwZ9+tjC79Wr/a5SIpnCRS20aZM1l0pIgEMP\ntWZXTz5poeLWW6F9e78rjCx3323zzaedBt9953c1obNmDRx1lB37/eGH9o5TpKratLFRzrlzbRtr\nv3529kynTrbraNas2j+tKFWncFGLrFwJd95pq77POcfO3PjgA/jhB+uo2aSJ3xVGpnr17LCuXr3s\nrJJgLdhA8scfdsDb/PkwY0bNdVGV2svz7M3KtGm28Pvmmy1s7L+/teB/8knbeSICChe1wty5Nk/a\nuTPccYdNg/z0k/V2OPJIzY/uiubNrQ+GczaKsX693xVVX26uvQgsXQqffFI7thNLZOnQAW680UZF\n//c/G8W46CJrNX7ppfaGRuo2hYso5ZwNdR97LOy7r70w3nyzvaA8/ridwSFV07GjfR/nzrW24tG4\nDW/5cpsXX7kSPvtM62okvOrXt1b7b74JCxfaouj0dPt/l5xsI4I646duUriIErkFuSRPTiZuUjw9\nJiSz7+A8jjqq9PCuRYvgppusra9U38CBMHWq9YK48Ua/q6maJUtsXUVhoQULBUypSV26wO232//D\n//7X+micdRbstZctmp4/3+8KpSYpXESJlGmpZC7NZOHaHH7dlEneoSl88glkZdkPcKNGfldYe4wY\nYQe03XOP7aqJBjk5FiyKi+Hzz6FHD78rkrqqYUPbtvrhh/DLL9Ydd/Jk6NnTpmlffdW2w0vtpnAR\nJRavKr/KsFXnIMOGaT1FuFxxha1jufBCW7cQyX7+2YJFo0YWLLp187siEdOjh52gvGyZHSewfr0d\nNdC1q/XOUKvx2kvhIkq0aahjsWuS58HDD9vCyNRUewGPRD/+CIccAi1b2lTIXnv5XZHIXzVpAmee\nCZmZtt375JNh0iQLwiecAO++C1u3+l2lhJLCRZQYt086LE6ic0wcSZ2TSE9L97ukWq9hQ2s81qED\nHH+8be+MJLNn2+LNDh2sWVpAeVOiQP/+do7R8uW2+Py332xRaPfuMGECzFlo68viH4oneXIyeYXq\nPR6NFC6iRDMXC1My+PK0BWSMyiC2eazfJdUJrVpZ5861a+3dVlGR3xWZr7+2UZVu3azXQLt2flck\nUjUxMdZ6PyvLOsgeeqgtCO1/p60vy1mdQ+bSTFKmpfhdqlSDwkWUKHlRa9zY3zrqorg4eOMN60Q4\nerT/3QgzM61BVu/e8NFH1kFRJFp5HgwZYos+ly2DNl3Lry/79fegpkyikMJFlFC48NeBB9ovv+ee\ns6Fbv3z6KRx9tLV2f/99W2shUlu0bg37dCo/v5f7a4DeveGppyJn5FB2TuEiSihc+O/00+1slptu\nsrUYNe2DD6xp2oEHWvfVFi1qvgaRcEtPSyepcxJxrW192Xuj0unf30YN997bdp/k5/tdpeyMjlyP\nEiVd7ho29LeOuu6WW2zv/jnn2Ha6/fevmc/79tu2a+XII61Bkc6JkdoqtnksGaMyyl13dLLt2Lrv\nPmtud+edcMkl1hE0VsvPIpJGLqJEUZGNWqivhb88D55+2qYlTjwRFi8O/+d89VVbTDp8uLVWVrCQ\nuqhXL5saWbjQDmJ86CEL+JdeatdJZFG4iBIl4UL816QJvP46NGtmL/jhHKJ96SVIS4NTT7XTKNWJ\nVeq6Tp1samTJEjtP6ZVXrFnXGWfA99/7XZ2UULiIEgoXkaVdO5uqWLoUTjsNtmwJ/eeYMsUaD511\nlnU3bKBJTJE/tW5t658WLYIHH7RdVAMGWM+Mzz/3f1dXXadwESUULiLPvvva+ocPPoCxY0P73I8/\nDqNG2SLNcPlsAAAgAElEQVS2p5+20ydF5K+aNbOpkfnz7RTWpUuta21Skp3WGo2nG9cGChdRoqhI\nc+2R6Igj4NFH4ZFHrF14KEyaBBdfbIvVHnsM6umnVGSnGjYsnRr53/9spO/EE6FfP9tCrsPSapZ+\nbUUJjVxErtGj4cor7bCzd97Zvee6+24bBbnuOpg4UQt4RarK80qnRjIyID7ednfFx9v0SWGh3xXW\nDQoXUULhIrLdey8ccVIuI15LpvN9VT8TwTnroXHDDTBunDXqUrAQ2T0lUyM//GDn8Fx1le0wGT8+\n8s4Kqm0ULqKEwkVkq18f8o9JpXivTH5bX7UzEZyzvfvjx8Ndd1nIULAQCZ2+fW1qZMECa4Z3zz3Q\npYuNEurY9/BQuIgSCheRL29D+TMRlq8LVnJPk1tgpz+2GhfP3bnJ3PZAHjfcEM4KReq2rl2tP8bi\nxXD11fDss3Z20N//Dj/95Hd1tYvCRZTYuFHhItIFYsqfibA+d8dnoKdOt9Mf8+vnQNdM3m+l0x9F\nakK7djZSuGSJTWl+9BH06QMnnQQzZ/pdXe2gcBElNHIR+cqeidC9URK5D6Xz2GOV33/xqvIjG8GC\nHY90iEhoxcTY1EhOjh1MOG8eDB1q6zPefVe9MnZH2MOF53mXep630PO8DZ7nzfQ8b/BO7j/M87ws\nz/M2ep73i+d554S7xmigcBH5Ss5EWDBmAfNvyGDMebGMGWMnmVa0aROsXlJ+ZKPiyIeI1IxGjeDc\nc21qJD0dNmywHScDB8LUqeFpkrcjJVOm8Q9tf3H4jm7f2WNriufCGM08z0sDngVGA7OAscCpQE/n\n3Mrt3L8bMAd4FHgaOAKYBBznnPuwks+RAGRlZWWRkJAQhq8iMhx8sM0XPv+835XIrtqyBY45BmbP\nhq+/thMdS/zznzDhoTz63pHCOhckEBMgPS2d2OY6hUnEb87BZ5/Z1vD337ef3WuusbUZTZuG9/Mu\nWwbHTE/mx3WZf17fYVMSx+dmsHGjTZF/3C2Z1S1Kbw9sTuL8ehm0bg2PFSUzv6j0tsEdB9OofiOC\nBdv/PZOdnU1iYiJAonMuO1RfS7jDxUzgK+fc5ds+9oClwEPOuXu3c/97gGOdc/3LXDcVaOmcO66S\nz1EnwsX++1szmKee8rsSqYo//oAhQ6B5c/jiCxuGnTXLjk0fN85ChohErm+/td0lr7wCbTrn0mxU\nKvVaBOnUsmpvCHILckmdnkqwIEj7ZgHuHJDOioWxzJtnJ77Om2cnLhcUAGPioU3On49tVBBHQsYC\nmjSxZoqfD4xnfePS2xsWxNH+5QWsXg2F55V/rLelMa5BUelz1WtEw/oN2VJswzHdi7rz4x0/QojD\nRdimRTzPawgkAjNKrnOWZD4ChlbysAO23V7W+zu4f52haZHotOeets9+4UJr5LN+vf2ZkIB2hohE\ngYED4eWX7YW/wRmpLHGZLMq37eaJ96bw8cc2mrAjubmQ/LAt4M5ZncOXyzI57LEU0tKss++SJfZ5\nxo2Dt96ChJ7lp0gH9w7w5ZfwySe2FmRg9/K3D+kdYOlSCyZD+5e/rV6FowM2FW+icHMhRVuLKNpa\nxI95P1b7e7Mj4TwKqS1QH8itcH0u0KuSx3So5P57eJ7X2DlXtJ3H1AkKF9GrTx948UVbid68uf07\nfvutDiITiSbx8dCsfRBWl163fF2Qww+3n+mkJDj8cLskJtobitdft8uXX4K7LAhtSh/bcZ8g36+0\nNyAVDSlMJ2VaSrmpjLLS0yq//fXTyt+2uXgzs5bNCvW3Y6dqza+3sWPH0rJly3LXjRw5kpEjR/pU\nUWjpbJHodsIJdg7Jhx/C4MHQu7ffFYlIVQViAuSsLp1yOKBvgMe+gxkzbDvrTTfZpaz4eDt88PFN\nAWb9XvrYvdsGthssoHRxeGV2dHvF2/IK8xj+0nC+Wf4N7gcHP1R4wE5GXaornOFiJbAVaF/h+vbA\n75U85vdK7p+/s1GLiRMn1uo1Fxq5iG7r1tmwKkBWlrUj7tfP35pEpGoqjhi8fHI6P86yn+3Zs7f/\nmAUL4Prr4cAj0+neJ4WiRkG6tP7raES4xDaPpWG9Rjgc9MMuZS0H/hP6zxu2cOGc2+x5XhZwOPAm\n/Lmg83DgoUoe9iVwbIXrjtp2fZ2mJlrR7aqrYOVKmDMHzjzTRjK+/hratvW7MhHZFSULMpetDdJo\nc4C2X6XTZ2ws+fm2myQtzU5hTUoqnfJcvx4yM+Hjj2HGjFhypmZQXAwN9oabvoTDDoN+B+Ry0Sep\nle7mqE6NJc/1+PDHGf3WRcxa9lUIvxO7Jty7Rf4GPANcROlW1FOAfZxzKzzPmwB0dM6ds+3+3bBB\nm0eByVgQKdmKWnGhZ8nnqBO7RWJi4PbbreGLRJd33oHjj4cnnrATVJcsgUGDbC3GBx/YUdEiErmc\ng/6TkpmTX7rFs/kfSVzXLoOTTrKzS3blPKA1a6zvzYwZFjh++gk4Nxm6lj5vxy1JXNLEtpW2aQOt\nW5e/FDXIJS09leXrgsQ2DfBwcjoNimJZswYu+SaZeetLn6vephYUN1q346JKRy5CulskrGsunHPT\nPc9rC9yGTW/MBo52zq3YdpcOQOcy91/ked7xwERgDPAbcF5lwaIu0bRIdFq1Cs4/3/pdXHCBXdel\nC7z6qi38GjsWHnnE3xpFZPucszcHd90FcwaVX5DZvnuQf46p2vO1amULu086yT4OBmG/yUHyyjTp\nWrEhyINP2O+OrVu38yTnpv4ZRhauyWHI/SkwZdsaizHla3QNN+y0pp579uQXfqnaF7ILwt6h0zn3\nqHOum3OuqXNuqHPumzK3neucO6zC/T93ziVuu38P51ydbxtVXGwNmRQuos8//mHd/p56qvw7m4MO\nslDx73/Df8Iw3yki1bd1K0ybZttDhw+36/btEvqOuoEA9AhU2Fa6b4C8PNi8GfLz7ZC12bNtG2p6\nOrSNK39MQKBXkG++gV9/hSF9yj9XTOPyHb9aNGpBXOs4kjonkXt1Lm6cY+opU3f769ieWrNbpDYr\n2raUVeEiurzyirUOfvFF6NTpr7ePHg3ffQeXXmq7Rw46qOZrFJFSmzbBCy9YZ8758+Goo2wa4+CD\nYcX6HW8Pra7KtpV6HrRoYZcuXUrv/8DqACuXlu46iWsXwBpswlsdyj/XE8Of4MK3LwzJeo6qCuua\ni5pQF9ZcrFljc23Tp8Opp/pdjeyK33+3edhhwyxkVDYfu3mz/QL78Udb4Nm1a42WKSLYwsunnoL7\n7oPffoOUFGtyN2iQ35X9VV5h3l/CyO4EhnC1/9bIRRTQyEV0cQ4uvBDq14fHHtvxQq+GDS18DB5s\n87AZGdZoS0TCb+1aePRRmDjR1jicfrptG913X78rq9zOemBECh25HgUULqLLs89ay+8nn4R27XZ+\n/7Zt7f7z59vJjFE+mCgS8VasgJtvtpHCW2+F1FT7+XvuucgOFtFEIxdRQOEieixZApdfbueHnHDC\nrj+uXz878TYlBQYM+GuXPxHZfb/9Bvffb4uo69WDiy+GK6+0hZUSWgoXUUDhIjoUF8OoUdCyJTz4\nYNUff/LJ9i7q5pstbFQlnIhI5X791U42ffZZ6xl07bVw2WXbP9dDQkPhIgqUnLins0Ui22OPWXOc\nDz+0gFEd//yntQY/4wyYOdMabYlI9Xz/PUyYYIvhY2OtX8WFF9oODAkvrbmIAhq5iHzz58M119i2\n0iOOqP7z1KsHzzwDcXHWSnjVqpCVKFJnzJxpI38DBtjfH3nETim9+moFi5qicBEFFC4i29attsai\nUycbet1dMTF2TPOaNfC3v1kDNRHZMeds5PCww2DoUJsKee45O1Ts4os18lvTFC6igMJFZLv/fvjq\nK5vPDdU20r33hv/+Fz77zA49E5HtKy6GN96AAw6wUcP8fGuvP2cOnHWWzu7xi8JFFFC4iFw//AC3\n3GLDrQceGNrnHjbMFoY+9BBMnhza5xaJdlu2wEsv2dTHSSfZyMT771szupQUm2IU/+jbHwUULiLT\npk1w9tnQowfcdlt4PsfFF1ub8Isugi++CM/nEIkmRUW2lbRXL1v43KUL/N//2SjfUUft2umkEn7a\nLRIFFC4i0+2329DrrFnh+7fxPHj4YZg7196Nff01dO6888eJ1DaFhRYq7r/fThM95RSbOhw40O/K\nZHs0chEFSsJFo0b+1iGlZs2yLW633BL+X26NGtkv0caNrRfGhp2foixSa6xebUG+a1frT3H00Ra2\np09XsIhkGrmIAkVF9gKj4T7/5RbkcvLLqXzzc5CmlwQYdVk6EP5TBmNjbdFaUhKcd56dtKr/D1Ib\n5Rbkkjo9ld/WBtmyOsCaJ9LZmh/L+efbdu+yJ4RK5NLIRRQoKtKUSKRInZ7Kl8sy2RyTQ8GemaSl\np9TY595vP+uBMXUq3HtvjX1akRp17LOpZC7NZHF+DsvqZ9LqwhQWLbLpQQWL6KFwEeFyC3J5qCCZ\nwvPjSZ6cTF5hnt8l1VnOwU9Lg+WuCxYEK7l3eJx6qrUHv+EG+N//avRTi4TNxo02GnfwwfDtL+V/\nphq3DdK+vU+FSbUpXES41OmpLCGT4pY5ZC7NJGVazb1TllLFxXYg2erF5U84CsTU/IlH48db98HT\nT7e5Z5FoNXcujB1rDejOPNN6UvTq5P/PmOw+rbmIcBXfGf+8PMiWLdBA/3I1ZtMmOwp96lS455F0\n3mySQrAgSCAmQHpaeo3XU6+enaA6eFguCY+k0qF7kE4trZbY5uFf/yGyOzZutCZXTzxhW0jbtrV1\nROefDz17Ql5hOinT/P0Zk92nl6gIF4gJkLM658+PV+YE2G8/m3M/9lgt6gu3wkLb8jZjBkybBqee\nGsu1ZPhdFi1aQMx5qWxckcmifFiUn0PKtBQyRvlfm8j2zJ1rW0mfe87OzDnsMHj5ZWuAVXZNWWzz\nWP0/rgU0LRLh0tPSSeqcRFzrOJI6J/H+eem0awfHH2+tbr/91u8Ka69Vq+DII+3d1Tvv2HqHSLJ6\nS/lRrXnLghQX+1SMyHZs3AgvvGBrKfbd1/5+3nnw888W2NPStFi9tlK4iHAlKX7BmAVkjMrgqKRY\nPv4Y3noLli+HxETrErlkid+V1i7LltkvxF9+gU8+2b2TTsOl4lz0HwsDHHus/b8Q8VPJWoqOHUvP\n93j5ZfjtNxt17dnT7wol3BQuopDnwfDhdq7FY49ZP/2ePW0Hwdq1flcX/ebPt34Sa9dCRgYMHux3\nRdtXcVRrWmo6P/wA/fvDa6/5XZ3UNdsbpTj/fI1S1FUKF1GsQQO48EI7Wvjaa+2Qq+7d4ZFHYPNm\nv6uLTt9+a8GiaVM7y2OfffyuqHIVR7X+dnws339vv9xTUmz4uaDA7yqlttMohWyPwkUt0KKFHZw1\nf75tURwzBvr0sXevzvldXfT49FM45BA77vz//i86z/Bo29ZW4j/9tC1A3W8/mDnT76qkttmwQaMU\nsmMKF7VIp072ojJ7NsTF2bvXgw6Cr77yu7LI9/rrcMwxcMAB9suxbVu/K6o+z4NRo+C776BdO0hO\nhltvtSOqRXbHTz+V9qXQKIXsiMJFLdS/P7z3HnzwAaxbZy+YaWmQk7Pzx9ZFkydDaqqN+rz1FsTE\n+F1RaMTH2wjMP/8Jd9xhIePXX/2uSqJNySjFQQfZiKhGKWRXKFzUYkceCdnZMGUKZGba+oErr7Qt\nlmLuvdfWJowebU2yatsvyQYNYNw4CxkrV9o0yeTJmi6Tnas4StGokUYpZNcpXNRy9evD3/9uWyrH\njYMnn7R3tPffb6u76yrn7ITF666zd/aPPmrfq9pq6FBbrHraaRamTjkF/vjD76ok0lQ2SvHLLxql\nkKpRuKgjmjWDm26yYfGRI+H666F3b3u3XtcaL23ZYi+w999vO2xuu61udDpt0QKeesoWfH76KfTr\nZ1NnIj/9BFdcUfkoRY8eflco0Ubhoo5p397epc+ZY2szTj/d1mR8/rnfldWMDRvsXfvzz9u7sjFj\n/K6o5qWkWI+Ufv3g6KPtQLYNG/yuSmrahg32c1AySvHiixqlkNBRuKij9tkH3njD3sGCbcE88USY\nN8/XssJq7Vo7j+WDD+DNN+GMM/yuyD8dO8K779rIzRNPWKOw777zuyqpCWVHKc4+W6MUEh4KF3Xc\nIYdYH4SXXrIXl7594ZJLIC/P78pCKzcXhg2zr/Gjjyxk1HX16tnIzTff2HqTIUPggQfq3jRZXaBR\nCqlpChdCvXq2DmPePLj7bluH0b073HknrF/vd3W7b+FC24aZm2vTPwce6HdFkaVvX5g1y4LG1Vfb\nLqPffvO7KgkFjVKIXxQu5E9NmtiLy6+/2oLH8eNtu9kzz8DWrX5XVzW5BbkkT06m8/3x9LonmS2N\n8/jiC1tnIH/VuDHcd5+9i/35Z/s+TZ/ud1VSHRqlkEigcCF/seeeMHGinRmQlATnngsJCfDhh35X\ntmNr1kBWFrzyCuz/r1Qyl2byW2EOmwOZxF6WQrduflcY+Q47DL7/Ho46yl6Ezj4b8vP9rkp2pGyQ\n7nhzMoHueRqlEN818LsAiVzx8XY+xRVX2IjGUUfZ7oJ777WdJjVtyxb7RZmTAwsW2J9l/756del9\nvSuC0Kr045VFwZovOEq1aWMvSscfD//4hzXgev55m1qSyBAMWoO87GyYlJ/KqphMu6FhDp3OS+Hr\nszIUJsRXCheyU0OH2tHjr71mTaf2289GM267zeZyQyk/vzQ0VAwRixaVno9Rr54dLBYXZ/Wkptrf\n4+IsFJ3weoDMpaX9zgMxgdAWWst5no1aHHSQ9T045BC44QZrxNawod/V1R3OwdKlpUEiO9tG537/\n3W5v3Ro2jC4fnBu3DSpYiO8ULmSXeJ71Rxg+3LYujh9vCz+vvto6XbZosWvPU1wMy5ZVPvqwcmXp\nfZs3t6AQF2fbZEuCQ1wcdO1qw76VSU9LJ2VaCsGCIIGYAOlp6bv3Daij9t4bPvsM7rnHgsX779sc\nvlo/h55z9jOQlVU+TJR0Uo2NhcREWw+VkGCXrl3hoCkK0hJ5PBflhwx4npcAZGVlZZGQkOB3OXXG\n2rUwYQJMmgQtArm0vCCV4uZBOrYI8PzwdApyY7c7+rBwIWzaVPo8nTqVBoay4SEuzk70rAudM6PF\n11/DmWfa1NTEiXDBBfr3qa6tW22BZdkQ8e239nMFsNdepQEiMdH+DAS2//3OK8z7S5CObR5bs1+Q\nRK3s7GwSExMBEp1z2aF6XoUL2S1LlsCgR5NZ0TSz9MrFSTAlA4CmTbcfHOLjoVs326Ei0aOwEK66\nykavRoywduKxeh3boc2bbXF0yZRGdjbMnl26zTsurjRIJCTAwIH6nkrNCVe4CNu0iOd5rYFHgOFA\nMfAqcLlzrnAHj5kCnFPh6vecc8eFq07ZPV26QIuOQVaUWUzZLj7Iaxn2S7NDB727rU2aN4fHH7fF\nnuedZ1tWp0yB4/QTCthhgHPmlB+R+P57KCqyn4OePS1AnHxyaZBo3drvqkVCL5xrLl4C2gOHA42A\nZ4AngDN38rh3gb8DJS9JReEpT0IlEBMgZ3XpnG/PQICkJB8LkrAbMcLOJxk1yoLGJZdYn4xmzfyu\nrOYUFlpwKLvQ8scfbdFx/fqw774WIM480/4cMGDX1yaJRLuwhAvP8/YBjsaGWb7ddt1lwP88z7va\nOff7Dh5e5JxbEY66JDy0eLJuat8e3n4bHnvMpko+/tgWe9bG2cn8fFsTUXZEYt48W6DcsKGN4Awe\nDBdeaF9///42JShSV4Vr5GIosLokWGzzEeCA/YE3dvDYYZ7n5QKrgY+Bm51zq8JUp4RAbPNYMkZl\n+F2G+MDzbNTi0EPtHfoBB8Dtt9suovr1/a6uev74469BYv58u61JE9v6PGwYXHmlBYk+fXa8c0mk\nLgpXuOgAlDv6yjm31fO8Vdtuq8y72NqMhUA8MAF4x/O8oS7aV56K1GK9e8OXX9p21RtugHfegeee\ns62SkSw3t/xCy+xsWLzYbouJsTURxx1Xuthyn32ggTbwi+xUlX5MPM+bAFy3g7s4oHd1i3HOlT3N\n4EfP834AFgDDgE929NixY8fSsmXLcteNHDmSkSNHVrccEamCRo1se/Kxx1rjrQED4NFH4fTTa76W\n3IJcUqen/jlV9+rf0tm0OrbcaER2Nixfbvdv3drCw9/+Vhokune3Zm0itcXUqVOZOnVquevWlux/\nDrEqbUX1PG9PYM+d3C0HOAu43zn35309z6sPbAROcc7taFqk4ufMA25yzj1Zye3aiioSYdassdbh\nL75oJ+4++ii0arXzx+2O4mKb0ggG4fSPkvlxXen26AbLk9jyH5u6a9eutHdEyaVbN+1qkropIrai\nOuf+AP7Y2f08z/sSaOV53sAy6y4Ox3aAfLWrn8/zvL2wMKODIUSiSKtW8MILtpPk4ottgePzz1sb\n8aravNmmL4LB0svvv5f/uOS6kvbwjAlCm9LnaNExyLNvWpDo2FFBQiTcwjJ76Jyb53ne+8CTnudd\njG1FfRiYWnaniOd584DrnHNveJ7XHBiHrbn4HegO3AP8ArwfjjpFJLxGjrSTdc8+2xZ9XnptLlk9\nUsktDNK+WYBJQ226YntBoeTvK1daa+wS9epZk6lAwC79+9uBeiUfBwLwj28DZOWVbo/et3OAESN8\n+AaI1FHhXJp0OtZE6yOsidZ/gcsr3KcHULJQYivQHzgbO89yORYqbnHObQ5jnSISIs7ZlEjFsLDf\nfnZGySN5qbCtm2vO6hz2fyDlz26ujRtbMOjQwf5MTi4fGEou7drtfCfKO/20PVrET2ELF865Neyk\nYZZzrn6Zv28EjglXPSJSfVu3Ql5e5VMSZUcciiq0vdtjDwsFw4ZB5p5Byr5TaN8jyMc/2u2tWoVu\nukLbo0X8pU1VInVQyW6K5flB2jQOML53OkWrYisNDnl5tmCyrHbtSkcT9tnHpj1KRh3KXsp27Uye\nXP4EzzVLAyxebN0sRaT2ULgQqSOKiuC772DmTLh9WSorm9n0xEJyGP61TU80aFA+IAwZUj4olNzW\nvr11pqyqst1cW9UP0HxmOscdB6mpdtJq584h/qJFxBcKFyK1kHPWDGrmTPjqK/szO9uOu2/UCLwr\nym/A2mvfINl5sOee4e3tUHG6wl0C06bB2LHWiGv8eBgzpnrBRUQih1rEiNQC69bZ2R4TJsCJJ9oI\nw957226NN9+0E2rvu8+CRn4+DOoVKPf4rm0CtGtX802jPA9OO83O6TjvPLj2WutBkZm588eKSOTS\nyIVIlNm61V6MZ84sHZmYM8dGK1q0sKmMCy6wcz6GDLFtmxVF2mFzLVvCgw/COedYX4zkZDj3XLj3\nXmjb1tfSRKQaFC5EIlxeXunUxldfwaxZNlJRr54dmnXAAXD55fbnPvvs2oFhkbqbIiEBvvgCnnoK\nrr8e3ngD7rnHjnZXK26R6KFwIRJBiopg9uzyayUWLrTb2re3AHHDDfbnoEE2UlHb1K9vR5effLJN\nk1xwAUyebEe7Dxjgd3UisisULkRqyPYO01q/IvbPEDFzph31XbLoMiHB1k8ccADsv7+dMFqX2lbH\nxsIzz9ioxcUX21qMMWNs0WdtDFUitYnChUgNSZ2eSubS0u6Una9KYfO2w7Ti4ixEnHGGBYkBA6xj\npcDBB1vomjTJgsW0afb3U06pW2FLJJooXIjUkIUry2//bB4I8vxbFibatfOpqCjRqJFNkaSl2fqS\nv/3NzhN55BE7Gl1EIouWSInUgMJCWLW4/PbPPl0CDB+uYFEVXbvC66/b9tp586BvXxvN2LjR78pE\npCyFC5EacO21wLR0EtomEdc6jqTOSb5v/4xmI0bATz/BlVfCnXdCv37wwQd+VyUiJRQuRMLsvffg\n0Ufhgdtiybo0gwVjFpAxKoPY5ttpQCG7rFkzuOsua2m+1142TZKWBsuW+V2ZiChciITRH3/Yboej\nj7YdDxJ6vXtbd9IXXoBPP7WPJ02CLVv8rkyk7lK4EAkT5yxQbNxofRq0syF8PM922vz8M5x1lk2X\nDB5s23tFpOYpXIiEyUsvwSuvwOOPQ8eOfldTN7RqBf/+tzUgq18fhg6F0aNh1Sq/KxOpWxQuRMJg\n6VK49FI4/XTbNik1a/BgCxiPPALTp0OvXjBlChQX+12ZSN2gcCESYsXFduhWixb24ib+qF/fAt68\nebbmZdQoOOQQO+RNRMJL4UIkxB5+GGbMsNbVrVv7XY106GCLPWfMgBUrYL/94JproKDA78pEai+F\nC5EQmjvXTvMcMwYOP9zvaqSsww6zbau33WYjSr0Scul9XzLxD8aTPDmZvMI8v0sUqTUULkRCZNMm\nOPNM6NYN7r7b72pkexo3hhtvtAZcG05IZd76THLW5JC5NJOUaSl+lydSa+hsEZEQuf12+P572/7Y\ntKnf1ciO7L03tO4cZPWa0uuCBcHKHyAiVaKRC5EQmDnTukXecosdDS6RL9Ci/FkvgZhAJfcUkapS\nuBDZTYWF1rhp8GC44Qa/q5FdlZ6WTlLnJAJN4mBxEmnorBeRUNG0iMhuuvpqWL4c3nkHGugnKmrE\nNo8lY1QGYN09b7sOzjgJ2rTxuTCRWkAjFyK74d13rQPn/fdDjx5+VyPV9cADtiBXI08ioaFwIVJN\nJYeSHXMMXHSR39XI7ujQwY5u/89/dB6JSCgoXIhUg3MWKDZt0qFktcXFF9ti3Isu0omqIrtL4UKk\nCnILckmenEz7u+L5b4tk7nkkj4A2GdQK9evbFNf339vhZyJSfQoXIlWQOj2VzKWZrNiSA10zeWaD\nGi/VJoMG2QjGP/9pi3RFpHoULkSqoGKjpW9+DvJ//+dTMRIWd95pTdDGjvW7EpHopXAhUgUVGy01\n2BDg4INh+HD44QefipKQatXKdo9Mnw4ffOB3NSLRSeFCpApKGi/FtY4jqXMSv96ZztSpdqz3gAFw\n9pNyp1EAAA55SURBVNmwaJHfVcruOuMMOPRQO7J940a/qxGJPgoXIlVQ0nhpwZgFZIzKoEOLWE47\nzU5D/fe/7Z1uz55wxRV2vLdEJ8+zf8/Fi+Gee/yuRiT6KFyIhEDDhrYQ8NdfYdw4mDIF4uJg/HhY\nt87v6qQ6eve27qsTJti/q4jsOoULkRCKiYGbboIFC2D0aDvMLD4eHn7YemJIdLn5ZggEbHrEOb+r\nEYkeChciYdC2rS0KnD8fjj/epkn22QdefBGKi/2uTnZVs2YWDD/4AP77X7+rEYkeChciYdSli02R\nfP899OsHZ54JCQl2JoneCUeH4cPhpJMsIObn+12NSHRQuBCpAX36wBtvQEYGtGgBxx1nuxF0jkV0\nePBBWL0pl973JhP/UDzJk5PJK8zzuyyRiKVwIVKDkpLg88/hrbfs4LOhQyElxbaySuTq0gXaXZbK\n8oaZ5KzOIXNpJinT1J1VpDJhCxee593oeV6m53mFnuetqsLjbvM8b7nnees9z/vQ87zu4apRxA+e\nZ0Pts2fDs89CdraNbJx/Pvz2m9/Vyfbk5EAwv3x31ordWkWkVDhHLhoC04HHdvUBnuddB/wDGA0M\nAQqB9z3PaxSWCkV8VL++Nd36+Wdb/Pn669CjB1x7Laza5Tgu4TZrFhxwANTfUL47a8VurSJSKmzh\nwjk33jn3IFCVpsiXA7c75952zs0BzgY6AieFo0aRSNC4sS0WzMmBa66BRx+17at33w3r1/tdXd32\n5pswbBh07w7Z15Xvzpqelu53eSIRK2LWXHietzfQAZhRcp1zLh/4ChjqV10iNWWPPeC226xHxhln\n2MmcPXrAf/4DW7b4XV3d8+ijcPLJcMwxMGMG9O5SvjtrbPNYv0sUiVgREy6wYOGA3ArX5267TaRO\naN8eHnnEFnkecghceKGtyfjvf7V9tSYUF8N111njrMsug1desVNSRWTXVSlceJ43wfO84h1ctnqe\n1zNcxYrUJfHx8NJLtuBz773h1FNh//3h44/9rqz22rgRTj8d7rsPJk6ESZNsbYyIVE2DKt7/fmDK\nTu6TU81afgc8oD3lRy/aA9/u7MFjx46lZcuW5a4bOXIkI0eOrGY5IpFh4EB47z345BO4/no4/HA4\n6ig78yIhwe/qao9Vq2waZNYsG61ITfW7IpHQmjp1KlOnTi133dq1a8PyuTwX5nFWz/POASY659rs\nwn2XA/c55yZu+3gPLGic7Zx7pZLHJABZWVlZJOg3rdRyzsFrr8GNN9ouk7Q0uOMOW3Ao1bdoERx7\nrJ1k++abcOCBflckUjOys7NJTEwESHTOZYfqecPZ56Kz53kDgK5Afc/zBmy7NC9zn3me551Y5mGT\ngJs9zxvheV4/4DngN+CNcNUpEk08z5puzZkDTz5pHT9797b1Ab//7nd10Skry7aabtoEX3yhYCES\nCuFc0HkbkA2MA2K2/T0bSCxznx7An3MZzrl7gYeBJ7BdIk2BY51zOk9SpIwGDazp1vz5cOedtjYj\nPt5O8QzTKGet9M47tmi2a1f48kvoqRVjIiERzj4X5zrn6m/n8nmZ+9R3zj1X4XG3Ouc6OueaOeeO\nds79Gq4aRaJd06bWdCsnx3Y2PPCAhYx//csWJ0rl/vMfOOEEOOIIW88Sq52lIiETSVtRRaSaWre2\nplu//mrTJtdeC716wTPPwNatflcXWZyDm26yLb4XXQSvvmpHq4tI6ChciNQinTrZO/Iff4TBg+Hc\nc2HAAFukWNd7ZOQW5HLg08nscXM8dwWTGXdvHg8/rK2mIuFQ1a2oIhIFevWypluzZtn21RNPtBNZ\n774bkpP9ri60nLM26atX22XNmtK/l7283CyVFU0zoRHQNYeP9kzhVi/D7/JFaiWFC5FabMgQa139\nwQcWMg46yE5kvesu6NfP7+pKOQfr1pUPA5WFhO3dtnnz9p+3WTObMmrdGtYdp1NNRWqKwoVILed5\ncPTRcOSRMG2a7SgZMADOOsvOMunaNTSfp7jYdqrsShioeP2aNZWvDdljD2jVqjQktG4NHTuW/r3i\nbWWvb1TmPOXkyQEyl5b2+NOppiLho3AhUkfUqwcjR1rnySeftGDx8stwzqW5fNcrlZVFQTo0D/DU\nkek0KIqtUkhYvRry87e/rsPz7IW+Ygjo1u2vYaBiQGjZ0rbdhkJ6Wjop01IIFgQJxAR0qqlIGIW9\nQ2e4qUOnSPUUFNjZGbcuTmbrXpmlNyxOginl1yLUr7/jUYIdBYQ99rBgIyKRJ1wdOjVyIVJHxcTY\nFMlTE4Mszi+9vn33IC9/Uj4gxMTYCISIyK5QuBCp4/ZqGWBxfulahO4dAgwb5l89IhL9FC5E6jit\nRRCRUFO4EKnjYpvHkjFK/R5EJHS0zEpERERCSuFCREREQkrhQkREREJK4UJERERCSuFCREREQkrh\nQkREREJK4UJERERCSuFCREREQkrhQkREREJK4UJERERCSuFCREREQkrhQkREREJK4UJERERCSuFC\nREREQkrhQkREREJK4UJERERCSuFCREREQkrhQkREREJK4UJERERCSuFCREREQkrhQkREREJK4UJE\nRERCSuFCREREQkrhQkREREJK4UJERERCSuFCREREQkrhQkREREJK4UJERERCSuFCREREQkrhQqpl\n6tSpfpdQ5+h7XvP0Pa95+p7XDmELF57n3eh5XqbneYWe563axcdM8TyvuMLlnXDVKNWnXwA1T9/z\nmqfvec3T97x2aBDG524ITAe+BEZV4XHvAn8HvG0fF4W2LBEREQmnsIUL59x4AM/zzqniQ4uccyvC\nUJKIiIjUgEhcczHM87xcz/PmeZ73qOd5bfwuSERERHZdOKdFquNd4FVgIRAPTADe8TxvqHPOVfKY\nJgBz586tmQoFgLVr15Kdne13GXWKvuc1T9/zmqfvec0q89rZJJTP61X+mr2dO3veBOC6HdzFAb2d\nc7+Uecw58P/t3VuoVFUcx/HvD7GiA12tc4JEI6Mi8ihmdiG1ThcssjAoSrK3Agush+ytopfIii4Q\nEUUSlUFPId3sYtFDaSRaUVlSCZEeS5MKLyX672EtbTzOnDMz7jlz9pzfBxbM7Fl7ztp//sz8Z5+1\n9+KJiGj4DISk04Afgb6I+KhGn1uAVxt9bzMzMztgfkQsK+rNGj1z8RiwdIg+PzU5lkNExM+StgKT\ngKrFBbACmA9sBHYX9bfNzMxGgaOAiaTv0sI0VFxExDZgW5EDGIykU4ETgc1DjKmwasvMzGyU+bTo\nN2zlfS7GS+oFJgBjJPXm1lXRZ72k6/LjLklLJM2QNEFSH/AG8AMFV1RmZmbWOq2c0PkQsKDi+f4Z\nOpcCn+THZwDH5sd7gcl5n+OATaSi4v6I2NPCcZqZmVmBGprQaWZmZjaUkXifCzMzMysxFxdmZmZW\nqFIWF14Ubfg1E/O830OSNknaKel9SZNaOc5OIul4Sa9K+lPSdkkvVE6IrrGP87wBku6U9LOkXZJW\nSZo+RP/ZktZI2i3phyaWNxj1Gom5pFlV8nmvpJOHc8xlJukSScsl/ZrjN7eOfQ47z0tZXPD/omjP\nNrjfO0A30JPbzQWPq5M1HHNJ9wF3AbcD5wM7gBWSjmjJCDvPMuBsoA+4BpgJPFfHfs7zOki6CXgc\neACYCnxJys9xNfpPBN4EPgR6gaeAFyRdMRzj7QSNxjwL0uT//fl8SkT81uqxdpAuYB2wkBTLQRWV\n56We0NnI3T8lLQWOjYh5rR9Z52ow5puARyPiifz8GGALcFtEvN7akZabpLOAb4FpEbE2b7sKeAs4\nNSL6a+znPK+TpFXA6ohYlJ8L+AV4OiKWVOn/CDAnIiZXbHuNFO+rh2nYpdZEzGcBK4HjI+KvYR1s\nB5K0D7g+IpYP0qeQPC/rmYtmeVG0YZJv3d5Dqn4ByB8Oq4EL2zWuErkQ2L6/sMg+IP3ymDHEvs7z\nIUgaC0zj4PwMUoxr5ecF+fVKKwbpbxWajDmAgHX536vvSbqotSMd9QrJ89FUXLxDuofGZcBiYBZp\nUTS1dVSdq4f0RbhlwPYt+TUbXA9w0KnfiNgL/MHg8XOe12ccMIbG8rOnRv9jJB1Z7PA6UjMx3wzc\nAdwAzCOd5fhY0pRWDdKKyfMRsyqqmlgUrREDTsN/I+lr0qJos6m9bklHa3XM7VD1xrzZ93eeWyfJ\nnz2Vnz+rJJ0O3AN4Mu0INmKKC0bmomidrpUx7yedzuzm4Cq4G1hbdY/Rod6Y9wMHzYiXNAY4Ib9W\nF+d5TVtJdwXuHrC9m9rx7a/R/6+I+KfY4XWkZmJezefAxUUNyg5RSJ6PmOJiJC6K1ulaGfP8pdZP\nutLhKzgwoXMG8Ewr/mYZ1BtzSZ8Bx0maWjHvoo9UsK2u9+85z6uLiD2S1pBiuhwOTC7sA56usdtn\nwJwB267M220ITca8mik4n1upmDyPiNI1YDzpEpn7gT/z416gq6LPeuC6/LgLWEL6YptASuYvgO+A\nse0+njK0RmOeny8mfZFeC5xLWohuA3BEu4+nDA14O+fpdNIvte+Blwf0cZ43H98bgZ2kOSpnkS7z\n3QaclF9/GHipov9E4G/gEeBM0qV9/wKXt/tYytKaiPkiYC5wOnAO8CSwB5jd7mMpS8ufC72komwf\ncHd+Pr5GzAvJ87YfeJPBWko6vTawzazosxdYkB8fBbxLOt2zm3Ta+dn9Ce1WfMwrtj1IWoRuJ2nG\n8aR2H0tZGmkBv1dIxdx24Hng6AF9nOeHF+OFwEZgF+mX2XkVry0FVg7oPxNYk/tvAG5t9zGUrTUS\nc+DeHOcdwO+kK01mDveYy9xIk7r3VfnsfrFazPO2w87zUt/nwszMzEae0XQpqpmZmQ0DFxdmZmZW\nKBcXZmZmVigXF2ZmZlYoFxdmZmZWKBcXZmZmVigXF2ZmZlYoFxdmZmZWKBcXZmZmVigXF2ZmZlYo\nFxdmZmZWqP8A9WIckNBeej4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a0bc240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ys[:, 0], ys[:, 1])\n",
    "plt.plot(res.y[0], res.y[1], '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now time:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Python `solve_ivp`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 562 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit res = solve_ivp(fun, [t0, 5000], y0, method='RK45')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Cython `solve_ivp` with Python level callable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 328 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit ts, ys = solve_ivp_cython(fun_python, [t0, 5000], y0, method='RK45')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Cython `solve_ivp` with Cython level callable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 loops, best of 3: 26.5 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit ts, ys = solve_ivp_cython(fun_cython, [t0, 5000], y0, method='RK45')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we have ~2x speed-up by using Cython version with Python callable and ~20x speed-up with Cython callable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}