examples/stats/custom-model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Writing a model\n",
    "\n",
    "This example shows you how to write models (or model wrappers) that can be used with Pints.\n",
    "\n",
    "Pints is intended to work with a wide range of models, and assumes as little as possible about the model's form.\n",
    "Specifically, a \"model\" in Pints is anything that implements the [ForwardModel interface](http://pints.readthedocs.io/en/latest/core_classes_and_methods.html#forward-model).\n",
    "In other words, it's anything that can take a parameter vector $(\\boldsymbol{\\theta})$ and a sequence of times $(\\mathbf{t})$ as an input, and then return a vector of simulated values $(\\mathbf{y})$:\n",
    "\n",
    "$$f(\\boldsymbol{\\theta}, \\mathbf{t}) \\rightarrow \\mathbf{y}$$\n",
    "\n",
    "This might not match well with other ideas of a \"model\".\n",
    "For example, if you're modelling something using ODEs, Pints will regard the combination of the ODEs _and an ODE solver_ as a \"forward model\".\n",
    "\n",
    "In the example below, we define a system of ODEs (modelling a simple chemical reaction) and use [SciPy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html) to solve it.\n",
    "We then wrap everything in a `pints.ForwardModel` class, and use a Pints optimisation to find the best matching parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving a system of ODEs\n",
    "\n",
    "First, we define an ODE. In this example we'll use a model of a reversible chemical reaction:\n",
    "\n",
    "$$\\dot{y}(t) = k_1 (1 - y) - k_2 y,$$\n",
    "\n",
    "where $k_1$ represents a forward reaction rate, $k_2$ is the backward reaction rate, and $y$ represents the concentration of a chemical solute.\n",
    "\n",
    "Next, we write a Python function `r` that represents the right-hand side of this ODE, choose a set of parameters and initial conditions, and solve using SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XXWd//HXJzd72qZL2kLXFGiR0EKBtCwiMMiqUlSQAZcZFK064saMM8zGD9HHbxxxGcbh8dPqMCIuKKLSEbDIIiAWSEo3uqd7uidt0nTJcu/9/P64t+k1pMltm5Nzl/fz8biPnO3e+zlJe9/3fL/nnK+5OyIiIgAFYRcgIiKZQ6EgIiLdFAoiItJNoSAiIt0UCiIi0k2hICIi3RQKIiLSTaEgIiLdFAoiItKtMOwCjldVVZVXV1eHXYaISFZZtGhRk7uP7m+7rAuF6upq6uvrwy5DRCSrmNnmdLZT85GIiHRTKIiISDeFgoiIdFMoiIhIN4WCiIh0UyiIiEg3hYKIiHTLuusURCT7uDtxh1jcibsTizvRuBOPOzE/+jMaO7o+7oltYnEnHoeYO7F4nLhDPO44EHcHh7iDk3gPd8ePzMfp3s6PrEuZjyeHIz463+P5Tvd2ztF1PZ+f2MfkTzxluvflqb+X/rb1lBd/51ljOXfi8AH7u/RGoSCSI+Jxpz0ao70rTntXLPmIJ5fF6IjG6ejqsT6amO6KxYnGnM5YnK5YnK6o0xWP0xVzuqJxovE4ncnprlicrvjR6Wjc6UyZjsYSH9zReDzlw1xjwZ8sMxgzrFShIJLL3J1DnTFaDnfRcqiT1sNd7D/cRVt7lIMdUQ52xhI/O6Ic6EhOdybXdcQ4kJw/1BGjMxY/4ToiBUZRxCiKFCQfieniSAGFPZYXFxZQ0WO7I9OFEaOwoIACS0wXmBEpgIgZkYICIgVQUGDJ+aOPAjMKC+yY6yIFR9cXGBSYYQZGYt7syM/EtHF0mwIz6DGf+An0fD7J17Ue25N83ZT5gkQBJF++e5uj00eWW/c0PZYfmbfudZYyffT1BpNCQWQAtXfFaDrQQdOBTpoPdHRPNx3ooOXQ0Q/+luSHf8uhLqL9fIuOFBgVxRGGlBRSkXwMKSlk9NASKooT8+UlEcqKIpQWRSgpLKC0KEJpUQGlhcllRcllhcnlyW1LiwooKYwQKRj8Dx/JTAoFkTTE407zwU52tB5me0s7O1oPs6O1ne0th9m1vz3xwd/WQVtHtNfnDykpZHh5EcPLi6gsK+LUyjKGlSXmh5cllg0vL2JYcnpoSREVJREqSgopKSwI5Ruj5CeFgkjSgY4om5sPsrn5UPJxkE3NB9ne0s7O1va3NM8UFxYwrrKUscNKqRk3jNFDSqgaUkzVkJLEY+jR+dKiSEh7JXJ8FAqSV9ydXfs7WLurrfuxYc9BNjUfoulAx59tWzWkhMmjypk5cTinzihlXGUZp1aWMm544ufIimJ9g5ecE2gomNl1wANABPiBu3+tl21uAe4lcebVUnf/YJA1Sf5o74qxemcbyxtbWLmjjXW72lizq4229qNNPFVDSjh9dAVXnTWGSaPKqR5VweRR5UweVcGQEn1nkvwT2L96M4sADwJXA41AnZnNd/eVKdtMBf4ReLu77zOzMUHVI7ktGouzemcbSxtbWN7YyvJtrazZ2dbdiVtZVsSZY4dy48xxTBs7tPsxsqI45MpFMkuQX4VmAw3uvgHAzB4FbgRWpmzzCeBBd98H4O67A6xHckh7V4zFW1qo27SXuk17eWPzPg52xgAYXl7EjPGVzL3sNM6ZUMmMCcMZV1mqph6RNAQZCuOBrSnzjcCFPbaZBmBmr5BoYrrX3X8XYE2SpaKxOEsbW3lp7R7+2NDEssYWumKOGZw5dijvP38CtdUjOH/SCCaMKFMAiJygIEOht/+VPU/ILgSmAlcAE4CXzWy6u7f82QuZzQXmAkyaNGngK5WMtGt/Oy+s3s1L6/bwx3VN7G+PYgbnTBjOxy6dwuzqkdROHklleVHYpYrkjCBDoRGYmDI/AdjeyzavunsXsNHM1pAIibrUjdx9HjAPoLa2VtfL57BNTQf53YqdLFixk8VbEt8Nxg4r4dqzT+HyM0dz6RlVDC9XP4BIUIIMhTpgqplNAbYBtwI9zyz6DXAb8EMzqyLRnLQhwJokA21uPshvFm/n6Td3sHpnGwAzxlfyd9dM46qasZw5dqiag0QGSWCh4O5RM7sTWECiv+Ahd19hZvcB9e4+P7nuGjNbCcSAL7l7c1A1SeZoPdTFb5dv59dvbKN+8z7MoHbyCP71PTVce/ZYJowoD7tEkbxk3vNerhmutrbW6+vrwy5DToC78+qGvfz4tc38fsUuOmNxpo4ZwvvOH897Z45n3PCysEsUyVlmtsjda/vbTlfnSOD2t3fxq0WN/Pi1LTTsPkBlWREfvHASN50/genjh6lpSCSDKBQkMI37DvGDlzfy87qtHO6Kce6ESu6/+RxuOHec7gUkkqEUCjLgVu/cz/de3MD8pdsxYM7Mcdx+STXnTAh2cBAROXkKBRkwK7a38q1n1vLc6t2UF0e4/ZJq7rh0ivoKRLKIQkFO2vo9B/jW79fy5LIdVJYVcdfV0/iriyfregKRLKRQkBO2a38731iwhsffaKS0KMJnrzyDj7/jNCrLdIWxSLZSKMhx64jG+O8/buS/nm8gGnM++vYpfPqK06kaUhJ2aSJykhQKclyeXbmLrzy5ks3Nh7i6Ziz/8u6zmDyqIuyyRGSAKBQkLbv3t3PPEyv43YqdnDFmCI/cMZt3TB0ddlkiMsAUCtInd+ex+ka++uRK2qNx/v66M/nEO06jKFIQdmkiEgCFghzTjtbDfOmxZfyxoYnZU0bytffP4LTRQ8IuS0QCpFCQXj29fAd3/2o5XbE4X33vdD44exIFBbodhUiuUyjInznUGeXL81fy8/qtnDOhkgduPY8pVepIFskXCgXptnZXG5/68SI2Nh3kb644nS9ePU19ByJ5RqEgADy1fAd/99hSyosL+cnHL+SS06vCLklEQqBQyHOxuHP/gjV898X1nDdpOP/vQxdwSmVp2GWJSEgUCnnsQEeUz/zkDV5cu4cPXTiJe26ooaRQt7QWyWcKhTy1s7Wdj/6wjrW72vi/75vBBy+cFHZJIpIBFAp5aPXO/Xz0f+rYf7iLh26fxeXTdGWyiCQoFPLMwvXNzP1RPeUlEX7xqYs5e1xl2CWJSAZRKOSRP6zZzScfWcSkkeU8/LHZGvxGRN5CoZAnFqzYyZ0/fYOpY4byyB2zGaXbXItILxQKeWD+0u188edLmDG+koc/OpvKcg2CIyK9UyjkuCeX7eALjy6mtnokD90+iyEl+pOLyLHpEyKHvbB6N59/dDHnTxrBDz86i/Ji/blFpG+6sU2O+tP6Jj7140WcdeowHlIgiEiaFAo5aOnWFj7+cH33WUbDStWHICLpUSjkmK17D3HHw3WMrCjmJx+/kJEVxWGXJCJZRG0KOaT1UBe3/8/rdMWcR+fOZsww3dhORI5PoEcKZnadma0xswYzu7uX9beb2R4zW5J8fDzIenJZRzTG3Efq2br3MPM+cgFnjNGwmSJy/AI7UjCzCPAgcDXQCNSZ2Xx3X9lj05+7+51B1ZEP3J1/fHw5r23cywO3zuTC00aFXZKIZKkgjxRmAw3uvsHdO4FHgRsDfL+89cM/beJXi7fxxaumcePM8WGXIyJZLMhQGA9sTZlvTC7r6SYzW2ZmvzSzib29kJnNNbN6M6vfs2dPELVmrVc3NPPVJ1dxdc1YPnvlGWGXIyJZLshQsF6WeY/5/wWq3f0c4Fng4d5eyN3nuXutu9eOHq3bPB+xveUwn/nJG0weVc63bjmXgoLefuUiIukLMhQagdRv/hOA7akbuHuzu3ckZ78PXBBgPTmlIxrj0z9eREc0zryP1DJU1yKIyAAIMhTqgKlmNsXMioFbgfmpG5jZqSmzc4BVAdaTU/796TUsbWzlm7ecqzONRGTABHb2kbtHzexOYAEQAR5y9xVmdh9Q7+7zgc+Z2RwgCuwFbg+qnlzy/OpdPPTKRm6/pJprzz4l7HJEJIeYe89m/sxWW1vr9fX1YZcRml3727n+gZcZO6yUX//NJZQWRcIuSUSygJktcvfa/rbTbS6ySCzufOHRJRzujPGd285TIIjIgNNtLrLID17ewMINzXz9pnPUjyAigdCRQpZo2N3GN3+/lmvPHssHaieEXY6I5CiFQhaIxuL87WPLqCiO8NX3zsBM1yOISDDUfJQFvv/yRpZubeE/bzuP0UNLwi5HRHKYjhQy3LpdbXz792u57uxTuOGcU/t/gojISVAoZLB43Pn7x5dRURLhK++drmYjEQmcQiGD/axuC4u3tPCv76lRs5GIDAqFQoba09bBvz+9motPG8X7ztPtsEVkcCgUMtS/PbWKw10xNRuJyKBSKGSgheub+dXibXzystN1kZqIDCqFQoaJxuLc88SbTBxZxp0aNEdEBplCIcP87PUtrNt9gH95d43ubSQig06hkEFaD3Xxrd+v5aLTRnJNzdiwyxGRPJTWFc1mNhr4BFCd+hx3/1gwZeWn7zy/jpbDXfzre2rUuSwioUj3NhdPAC+TGEc5Flw5+Wtj00EeXriJWy6YyNnjKsMuR0TyVLqhUO7u/xBoJXnua0+vojhSwN9eOy3sUkQkj6Xbp/BbM3tXoJXkscVb9rFgxS4+efnpjBlaGnY5IpLH0g2Fz5MIhnYza0s+9gdZWD75xjNrGFVRzB2XTgm7FBHJc2k1H7n70KALyVevNDTxSkMz97ynhooS3clcRMKV9qeQmc0BLkvO/sHdfxtMSfnD3bl/wRrGVZbywQsnhV2OiEh6zUdm9jUSTUgrk4/PJ5fJSXh21W6WbG3h81dN1YVqIpIR0j1SeBcw093jAGb2MLAYuDuownJdPO5885k1TKmq4KbzNeayiGSG47mieXjKtE6kP0m/X7WL1Tvb+Nw7z6AwogvLRSQzpHuk8G/AYjN7ATASfQv/GFhVOc7d+a/nG5g0spwbzhkXdjkiIt3SPfvoZ2b2B2AWiVD4B3ffGWRhueyldU0s39bK194/Q0cJIpJR+vxEMrO3JX+eD5wKNAJbgXHJZXKc3J3vPLeOUytLeb/6EkQkw/R3pHAXMBf4Zi/rHLhywCvKca9t3Ev95n3ce0MNxYU6ShCRzNJnKLj73OTk9e7enrrOzPq9H4OZXQc8AESAH7h7r6exmtnNwGPALHevT6fwbPXgCw1UDSnm1tm6LkFEMk+6X1X/lOaybmYWAR4ErgdqgNvMrKaX7YYCnwNeS7OWrLVieysvr2viY5dO0XUJIpKR+jxSMLNTgPFAmZmdR6KTGWAYUN7Pa88GGtx9Q/K1HgVuJHHxW6qvAF8H/u74Ss8+//3yRsqLI3xo9uSwSxER6VV/fQrXArcDE4BvpSxvA/6pn+eOJ9EpfUQjcGHqBsmgmejuvzWznA6Fna3tzF+6nQ9fNJnK8qKwyxER6VV/fQoPAw+b2U3u/vhxvnZvQ4d590qzAuDbJEKn7xcym0uiw5tJk7KzLf5HCzcRc+djb9edUEUkc6V7ncLjZvZu4GygNGX5fX08rRGYmDI/AdieMj8UmA78ITn05CnAfDOb07Oz2d3nAfMAamtrnSxzqDPKT17bwrU1pzBpVH+tbiIi4Un3hnjfBf4S+CyJI4APAP01jNcBU81sipkVA7cC84+sdPdWd69y92p3rwZeBd4SCLng8UWNtB7u4uPv0FGCiGS2dM8+usTd/wrY5+5fBi7mz48C3sLdo8CdwAJgFfALd19hZvclb8OdF+Jx56FXNjFz4nAumDwi7HJERPqU7r2PjlyjcMjMxgHNQL9fe939KeCpHsvuOca2V6RZS1Z5uaGJjU0HeeDWmSSbyUREMla6ofC/ZjYcuB94g0SH8fcDqyqHPLJwM1VDirl++qlhlyIi0q9+QyF5ltBz7t4CPG5mvwVK3b018Oqy3LaWwzy/ehefuvx03dJCRLJCv59UyYF1vpky36FASM+jr2/BQUNtikjWSPfr6zNmdpOpUTxtndE4P3t9K1eeOYYJI3Qaqohkh3T7FO4CKoCombWTOC3V3X1YYJVluWdW7qTpQAcfvki3tBCR7JHuxWtDgy4k1/z41c1MHFnGZdNGh12KiEja0r147bl0lknChj0HeHXDXm6bPYlIgVrcRCR79HeX1FISd0OtMrMR/PldUjW48DE8tqiRSIFxs0ZWE5Es01/z0SeBL5AIgEUcDYX9JMZKkB5icedXbzRyxbTRjBnW7zhEIiIZpb+7pD4APGBmn3X37wxSTVntpXV72LW/gy/P0VGCiGSfdDuav2NmlwDVqc9x9x8FVFfWeqx+KyMrirnybWPDLkVE5LilFQpm9ghwOrAEiCUXO6BQSLHvYCfPrtzNhy6apCuYRSQrpXudQi1Q4+5ZN5bBYHpiyTY6Y3E+cEGfN5AVEclY6X6dfZPEIDjSh8cWNTJ9/DBqxumaPhHJTukeKVQBK83sdaDjyEJ3z5txEfqzeud+Vmzfz7031IRdiojICUs3FO4NsohcMH/JdiIFxnvO1eUbIpK90j376EUzmwxMdfdnzawciARbWvZwd55Ysp1Lz6iiakhJ2OWIiJywdG9z8Qngl8D3kovGA78Jqqhs88aWfWxrOcyNM3WUICLZLd2O5s8AbydxJTPuvg4YE1RR2eaJJdspKSzgmrPVFy8i2S3dUOhw984jM2ZWSOI6hbwXjcV5ctkOrjprLENK0u2iERHJTOmGwotm9k9AmZldDTwG/G9wZWWPV9Y303ywkzlqOhKRHJBuKNwN7AGWk7hJ3lPAvwRVVDZ5Ysk2hpYWcsWZGjdBRLJfuu0dZcBD7v59ADOLJJcdCqqwbNDeFeOZFbt414xTKCnUyVgikv3SPVJ4jkQIHFEGPDvw5WSXF9fu4UBHlBt0bYKI5Ih0Q6HU3Q8cmUlO5/1o9L97cyeVZUVcdNqosEsRERkQ6YbCQTM7/8iMmV0AHA6mpOzQGY3z7KpdXF0zlqKI7ogqIrkh3T6FLwCPmdn25PypwF8GU1J2eGV9E23tUa6frmsTRCR3pHubizozextwJokhOVe7e1eglWW43y3fyZCSQi6dWhV2KSIiA+Z4rraaxdGR184zs7wdeS0ai/PMyp1c+bYxOutIRHJKuvc+egT4BnApiXCYRWLgnf6ed52ZrTGzBjO7u5f1nzKz5Wa2xMz+aGZZcd/p1zfuZd+hLjUdiUjOCWzkteS1DA8CVwONQJ2ZzXf3lSmb/dTdv5vcfg7wLeC6dN8jLE+/uZPSogIu1wVrIpJjghx5bTbQ4O4bkvdNehS4MXUDd9+fMltBFtxPKR53FqzYyRXTxlBerHsdiUhuCXLktfHA1pT5RuDCnhuZ2WeAu4Bi4Mo06wnNksYWdrd1cJ2ajkQkBwU58pr1suwtRwLu/iDwoJl9kMT9lP76LS9kNheYCzBp0qQTKGXgPLdqF5EC4y/O1J3DRST3pNV85O4vAquBocnHquSyvjQCE1PmJwDbj7EtJJqX3nuM95/n7rXuXjt6dLjt+M+t2k3t5BFUlheFWoeISBDSPfvoFuB14APALcBrZnZzP0+rA6aa2RQzKwZuBeb3eN2pKbPvBtalW3gYGvcdYvXONq46a2zYpYiIBCLd5qN/Bma5+24AMxtN4oZ4vzzWE9w9amZ3AgtIjOf8kLuvMLP7gHp3nw/caWZXAV3APnppOsokz6/eDcCVZ6npSERyU7qhUHAkEJKaSeMow92fIjH2Quqye1KmP5/m+2eE51btZkpVBaePHhJ2KSIigUg3FH5nZguAnyXn/5IeH/a57mBHlIXrm/nIxZPDLkVEJDB9hoKZnQGMdfcvmdn7SVzRbMBC4CeDUF/G+GNDE52xOO9U05GI5LD+moD+A2gDcPdfuftd7v5FEkcJ/xF0cZnk+VW7GVpayKzqkWGXIiISmP5Codrdl/Vc6O71JG6Olxficee51bu5fNpojZ0gIjmtv0+40j7WlfWxLqcs39ZK04EONR2JSM7rLxTqzOwTPRea2R3AomBKyjwvrd2DGVw2VTfAE5Hc1t/ZR18Afm1mH+JoCNSSuE/R+4IsLJO8vK6J6eMqGTWkJOxSREQC1WcouPsu4BIz+wtgenLxk+7+fOCVZYi29i7e2LKPuZedFnYpIiKBS3c4zheAFwKuJSP9aX0z0bhz2TQ1HYlI7tOpNP14ed0eKoojnD9pRNiliIgETqHQj5fWNnHx6aMoLtSvSkRynz7p+rC5+SBb9h5S05GI5A2FQh9eWrsHgHfoVFQRyRMKhT68uLaJiSPLqB5VHnYpIiKDQqFwDF2xOAvXN3HZ1NGY9TayqIhI7lEoHMPiLS0c7Iyp6UhE8opC4Rj+tL6JAoOLTx8VdikiIoNGoXAMf1rfzNnjKqksKwq7FBGRQaNQ6EV7V4wlW1p0lCAieUeh0ItFm/fRGYtz8WkKBRHJLwqFXixc30ykwJg1RaOsiUh+USj0YuGGZmaMr2RISVr3CxQRyRkKhR4OdkRZulX9CSKSnxQKPdRt2ks07lyiUBCRPKRQ6GHhhmaKIkbtZPUniEj+USj08Or6ZmZOHE5ZcSTsUkREBp1CIcX+9i6Wb2vVqagikrcUCinqNu4l7nCR+hNEJE8pFFLUbdpHUcQ09KaI5K1AQ8HMrjOzNWbWYGZ397L+LjNbaWbLzOw5M5scZD39qd+0lxnjKyktUn+CiOSnwELBzCLAg8D1QA1wm5nV9NhsMVDr7ucAvwS+HlQ9/WnvirGssZVZ1TrrSETyV5BHCrOBBnff4O6dwKPAjakbuPsL7n4oOfsqMCHAevq0fFsrnbE4tQoFEcljQYbCeGBrynxjctmx3AE83dsKM5trZvVmVr9nz54BLPGo1zfuBeCCyepPEJH8FWQo9DaGpfe6odmHgVrg/t7Wu/s8d69199rRo4MZCa1+017OGDOEkRXFgby+iEg2CDIUGoGJKfMTgO09NzKzq4B/Bua4e0eA9RxTPO7Ub97HrGodJYhIfgsyFOqAqWY2xcyKgVuB+akbmNl5wPdIBMLuAGvp09rdbbS1R9XJLCJ5L7BQcPcocCewAFgF/MLdV5jZfWY2J7nZ/cAQ4DEzW2Jm84/xcoGq27QPQKEgInkv0AED3P0p4Kkey+5Jmb4qyPdPV/2mvYwdVsKEEWVhlyIiEipd0QzUb9pHbfVIzHrrGxcRyR95HwrbWg6zreUws3QqqoiIQqF+U+L6BF20JiKiUGDxlhbKiiK87ZShYZciIhK6vA+FJVtbOGdCJYWRvP9ViIjkdyh0RGOs3L6fmZOGh12KiEhGyOtQWLF9P52xOOdNVCeziAjkeSgs2dICwHk6UhARAfI8FBZvbWFcZSljh5WGXYqISEbI61BYsnWf+hNERFLkbSg0Hehg697D6k8QEUmRt6FwpD9BRwoiIkflbSgs3rqPwgJj+rjKsEsREckYeRsKS7a28LZTh1JWHAm7FBGRjJGXoRCLO0u3tqo/QUSkh7wMhYbdBzjQEWXmRPUniIikystQWLI1MdKaOplFRP5cXobCssZWhpYUMmVURdiliIhklLwMhTe3tTJ9fCUFBRppTUQkVd6FQmc0zqqdbcyYoFNRRUR6yrtQWLurjc5onOnjFQoiIj3lXSi8ua0VgBkKBRGRt8i7UFi+LdHJPHlkediliIhknLwLBXUyi4gcW16FQmc0zqod6mQWETmWvAqFtbva6Iypk1lE5FjyKhTUySwi0re8CoXl21oZWqpOZhGRYwk0FMzsOjNbY2YNZnZ3L+svM7M3zCxqZjcHWQskQmH6OHUyi4gcS2ChYGYR4EHgeqAGuM3ManpstgW4HfhpUHUc0RmNs1qdzCIifSoM8LVnAw3uvgHAzB4FbgRWHtnA3Tcl18UDrANQJ7OISDqCbD4aD2xNmW9MLguFOplFRPoXZCj01nDvJ/RCZnPNrN7M6vfs2XNCxYysKObqmrHqZBYR6UOQzUeNwMSU+QnA9hN5IXefB8wDqK2tPaFguebsU7jm7FNO5KkiInkjyCOFOmCqmU0xs2LgVmB+gO8nIiInKbBQcPcocCewAFgF/MLdV5jZfWY2B8DMZplZI/AB4HtmtiKoekREpH9BNh/h7k8BT/VYdk/KdB2JZiUREckAeXVFs4iI9E2hICIi3RQKIiLSTaEgIiLdFAoiItLN3E/oWrDQmNkeYPMJPr0KaBrAcrKB9jk/aJ/zw8ns82R3H93fRlkXCifDzOrdvTbsOgaT9jk/aJ/zw2Dss5qPRESkm0JBRES65VsozAu7gBBon/OD9jk/BL7PedWnICIifcu3IwUREelDToaCmV1nZmvMrMHM7u5lfYmZ/Ty5/jUzqx78KgdWGvt8l5mtNLNlZvacmU0Oo86B1N8+p2x3s5m5mWX9mSrp7LOZ3ZL8W68ws8DHPw9aGv+2J5nZC2a2OPnv+11h1DlQzOwhM9ttZm8eY72Z2X8mfx/LzOz8AS3A3XPqAUSA9cBpQDGwFKjpsc3fAN9NTt8K/Dzsugdhn/8CKE9Ofzof9jm53VDgJeBVoDbsugfh7zwVWAyMSM6PCbvuQdjnecCnk9M1wKaw6z7Jfb4MOB948xjr3wU8TWJ0y4uA1wby/XPxSGE20ODuG9y9E3gUuLHHNjcCDyenfwm808x6Gz40W/S7z+7+grsfSs6+SvbfsjydvzPAV4CvA+2DWVxA0tnnTwAPuvs+AHffPcg1DrR09tmBYcnpSk5whMdM4e4vAXv72ORG4Eee8Cow3MxOHaj3z8VQGA9sTZlvTC7rdRtPDAbUCowalOqCkc4+p7qDxDeNbNbvPpvZecBEd//tYBYWoHT+ztOAaWb2ipm9ambXDVp1wUhnn+8FPpwcsOsp4LODU1pojvf/+3EJdJCdkPT2jb/nKVbpbJNN0t4fM/swUAtcHmhFwetzn82sAPg2cPtgFTQI0vk7F5JoQrqCxNHgy2Y23d1bAq4tKOns823AD939m2Z2MfBIcp/jwZcXikA/v3LxSKERmJgyP4G3Hk52b2N6OmkdAAAClElEQVRmhSQOOfs6XMt06ewzZnYV8M/AHHfvGKTagtLfPg8FpgN/MLNNJNpe52d5Z3O6/7afcPcud98IrCEREtkqnX2+A/gFgLsvBEpJ3CMoV6X1//1E5WIo1AFTzWyKmRWT6Eie32Ob+cBfJ6dvBp73ZA9Olup3n5NNKd8jEQjZ3s4M/eyzu7e6e5W7V7t7NYl+lDnuXh9OuQMinX/bvyFxUgFmVkWiOWnDoFY5sNLZ5y3AOwHM7CwSobBnUKscXPOBv0qehXQR0OruOwbqxXOu+cjdo2Z2J7CAxJkLD7n7CjO7D6h39/nAf5M4xGwgcYRwa3gVn7w09/l+YAjwWLJPfYu7zwmt6JOU5j7nlDT3eQFwjZmtBGLAl9y9ObyqT06a+/y3wPfN7IskmlFuz+YveWb2MxLNf1XJfpL/AxQBuPt3SfSbvAtoAA4BHx3Q98/i352IiAywXGw+EhGRE6RQEBGRbgoFERHpplAQEZFuCgUREemWc6ekigwUMxsFPJecPYXEKZ5Hzn8/5O6XhFKYSIB0SqpIGszsXuCAu38j7FpEgqTmI5ETYGYHkj+vMLMXzewXZrbWzL5mZh8ys9fNbLmZnZ7cbrSZPW5mdcnH28PdA5HeKRRETt65wOeBGcBHgGnuPhv4AUfv2PkA8G13nwXclFwnknHUpyBy8uqO3HvGzNYDzySXLyd5HyLgKqAmZdiOYWY21N3bBrVSkX4oFEROXuodZ+Mp83GO/h8rAC5298ODWZjI8VLzkcjgeAa488iMmc0MsRaRY1IoiAyOzwG1yYHWVwKfCrsgkd7olFQREemmIwUREemmUBARkW4KBRER6aZQEBGRbgoFERHpplAQEZFuCgUREemmUBARkW7/HwSjgPwb7ThPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb05743db38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# Define the right-hand side of a system of ODEs\n",
    "def r(y, t, p):\n",
    "    k1 = p[0] # Forward reaction rate\n",
    "    k2 = p[1] # Backward reaction rate\n",
    "    dydt = k1 * (1 - y) - k2 * y\n",
    "    return dydt\n",
    "\n",
    "# Run an example simulation\n",
    "p = [5, 3]    # parameters\n",
    "y0 = 0.1      # initial conditions\n",
    "\n",
    "# Call odeint, with the parameters wrapped in a tuple\n",
    "times = np.linspace(0, 1, 1000)\n",
    "values = odeint(r, y0, times, (p,))\n",
    "\n",
    "# Reshape odeint's (1000x1) output to (1000,)\n",
    "values = values.reshape((1000,))\n",
    "\n",
    "# Plot the results\n",
    "plt.figure()\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Concentration')\n",
    "plt.plot(times, values)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing a wrapper class for Pints\n",
    "\n",
    "Now we'll wrap the model in a class that extends `pints.ForwardModel`.\n",
    "\n",
    "It should have two methods: one that runs simulations, and one that tells Pints what the dimension of the model's parameter vector is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pints\n",
    "\n",
    "class ExampleModel(pints.ForwardModel):\n",
    "    \n",
    "    def simulate(self, parameters, times):\n",
    "        # Run a simulation with the given parameters for the\n",
    "        # given times and return the simulated values\n",
    "        y0 = 0.1\n",
    "        def r(y, t, p):\n",
    "            dydt = (1 - y) * p[0] - y * p[1]\n",
    "            return dydt\n",
    "        return odeint(r, y0, times, (parameters,)).reshape(times.shape)\n",
    "    \n",
    "    def n_parameters(self):\n",
    "        # Return the dimension of the parameter vector\n",
    "        return 2\n",
    "\n",
    "\n",
    "# Then create an instance of our new model class\n",
    "model = ExampleModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use this model class to run simulations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XXWd//HXJzd72qZL2kLXFGiR0EKBtCwiMMiqUlSQAZcZFK064saMM8zGD9HHbxxxGcbh8dPqMCIuKKLSEbDIIiAWSEo3uqd7uidt0nTJcu/9/P64t+k1pMltm5Nzl/fz8biPnO3e+zlJe9/3fL/nnK+5OyIiIgAFYRcgIiKZQ6EgIiLdFAoiItJNoSAiIt0UCiIi0k2hICIi3RQKIiLSTaEgIiLdFAoiItKtMOwCjldVVZVXV1eHXYaISFZZtGhRk7uP7m+7rAuF6upq6uvrwy5DRCSrmNnmdLZT85GIiHRTKIiISDeFgoiIdFMoiIhIN4WCiIh0UyiIiEg3hYKIiHTLuusURCT7uDtxh1jcibsTizvRuBOPOzE/+jMaO7o+7oltYnEnHoeYO7F4nLhDPO44EHcHh7iDk3gPd8ePzMfp3s6PrEuZjyeHIz463+P5Tvd2ztF1PZ+f2MfkTzxluvflqb+X/rb1lBd/51ljOXfi8AH7u/RGoSCSI+Jxpz0ao70rTntXLPmIJ5fF6IjG6ejqsT6amO6KxYnGnM5YnK5YnK6o0xWP0xVzuqJxovE4ncnprlicrvjR6Wjc6UyZjsYSH9zReDzlw1xjwZ8sMxgzrFShIJLL3J1DnTFaDnfRcqiT1sNd7D/cRVt7lIMdUQ52xhI/O6Ic6EhOdybXdcQ4kJw/1BGjMxY/4ToiBUZRxCiKFCQfieniSAGFPZYXFxZQ0WO7I9OFEaOwoIACS0wXmBEpgIgZkYICIgVQUGDJ+aOPAjMKC+yY6yIFR9cXGBSYYQZGYt7syM/EtHF0mwIz6DGf+An0fD7J17Ue25N83ZT5gkQBJF++e5uj00eWW/c0PZYfmbfudZYyffT1BpNCQWQAtXfFaDrQQdOBTpoPdHRPNx3ooOXQ0Q/+luSHf8uhLqL9fIuOFBgVxRGGlBRSkXwMKSlk9NASKooT8+UlEcqKIpQWRSgpLKC0KEJpUQGlhcllRcllhcnlyW1LiwooKYwQKRj8Dx/JTAoFkTTE407zwU52tB5me0s7O1oPs6O1ne0th9m1vz3xwd/WQVtHtNfnDykpZHh5EcPLi6gsK+LUyjKGlSXmh5cllg0vL2JYcnpoSREVJREqSgopKSwI5Ruj5CeFgkjSgY4om5sPsrn5UPJxkE3NB9ne0s7O1va3NM8UFxYwrrKUscNKqRk3jNFDSqgaUkzVkJLEY+jR+dKiSEh7JXJ8FAqSV9ydXfs7WLurrfuxYc9BNjUfoulAx59tWzWkhMmjypk5cTinzihlXGUZp1aWMm544ufIimJ9g5ecE2gomNl1wANABPiBu3+tl21uAe4lcebVUnf/YJA1Sf5o74qxemcbyxtbWLmjjXW72lizq4229qNNPFVDSjh9dAVXnTWGSaPKqR5VweRR5UweVcGQEn1nkvwT2L96M4sADwJXA41AnZnNd/eVKdtMBf4ReLu77zOzMUHVI7ktGouzemcbSxtbWN7YyvJtrazZ2dbdiVtZVsSZY4dy48xxTBs7tPsxsqI45MpFMkuQX4VmAw3uvgHAzB4FbgRWpmzzCeBBd98H4O67A6xHckh7V4zFW1qo27SXuk17eWPzPg52xgAYXl7EjPGVzL3sNM6ZUMmMCcMZV1mqph6RNAQZCuOBrSnzjcCFPbaZBmBmr5BoYrrX3X8XYE2SpaKxOEsbW3lp7R7+2NDEssYWumKOGZw5dijvP38CtdUjOH/SCCaMKFMAiJygIEOht/+VPU/ILgSmAlcAE4CXzWy6u7f82QuZzQXmAkyaNGngK5WMtGt/Oy+s3s1L6/bwx3VN7G+PYgbnTBjOxy6dwuzqkdROHklleVHYpYrkjCBDoRGYmDI/AdjeyzavunsXsNHM1pAIibrUjdx9HjAPoLa2VtfL57BNTQf53YqdLFixk8VbEt8Nxg4r4dqzT+HyM0dz6RlVDC9XP4BIUIIMhTpgqplNAbYBtwI9zyz6DXAb8EMzqyLRnLQhwJokA21uPshvFm/n6Td3sHpnGwAzxlfyd9dM46qasZw5dqiag0QGSWCh4O5RM7sTWECiv+Ahd19hZvcB9e4+P7nuGjNbCcSAL7l7c1A1SeZoPdTFb5dv59dvbKN+8z7MoHbyCP71PTVce/ZYJowoD7tEkbxk3vNerhmutrbW6+vrwy5DToC78+qGvfz4tc38fsUuOmNxpo4ZwvvOH897Z45n3PCysEsUyVlmtsjda/vbTlfnSOD2t3fxq0WN/Pi1LTTsPkBlWREfvHASN50/genjh6lpSCSDKBQkMI37DvGDlzfy87qtHO6Kce6ESu6/+RxuOHec7gUkkqEUCjLgVu/cz/de3MD8pdsxYM7Mcdx+STXnTAh2cBAROXkKBRkwK7a38q1n1vLc6t2UF0e4/ZJq7rh0ivoKRLKIQkFO2vo9B/jW79fy5LIdVJYVcdfV0/iriyfregKRLKRQkBO2a38731iwhsffaKS0KMJnrzyDj7/jNCrLdIWxSLZSKMhx64jG+O8/buS/nm8gGnM++vYpfPqK06kaUhJ2aSJykhQKclyeXbmLrzy5ks3Nh7i6Ziz/8u6zmDyqIuyyRGSAKBQkLbv3t3PPEyv43YqdnDFmCI/cMZt3TB0ddlkiMsAUCtInd+ex+ka++uRK2qNx/v66M/nEO06jKFIQdmkiEgCFghzTjtbDfOmxZfyxoYnZU0bytffP4LTRQ8IuS0QCpFCQXj29fAd3/2o5XbE4X33vdD44exIFBbodhUiuUyjInznUGeXL81fy8/qtnDOhkgduPY8pVepIFskXCgXptnZXG5/68SI2Nh3kb644nS9ePU19ByJ5RqEgADy1fAd/99hSyosL+cnHL+SS06vCLklEQqBQyHOxuHP/gjV898X1nDdpOP/vQxdwSmVp2GWJSEgUCnnsQEeUz/zkDV5cu4cPXTiJe26ooaRQt7QWyWcKhTy1s7Wdj/6wjrW72vi/75vBBy+cFHZJIpIBFAp5aPXO/Xz0f+rYf7iLh26fxeXTdGWyiCQoFPLMwvXNzP1RPeUlEX7xqYs5e1xl2CWJSAZRKOSRP6zZzScfWcSkkeU8/LHZGvxGRN5CoZAnFqzYyZ0/fYOpY4byyB2zGaXbXItILxQKeWD+0u188edLmDG+koc/OpvKcg2CIyK9UyjkuCeX7eALjy6mtnokD90+iyEl+pOLyLHpEyKHvbB6N59/dDHnTxrBDz86i/Ji/blFpG+6sU2O+tP6Jj7140WcdeowHlIgiEiaFAo5aOnWFj7+cH33WUbDStWHICLpUSjkmK17D3HHw3WMrCjmJx+/kJEVxWGXJCJZRG0KOaT1UBe3/8/rdMWcR+fOZsww3dhORI5PoEcKZnadma0xswYzu7uX9beb2R4zW5J8fDzIenJZRzTG3Efq2br3MPM+cgFnjNGwmSJy/AI7UjCzCPAgcDXQCNSZ2Xx3X9lj05+7+51B1ZEP3J1/fHw5r23cywO3zuTC00aFXZKIZKkgjxRmAw3uvsHdO4FHgRsDfL+89cM/beJXi7fxxaumcePM8WGXIyJZLMhQGA9sTZlvTC7r6SYzW2ZmvzSzib29kJnNNbN6M6vfs2dPELVmrVc3NPPVJ1dxdc1YPnvlGWGXIyJZLshQsF6WeY/5/wWq3f0c4Fng4d5eyN3nuXutu9eOHq3bPB+xveUwn/nJG0weVc63bjmXgoLefuUiIukLMhQagdRv/hOA7akbuHuzu3ckZ78PXBBgPTmlIxrj0z9eREc0zryP1DJU1yKIyAAIMhTqgKlmNsXMioFbgfmpG5jZqSmzc4BVAdaTU/796TUsbWzlm7ecqzONRGTABHb2kbtHzexOYAEQAR5y9xVmdh9Q7+7zgc+Z2RwgCuwFbg+qnlzy/OpdPPTKRm6/pJprzz4l7HJEJIeYe89m/sxWW1vr9fX1YZcRml3727n+gZcZO6yUX//NJZQWRcIuSUSygJktcvfa/rbTbS6ySCzufOHRJRzujPGd285TIIjIgNNtLrLID17ewMINzXz9pnPUjyAigdCRQpZo2N3GN3+/lmvPHssHaieEXY6I5CiFQhaIxuL87WPLqCiO8NX3zsBM1yOISDDUfJQFvv/yRpZubeE/bzuP0UNLwi5HRHKYjhQy3LpdbXz792u57uxTuOGcU/t/gojISVAoZLB43Pn7x5dRURLhK++drmYjEQmcQiGD/axuC4u3tPCv76lRs5GIDAqFQoba09bBvz+9motPG8X7ztPtsEVkcCgUMtS/PbWKw10xNRuJyKBSKGSgheub+dXibXzystN1kZqIDCqFQoaJxuLc88SbTBxZxp0aNEdEBplCIcP87PUtrNt9gH95d43ubSQig06hkEFaD3Xxrd+v5aLTRnJNzdiwyxGRPJTWFc1mNhr4BFCd+hx3/1gwZeWn7zy/jpbDXfzre2rUuSwioUj3NhdPAC+TGEc5Flw5+Wtj00EeXriJWy6YyNnjKsMuR0TyVLqhUO7u/xBoJXnua0+vojhSwN9eOy3sUkQkj6Xbp/BbM3tXoJXkscVb9rFgxS4+efnpjBlaGnY5IpLH0g2Fz5MIhnYza0s+9gdZWD75xjNrGFVRzB2XTgm7FBHJc2k1H7n70KALyVevNDTxSkMz97ynhooS3clcRMKV9qeQmc0BLkvO/sHdfxtMSfnD3bl/wRrGVZbywQsnhV2OiEh6zUdm9jUSTUgrk4/PJ5fJSXh21W6WbG3h81dN1YVqIpIR0j1SeBcw093jAGb2MLAYuDuownJdPO5885k1TKmq4KbzNeayiGSG47mieXjKtE6kP0m/X7WL1Tvb+Nw7z6AwogvLRSQzpHuk8G/AYjN7ATASfQv/GFhVOc7d+a/nG5g0spwbzhkXdjkiIt3SPfvoZ2b2B2AWiVD4B3ffGWRhueyldU0s39bK194/Q0cJIpJR+vxEMrO3JX+eD5wKNAJbgXHJZXKc3J3vPLeOUytLeb/6EkQkw/R3pHAXMBf4Zi/rHLhywCvKca9t3Ev95n3ce0MNxYU6ShCRzNJnKLj73OTk9e7enrrOzPq9H4OZXQc8AESAH7h7r6exmtnNwGPALHevT6fwbPXgCw1UDSnm1tm6LkFEMk+6X1X/lOaybmYWAR4ErgdqgNvMrKaX7YYCnwNeS7OWrLVieysvr2viY5dO0XUJIpKR+jxSMLNTgPFAmZmdR6KTGWAYUN7Pa88GGtx9Q/K1HgVuJHHxW6qvAF8H/u74Ss8+//3yRsqLI3xo9uSwSxER6VV/fQrXArcDE4BvpSxvA/6pn+eOJ9EpfUQjcGHqBsmgmejuvzWznA6Fna3tzF+6nQ9fNJnK8qKwyxER6VV/fQoPAw+b2U3u/vhxvnZvQ4d590qzAuDbJEKn7xcym0uiw5tJk7KzLf5HCzcRc+djb9edUEUkc6V7ncLjZvZu4GygNGX5fX08rRGYmDI/AdieMj8UmA78ITn05CnAfDOb07Oz2d3nAfMAamtrnSxzqDPKT17bwrU1pzBpVH+tbiIi4Un3hnjfBf4S+CyJI4APAP01jNcBU81sipkVA7cC84+sdPdWd69y92p3rwZeBd4SCLng8UWNtB7u4uPv0FGCiGS2dM8+usTd/wrY5+5fBi7mz48C3sLdo8CdwAJgFfALd19hZvclb8OdF+Jx56FXNjFz4nAumDwi7HJERPqU7r2PjlyjcMjMxgHNQL9fe939KeCpHsvuOca2V6RZS1Z5uaGJjU0HeeDWmSSbyUREMla6ofC/ZjYcuB94g0SH8fcDqyqHPLJwM1VDirl++qlhlyIi0q9+QyF5ltBz7t4CPG5mvwVK3b018Oqy3LaWwzy/ehefuvx03dJCRLJCv59UyYF1vpky36FASM+jr2/BQUNtikjWSPfr6zNmdpOpUTxtndE4P3t9K1eeOYYJI3Qaqohkh3T7FO4CKoCombWTOC3V3X1YYJVluWdW7qTpQAcfvki3tBCR7JHuxWtDgy4k1/z41c1MHFnGZdNGh12KiEja0r147bl0lknChj0HeHXDXm6bPYlIgVrcRCR79HeX1FISd0OtMrMR/PldUjW48DE8tqiRSIFxs0ZWE5Es01/z0SeBL5AIgEUcDYX9JMZKkB5icedXbzRyxbTRjBnW7zhEIiIZpb+7pD4APGBmn3X37wxSTVntpXV72LW/gy/P0VGCiGSfdDuav2NmlwDVqc9x9x8FVFfWeqx+KyMrirnybWPDLkVE5LilFQpm9ghwOrAEiCUXO6BQSLHvYCfPrtzNhy6apCuYRSQrpXudQi1Q4+5ZN5bBYHpiyTY6Y3E+cEGfN5AVEclY6X6dfZPEIDjSh8cWNTJ9/DBqxumaPhHJTukeKVQBK83sdaDjyEJ3z5txEfqzeud+Vmzfz7031IRdiojICUs3FO4NsohcMH/JdiIFxnvO1eUbIpK90j376EUzmwxMdfdnzawciARbWvZwd55Ysp1Lz6iiakhJ2OWIiJywdG9z8Qngl8D3kovGA78Jqqhs88aWfWxrOcyNM3WUICLZLd2O5s8AbydxJTPuvg4YE1RR2eaJJdspKSzgmrPVFy8i2S3dUOhw984jM2ZWSOI6hbwXjcV5ctkOrjprLENK0u2iERHJTOmGwotm9k9AmZldDTwG/G9wZWWPV9Y303ywkzlqOhKRHJBuKNwN7AGWk7hJ3lPAvwRVVDZ5Ysk2hpYWcsWZGjdBRLJfuu0dZcBD7v59ADOLJJcdCqqwbNDeFeOZFbt414xTKCnUyVgikv3SPVJ4jkQIHFEGPDvw5WSXF9fu4UBHlBt0bYKI5Ih0Q6HU3Q8cmUlO5/1o9L97cyeVZUVcdNqosEsRERkQ6YbCQTM7/8iMmV0AHA6mpOzQGY3z7KpdXF0zlqKI7ogqIrkh3T6FLwCPmdn25PypwF8GU1J2eGV9E23tUa6frmsTRCR3pHubizozextwJokhOVe7e1eglWW43y3fyZCSQi6dWhV2KSIiA+Z4rraaxdGR184zs7wdeS0ai/PMyp1c+bYxOutIRHJKuvc+egT4BnApiXCYRWLgnf6ed52ZrTGzBjO7u5f1nzKz5Wa2xMz+aGZZcd/p1zfuZd+hLjUdiUjOCWzkteS1DA8CVwONQJ2ZzXf3lSmb/dTdv5vcfg7wLeC6dN8jLE+/uZPSogIu1wVrIpJjghx5bTbQ4O4bkvdNehS4MXUDd9+fMltBFtxPKR53FqzYyRXTxlBerHsdiUhuCXLktfHA1pT5RuDCnhuZ2WeAu4Bi4Mo06wnNksYWdrd1cJ2ajkQkBwU58pr1suwtRwLu/iDwoJl9kMT9lP76LS9kNheYCzBp0qQTKGXgPLdqF5EC4y/O1J3DRST3pNV85O4vAquBocnHquSyvjQCE1PmJwDbj7EtJJqX3nuM95/n7rXuXjt6dLjt+M+t2k3t5BFUlheFWoeISBDSPfvoFuB14APALcBrZnZzP0+rA6aa2RQzKwZuBeb3eN2pKbPvBtalW3gYGvcdYvXONq46a2zYpYiIBCLd5qN/Bma5+24AMxtN4oZ4vzzWE9w9amZ3AgtIjOf8kLuvMLP7gHp3nw/caWZXAV3APnppOsokz6/eDcCVZ6npSERyU7qhUHAkEJKaSeMow92fIjH2Quqye1KmP5/m+2eE51btZkpVBaePHhJ2KSIigUg3FH5nZguAnyXn/5IeH/a57mBHlIXrm/nIxZPDLkVEJDB9hoKZnQGMdfcvmdn7SVzRbMBC4CeDUF/G+GNDE52xOO9U05GI5LD+moD+A2gDcPdfuftd7v5FEkcJ/xF0cZnk+VW7GVpayKzqkWGXIiISmP5Codrdl/Vc6O71JG6Olxficee51bu5fNpojZ0gIjmtv0+40j7WlfWxLqcs39ZK04EONR2JSM7rLxTqzOwTPRea2R3AomBKyjwvrd2DGVw2VTfAE5Hc1t/ZR18Afm1mH+JoCNSSuE/R+4IsLJO8vK6J6eMqGTWkJOxSREQC1WcouPsu4BIz+wtgenLxk+7+fOCVZYi29i7e2LKPuZedFnYpIiKBS3c4zheAFwKuJSP9aX0z0bhz2TQ1HYlI7tOpNP14ed0eKoojnD9pRNiliIgETqHQj5fWNnHx6aMoLtSvSkRynz7p+rC5+SBb9h5S05GI5A2FQh9eWrsHgHfoVFQRyRMKhT68uLaJiSPLqB5VHnYpIiKDQqFwDF2xOAvXN3HZ1NGY9TayqIhI7lEoHMPiLS0c7Iyp6UhE8opC4Rj+tL6JAoOLTx8VdikiIoNGoXAMf1rfzNnjKqksKwq7FBGRQaNQ6EV7V4wlW1p0lCAieUeh0ItFm/fRGYtz8WkKBRHJLwqFXixc30ykwJg1RaOsiUh+USj0YuGGZmaMr2RISVr3CxQRyRkKhR4OdkRZulX9CSKSnxQKPdRt2ks07lyiUBCRPKRQ6GHhhmaKIkbtZPUniEj+USj08Or6ZmZOHE5ZcSTsUkREBp1CIcX+9i6Wb2vVqagikrcUCinqNu4l7nCR+hNEJE8pFFLUbdpHUcQ09KaI5K1AQ8HMrjOzNWbWYGZ397L+LjNbaWbLzOw5M5scZD39qd+0lxnjKyktUn+CiOSnwELBzCLAg8D1QA1wm5nV9NhsMVDr7ucAvwS+HlQ9/WnvirGssZVZ1TrrSETyV5BHCrOBBnff4O6dwKPAjakbuPsL7n4oOfsqMCHAevq0fFsrnbE4tQoFEcljQYbCeGBrynxjctmx3AE83dsKM5trZvVmVr9nz54BLPGo1zfuBeCCyepPEJH8FWQo9DaGpfe6odmHgVrg/t7Wu/s8d69199rRo4MZCa1+017OGDOEkRXFgby+iEg2CDIUGoGJKfMTgO09NzKzq4B/Bua4e0eA9RxTPO7Ub97HrGodJYhIfgsyFOqAqWY2xcyKgVuB+akbmNl5wPdIBMLuAGvp09rdbbS1R9XJLCJ5L7BQcPcocCewAFgF/MLdV5jZfWY2J7nZ/cAQ4DEzW2Jm84/xcoGq27QPQKEgInkv0AED3P0p4Kkey+5Jmb4qyPdPV/2mvYwdVsKEEWVhlyIiEipd0QzUb9pHbfVIzHrrGxcRyR95HwrbWg6zreUws3QqqoiIQqF+U+L6BF20JiKiUGDxlhbKiiK87ZShYZciIhK6vA+FJVtbOGdCJYWRvP9ViIjkdyh0RGOs3L6fmZOGh12KiEhGyOtQWLF9P52xOOdNVCeziAjkeSgs2dICwHk6UhARAfI8FBZvbWFcZSljh5WGXYqISEbI61BYsnWf+hNERFLkbSg0Hehg697D6k8QEUmRt6FwpD9BRwoiIkflbSgs3rqPwgJj+rjKsEsREckYeRsKS7a28LZTh1JWHAm7FBGRjJGXoRCLO0u3tqo/QUSkh7wMhYbdBzjQEWXmRPUniIikystQWLI1MdKaOplFRP5cXobCssZWhpYUMmVURdiliIhklLwMhTe3tTJ9fCUFBRppTUQkVd6FQmc0zqqdbcyYoFNRRUR6yrtQWLurjc5onOnjFQoiIj3lXSi8ua0VgBkKBRGRt8i7UFi+LdHJPHlkediliIhknLwLBXUyi4gcW16FQmc0zqod6mQWETmWvAqFtbva6Iypk1lE5FjyKhTUySwi0re8CoXl21oZWqpOZhGRYwk0FMzsOjNbY2YNZnZ3L+svM7M3zCxqZjcHWQskQmH6OHUyi4gcS2ChYGYR4EHgeqAGuM3ManpstgW4HfhpUHUc0RmNs1qdzCIifSoM8LVnAw3uvgHAzB4FbgRWHtnA3Tcl18UDrANQJ7OISDqCbD4aD2xNmW9MLguFOplFRPoXZCj01nDvJ/RCZnPNrN7M6vfs2XNCxYysKObqmrHqZBYR6UOQzUeNwMSU+QnA9hN5IXefB8wDqK2tPaFguebsU7jm7FNO5KkiInkjyCOFOmCqmU0xs2LgVmB+gO8nIiInKbBQcPcocCewAFgF/MLdV5jZfWY2B8DMZplZI/AB4HtmtiKoekREpH9BNh/h7k8BT/VYdk/KdB2JZiUREckAeXVFs4iI9E2hICIi3RQKIiLSTaEgIiLdFAoiItLN3E/oWrDQmNkeYPMJPr0KaBrAcrKB9jk/aJ/zw8ns82R3H93fRlkXCifDzOrdvTbsOgaT9jk/aJ/zw2Dss5qPRESkm0JBRES65VsozAu7gBBon/OD9jk/BL7PedWnICIifcu3IwUREelDToaCmV1nZmvMrMHM7u5lfYmZ/Ty5/jUzqx78KgdWGvt8l5mtNLNlZvacmU0Oo86B1N8+p2x3s5m5mWX9mSrp7LOZ3ZL8W68ws8DHPw9aGv+2J5nZC2a2OPnv+11h1DlQzOwhM9ttZm8eY72Z2X8mfx/LzOz8AS3A3XPqAUSA9cBpQDGwFKjpsc3fAN9NTt8K/Dzsugdhn/8CKE9Ofzof9jm53VDgJeBVoDbsugfh7zwVWAyMSM6PCbvuQdjnecCnk9M1wKaw6z7Jfb4MOB948xjr3wU8TWJ0y4uA1wby/XPxSGE20ODuG9y9E3gUuLHHNjcCDyenfwm808x6Gz40W/S7z+7+grsfSs6+SvbfsjydvzPAV4CvA+2DWVxA0tnnTwAPuvs+AHffPcg1DrR09tmBYcnpSk5whMdM4e4vAXv72ORG4Eee8Cow3MxOHaj3z8VQGA9sTZlvTC7rdRtPDAbUCowalOqCkc4+p7qDxDeNbNbvPpvZecBEd//tYBYWoHT+ztOAaWb2ipm9ambXDVp1wUhnn+8FPpwcsOsp4LODU1pojvf/+3EJdJCdkPT2jb/nKVbpbJNN0t4fM/swUAtcHmhFwetzn82sAPg2cPtgFTQI0vk7F5JoQrqCxNHgy2Y23d1bAq4tKOns823AD939m2Z2MfBIcp/jwZcXikA/v3LxSKERmJgyP4G3Hk52b2N6OmkdAAAClElEQVRmhSQOOfs6XMt06ewzZnYV8M/AHHfvGKTagtLfPg8FpgN/MLNNJNpe52d5Z3O6/7afcPcud98IrCEREtkqnX2+A/gFgLsvBEpJ3CMoV6X1//1E5WIo1AFTzWyKmRWT6Eie32Ob+cBfJ6dvBp73ZA9Olup3n5NNKd8jEQjZ3s4M/eyzu7e6e5W7V7t7NYl+lDnuXh9OuQMinX/bvyFxUgFmVkWiOWnDoFY5sNLZ5y3AOwHM7CwSobBnUKscXPOBv0qehXQR0OruOwbqxXOu+cjdo2Z2J7CAxJkLD7n7CjO7D6h39/nAf5M4xGwgcYRwa3gVn7w09/l+YAjwWLJPfYu7zwmt6JOU5j7nlDT3eQFwjZmtBGLAl9y9ObyqT06a+/y3wPfN7IskmlFuz+YveWb2MxLNf1XJfpL/AxQBuPt3SfSbvAtoAA4BHx3Q98/i352IiAywXGw+EhGRE6RQEBGRbgoFERHpplAQEZFuCgUREemWc6ekigwUMxsFPJecPYXEKZ5Hzn8/5O6XhFKYSIB0SqpIGszsXuCAu38j7FpEgqTmI5ETYGYHkj+vMLMXzewXZrbWzL5mZh8ys9fNbLmZnZ7cbrSZPW5mdcnH28PdA5HeKRRETt65wOeBGcBHgGnuPhv4AUfv2PkA8G13nwXclFwnknHUpyBy8uqO3HvGzNYDzySXLyd5HyLgKqAmZdiOYWY21N3bBrVSkX4oFEROXuodZ+Mp83GO/h8rAC5298ODWZjI8VLzkcjgeAa488iMmc0MsRaRY1IoiAyOzwG1yYHWVwKfCrsgkd7olFQREemmIwUREemmUBARkW4KBRER6aZQEBGRbgoFERHpplAQEZFuCgUREemmUBARkW7/HwSjgPwb7ThPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb0805195f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run the same simulation using our new model wrapper\n",
    "values = model.simulate([5, 3], times)\n",
    "\n",
    "# Plot the results\n",
    "plt.figure()\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Concentration')\n",
    "plt.plot(times, values)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running an optimisation problem\n",
    "\n",
    "Now that our model implements the `pints.ForwardModel` interface, we can use it with Pints tools such as [optimisers](http://pints.readthedocs.io/en/latest/optimisers/index.html) or [MCMC](http://pints.readthedocs.io/en/latest/mcmc/index.html).\n",
    "\n",
    "In the example below, we use the model to generate test data, add some noise, and then define a score function that characterises the mismatch between model predictions and data. We then use the [SNES](http://pints.readthedocs.io/en/latest/mcmc/index.html) optimiser to estimate the model parameters from the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimising error measure\n",
      "using Seperable Natural Evolution Strategy (SNES)\n",
      "Running in sequential mode.\n",
      "Iter. Eval. Best      Time m:s\n",
      "0     6      11.12974   0:00.0\n",
      "1     12     1.777777   0:00.0\n",
      "2     18     1.777777   0:00.0\n",
      "3     24     0.923      0:00.0\n",
      "20    126    0.412      0:00.1\n",
      "40    246    0.408      0:00.2\n",
      "60    366    0.408      0:00.4\n",
      "80    486    0.408      0:00.5\n",
      "100   606    0.408      0:00.6\n",
      "120   726    0.408      0:00.7\n",
      "140   846    0.408      0:00.8\n",
      "160   966    0.408      0:01.0\n",
      "180   1086   0.408      0:01.1\n",
      "200   1206   0.408      0:01.2\n",
      "220   1326   0.408      0:01.3\n",
      "240   1446   0.408      0:01.4\n",
      "260   1566   0.408      0:01.6\n",
      "280   1686   0.408      0:01.7\n",
      "300   1806   0.408      0:01.8\n",
      "309   1854   0.408      0:01.9\n",
      "Halting: No significant change for 200 iterations.\n",
      "Score at true solution:\n",
      "0.408147155296\n",
      "Found solution:          True parameters:\n",
      " 4.97787481295077772e+00     5.00000000000000000e+00\n",
      " 2.97398166176784784e+00     3.00000000000000000e+00\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4nNWZ8P/vmRmNRr1bVpfcey8YgyEE05IACWAgySaQBFI2pOxu3pdkCQlk8/ttsruBLEk2wIYSlsUQqgkmgAHTjCs2Nu6yLUuyLKt3jTTlvH9M0XSNZI1kae7Pdfmy5plnnjmjcu7ntPsorTVCCCEEgGGsCyCEEOLcIUFBCCGElwQFIYQQXhIUhBBCeElQEEII4SVBQQghhJcEBSGEEF4SFIQQQnhJUBBCCOFlGusCDFVubq4uLy8f62IIIcS4smvXriatdd5g5427oFBeXs7OnTvHuhhCCDGuKKVORnNeTLuPlFJXKKUOK6UqlVJ3hnj+PqXUHve/I0qptliWRwghRGQxaykopYzA74G1QC2wQym1QWt9wHOO1vqHPuffASyOVXmEEEIMLpYthRVApdb6uNa6H1gPXBPh/JuBp2JYHiGEEIOIZVAoAmp8Hte6jwVRSpUBFcBbMSyPEEKIQcQyKKgQx8Jt3nAT8KzW2hHyQkrdrpTaqZTa2djYOGIFFEII4S+WQaEWKPF5XAzUhTn3JiJ0HWmtH9JaL9NaL8vLG3RGlRBCiGGKZVDYAUxXSlUopcy4Kv4NgScppWYCWcCHMSyLEEKIKMQsKGit7cB3gdeAg8AzWuv9Sql7lVJX+5x6M7Bey76gQohhON7YRafVNtbFmDBiunhNa70R2Bhw7O6Axz+PZRmEEBPbS3vqSE00cduaKWNdlAlBch8JIcYtTwdDV599jEsycUhQEEKMW9LpPPIkKAgh6Lc7x7oIw+KUqDDiJCgIEefOdFj5/duVVDZ0jtp7Hm/sorGzD4D2Xhu/3XTU+7i1u5+jZ6IrizOKmHDwdAcHT3cMu6zxRoKCEHHuTIcVgKqmHrTWbD/RgtUWch3piHlpTx3/s9WVtPN4YxdOrfnkVDsAj22p4q97T0d1He2zHra3P3SZ//ZJPX/7pP4sSxw/JCgIMUE0dfVx3xtHONXWO6TXGZQr+YBTa6qae/igsonNh0cvc4Dn/XXYhAfh+fYe/fGdY8PuBjvV1jsi01odTs3xxq5Bz9t+osUbFH396f0TPLOzJsQrRo8EBSHC6LTa2FLZxHhYQrO3to1X3HfXh+uDu0oaOqxsPtwQMWBo8FaMNsfIjzH02R0hZwm5YwLOYbxl4JiCI5r+pBCe2VHDk9uqh/VaXx9UNvHSnrpBA/MHlU3e7jJwzaLSWtPRa+NUqytAvXnwzLA/z9mQoCBEGH/7pJ5tJ1po8PnjPZd8VN3KfW8cwenUvHmwgZbufgBC1edPbqtmd3Ubz+wIvgv13qlr13Vi5X+2VvPwu8fDvn9gBR8qGFttDpq6Bn4egXVmZUOX9/sQLbv7Gxau+2komrtdZbMNocXy5sEz3L/pKLtOtvoca2BvbTs1LT1nXaahkqAgJrSalh5vX/VQee7SzmaGS2t3P332keuf/+RUO0+4ux3eP9oEgCPobnlot9yeO/VIH9Pp1N6xh8G099i8Fa2vjl5XK6SyIXT3SuDbh7pJfnpHDU98ONDtEhg4Nh08w+NbqkJef9vx5pAtoN4RHD+x2V3lMRpC5QMNbW+t6/fzyJng78tYzK6SoCAmtGd31fLGgTPDem00leVgHttSxTM7a4d/gQBvHDhDU2cfWmtvGuLA8g2152fgTn3gmFLQ1tPP41uqqGvrZcuxZv53W7XfXXooTqfmkQ9OsDHCwO7LH/vnxRz4Pg/eFRTYChhK78qWY83srGoNOu5pIQylIvd48+AZDtV30N5jY09NGzZ3QNba9XlOt0c/vqNU8Ndj0XEpQUHErV0nW7wVVHefPahSUngGQM9O0wh0P+2savG7U/etMAPvJu0RWgoqRL03UAH5X2d/XQct3f1UNnTR0Ol67y5r8JhAh7v/2+5westyvLGLhg7roIO3Tqemo9d1zcDg1t1npzvMSuXmrr5h9bmHain0uIOC2TRQHdodTn731lEOhRif8bW3tp1X99Xz7Ee1vH2ogYYO18/a7nSy71Q767fXcLyxi00HzrCnJvJuw6E+y1gMZ8U095EQvg7Vd5CRlEBBRtKov7fWGhVQI757xNX9UtnQxcsf1/H5xUWU56YMnOAdAB2dv8xjjV3kpSWSbkkAXHfFu6tb+dTMSbzn7irysEcICoGPX903ML3TECIqeI4EfkzPZXz7uoPey6lZv72a7j4HkzMszMhP877WM3D7w7UzQn1cAF7bX8+h+k7v+/t2tT3m7ga6an4BU/JSSDAOVNpvH26kpqWHvLTEkNfddryZFRXZQT/zUDzdR4k+QaHDasfm0Hx4rJlZk9MHvUZg8HppTx3FWUnerz0WlWT6ndfhEzR9Z04NlFtzoK6D3DQzk9Isg5ZjJEhLQYyaV/fVs367/0DnyeZuevpjn7fGHqFi99wNdvfbae7qY9fJFnr7HSF3iQr05sEz7KsNPWYx1FlLG/bUsX77wAyYtw65BhvrQnRBOJx6YNZOwNtUNfXQ7NPN46l0Pa97ZkcNu6sHKnrPy33Lq1Ahp4gGBoW9p9rp7nNVqm09tpB94M1dfbx3NPQUV9+yabRfBeqxcd9pfvdWpd9n8gzAhhtU3nKsmUP1neyv8//ZhPqJeFoKdofmRFO3+5jrdzLZbPQ7t6HDyn1vHKGmpccvEIS6y69tHbzrqL59oPWXENBSAXj549O8tr+eJ7ee/cyoaElLQYwZp1Pz/EenyE0183erykf8+r6VnM3h9LvT9D/P9f+OEy209rju3KqaerzdNZHqdtcgYTvzizOCnts3hAFuT2uku8/B24ca2FPTxrRJqUDoZG92R/iWAsDpdis5qYkhF22dauvlVFsveWmJmI2GsIOZoQ6f6ehjal4qSimcTk2DT5eWqx89+DWbDp6hrm3wQWqt4VSEirStN7grymoL31UW6rOHCtSeMYWuPjsv7j7FbWumeANFktlVRbb19GNzaI66B8mf3TW8caK2nn4yk83ex77dWWbjwG1IzwjMhBouCQpizHhmzTR1DW0KoS+rzUGiyRCym8D35s1m12AOOgUYuHv0BASAurZeb+si2tlD/XYn7b026tp6WVCcMaTpnb4ziDx9z5671Cr33auvnSdbsLkDgw5RL3q6iSKld/iLewD88rmTgehmumw/0UJHr41D9Z2YTYagxWKhrhHtz3ew9w+1MG0k5vEHtlT77U7vexmV4vX99eyvc30fs1PC/BJF6bX99dy4vNT72OYT3H0DXOMYToOWoCDGzNn+QTd39fHnD0+ydk4+84qC79R9r9/vviNzOjUv760jIykh4rV9i/bXvaf5/qdTMQTMTvFNBeF0ajZ8XOft1khMCG6VfFDZxJEzndy6usLv+B82V5KaGP5P8eDp4DxAe326rAKnpIJroDPa7itPZezb+tAR1hd7unwCK+mTLd0sKs0MOj/aVcZ9Ee76IfTd81AX2TW6Z2619tg40dRFYWZS0JTUPrvDO2PtSEAOpqGugQgU2LLxnbo72LX31bZTmGkhJzX0OMpIkaAgxkykfv5oeP6ITjR1U56bwuH6TpaUZnpbDb53np4A0WNzcLzR/847VOUZtFJWazbtP8PCkkzy0y309Nt58J2BhViPbami3ad7wzNf3aPTamP7iZaQn6PP5qTPFlwhnGyObuHSmweDp9y+ebAh6srY0z3lG+SOnulicYgKPpKGjj5qW4e/2GqwVcCh1psMNa1FbWsv71c2+U1NLcjwH8D9uGZ461qiEXgj5NtSGMymg2f49OxJMQ8KMtAsxoxjCH8Qg3n54zrePdLond74UXUrH/kMpnruKIebsqLLamd/XQcv7j4F4H0fj/aA/u6qZv/A49sHHe3K2cBrhhNuQPN4iG6nUD481gy47sR9WyzD+U4dro9dptVQd9J9w0jHsafaf2ro6Xb/8Y5YZlR1au0XfOs7/H92aZbI9+kmQ+yrbGkpiDETqttjOJQamD+v3H8z7wQkdPOuTg5Rh0RTjPcrB6aEHq7v9Ju+GErgqt02n/GKJ7ZWcfuaqQAxzUYaadA2FN/uGaVCTMXVmtRERa+1D4N2YNAOlHZg1HaUtmPQDnSTgeye3rDPe4/jBK1RODFop99jpZ0YcIJ2otCux9oJDDxW7tcorSnH97HPOSEeg2ucwOl0Aho1MPfK9bm19juu9MDzrmMaNANfM/C10gRf07vI0HXcaFA0HFB8pt/hPTaPgRH6nBQzLd19rmto32u5JGbcBoXXDunnOlQSFMSYibTIaii0xm8laej30mHfM5qxDU8l39PvYOO+06QPMiYRiWcK56v7TvtNyRyM0nYSHFYSnFZMjl4SnFYSHL2YvP/3keDsJcFhxeS0YnT2Y9Q27/8mZz9Gpw2j9v3f9bXJ2R90PAEHFzjt3krdqMduRsxQaRROZUBjAGVwV9sGQKF9JiW4qmwFSqEJOB5wzFX9+x/3Pu895j7ufY+B4xpQyoBGkxJ0ruu8JG0io8+BVq7XKqVwap9rOWM/AC1BQUTt6R3V9PY7uCVgoHS4RjIDpGeQNNxCM897/dknb45H9TCSjnVE2bUDgNaYHd0k2duw2DpIdHRS9d5ODCfrWOzoItHejdnRhdneTaKji0T312ZHF4mObkwOKwnOXox66Os5HMqIQ5lxGMw4VILrf0OC+5jnsZk+U6r3eadKwGFIIDUlhRarRisjUydnYjIlsO90N05lxKlMOJURszmR4pw0DjX0opURgymBfqfBdQ5GFpTlcfBMD902cBpcr0lPTqSt1+m6T1YGn/8NrJiSS0luKs/sPIXG4D2ulRp47Dnm8/rV0yfx3rEWb8Ufcun2OSDUjC1fa+fk+6VlSU00+U1Jvm56cUzLBxIUxBBEM9d8KOwjOKYwWPK6EU0F7a7kU/qbXP9sTST3t5BkbyfJ1uau/NtJsreRZGvDYm8PrtD3Q7nPw35DEv2mVPqMKfQbU+kzpdKZmE+/MQWbMQmbwYLdaMFmSMJusHiP2YwW7IYkbJ7njBbXuQYLDkMCWvkvvhquGRdWkGZJYNsbR/yOT86wcP6KUp5zH1fKv7U2rSyfK1am89s3jw58+9ItNBhD/y5NnlJOWoqZ3JKUoAkBkZjMFrQ696uzwW6EBsu/ZDTGPtid+99FMSGEGuAd7A/k6JlO/rr3NMVZSUzPTwtKEeDh239f3dITcjDO4dR+ZbAkGEP352uNxd5Oet9p0q2nSe+rI62vgRRb00AQ6G/C7Azur3dixJqQTq8pk96ETFotJZxOm4cpNZcGRypWUwbWhAz6jGn0mQYq/35j8jlfoXnWPVw0My9ovAbgklmTeOtQQ1D3ndGggqbypiSGD1TDSUoHZ79+YLQ4nNqvtTC7IJ1ks9GbSiTw8wc2eBKG+f0ZinP7N1FMGKFu4EMNNNe3W9l2opnPLihkp/sPpba1l9rWXpxa02dzsmpqDhB6ds3mw40hdw2zO7V3PrrBaSen7xSJbcfI7K0mw1rnCgJ9daRb6zE7/buT+g1JdJvz6DbncCZ1Nt3mXLoTcl3/e/4lZNNnShsY6fYxIz8taL77WMlJNdM8jMWCRnftNHtyul9Q8PwIF5Zk8u6RxqBpxsYQ3ThT81LDtgIS3HfC0eQs8hVqXci5qt/uxKAUTq1JMhtZMyPPGxQMSpGbag674G+4QXMoYhoUlFJXAL8FjMB/a63/NcQ564Cf4xrK/1hr/cVYlkmMjVABwLel4HBqPq5tY091G+29tpDTMT2Vkdaa86flcqAu8tRBk8NKTs8xcnsqsZyq4UxnFV/tqSLDWus3aGo1ptFhKaDNUkp1xko6LAV0JLr/WQpdlf1ZsAyjwgpVeScYVdh57SaDayBzsNbXNQuL+MuuGjpDZDuNxFNHmyJ0X/gGhPOn5tDc3U9hZnDyw0j1/WBTLldPy+WDyqag457PnZNqZlFJZtBq8pxUMxdMyw2ZWymUitwUmrr6/L5PRoPyvs/MyWlhp99W5Kag0VQ1hR+r8nRzmgNSrygF65aX8Ie3j4V8nSlMqpaRFLOgoJQyAr8H1gK1wA6l1Aat9QGfc6YDPwZWa61blVKTYlUeMbi2nn6SzEYSTSPTD+1rU4g9DXzHFHZXt/plAo00BrDtRAsrKrL9jqX0NTCp+zB53UfJ6z5KbvcRsqw13il9dpVAW1IJzclTqMz5FP0ZU6kxFtFqKaUvIXg19FBlp5jDrkgdzvcz3ZIQFBTSLAkh3+PW1eWkWRLotTlo7urj+Y9OhbxmSXYyGckJTMlLCVqgFTjAGcjTfWSKcKeam5boTRNemJnEyik5Qed8ZVVZxKnInuvPLUznWENXUOWbnmRien4qR90b0uSkminNTvY+X5CRxILi4KCgcHUZ+vK9TqBrFhVitTn54zsDlbPWcN2SYpITjWFvSNKTEvjcwkKcWvO7tyrDfk6PInfQ9HzvtNYkmoysrMgmPSmBbQELHgebCj0SYtlSWAFUaq2PAyil1gPXAAd8zrkN+L3WuhVAax27vQDFoB79oIq8tES+fF5ZxPPae2xsPtLA5XMnB/2hhRNq6qXvXW1fwIwMq80RdnqpydHL8y88zbLO/Uzu3EdB135S+we6NNosRTQmT+dw3uW0p8+kLnEqHZYCv0HXSE304bhpRQkb9tR5F5JNyRsYKM1PH/oK1Bn5ad6MnR7h6mNPgrXURJPf4rO8tMSgfYBd5wVPp50UJgW1h6fbIrBbx/fhFxYX8VCI7TZ9pSclkGA08O2Lp/LktuqgWVye8YepeanelNs2h9P7vQxsSXzFnUhRa80V8yZ7kwgGcurgrpdI61OUUiSZjX6BDqA0xxWAQqUgB1er0GhQGH2msi4py+Kjk/6b+1w1v4DCTAtp7jTpmUkJ7qDgev78abkAbDnm3yoa70GhCPDNk1wLrAw4ZwaAUuoDXF1MP9da/y2GZRJheO6IoknEtfGT09S3WznTYaUsJyXseU9sPUleqpkr5hX4Ha9p6aEkO9l7xxjqD6yuzerNUmqxtVHc8RHF7bso7thDTnela3ET0GopoSZjKfWp8ziTOovm5Kn0mwYqhjSLKairJCXRSGayf1AozLQMeXbV3MJ0b6K0RJORK+cXePcg/sz8An73diVaQ36GhVtXl/PoB1VRXzs3LXjgNHDAdjA3ryjlP31m/XgG6qdPSvV2wcwpTCct0RR2XwLve/u89cqKbFp6+jl6psuvYk1JNFGem0xVU0/YbizPGIMlwUi0E2muWVTEC7trqWrqCdunrpRidkH4fQ+cWgf9noWLCb6rih1hWqz9jtBrNkL9Ll80Iy8oKMyc7N8lGWr3OwhOgzHUsZbhiGVQCFX6wJ+DCZgOXAwUA+8ppeZprf3WoSulbgduBygtLUWMvKFsWu5Z+WpzOKls6Aq6O2vt7qff4aSps4+mzr6goPDsrlp+uHaGdy/hwN9zk6OXMztfZE37Tkrad5DXXYlCYzNYOJOxkO1Zt3I6bT71aXOxJgTPSCrLSfbmDUo0GQhso6yakkt9wH7DxiGmD7hhWTFFmUlU5KZwzH0X69s/bDIavBVmosng14UUrl/cY3l5NrkpiaQnJfjdSYeqcCpywwdl3/rTd6ObrBQz50/NYcuxZrKSzd6uuC+tLOXJbdVBd8fgXxmdPy2X+nZryK4Xz518uIWJvh/B6P5+XbOocNC+fk/d7Dtw/ZkFBWHODubUwV1foWbEfXr2JOYWDnQn+o6T+KYI7O0P8/miKMunZwf3kofb/c4zS6kiN2XUll7EMijUAiU+j4uBwJ98LbBVa20DTiilDuMKEjt8T9JaPwQ8BLBs2bKRm9wuhsXzx/TXvafRGm5fM4UUn24Lz45ZHoH94J5tDz1jCgqwdFazqO4VKlrfp7j9I0y6H7syU5e+gA9Lb6cmYxn1qXNJS0mOmBMoNdHk98djMhqYOimVYz7TVjU6aPB3sPG7aZNS/aa+Fme5uhGm56cx3b3bWEKYW9/AwUSTUfm1MjwumpmH3aFZXp6FUoqvX1DBbzcd9Q5KTkpL9NuU5bMLCrzvHYqnIg+VT2dJWRZ9dqdf0jtvxa81pdnJw1rUl+Xuygrs5vG0IHyDi2+r4ZsXTYm4N8LC4gxqWnpCtqBC+fJ5ZVQ1dzNrchpP76jhohl5QS2tUN1HBqX8WiP9YVoK4YKeb+A+b0oOk9xdh55pqGtm5DGvMHgMa/XUXHr7HUFB3nNjcM2iwlFpJUBsg8IOYLpSqgI4BdwEBM4sehG4GXhMKZWLqzspcqekGHOe7hjPH5XV5sBqc5CelMDD7wX/+F795LTf4367k7cPniGt/SCrT77A1OZ3yOk9AUCLpZSPC67nRNZq6tIX4jD4d2uEq3g9AmfHGA2KqxcWAq7Bbs/GN0WZSeyk1ee8yFHB9+8xsOk/cI7/e9+wrJiall7v8UWlmd5kbKU5yUFBIS81kRKfQVOAOy6ZBkBnn53GTqtfyuxo7o7WLS8hM0RKjgSjgTUz8vyOeepCp448yyiSVVNzyEtLpCzH/3N8bkFhUIpqz527U2uSzSaSI9T30/PT+OFa1/c9xb3xTbhNk8A1nuLpEvvGhVOA4M2KSrKTgsZtAlOqB2a79Qi7Ktnn2+aZOg2uINXc1ceUvNBjHhnJCVy3NHi18rplxTR39Y9aQIAYBgWttV0p9V3gNVzjBY9orfcrpe4FdmqtN7ifu0wpdQBwAD/SWjfHqkwiNnr6HTy7q5YpeSkhc+IrBu5Ac3qOMaPpDWbs2kS2tRonRmozlrBv8rWcyLqAtqTI3YORKgJwBQHfO0C/zDUDxWBKXiq3ri5nw8d1NHf1+3UlhJr6qVB848IKNHj3UA5lRn6adzZMcVayt0URaNbkdLJTzH7bLIYatPfc3WYkJQSN90SzKU5RiCmh4XhacNkpZu8GP+dPzQk7eBuK0aBCBk2T0UBawM/O89mGmgLrgum55KYmUp4T+nsbTmD30ZLSLBKMBt482EBZTjKrpuYE7R+ebDZ6g4nvt9vz++E7oQDCdx9lJCUMuodHKGmWBO9g9GiJ6ToFrfVGYGPAsbt9vtbAP7j/iXHKs/l4uAVJHc2nWdKwkbkNL5PbcwwnBmozlvJR0Zc4mnNJyHGBcAabpx0YFHyrzYF+W5fMZHPIithsMnD+tGy/RVpKEdUfZ8R+bvcbeyqOwI3Yk8yRZ3IF3iyOUJJZrzRLAtcvLWZSeiIKRU5qIguLM2J2l+rpwgu/nU9oCUZDyO1PB+PbtWNJMKKUosQdtMtyUoICAsC6ZSXUtvXw+n7/6bqe7qMlpVl+v/fhZiWNJ7KiWYTV1NXHEx+e5OKZeSwuzQp63pPn5t0jwYOmSjsoa93KvIYNTGl5F6O2czp1Hm9O+T9U5lxCjzl4Dns0Bus+mluY4Td+EGpFra+r5hdwoK6DPvvA5jtLSrNYUprlFxRyYpRG4epFhWxwD7JaBpluGPhJwgWFaZNSh7w4zcO3+ypcWhHwWcx2FitsPz0rn+xks7dijjVPWctykvmcu0sxK8XM7WumeFtGgTKSE0i1pAcFhcvnTmb7iRZyUv1/L0Zhu4OYk6AggjR39ZGTmsgT7oyiu062hgwKCUbX4JlvDqFEewfzzrzEotPPkN5XT48pkz0F69iffzXNyVPPumyRuo9SEo0sLM7wGxD2/SP1dGP5dhWlJppYUZHNyeZudla1cuPykqBVuOuWl1CQ7n9XP1Km5qXy9QsrqG+3DtoK8tyFFmRYyEtLZEZ+6G4dT4UXS5PSEllRkT2sO3aPJLPROx9/NBgMiq+triAl0ej3vU6JsBUqDIy1+K4ByU+38LmFhUH7d6uo5h+d2yQoiCCN7qDgEW7bTN+b8KyeKhadfpq5DX8lwWmlJn0p75b/gGPZa3AaBrpdQq0b8DVrclrEPQY8d3tFWUl8auYk3q9s9KYTUO7889MnpXr3Sp6Z7zN3PaD7yFdZTgp3XDItZMU8lH754Ui3JEQcp/AoyU5mYUkGy8uzR72fOZBSitWjWKGPlIzkoX/flFJcOX9yyO6lwCAwAXqPJCiIYEfPdPnlbQm3EMlm1+R1HWZl7Z+Y3vw2dpXAobwr2FNwI42pM0O+piwnxW+v3cBFY5fPnUx+hoV3DjeSkZQQNP3Uc7dsNhrIS0ukMCPJW1bPwOuC4gzmFKaj8B+DSHA3G8J1KY1GXpmzYTQoLpmVP9bFiEuzJodeGBfYezaas4RiRYKCCBK4lWS/3RkUGPI793NezX8zpfV9rMZUthZ/nY8Lbhh0rCBwwdDqabn8Zadr/+LPLSzEYFDeyjvUtEhPq8XznO/VPF8rpUKOPXgWac0tDL/yNZY8A6oToeIQLoE/yyWl0U+aOFdJUIhz2443B83fDuXpHa6MJTk9x7ig6ndMaX2fXlMGH5R+iz0FN/qllogkNy1w3YGBqxcV8smpdqbmuRbueCp8h1Mzc3IaiSaDd36+Z9aHZ3FURW6Kd+P5wWbjmE0GLpg+/ro8xLkr3Krx8UyCQhz46946OnrtfHFl8BqALceiWxbS1VjDpTUPMvfMy9iMSbxf9h32TF6HzRQ+zUIgpWBxSabfrB6jQTE1L5WpPot6EnyCwlXzXVM8PUHB04frmS2Sn27hu5dM43dvVVKUNbJ9/0vLsuiOImBGyzUFtZ3sSKu0xLjiaSmYRyFR3WiRoBAHAnPUePYxjibBmtFhZfmpP7Ps1BMYtJ09BTeyreRrEdcWhMv7n5eWGNTcDjWl0dMK8F2c5RmgnpGfSl5aol96hgSjgb9bVTasxUGRBK74PVtzC9MpzEwaN7uEieh8evaksIsUxyMJCnGkuauPJ7dV43Bq0iwm7/L/cMpbPuCS478mo6+OwzmX8kH539NuGXzj8JmT0znR1EV33+BJ9kJlvfR0H/kOY+SlJdJptaOUCtpLASA3dejpqUebUkoCwgS0oHj8jyP4kqAQR/aeavcOGEeaFprWV89Fx/+D6S2baU7C6tP+AAAgAElEQVQq5y9z/4vazGVRv4/WmhUVObx9aPDtMULttOVpitt88st4ksqNwm6EQsQ1CQpxJDAvUWAyMLRm3pkXuejE/SgcvF/29+wq/JLfOoNoJJgMLCrJZFFJJve9cSTiuaFaCknu1BO+6yM+NWsSGckJlEfYv0EIcfYkKExwTV0DSdQOnvbPyvni7oFtG1P7znBp5S+paPuQ6oxlvDHtLjosRYArRfOB0x1hcxv5KspKYvXU6Gf4hBpTSAqRj8iSYOT8IVxXCDE8E2fIPM5ZbQ56+l1dQg6nptq9ycyugB2fQpnV8Cp/t/smijt289aUH/Hc3N97AwK4KuTL506OqhyrpuSEnYmxKsSevaEGu8/1RWRCTGTSUpggHnznOE6t+eHaGbxf2cRHJ1u5eUVpyN2lPEwOK586/mvmNbzMqbSFvDb9Z7QnlYQ815Jg5FsXTeXA6Q4sCQZvgjCTQfl18wRW8lfOn8yr++oBgnLJrwwxYOxxxbzJMUtCJ4QIT4LCOFHb2sNfdtby9QsrQubJ8UzftDmc3v1gn9pezawwG8Jk9xzns4fuJLu3iq3FX2Nr6W1oFfrXwZPSOclsZGlZFnaH0xsULAlGv8VvgSkkZk1OJynBiM1nByuzyYDWOmIytEj77QohYkeCwjixz71461RrL+kF4Qd+d1b5dxeFSls0s/FvrK38Jf3GZJ6f+wDVmSvDXq8oKylouqfvWoOsFLNfUAiVwaEsYHD4m2siT4UVQowdCQrjRLTpcgL3jvXdwERpB+ef/C9WnHqc2vTFbJz5/9Ftjjx4GypDqG8P0WcXFNDea2PDnjq6+uxRJQ6WMQMhzl0SFMaZwfL7BCau8zw027u48shPmdL6Pnvzv8DbU/4pqqmmoSp5FbCDlSXBSJLZGFUOJSHEuU2CwgSz270xvMexhi7S+uq59sD3ye45yZtT/g97C24I+/o5hekcqe/0Dh5HyuiZ55PczhuKZHGZEOOatOPHIavNwRsHztBvH3zH85zuSm7c+3VS+xp4fu4DEQMCuDYYv3L+wPTTcDHhSytLuX7pQMoLz0whs3QNCTGuSUth3BionT862conp9rJTE5geXn4aZ1F7bu5+uA/YjdaeGb+wzSnTBv0XcpzUshPTyQ7xUxLd3/YG/9JAdtTXjo7n7mF6WRKBlAhxjUJCuOMRnu7dDy5gULtjDal+R0+c+Sf6UiczPNzHqDTUjDotX3zwU/NS6WluyWqTKrgmmYaOMtICDH+SFAYhzz7DdjcwWDrcf89EaY1vcVVR35CQ8osXpxzf8Q01+F41j3IEIEQ8UU6gMeh9442AWB3OKls6GL7iRbvczOa3uAzh3/CmdS5PD/3d34BYWFJBmU5/nnfv33x1JDvMbC15ciWXQhxbotpUFBKXaGUOqyUqlRK3Rni+VuUUo1KqT3uf9+IZXnGM0/l7Dsl1e7UvPxxnffxjMbXufLwTzmdPp/n5/xnyC0yAyt5S4jkc673kf2EhYhHMes+UkoZgd8Da4FaYIdSaoPW+kDAqU9rrb8bq3LEi6nNm7nyyE+pS1/Ii3Pux2YMvROUirJDyBN8JCQIEV9i2VJYAVRqrY9rrfuB9cA1MXy/uOCbfM5TYRe37eSqwz+hPm1uxIAAkJEc3d4IntQWslOYEPEllgPNRUCNz+NaIFSSneuUUmuAI8APtdY1gScopW4HbgcoLQ3efD4e1LX1AnCiaWC/ZYNS5Hce4JqD/0hbUgkvzr4vYkAAuHBaLiVZSSSajCQmhL8nmFeUTn5GonuzeSFEvIhlSyFUz0Pg3MmXgXKt9QJgE/B4qAtprR/SWi/TWi/LyxvZzdTHi7YeGwBVTT3eY0kdx7n2wPfpTcjk+TkPkJg++CY0JqOBaZPSKMlOjljhK6UkIAgRh2LZUqgFfJPzFwN1vidorX3nUj4M/CqG5Rm32nttQccstjbmbv46KMXzcx+gO3ESX1tSTE1rDymJJkwGRbLZyJ8/PDkGJRZCjFexDAo7gOlKqQrgFHAT8EXfE5RSBVrr0+6HVwMHY1iecUdrzfYTLWw55r8Owejs5+pDPyK1r5Fn5/0XbUmlzClMJyM5gYzkDO95np3YhBAiWjELClpru1Lqu8BrgBF4RGu9Xyl1L7BTa70B+J5S6mrADrQAt8SqPONRV589KCCgNZdW/pKijj28MuOXnE5fAARvbgOuMYeAlwohREQxXdGstd4IbAw4drfP1z8GfhzLMoxnodJXLK99jDmNG9lScjtH8i7zHjeGSEcRGBR8dz/zdcOy4pCvF0LEn6iCglIqD7gNKPd9jdb6a7EploDgsYSy1i2srv4vDuVezrYS/3V+ZlPwnIHAir7fEbqpUJwVecaSECJ+RNtSeAl4D9cMIUfsiiN8Pf/RKe/XadbTXHnkbpqSp/LGtLv8liavqMgOmS018OY/mlTbQoj4Fm1QSNZa/9+YlkSE9YX5OST+z1cwaDt/nfUrCvKyqWlxTU39xoUVpFlCL0gLTFERrvtICCE8ol2n8Fel1FUxLYnw+uRUOw+9e8z7ePKWnzO56yB/m34PbUmlXDY33/tcuNxFoSwsHnq2VCFEfIm2pfB94CdKqX7A09GttdbpsSlWfHvjwBkAks1G5rZuInHvn9le9FWO51wEQLpPy8AU5QCx714JQggRTlRBQWudFuuCiAFKuaaPGjtPsXL/L6B4BVtKvxXm3MGDQqR0FkII4SvqKanu9QRr3A83a63/Gpsixa/6ditpFpMrk6m2c8WRu1Fo+MJD6F2uBtqM/KHF5y+uLCUlUfZSEkJEJ9opqf8KLAeedB/6vlLqAq110B4JYvie2l4NuKaSLq39M8Udu6la8xvKsytw5QuEy33GE6KRny75i4QQ0Yv2FvIqYJHW2gmglHoc2A1IUBgh2me5cW77J6yqeZBDuZfhmHmd33kmo6sr6LMLCuiwShoLIcTIGkq/QiauVBQAGZFOFEP3oTudhdHZz+VH76EnIYe3pt7J5ebQP6LpQ+xGEkKIaEQbFP5/YLdS6m1cKbHXIOkpRoTWmn6Hk23ufZZX1D5KTu8JXphzP32mNBJDrFQWQohYiXb20VNKqc24xhUU8H+11vWxLFi82FPTxubDjQDkdh9lee2jHMi7iqqs1QAkDWEdghBCnK2IQUEpNUtrfUgptcR9qNb9f6FSqlBr/VFsizdxdfXZMSrFscZuAJS2c9nRX2A1ZfBOxQ+95yWZXUHhaxdUSCpsIUTMDdZS+Adc22D+R4jnNHDJiJdognM6Nafaenl2Vy1KDSSjW3rqf8nvPshfZ/4r1oSBlccWkysoZCQlkJEU3f7KQggxXBGDgtb6dveXV2qtrb7PKaVkruMwbD3e7B0/8Ew4SrOe5ryahziafTFHcz/td75BUloLIUZRtAPNW4AlURwTg2js6vN7rICLTvwGUGye8o/e42vn5FOQIXFXCDG6BhtTmAwUAUlKqcW46jCAdECS8A9D4MY5aTVvM71lM++X/T1diZO9x1MTTeSkJo528YQQcW6wlsLluLbILAZ+43O8E/hJjMo0ofkGBaOzn+WHfkVLUhm7Cr/kd57shCaEGAuDjSk8DjyulLpOa/3cKJVpwnI6tV9QWHrqf8iy1vDc3N/hNPgPIstYghBiLES7TuE5pdRngLmAxef4vbEq2ETT0t3P41uqvI9T+hpYUfsoR3MuoTpzJQDnTclh63H3yuYosp8KIcRIi2q5rFLqj8CNwB24xhVuAMpiWK4Jp7a1x+/x+dUPorSDd8u/5z22amqOd69lgyxkFkKMgWirnvO11l8BWrXW9wCrgJLYFWvi6ekf2No6t/socxteZk/BOjosRQDkpJqBgU1zpKUghBgL0QYFzxqFHqVUIa7d1ypiU6SJqd8+sD/yhVX/SZ8pje0lX/Meu35pMTAwwCwDzUKIsRBtUHhZKZUJ/BvwEVAFPDXYi5RSVyilDiulKpVSYdNsK6WuV0pppdSyKMszrvT0270DzGWtH1LetpWtxV+nzzSwm2myOxuqJxhoHXwdIYSItUEHmpVSBuBNrXUb8JxS6q+ARWvdPsjrjMDvgbW4cibtUEpt0FofCDgvDfgesG2Yn+GcprXmwXeOux84ubDqAdosRewtuN57ziWzJnm/rshNYXd1m3dsQQghRtOgNY97Y53/8HncN1hAcFsBVGqtj2ut+4H1wDUhzvsF8GsGuqjGvT77wPiB3WcK6oymTeT1HGVL6bdxGMze476NgjXT8/jaBRWyhaYQYkxEezv6ulLqOhXNLvEDioAan8e17mNe7lXSJRNpv+eTzd384e1j1LS4ZhvZHK6xBKXtrKp5iKbkKRzOXRv29QaDksR3QogxE21Q+AfgL0CfUqpDKdWplOoY5DWhAoj3ptjdLXUf8I8hzvO/kFK3K6V2KqV2NjY2RlnksVHb2gvAqTbX/za76yPPanyN7N6TfFjyTVDSNSSEODdFVTtprdO01gattVlrne5+nD7Iy2rxn7ZaDNT5PE4D5gGblVJVwHnAhlCDzVrrh7TWy7TWy/Ly8qIp8phTwPHGLmpaezA47ZxX8zANKTNIW3xt0LmzC2RrTSHEuSHaxWtvRnMswA5gulKqQillBm4CNnie1Fq3a61ztdblWutyYCtwtdZ6Z9SlP4e1dPfz0p463jhwhtmNr5BpPYV57V0sr8j1O2/V1BwSTbK7mhDi3DBYllQLrmyouUqpLPyzpBZGeq3W2q6U+i7wGmAEHtFa71dK3Qvs1FpviPT68cozlbTX5hpsNjhtrKz5E/35i8lceDW9toH1ClfOn8z0SdJKEEKcOwab4vJN4Ae4AsAuBoJCB67pphFprTcCGwOO3R3m3IsHu9544hmSn9X0Ghl9p7Fe9BtQynvcoBSzJg/WAyeEEKNrsCypvwV+q5S6Q2v9wCiVaVzT7rH0PpsTtJNltX+mKWUGubOvBFzBAMAoY81CiHNQtFlSH1BKnQ+U+75Ga/3nGJVr3DvdbmVqyzvk9J5g78rfkOsOBp7sFdJKEEKci6IKCkqpJ4CpwB7AszJLAxIUAnjXqmnNitrHaLMU01hyhfd5k9HAty6aSqKsWBZCnIOiXTa7DJijtWTkCae338HhM5043VGhpH0Hk7sOsGnqj7Ek+m+rmWSW2UZCiHNTtEHhE2AycDqGZRm3tNb88Z1jAExKdwWAFbWP0ZWQy4FJn+XbU7LHsnhCCBG1aINCLnBAKbUd6PMc1FpfHZNSjTO+eyX09jvI6zpMafsO3iu7A4fBTIKMKgshxolog8LPY1mI8czh1Lxx4Iz3cZ/dyarT6+k3JLFv8ufHsGRCCDF00c4+ekcpVQZM11pvUkol41qQFvcOnu7gRFO397Gpp5GZja/xSf619JlkYZoQYnyJNs3FbcCzwIPuQ0XAi7Eq1HgSOPQ+v/55TNrG7sIbx6ZAQghxFqLt7P57YDWulcxorY8CkyK+Ik74bptpdPazsP5Zjmetpi2pbAxLJYQQwxNtUOhzb5QDgFLKhP/eMHHLNyjMaHqDFFsLuwtvHsMSCSHE8EUbFN5RSv0ESFJKrcW1t8LLsSvW+OENClqzuO4pmpMqqM5YMbaFEkKIYYo2KNwJNAL7cCXJ2wjcFatCjQc2h5P73jjC7upWAAo695HffZg9BTcOZMMTQohxJtopqUm4Ul8/DKCUMrqP9cSqYOe6LqsdGNhpbf6Z5+k3JHMw74pILxNCiHNatC2FN3EFAY8kYNPIF2f8sNoHFqwl2juY3byJg5OuxGZKASAzWfZZFkKMP9G2FCxa6y7PA611l3utQtyy+myWM7thIwZHH/vyBxarfW5hIR9UNjGnQLKhCiHGj2iDQrdSaonW+iMApdRSoDd2xTr39fS7uo/QmgX1z2MvWEJj6kzv87mpiVyzqGiMSieEEMMTbVD4AfAXpVSd+3EBENers3rd+Y6KOvaQ03uCnot/4F7FIYQQ41e0aS52KKVmATNxbcl5SGtti2nJznGePZjnn3mePmMKjtnXwrYGAD6zoGAsiyaEEMMWbUsBYDkDO68tVkrF9c5rPf0OEm3tTG96i0/yr2F6YiqWhGasNgfTJ6WOdfGEEGJYZOe1YbI5nMxqeg2T7ueT/GuYqRRfXFGKMoCSdQpCiHFKdl4bhkP1HRw908XNDa/QmDydxtSZKAUZMg1VCDHORbtOwbPzmgBe3VdPds9xJncdYP+kzwJgkNaBEGICkJ3XoqS15oXdp1hQnAnAnIZXcGLkcN7lABgkJgghJoCY7rymlLoC+C2uDXn+W2v9rwHPfwtXWm4H0AXcrrU+MJz3ijWnhpPNPZxs7kFpB7MbX6UqaxU95hxAWgpCiIkhqu4jrfU7wCEgzf3voPtYWO78SL8HrgTmADcrpeYEnPa/Wuv5WutFwK+B3wyx/KPG7hxYwVzatp3U/kYOuLuOQHLgCSEmhmh3XlsHbAduANYB25RS1w/yshVApdb6uHsvhvXANb4naK19l3ulcA7v0eATE5jT8ApWUzrHsy/0HpMZR0KIiSDa7qN/BpZrrRsAlFJ5uBLiPRvhNUVAjc/jWmBl4ElKqb8H/gEwA5dEWZ5RVdPSQ7rFNbPIbO9iWstm9k/6LA6DeYxLJoQQIyva2UcGT0Bwa47itaFunYNaAlrr32utpwL/lzB7NCilbldK7VRK7WxsbIyyyCOjocPKs7tq2XzE9fGntWzG5OzjwKTPeM/JTUsc1TIJIUSsRBsU/qaUek0pdYtS6hbgFVwb7URSC5T4PC4G6sKcC67upWtDPaG1fkhrvUxrvSwvLy/KIo8MTzbUxk7XpKsZja/TnlhIfeo8AIqykrhsTv6olkkIIWIlYveRUmoakK+1/pFS6gvABbhaAB8CTw5y7R3AdKVUBXAKuAn4YsD1p2utj7offgY4yjnGM1TQabVjsbVR1radnUVf9j6xbllJhFcLIcT4MtiYwv3ATwC01s8DzwMopZa5n/tcuBdqre1Kqe8Cr+GakvqI1nq/UupeYKfWegPwXaXUpYANaAW+epafZ8T5jh9Pb34LAw4O510GwLrlEhCEEBPLYEGhXGu9N/Cg1nqnUqp8sItrrTcS0M2ktb7b5+vvR1fMseM7q2hm4+s0J5XTlDwdgKLMpHAvE0KIcWmwMQVLhOfiokb0hISUvkaKOz7iSO5lWMwmfnDp9DEtlxBCxMJgQWGHUuq2wINKqa8Du2JTpHPTjOZNKDSHc9ditTlkXYIQYkIarPvoB8ALSqkvMRAEluFaU/D5sK+aQJzuxLAzmt6gIWUGrcnlY1sgIYSIoYhBQWt9BjhfKfUpYJ778Cta67diXrJzhNaQbq2jsHMfjkt+htlhoDgrLnrOhBBxKNrtON8G3o5xWc45NoeTZ3fVsqTZ9dGN867l25kVkudICDFhRbt4LS71210L16a1bKYxeTpkT8FgUDKeIISYsCQoRKCB5P5mCjs+pjLn4rEujhBCxJwEhQgcTs2UlndRaFpLLx/r4gghRMxJUAihvdfGfW8coaalh2nNb9NmKWL5ygsHf6EQQoxzEhRCqGnpAeCdvccobd9BZfbFGGS/TSFEHJCgEEFF6wcYtZ1jOZ+S7TaFEHEh2k124sZzu2qpdrcUpjVvpjshm7q0+TINVQgRF6SlEMATEIzOPspbt3As+yJQBpmGKoSICxIUfDicAxvDlbbtwOzs8U5FlSEFIUQ8kKDgw5PnCKCi9X36DUnUZiwDkDEFIURckKDgwxsTtKai9QOqM1fgMJgBZExBCBEXJCj4ONHUDUBOzzHS++o5kbXa+5xPI0IIISYsmX3kY+O+04BrKipA0pwrWJObx4mmbpISjGNZNCGEGBUSFEKoaP2AhpQZXLB0IQBLy7LGuERCCDE6pPsoQKK9g8KOvZzIumCsiyKEEKNOgoJbb78DgLLWrRhw+I0nCCFEvJDuI6Clu5/Ht1QBMKX1fXpMmdSnzR3bQgkhxBiI65ZCY2cfWmtae/oBUNpBeesWTmatQisZWBZCxJ+4DQo1LT38z9aTfFzbjnbPN53cuZ8kezvHZTxBCBGnYhoUlFJXKKUOK6UqlVJ3hnj+H5RSB5RSe5VSbyqlymJZHl/tvTYAGjqs3jUIZW1b0ShOZq4crWIIIcQ5JWZBQSllBH4PXAnMAW5WSs0JOG03sExrvQB4Fvh1rMoTjlIKz7q0srZtnEmdTV9CxmgXQwghzgmxbCmsACq11se11v3AeuAa3xO01m9rrXvcD7cCxTEsjx/fFcpOrTHbu5jcuV9aCUKIuBbLoFAE1Pg8rnUfC+frwKsxLE9YWkNx+y4MOKiWoCCEiGOxnJIaKoVcyAxCSqkvA8uAi8I8fztwO0BpaelIlc/LqTVlbVvpNyRxOm3+iF9fCCHGi1i2FGqBEp/HxUBd4ElKqUuBfwau1lr3hbqQ1vohrfUyrfWyvLy8ES+o1lDWtp3ajCXerKhCCBGPYhkUdgDTlVIVSikzcBOwwfcEpdRi4EFcAaEhhmUJot2NFgWYOmrJslZTnbliNIsghBDnnJgFBa21Hfgu8BpwEHhGa71fKXWvUupq92n/BqQCf1FK7VFKbQhzuZhp6emnZucrADKeIISIezFNc6G13ghsDDh2t8/Xl8by/aNxqrWXhW3b6DLn0Zw0BYDZBeksKJZpqUKI+BO3K5o9lHZQ2r6Dk5krvNurlWYnU5iZNMYlE0KI0Re3QcGzTmFS12GS7O1UZwx0HWUkJ4xRqYQQYmzFbVDwKGnfCUB15nIAlpdnUyStBCFEnIr7oFDcsYvmpAp6zLkAJCdKdlQhRPyK66CgtJ2i9j3UZizxHtMhl9cJIUR8iOtNdiZ1Hcbs7CF91qfAlTSVhTLrSMQJm81GbW0tVqt1rIsiRpDFYqG4uJiEhOGNjcZtUNBASfsuAHoKz4OTriaCyRjXjScRR2pra0lLS6O8vBylQmWlEeON1prm5mZqa2upqKgY1jXiugb0jCfolEljXRQhRp3VaiUnJ0cCwgSilCInJ+esWn9xGxScDlvQeIIQ8UYCwsRztj/TuA0KSU2fYHb2UJuxdKyLIoSIwt13382mTZtG/Lo7d+7ke9/73ohf9+c//zn//u//PuLXjbW4HVNIO/0hALXpSyjVML8og6wUWbQmxLnq3nvvjcl1ly1bxrJly2Jy7fEoblsKqfXb3OsTctBoLp2Tz9Ky7LEulhBxoaqqitmzZ3Pbbbcxd+5cLrvsMnp7ewHYs2cP5513HgsWLODzn/88ra2tANxyyy08++yzANx5553MmTOHBQsW8E//9E90dnZSUVGBzeaaRtjR0UF5ebn3scdf/vIX5s2bx8KFC1mzZg0Amzdv5rOf/SwAjY2NrF27liVLlvDNb36TsrIympqaIpb34YcfZvny5SxcuJDrrruOnp4exrP4bCk47KSd2cEnuVcAsjZBiM2HG2jsDLmdybDlpSVy8czwkziOHj3KU089xcMPP8y6det47rnn+PKXv8xXvvIVHnjgAS666CLuvvtu7rnnHu6//37v61paWnjhhRc4dOgQSina2tpIS0vj4osv5pVXXuHaa69l/fr1XHfddUHTMu+9915ee+01ioqKaGtrCyrTPffcwyWXXMKPf/xj/va3v/HQQw8NWt4vfOEL3HbbbQDcdddd/OlPf+KOO+4422/fmInPlsLpjzHau73jCRIThBh9FRUVLFq0CIClS5dSVVVFe3s7bW1tXHSRaxPGr371q7z77rt+r0tPT8disfCNb3yD559/nuTkZAC+8Y1v8OijjwLw6KOPcuuttwa95+rVq7nlllt4+OGHcTgcQc+///773HTTTQBcccUVZGVlRSwvwCeffMKFF17I/PnzefLJJ9m/f//ZfFvGXNy1FGpaejDt3UQBUJu+GHDN7RUinkW6o4+VxMRE79dGo9HbHTMYk8nE9u3befPNN1m/fj2/+93veOutt1i9ejVVVVW88847OBwO5s2bF/TaP/7xj2zbto1XXnmFRYsWsWfPHr/nI9UF4cp7yy238OKLL7Jw4UIee+wxNm/eHNXnOFfFXUvh2V21dFd+QEdSiTffkRDi3JCRkUFWVhbvvfceAE888YS31eDR1dVFe3s7V111Fffff79fxf6Vr3yFm2++OWQrAeDYsWOsXLmSe++9l9zcXGpqavyev+CCC3jmmWcAeP31173jGZF0dnZSUFCAzWbjySefHNLnPRfFXUsBrSns3EttzgXeQyZD3MVGIc5Zjz/+ON/61rfo6elhypQp3i4hj87OTq655hqsVitaa+677z7vc1/60pe46667uPnmm0Ne+0c/+hFHjx5Fa82nP/1pFi5cyDvvvON9/mc/+xk333wzTz/9NBdddBEFBQWkpaXR1dUVtry/+MUvWLlyJWVlZcyfP5/Ozs6z/A6MLTXeuk6WLVumd+7cOezXP7rhDW796HrenvHP7Mm7FoDvfXo6RoMs4hHx5eDBg8yePXusizGinn32WV566SWeeOKJYb2+r68Po9GIyWTiww8/5Nvf/nZQF9N4EOpnq5TapbUedO5tXLUU2nr6Kez4GIDqlAWcNyWHVVNzxrhUQoiRcMcdd/Dqq6+ycePGwU8Oo7q6mnXr1uF0OjGbzTz88MMjWMLxIa6CwqMfVHFpx16spnRaksoxm6TbSIiJ4oEHHjjra0yfPp3du3ePQGnGr7irFQs7P6YubQEoA2bJiCqEEH7iqla02NrI6a2iLn0BgLQUhBAiQFzVioWdewGoS1sISFAQQohAcVUrFnZ8jEOZqE+dA0hQEEKIQDGtFZVSVyilDiulKpVSd4Z4fo1S6iOllF0pdX0sywKulkJDyiwcRgsACUaZhiqEGLqhpMUeb6m5Yzb7SCllBH4PrAVqgR1KqQ1a6wM+p1UDtwD/FKtyeNn7yO88wMcFA7FHIUFBiLGmtUZrjWGMF5Ha7XZMppGvEqhlJRoAAAq8SURBVMdbau5Y/hRWAJVa6+Na635gPXCN7wla6yqt9V7AGcNyuN6rbg8m3U9d+kLvsaxk2T9BiLHgSUX9ne98hyVLllBTU8Prr7/OqlWrWLJkCTfccIN3FfGOHTs4//zzWbhwIStWrKCzsxOr1cqtt97K/PnzWbx4MW+//TYAK1eu9EtId/HFF7Nr1y66u7v52te+xvLly1m8eDEvvfQSAI899hg33HADn/vc57jssssA+Ld/+zeWL1/OggUL+NnPfua91i9/+UtmzpzJpZdeyuHDh0N+romQmjuW6xSKAN/EIrXAyhi+X0TO6q0YGRhkXjsnH5NMSRXC5dU7oX7fyF5z8ny48l/DPn348GEeffRR/vCHP9DU1MS//Mu/sGnTJlJSUvjVr37Fb37zG+68805uvPFGnn76aZYvX05HRwdJSUn89re/BWDfvn0cOnSIyy67jCNHjnDTTTfxzDPPcM8993D69Gnq6upYunQpP/nJT7jkkkt45JFHaGtrY8WKFVx66aUAfPjhh+zdu5fs7Gxef/11jh49yvbt29Fac/XVV/Puu++SkpLC+vXr2b17N3a7nSVLlrB0afCujRMhNXcsg0Kovplh5dRQSt0O3A5QWlo6rMLYp1/Jm1U2esyyglmIc0FZWRnnnXceAFu3buXAgQOsXr0agP7+flatWsXhw4cpKChg+fLlgCttNrhSXHsqxlmzZlFWVsaRI0dYt24da9eu5Z577uGZZ57hhhtuAFzJ7TZs2ODtg7darVRXVwOwdu1asrOzvee9/vrrLF7syqDc1dXF0aNH6ezs5POf/7w3TffVV18d8jN5UnOvW7eOL3zhC0HPv//++7zwwgvA0FJz33XXXbS1tdHV1cXll18+tG/0EMUyKNQCJT6Pi4G64VxIa/0Q8BC4ch8N5xq2zCnszx/4Qcp+5UL4iHBHHyspKSner7XWrF27lqeeesrvnL1794bciD5czraioiJycnLYu3cvTz/9NA8++KD3/Oeee46ZM2f6nb9t27agcvz4xz/mm9/8pt95999/f8hyBJoIqblj2X+yA5iulKpQSpmBm4ANMXy/iLZUNvk9lkFmIc4d5513Hh988AGVlZUA9PT0cOTIEWbNmkVdXR07duwAXBlS7XY7a9as8aapPnLkCNXV1d4K/6abbuLXv/417e3tzJ8/H4DLL7+cBx54wFsph0tlcfnll/PII494xzNOnTpFQ0MDa9as4YUXXqC3t5fOzk5efvnlkK+fCKm5YxYUtNZ24LvAa8BB4Bmt9X6l1L1KqasBlFLLlVK1wA3Ag0qpmG1ZVJSV5PdYWgpCnDvy8vJ47LHHuPnmm1mwYAHnnXcehw4dwmw28/TTT3PHHXewcOFC1q5di9Vq5Tvf+Q4Oh4P58+dz44038thjj3nvtK+//nrWr1/PunXrvNf/6U9/is1mY8GCBcybN4+f/vSnIctx2WWX8cUvfpFVq1Yxf/58rr/+ejo7O1myZAk33ngjixYt4rrrruPCCy8M+fof/ehHzJ8/n3nz5rFmzRoWLlzo9/zPfvYzXn/9dZYsWcKrr77qTc0diSc199q1a5k1a9ZQvq3DEjeps51Oze6aNuraeqls6OIzCwqYkR/5hyHERDYRU2ef60YrNbekzo6CwaBYWpbF3MJ0spLNTMtLHesiCSHizHhIzR03QcHDkmDkgumyDacQYvSNh9TcMlFfCCGElwQFIeLYeBtTFIM725+pBAUh4pTFYqG5uVkCwwSitaa5uRmLxTLsa8TdmIIQwqW4uJja2loaGxvHuihiBFksFoqLi4f9egkKQsSphIQEKioqxroY4hwj3UdCCCG8JCgIIYTwkqAghBDCa9yluVBKNQInh/nyXKBp0LMmFvnM8UE+c3w4m89cprXOG+ykcRcUzoZSamc0uT8mEvnM8UE+c3wYjc8s3UdCCCG8JCgIIYTwireg8NDgp0w48pnjg3zm+BDzzxxXYwpCCCEii7eWghBCiAgmZFBQSl2hlDqslKpUSt0Z4vlEpdTT7ue3KaXKR7+UIyuKz/wPSqkDSqm9Sqk3lVJlY1HOkTTYZ/Y573qllFZKjfuZKtF8ZqXUOvfPer9S6n9Hu4wjLYrf7VKl1NtKqd3u3++rxqKcI0Up9YhSqkEp9UmY55VS6j/d34+9SqklI1oArfWE+gcYgWPAFMAMfAzMCTjnO8Af3V/fBDw91uUehc/8KSDZ/fW34+Ezu89LA94FtgLLxrrco/Bzng7s/n/t3VuIVWUYxvH/IxZdpBVaFCUYkZAYHTDRghKUCC/0oi6KrAwJCqyI6KqLpG6iA9JF0JkOF4UV1BCFQmciS+80oTATGwqCDpLYyXy6WMvlYKOznL3XWu7N84Nh1rfnY/b7zt573v0d+DZwWtk+o+u4W8j5GeCO8nousKvruHvM+UrgUmDbEX6+DHgPELAQ+KKf9z+MI4UFwA7bO23/DbwGrDiszwrgpfL6DWCJJLUYY79NmLPtD23vK5ubgMkfo3h8qPM4AzwEPAL82WZwDamT823Ak7Z/BbD9U8sx9ludnA1ML69PAX5oMb6+s/0J8MtRuqwAXnZhE3CqpLP6df/DWBTOBr4f0x4tbxu3j+39wB5gRivRNaNOzmOtpninMcgmzFnSJcAs2++0GViD6jzOc4A5kj6TtEnSNa1F14w6Oa8FVkoaBd4F7mwntM4c6+v9mAzj0dnjveM/fItVnT6DpHY+klYC84GrGo2oeUfNWdIUYB2wqq2AWlDncZ5KMYW0mGI0+KmkebZ/azi2ptTJ+QbgRduPS1oEvFLmfKD58DrR6P+vYRwpjAKzxrTP4f/DyaqPpKkUQ86jDdeOd3VyRtJS4H5gue2/WoqtKRPlPA2YB3wkaRfF3OvIgC82131uv237H9vfAV9TFIlBVSfn1cB6ANufAydRnBE0rGq93idrGIvCZuB8SedKOpFiIXnksD4jwC3l9XXABy5XcAbUhDmXUylPUxSEQZ9nhglytr3H9kzbs23PplhHWW57Szfh9kWd5/ZbFJsKkDSTYjppZ6tR9lednHcDSwAkXUBRFIb54+RGgJvLXUgLgT22f+zXLx+66SPb+yWtATZQ7Fx4wfZXkh4EttgeAZ6nGGLuoBghXN9dxL2rmfOjwMnA6+Wa+m7byzsLukc1cx4qNXPeAFwtaTvwL3Cf7Z+7i7o3NXO+F3hW0j0U0yirBvlNnqRXKab/ZpbrJA8AJwDYfopi3WQZsAPYB9za1/sf4L9dRET02TBOH0VExCSlKERERCVFISIiKikKERFRSVGIiIjK0G1JjegXSTOA98vmmRRbPA/uf99n+/JOAotoULakRtQgaS2w1/ZjXccS0aRMH0VMgqS95ffFkj6WtF7SN5IelnSjpC8lbZV0XtnvdElvStpcfl3RbQYR40tRiOjdRcDdwIXATcAc2wuA5zh0YucTwDrblwHXlj+LOO5kTSGid5sPnj0j6VtgY3n7VspziIClwNwxH9sxXdI027+3GmnEBFIUIno39sTZA2PaBzj0GpsCLLL9R5uBRRyrTB9FtGMjsOZgQ9LFHcYScUQpChHtuAuYX37Q+nbg9q4DihhPtqRGREQlI4WIiKikKERERCVFISIiKikKERFRSVGIiIhKikJERFRSFCIiopKiEBERlf8A97gmCp/M2i4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb0573f4b70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define the 'true' parameters\n",
    "true_parameters = [5, 3]\n",
    "\n",
    "# Run a simulation to get test data\n",
    "values = model.simulate(true_parameters, times)\n",
    "\n",
    "# Add some noise\n",
    "values += np.random.normal(0, 0.02, values.shape)\n",
    "\n",
    "# Create an object with links to the model and time series\n",
    "problem = pints.SingleOutputProblem(model, times, values)\n",
    "\n",
    "# Select a score function\n",
    "score = pints.SumOfSquaresError(problem)\n",
    "\n",
    "# Select some boundaries\n",
    "boundaries = pints.RectangularBoundaries([0.1, 0.1], [10, 10])\n",
    "\n",
    "# Select a starting point\n",
    "x0 = [1, 1]\n",
    "\n",
    "# Perform an optimization using SNES (see docs linked above). \n",
    "found_parameters, found_value = pints.optimise(score, x0, boundaries=boundaries, method=pints.SNES)\n",
    "print('Score at true solution:')\n",
    "print(score(true_parameters))\n",
    "\n",
    "print('Found solution:          True parameters:' )\n",
    "for k, x in enumerate(found_parameters):\n",
    "    print(pints.strfloat(x) + '    ' + pints.strfloat(true_parameters[k]))\n",
    "\n",
    "# Plot the results\n",
    "plt.figure()\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Concentration')\n",
    "plt.plot(times, values, alpha=0.5, label='noisy signal')\n",
    "plt.plot(times, problem.evaluate(found_parameters), label='recovered signal')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}