Simulating OU Process with Cython

Simulating OU Process with cython.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:61ae7c8a36316f95faf22b4a297dca42bea4229351f8f4d74d669157de833109"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic\n",
      "%load_ext memory_profiler\n",
      "%load_ext autoreload\n",
      "\n",
      "\n",
      "%autoreload 2\n",
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Pure Python"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from math import sqrt\n",
      "from numpy.random import randn\n",
      "import numpy as np\n",
      "\n",
      "def sdeevolvepy(x, tout, A, B, dt=.01):\n",
      "    t = 0.0\n",
      "    while t < tout:\n",
      "        dt = min(dt, tout-t)\n",
      "\n",
      "        ax0 = A(x)\n",
      "        bx0 = B(x)\n",
      "        \n",
      "        x += dt * ax0 + sqrt(dt) * bx0 * randn()\n",
      "        t += dt\n",
      "        \n",
      "    return x\n",
      "    \n",
      "def sdepy(x, tout, A, B, dt=.01):\n",
      "    t = tout[0]\n",
      "    \n",
      "    xout = np.empty_like(tout)\n",
      "    dt = min(dt, np.diff(tout).min())\n",
      "    \n",
      "    for i, (t0, tNext) in enumerate(zip(tout[:-1], tout[1:])):\n",
      "        x = sdeevolvepy(x, tNext-t0, A, B, dt=dt)\n",
      "        xout[i] = x\n",
      "        \n",
      "    return tout, xout\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Cython"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython -lgsl\n",
      "#cython: boundscheck=False\n",
      "#cython: wraparound=False\n",
      "##cython: profile=True\n",
      "cimport cython\n",
      "from cython.view cimport array as cvarray\n",
      "cimport numpy as np\n",
      "import numpy as np\n",
      "from libc.math cimport sqrt, fmin\n",
      "\n",
      "\n",
      "\n",
      "cdef extern from \"gsl/gsl_rng.h\":\n",
      "    ctypedef struct gsl_rng_type:\n",
      "        pass\n",
      "    ctypedef struct gsl_rng:\n",
      "        pass\n",
      "    gsl_rng_type *gsl_rng_mt19937\n",
      "    gsl_rng *gsl_rng_alloc(gsl_rng_type * T)\n",
      " \n",
      "cdef extern from \"gsl/gsl_randist.h\":\n",
      "    double gamma \"gsl_ran_gamma\"(gsl_rng * r,double,double)\n",
      "    double gaussian \"gsl_ran_gaussian\"(gsl_rng * r,double)\n",
      " \n",
      "cdef gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937)\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "cdef class RandStream:\n",
      "    cdef gsl_rng * _strm \n",
      "\n",
      "    def __cinit__(self):\n",
      "        self._strm = gsl_rng_alloc(gsl_rng_mt19937)\n",
      "\n",
      "    def __call__(self, double sig):\n",
      "        return gaussian(self._strm, sig)\n",
      "    \n",
      "    cdef double normal(self, double sig):\n",
      "        return gaussian(self._strm, sig)\n",
      "\n",
      "    \n",
      "cdef RandStream rv = RandStream()\n",
      "\n",
      "cdef double A(double x):\n",
      "    return -x\n",
      "cdef double B(double x):\n",
      "    return 1.0\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "cdef double sdeevolve(double x, double tout, double dt ):\n",
      "    cdef double t, ax0, bx0\n",
      "#     cdef np.ndarray[np.float64_t] rands = np.random.randn(int(tout/dt) + 1)\n",
      "    cdef int i = 0\n",
      "    cdef bint dobreak\n",
      "    \n",
      "    t = 0.0\n",
      "    \n",
      "    if tout < 1e-9 : t  = tout +1\n",
      "    \n",
      "    while t < tout:\n",
      "            \n",
      "        if t + dt > tout:\n",
      "            dt = tout - t\n",
      "            dobreak = True\n",
      "        \n",
      "\n",
      "        ax0 = A(x)\n",
      "        bx0 = B(x)\n",
      "        \n",
      "#         x += dt * ax0 + sqrt(dt) * bx0 * rands[i]\n",
      "        x += dt * ax0 + sqrt(dt) * bx0 * rv.normal(sqrt(dt))\n",
      "        \n",
      "        if dobreak:\n",
      "            break\n",
      "\n",
      "        t += dt\n",
      "        i+=1\n",
      "        \n",
      "    return x\n",
      "\n",
      "\n",
      "def sde(double x, np.ndarray[np.float64_t] tout, double dt=.01):\n",
      "    \n",
      "    cdef double t, tNext, t0\n",
      "\n",
      "    \n",
      "    t = tout[0]\n",
      "    \n",
      "    cdef np.ndarray[np.float64_t] xout = np.empty_like(tout)\n",
      "    dt = min(dt, np.diff(tout).min())\n",
      "    \n",
      "    cdef int i\n",
      "\n",
      "    xout[0] = x\n",
      "    for i in xrange(tout.shape[0]-1):\n",
      "        t0 = tout[i]\n",
      "        tNext = tout[i+1]\n",
      "        x= sdeevolve(x, tNext-t0, dt)\n",
      "        xout[i+1] = x\n",
      "        \n",
      "    return tout, xout"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x0 = 0.0\n",
      "tout = np.mgrid[0:100:.1]\n",
      "tout ,xout = sde(x0, tout, dt=.0001)\n",
      "plot(tout, xout)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10c49d9d0>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEACAYAAACgS0HpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl8VOXVx3+HBAgkkEASwEAKKKigKCAqWlvSUhXBghsi\nrdVXqVvL28W+rq1V7KK1rfpa1NK6VGtdWrcXK6KIxg2LAmUR2YJG1ixANiYLCXneP848fe7c2e7M\nnUlm7pzv58PnLnPvzJPLzP3dc56zkFIKgiAIghCNHt09AEEQBCE9EMEQBEEQHCGCIQiCIDhCBEMQ\nBEFwhAiGIAiC4AgRDEEQBMERrgWDiKYR0WYi2kZEN4U55gH/6+uIaILTc4noJ0TUSUQD3Y5TEARB\ncIcrwSCiLAALAUwDMBbAXCIaYztmOoBRSqnRAK4G8LCTc4moFMCZAL5wM0ZBEAQhMbi1ME4BUKGU\nqlRKtQN4FsAs2zEzATwBAEqplQAKiGiIg3PvBXCjy/EJgiAICcKtYAwFsNOyvcu/z8kxJeHOJaJZ\nAHYppda7HJ8gCIKQILJdnu+0rgg5fUMi6gPgVrA7KubzBUEQhOTgVjB2Ayi1bJeCLYVIxwzzH9Mz\nzLlHARgBYB0R6eNXE9EpSqkafTARSREsQRCEOFBKxfcQrpSK+x9YcLaDb/C9AKwFMMZ2zHQAS/zr\nkwH8y+m5/uM+BzAwxH4lMLfffnt3DyFlkGthkGthkGth8N8747rnu7IwlFIdRDQfwOsAsgA8qpTa\nRETX+F9fpJRaQkTTiagCgA/AFZHODfUxbsYoCIIgJAa3LikopV4D8Jpt3yLb9nyn54Y45ki3YxQE\nQRDcI5neHqCsrKy7h5AyyLUwyLUwyLVIDKTStIESEal0HbsgCEJ3QURxT3qLhSEIgiA4QgRDEARB\ncIQIhiAIguAIEQxBEATBESIYguARli8HNm7s7lEIXkaipATBIxABRx4JbN/e3SMRUhmJkhKEDEc/\nO/Xq1b3jELyNCIYgeIDWVl52dnbvOARvI4IhCB5AC0ZNTeTjBMENIhiC4AFaW4FBg4CDB4FDh7p7\nNIJXEcEQBA/Q1gb06QMUFwO1td09GsGriGAIggdobQV692YrQ9xSQrIQwRAED9DWBuTkiGAIyUUE\nQxA8gFgYQlcggiEIacK4ccBrYdqNiYUhdAUiGIKQJnzyCfDGG6FfEwtD6ApEMAQhRfnJT4JrQ3V0\nBB93+DBw9tliYQjJRwRDEFKM1lbg7ruBe+8FFi4MfG3hQmDnTuDaa4HZs4EXXwR8Pn7t0CGgsBDY\nv7/rxyxkBtndPQBBEAJ5+WXgllt4ffduXlrrbH74IbBoEa8//zxQVcXr9fXAgAFAXV3XjVXILMTC\nEIQUY+5cs15TAzQ1AT38v9QzzwQ+/xzo25cnwfv1A1pa+LVRo1gwDhzo+jELmYEIhiCkGCeeaNZ9\nPuNiuvBCYPx44OabgeZm4JVXgPx8IxgPPigWhpBcRDAEIcVoaQEKCni9udkIwqBBwMCB5rghQ4Dq\naj5m/Hg+Z+BAtjAOH+76cQveRwRDEFIMnw8YOdKsHzzI69nZPKmt6d0b6N+f5zlycnhf377A2LHA\nBx907ZiFzEAEQxBSDKtg1NcDjY28npXF8xv33gvs28f7Bg0C1q/nwoOaiROBKVOAv/+9a8cteB8R\nDEFIMXw+oLSU10tLgX/9i9eHDAHy8oAf/9hYGps2AbfdFigYJSW8fPXVrhuzkBmIYAhCCtHezvMP\nxcW8PX48sHIlcMYZwPXXBx+/aRMvtUsKMIKRnQ288AJwzjnJHbOQOYhgCEIK0dwM5OYC06bx9siR\nwNq1wHHHAT17Bh9/zDG8bG83+7Q7q2dPjqRaujS5YxYyBxEMQUghfD4WjJNO4ozvESM4s3vYsNDH\nE/Gyd2+z7/jjeZmdzWG3gpAoRDAEIYXQggGwCGhrQc9phGPCBLNeUgJ885ssOFlZvK+tLfFjFTIP\nKQ0iCCmEdklpRozgpXY9haKhgSfDNUTAvHnAww8DgwfzvqamQCtEEOJBLAxBSCF8Ps6l0IwaBdx5\nJ3DqqeHP6d/flA7RnH46sGKFyRLXobmC4AYRDEFIIawuKYAnrm+7zcxVOKW4mAsWfv45b7sVjJkz\nzXsJmYsIhiCkEHbBcENxMbBtGy/dCsYrr0gioCCCIQgphX0Oww1FRRxuO3KkO8HQdalaWxMzLiF9\ncS0YRDSNiDYT0TYiuinMMQ/4X19HRBOinUtEv/Afu5aIlhNRlBgRQfAG9jkMNxQXsyvrS1/iSe94\naW7mZUNDYsYlpC+uBIOIsgAsBDANwFgAc4lojO2Y6QBGKaVGA7gawMMOzr1HKXWiUmo8gJcB3O5m\nnIKQLiTSJVVUxBVsc3PdWQe6+GF9fWLGJaQvbi2MUwBUKKUqlVLtAJ4FMMt2zEwATwCAUmolgAIi\nGhLpXKWU9XkoD8A+l+MUhLQgkS6p4mKuOdWnjymRHomzz+aeGnbeeYeXIhiCW8EYCmCnZXuXf5+T\nY0oinUtEvyKiHQAuB3C3y3EKQlqQaAujqMi5YLzxBvcIt6M7AO6Tx7aMx23inop+CAAgxqBAQCn1\nUwA/JaKbAdwH4Ar7MXfcccd/1svKylBWVhbrxwhCynD4MPD228DFFyfm/aJZGHV17LKyhuz26hX+\n/bZuTcy4hK6lvLwc5eXlCXkvt4KxG4B1QroUbClEOmaY/5ieDs4FgKcBLAn14VbBEIRkUFXFN1Fr\np7tk8fjjwIcfAv/1X4l5vzPP5JIiK1cGz2G89RYwdSrw7LPAnDlmf6gChwBw1VXA00/zxLfUp0ov\n7A/TCxYsiPu93LqkVgEYTUQjiKgXgDkAFtuOWQzgMgAgoskA6pVS1ZHOJaLRlvNnAfi3y3EKQlyM\nGgV0leG6ejUvsxNUsOdLXwLOOivYwrj+ehaLoiJu52rFbmHs2MGl0++/n4VCT4ALmYmrr6ZSqoOI\n5gN4HUAWgEeVUpuI6Br/64uUUkuIaDoRVQDwwe9aCneu/63vIqJjABwGsB3AdW7GKQjx4vNx8ltX\noHtbOJlviAW7YNx3Hy9nz+a/D+D5CyCw3lRNDTB8OK/37cuvSS5GZuP6WUYp9RqA12z7Ftm25zs9\n17//IrfjEoRE0doKbNkSuQCgW+rqgI0beV3fxBNFnz6BN/qvfQ344Q/ZovH5eA7jm9/k16zWzYcf\n8lJbHTk5IhiZjmR6C4IDdBG/Q4e4qF+iOe88jkK65BJ+8k8kOTmBFobPx1Vsc3NNMt5nnwEDBhix\nUgrYs4fXdSKhCIYg5c0FwQG6tMYLLwDf+hbfUGMh2mSxnkv4wx94biGR2F1SBw9yOfTcXGD3bt63\ndy9w1FHAp58CnZ08ya/7hOvS6Tk50lcj0xELQxAioHtlawujo4OX1pao0fj4Yw5fdfI5/frFNj4n\n2AWjqYk/xyoYBw6wW2r3bqC2lgWuqopfszZ0EgsjsxHBEIQInH46cOyxJmlNL3eFCgC3sXYtC4WT\nY0O1Wk0U9jkMbWH07Rs4tvx83v/ZZ2bfaacBF/lnFMUlJYhgCEIIiICXXuJ8hUGDTPE+/dTtpPrr\nc8/xk7rOc1gSMpvIfF6ysFoYjY38t2iXlFUw+vdnEamo4E5/ixfzfM0vf8mvi0tKkDkMQQjDc8/x\nsrDQVGzVcw1O8hF0FzztvpoxI3ju47rrgA0buAT5N77hfsyh0JPehw+zFaFDZHNzeV/v3iwE/frx\na59+Cpx0komc0ohLShALQxDC8PbbvLRGD9XV8dJJ6Kt9nuOII4KPefVV4IMPWJBOOSX+sUZCu6T0\n2HU2t56b0O1ftYWxejUwZkzw++TksEtOFyMUMg8RDEGwcegQL2tqeNmrl7Ew6ur4xurEwtDumx/9\nCFi+HBg9Ovyxzc2J64NhR7uk9MS9zrXQgjHB36FGC8b69aHHmpMD/Pa3XZf5LqQeIhiCZ9A3+lh4\n8UXTUU5jtx7GjAm0MIYNcyYYra2cJPeLXwBDhgA7d3LIaij0vEIy0ILx+OO8ffzxvNSCcfTRvNSC\nUV3NNajs5OSYuZumJuD115MzXiF1EcEQPEFnZ+w+9poa4MILgS++CNzv85mb97hxfPPUFkZVFU8I\nO7UwvvMdfq9jjuEb9/vvBx6j5zQOHEhegcOcHBa63/yGt3Vvbi0Y2pro3x9YtozXh9qbFICvr548\nf+ghYNq05IxXSF1EMARPoCu8xiIY69bxUvv2NT6fmW/QDY2amzm57dAhYOxYFgylgMrK0O89axbw\nxBMmvyIrC/jqV81narR1c+AAT64ngz59jDC98gpHfQFGMEaO5GVenhlDOAtDW3Ey+Z2ZiGAInuCv\nf+VlLIX7tFBo375GNzE6+mj27/fty/sWL+Z9BQXAU08Bq1bxzXb79sDzOzr4WCAwr2L06OBjddRV\nQwNPrieDHpZf+cSJZl1nchcUcOHDggIOqW1vN69Z0eIHBFe5FTIDEQzBU8Ty5BtOMPR8wiefAM88\nw66i/fuBa69lC6GwkENhddjtJ58Enj92rFm33mQHDDC1mwBueaojlnr1Yisk2RQXm3Ui4NFH+e85\n9lje16dP+PLqugjhcceZaxZpIl/wHiIYQtqzaZOZuI3FwtA9qu1P/fX1fHPv2ZNvnqWlPGEN8NO3\nrgmlBcda/rylJXDbamHk53PV2ylTeHv/fnMDP/lk5+OOh3vvZYvI3iDpyiudJw1qt9bgwTwxDrBF\n8umniRunkNqIYAhpTWcnP9Frf7xTC6O9Hbj5Zs5BeOYZ3nfbbXx+fX1g7ad+/cyNv0cPM+/w2GNs\nFVhvmNu2cVTV3f4u9FbBKCjgHIdVq3jbGko7frzzvzkefvxjTsZzg47wys83pc8B4J//dPe+Qvog\ngiGkNdqiqK8PjOKJxvLlvPze97jg3r59XAKjoiJYMADuXgcAJ5xgrBkAmDSJw1XXr+ft2lp+Aj/z\nTN62uqQKCnjSuLmZ//3jH0Z8Tj/d+d/cXVjDj62hx1Y3W7rS2cmRYURcMVgIjQiGkNZYE+oGDHBu\nYWzYwE2ELruMy3rr5kVEoQVj2DBevvUWZ2RrN86LL/JSWxm1texmOuoo4NJLAzOmre9ZW8v5GZs3\ns4hY+2qnKlow/vGPwP1eEAyfz/T/2LKFl1/+MvDII903plREBENIa/STbn09T047tTAqKoylMHy4\n6afd3MzvZe9dMWQIL3V+xnXXcQZ3SQlw661m3mLfPhaM/HyO3NKuMut7AMC/LV3q7fMKqYp2SWVl\nGaEGgAcfBP7v/7pnTInCmlezZQu7FFesAJ5+uvvGlIpI8UEhrbG6RmKxMKwTzqWlwJo1vK7dRfas\na3ufigcfNOtHHGEsjJUrw89HWEt/pONEsfUa2MNuzzuPRfioo7p2TIlg40ZTxiU7G3jzTf4HcD0x\npZJbTTidEAtDSGvsguHUwtAuLICFY9Mm836h6jpZ5yLsWKOGPvwQmDkz+udv3uxsnKnE/PmRx/3u\nu+mZn3H88cAPfgCceCJw553Br+toOkEEw/O8+y67XLyKVTAGDnRuYVhLcRQXmyd+bWHYBSNUIpvG\nKhgNDZG76737LnD11ekpGD17cokTTY8e3GDpkkt4+8or07dcyP79bFXaLcnBg83chiCC4XmmTAF2\n7OjuUSQPLRjjxnEtJKcWhl0wtNCEszCiCUZNDbsuGhp4HOH4yld4Al0n+512mrPxpiJVVcDSpYFl\n2e1JkOlCWxvPN9kFY+zY8OVfMhERDA9jb9bTXezdm7z3rq3lZV5ecCvScCjFN3g9h6EjoAAe69Kl\n8VkYra381B2tzeqZZ7Kw/fKXPLGarhQXszh+5StmXzJazHYFWjAmTACmT+d9mzdzKfff/c4clyq/\nqe5CBMPDOGkjmmx27uRIonhKjzt9f4Bv6LqzXDRqakyLUsDUVzrtNOBXv+J1u2Bcdln4iJn8fL7W\nffsGR1eFYtIkXiar2GBXM2mSmbuINNeTyjQ2AkVFnGfz6qvA1q3sfrv5Zk60rK/nKLEePQIjxDIN\nEQwPs2+fWe+uJyPdM3rt2uS8v3a3XXZZeAujpsb0ggDYxTBihNnOzuYcg6uuMuGVdsHo3x+YOzf0\nGKwRNJHmL6yfB3CRQq+g/+506/mtfxcHDwaGPesaWb168Wvf/775PaWr2y0RiGB4mG9/26x311OR\nFoxkRZrs3AksWQJcfnl4C+Odd3hCVrdctQsGwE+O1pDQcI2OouG0xMcf/xhegNIRIs5WHzWKI9AW\nLeruETnDmnQYyTp6+mnzXdY9QzIREQyPohTnBGi6yz2lE+IS8eRZWxvs2tq50/RuCGdh6HNuuomz\nlDduDBYMwAjG8OGmemus6G520bjmGu+4pDS33MJW07PPcmXfdECLABC6jznA3+Hx47mEDADMm5f8\ncaUqIhgexd4UqLsE4/33ORopEYIxaBCwYIHZXrOG8ye0YOTkcMbxSSdx6Q9NUxMvP/4YuPhiLskR\nqiy3bpp0//3xZ18nq2teOtC3L1f+TadSIdbCieefH/qYgQN5jkaXDMlkRDA8yt69PGlXXs5+2O4S\njLa22PIjovG//wv8/OdsQT31FDB7tplobm7mUM81awKrqV53XfD7XHVV8L4ePbg21IwZsY9r1Che\nOrUwvEi/fnxTveWW7h6Jcx5/nIMeHnss/DEDB/L36oYbONnTyTyVVxHBSCP0k7IT9u7lJ+YpU4Az\nzuhewejfP3GToT4fWwhNTcDnnwcW7bOGdOrJ6/Z2XuqJ6XPOCSweaOf88+OzLlat4sl13esiE9EV\nfTUdHezyCSXYqcLAgcDChcAVV4Q/pl8/49a88MLIIdZeRwQjjejfP7A5TyT27OFyzfq87hSM/Hz3\ngmF/qquuZvGwJlpdeaVZtxYlBIwL6u9/D5zbSRT5+YHd7DIR3Stc09jILp8//pF7nKcioSoT29EP\nF+PHA7//fXq53BKNCEaaoJ+U33rL2fF79nD+A9D9gpGXxzfqRFJTw1aEtUig1WrQgqHncs44g91Y\n9qKCQuIg4oZUmspKUxJd9zhPNZwIhiY/31gb27alRp5TVyOCkSZod9TWrWZfZyffkNesAe67z+yv\nqWG3jRaM/Pzuix0/dIijS955x937HHUUcOSRZnvvXhYMa/lwgCOgbr/duKTq6ti6eOghd58vOOO/\n/9usL15sHnRSkfff59+FU8Ho0YNF8fjjuTT+N76R3PGlIiIY3UR7O9/YnaKfZioqzL6xYzky6Fvf\nAq6/3rzfY4/xBLAWjNNPd26ZJJq2NjMnEG8uyMMPsyjq/hXz5nF0lM8XbDGMHcu+dKuFMWJE+pas\nSDf0/8e553IXQnu0Xqpw4ACXNGlvd56dri1YnT/z8cfJGVsqI4LRTTz8MNcgsjZuiYS2MKw/QB3m\nt2ULZw/rgnY6S1kLxtFHd1/FzbY2k9Ecj5XT2cltVAFOBtuzh+v7bNwY7JLS5OYGzmFkclRLV6Mt\nvlmzWNStjaJSCV0h4Ktfdd7rokcPc06m4lowiGgaEW0mom1EdFOYYx7wv76OiCZEO5eIfktEm/zH\nv0hEDir0pBf6h6QbtUSjsZGfhKx+0xNOMF/eb34TeOUVdtXocgf6Rtm3b3IzvQ8cCCw9snixmRhs\nazNuCWupEqdYBXXwYI78KinhMMdQLimARcTqktJ9L4TkowV86lQu3vfKK8A99/BNOd7s+WSwcydH\nzMXiKtXCYp/czyRcCQYRZQFYCGAagLEA5hLRGNsx0wGMUkqNBnA1gIcdnPsGgOOUUicC2AogjSK7\nnVFdzTeyWCyMkpLA0NqmJvbNf/wx5wHcfz/3Jti/n/sSHHccH5dswSgsDJzUnjULuPdevkF0dJjq\nn/EIRkMDi8QNNxi3UnExsG4db9trPgGBFoYIRtdy6qn8nbZWAJ44kR92EpWLkwg2bIg9m19bGNZo\nuEyrXuvWwjgFQIVSqlIp1Q7gWQD2ALqZAJ4AAKXUSgAFRDQk0rlKqWVKKf08shLAMHiMgwe5OqbT\nH1FzMxdHs1oY+/dz6OykSSa8NDub90+fbp6IkikY27fz0t4zYP9+nvDu1Qu46y7Ol4hHMBob2VK6\n5x6zb9AgdjXNnGl+xFZEMLoPIr7+1lyW7OzkP7TEyqpVwOTJzo9/+23gT3/idW1FZWVlXqSUW8EY\nCmCnZXuXf5+TY0ocnAsAVwJY4nKcKUc0wejsDHyttZVvlNYvaHOzcclowRgwgN011sqbyfqx/v73\nJsPZXvSvro7dUdoqKCqKXzDsJcN1DaZwyVbikkoN5s0DnnsO+PKXU08w6utjK+NSVmYSE4nYsigs\nzLxChG4Fw6lBFlcLdSL6KYBDSqkwnQjSF5+Pv3DhBGPhwsCM0pYWvvHpUNr2dv7S6ic5/dQzYIDJ\n8tZod0AifchKAf/zP2Zbi4HVRHcrGO3tHOH10UeB+7Oz2XoJl1UtFkZq8MgjXLsrO5u/g2+80d0j\nMlgftuKlpoZL0zjh6acD61alK9kuz98NoNSyXQq2FCIdM8x/TM9I5xLRfwGYDmBquA+/4447/rNe\nVlaGsrKyGIbevUSzMHSPaE1rKwuILoSWmxsoKDqTurMzWDB69DCiEcrnHw+7dvFn6G56OgJKP0W2\ntXHOiP5RFhVx1Ews6PcOJXSRyneIYKQe27Zx/a7vfje5n/PWW1x8Mlojq0QIRizoVgOJnPPYtInr\nxYVyy1opLy9HeXl5Qj7TrWCsAjCaiEYA2ANgDgB7lf/FAOYDeJaIJgOoV0pVE9H+cOcS0TQANwCY\nopQK6+W3Cka6EU0w9JdZh462tPBN/4gj+EcxblygYOgSGD4f37yLigLfT7sEEiUY9fVsIe3dy2Gv\numSJjo5qaAhs3ZmfH3tJBd3v+osvYjvP6pKSsNrUorU1uV35pk4Fbr3VdE4Mh8+XuN9CNHQdqt69\nTde+eNBthfXc5NixXK7/oosin2d/mF5gLfkcI65cUkqpDrAYvA7gUwDPKaU2EdE1RHSN/5glAD4j\nogoAiwB8L9K5/rf+A4A8AMuI6N9E5Kk8XaX4hhbJJaVbXuoa/C0tLBBDhgCXXsp1baxf+Guv5SZC\nTU2hE9pGjgT+9a/E/Q06pLWujucStLtJi4I9XFi3MY0FnTtiL2oXjV69+If55JMcviwWRuqggySS\nya9/bUqShCMRFsaGDXzTjsbOndxj5cgj3bmlBg8GPvggcF8sBUkTges8DKXUa0qpY5RSo5RSd/n3\nLVJKLbIcM9//+olKqTWRzvXvH62UGq6UmuD/9z2340wldPZzXl5owViyhCeUAdPgRbukdHkMpQIt\njMJCFow9e/iHkJUV+J5f/zonuyUKbfkUFLA1U10NvPwyC0aoJ7f+/TkmPxbRGD8+vklFIi4Hcvnl\nvJ3JPSpSjauvTt5NzioS//hH+OOamjgwxK2F0acPP8itXh3Z1fTFFywY554LfPppfJ+lg0r036iv\nodOw/EQhmd5dTEcHf9H69g0fm677X3/5y+YpW7ukLrnEHGcvszx4MM8bhPLf5uQktt+y1YopKuJx\nnn8+96EYFiIIWo/pwQedf0ZLS+j3coJ+8vvOd8TCSAX+9S8uALliRWCvkkRy9928vPHG4An2/fv5\nt3HNNcBZZ/E+txaGFoxJk4ClS8MfpwUjL8/MrcWKDlvXwlFVFbjsKkQwuhh9025uDhaM2bOB55/n\n1372M37CfvNNfmrXFoYu9wEEh7KWlPC5oQSjd+/EJk5Zs6ytPzw9EQdwo6N163hdh/3eeivXGHJC\nYyNbJvHwwAPACy+wW0rofk49FbjgAl5P1tzBz35mPsvaQ76lhR9qxo3jXArtms12OYPbpw+3DQYi\ni+BPf8pzD7m58VsE+u+xlwiy/p1dgQhGF6MFo62Nv0BW8/z554EnnmA31IgR/PqTT7LbZ+9e/qEN\nHmyOt5v2Wijs7igg8Zm21jpO1lo8u3aZ0gnFxVy+BAh0C33+uVmPZMq7EYyhQ80NSkgN9Pcl2fkY\nBQXASy9xd8Zf/5pdvIDJ7C4uBm6+2f3n5OQYF5Gu4xaKvXu5i6MbC0P/1vVSC4UIhsexuoWKithU\nXrkSuPNO3vfPfwLvvcflFKw+2Y8/5id362S23SWlb9y62J+V3r2DXVIXXWQm12MlXB2nZcv4Zn3E\nEWwhaQYO5E53gJnHOHAgfMTI4cP8ZNiVoY9CctHRPPHeNCOxyxLMr12QP/oRP90/8gi7JrVFXlsb\nv6vTio726tEjcH7wnnu41S/An9mrF/C1rwWGe8dKKMHo3bvrmzmJYHQxdsHYt48nuG+/3ez/7DN+\nMrf2Eti9m/cRmRuuzrK2EypvIZSF8cILkZ+MImHPwH7vPRaEjg6+MezZwz5rK7fdxkv7lz/UeJua\nWBydVhIVUp8BAzgfIRkTtfPm8bKuLjiMeulSfgCzVmxOxLyW/m5OmMDzFG1tHPxx0008VwLwA2Fh\nIR9rDfeOFf1baWjgsiY33sheCLEwPI6OyQaMhWFn5Eh2K1mPHTs2sAzIjh3As886/1y7YOhqufH2\nJ7YnxJ1xhvkBjRsX+pxzz+UigvrLr5+27HMxAB8TrztKSF169kxOxrO2vAsKQkfFDRsW+HSvy8sk\ngpwc/k3W1xurWoeZa8EA3FkYu3fzZ3zyCfCHP7BAiWBkANrCOPts/mIrFRwCqL9g555r9vXqFXhM\naWnoye1//hN46qng/dZJb6X4icsNoRLinHx5+/ULtjBChVm6mb8QUpcnnkh8u16leF5k/nzezs83\n7k+N3QVlrbWWCHJzuS+NjlrSD08HDhgBczPp/bOf8e9k/Xrjxj3hhOCKEMlGBKOLaWvjm/XSpWxF\nrFvH7UfnzDHH6BvljBnAddfxeqRSGFZmzODSCHZycvhJ/uWXAxP44p0Ir68PNuvnz2e/cSSsgqFd\na6F+RCIY3iTa9wNgl6w1uCMaS5fyvxNPNPt0VJ5GC4Z+yLKWzkkEeXnAo4+abT2/YbWU3Ux6A9x0\nra7OPJgdeSS/X1fmYohgdDHWgnwATxBXVHA01EcfAbfcEljWQK87FYxw9O4NvPYa50qcfnrgeOKh\nujrYwpgSUgzXAAAbcklEQVQ/P7C3eCiKivjvHT8+cvJRY2Pwj15If379a7aWd+4Mf8zbb8fWvvix\nx3hpdeE+8ww/mGlLWguQ/u4n0iUFsPVQUWF6mmvBsEYTunFJ9evHUX/19caN/ZWvcBWE9evj62YZ\nDyIYXcyhQ6H7S/fqBZx8Mv+grHX69VO8Nf8iHuy9lXUCYDyCsX07/ws3VxGJ4cM5zHHdOvMl37o1\n+DixMLxJTg7/BsKVe2loiCwmQLCY6DBtq2AUF3MG9vHH87Z+4Cot5SS4UKHn8fDII8Bvf8tisGIF\nF1jcssWIhF0w4rEGlDLVrXv04Iiw9eu5QVp+Pif4WiMSk4kIRhdjtzCcsGOHeYqKl7PPBv7yF7Yy\nAC7RcN558QnGli0sbvEU9Rs+3KxffTX/mKzuOE1VVWa3wsxE3nmHv1M6OjBUjs62bcHuqqYmLh+u\nI5OsWJ/oTzyRG4tZv4NumTePi2TqeYXjj+e/QQdyNDUZSzlel5TPx/eMrCx+788/NxaStmR22WuE\nJwm31WqFGIlHMEpLox8Tjfx8U5wQ4JyOULkZTvjsM1PTKlasf0tWFlsZoTJ/t2/nuR3Bu9h/C2v8\nVeb0zbalJfi7oXOMOjpMpnZTE/dGCRXxZ32i1yV3koH+XRGZkiH687WF0bcv74+1Yu3Bg8EN0vRE\nujUXpCsQC6MLeOopUw8/HsFIJP368ZNbSUn8glFVFb+LjIgnKInYz9yzJ4/HnmxYWcnhxYJ30WU1\n9u1jq/X993lbZ4KHct9ot5O1GVek+a6uslKt1pAWDF2VWo+tRw9+LdZMd6uVokPrtVDopdPmaHqO\nJV5EMJKMz8dZpk8/zTfJOXNSpz9DvOVC9u93VwH27LP5Cz57dvATmWbfPvZDC97lvfd4eeqpXLZD\nZ0drt00owdDfVx1OGq0iwKJF8bUGjpWXXwZef53Xs7NZHDo6eO7QGv4ez8T32rXG2rYHv1h7i0Rr\nztTUxGHNbhDBSDK6a5yVo4/u+nGEIl4L48CBxCc+2YXLGr8ueBNdmNJecVULhVUwLrgA2LzZPFjo\nSV5doiacS6ZPn8RHRIWipMRUwdWf29LCCXfWHJB4Jr5XrzbtiK3Z6oARjJ49o7cOqK52/xAmgpFk\nQpUf1p3kuht7BdvPPnPW3c6thWEnlIVRVyeC4XX0Dc6alEpkIvqsE7kvvQQsXhz6wSIVy9cXFHD0\n3+bNHDqviWfiu7HReCWeeQZYvjz4mJEjo098NzREb10bDRGMJFNVxbkPAE86t7VxGFwqYP/yjhvn\nLDyvqiq4BawbrBbG2rU8+VlVlZo3AiFxaMGwuln69eMHkkmTuKOdFZ+PHyxmzODtw4f5e5LorO1E\ncMQRHElYWcklPDSxuqQOHAgMMT/tNG6GptEegrPOil4qaNIkUxIoXkQwkoy1lkxhYXCJj+7EmnUN\n8I8xWnmPpiYO6zvuuMSNQ1sYNTVcyO3VV3l/Mns/C93H/fcDZ54Z2sJobeWeKpMnBxfG9Pn49X79\n+AaqO+elomBYq0Db5zCcuqReeYXvGZEsA30NL7gAWLjQvSBEQwQjyej/7PXrgTvu6O7RBKIFQyl+\nEnIiZrt2cWhsIoVPWxi6f/mhQ8D3v5+49xdSix/+kCsa1NXxd09bGB98wP/3hw+zoGgL4/rrednU\nZDpPasGork5NwdChtCtWBO93amHoOZ5ISawVFbycNIkf9qw14jo6uCI1kLheOCIYSUaXAR83LvVK\nXfTvz+N7+WX2gepw30hf6ObmxHdM0xaG9sFWVop14XX69+cw2pISM2Ft9fV/7Ws8B9DRYcrNVFfz\nTbFPH/4tvflm6loYugOffR4uFpeUtlIiCcaCBVy9NtS9ZcMGbjXQ2WmCb3RDs3gRwXDAli3xn5uI\niaZkoS0M7ZbSghHpKShUMpVbdGy6LvlQWRl/2XUhPdD5EVVVnAQ6e3ZgBna/fiwm+gka4Jve00/z\nw9f27cCVV3LUUCoKRk4O8L3vBWeVx9ITQ+cmVVeHF4y5c02V3uuuM/M7gJnf2LOHXeMTJrhPXhTB\niMKaNRwjHqqLXTTWrk3tmkjarNd1dawJhQsWmPWdO004XjIsjAED2D2hzeY9e8TC8DrWxE8i4Lvf\n5XWlzIPDsGH8Xejbl6ONtmzh3+Nll5m+8R99lJqCAQAPPhj8Pe7XL3r4q0Y/tO3d6+weMmtWYD0t\nbaFUVppQeLcNyUQwoqCTfmKtFltZyYqe6hZGY6NxCVgFw1oCffduk4179tlczyeRFBbyE5AOrd27\nVywMr2Mt/tfSEvjd0w8nhYVsgbS1ce6SrgqQmwu88QYn/K1bl/hS5cnE3jSNiCf57bz/fmD9OCeC\noV3MGl3Esa4ucaHwIhhRCJVH4QQ9YffZZ12TOBQP/frxOHVBQmvy06pVbFX9/OdmAvovf+FltGqi\nsVJYyL3AtYXR0iKCkQlUVfFNrrY2dLmcwkKOyMvP5xur9ZghQ7ipEJBeNceKikxJFM2UKcFNxP78\n58BtJ+WErFGP778PXHstr3/4IUceJuI+JIIRhXg7WulJpg0bUrfEhX5q+etfeWktG93Swr3Gf/EL\nUxTu7rt5mehaWJWV3CnQKhTikvI+gwezGFRXh/5ONTayKGgL3V76QhfATGROULIpKmLL4Qc/MDWl\namuBjRt5fckSrtZrDct1ytChfN/Zs4fnejR33QX87W+Jqc0mghEF6+RvuFotofZbC4ylqmDYIysa\nGsx6797AzTebbavPOdbiadG45RZe7ttn/NFiYWQGAwZw5FOoeTHdolj73e3CMHYsWyBu/fJdia7W\nXFERWOOqZ0+26GfMYGvbScUFOwMGcLLgunXAypUc+n7OOeb10aPdjR0QwYhKczMnxQCh6y41NLAr\nx14t0npTTdUSF6Ge6pTiYon2VprWgomJroV17LFmglP7o8XCyAz0A4LOW7By6aW81OHWb7xhmiVp\nrFnU6cCkScDFFwOjRgUKhs9nojFnzAjOcnfKsGHsMj54kN9j6lTe/8gjwLRp7sYOiGBExefjmPD8\n/NDJL9pnaJ8I9vn4CX3Fiq6rVZ8onnySXQF33WX2/f73wJgx/CNevTrxn5mTw24pfQNJ1UABIbHo\nB4RQgqHRkT+DB6efQIRi2jR2t1lLn/t8gQIydSpbCbFWti4tNYKRm2vOnz49MQ9haXYr63p8Pr7w\n4UqB68gee6RDczO7olKl0GA0cnNNSQ6Af8A6wxbgL/mnn/J8R6Qfd7zoeHvtmkjVQAEhsWjBCFee\n3IuWZn4+eyaam7nH+AUX8H3Gmp/xne8Ap5wSe2XbggIjRnl5RjAS9ZuVjntR8Pn4JhZNMOxlzJOR\nr5AMDh1iX+eoUfwUYkWX/3jqqeSP48MPWVz1NU5VN56QWHSByewwdyLdg9tL5OfzvI2+R+hEWWuu\nl7Y8brwxNg+FrpqgLQxtqSfqXiSCEYXmZr7wffo4F4yODuBPf2LXTqqj80usEVJWtm5lMUk2kyfz\nhN055wDPPy+CkSlEK5fz3nvx9WxJZYqLOTJKC0ZuLkeKvfMOC2dHh7EIfvWr2N67Tx+eu+js5PfS\nD33WvBc3iGBEweqSsvdsAMzktlUwdD3/aJVfU4Wnngr/w01EZIVTPvqIr+e8eRIllSlEc5UMHtw1\n4+hKSkt5Il8LRv/+3Odi7VozBxFv3bmcnMCk20R7OUQwohBNMPQ+a4alDsW115FJVXS/8VSgb9/o\nrSYF7zB1Kpe0yCQKCjjXorqav+8lJaZSs574jlcw7A9aJ53ED2KJQia9o6AFo6gotNtGl1u25jD4\nfBxRNHNm141TENKRkhKulpxJELHlvmYNC8awYSb7+8ILeWmt3BsLWjB04h4Ru3oThVgYUdBzGMOH\nh06maWlhs9lqYezaJS4VQRDCM3EiN0g66STT83vBAq44W1QUf7dJfd/R1YATjQhGFLSF8aUvATt2\nBL/e0sK5A1bBSESCjCAI3qW0lN1P/fsbwRg2jCfEdf+PeNCCYa3MkEhcu6SIaBoRbSaibUR0U5hj\nHvC/vo6IJkQ7l4hmE9FGIjpMRBNDvWdXocNqdSlwO83NLBjaJRVPGXRBEDILnWc0aZKZ2E/ETV5H\nRSWiblQoXAkGEWUBWAhgGoCxAOYS0RjbMdMBjFJKjQZwNYCHHZy7AcD5AN51Mz63dHZyKG2fPia+\n2Y7dwtBFxARBEMKhXU4nnmhCXhORxa4TIJOV8OjWwjgFQIVSqlIp1Q7gWQD2mIeZAJ4AAKXUSgAF\nRDQk0rlKqc1Kqa0ux+YaXT21R4/IgjFoEAsGkenDm+h6S4IgeAed06Vv8EpxTTW3jByZ3ChDt4Ix\nFIC1O8Iu/z4nx5Q4OLdb0fMXQGTBsMaS19Vxa0Y3bV0FQfA2J5/M1kW64VYwnGpZGhUgNrS2GtMu\nkmBYI6Kqq71Z/0YQhMQxfrz7/trdgdsoqd0ASi3bpWBLIdIxw/zH9HRwbkTuuOOO/6yXlZWhrKws\nltOjcuiQKQEerjRIc3OgYFRVpVfLSEEQvE15eTnKy8sT8l5uBWMVgNFENALAHgBzAMy1HbMYwHwA\nzxLRZAD1SqlqItrv4FwggnViFYxkcOiQqbXk1MKoqkpehIIgCEKs2B+mFyxYEPd7uRIMpVQHEc0H\n8DqALACPKqU2EdE1/tcXKaWWENF0IqoA4ANwRaRzAYCIzgfwAIAiAK8S0b+VUucEDSDJ6EquQGTB\nsNZrqaqSpD1BELyJ68Q9pdRrAF6z7Vtk257v9Fz//pcAvOR2bG5pb48sGJ9/Drz7LnDNNZzKv3s3\nC4bMYQiC4EWkllQErC6p3NzgZiZz53Kz9kGDgNdeY9EQwRAEwauIYETA6pIaMICzuQ8fNq+3t/Py\n2GOBceO4V69SIhiCIHgTEYwIWAUjO5tLDlt7XOTlAVOmmDkM3Q5RBEMQBC8ighEBq2AAXP9l/36z\nffAg8LvfmW0RDEEQvIwIRgTa280cBsBlh3WDE4BdVLpnLsAliwFOyhEEQfAaIhgRsFsYI0cCFRVm\n2y4YZ5zBvTB0uWJBEAQvIYIRAbtgjB0LbNrE60oFCwYQf6csQRCEVEcEIwLWsFqAu+7t8hcvaW3l\nKra6dIggCILXEcGIgN3CGDyYiwsC3I9XSoAIgpBJiGBEwG5hDBkCLFsGfPIJcM89wA9+0H1jEwRB\n6GpEMCLwxReBcxKl/tq6b74J7N1roqIEQRAyARGMCKxfHxgiO3AgN2jfvp3LmuvmSoIgCJmACEYE\namvZDWWloABoamLBsFapFQRB8Dquq9V6kdparkzb1MTlQKzoIoQiGIIgZBoiGCGYNo2joIBgwcjL\n417fPp8IhiAImYW4pEJgLTAYycKQRkmCIGQSIhghIEtTWHtiXl4e98DIygoMuRUEQfA6Ihg2Nm/m\nKCgN2TqK5+YCNTXijhIEIfOQOQwbe/bwcvLk0D28c3O5Yq3OyRAEQcgUxMKwsXcvL+fNA9auDX59\nwABe6t4XgiAImYIIBngSu6GB1/fuBUaMAObMCX2sdkVZa0wJgiBkAiIY4DarEybwemMjcMUVwdFR\nGj2nYe3tLQiCkAnIHAZMzgXAyXpO5ieUSt54BEEQUhERDD95ebwMld1t589/lq56giBkHmktGPfd\nxz0qvvWt+N9j/35eDhzIy4MHjXiE47vfjf/zBEEQ0pW0Fozrr+fy424EY+dOXh59NC+dWBiCIAiZ\nSNpPervNtm5sBLKzWSgAXkazMARBEDKRtBeMbJc2UmMjUFICrFwJvPACZ3EXFSVmbIIgCF5CBKPR\nuKMuugjYscNsC4IgCIa0F4xEuKSOOgo44QTezsuTpDxBEIRQpL1g2KvJxsq2bVzmY9063q6tdT8m\nQRAEL5L2gnHokLvzFy8GZs/m9aws9+MRBEHwKmktGP37A62t8Z+vFLBrFzBmDG+LYAiCIIQnrQVj\n0CB3grFvHxcT1AUFRTAEQRDCk9aCUVjoTjB27Qos8dEjra+GIAhCcknrW2RRUfyCsWcPMHFiYKFB\nEQxBEITwuL5FEtE0ItpMRNuI6KYwxzzgf30dEU2Idi4RDSSiZUS0lYjeIKKQ7YrcWBi7dvEyJ8fs\nGzlSyoIIgiCEw5VgEFEWgIUApgEYC2AuEY2xHTMdwCil1GgAVwN42MG5NwNYppQ6GsBy/3YQAwcC\n7e1AZ2ds41bKhNFef73ZX14OVFTE9l6CIAiZglsL4xQAFUqpSqVUO4BnAcyyHTMTwBMAoJRaCaCA\niIZEOfc/5/iX54X68Lw8zsNoa4tt0KtXA1dfDcyaBZx+utk/YABPpAuCIAjBuBWMoQB2WrZ3+fc5\nOaYkwrmDlVLV/vVqAINDfXhzM7uUamtja2jU0cFLXdpcEARBiI7b8uZOb9Pk8Jig91NKKSIK+Tmd\nnXegvR0YPhy47royPPRQmaPBHDzIy02bHB0uCIKQtpSXl6O8vDwh7+VWMHYDsDY0LQVbCpGOGeY/\npmeI/bv969VENEQpVUVERwCoCfXh9913B5YsAbZuDZy8joYWjLvucn6OIAhCOlJWVoaysrL/bC9Y\nsCDu93LrkloFYDQRjSCiXgDmAFhsO2YxgMsAgIgmA6j3u5sinbsYwOX+9csBvBxuAP378zKWEiEH\nDwLf/jZw1VXOzxEEQch0XFkYSqkOIpoP4HUAWQAeVUptIqJr/K8vUkotIaLpRFQBwAfgikjn+t/6\nbgB/J6J5ACoBXBxuDLr4YEOD83FLkyRBEITYIRXLbHEKQURKKYUnnwTmz+cmSEOHAsuXRz/3t78F\nqquB3/0u+eMUBEFIJYgISikn88pBpHVPbwC47DLO1v7614EtW5ydc/CgWBiCIAix4oliGIWFZr2l\nJfrxTU2S0S0IghArnhAMawHBL76IfOzhw+y2EgtDEAQhNjwhGAMGmPXKysjHLl8OrF8vgiEIghAr\nnhAMIg6RPfdc4MknIx+re4CLS0oQBCE2PCEYAPCnPwFXXAEsWxb+GKVMz26KK0ZAEAQhc/GMYADA\nWWdxfalwkcIffADMmcPrTibHBUEQBIOnBCMvj9usNjaGfl0XHQQ4DFcQBEFwjqcEAwByc4E77wz9\nmu6bMXUqd+sTBEEQnOM5waiqAu69l8Nn7fh8vJSy5oIgCLHjOcHQfPJJ8L7mZl6KYAiCIMSO5wRj\ngr9j+IoVwa+JhSEIghA/aV9Lys6HHwK/+U3oBD4tGNY+3oIgCIIzPGdh9O7NlWt1voUVnw+44Qbg\nF7/o+nEJgiCkO54TDAAoLjaCoSOjAKlSKwiC4AZPC8a2bcCYMUY0GhoC604JgiAIzvGkYAwdCuza\nxU2Stm4F3n2X99fVAQUF3Ts2QRCEdMWTgjFsGLB7N/DMM7z9+OPAgQPA3/4mgiEIghAvnhSMrCxe\nPvQQcMwxwEcfAevW8T6pUisIghAfnhQMgCe4AWDwYKCtDejh/0vz87tvTIIgCOmMZwUjN5eX2dlA\nayu3ZT37bJPYJwiCIMSGZwVDU1/PFsbBg2JdCIIguMFzmd5WPvoIaG8HzjxTcjAEQRDc4mnBOPlk\nFgxtYYhgCIIgxI/nXVLZ2dyBr75eBEMQBMENnhcMIq4vtWMHZ4ALgiAI8eF5wQC4f/fjj3MGuCAI\nghAfGSEYmsLC7h6BIAhC+pIxgjFxInDqqd09CkEQhPSFlFLdPYa4ICLldOxEwOrVLBqCIAiZDBFB\nKUXxnOvpsFpNmmqiIAhCSpExLilBEATBHSIYgiAIgiNEMARBEARHiGAIgiAIjohbMIhoIBEtI6Kt\nRPQGEYXsZUdE04hoMxFtI6Kbop3v3/82ETUR0R/iHZ8gCIKQWNxYGDcDWKaUOhrAcv92AESUBWAh\ngGkAxgKYS0RjopzfCuBnAP7HxdgyivLy8u4eQsog18Ig18Ig1yIxuBGMmQCe8K8/AeC8EMecAqBC\nKVWplGoH8CyAWZHOV0o1K6U+ANDmYmwZhfwYDHItDHItDHItEoMbwRislKr2r1cDGBzimKEAdlq2\nd/n3OTlfsicEQRBSiIiJe0S0DMCQEC/91LqhlFJEFOoGb99HIfZFOl8QBEFIFZRScf0DsBnAEP/6\nEQA2hzhmMoCllu1bANzk5HwAlwP4Q4TPV/JP/sk/+Sf/Yv8X733fTWmQxeCb+m/8y5dDHLMKwGgi\nGgFgD4A5AOY6PD9irZN4a6EIgiAI8RF38UEiGgjg7wC+BKASwMVKqXoiKgHwZ6XUDP9x5wC4H0AW\ngEeVUndFOt//WiWAfgB6AagDcJZSanN8f6IgCIKQCNK2Wq0gCILQtaRlpne4ZMBMgIhK/YmNG4no\nEyL6gX+/o0RKr0FEWUT0byJ6xb+dkdcBAIiogIieJ6JNRPQpEZ2aideDiG7x/z42ENHTRNQ7U64D\nET1GRNVEtMGyL+zf7r9W2/z307OivX/aCUaUZMBMoB3Aj5VSx4GDCr7v//ujJlJ6lB8C+BQ8mQdk\n7nUAgP8FsEQpNQbACeDAkoy6Hv750qsATFRKjQO7wi9B5lyHx8H3Rish/3YiGgueVx7rP+chIoqo\nCWknGIicDOh5lFJVSqm1/vWDADaBc1ucJFJ6CiIaBmA6gEdggiQy7joAABHlA/iKUuoxAFBKdSil\nGpB516MR/FDVl4iyAfQFB9xkxHVQSr0Hnve1Eu5vnwXgGaVUu1KqEkAF+P4alnQUjEjJgBmF/2lq\nAoCVcJZI6TXuA3ADgE7Lvky8DgAwEkAtET1ORGuI6M9ElIsMux5KqQMAfg9gB1go6pVSy5Bh18FG\nuL+9BHz/1ES9l6ajYMgsPQAiygPwAoAfKqWarK/5e9d6+joR0bkAapRS/0aYEOxMuA4WsgFMBPCQ\nUmoiAB9sbpdMuB5EdBSAHwEYAb4h5hHRpdZjMuE6hMPB3x7xuqSjYOwGUGrZLkWgSnoeIuoJFou/\nKqV0/ko1EQ3xv34EgJruGl8XcTqAmUT0OYBnAHydiP6KzLsOml0AdimlPvZvPw8WkKoMux6TAKxQ\nSu1XSnUAeBHAaci862Al3G/Cfi8d5t8XlnQUjP8kAxJRL/CkzeJuHlOXQUQE4FEAnyql7re8pBMh\ngfCJlJ5BKXWrUqpUKTUSPKn5llLqO8iw66BRSlUB2ElER/t3fQPARgCvILOux2YAk4moj/+38g1w\nUESmXQcr4X4TiwFcQkS9iGgkgNEAPor0RmmZhxEuGTATIKIzALwLYD2M+XgL+D86ZCKk1yGiKQB+\nopSaGSkh1OsQ0YngAIBeALYDuAL8G8mo60FEN4JvjJ0A1gD4LjgR2PPXgYieATAFQBF4vuLnAP4P\n4ZOkbwVwJYAOsHv79Yjvn46CIQiCIHQ96eiSEgRBELoBEQxBEATBESIYgiAIgiNEMARBEARHiGAI\ngiAIjhDBEARBEBwhgiEIgiA4QgRDEARBcMT/A1j2txjRZaLMAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c453050>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Using MKL"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file sdemkl.pyx\n",
      "#cython: boundscheck=False\n",
      "#cython: wraparound=False\n",
      "##cython: profile=True\n",
      "cimport cython\n",
      "from cython.view cimport array as cvarray\n",
      "cimport numpy as np\n",
      "import numpy as np\n",
      "from libc.math cimport sqrt, fmin\n",
      "\n",
      "ctypedef void * mklRandomStream\n",
      "cdef extern from \"mkl.h\":\n",
      "    int VSL_BRNG_MT2203\n",
      "    int VSL_RNG_METHOD_GAUSSIAN_BOXMULLER\n",
      "    int mkstream \"vslNewStream\"(mklRandomStream *, const int  , const int)\n",
      "    int mklgaussian \"vdRngGaussian\"(const int  ,mklRandomStream  , const int  , double [], const double  , const double  )\n",
      "\n",
      "cdef  mklRandomStream mklstream\n",
      "\n",
      "\n",
      "cdef double A(double x):\n",
      "    return -x\n",
      "cdef double B(double x):\n",
      "    return 1.0\n",
      "\n",
      "cdef  double a = 0\n",
      "\n",
      "\n",
      "\n",
      "cdef class RandN:\n",
      "    cdef mklRandomStream _strm \n",
      "\n",
      "    def __cinit__(self):\n",
      "        mkstream(&self._strm, VSL_BRNG_MT2203, 111)\n",
      "\n",
      "    def __call__(self, double sig):\n",
      "        cdef double r\n",
      "        mklgaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, self._strm, 1, &r, a, sig )\n",
      "        return r\n",
      "\n",
      "    cdef double n (self, double sig):\n",
      "        cdef double r\n",
      "        mklgaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, self._strm, 1, &r, a, sig )\n",
      "        return r\n",
      "\n",
      "\n",
      "\n",
      "cdef RandN rv = RandN()\n",
      "\n",
      "cdef double sdeevolve(double x, double tout, double dt):\n",
      "    cdef double t, ax0, bx0\n",
      "#     cdef np.ndarray[np.float64_t] rands = np.random.randn(int(tout/dt) + 1)\n",
      "    cdef int i = 0\n",
      "    cdef bint dobreak\n",
      "\n",
      "    \n",
      "    t = 0.0\n",
      "    \n",
      "    if tout < 1e-9 : t  = tout +1\n",
      "    \n",
      "    while t < tout:\n",
      "            \n",
      "        if t + dt > tout:\n",
      "            dt = tout - t\n",
      "            dobreak = True\n",
      "        \n",
      "\n",
      "        ax0 = A(x)\n",
      "        bx0 = B(x)\n",
      "        \n",
      "#         x += dt * ax0 + sqrt(dt) * bx0 * rands[i]\n",
      "        x += dt * ax0 +  bx0 * rv.n(sqrt(dt)) \n",
      "        \n",
      "        if dobreak:\n",
      "            break\n",
      "\n",
      "        t += dt\n",
      "        i+=1\n",
      "        \n",
      "    return x\n",
      "\n",
      "\n",
      "def sde(double x, np.ndarray[np.float64_t] tout, double dt=.01):\n",
      "    \n",
      "    cdef double t, tNext, t0\n",
      "    \n",
      "    t = tout[0]\n",
      "    \n",
      "    cdef np.ndarray[np.float64_t] xout = np.empty_like(tout)\n",
      "    dt = min(dt, np.diff(tout).min())\n",
      "    \n",
      "    cdef int i\n",
      "\n",
      "    xout[0] = x\n",
      "    for i in xrange(tout.shape[0]-1):\n",
      "        t0 = tout[i]\n",
      "        tNext = tout[i+1]\n",
      "        x= sdeevolve(x, tNext-t0, dt)\n",
      "        xout[i+1] = x\n",
      "        \n",
      "    return tout, xout\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Overwriting sdemkl.pyx\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%file setup.py\n",
      "\n",
      "from distutils.core import setup\n",
      "from distutils.extension import Extension\n",
      "from Cython.Distutils import build_ext\n",
      "import numpy as np                           # <---- New line\n",
      "\n",
      "ext_modules = [Extension(\"sde\",\n",
      "    sources = [\"sdemkl.pyx\"],\n",
      "    libraries = ['mkl_intel_lp64', 'mkl_sequential', 'mkl_core', 'pthread', 'm'],\n",
      "    )]\n",
      "\n",
      "setup(\n",
      "    name = 'sde',\n",
      "    cmdclass = {'build_ext': build_ext},\n",
      "    include_dirs = [np.get_include(), '/opt/intel/mkl/include'],         # <---- New line\n",
      "    library_dirs=['/opt/intel/mkl/lib'],\n",
      "    ext_modules = ext_modules\n",
      ")\n",
      "# -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Overwriting setup.py\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!python setup.py build_ext --inplace"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "/Applications/Canopy.app/appdata/canopy-1.4.1.1975.macosx-x86_64/Canopy.app/Contents/lib/python2.7/distutils/dist.py:267: UserWarning: Unknown distribution option: 'library_dirs'\r\n",
        "  warnings.warn(msg)\r\n",
        "running build_ext\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "cythoning sdemkl.pyx to sdemkl.c\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "building 'sde' extension\r\n",
        "gcc -fno-strict-aliasing -fno-common -dynamic -arch x86_64 -DNDEBUG -g -O3 -arch x86_64 -I/Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include -I/opt/intel/mkl/include -I/Applications/Canopy.app/appdata/canopy-1.4.1.1975.macosx-x86_64/Canopy.app/Contents/include/python2.7 -c sdemkl.c -o build/temp.macosx-10.6-x86_64-2.7/sdemkl.o\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "In file included from sdemkl.c:314:\r\n",
        "In file included from /Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4:\r\n",
        "In file included from /Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17:\r\n",
        "In file included from /Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1761:\r\n",
        "\u001b[1m/Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1m\r\n",
        "      \"Using deprecated NumPy API, disable it by \"          \"#defining\r\n",
        "      NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\" [-W#warnings]\u001b[0m\r\n",
        "#warning \"Using deprecated NumPy API, disable it by \" \\\r\n",
        "\u001b[0;1;32m ^\r\n",
        "\u001b[0m"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 warning generated.\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "gcc -bundle -undefined dynamic_lookup -g -arch x86_64 -headerpad_max_install_names -arch x86_64 build/temp.macosx-10.6-x86_64-2.7/sdemkl.o -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -o /Users/noah/workspace/ergnum/results/sde.so\r\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sde as sdemkl\n",
      "tout, xout = sdemkl.sde(x0, tout, dt=.01)\n",
      "plot(tout, xout)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10c4c61d0>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXeYHNWV9t+jyaMsJI0CEgIlJIIRYIHBhiGasAaDYQEH\nMA7wERYv2MBnbLDEOu4aDCyYZdcYs/gDr5cMBgMGxoAMIgqEgomy8kijkSYnae73x5njul1dHaqr\nq6u76/yeZ56qrq7pqq6uuu894Z5LxhgoiqIo8WNY1CegKIqiRIMKgKIoSkxRAVAURYkpKgCKoigx\nRQVAURQlpqgAKIqixJRAAkBE04joeSJaQUTvEtFlHvs0ElEbEb019Pf9IMdUFEVR8kNlwP8fAHC5\nMWYZEY0A8AYRPWOMWeXa78/GmFMCHktRFEXJI4EsAGPMZmPMsqH1TgCrAEzx2JWCHEdRFEXJP3mL\nARDRDAALACx1vWUAHEZEbxPRE0Q0P1/HVBRFUXInqAsIADDk/rkfwLeGLAGbNwFMM8Z0E9GJAB4G\nMCcfx1UURVFyh4LWAiKiKgCPA3jSGHNTFvt/DOAgY0yra7sWJVIURckBY0xObvagWUAE4E4AK1M1\n/kTUMLQfiGghWHRavfY1xuifMfjBD34Q+TkUy59eC70Wei3S/wUhqAvocABfBvAOEb01tO0aANMB\nwBhzB4AzAFxERDsBdAM4O+AxFUVRlDwQSACMMS8hgxVhjLkNwG1BjqMoiqLkHx0JXIQ0NjZGfQpF\ng14LB70WDnot8kPgIHC+ICJTLOeiKIpSKhARTBRBYEVRFKV0UQFQFEWJKSoAiqIoMUUFQFEUJaao\nACiKosQUFQBFUZSYogIQA4zhP0VRFBsVgBhw5ZXA7rtHfRaKohQbKgAx4N13gY0boz4LRVGKDRWA\nGDBhQtRnoChKMaICEAPGj4/6DBRFKUZUAGJATU3UZ6AoSjGiAhADhg39ypoJpCiKjQpADBgc5OXO\nndGeh6IoxYUKQAwYGOBlb2+056EoSnGhAhAD+vt5qQKgKIqNCkAMUAtAURQvVABigFoAiqJ4oQIQ\nA8QC6OoC7ror2nNRFKV4UAEoc9rbgb4+Xn/sMeBrX4v2fBRFKR4qoz4BJVxGj3bW33knuvNQFKX4\nCGQBENE0InqeiFYQ0btEdFmK/W4hoveJ6G0iWhDkmEru3H8/L3ftivY8FEUpDoK6gAYAXG6M2QfA\noQAuIaJ59g5EdBKAWcaY2QAuAHB7wGMqARGXkKIo8SaQABhjNhtjlg2tdwJYBWCKa7dTANw9tM9S\nAGOIqCHIcZXsuOEGXt5+OzBrlrNdsoIURYk3eQsCE9EMAAsALHW9NRXAOuv1egA6PUkB+M53eHne\necCYMc52tQAURQHyFAQmohEA7gfwrSFLIGkX12vPsmSLFi36+3pjYyMaGxvzcXqxp7YWuO464JRT\n+LUKgKKULk1NTWhqasrLZ5EJWCKSiKoAPA7gSWPMTR7v/weAJmPM74ZerwZwpDGm2bWfCXouikN7\nO2cADRvmBH1vuw247DJg9Wpg9uxoz09RlPxARDDGuDvZWRE0C4gA3AlgpVfjP8SjAM4d2v9QADvc\njb+Sf1avBg48MDHj55JLgLlzNQagKAoT1AV0OIAvA3iHiN4a2nYNgOkAYIy5wxjzBBGdREQfAOgC\ncH7AYypZsGIFsM8+3u+tX5/6PUVR4kMgATDGvIQsrAhjzKVBjqP4p7XVeyrIVauAE07QyWEURdFS\nEGVLVxcwfHjq96U+kKIo8UUFoEzJJADbthXuXBRFKU5UAMqUVAKwYQMHgrduLfw5KYpSXKgAlCld\nXcCIEcnbp0wBJk8GNm8u/DkpilJcqACUKelcQIcdBjz7bGHPR1GU4kMFoExJJwAHHgi8915hz0dR\nlOJDBaBM2bEDGDXK+72GBmDLlsKej6IoxYcKQJny0UfAnnt6vzdxItCsY7GVEOnp4dnnWlujPhMl\nHSoAZchLL3GQd4q7MPcQEyeqBaDkj/feSy4wuGQJzz+9227A4GA056VkRgWgzFi+HDjjDC4BPSzF\nrzt6NNDby/MEnH56Yc9PKT/mzgV+/nPn9W23Accd57zevr3w56Rkh84JXEZ8/DFwzTXs3jn55NT7\nEbEVcMMNwIcfFu78lPLFbuTdlYqbm9kSUIoPtQBKnK1buUEHgIsuAh5/nCeC+cIX0v+fuoGUfGK7\neaZN4+W11/JSx5wULyoAJY6UdPif/wGeeorX9947tftHmDwZ6OhI3KamupIrdtlxgF1C118PnHYa\nZ6QpxYkKQInTOTT/2tlnO9uyKfT26U8nvt66FRg3Ln/npcQLWwDa250U5OHDnXtUKT40BlDiuHvx\nd90FnHNO5v+bOTPx9fLlvDTGcSkpSrbYArB9OycaAFyOpKsrmnNSMqMWQInjFoCTTwZqajL/39ix\nia/ffpuXWiZayQURgJ4e4MEHgXnz+LVaAMWNCkCJ4xYAr0lgvHALwLJlvFQBUHJBBKCtjUea77cf\nv1YLoLhRAShxbAE466zs3TduAVizhpcqAEouiAB0dQH19c52tQCKGxWAEqejAxg5ktevvjr7/3MH\nfKU0hE4Yr+SCCEB3d2IRQrUAihsVgBKno4NTOgHv+v+pkCAdwDncmzcDtbVqASj+kPz/nTt56bYA\ndttN604VMyoAJU6uAmC7iioquOc2YYIKgOIPuV96e9kKcJchX7gQePnlaM5NyYymgZY4HR1O0Tc/\nAuCmpob/1AWk+EGKwK1aBVRWch0q2wKYNo0tAE0vLk7UAihxOjqASZN43X7wsuHDD526LbvvDlRV\nqQWg+EP8+3/9Ky+XLUu0ACoqWBi0Y1GcBBYAIvo1ETUT0fIU7zcSURsRvTX09/2gx1QcOjqcUZcV\nFf7+d6+9gCOP5DEATz+tAqD4x10+ZPPm5Jno6up4fIBSfOTDArgLwAkZ9vmzMWbB0N8P83BMZYiO\nDu69B2H//dlUr67Wnprij9ZWJwtt/nxO+ZwxI3GfQgjAzp3A+vXhHqMcCSwAxpgXAWQqI6bev5Do\n6OA5fvMx85JaAIpftm8Hpk7l9X324eXcuYn7DA4CH3wQ7nk89JATb/DD2rXA+eeHc06lQCFiAAbA\nYUT0NhE9QUTzC3DM2CDjANwDu3JBBUDxy7ZtThLCfvsBn/wkcILLH7B1K3DEEeGeR20tL++919//\nPfkk8Jvf5P10SoZCZAG9CWCaMaabiE4E8DCAOV47Llq06O/rjY2NaGxsLMDplTb2QLCgqAtI8UNz\nM/eer7kGeO457oS8+mo05yLZSG1t/v5PBrCVUpZSU1MTmtyz7uRI6AJgjOmw1p8kol8S0ThjTJLT\nwhYAJTvyKQBqASh+eOMNdrtcdx3w4x/7T0LIJyIAfuceEMHo7MzfcxQ27s7x4sWLc/6s0F1ARNRA\nxNpKRAsBkFfjr/inp4cH4PhN/0yFWgCKH95/Hzj1VKf6bG9vdOfS28s9eL8WgEyoFNfZ8QJbAER0\nH4AjAYwnonUAfgCgCgCMMXcAOAPARUS0E0A3gLNTfZbiD6n7ny/TVdP1FD+0tPDocSHVvXPsscCf\n/hTuufT28ngYvxaAjGOIa8G6wAJgjEk7/Ygx5jYAtwU9jpJMvm9aFQDFD9u2OZk/QOp754kn+N7y\nYsMGHkR29NHBzqW3l8tQ79jBAeENG7KbiL67O3EZN7QURAljjDMPcD6or4/vg6D4p6XFaWSvvz5x\nWlKbykoOtg4OJs9VfcQRwEcf8b0cBLEAmps5HtDcnJ0AiAUQ14qlWgqihNm40cnBzgdqASh+2LbN\nmYDo2muB2bO99yPi+JJXgsGmTfk5l+9/n0cgb97Mr92T1Keiq4uD1yoASsmxdSswcWL+Pk8FQPHD\ntm3Z9bKB1AkG+cgckgyg7dsdQcm2Qe/q4mcorpavCkCJsnMnZzxk+wBmg7qAFD+0tGQ/BWmYArBh\nAy9lTgIA+NSnsvtfKYOuFoBSUrS08Kxebp9qENQCULLFGH8WQE2N01O3ycf9u34937v//d/+/1cs\nABUApaRwT7yRD1QAlGzp7mbffrZjULwsgEMPTa4mmuu5fPrTwB57+P/fri7OHrLn1o4TKgAlSl+f\nMwAnX6gLSMmWFSuAWbOy399LAJYuddaDZAENDPDnu8mmV9/dzd9j48bcj1/KqACUKGEIgFgAzzyT\n389VyoudO4FvfIN78NmSaZS5l3soWwYGuIyJmyefTN62fTtw553O664uYO+9gXXrcj9+KaMCUKKE\nJQCtrcDxxycG1BTFZsMGYPlyfynIXgIgFTzr64O5YGwB+PKXgRtu4PUzz+TxACtWAH/5C/DOO1wt\n9Bvf4Pd37eLnaJ99gJUrg49FKEVUAEqUsFxAUhOlvT2/n62UD5VDw0f9ZPB4CYCMDh4/PtiodlsA\n7rkHuOIK571Jk4B99wUOP5xHG9tB554evuf33Zefpw8/zP0cShUVgBIlLAtg1SpeVwFQUiHWoR8r\n0UsAurrYBz9yZP4EIB0dHYl1s7q6WACIOA6wdm3u51CqqACUKGEJgOC3qqISH6Thv/ji7P/HLQC9\nvexyqa1lAciXCygVFRV8fLcASCbd1KnxDASrAJQoYbmABBUAJRW7dnHZBz+j0GtqEgWgrQ0YM4Yb\n5BEj8m8BLFiQ+HrECF6KC4gIWLbMEYApU1QAlBIibAsgH3MMK+XJzp1OHCBbamoSx5i0tQGjR/N6\nrhbATTcBv/gFcOmlyQLw5pvAW285ryXAa7utNm92Oj2jRsVzLIAKQIkStgD88Y/5/WylfNi5038J\nh+HDE8eY2AKQqwVw+eVOwNdrHMDMmc66xLRsEbrkEud1XMfAqACUKGG6gC6+OJ7msJIduVgAw4cn\nDszasSNRANavD3ZOXjEAryke3SPdJf+/vj6eo+BVAEqUMASgthZYswb4whfiO0OSkpldu4ILQHu7\nIwB1dVzOOcjYk1RB4LY2HhsguBt5iQ3U1akFoJQQYQgAwPVUggbllPLl4Ye5DLlfF1B9faIA9PQ4\nLkcJDgcZDZxKAEaN4onr7ePaDA4656cCoJQMYQkAwAIQx4CYkpnTTuPga1ALwK7fIw1/kEnl01UV\nHTvWWe/pAS680HktE8eoACglRZgCEHRgjlLe5MMF1N/v9NrzIQDpBi66BWDhQkd8bAHQGIBSMoRt\nAagAKKnIRxDYtgB2352XQRrgDz5I/Z4tAN3d7Hr6xCf4tbiA4hoD0EnhSxR1ASlRsWuXU8gtW9xz\nTdgWwOLFwN13524BnHgicOSRqd+3BaC9nc9lyRKOZ9iZSHG851UASpQwBUAezFx6ekr5k4sLqLo6\nMchrWwCVlVy0LVcBeOKJ9O9/6lPA9dcDP/oRl4Ouq+N7/MwznX0mTcrfBPWlRGAXEBH9moiaiWh5\nmn1uIaL3iehtIlqQaj8le3p7/ffC/FBbG8wnq5Qfkq3T05PbSGBbAGwLAAj3fqurA669lv38S5cm\nDngUxo9nF1Xc7vl8xADuAnBCqjeJ6CQAs4wxswFcAOD2PBwz9oRpAQD8kMTtYVDSI3Ghtjb/aaDu\nWkD9/YmjdwvR4ZDpJ70EgIgLwq1eHe45FBuBBcAY8yKAdDN7ngLg7qF9lwIYQ0QNQY8bd8IWgNra\neGZFKKnZsYOXbW3BLIAVK4Df/Ca4BTBhQm4j1r0EAADOOAO4/37/n1fKFCILaCoAe8K19QB2L8Bx\ny5pCCIBaAIpNSwsvgwrAt77FpR9sC8AdJM6Gvr7UjXk6Uk1Nuc8+PBI+ThQqxEeu156Try1atOjv\n642NjWhsbAzvjEocFQCl0IgADA7m5gISAZCcfdsCGDvWfwXa7m6nnLMfJAXUzR57lIYANDU1oamp\nKS+fVQgB2ADAGoyN3Ye2JWELgJIeFQCl0LS08BwAW7YEswBEAGwLQD43WwYGeJnNTGA2kyenFq8J\nE4Bt2/x9XhS4O8eLFy/O+bMK4QJ6FMC5AEBEhwLYYYxpLsBxy5qennCzgHIxyZXyZtMmp8RyvgWg\noYEncM+W7u7ECYyyJd3E73G85/ORBnofgL8AmEtE64joa0R0IRFdCADGmCcAfEREHwC4A4CPieSU\nVHR1eZe7zRdqAShu3n8fOOAAXs/VBXTrrU6+vd1792sByHy++SSOAhDYBWSMOSeLfS4Nehwlkc7O\n3Pyf2aICoLj58EPg7LN53WsClnSIACxb5mwLagH4vf9nzADmzUv9fhwFQGsBlShdXU4t8zBQAVDc\nNDcDs2bxut/GVwRg8mRnm92D92sB5OICWrmSyz+kQuoB3Xdf6kyhckMFoATZuJFT8cK2AOLWG1LS\n09IC7Lknr/tNv5QOhV2XZ3crGXziRH8WQC4uoLq69JZLVRWXufjiF4FVq/x9dqmiAlCCnHgiL/1m\nQPhBRwIrNoODnKYpjbbfDLSqKr6nWlt5Ll4AmD7deX/cOC7Glm3Pe9my/FvAZCWrx6UargpACRJk\n6rxsUReQYtPaykkHfn3/NuPGcQB4zBjgmGP4tTBsGNfj2bo1u89autTpCIWBCoBStOy7b/jHUAFQ\nbFavBubMcV5LHX0/bNgA/PrXfG/96U+JPW6AA8HZxgF6e4EpU/yfQ7aoAChFy8SJwM03h3sMFQDF\n5t13gf32c17nIgAy+1Yq95GfOEBfX7jjYFQAlKLFnlA7LCQIbAzw1a+mH0CjlD9tbYkB3FwE4NRT\neZnKjbTbbtmPxA2rHPrgIHDRRSoAShFTKAH42c94yP3ddyfWclfiR29v4j2XyyDEq67iZSoLoK4u\n+/ssLAEgiteUqCoAJUghBKCtjZcybL+7m9PjChGAVooP+57761+dTB4/iGikEoCamuzdjmFOiDRy\npAqAUsQUQgA2b+almORtbTxAxh7JGWfuuAOIU7Fa+56bMye3FORMAuAn7hT2nNgqAErRUggBkICf\nCICUArYHyGSbsleO/P73wJ//HPVZFA63CygXRABS9dz9CECYFoAKgFIQ7r8/N996a2u4ZSAANvHn\nznUE4PnneWmPDp440RGGcuWFF4Bf/jJxmzHAc89Fcz5RkY/qs6NG8TKV5VRbC3zve8C//mvmz1IB\nyA8qABFy5pnAK6/4+5/mZmDdOqcqY5gMH+4IwNVX81IEQLKC/va38M8jSq68MtnfHZfGwSYfVmdV\nFd839gAwG3HpPPNM5s8KMw1UBUAJHZlfVXpF2dLayj3vMMtACLYACCIAMiHHxx+Hfx5R4hX0jmON\npHy4gDIhDXo25SB6ezUGkA9UACLio494KQ1ptvT05L8OeipGjEh28UjjJw+pZAuVK6kEYMoUThmM\ny/iIQqUeA5mfCWM4Oy2s+TBUAJTQEQHwW3a2uzv8B1GwLQARne5uXkpZ4HIfLezVGPX28rUZNswZ\n3VruFOK+EwHIFBdrb+f7MUhdonSoACih8+GHvPQbBC6kBSACcPLJwEknOccHnCH7cRGAhx92ZrKS\n3nB1dXzqxre3A6NHh3sMcelkuqYtLTxqOCxUAJTQaW3lpd8GpBCmuCAuIMnOmDcv2f9d7gIgLqDT\nTgNuugl47LF4CkBbm/94lV+k1EQ2AjB+fHjnoQKghI40nLm4gAptAdTWctbR974XPwGwXUADA8Ap\np8RTAAphAUiV20zXdPt2LikdFsOHswDEIb6jAhARuQpAIS0AWwAA7zlTy1kANmzgP2HdOl7+4Q/c\nOMRFAAYHuUEMe+zJ1KnAeedlvqfCzkiqrub4Thx+WxWAiBDffzELgO0CApw5U//t35x9ylkA3Pno\nb77Jyxtv5FHAIgCLFvkfz1FKdHZyZ6CiItzjEAG3385lSB56KPV+YQ4CE/r7geuvD/cYxYAKQERI\nHnOxZwHt3JlsAUhVR6C8BaCzk0et1tYClZVO5pYgArB4MfAv/xLJKRaEQvj/Bbm3040GLoQAAMBP\nfxr+MaJGBSAienv5ofKbBdTVFV7+sxtJs5NlfT0f36acBUB+o95eYObMxPeeeoqvi/x+fiY0LzXa\n2sL3/9ssWZL+/UIJwHHHhX+MqAksAER0AhGtJqL3iehqj/cbiaiNiN4a+vt+0GOWA319/FBlsgDu\nuYcnyxbEHC8Ep5/OS8mEqatLHvhVzgLQ1+eIrdvqmjkzMQZQzjWR2tsLZwEA/FykG2B40UXhn8ON\nN3LWW1BOOKG4K+gGEgAiqgBwK4ATAMwHcA4ReV22PxtjFgz9/TDIMcsF6V2mE4DmZuDccxNvoEIE\n4wTp9ck5xk0AentTV7AcMSI+QeBCWwCZBADIfuawXKmqCv7bbtzIluKLL+bnnMIgqAWwEMAHxpg1\nxpgBAL8DcKrHfuSxLdZI45LuJlu7lpf29HuFFADBLQBjxzoVMuMiAO66MyNHJsZwhpWxMzUKC2Dj\nRi6WmCoVM+wUzaoq/2Va3Hzxi7wMc8xCUILetlMBrLNerx/aZmMAHEZEbxPRE0Q0P+Axy4JsLADp\n5diNbBQCIH7u+np+KKqq2Ax/4YXyFwDp+YoFsPfevKyrA5qanIA4lXEXp9AWgNzf99/vFE0UpPRG\n0MY5E/mw7kSkCv28+qEy4P9no8NvAphmjOkmohMBPAxgjteOixYt+vt6Y2MjGst4yqW+vuwFwM69\n7+qK1gIAnKCwnwk8So2lS3kuZEl57evj4OS8ecDLLzsN/uuv87KcLYDHHwc++cnCHY8IOPponnNh\n27bEyeil4Q/7vsuHBTB7NneS8j2NalNTE5qamvLyWUEFYAOAadbraWAr4O8YYzqs9SeJ6JdENM4Y\n0+r+MFsAypm//Y1n05o1KzHAa7NrF7BlC6/bAlDIIDDAM4MdfTSvu90g5SwA557Lv430/LduBQ47\njNelLtL3vsdWwJIl5SsAr77K5S9+9KPCHvc3vwGmT2cBkMKDgNMZCfu+q64OJgBr1gB33snr+bZW\n3J3jxYsX5/xZQQXgdQCziWgGgI0AzgJwjr0DETUA2GKMMUS0EAB5Nf5xYs0aYP58LmiVKph1+unA\no4/yulsACmkBvPOOsy693vVDEl/OAiBZPbW1wB578BwMbo44ghtIoHxdQC+9BFx6qTNFaKGYNo0z\naNzPhwjAlVeGe/ygQWAZNAiE764KQqB+izFmJ4BLATwFYCWA/zHGrCKiC4nowqHdzgCwnIiWAbgJ\nwNlBjlkOtLbyrEjp/Ixvv+2sRykAbmbPdtbLQQA+/phdOm6kWF9tLbB6NVcDdVNb6xQNK9e6MW1t\nqWfwCpuGBg4G2/T3A5MnA5/7XLjHrq4GHnkEePDB3P5fZspbuLC4BSCoBQBjzJMAnnRtu8Navw3A\nbUGPU060trJfM50ATJzIN9GkScUlAHY2SDkIwA03ALfdltiA2w/sqFH8Pb0GHtXWOpZCLnM7lwId\nHTz5TRQccEBiTxrg5yWseQBsZMa9665zxsP4YWAA+M53uHBdMQtAUXkud+wo356UTTYWwNShXKrd\nd09sZKMIAtvYDWEqAXjoodL5HSW7Zdcuvv/q6jj2Mnky9z7T5SHU1HBsYNiw8hWAzs7CjTx3M3t2\n8pSjhRIAOUZDQ27/v3Mnlw/JRzA5TIpKAMaO5dQvPxjDhblKiU2bgAkT0tcCknK3Y8Y4FsCuXYUt\nBudFJgEYHOQe0+bNhT2vXJEZzrZtY4urt5eFYMwYFoF0vv3aWt53woTyFID+fuC//is6AairS76u\nAwOFtQBy7chIurQKgE+2bvW3/44d3EvLd6pVmPzlL8AhhyTWknHT0wPccgsHwkQAZC6AKDNObAGo\nrOQG3772krddKiNk29t5uWmTM7ituTkx9TAVkhVVKgLwu98l59Wn4733eFmIujte1NYmlx/v7Ay/\nKingjDd4/vnUmXrpEAGorCzutqnoBMCdapgJubilVIxr7Vo2b9O5gHp7ORNCKnAScW8s6obVtj6I\nuHcojSjgCLi7aFyxIud+wAHAf/4nr69dm92EI9IwigDYI7aLjeXLgXPOcUQuG6Tn2hpRzp6XhXno\noYmZaWFhzwi2cqX//1cLIEf8mnfSINoTdxQ74mLwEoCBAf7r6eEHwJ6E5dvfBg48sPDna+MukDVl\nSmKmRqkJQEcH8ENXdSq/AjBuHAftRUCKkccf56W4vLJBxPEf/zH/55MNUSYZTJ/OyyOOSLZCsmHn\nThWAnJCe/DvvZFdEqdQEYO1aZ2J3LwE45hjg8MMdX39dXeJD8PWvF/Z83fzgB4mFukpdANrbgaOO\nStz2gx/4cwGNH88NxkUXFa+539bG1pqf8+vo4EFvhawDZBOlAOy7L/v/c50feGBAg8A5IQM8PvUp\nVl8779yLUhMAqSVD5B0Efukl4LXXnGnv3H7QQvg/01FRkdggTJ6cGPAVASiVSbWl0NkDDyRu92MB\n1NY616RYha+/35ngJ1sKXQTOjbvzIxRyJPyIEf5+0/5+TitWF1BAJkzg5QcfpI/El5oA2I15TU2y\neSnf1XYByYhgIPwBMH4ZPtxxK9x4o5PFVUwN4bXXOpUZhdZWHtEsjZw71zsbC0AyhLq7ncB8MX1v\nm74+tjoluJkN77wTfcaZ+/mYN4/rNBUKmSA+W959l0dOqwAExA6opctcEAFYty71PsXEMcc46zNm\n8DSDXg+lWADuB1CEsVior3cE4NvfdubRjdoCsI//298C992X+P5JJ3GQ3d3LXbiQl9lYAAC76044\noTQEwK8FcOedPBAxKrxcQIW2SvxaAGKh9/aqAATCK7PECxGANWtCPZ28MXMm8B//weujRvFAkw8/\nTNyHKNECEBYsKNx5ZosIgF1OoaEh+nEAI0c6gc/6el6+9przvsQtOjoS89xlYFg2FgDALrvTT3cE\nIGrhS4UIgB8LYNIk4OwIC7fU1PB5i1X82c+ypV/IcQl+YwDSGdqyhRv/UaO4WGCxpgkXpQAMDibm\n3mYSgD33dCbs3r7d301eKLZu5XPr73cGmQDsQ5eqn8KwYYkxAIDdQHYDViyIALz1lrNtzhynYFyU\niAiJAKxY4bwn90hdXWJcRQTAq/hbOkrBAqiv92cByD0YFcOGOSIAAE8/zctCWgDuNOdMXHABLzdv\n5iDwXnuJWex5AAAgAElEQVQBb7zhVAYtNopKAI44gpfbtzsPLZA+da2/n0355ma+UcaNA37843DP\nMxcmTgS+9KXkkYzjxyfPJ1tR4WQBSUM0fXr0AWAv6uqAn/wEuPdeZ9vs2cUhAHK9pHG2e47SELob\nEwnS+7W25HcqVgGQILCfzpFYoVHi5QYq5EDIsWO5PcqGlhaOAQA87qKqykliueSS3I4/OBiu9VBU\nAlBfz38tLfxgfvObvD3dBejvZzNtcNC5Wd0FpIqFESOSLQAvAbBdQA0NHAP5xCcKe67ZIg3pkiXO\ntj33LI6BeSIAUs/FbpzlvO0GcfNmYNEivn/8ZprcfDNw0EHF7QIqNQsA8A4EF5Jx47IfCOeOVVZV\ncU2vK67I/fjXXBPub1BUAlBby73jrVvZFP/P/wQ+//n0o1+9ikMVolZILowc6W0BuF1cUkhKGrBC\nTsfnFy/zeI89/Jf0CINhw7jGz2OP8bW33YrS8Nvut4YGvua5xFqGD2cBKMT39luf5jvfAf74R/9Z\nQMVqARQSPwLgPk/p6P3bv3GnLhfX9GuvhVtYsegEQCosSqNn+wC9cAvAnnsWX8BFMppkkgnbAvDK\nMti1K/qeV7bYAvCZz/Byjz2S4xqFRK73++87E3JPnJjYOw9jwFYhgt+bNrGw+bm+v/0tL/1YAIOD\nxWEB2GMB6uu50mwh8eMCkvMUV7Y858OG5T6gLGzxKyoBqKnhxvz00zn/X7Zdc03qB8stAPX1xVej\nXhrJ3t5kC8Au9SBUV0ff88qWSy7h2MbVVwM//zlvmzCBRcxP2YF8ItfTTv0cPTrRAgjDrdDQEL7r\na/VqXvo5zl578TLbNNDeXraEBgejjzuJBWAM31PHHVfY42cKAu/a5aQ+9/RwWvBNN/HrSmu2FbcF\nmi2xEgBxAQGOX3zHDk7xTOXX7+1NLCCXrsJmmBiTeryCbO/uTrYAbAEQS6CiIvqeV7bMmcM9zJ/+\n1AmESvA6KjeQCM+GDcB55/E6UeID6LfoYDZMmJAcz8k30iD4qeq522687O/Pzg2xdq3/8woLiQHI\nYLtCjgKW46drT15+GTj+eF4Xi0kSC0R4Af/ZREKsBKCmJvnBlAuQajBFW1uij7y6OhoL4NVXgVNP\n9X5PHtb//V+e4MK2AOwg14YNzuuoarAHQQap1dXxelRuIDG1u7udOj9uC6CqCvjKV7gsd76wR0WH\nhdwrdj2mTLS3A9/4BlfSzMYCKKZAtlgALS2OO6+QZHJB2wNWe3v5fGfM4ClEp01z3hs1yr8FYExu\nlUj9UFQCMHq00zjeeCMvpeGXG9/dO5HKmqK2UQlAR0fqHq/dKKxdm2wByPl2dgKzZvH6jBmhnGao\nSO+svp4FICoLQExygHu/H37IDWBHB4/0ffllvtZ33QW88kr+jmuPig4LeQ78WAA7dgAXX8ydimws\ngI4Ofp7OPTe3c8wnu3YBr7/O91IUo+AzBaElJXXnTuCee3j/iorkzuDIkewi9fMdJPhsp8Tnm6IS\ngHHjHEX953/mpS0AGzdygNFGBEBG00YlAH19yb2yDRs4aOXuQaSKAfT2OtM9ZjsStdgwhm/2iROj\nswB+/3tHSBsauDETC+C117hWSxj+7UIKQLaBSYDPacQI/r7ZWgBz5wJ3353bOeaTF1/k+NKPf5z7\n9IxB8LIABgeBF15w1gEeEf6//5s6djdqFA9WbWnJPqtHnp/dd3e2dXfnNyuoqARg7FhHUaXQli0A\n8sUvvhg47DBeFwEQohKA/v5kH9//+T8c0O7rA/bfP9FKEdwCUFfH6YS5DhwpFqK0ANaudX6Lgw7i\npR2EC2ucSH09uwJvuCGczwece8VP7StJlKiszN4CiHLeaS8eesiJ5xSSmhpug+xGd9Uq4MgjucPn\ntshSxe7sGEA2LjZjgH//dy7HYQvQ3nsDp53m/3ukojLzLoVj3LjklEgZA9DT46jt7bc777e1OQJw\n2WXAscc6w7ELSV8f/7C7djk9S+kN9PWx5fKTn3Cdcduks11Aknf93HOFPfcwiDIIvGMH3wtdXYmj\ngDs7nR7dr36V/+PK77p8ef4/W+jtZfegZMllgwiAHwugWGJQr7zCsYtzzwXOOKPwxydil21fX+Lz\nDPBzKh1Vsf5TJReMHOlkbrW0ZL6+3d3czlVUcJtCxBbBunX5HWVfVBbA+PE8b+kf/+hssy0Ar0Dw\njh1OEPjmmzkNa/NmZz7TQiFCZVsB0mtobuYbY599+Mez/ft2ELgY8q7zRSGDwD/7mdMrkmysK69k\nwRXEApg4ke+XMArriQDkw2e7cycnDLjp6fFfa6mvr3QtgEMOAY4+mqdDjSol1e0GEjefzNsBOPd6\nqkGodhC4pYXHc6Sb2rKvj3+v1193tk2cyEJQVC4gIjqBiFYT0ftEdHWKfW4Zev9tIkr56B1yCE95\n+NnPOtukN9/Tk9h7OfpoXrpdQPLwzZ2b2/fJFREAOw4gYvDaa07PYOrUxP9zu4BKJf8/E4V0Ab38\nMmddAHwtKyqSe2IiAH19nEsfxtSa+RSAK65ITCMUenr42mbjRvjgA+49l7IFAADPPhvt6H53IFi8\nFLYLSKoRpzpP+3p+9as8diZdeZfeXu4QSyxLyHf6ciABIKIKALcCOAHAfADnENE81z4nAZhljJkN\n4AIAtyd9kJyMx9n80z/xUGrbApg1iwOsRMkCEFUDKj0ELwH41a9S/3C2AEgBuHLAbx31INgPnft+\nEGwBCGMMAOA0/Pn4DZ9/3nt7dzcLQDbX9oILePKU/n7+zqVoARQDbgtArn1np/Ps/vKXvEwlAPY9\nuXJl5t9BOoPSns2cyZ1e+/PtCry5EtQCWAjgA2PMGmPMAIDfAXBnw58C4G4AMMYsBTCGiHzF8ydM\n4Cwf6b3cfrtT/nnr1uQH/u23gf328/tVguFlAdjr6QRAehflZAGEOSDPGGdELOB0HIxh0fXqvdbW\n8v3T2Rleb7KqCjj44HA+W2hpYRdiJgvAGEdEZHaqysrsspSKzQKIGrcF0N3N91BXl/cofi/cM85l\nejakoyKjiceNS6x60NPDVmzQrLOgAjAVgJ2PsH5oW6Z9docPPvc5DrgMDPCXnjXLsQb6+pJ7K+PH\ns8oWcoJu+UHtGEBbG6ccAqkFwI4BlJMF4DXfcb544gmeGlCQ4/T18QPhNVqUiO+TXbvCswAA4Mwz\n81NmQoKLbrZs4XpXmQTAq8JsRQU/F5myoDo71QKwqa1lq0h67V/5Ct9zHR3JheJSuf/cM6tlGhTm\n7gxWVvI2SYT52994+f3vZ/cdUhE0CyjbcIT7dvb8v0WLFv19vbGxEY2NjQA4PbSzk/2iVVXJN6fb\ndTR+PP9YTz0FnHxylmeYIy+9BFx4IXDWWfxaev3GcCXKyy8Hbr01OxfQli3hDvooJJlGUAbB3euR\nxrCjI7UAAHyt29oSa7TkG6/aTvlkyxYeYbpzp9Oz92LdOs446+lJnnEuU0kC9yxpcWfsWOD88zlw\na8e1nn6a/666ilNUt2zJPLJ84UJOFbYH8i1bBhxwQOJ+dtYRwPdsR4fcW034l39pAsBJM0EI+ihs\nAGANeMY0cA8/3T67D21LwhYAG+kNvfgiZ/lkKo9cXc0pY4XwQT/7LPeq3DGAzk7+0WRAVyYX0OAg\n8MgjyXPXliphCoBb8EUA2tv5N08lomH2/IX6+nAFoKWFXaKDg8Bf/sL56F5IOZH99nMEYP/9uVJl\npkFkagEkMn48tz2A0xZ99rPcwQSAT34SmD+f/9IxOMj/f8opznSlbW2ckdbZmdhxcdc4E/ddby8w\nfHgjjjyyEffey8fctGlxzt8tqAvodQCziWgGEVUDOAvAo659HgVwLgAQ0aEAdhhjcq6ZKBMtZ6K6\n2nETuXtA+URMsv5+NrFFALZt4zIE4tJJ1fgMG8bfp6KCyxe7o/6lSnV1eC4geQiNAX7zG84CImIB\nSGcBZHPfBCUsC0B6jHbtq3Qz34k70Z11ttdemQVALYBEpJiezVVX+f8cuW8nTHBSOSW28Pbbifu6\nXUAVFexukpL527bx9h/+0P952ASyAIwxO4noUgBPAagAcKcxZhURXTj0/h3GmCeI6CQi+gBAF4Dz\nA52wdcazZvFoWy+kAerv5/1EffONCEBvL/+wIgCtrXzjSMOfrkdlN5SFnO80TMK0ACS209fHpjkA\nTJniCECUbrSwBGDsWODPf+bPHj4cuPZa/v5LlgCf/nRybrgIgNtazqa+vVoAiXgVodtvP260P/jA\nf4lquzMosQB3yrQ7W01mBhw/nttAcS8feqi/Y7sJ7A01xjwJ4EnXtjtcry8NepynnmKzy+7F/cM/\nAN/+tvf+MvmKXU3UT/ZHeztnmixcmH4/CQw1N/P8nyIA0hCJ6GTbsIchUlFQiBiA7eITAUjnAioE\ndlZXvnnkEW78hw0DJk/mBkgGE61dy/NGCzKq3J0hN2ZM+kJy110HrFihFoCNV12u2lp2qe2/v//P\nk4KWdoVa9z1jWwAff8z39P/7f7w0huMIMhYqCEU1EjgdYobZFsAnP5l6f3EB2YOs/LBoEQd0Hn+c\nsy5SIRbA8uVcp0MCbG7BidsDVQgBsDNhpkzJHAQOc2o9IV8WgFdHYNs25z4aPZrvNRkde8wxifvK\nqPKZMxO3ZypYd+utvFQLwMHr2Q0ST7r5Zo4p1NQkThZl093tuI9nzHBESAb0SUwiKCUjAPJQiwVg\nDPDFL6beX1xAdo69H6R3uWSJM8rPC2lUVq9mARALwD3xy7hx/o5f6oQZA5AG1haAhgZ+mLZt8x4I\nBpSWAHjR2upYkqNH870mAmDXpV++nDNL6uqAE09MdC9kClKLRRu3Dks6vK5FkHjSzJnstqutBb7+\ndd7m/k3sGmdyvMpKdjnJc5WP2edKRgCkR5LthXcLgN/eqDwIYnGkEhD7wZs3z/Gv2hbA0qXAUEZr\nbKis5GuTy0TYmZAe7KuvOttGjWIBeOcdTn/0opQEQK7b17+eOIG9bQFs3ux9X++/P/CLXzg9SNuH\nXVeX3gKQc4+y9EKx4dXm5MNVW1npJKi425ft25NdT3Pm8FIEYPbs4OdQMgIgFkC2OdwSA/DrAhoY\n4EEWEmiUVK977vFuQGwB2Gcfx1qwLYCFC7O7YfI5O1XUEOXfDbR2LT8EMvjma19z3pNyux9+mJ8H\nI1fyJQCSwfbrXzsP/NKlzgC4gw/me23VKn4tlsBrrzmf4TXgSyyAlSsT7133ccslFpUPvErU5AOp\nZgAk3zNeJU1EEOSZuuKK4OdQtgIgMQC/LqD77mOfm1zk117jBuWCC5InpidyRuQBPHFDayu7j/wG\nnYmcSSbKBRGABx7IT4OyeTOnyt57b/J7tbXOnAyl6AI66qjEe9R2n9kiKlVMa2u5Z79pE19nKX1i\n14exGxj7/Lq7ubPy0EPJ72vPPxkR1zBF0d0+eQmA2wLIByUjAHJjZvsQ5xoDEB/+BmuompjSXj0m\n6YEB3FOYNIl9c+4YQDoeeIAfxnJ7+Cor2ZLKV2lu6Z3aZabHj+frLb+3e45om2IWgKamRJ+uXfrc\nfuDtRqGuDnjjDeALX3BiVnb5E68OhR0DsM/zjTfYwthjj8R7WuFUy2OOCbcctbt92rIlefzBrbdy\nhzOfVnXJCIDg1Qh7UVXFN7hMCp6tAEgv3+49iQ/2jDOSL77st2QJL6XqpB8L4PTTU08oX8pUVXGD\nlK+AoteN39PD11lmgou6kJkfAXjppcREBrnPBgcT/fT297azc+rr2Q30uc8llsMQ3BlAcn7y2XaP\n9uCD+T7s6SmfciT5YsoU4E9/Sqw/lW/se8YYngfAPWdFfb2T6puvciYlJwDZpl9VVvKUaoJfAdiw\nwQmeyTzEr7zC8xLbyOfKDyLBSD8WQLkiFoAIoZTMzRX7NxS/bFcX3xPV1ex+q69P3VMrhAVQW8sN\ndjbH+ugjbsCl4Zee/ssvc56/YA/csgVA7rkFCxyL9f330x/TtgDcLo3W1vIqSJhvXnqJxx7lay6J\ntjanfISM7AVYzLu7E+cCdpOvNN2yFQB3itSzz2Y3Q9X27c4ADKmz8vvfOw2O16xkgNPYiwD4jQGU\nI5WVfB2kxxk0bU16wnfdlVh3RSyAlpb0A+4KIQAS/M6mw7F1KzcC8r2kYd6wgX30Us5Z5r8GEq0b\ncVfOnMkdl+5uduPccguLiBdeWUByT48ZowKQjlGjgEcfTQy0B/08aTcefNBx33V1ZW7g041N8kPJ\nCUC29fJtHz7A0wZefHHm/2tvd7JI/vVfOfNk+HBn29y5LBJ2euNllzkjAsUFpBaAYwGIfzqo77Kv\nD/jHf+QZlexeckWFYwGke3AKIQBA9m6glha+39wC0NLC1mdjI2eG2SOe7e8nlkFlJZc7+etfOXj4\n+c+nLhFQX++MBO7v52v27rssXIODKgCZIMpvVpBYqxMmOCOE0w1mBPieuPDC/By/5ATAHu6eDpkj\n1sadxeNFW5vjO504kUvvAok9/48/Tnx9wQXOD6kWgIPEALq62I9qT5CTC3Z9lBEjEoNk1dXcsBWD\n/zpbAUhlAYgAAMmCZvt+7Vr08+Zx8DbTjGd1dc5cAT09wLHHsktj+nT+nSoro5t7N46I2M6a5YwJ\n6OpKLwCvvBJTAdi0yZlgJRN2UPXnP+el7WdLRXu703DbD5+dXbF2beJr2+0wYgT78NQCSLQApk7N\nrwAAievV1fz56Xqv110XfAKNbPBjAXR0OO4iWwBE3Oz6V/vum1jds6vL6Y3OnMkdE3cZYTf19U4i\nRU+P80yIAGjvv7D8/vc8XmO33RyLrpD1rEKcGiP/uGfVyYa993aq9WUrAAcdxD+Mjd3j37gx8fXE\nic665L6rBeDEANra2JIKKgDuErl+BeCb3wx2/GzJJACXX8732SOP8Gu5L+V/7HLMV17Jf15MnepY\nxGPHskWRjQUg9PQ4nRf5nGKwoOLElCn8J54DILMFkE9KygLIhYoK9s+/+io/aJlKE7S1AZdckuwv\ntuvOtLcnWgDuhqi1lR/uuFsA4gLavJkHsWSaBi8TmSwAY4qjAUsnALt2ATfdxCN8BUlOeOMNXkoh\nt0ysWAE88wyvjxnDPUh73lgvbP/19u1OQzNlinPuSuGxBSBTDCCfxEIAAK4cWlfH+bWpGBjgB8ir\nEZHG6+CDHR+/FCCzqanhSeuXLy/MDFTFjLiANm/miUjs7JNVq5ID9Zm46qrE38+2BqTRK4YGbNs2\n4OGHvd9zl2IeP94p1vbgg7zMNhA7erTTUIwezZ9TVZVdkLKxkceuSCqouJyK4frFkVGjnDZGLYA8\nYjfmXV2cHfHxx977trfzD+E15HvFCo5BnHdeYpDXPeiopsaxDrzqiMeJdAIwfz5XqvSLbXnZDV0x\nCcBHHwE/+Yn3e+LnHTeOBxfNmuU03HIvSS1/P4wZwymKqdKU3Rx4II8qFUtXRk/H3WqNipEjgeuv\n57bnjTcKZ8mWvQB4KenSpd77igB4MX8+xyDEVPvKV3jibTfiigC8p5KLE5WV3GNvbub6Su45mjNN\nTu6muhp47DHn9fXXA88957wHFIcApEMEYPJkLi8gPffx453rkUsqptesVemYPp0tFRFlafhVAKJh\n1CgnOH///WoB5A0vJU0VjGxryzxzlwhAqgkZamqcXlimyevLnaoq4Dvf4fVRo5IHIPnJyx8c5N6/\nfU0PPNAp9SHpkSIIUXLLLan98CIA553Hy9GjOT61xx7BBGDBgvQjR90cfji7R92p0fkqMaD4w05w\nqapSAcgbXhcyVc+zvT1zoy3BtlTYFkDc86ntxsRrJqp0denddHbyZ6Tyb8+dy43aOef4P898c9ZZ\nqTsS27dzTSnJ7Bk1ikf8nn8+3393351bOuawYdnnhhvDsaxx45LnolULIBqk3AzAFnOhBKCs9f6n\nP+U8ajepBCAbC2DixPQlJezAby7zhZYTIgCjRzuZMcawn1MGJLlTO1Nx/vmJmVhuKiq4VksxIBOv\ny3e1cU/0IR2OPfbgLKevfpVf+40BAP79xuPGcSD+E59g8QTUAogKEYDaWm6fdBxAHrj6au/tqVxA\n27Zl9qU2NKQfUSym/49+pNPqSWPy8MPcQ62pcapNErE1tXWrM9paaGriAXUHH+xskwyZUkAycbzG\ngrgFQDocEyY4kxgBucUy9t7b3/7Sy1y2zNmmFkA0SJnv2lruFKkFEAKLF/MgrlQWgD0EPxVjxyYH\nM23sUgVxRxoTuZnFDVRfzw3dbrt51wc66qjE9EiASx77TRuNEmnMvQTATg4QC2D8eJ4nVmIYuTQA\nJ53kby4Cr96+CkC0SAdAYwAhcN11wPHHewvAunXZCcCwYdxbS4VXGYm4Ig2MXAsZJS1zBQ8fnjy7\nkcRP3L7+qirgu98N93zziUxQY7N0KY9/sC0A6elPnAj86ldOMDDVrGaZyMV1ZKMCEB2HHsrlpoES\nEAAiGkdEzxDRe0T0NBF53rJEtIaI3iGit4joVa99Cok94s5m+nSeuSqbdLqGBl4++mjye2oBOIgA\nyM0sAtDXxw1kTU1yI/nAA7x0+87TpegWIzIl6d13O/7dQw8FHn888d6QEtm1tVziV2pdFWJO3vXr\nk7fFPXMtSl5+GTj3XF73m9abK0EsgP8L4BljzBwAzw699sIAaDTGLDDGLAxwvLwwalRyDEB6nW+/\n7U8ATj45+T0J3hRKwYsZ6d1KLEQEoL/fmcTF7QJas4aXzc3A737nbLfr45QCYgG88AIXD7RdM3Yv\n/eyzOV4kXHFFok8+TBYtAm68MXFboRoexRspyWEX/QuTIAJwCoC7h9bvBvD5NPsWoD+THV4WgOTt\nf/RRdg/ApEncwHulJKYq4xtH5BqIO0MEQEpoeLlJ7Nd2Smd7e2kKgFhBdvzCzhSbPRu45hrndV0d\nZ+UUgvPO48J0Nuncm0r4iAuwFASgwRgjczw1A2hIsZ8B8Cciep2IClSPMTWjR6cWACB7CyCVO0Ie\nIBUAxxoSd4btArJf27iDmFK8r5D1UfKBWwDWrXPiQ0H99GHxT//EFokSHZMmAb/9beHukbRZQET0\nDACvIszfs18YYwwRpRrXebgxZhMRTQDwDBGtNsZ4jqNdtGjR39cbGxvR2NiY7vRywssCsHud2fSA\nGhpS90al4Y/7IDCARwGfdJLz2t3ge1kA9gjr2loOmu67b+nNVCVZQCIAK1cC++3HdV6KVQBuuSXq\nM1CIgC99Kf0+TU1NaGpqysvx0gqAMea4VO8RUTMRTTLGbCaiyQA8h0cZYzYNLbcS0UMAFgLIKABh\nUV/v1OuXjAfbAsimgFs6C0B6u2pKsxgutKI+bgHwCgKvXOmsf+lLnBZZigLgtgDeegv4zGdYAOJe\nJVYJhrtzvHjx4pw/K4gL6FEAQxVNcB6ApAK4RFRPRCOH1ocDOB7A8gDHDAyRM2+v8PWvO+vZ9Nzn\nznXmCPbCmML58EqJmhpnOkLAOwjc3g7cey+vz5nDFSuB0hUACbJu2eJM5F2sFoASP4IIwE8BHEdE\n7wE4eug1iGgKEf1haJ9JAF4komUAlgJ43BjzdJATzge2G8gYTs0Dsh9tetBBTiOlZM8TTwBf/KLz\n2m0ByBiBs87i+VEnTeJsoB07+HcqpRz16urEWkft7U72mAqAUizkPBLYGNMK4FiP7RsBnDy0/hGA\nA3I+u5CwA8H2DGFeE8kr4eG2AKQY37BhPH/A+++zAIhbrhC58flCJqkXduxwBEBdQEqxEKuRwII9\nFsDtg1bCZ+pU7tH/8Y+J8/S6i/FNmJBcrbJUqK5OrBq7Zo1aAErxEUsBWLIEOOIIXlcBKDzSA3Zb\nXO5y3CNGpK+7VMxUVSVaAG1tTg0grbipFAuxFIC5c511FYDC8corvJSZj775zcSS2W4LYPjw0hWA\n2lquLmszejTwwx9quQWleIilADzxhFOfRQWgcBxyCC+lUfeKAbgFQGI1peY3nzjRKWsh1NQA3/te\ndpO2K0ohiOWtOHkyT/BujApAFEh2jHtcQFtbYu9YLID6eq6RXkpMmsSZTIpSzMRSAOrquBfW3a0C\nEAXiAnKPBHZbAJL2WYqjqidN4tpSilLMxFIAAB4M1tmpAhAF7uqggtsCEPxMHl8sjBvHWUAzZgDH\npRxPryjREtt8hBEjeDSwCMDnPhft+cQJqQ7qZQF4leIQi6GUkFTPo45KHGuiKMVEbAVALIBTTuHX\nXpO7KPnn9dedRt7LApgxI/l/7BG1pYKUrRg+HDjmmOSAsKIUA7EVgHfeARYsiPos4sdBBznrYgEY\nw6N8U836VYpZMyIA9fXA5z/Pf4pSbMRWAJToGTaMB0UNDLAYeMUALr20NOdWEBeQzImgKMWICoAS\nKeIGqq72tgD+/d+jOa+g2C4gRSlWStC4zg/2vMD2pCVKYbEDwamygEoRtQCUUiC2AjBqFHDxxbz+\nhz+k31cJj9paZ5BXqc37mw61AJRSILYCAJRmemG5sdtuwCOPAO++y6N+S9Hf74UIQEOqmbIVpQiI\ndQxABSB6Jk0CLrkEmDat9Gb9Soe4gOxid4pSbMTaAijFEablxqRJvCQqLwEYNoxdi5MnR30mipIa\ntQCUSJk1i5cDA1zzp5xq5WtygVLsxNoCUAGInuOP5+WmTSwCiqIUDhUAJVKmT4/6DBQlvsRaADQG\nED2aJaMo0RFrAZg/P+ozUMrJ568opUbOAkBEZxLRCiLaRUQHptnvBCJaTUTvE9HVuR4vDK66KrEa\npRINzz0X9RkoSjwJYgEsB3AagBdS7UBEFQBuBXACgPkAziGieQGOmVeIuBRBsdHU1BT1KRSUo47i\nALDXFIpxuxbp0GvhoNciP+QsAMaY1caY9zLsthDAB8aYNcaYAQC/A3BqrseMC3G8uSsrgb32St4e\nx2uRCr0WDnot8kPYMYCpANZZr9cPbVMURVEiJm0IjoieATDJ461rjDGPZfH5mmejKIpSpJAJmAtJ\nRLfjOm4AAAPfSURBVM8D+LYx5k2P9w4FsMgYc8LQ6+8CGDTG/MxjXxULRVGUHDDGUC7/l68kvFQH\nfx3AbCKaAWAjgLMAnOO1Y65fQFEURcmNIGmgpxHROgCHAvgDET05tH0KEf0BAIwxOwFcCuApACsB\n/I8xZlXw01YURVGCEtgFpCiKopQmkY8ELuaBYmFDRNOI6PmhAXXvEtFlQ9vHEdEzRPQeET1NRGOi\nPtdCQUQVRPQWET029DqW14KIxhDR/US0iohWEtEhMb4W3x16RpYT0b1EVBOXa0FEvyaiZiJabm1L\n+d2HrtX7Q23q8Zk+P1IBKPaBYgVgAMDlxph9wK60S4a+//8F8IwxZg6AZ4dex4Vvgd2FYprG9Vrc\nDOAJY8w8APsDWI0YXouh+OE3ARxojNkPQAWAsxGfa3EXuH208fzuRDQfHGedP/Q/vySitG181BZA\nrAeKGWM2G2OWDa13AlgFHidxCoC7h3a7G8DnoznDwkJEuwM4CcCv4CQWxO5aENFoAJ8xxvwa4Fia\nMaYNMbwWANrBHaV6IqoEUA9OKInFtTDGvAhgu2tzqu9+KoD7jDEDxpg1AD4At7EpiVoAdKDYEEM9\nnQUAlgJoMMY0D73VDCAuNTN/AeBKAHah7jheiz0BbCWiu4joTSL6LyIajhheC2NMK4AbAKwFN/w7\njDHPIIbXwiLVd58CbkOFjO1p1AKgEWgARDQCwAMAvmWM6bDfMxylL/vrRET/AGCLMeYtpEgrjsu1\nAKdnHwjgl8aYAwF0weXiiMu1IKKZAP4ZwAxwAzeCiL5s7xOXa+FFFt897XWJWgA2AJhmvZ6GRAUr\ne4ioCtz432OMeXhoczMRTRp6fzKALVGdXwE5DMApRPQxgPsAHE1E9yCe12I9gPXGmNeGXt8PFoTN\nMbwWBwP4izFm21Ba+YMAPoV4Xgsh1TPhbk93H9qWkqgF4O8DxYioGhzAeDTicyoYREQA7gSw0hhz\nk/XWowDOG1o/D8DD7v8tN4wx1xhjphlj9gQH+Z4zxnwF8bwWmwGsI6I5Q5uOBbACwGOI2bUAB78P\nJaK6oeflWHCSQByvhZDqmXgUwNlEVE1EewKYDeDVtJ9kjIn0D8CJAP4KDlh8N+rzKfB3/zTY370M\nwFtDfycAGAfgTwDeA/A0gDFRn2uBr8uRAB4dWo/ltQDwCQCvAXgb3OsdHeNrcRVYAJeDg55VcbkW\nYGt4I4B+cLz0/HTfHcA1Q23pagCfzfT5OhBMURQlpkTtAlIURVEiQgVAURQlpqgAKIqixBQVAEVR\nlJiiAqAoihJTVAAURVFiigqAoihKTFEBUBRFiSn/H/dCl1gKBXasAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c46a190>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Comparison"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tout = np.mgrid[0:100.0:.1]\n",
      "x0 = 0\n",
      "dtmax= .001\n",
      "\n",
      "%timeit sdemkl.sde(x0, tout, dt=dtmax)\n",
      "%timeit sde(x0, tout, dt=dtmax)\n",
      "%timeit sdepy(x0, tout, lambda x:-x, lambda x:1.0, dt=dtmax)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 988 \u00b5s per loop\n",
        "10000 loops, best of 3: 171 \u00b5s per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1 loops, best of 3: 445 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "GSL is far and away the winner. It is around 3000x faster than the pure python code."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Theano way to calculate SDE"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import theano\n",
      "import theano.tensor as T\n",
      "from theano.tensor import sqrt\n",
      "\n",
      "from theano.tensor.shared_randomstreams import RandomStreams\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "def A(x):\n",
      "    return -x\n",
      "def B(x):\n",
      "    return 1.0\n",
      "\n",
      "def onestep(x, dt):\n",
      "    return x + A(x) * dt + sqrt(dt) * B(x) * rv\n",
      "\n",
      "srng = RandomStreams(seed=31415)\n",
      "\n",
      "x = T.dscalar(\"x\")\n",
      "dt = T.dscalar(\"dt\")\n",
      "\n",
      "rv = srng.normal()\n",
      "\n",
      "xout, updates = theano.scan(onestep, non_sequences=[dt], n_steps=100, outputs_info=[x,])\n",
      "\n",
      "sde = theano.function(inputs=[x, dt], outputs=xout, updates=updates, allow_input_downcast=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/noah/epd/lib/python2.7/site-packages/theano/scan_module/scan_perform_ext.py:85: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility\n",
        "  from scan_perform.scan_perform import *\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sde(0.0, .01)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "array([ 0.09437863,  0.24098799,  0.16236896,  0.33597058,  0.39227402,\n",
        "        0.42859804,  0.4297515 ,  0.41065422,  0.42676154,  0.46194971,\n",
        "        0.23951774,  0.22655675,  0.22447993,  0.2832613 ,  0.30663556,\n",
        "        0.2521183 ,  0.15739451,  0.20163056,  0.44691499,  0.49963293,\n",
        "        0.44309233,  0.45447732,  0.76954386,  0.99497569,  0.94931493,\n",
        "        1.00052418,  0.94978512,  0.87785128,  0.94674234,  0.99830616,\n",
        "        1.01836497,  1.13322821,  1.24859349,  1.10906753,  1.1890498 ,\n",
        "        1.14122232,  1.09423422,  1.08703787,  1.01070417,  1.0176151 ,\n",
        "        0.93632021,  0.97563049,  0.87088073,  0.697739  ,  0.73428342,\n",
        "        0.71105801,  0.7295656 ,  0.77290335,  0.87487668,  0.70402478,\n",
        "        0.84913249,  0.92544088,  0.86158857,  0.85057093,  0.9961133 ,\n",
        "        1.07968195,  0.97135231,  0.86316631,  0.7943823 ,  0.72409596,\n",
        "        0.92205431,  0.8814426 ,  1.0044605 ,  1.00126224,  0.90335298,\n",
        "        0.82920857,  0.70217597,  0.7704939 ,  0.69542109,  0.56746137,\n",
        "        0.60785471,  0.58045189,  0.6764967 ,  0.63270697,  0.66330246,\n",
        "        0.72481229,  0.64581193,  0.7229394 ,  0.72425885,  0.82994694,\n",
        "        0.78466216,  0.76627463,  0.83304857,  0.68643934,  0.58504504,\n",
        "        0.42661802,  0.51153621,  0.4141255 ,  0.37583259,  0.35615542,\n",
        "        0.24124443,  0.11700228,  0.0925223 ,  0.01597322, -0.04461747,\n",
        "       -0.01722394,  0.0083655 ,  0.17626539,  0.09930939,  0.13939293])"
       ]
      }
     ],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}