{
 "metadata": {
  "name": "",
  "signature": "sha256:9b691bbd41c05ac28e17734370358c3d0ae3bba31bb28a91f38cd5bd6bfcf5c3"
 },
 "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 0x10d002a90>]"
       ]
      },
      {
       "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 0x10c458110>"
       ]
      }
     ],
     "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_MT19937 \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_MT19937, 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/Dropbox/gnl/ipynb/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 0x10d02b290>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmcXVWV738rqSSVVCWVVCWpzCQhYQhDmAwIImUEjCIg\njbTQYAMtgwpPFFsQJxJ92tDiQxDpBxp9YCugCBiGAGEoiC1DkJkkZGLIPFZqnlK13x/rLvc+5547\nnnOHund9P5/6nOGee86pc+/dv73WXnstMsZAURRFKT8GFfoGFEVRlMKgAqAoilKmqAAoiqKUKSoA\niqIoZYoKgKIoSpmiAqAoilKmhBYAIvoNEW0norcSvN5ARM1E9Frs73thr6koiqKEpyKCc/wWwC8A\n3J3kmOeMMWdEcC1FURQlIkJbAMaY5QCaUhxGYa+jKIqiREs+xgAMgOOJ6A0ieoyI5uThmoqiKEoK\nonABpeJVAFONMR1E9GkADwE4IA/XVRRFUZKQcwEwxrQ660uJ6HYiqjXG7HGPIyJNSqQoipIFxpis\n3Ow5dwERUT0RUWx9HgDyN/6CMUb/jMH1119f8Hsolj99Fvos9Fkk/wtDaAuAiO4BcBKAsUS0EcD1\nAIbEGvQ7AHwewFeIaB+ADgDnhr2moiiKEp7QAmCMOS/F678E8Muw11EURVGiRWcCFyENDQ2FvoWi\nQZ+FRZ+FRZ9FNFBYH1JUEJEplntRFEUZKBARTLEOAiuKoijFiQqAoihKCo44AlizptB3ET0qAIqi\nKCl44w3ghRcKfRfRowKgKIqSBqU4RKkCoCiKkgb9/YW+g+hRAVAURUkDFQBFUZQyRQVAURSlRNm1\nC/jBD+L39/Z6l6WECoCiKAqARx4BfvSj+P2dnbxsb8/v/eQDFQBFURQk7uGrACiKopQ4iQSgq4uX\newKT2A9sVAAURVEA7NsXvF8sgG3b8ncv+UIFQFEUBaldQNu35+9e8oUKgKIoCoIF4J13gKOOAsaM\nibcAtm8HiIDGxvTO39ER+hYjRwVAURQFwQKwYwcvp0+PF4D163n5wAOpz715M1BVFer2coIKgKIo\nCoDnnuOlm/OHYln2hw4Fenq8vfitW3mZTnSQe0xvb/FYAyoAiqIoAFau5KU7GNzWxsvdu4H6euD4\n44ENG3jfhg3sGpJjkjEo1tL29wOXXMLnKgZUABRFUQAMHszLnh67T3runZ1AZSWnhV66lPf96U/A\n2WenZwHIOTs7gddfT0808oEKgKIoCuwYgCsA0lB3dLAVAADvvcfLlhYeIE6nMe/utufp64vmfqNA\nBUBRFAW24W9ttftcC0CSwbW08LKjAxg/Pj0BkHN3dHBkUbGgAqAoigLbSJ9wgt3X1MTLrq5gARg3\nLj0XkFgAjz9u9xVDgZnQAkBEvyGi7UT0VpJjbiWitUT0BhEdGfaaiqIoUSMuoE2bbITPAw8AX/4y\n8Mc/2ga7uRlYtYothQkTrCAkQwTgy1+2+773vejuPVuisAB+C2BBoheJ6DMAZhljZgO4DMB/RXBN\nRVGUyDCGLYDKSt6eNImXu3YB3/0ucM45VgBaWoA5c/j4SZOS5wjq7WXrQQTA5Sc/AZ56Ktr/I1NC\nC4AxZjmApiSHnAHgrtixLwEYTURFEgSlKIrCoZ+DB3Pv3qW5Gaip4XW/CwjgyV3G2HQRfi64gCeR\nuQPL8+bxcsIE4JRTgGnTIvkXsiIfYwCTAWx0tjcBmJKH6yqKoqRFby9P9ho61O7r7uaGXWbwTp/O\nS3dGMBFQW2sjhPysWMEpI7q7gYkTed9BBwGrVwNbtvD2xo3B780HFXm6Dvm2A4c/Fi5c+I/1hoYG\nNDQ05O6OFEVRYvT0AEOGePdVVnLvXyZx/fWv7BKaPdt7XF0dC8CUgG6tTCrr7ubjtm7lyWMHHug9\nbudOHlBOh8bGRjSmm4AoBfkQgM0ApjrbU2L74nAFQFEUJV+IBeBH3D8AMHo0UF1tt489lpejRiUO\nBRUBePZZFgCALQY/48enHxXk7xwvWrQovTcGkA8X0BIA/woARHQcgL3GmBJMrKooykBl2TLu3fsZ\nP967XeF0mV94gZeVlYnHACSy6O67gcmTeT1IANxj80loC4CI7gFwEoCxRLQRwPUAhgCAMeYOY8xj\nRPQZIloHoB3AxWGvqSiKEiXnnx+8f9asxO+RRHGVlbZqmB+3UT/sMHu88MYbwNy5vN7amlgcckVo\nATDGnJfGMVeGvY6iKEq+kXDQZAwfntgCcNM+7LcfL92IoMMPt+stLfkXAJ0JrCiKkgCJ3HH55je9\n28ksAHLCXyZOBNas4WygQaQzoSxqVAAURSl7FiywUUDDhtn9hx4af+xNN3kHbJNZAIOcFnbiRI4g\nChpsBlQAFEVRsqavz8bWZ8rYscDixbwudQFuvpmFIRXJLAB30DjImgCAo4/m5Zo16d1rlKgAKIpS\nEtx7r420yZTubtsznzmTl+nG5Q8fnlgA3DKQI0cGH/P008DChRyJlG9UABRFKQmkgU3kjklGd7fX\n9QOkX8M3WRiozAO4/HLveIBLTQ2nh7j3XuB3v0vvmlGhAqAoSkkg/vZkydlc7rkHePhhXu/p8frm\nX3wROP309M6TzAUkSeBmzEh+Don+Wb48vWtGhQqAoiglgcTcB2XeDOJf/gX44hfte1wL4NhjbYnI\nVCQbBO7uBn7/e+Cyy5KfY8wYXv7qVzw3IF+oACiKUhKIACTqja9cCZx8snefRPMEuYDSJZUFcPbZ\ntoFPhJty4ogjsruPbMhXMjhFUZSc4gpAXx+weTOnWr7jDmDdOmDqVB5wdREB8LuAMiGRBdDfnzjH\nkJ/x49lSSDQjOVeoBaAoSknguoB+/Wueefvcc1yF66abgt+TSwtAIosSDf66ELFL6oUXbJK5fKAC\noChKSeBaAJLYbefO9N7b05O9ACQKA/37373zANJh9Ghg797s7iMb1AWkKEpJ4AqA9Oyl951oQNe1\nALJ1AQWFgW7eDJx4YubnyrcAqAWgKEpJ4LqApGFvbwfOPdfW/JXjpCBLFC6gIAvgww+zO1dNTbAA\n9PWxmElZyu9+F3j++eyu4aIWgKIoJYFrAUhD2dbGETiVlVyaUfZJ2oWODuDqq8MNAgdZAMcfz8uL\nLsr8XPv28f/iVijr6ODl4MEsBj/5CbB2LfDxj2d3z4JaAIqilARB8wDa2riKV02NzRPkn227dGlu\nBoGrqoDf/jazcxHx+9rb7b4VK4BNm+y2WAj+EpbZoAKgKEpJIAJw0UXWtSMCUFVlLYCrrvK+b/Lk\n8C4gvwUwZw7wt79ld77qaltisrWV00RccIF9XQa4Mx1gDkIFQFGUkqC31zbi4gJqbeXGf8QIYMeO\n+ILuAHDQQfEul0wIsgB6eryVvzLBFYArY6W0Xn3Vvi4CoBaAoihKjN5eHvA94ADbIO/Zww2qCMC0\nad73/Oxn7G4ZMsSbuz8Thg6Nr+fb1RWNAOzeHf/6jh28VAFQFEWJ0dvLA75tbdYls3GjFYDt260A\n3HYbsHUrH9/cnL37B2ABcMs8AuxSikIAgjjrLF6qACiKosQQAWhttQLw7LNWAABOBwFwGOiECdxI\nt7RkHwEEBAtAV1f2ojJyJIuSnxkzgDfftNsqAIqiKDG6uriyV3s7/0kCNhkEBoDp03kp2yIAxWQB\nzJgBbNgQv7+jAzjsMODtt3lbBUBRFCXGtm3ApEncmO/Zw5E4AAuANMZHHslLvwCEsQAqKjh2Xwae\n+/vDzSs48EDg3Xd5XaKZ6us5GggADjkEuOEGng8QFhUARVEGHNdcEz/wum0bu3WqqjgHUH0976+u\ntrNm587lpQz4VlaGHwMg8g4ES0hpOknggnAFoKODM5iuWQM89JA9prIy3urIhtACQEQLiGg1Ea0l\nomsDXm8gomYiei32972w11QUpbz56U95gNdl61YWgOHDgaYmXgdYEO68k8cDiPh9hx7Kr1VWsnCE\nEQDA6wYKM6cA8ApAczPnBxo1yhulNHgwcOut2ZW/dAk1lYCIBgO4DcDJADYDWEFES4wxq3yHPmeM\nOSPMtRRFUVzcHvAll/Bs2bFjbdqH4cP5tYkTre8fAKZMsetiJaSbNTQRrgCECQEFeGLatm3s/tm7\nlwXAT0sLL1f5W9oMCWsBzAOwzhjzvjGmF8C9AM4MOC5LY0hRlHKjoyO5f1t87e7kq8WLeVlVZf36\nRx8NnHBC8t74zJm8DNuTjtICGDzY5hfau9dbLUyQuseS0yhbwgrAZACuIbYpts/FADieiN4goseI\naE7IayqKUkKsWcMF2oWqKuDb3058vPja3Xw5ApHtfX/kI8Bf/5r6+ldcwX9hiNICAHjcorWVhSxI\nAJqaeJlsvkA6hM0mYdI45lUAU40xHUT0aQAPATgg6MCFCxf+Y72hoQENDQ0hb09RlGLn6quBRx8F\nzjsPeP113rduHS+NYd93R4d16YgAtLbyUtwhgvS+JdInFbfdlv29C34LIAoB2LmTI4z8OX8aGxtR\nXd2Iqirg/vvDXSesAGwGMNXZngq2Av6BMabVWV9KRLcTUa0xZo//ZK4AKIpSHriTnm65hZcSQim9\n/NZWKwDS0Erv98UXved75x1eyiBwPujtBRYtAv7jP7gMZdhB5epqTgMRdB7pHI8YwZPGnnhiUdbX\nCSsArwCYTUTTAWwB8AUA57kHEFE9gB3GGENE8wBQUOOvKEp5Iu4MwE5uEgGQxGdtbVw4HbACIBZA\nVxdPnvrkJ73nyzYMMxs++ID/2ts5XPNjHwt3vupq9vMnExJ/2uhsCCUAxph9RHQlgCcADAaw2Biz\nioguj71+B4DPA/gKEe0D0AHg3HC3rChKKeH6sUUAZCnROe4xIgBSdauzEzjmGOBXv7LHiLWQb8SF\nFdYC6OjgsNVUAiACmS2hM0obY5YCWOrbd4ez/ksAvwx7HUVRShM3nDOZBSDIGMDKlbzs6opv8CX3\nT7649lrgxhttMZqwYwD77w/88Y9sCSQiCgtAZwIrilJQpNE8/XQb/vnWW7xMZgFIDPxFF8UnT8u3\nANxwAw/WyoB02Dw9l1zC58u1C0gFQFGUgiIC8Mgj7NefO5eraXV2WgugtdUe39PDPeS1azkHD+CN\n4//Wt4Dvfz8/9+6yb59tkMOOP4waxeKXTAAmTuRJY2FQAVAUpWAY463h+/LLtlFbutRaACtW2GN6\nenh2bE0ND7wC3gb3P/8TuPTS3N53Kkw6AfJJqKnh/zOZAMyfzykxwhBBVUlFUZTs2LmTGzupfLV2\nLXDKKbx+9tncAJ55prckomTaHDfOhnxGkRgtSsIKwKhRvAw7mJwKtQAURSkYGzfaIi0A5/Cf4+QK\n6O7m5Ghu2ofeXhaAt99mcQCKRwCujaXDDDsGoAKgKErJs2mT149dWwtcfrn3mIMP9gpAd3d8rv1i\nEQD5X2prw51HBOCpp8KdJxUqAIqiFIz2dp7NOmsWbwf1eA880DvI29TEloKbv+cLX8jtfabDm28C\nX/4yr0s1smwZNIifiSS+yxUqAIqiFIzeXnaXrF3L23sCcgTU1XktgF27eN/XvsbbV1wBfPObub/X\nVBx2mHX9RJHG7JprgDNynERfB4EVRSkYIgAAjwWcfXb8McOHewVg927O+y+Tv447Lvf3mQlhB4CF\nSy/NfTSTCoCiKDljwwYO57zoouDsnK4ArF/vzXw5YwZw3XU2N76waxe7R2SyV1C6ZCU9VAAURckZ\n//IvwEsvcYnDW2+Nf90VAH/kzKxZ3ANubbUWwIsv8mzbUaOsAARVzFLSQ8cAFEXJGZKZc8eO+NfW\nrOHZs4lCJmVyl7iAjAE++lHg/fd5n+TbkYgZJXNUABRFyQk33mhLFooLp7ubG/auLhvfHyQABx5o\nJ4RVVHCOoGuu4e3du7nxF4HIdax8KaMuIEVRcoJb1lFy5GzZwsv33uNlU1N8TD8ArF4dv++mm3i5\ndavt/Z98MjBtWjT3W46oBaAoSs7p6ODlpli9QOnd796d3qzZP//ZDiLv2WMFYNmy/Gf+LCVUABRF\nyQluLnsRgK1bebl5My/TFYBp07ypjwtV8KXUUAFQFCUn7Ntn0znLsqODG/z6et7esyc9AfA3+GEL\nriiMCoCiKJHT1cVpDGpreQB4xw7e7unhfDl79/Jx6VoA0uB/5CPebSUcKgCKokSO5Osh4sa6qor3\n9fRw2KbUANixIzMLQAZ8VQCioagEQIo/KIoyMHn0UeCII6wACBMmANu3swCMHGn3p2sBiABI6mgN\n/YyGohKA8eMLfQeKooThzTeBN95g374rAPX1wLZt3PN3BQDIzAU0diwvg0JHlcwpKgFQFGVgI7V7\nX3nFGwUUZAFIbz4TARg1iq8RNt++whTdRDA3N4iiKAMLKe34jW9494sF4ArAYYdxRbB0fu8y63f2\nbK+wKOEIbQEQ0QIiWk1Ea4no2gTH3Bp7/Q0iOjLZ+aQHoSjKwOPRR3k5eDDwxBN2v1gA3d02d8/h\nh/Ny+vT0zv3aa8CnPhXZrSoIKQBENBjAbQAWAJgD4DwiOth3zGcAzDLGzAZwGYD/SnbO5uYwd6Qo\nShASdZNLtm61E7z6+ry1fidP5txArgDMmsU1dOfNS+/8RxxhLQElGsJaAPMArDPGvG+M6QVwL4Az\nfcecAeAuADDGvARgNBHVJzphS0vIO1IUxcOqVexD37Ytt9f58EPgmGPstpulU0o2bt1qXUD19cAN\nN3D5Q6UwhH30kwFsdLY3xfalOmZKohOqC0hRomXVKl5++GFur7NxI/f6f/xj3nYLtUjY5p//bIXh\nhBNyez9KasIOAqdb/MxvuCV430LccQfw1FNAQ0MDGqIorKkoZY5k3sx152rHDu7Vn38+8N3vBlcA\nAzh5W3Oz5vHPlsbGRjQ2NkZyrrACsBmA4+nDVHAPP9kxU2L7AliIc87JfSFkRSknxPUTpXvVGE7t\nMHiw3dfayu6dyZOBq66K99c3NnKx9KFDtfEPg79zvGjRoqzPFdYF9AqA2UQ0nYiGAvgCgCW+Y5YA\n+FcAIKLjAOw1xmwPvJlB3uLPiqKER/LuRGkBXH+9TegmiABUVAA//3n8e+R4bfyLh1ACYIzZB+BK\nAE8AWAngPmPMKiK6nIgujx3zGIANRLQOwB0AvprofOPGeYs/K9HT3Jx7X7BSXDQ18Sz7KC2Avj6O\n+ZeKXwALQLIYfZm8dfDBiY9R8kvoiWDGmKUAlvr23eHbvjKdc40dqwKQa/71X4ElS9iEV8qDpiZO\nohalAEgkzyuvAAccwOttbfFpHlxqa3k+gBseqhSWogrAmjNHBSDXbA90vimlzN693Oi6BVXC0tPD\nSxlgBqwLKBEVFRwGqmGfxUNRfRQzZ+oYQC7p6PAKQHu7TqwpB1pb2QUUpQB0dwNTpgDvv++9TjIB\nUIqPohKA4cPVAsgll17q/cFKgyC9OaU06eoC6upsWcYo6OkBDjwQeOABYOVK4JJLUo8BKMWHCkAZ\nIYm6BGkQdPJdadLTw8kVRQCitgBmz+a0z5//PLB4sVoAA5GiE4DWVuChhwp9J6WJv3cmDYKm3yhN\n5s8HTjqJG+va2mgtgO5uHrMDuPYvwN8jFYCBRdEJwP33A2edVeg7KU38AiANggpAabJqFfDCC2wB\nRC0APT080/e88+xA8M6dKgADjaISgMrKeDeFEh2uACxZYicIqQCUJhMm8LK3Fxg9OnoX0LBhnNZB\nLID2dhWAgUZRFYSRup8ATzPXcLFoGT3arp95JnDccbwuQlBINm/mDkBdXaHvpHQYN46Xw4Zxbz1q\nF5AIgIv7G1aKn6JqYt0vjw4GR8/gwcAnP2m3N8WyNmU6N6C1Fbj55mgnk02ZwvljlOiQMopA9gJw\n663AaafF7+/p4Zw+fgHQsOKBRVEJgPuFDTJXR4xg14WSHb29XFRDEAHYnCA1XyKWLgWuvpojQKIU\ngb6+6M6l2MLpXV3828nGBfTqq8Bjj8XvT2QBKAOLohIA1wII+rJ2dgIvvpi/+yk1ent5QtCXvmT3\nnXNO5gKwZw8vx44Frrsu/H2JiETpolCsAADZWwDuOVyCBECTvA08BpQAAHbAScmc3l4uwC2+f4DD\nBNMdA9i0iRuRpia773e/C39fUq6wGMYiSokKZ4RvxIjsBKC3N3h/ezsHFYgAzJ4NfPazmZ9fKSxF\nOwic6MuqboLs2bePG4UZM+y+/fZL3zUwdSonk3PTAG/ZEv6+JAqprS38uRSLdJbOOIPdq93d/Ptx\nc/inIpEAtLWxVSFzAd55xys4ysCgaC2AY4/1viZfZhWA7BELYPp0u6+62ja8mzenLh6+ciW7gI48\nMrr7EgHQPFDR0dtrG+8LLuDB2REjMg+ukM/EFZPXXuPvTHU1W5DPPMPfKx0AHngUlQBI2FoQl13G\ny/7+/NxLKSICMHOm9ddWVXF+oN5ejsS5/vrk59i9m0v/nX56dPfV0sKNiQpAeH70I2DRIvbdv/IK\n75N6vNkMBMvx0jF4+GEOAhAXEBHwiU9Ec+9K/ikqAUgWQ7xyJS/VAsiefftsT00qylVXc4GYm27i\n7R07kp+juxv44APgmGPsvrCRQO3tHP+fyvpQkvOjH7GAL1zI21u2cG1eCePMZiBY8kS5n01FhXUB\nKQObohIAP27DIrNYb79di5lkS2+v9dPKD1p+xNIwpDLj+/tZMI4+2u4LO2ejuxuoqVELICw/+EH8\nb2P+fOvzHzkys1nfxgCrV/N6XR2waxev9/fz9yRRhJAycCg6AWhtBWbN4nU3TbGYsYAtcq1khriA\nAPtsRQDE/ZZq9nVvL39GEyfyJKHq6sQDhemiAhA9l1zCS/m8AR7Ez6Qc6KZN3oHdv/6Vl3v3atrn\nUqHoBKC62vZi3J6lNEwf+xjw7rv5v69SwBWAxYt58E7GAkRgU1kAe/dyY00E/K//xb3AKARg1Ch1\nAYXBDY9essR2olwBmD7dWw8iFRs22PMANknj448Dp56a7Z0qxUTRCQBgBcDtEU6bxsuqKk0TkS0S\nBgpwKOgnPsENxJVX2sY3kQDI2Etfnzen0JAh4edmdHez8Pf16RhPtriunaoq20MPIwDvvecNGRaa\nm4FJk7K5S6XYKEoBkN6o3wK45Rb+QoftcZYrrgXgUllpxfbOO4Pf295uxaOmxu6vqIjGAhg2zMaq\nK5njCkB1dTQCsGsXzxz309ys/v9SoSgF4PHHeelaANJIqABkT3e3N9+S4G94g55vWxunfhg8ON4C\niEoAhg2zn/nChTYySUkOEbBiBYf3AtxZChKA/fbzFnFPRU+P7Yy535vmZu+YnDJwKUoBmDaNk5a5\nFoBfAO68kyOClPTp6gr+4Q4b5g0PDCoR2dbGUSQjR3otgKhcQGIBiAD88IfAt74V7rzlgETpvPAC\nMHkycMMNPDtXBMCdWzN2rDeNRyok4yfAn8v48ez77+9XASgVshYAIqolomVEtIaIniSi0QmOe5+I\n3iSi14jo5XTPX1mZXAAuvzyaRGTlwgcf8I84kQXg9gznzYt/tlLwe+RIrwUweHD4FA7y2dbWclUp\ngBszJTnGAAcfzOurV3Njf+21POFLrDK3vkJVVWYTwUQAZH7O9u128FddQKVBGAvg2wCWGWMOAPB0\nbDsIA6DBGHOkMWZeuicfNcrr15Qvo+ty0FS06bFtG/t/E7mAhg0D7rvPbq9fzz1JF7EARo3yCsDq\n1d45AdkgAjBzphUiHWRMjYglALz1lre3LwLqDuq7aT/SwS8AgP3NqQVQGoQRgDMA3BVbvwvA55Ic\nm3GWkNpar7mqYwDZI5E1bW3BP1wRhf/9vxOfQ3K/+F1AURAkABKeqhlCE7NqlV3ftIldPMKRR8ZP\nChs+nK3AdCOtRAAOPdS6lORzUQugNAgjAPXGGKkltR1AfYLjDICniOgVIro03ZPX1tq884BtJIwB\nvvzlbG+5PJFJX1J20Y/08CQ9hGAM8LOfARs3coqIsWPjXUBRIJ/tjBkcew5YV8W6dcEFSRQWgLPP\nton5amuTHz9oUGZpoUUAHn2UXYiAFX+1AEqDpAlciWgZgAkBL33X3TDGGCJKlKDhBGPMViIaB2AZ\nEa02xiwPOnChJDEB0NLSgKamhn9sSyPhuoU0JUR6JJpRLUiDPnIkcOON7EcGuPf97/8OPP88v+/Q\nQ7lhDrIAMk0z7CJjEzNmAMuWccK5F17g6/z858Dvf6+fdRDr13NOplNO4U5ROgVZZBwgWfF2CRfu\n7mYBcGf9qgAUnsbGRjQ2NkZyrqQCYIw5JdFrRLSdiCYYY7YR0UQAgWnEjDFbY8udRPQggHkAUgrA\n//k/3mnr0htx3T/NzdwwaBra5LghnkEWgCsAp59uBWDjRh5ElDKcF10E/NM/Bfvnr7mGrYVsEPfS\npEk8XvH667z/wAMzS11QbkipR0nnkY4ApBoHePdd4KCD+HflRgEJIgDqAiocDQ0NaHDM9UWLFmV9\nrjAuoCUALoytXwjgIf8BRDSCiEbG1qsAnArgrXROPm0aT/z6xjd4WywACTncu5e/8PnKC7RgAYvS\nQMQVgKCe25gxvBw50pvhcdMmm/XznHOAz3yGrYAgV4OEI2aDCMDw4Rz51dMDfOpTwCGH2LrFSjy9\nvdwQZyIAqSKB3DGXZAKgFkBpEEYAbgBwChGtATA/tg0imkREj8aOmQBgORG9DuAlAI8YY55M5+Qy\nqeW//ouXIgBiAdTUsMsgX3mBnngCuP/+/FwralwBCHLTuALgRlbt3m3DCFO5d8I0CCIAMg9Atuvq\nMq9XXE5IAy0CkMytI7hzLYKQz3H9+mABmDiRl5qOpTTIuoibMWYPgJMD9m8BcFpsfQOAI7I5/0EH\ncY/w5NgV/AIAcI/0/vvjBy9zRVAahYFAdzew//78ow5CXEDDhnkjRNrbbeOSarKXO2CfKRJiKhaA\nKwDu+IXipaeHv5OZWAASCZQISbo4axaPLfgFYOhQHnQOyhGkDDyKciYwwD3RP/zB+veDxgBmzMhv\nmGC+ap5edlk0tXaF7m4u2v1WAufbiBF2kHX4cODb32bBkKpPS5cm9+9PmGCjRLLB7wJqaWFBcCcx\nKfHIb+Lr0CJHAAAgAElEQVSQQ3ieRzqT51JZAO7v6803gy27V1/1FgRSBi5FKwCAt4KRfwwA4O18\nJA9bHhuyzpcF8OyzmSXtSoU8u0MPTX0sEUeU9PZaC2DBApuNNYj99mNXTTaZPJubOcRUBKCrC/j6\n13lMQQUgOTIGMHo0z59IFQYK2Bn2fX02v7//nML27TrYW+oUtQC4MctBLqB8CcAzz/AyXxZARwc3\njGHp7ubetITzpcvQody7dF1AiRg1Cvj4x7lhCcohlIpZs1g8qqv58xSXz5o1XgHQMNB4xAWUCSKy\njz4KnHhi/Ot+V19Q1JhSOgwYAZDMhB/9KPc4AW+DkUukIcqnAGRSui8RV1zBk7fcrI7pIK62dASg\nqYnnDoj7JlOkzOCYMWx9yH0+9JB3vkHYhHMDif5+4MILgd/+NvlxQYO0qRAXUKIxG/8se7XCSpsB\nIQBuTPItt1j3SL4sgO3b7fXygfjBw7J5M/+gxXpKl0wsgEGDuOFOJgDbtvF4DsDzBRYvjj9Gepoy\nK3n2bO89l9Ng8J13AnffzZPgEtHTw887UwEQCyCRteYXADe9hFJ6FLUASMxyby+HIfrr1eZaAMSn\nLTly8jEG0NfH/9P//E/46lhiPa1fn1nivEwEQEgmAD/5CXD++bz+058CP/5x/DEy2C8CIOmhhXIq\nFLNxIy+T1WeeO5e/I5l+J2UMIFEHY98+rhQ3fz5vpxNaqgxciloAamo4yidRD3bo0Nw1DCeeyC6f\nJ56w+/LRCEkjetddwIMPhjuXmPk33phZEe8hQ1h81q5NPvjrEiQAxx7LhV0kqd+vfsXLZFEo0uBU\nVJSvBSAkm3shE+8ytQD27WPrIpkFIOU5AZ1lX+oUtQDIF/HllxOnMY6iUX7ySa4v4CIREk8609aS\nNVxR4Tai6SbtSsTevdaET7cnD3DPs7aWUzIceGB67wkSgJdfBh5+2PZoL7uMl+kIgDseAJSXBSCN\nrl8AguZxZCoAGzbw55rIApA8QOUouOVIUQsAEVchWrwYuPji+NejEoDvfCdxLdwXX+TlpEnZRblk\nSthG36WpCXjkEV7PxAIAbEOcTmghkNgFZAwLwFFH2X2uAJx6KkekCO5kJlf0y7FB8gvArFnxKTcy\nFYA77+R5GzJ/pq+PB53//nfe3rePra9yGnQvZ4paAAAWgGeeYXeCn6iigLZutevLfWnq/vY3Xj7+\neLjJTukSlQDI4K9UjMrEAgAyD/9LJgA7dnCJT8Ed2+ju9hYccX3O5eoCkpBX1/3iD4OV6Jxsw0Al\n+qq7G3j6aZ7Y9d57wHnnRVPmUxkYFL0AjBnDUThBU8+jsgDcwbaPf9xG/QgjR7IrZMuW3BejcQUg\njBg0N3NvWnrUmTagUQlAXx8PJk+dave54bT+8R3XAvAf53L33fb57NgB/PM/Z3a/xYw8R9dSkgye\nIp7ynLL9nHbv5u2eHltZTDo4KgDlQ9ELgPQIx4+Pfy1qAZD0C/LjEIyxSbeiCM9MhtuIhnE5SRoH\nIZNi4EB2DUuQb7+lhZ+bJJwDvMf5Y9kTRZ34BezCC+0A/YsvAn/6U2b3W8xIOux6p8SSfH7y7Do7\nOVPqlCmZnbuy0kZ4yfkkQksmHw4erAJQLgwYAXDrnQqS2zzsLFExtWUg2G8ByPlzVY5S3CQA92rl\nfw4jNh0dNvTz/PM5cV4mRGUBtLR46whfcYWdyAfEWwCJBCBI6OV98vmIu2sgYwyHd95/v/d7/VAs\n2bo8h337sgtLlsF1aezFxQkAn4sVde3qUgEoF4peAMTkDWqQhgzhL3Mmha6D8MdbS42Bb36Tl/39\n9nq5EIDGRtvb6+jgOOyFC8NZAK4A/Pd/27KB6ZLpoHEiAWhutpk+AY4EcnvzyVxALu57pGEUy0G2\nw9QkKDTHHANcfz33zAcP5g6P+/n/4Ae8FAtAonWyYfhwKwBnnw189aveMaL2dhWAcqHoBSBVmuEx\nYzJ3b/iRBkREpqmJfyRSaEeiMaISgF27gP/8T7vt9m47O7nhnjw5OgHIhjvu8PYOU5FIAPr7uVGf\nOZOzVk6c6P1//S6gj3wk/hz+tNDivpBGSgR6IPP3v3PtY6nBMGmSNyPszJm8P6wFAMR/VuPHe59h\nZ6cKQLlQ9AIg0QqJiEIAOjqAww6zlsDu3dzTkl6rOzkpih/GH/9oyy4CtgHs6rIN98iRXBc324k4\nYQWgvp7zLqVLspnANTXcw3377fh0xO3t3vs85RSv62P7dp6U9/bbwOGH8z4JYZTzRBk6W0j6+/m7\nN3YsdwA2b7YN8+bNnB7DtQCyzU3lRl3JtvvZ/eIXbI1cf31251cGDnlKb5Y9X/kKcNxxiV8fMyZc\nMRKAG5AzzwR++Uve3rqVGyoRBHFLRGUB+GsYyI9vxw6+l+HDWQBWrcr+GmEFIFOGD09cm8EdzHQF\nwBh+jztA7Gf8eHYRLV9u6xnI2Iicxy1xmE2CtGKhv5+/y2PG8POsquLtsWP5OzJmTHQuoETba9fy\nfIPZs7M7tzKwKHoL4KKLOAFcIurr4wdtM6GpiX9448ZZS2LzZq9PVGbD5koApAe7e7d1AaVT3SkZ\nhRCARBaAlBEEuNfa389jO2+9xctUDZm/EJC4QYIEIFm922KntdV+/gCPw8j/09nJA+nd3dYqSFWm\nMxFBAuBWAlPKh6IXgFRMmcIzTbP1Ay9fDnzsY94G1xWALVuA3/+e16MWAHF1SMPZ2mob7kwnbvlp\nbY3/oeeSZALgWgBE/Fl96Uuc0CwdhgzxWkPihhMBcMdK9uzJT8qOXLB+Pbu65HOTbLg7dvD/XF3N\n/1sY9w9gBffoo3kZlGhRKQ8G/Mc+dSqnGE43aZmftWs5fNBtcF0BmDjRRsRENUFGLBa/D1sEYPjw\n8D2xZcsy8+GHRbJMuoweDZxzDvBP/xR//F13pX/uv/zF5hMyxoqwPD93JvcnP+mddTzQeOaZeAEQ\nAZUKeWHcP4Bt7J97Ln6fUl4M+I9dfJWbN2fXOG/YwBEW8iP7xS84DDTIfRKVBbBpEy8lfFUEoKXF\nOwgchg0bMg/9DIPfAtiyhRuVe+/l5xvEhAkcGZQKd2JeT0+8ALjRMh98ALz7bmb3XkysX+8VAPeZ\n1tVxUITk68kWcR25nR4VgPJkwH/srhshm0LqO3fyQOMBB/C2JC0LcsFUVGQuAJ/9rHeQesMG7s2O\nHRssAK4POAwSTpgv/ALw9tvcE0/WsGzbxgKRCplEBvAzk89A/ONbtthMo7mgu5stzP5+rtecSzZt\n8gqAO6Yxbhx/X6OyAFLtU0qfAf+xT53Kg4nHH8+VwuQHki7d3ey+qK9nMZGcQ0ECkI0F8OijnFpa\nxii2bGGxmTDBCoB/DCAK372EE+YLvwB0dqY3mSydVAavvWbXpUAQYMdSmpqAf/93dgXmgh07WLSf\neooLpeSyPnFvr/38hw8HFiywr4kAhJkDAMQPHhujAlCuZP2xE9E5RPQOEfUR0VFJjltARKuJaC0R\nXZvouDAceihPIHr+eWD//YEzzkj/vVIwnYhzsEivOZEAZONmOvNMW1WstZUbxupqXn/3Xa6YVVPj\ndQGFobOTo2sKGQXU2Zk6ncTJJ3t794mYPt2ut7RYAZCordZWdpm5KSbuvjut204LSZZ2//28zEWk\nUX09cOutvC7Pzf/5TZzIFkLYQWB3bsnPfsYpIFQAypMwH/tbAM4C8HyiA4hoMIDbACwAMAfAeUSU\nk4wt8+cDt93GjcGbb6b/Pn/BdIkhD4oqCjMGILl+pLGSPEYrVvD+CRPiXUDyI8+0NKS4f/JZzckv\nAF1dqQUgm3j9X/+aG+AhQ6wAtLXFC8CFF2Z+7kSIZSEVzRLNd8gWIg4MkFnQ8tz87rG5czkaqqUl\nnAXgTq68+mqu+TBtWvjQY2XgkbUAGGNWG2PWpDhsHoB1xpj3jTG9AO4FcGa210zG3Lk2uiaThi9R\nucmgRjdTAXBFRCwHaawqKzlSR/ZPmBDvApJeWaZWR779/0B2ApBJI7Z4MX+ut9zCUUH19SwA+/bx\nZzhihFcADj00s/tPxtNPe7clj04UiDtp+nRb8+Ltt3npLwo/YgQ31lu3hhMACUJwef55YN267M+p\nDExybfhNBrDR2d4U2xc506bZiKBMzNlMBEAGgS+91M4aToYbjy6NhriAHnkEuOkme536+ngX0Je+\nxMt9+3jwON25DsUuAJJ3PhML4N/+jTOJyrnHjmUBaGtjdx2RFYArr0wceZQNw4fbbJxAtAIgHYra\nWv4ffvMb4Otf531XXcVprvfts0IxbBhbQGFcQI89Fl/4aMyY4Iy7SmmT9GtERMsATAh46TvGmIfT\nOH9Gw2ULFy78x3pDQwMaGhrSfi8R54xZuzaTKyYWgKBet1gAv/41DxZLg5QIt0Hcu5d/xOICEqQB\nGDuWe2Dt7Xb84fbbgT/8gY/Zf3/gd78DLrgg9f+0a1dhBMAVvGQCIH7/TF1A8j91dfE51q/3Ps+R\nI7khmzjRpo0Ii8xaXrAAuOEG4NvfjtYFJBFg0mnxlz79/Oe928OG8ZhEmAR4xxyT/XuVwtPY2IjG\nxsZIzpVUAIwxp4Q8/2YATi0oTAVbAYG4ApAN0mjv2cM/kHQsgUxdQJKVMp1c/W6DuHs35zW64w7g\nxhvt/q98hZcjR3JjJvnzBTcBnYwjpGLlSu/s23zgtwAaG4PLeALWfZGpAMhz6e7mhr69PV5QTzkF\nmDMnugIxLS0sZMOGcQK/FSuiTT4n57o2zfCIoUNzG/KqFD/+zvEiSVucBVG5gBJ53V8BMJuIphPR\nUABfALAkomvGIQ3uzJlcwzcdJArIT5AFMGeOLZ7trxoWhNsgrlwJvPIKr++3nz2PUF3NjY2/QXMj\nj6Q+QSoeeww499z0jo0KqTQlPdMnngj2NQP2eWcqAGJRtLWxBRAkAPfdx5ldo3LTtLV5w1mTpbzI\nhvZ2tu78Pf1EBHVWFCVbwoSBnkVEGwEcB+BRIloa2z+JiB4FAGPMPgBXAngCwEoA9xljQuS4TI4k\nCZs3j2cGp4M/CggAzjrLlslzOfpobsgluVkqV4DbUEyfbu+prg446CBel7QFtbV8vs5ObwhqNpPP\nduzw1uDNB1JpqrPTikCiRlji0DNNZibHNzVxo9/by8/MP2u6piY6N43rkgNYAKKyAC67DPjpTzPL\n+zR0KHcg/t//i+YelPIm66EkY8yDAB4M2L8FwGnO9lIAS7O9TiaIe0YmzKRDkAvogQeCjx01inuc\ns2ZxJEZzszeO3RhvBJJYJNXVPEj9zjvAN77ByefEDTJmDOf9nzWLBy+rq72uq0xrEEh5yUIM6NXU\neBv9VIXoMw1TFQHYs4ef34gRHPkVJAAtLfGfRzYECUBUFoCElSZLd+5HyjnmM9GfUrqU1PQPaTjD\nCkAixE8flIb4scfixxykoXjkEZus7uqr2ZUhjdmIEdwA1NRwg+VvzDKdfNbWxo1e2Gyi2SDFeaSH\nvHhxtOeXyBeJg6+q4nQSQc9MomXC4hcAf36ebHEHcTOZsDdsGFs3mdZsVpQgSkoAFi/mtAGTJnGm\nw1RpgTdt4mPSFYBRozgcc+1aboxcV4DUEXbp7ORByZNOsi6Z8eO9x8gPWSwCf3x3RUXqnrTLK69E\nGwOfCSIA7e3spojaCpEwT6khUF3Nzz0o5YTfGsmWXLmA3HNkkrJDxk3UAlCioKQEoL6eferHHMNC\nkKqk3fLlwGmnpT+pRnqae/eyH9/9EUtVK399X/mhTpvGx/gHPv3i4+8NDhkCnH663U6Vh+bxx735\nY/JJba0VgHR6tZlOcPvEJ4Af/5jXhw1LbAEALACpyommg8wzEPr7bVqPMLjWiZvqIhXyfVELQImC\nkhIAYf/9gcsvTz0noLU1s3BJ90dXV+f9EctAret6cusKT50afC1/4+VvOCsrgffes9upGs0XXgA+\n/vHkx+QK1wJIxwWVKEoo1TUAfm7JBOCEE3jy2PvvZ34NF///Ihlnw1oBrrBkIwBqAShRUJICQMS5\nYB58EPja1xIf5w8hTOe8Qn29txEQv7DfAhDROOGE+MiN+fPj78/v7vHPCUjlDtq9O/9zAAR3DCAd\nAcgmnYE0gNXVyQXg058GXn3VZnfNlu5ub2MrAhzWumhv5/QfgDeFRSrEglQLQImCkhQAgMcBAC7w\nkoiWlswLrxjD75s82SsAMt6QyAU0ZEj8xKinn+a5BcL69d6UA4C3x19V5T1/EE1N6WXYzAWuBZDK\nBbRrF89szhRp+MQCCIoCArx1iMOkb+7p8QrVzTezKKczDyQZbW12jCSTkF35/NUCUKKgZAXAbQAS\nkakFIIwcGV+sI8gCSCchmsvMmfE5bNwiN8OGpbYAmpqsmyTfZOICqqtLr16AH78FsGdP8GfohlYO\nGsRRWtnQ2+sdt6mp4aydfgvAmMxmH7e32/8/EwtAxrXynepDKU1KVgCGDmUXS7LeVbYCAHBD5G/s\ngcQWQLa4Vob/mkJPD/CDH3BjmO86AC6ZjgFkgwiAWACy7mfQII6+Elavzu56fgsA4IbbPwbQ3g78\n8z/blM6pEKE2JrM0zDNn8nsKZeUppUXJCgAAfO97yX3Ab7yRvb/c3xv3C0B3N18/3RnJiaip8V6z\nu5vHE9xkZx9+CPzoR8CSJXx8PusAuNTWsghFUdQmEa4FID3oRCLupgPJppAPEG8BAMGT88QlJOk+\nUrFrV34rtilKECUtAIl6zACH861YkVn1MJehQ+N7+4DdJ4VesnFzuDzzjE09PXQoi87f/uZN5yuF\nUdatK1zvH/CGgebaAqirS24BAF73W7aFfIIsABGA9nY7ocut+5wOKgBKMVC2AtDTww1qtqXw/BaA\nXwB27uSC8Lfdlt35hfHjrRtr2DDrenDz6Ejemw0bChsdIhZALgVArJtJk6zv/IADUr8vagugt5fF\n/Yc/5H2ZDgqrACjFQNkLQFTnlkZY9u3axY13FO4YSYEwerT1ZUuv313/8MPCZousreX7W7s2d5aI\niO7gwZy07/vfT+9aubAAAGDRIvbJZyoAhYzWUhShrAUgTGMp7hhBolFcAYiqhyd5hOrquMEBvIOa\ne/eyJdPUVFgLQKKP7rkndxbA/Pl2gl9dne2BpyJbAUhnDKC7O3MBSHeuhKLkkrIVgG3bwvXO/efe\ns4dDT2Xfli3xeX+y5ZBDuJf5wQfs5gF4bEBoauJrr1kTzqoJi3vtXLqAZs3K/H2Z5FMSvv99nkfi\ntwD8Cfqam1kATjqJB+jTIZduMkVJl7IQgL4+bjjcCUGHHRZuNmeQBTBxoo0Geu45TvscJdLA1tba\nMYfNm7lM4cSJ3Cj9z/9Ee81MuflmXhZb49bWltnxfX02508iC0BCPlta+POvr0/f0shlpJSipEtZ\nCICUb8zWDZDs3FIFa+tWLvJyxRVczHv79swm+KTDX/7CkT5vvGGF5t57eZnOxLd8IBOUiqVxezhW\nuTpTAXDDd4PGAHp7+bMnshZAfX36g80qAEoxUNICUFnJDaUM0EZZym/oUOC//5uF4J57eFr/pZfy\na5mmmU6XujpOdFdf77UAgOJJDVBby8tisQDkPjIVgC1bgAMP5PVEg8A9Pfy5t7SwAEyYkH4nI92M\nqYqSS0paAIYN4x+lTM6Jspi327hfcAG7DI4+Grj4YmD2bLYOcjUgO2QIX2/FCmDjRt6XbThr1IgF\nUCwCILmWWlsze9+WLbZsp7/j4ArA5MmcdDAbC6BYnpFSvhRJs5EbiLg+8LPP8nbUFkDQ9ty53LsL\nG2WUDvPmcUM1eTKHRBYDxeYCqq9noUzXAliyhBvzXbtssrbt273HuAJw4YWcrjvTMQC1AJRiIOua\nwAOFKVOA11/n9c5O7plHUTDc37hLEZaqKj5/RUV+euXNzd6UB5/9bO6vmYxiswAAnrCVrgCceSZw\n7rkcMDB+PM838EcdVVRw7iWAx16yGQNQAVCKgZK2AACeMSox852dHNonedjD4KYfmDsXuOUWXpcM\nlfmakOVPaS2DnoVi9OjC1SRORH09h/26sfp9fcADDwQff++9wH33sQtv3z7gnHO8r1c43abx43nW\nd0sLWwxiAUhB+r6++PN3d/N5VQCUQlPyAjB5sjXhOzrC53EX3JTLxxxj10eMyK8AbNzI2SQPOQR4\n/vn8XDMZgwaxyBYqJXUQY8YAJ55oXYEAu85k0D4RiT5Dd1B43Djg3Xe5DrMIBmCjtILqUu/dy/dU\nqKR9iiKUvAAccohd7+y0ydlmz+aykdniNnANDXa9qopFJp8zckeO5MbkxBPzd81kLFrk7SUXA8cc\nwxXChO7u1GNCiT5DvwUA8PdMwkMBG3ocdI29ezUNhFIcZC0ARHQOEb1DRH1EdFSS494nojeJ6DUi\nejnb62XLwQfbdTc//+mnA//3/2Z/XvcHLKGPAMf+r13rzdWTa4qtsS1Gpk+3EVMA98w7O5NXC0tk\nAbjP2w17dWcIS9RRUORZIYv2KIpLGAvgLQBnAUjleDAAGowxRxpj5oW4Xla4P+KODjsYGGSaZ4Lb\nCLg++NmzOTonypDTIG66KbfnLzUmTuTJekKQi2bNGu970rEA3KysrgUgAhBkAWgiOKVYyFoAjDGr\njTFrUh8JACgKb2dnpzXN3XKOYfEPeH70o9GdOxHz5/Ny3brcX6sU8AuA5GySBnrfPjvxS0iUVyko\nuqu/32sBHBWziYM6Aps389iUohSafIwBGABPEdErRJRi2C23dHbanlmUudgPPdS7/bnP5a9mq7+G\nsBLMuHHe3E/S85cGeseO+Pckiul3XUlCX5+1AFy3UpAFsHGjzfCqKIUkqfeYiJYBCAqa/I4xJt2A\nwxOMMVuJaByAZUS02hizPOjAhQsX/mO9oaEBDe7oagR0dPAP/667ONY7Ck49Nb6n2NAQLtFcOkgv\nVCNJ0qOmhuP1Bb8F4DbaJ53EyfwSufGmTInf19/Pn8mgQV7hSCQAxx+f2f0ritDY2IjGxsZIzpVU\nAIwxp4S9gDFma2y5k4geBDAPQEoByAWdnfw3dWp0aZOzrTQVlsMPB95+uzDXHohUVnIj3dVlc0QB\nLNQ33QRcdx27ib72Nc6uSpRYAL7yFU4v8pvf2H0yt6SigucFjBnDjXzQOfbsyZ+FqJQe/s7xIikS\nkgVRuYAC+6FENIKIRsbWqwCcCh48ziu7dnFoYmcn//CjTJwWNNEnHxB5Q1yV5BB5rQARgOXLgTvv\n5F77iBHc+AuTJiU+3623WlfQ9u22MM2QITzprK6Oz+dPIwHwPeggsFIMhAkDPYuINgI4DsCjRLQ0\ntn8SET0aO2wCgOVE9DqAlwA8Yox5MuxNZ0pdHU+WEgsgyhj9QgmAkjmuAIgLaNUqXkr6DqG9Hbjk\nksTnqqqyrqDx461FWVHBjf7YsdzRuPTS+KJEe/fyvShKock6gtwY8yCABwP2bwFwWmx9A4Ajsr67\nCBk+nM1xdy5AFKgADBxGj463AN5/n5eXXead4ZttmgbXApBxmvXr2e0jBYJ0IphSLJT8TGBh+HDr\nAlILoDwJcgGJi+b11+Pz/mdDRQVHFI0ZYweWL77YO0tbBUApFspGACZP5rS9agGUL0EuoG3b7OtR\nzKgeMoTHnNyZvm7YaH8/z0UZNSr8tRQlLGUjAMcey5EbUY8BFCoKSMmcIAugqQk4/3xuuKOyAHbu\n5B6+WADuBLSHHmIR0PQdSjFQNgJQXc0/yPZ2tQDKlSABAHgwt7c3OgtABCCIs88Ofw1FiYqyEQDA\nTtCJqvf1858DN9wQzbmU3NPfz/H+gLdk59SpvIzSAnDHAIRkiecUpRCUnSFaWxvd7NmrrormPEp+\nGDHC+v67uriR3rrVpvKIQgASWQDjx6u7UCk+ysoCAGz+dqX8+OEPuYduDAuAiMHcubyMwjJ0xwDc\ngeAhQ6zbqZiqpSnlTdkJgJu6WSkvhg7lBloa/8suA/70J9sgR2UBdHSwAPzHfwCnxJKpdHXxX11d\n+vWJFSXXlJ0LKF+lGpXiZNQoDsPs6uIqYWefbX3zQWmeM+XFF3k5Zgy7nCTts4hOPivFKUoqys4C\nUAEob1wBkMZYxoSi9NHLGMBPfgI88oi1APT7pxQTZScAUWUBVQYmo0ZxKGhQbzzKkF6Z6DVxInDa\naWxltLWpBaAUF2UnANoDK2/q6oDdu4MnBEZlAXzsY/HupMpKFh4VAKWYKCsBmD0b+PSnC30XSiGZ\nMIHz/uzdG5+Oob8/mmtUV8fvq6zka2oHRCkmymoQ2F/0Wyk/6uttzn+/AEQRIVZRYbN+uqgFoBQj\nZSUAiuLW4vXn5Hfj9rOluzt4omFlJfDww5o6RCkuysoFpChf/arNx+Pv8UeRonnQoMQC0NgIfOIT\n4a+hKFGhAqCUFYMHAx/5iF0X/vIXbznIqKms5DTRX/xi7q6hKJlCpkgyVBGRKZZ7UUqbdeu4nrK/\nVGMuOe444KWXgJ6eaGYcK4pARDDGZJXhTC0ApeyYNSu/jT9g8wBp468UEyoAipIHNP+PUoyoAChK\nHmhvL/QdKEo8GgaqKHngkEOiCTNVlCjJehCYiH4K4LMAegCsB3CxMaY54LgFAH4OYDCAXxtjbkxw\nPh0EVkqWnh6eAxBlOVJFAQo3CPwkgEOMMXMBrAFwXcCNDQZwG4AFAOYAOI+IDg5xzbKgsbGx0LdQ\nNJTKsxg6NHzjXyrPIgr0WURD1gJgjFlmjJHsKS8BmBJw2DwA64wx7xtjegHcC+DMbK9ZLuiX26LP\nwqLPwqLPIhqiGgT+NwCPBeyfDGCjs70ptk9RFEUpMEkHgYloGYAJAS99xxjzcOyY7wLoMcb8IeA4\ndeoriqIUKaFmAhPRRQAuBfBJY0xXwOvHAVhojFkQ274OQH/QQDARqVgoiqJkQbaDwFmHgcaie74F\n4KSgxj/GKwBmE9F0AFsAfAHAeUEHZvsPKIqiKNkRZgzgFwCqASwjoteI6HYAIKJJRPQoABhj9gG4\nEnp3DkwAAAOMSURBVMATAFYCuM8YsyrkPSuKoigRUDTJ4BRFUZT8UvBUEES0gIhWE9FaIrq20PeT\nT4hoKhE9S0TvENHbRPS12P5aIlpGRGuI6EkiiiBT/cCAiAbHLEoJMijLZ0FEo4nofiJaRUQriejY\nMn4W18V+I28R0R+IaFi5PAsi+g0RbSeit5x9Cf/32LNaG2tTT011/oIKgE4UQy+AbxhjDgFwHIAr\nYv//twEsM8YcAODp2Ha5cBXYXSimabk+i1sAPGaMORjA4QBWowyfRWz88FIARxljDgNnFDgX5fMs\nfgtuH10C/3cimgMeZ50Te8/tRJS0jS+0BVDWE8WMMduMMa/H1tsArALPkzgDwF2xw+4C8LnC3GF+\nIaIpAD4D4NcAJCig7J4FEdUAONEY8xuAx9JiaVbK7lkAaAF3lEYQUQWAEeCAkrJ4FsaY5QCafLsT\n/e9nArjHGNNrjHkfwDpwG5uQQguAThSLEevpHAmeVV1vjNkee2k7gPoC3Va+uRkcWdbv7CvHZzED\nwE4i+i0RvUpEvyKiKpThszDG7AHwMwAfghv+vcaYZSjDZ+GQ6H+fBG5DhZTtaaEFQEegARBRNYA/\nA7jKGNPqvhbLkFfyz4mIPgtghzHmNdjev4dyeRbg8OyjANxujDkKQDt8Lo5yeRZEtD+ArwOYDm7g\nqonoAveYcnkWQaTxvyd9LoUWgM0ApjrbU+FVsJKHiIaAG//fGWMeiu3eTkQTYq9PBLCjUPeXR44H\ncAYRvQfgHgDzieh3KM9nsQnAJmPMitj2/WBB2FaGz+IYAH8zxuyOhZU/AOCjKM9nIST6Tfjb0ymx\nfQkptAD8Y6IYEQ0FD2AsKfA95Q0iIgCLAaw0xvzceWkJgAtj6xcCeMj/3lLDGPMdY8xUY8wM8CDf\nM8aYL6I8n8U2ABuJ6IDYrpMBvAPgYZTZswAPfh9HRMNjv5eTwUEC5fgshES/iSUAziWioUQ0A8Bs\nAC8nPZMxpqB/AD4N4F3wgMV1hb6fPP/vHwP7u18H8FrsbwGAWgBPgdNsPwlgdKHvNc/P5SQAS2Lr\nZfksAMwFsALAG+Beb00ZP4trwAL4FnjQc0i5PAuwNbwFXHdlI4CLk/3vAL4Ta0tXA/hUqvPrRDBF\nUZQypdAuIEVRFKVAqAAoiqKUKSoAiqIoZYoKgKIoSpmiAqAoilKmqAAoiqKUKSoAiqIoZYoKgKIo\nSpny/wGssKBxJclwOgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10c46ae10>"
       ]
      }
     ],
     "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: 1.04 ms per loop\n",
        "10000 loops, best of 3: 162 \u00b5s per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1 loops, best of 3: 428 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": {}
  }
 ]
}