{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Optimisation: Using Transformed Parameters\n", "\n", "This example shows you how to run a global optimisation with a transformed parameter space.\n", "\n", "For a more elaborate example of an optimisation, see: [basic optimisation example](./optimisation-first-example.ipynb).\n", "\n", "First we will create a toy model which implements the logistic model." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from __future__ import print_function\n", "import pints\n", "import pints.toy as toy\n", "import numpy as np\n", "import matplotlib.pyplot as pl\n", "import math\n", "\n", "model = toy.LogisticModel()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The parameters for the toy logistic model are $[r, K]$, where $r$ is the rate and $K$ is the carrying capacity. We will create a wrapper pints model that assumes that the parameters are given as $(r, \\log(K)]$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class TransformedModel(pints.ForwardModel):\n", " def __init__(self, model):\n", " self._model = model\n", " def n_parameters(self):\n", " return self._model.n_parameters()\n", " def simulate(self, parameters, times):\n", " transformed_parameters = [parameters[0], math.exp(parameters[1])]\n", " return self._model.simulate(transformed_parameters,times)\n", "\n", "transformed_model = TransformedModel(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have our transformed model, we can use it in combination with all of pints' optimisation and inference algorithms. For now, we will use it to fit simulated data using CMA-ES:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Minimising error measure\n", "using Covariance Matrix Adaptation Evolution Strategy (CMA-ES)\n", "Running in sequential mode.\n", "Population size: 6\n", "Iter. Eval. Best Time m:s\n", "0 6 3.47e+07 0:00.1\n", "1 12 3.46e+07 0:00.1\n", "2 18 3.4e+07 0:00.1\n", "3 24 3.39e+07 0:00.1\n", "20 126 4077706 0:00.1\n", "40 246 358305.2 0:00.2\n", "60 366 95402.77 0:00.2\n", "80 486 95402.77 0:00.2\n", "100 606 95402.77 0:00.3\n", "120 726 95402.77 0:00.3\n", "140 846 95402.77 0:00.4\n", "160 966 95402.77 0:00.4\n", "180 1086 95402.77 0:00.4\n", "200 1206 92522.28 0:00.5\n", "220 1326 92519.12 0:00.5\n", "240 1446 92519.12 0:00.6\n", "260 1566 92519.12 0:00.6\n", "280 1686 92519.12 0:00.6\n", "300 1806 92519.12 0:00.7\n", "320 1926 92519.12 0:00.7\n", "340 2046 92519.12 0:00.7\n", "360 2166 92519.12 0:00.8\n", "380 2286 92519.12 0:00.8\n", "400 2406 92519.12 0:00.9\n", "420 2526 92519.12 0:00.9\n", "440 2646 92519.12 0:01.0\n", "460 2766 92519.12 0:01.0\n", "480 2886 92519.12 0:01.1\n", "482 2892 92519.12 0:01.1\n", "Halting: No significant change for 200 iterations.\n", "Score at true solution: \n", "93162.3109415\n", "Found solution: True parameters:\n", " 1.50473610420630947e-02 1.49999999999999994e-02\n", " 5.99762659416961963e+00 6.00000000000000000e+00\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX6wPHvmZKeUJLQEiD03lEQEEEsgChWVEDBsthdXcuiu3Z3f7hiARsWVFQUOwqoIEgR6Yj0Fnro6YS0Kef3x9xMZpJJATKZlPfzPHky99xzb87NJPedU6/SWiOEEEIUZQp0AYQQQlRNEiCEEEL4JAFCCCGETxIghBBC+CQBQgghhE8SIIQQQvgkAUIIIYRPEiCEEEL4JAFCCCGET5ZAF+BcxMTE6ISEhEAXQwghqpX169cna61jy8pXrQNEQkIC69atC3QxhBCiWlFKHShPPmliEkII4ZMECCGEED5JgBBCCOFTte6D8MVms5GUlERubm6gi1IjhYSEEB8fj9VqDXRRhBB+VuMCRFJSEpGRkSQkJKCUCnRxahStNSkpKSQlJdGiRYtAF0cI4Wc1rokpNzeX6OhoCQ5+oJQiOjpaamdC1BI1LkAAEhz8SH63QtQeNTJACCGEv63Yk0ziiVOBLoZfSYDwA6UUjzzyiHt78uTJPPvss2d8nnXr1vHggw+edTkSEhJITk4uNc9///vfsz6/EP5gcziZs/EIWutAF6VUo99fzSWvLgt0MfxKAoQfBAcH891335V5cy5L7969mTp1agWVyjcJEKKqeWtxIg98sYH5W48HuiiV6lBqNh/8vjfQxfAiAcIPLBYLEyZM4LXXXiu278CBAwwZMoSuXbsyZMgQDh48CMDXX39N586d6datGwMHDgRgyZIljBgxAqfTSZs2bTh58iQATqeT1q1bFwtAKSkpXHbZZfTo0YO77rrL6xPY1VdfTa9evejUqRPvvfceABMnTiQnJ4fu3bszZsyYEvMJUdTh9Bw+WbnfL+c+kp4DQHp2vlf6R3/sY87GI375mVXBuI/W8OK87aRk5QW6KG41bpirp+fmbGXbkcwKPWfHJlE8c2WnMvPdd999dO3alccff9wr/f777+fWW29l3LhxfPjhhzz44IPMnj2b559/nvnz5xMXF0d6errXMSaTibFjxzJz5kweeughFi5cSLdu3YiJifHK99xzzzFgwACefvpp5s2b53WD//DDD6lfvz45OTmcd955XHfddUyaNIk333yTv/76q9R80dHRZ/OrEpVg78ksjmXk0q91TNmZK9D4D9ew+0QWw7s0JiYiuFzHaK3JyrMTGVLyHBpXEDgKgKnIgIjn5mwD4MpuTc6y1CVLz85n5uqD3HNRK0ymwAzESDvtCoiOKtS0JjUIP4mKiuLWW28t1kS0cuVKRo8eDcAtt9zC8uXLAejfvz/jx4/n/fffx+FwFDvf7bffzieffAK4buK33XZbsTzLli1j7NixAFxxxRXUq1fPvW/q1Kl069aNvn37cujQIXbv3u2z3OXNJ6qGi19ZyugPVpc7v9Yau8NZ4v53luzh+w1JZZ4nzfh07zyDm9kXaw7R5dkFHEg5XWKe5+ZsI8dm/P0b9+m00/m8+uuucv8cX/48mEZSWnaJ+5/8fjMvz9/Jqn0pAKzZl8riHSfYczKL5DI+0dtK+X36kpFtIzvfzvoDqazYU9gK4DR+lXtOuH4/B1OycToDGyxqdA2iPJ/0/emhhx6iZ8+ePm/mBQqGjU6bNo3Vq1czb948unfv7vWpHqBp06Y0bNiQ3377jdWrVzNz5sxSz+dpyZIlLFy4kJUrVxIWFsagQYN8zmUobz5Rfb29ZA8vz9/J1ucuJzy4+L//S7/sAOCaHvFlnMn1d3YmH3Z/3uKqGexPyaZ5dLg7febqA/yy5Rif3tHHK//zc7bRs1k9Js/fyS9bj5X/B/lw7dsrXD970hXF9q3Yk8xPm43za7jx3ZWs3pfq3h9iNbHjhWElnvtUrp364UE+9+XbnWw5kkGPpnXd/5vdnl9AfL1QktJcTWlfTujLliOZ7ibhm99fxQtXd+ap2Vu4rX9CQO9jUoPwo/r16zNq1CimT5/uTuvXrx+zZs0CYObMmQwYMACAPXv20KdPH55//nliYmI4dOhQsfPdeeedjB07llGjRmE2m4vtHzhwoDtw/Pzzz6SlpQGQkZFBvXr1CAsLY8eOHaxatcp9jNVqxWazlZlPVD/rD6S5bzoZ2TYST5zircWJAGTl2SvkZ/j69Ky15vsNSeTku2oCK/Ykk5PvwGF8GjYX+RDzr++38Pvu4gM6svLsfPD7XvYXqXE8/s1G/khM5tNVB1i66+Q5X8PSnR7nUHgFB4BcW/Fr9OwnyMyx8dbiRBImzmPL4Qxumb7afe23TF/NtW+vYPHOE2Tl2UmYOA/AHRwAbnxvFS/M3YZnrC0o00d/7GfQy4t5ce42TmRW/oc1CRB+9sgjj3h1Jk+dOpWPPvqIrl278umnnzJlyhQAHnvsMbp06ULnzp0ZOHAg3bp1K3auq666iqysrBJrJM888wzLli2jZ8+eLFiwgGbNmgEwdOhQ7HY7Xbt25amnnqJv377uYyZMmEDXrl0ZM2ZMqflE9fLrtuNc984KPl/jGgQxdrprSGa2cePK87jp3f3pel77dZfPm/3h9JxincUABfd4m0OzP/k0Q19fxk+bj7LxUDpbj2Ty8JcbeeK7Tew9mcXo91fz3Jyt7gBhUjB30xF6v7jQfcMsSXiwpVgTz1frkhjzwWqemr2FcR+uYfxHa7jo5cVorfl2fRL7k0/jdGp+2XK0fE00HvEq9XTxay2wZl8qOfkOnE5NrxcXutOvmPo7L8/fCcCIN5bz++5kVu1L4a5P17mDzd6Tp5n++75Si3EqtzBon/YI4PtTsvlg+T7O/+8iMrJtTF20m4SJ89x9Fv5Uo5uYAiUrK8v9umHDhmRnF7Z9JiQk8NtvvxU75rvvviuWNmjQIAYNGuTe3rhxI926daN9+/Y+f250dDQLFixwb3uOovr55599HvPSSy/x0ksvlZlPVLzXft3FkA4N6Bpf95zPNX/rMfq1inZ3ABe08yeecP0tbj6c4ZU/z17Yz/XL1mP8svUYY/s2d6clTJzHx7edx/iP1gIw7oLmPDeyc7Gfu/dkFqv2prDj2CnunfknAD/c1x+A2X8dob/ReT5r7SF6NXf1iSWezOLpH7YWO5evvpEcm8PrxunLEuPT9tRFiby20NVX8d9ruvDk95t5YWQnbrkgocRjNyWl8+7SwqGlxzJ8f0rv+99FHCvhE/zp/OJ9hi/9vIMdxwon0U1ZuJtTZ1BrW7k3xWf68Km/c9gY5bU3+TS9Smjaqih+r0EopcxKqQ1KqbnGdgul1Gql1G6l1JdKqSAjPdjYTjT2J/i7bNXJpEmTuO666/i///u/QBdFVACtNVMW7eaqN/84q+NP5do4lWtzb9/16Xq6PLuAR77ayNDXCydvKRSjpq0sdnye3XUz9hyqurvIrOCC4AAwY6X3A8gKPnTfMWMd7xf5ZDzyrcJreuybTe7XBc0uvoIDQL6PAPH56oPuspZl1tqD7tdPfr8ZgKMZuV4jGRMmzuNIeg7Ldp3kUGo2172zwuscL87b7vPcJQWHkngGB+CMgkNpCoIDnHnn+NmojBrE34HtQJSx/RLwmtZ6llJqGnAH8I7xPU1r3VopdZOR78ZKKF+1MHHiRCZOnBjoYogKYnOU3vSRdjqfPLuTRnVC3Glfrj1IZIiV4V0ac/5/FhWO9vHw7Z+uEUgFN7rMXBtr9qcWy5drc7C3yCf50e+XPhpqx7FMZqzYT9+W0Zw4deZj9bcdLX3IeZ6Ptv4zcTQjBxMahcaExoSTr1bu5qMlWwkDTDgxoVmwfgev/7qTqBAzUQ4bJjS4jzGOV55pheVS7u/aI837vTQrhUNrd7pnj4uv485kPx7789JbAv4dgu7XAKGUigeuAP4D/EO5uvEvBkYbWWYAz+IKECON1wDfAG8qpZSu6vPthTgLvj4te7ro5cVk5tq9Rt3881vXp+L9k67wGRx8+Wa97yGruTYn93z2ZzlL6zL09d8B13BVE07CySWCHMJVDhHkEqFyCMf1OljZCCafYGyuL2V8L0hzb9uwYseiHIR8+iqzg9Kw4MCMAwtO47sDs3Jidac7MON0fzfjNG7oJQgpsv07jA8pYV818ue+POjZwa8/w981iNeBx4FIYzsaSNdaF9S3koA443UccAhAa21XSmUY+c9tvQohqqD8MppNMo12d621X1bQzcixsfN4QTOIJopsGqlUGqlU6nOK+uoU9dQp6mN8V6eoTyZ11WnCySFcnVkNIk9byMNqfAWRp63ku7etOLSJtYfzcBCBHTMOzNgxub/bnRYcmIx9JnceBwpHwed8rdAonMYXxvfS01yt7E5MXukATl2Y5hl+tEc9ojCtUGRoEJk5NnfeqFALGTl2H8cX3z/6/GY0iw5n0s87vM4bXy+cQ8Y8joLjxzUbdEbvwdnwW4BQSo0ATmit1yulBhUk+8iqy7HP87wTgAmAe5SOENWNZ4C457P1vDO2l898p/LsRJUy87i8IskmQR2jhTpGgjpG25Vf8Jl1L42NoODrhm/XJtKIJFVHkkYkiTqODGc4WYSSpUPJIoTTnq91KNf168DbK46Rq4M8AoIVfYbdnW0aRDDngQE88d1mvt9w+JyvvyL8cF9/r/6VksSoIJKdhSOMhiY0KnUeR8eoKEJjzKw/kMYNbXvTomNDzAfXF87NAIbENmBRygmv4z44z/ffTEXyZw2iP3CVUmo4ropcFK4aRV2llMWoRcQDBYurJAFNgSSllAWoAxRrPNVavwe8B9C7d29pfhLVkmeA+HlLyTePYxm5fLryAEt2Ft4cZq05WGL+YPJpq5LoYDpAR3WADqaDtFRHiFXe7f9Hj9TnlKrPDt2Uxc7uHNX1OaajOabrkUIUqTqSU4R53dgjgi0lzp+oF2bl2as60bpBBAf/WF5i+ZpHh3EgpeQZzQXaNIwgxGqmY+Mod4CwmlWxvpvezeux7kBameerCN2aFh9tFh0eREqR4aZFG8VbxhZOChzSvgGLdnjf6GMjg8nOd/1erRbX7/vVUd35afMv7jyDPY5rEBnMs1dVzuQ5vwUIrfUTwBMARg3iUa31GKXU18D1wCxgHPCDcciPxvZKY/9v1bX/wWw206VLF/f27NmzSU5O5pNPPmHq1KksWbKEoKAg+vXrF8BSikDKL7Kcyrr9qfRqXo88u5NnPDqOH/h8g0dTkMvE7zYbrzRxJNPbtJPepl30Nu2ijUrColzBJ0uHsEM3Y6GjF/t1I/brRuzTjTigG5LHmQ+PLO3fsVOTOozsHuc1ssqXemFB5QoQzxtDaj37alxNbZogs8md/uxVnRjxRskB6UwoBRueupTuz/9a7mN8rZtUdPmRVrER7tetG0YUCxDDuzTi2/VGEDTWgQqxmomrG+oetdQoKoRnr+zIs3O2MaxzI4Z3aVzuMp6LQMyD+CcwSyn1IrABKJhmPB34VCmViKvmcFMAylYhQkNDiy2VkZCQQO/evQHXkhYRERESIGqQb9cnEWI1c0XX8v3jFh26ef20lcREBNGoTghbDhd+2i8aHKLIYqBpM4PMG+ln2kIT5apkn9KhbHC2ZqGzJ1udCWzXzTioG5xx045SpS+f8coN3fhx45FiM5jTc1yfoiNDrGx46lI2H87g1g/XFDs+yOJdHrNJ8dboHtzt0WF+8/nN3AsAena/mJVi1RNDCA828+7SvdzQO54Qa/EVBcryx8SLCbWa6fmCdyCY/9BA6oYF8cyVHZm/9RjX9IhzDwwY3y/B57kcPkajFZ2bd02POBxOzcgeTXhvqfdy3h/fdh4XtY3l2z9dAcKzv8lzIEKe3cm1veJZeyCNB4a0Kfe1nqtKCRBa6yXAEuP1XuB8H3lygRsqozyBsGTJEiZPnsybb77JtGnTMJvNfPbZZ7zxxhtceOGFgS6eOEePfL0RgCu6Fl/rp0DCxHncfVErJg5r77OTOjkrn+Ss4rNj49VJRphWcon5T3qo3ZiVJk1H8IezM9Oc7VjnbMcO3czd4XouFjw0kC1HMnj4y43utIcvactrC3dhUorresVzXa/4YjOgj6QXzhOoFx7EwLaxPs9/90UtWbMvlfMT6rNmfyqjz2/G0M6NWfuvSzjvP67ZybkeN8ZxFySw4+gpftx4BJPCPez30cvbAWc+F+CT288nrm6oz31tG7rG0tzWvwW39W9Brs3B6n2p/HNoexpG+R7uVLS2cGnHhqwqMsnNZFKMOq8pAGEe61+9cXMPBrVr4JVXe3S7FjQ7gWtxxKgQK2+N7lnWJVaomj2T+ueJcGxz2fnORKMuMGxSqVkKnrEA0KJFC77//nv3voSEBO6++24iIiJ49NFHK7ZsosqbtnRPiQHCUwwZXGlewZXmlfQ0udZP2uRswZuOq9kZ0Ydf0uN9BoRfHrrQPRz1vVt6MeHT9SX+jGt7xLExKR2LyeSuqbRuEEGdUO9O8QeHtCY9J59rPRbwmzisvXukjdWseM5Hm3jrBhEknsjiroEteXfZXga3i+Xi9g3ZP+kKpi3dw5r9qYQGuWoAsZGFS4bneMxMDg+28PzITkaAKD6OxWo28dfTl7J010kOpGTTOS6KWWsOsWDbcUzK+9N8ZIjFZ+C6c0ALvvDRrxNiNfPqqO4+f3eF57R6zaR+d2wv5mw6wt9n/eUzf7hxvTf0ivdattzXCJ2pN/Xgka82YndqBrXzHXD9rWYHiADx1cQkap/XF+5i2a6TfHdv/2IBwfc8CM0Fpm2MMS/ictNarMrBVmdzJtluYq6zL0na9Wnzo2vO4yePWc6e6oUV9i2U9twFgFdGdUMphdaa//t5B38kJqOUwmL2DjxKqWIrit59UStu6duclKx8mkWH+Tz/j/f3J9fmZPPhDN5dtpeWHm3xBZPigi3Fg1yu3bt/pqDZxTOIeKobFsTI7nHu7YKlM94c3dO9/AfAiK7ez5GY9+AAcm1OejWvx79HdPR57qKm3NSd9QfSWLT9BIfTcxjauREfr9jv3m8yKUZ2j3MHiCWPDvI6vqAGkW0reo3GC4+AdlmnRmx+rlG5yuUvNTtAlPFJXwh/en1h4bM0PNuTv1p7iMe/LVyCIggb15p/507zT7Q2HSFdhzPDcRmzHINJ1N7Lbl/SoQGRxk0mPMjM9b3i3ctg1A8P8mqzjwwp/Pd+Z0xPOsfVYebqg8xcdYCWDSLcN16lFE8OL5xwZTGXb95FeLDF55LhBcKCLIQFwcA2MUy+oRtXdivsnynopA/yCEaPD23H/37ZyXU9va+5TqiV/1zTuVhzTElOGrO8WzeIoFVsOHtOutal6hJXxytfpyZ1ih1blpHd4xjZPY7IkB28tdhVGxzXL4G5G4/QxEfTVUJMuNf2ZR0bcmW3JvxzaDuv9I6N67Bqbyr1/Ly20pmq2QGiioqMjCQzs2KfdCeqNs9mk4LgEEouo82/8TfLPBqpNDY5W/BI/t3MdfYtNsqoWf0w0rPzuX1AC3fHbJO6oTw3sjPPjexMTr4Dpbw/kUcYN+8besUzzBj1MnFYex6/3PvmVJTVVLFLtCmluL6X902/oEbl2Wl976DW3Duotc9zjOnT3Ge6L+P7J/D0D1tpWi+MuHph7Dl5munjenNx+/IFmPJ45NJ2PHBxG0KsZlrEhBfrOG4RE86+5OIPRgqxmnnj5h7F0icOa8+wLo3o0Diq2L5AkgARAFdeeSXXX389P/zwg3RS1wIv/bLD61OxGQejzEt42PItDVQ661VnHs27m+XOzvhujXbdcGbcPhiARGNRPYdHA3tBWz7gHh5pNik2P3sZYUHe/+ZlPVLT6lGDWP7PweW7yDNUECB8NTGdq1svSOBWYwXXKTd2Z9GOEwzp0LBCf4bJpAgxlTyC6rt7+nG0hJVhfQmymDgvoX5FFK1CSYDwA8/lvgt4Lt3dtm1bNm3aVCyPqP7y7A6CLd43jneW7HGPzBls2sCTls9pYzrMWmdb7s1/kNSYXuzNKfkxnFBYGwCwGJ/w7WU860CpsvshfDF7BJD4er77F86VzSh78FkMUz0T9cKDitVeKkO98KAq11x0NiRACFGB2v37F67pEcffizQ5/PTHn0yzzmCoeS17nY24K/9h5jt7A4pW5ThveHDhjbTgBu4oIUAU7D/baaZKKR67vB0XlTBUtSL849K22B1OrukRV3ZmETASIIQ4R8uLPC7z+w2H+SPRlaZwMsa8iH9aZmHFziTbTXzgGI7d419P41qq4ubzm/H2kj0+f0ZEcGFNoEFUMPXCrPzrCt8reX44vjez1hwivp7v8f7lcd9g330BFSUmIpj/XV/8qYmiaqmRAcJfK2CK0pdbqK3GTi/+HIUTp/KIJY1XrNMYaN7MMkcX/m2/nYO6oftpZ24aNjx9GeCaEPXFmkO0bRjBZ3f24XSeg9cX7uKBiwtv2MEWszu/L60bRJZ72KYQpalxASIkJISUlBSio6MlSFQwrTUpKSmEhFTjRfQr2PvL9vpMH2Jaz/+s7xFGHv+y3c5MxxAKOqDDgswlLmnx4tVduHdQa2Iigl0dz5Ew5abio16EqAw1LkDEx8eTlJTEyZMny84szlhISAjx8ZXf6VdV/ecn70dUmnDyiOUr7rP8yFZncx603c8e7d3ObjWbWPbYYJYnJvPEd5u91rQ3mxRN6/unY1iIM1XjAoTVaqVFixaBLoaohaI4zRTrmww2b+Rz+2CetY8nn+KjiKxmVxDo08I1rFGa7URVVfGDkIWowZxOzeIdJziRmevVj9BMHef7oKcZYNrCk7Y7eNL+N6/gMPPOPrRv5FoMzmrMHi5Y82hAm5hKvAIhyq/G1SCE8Kd1B9K47eO1dImrw+bDGQB0Uvv5OGgSZpyMzv8Xa3V7r2OeubIj/VvHuGcNFwSI6Ihglj02mMZ1pU9HVE0SIIQ4AxnGs4YLgsMFpq28Z32VTMK5Nf+fxfobAIZ2di24VhAYPNc6KmmhOyGqAmliEqIcMnNtzN5w2GvRvYtNf/Kx9SUO6xg+aPsud107jCHtG/D6jYVLRO/+zzAa13HNRyhYmM5ZxgxoIaoKqUEIUYqtRzKICrHyj6/+Yu3+NP52oWsAxCDTBt6xvs4O3Yxb8idyeVAso85r6n4wzIaDaRxMzXbXGqBwYbq8M3zIjRCBIgFCiFJcMdX7eceZOXYuMm3kXetr7NLx3JI/kUwiKDrl5jnjmcqeCgJEWQ8LEqKqkAAhxBnYv34BnwS9SqKOY2z+k2TiegiOKmEVVk8Th7UnM8fGgNYyaklUDxIghPBh8Y4TtGkY4ZXWVh3i/aBXOKgbMCb/STIo3F/WEtoArWIj+PKuCyq8rEL4iwQIIYr4et0hHvvGezn2RqTwcdBL5BDE+PzHSScSs0m5V1QNC/LvstVCBIKMYhKiiKLBIZJsPgr6H5HkcFv+4xzGtQy253MTwiVAiBpIAoQQpVA4ec36Fq3VEe6xPcQ2neDeN7xz4QPlQ4OkMi5qHgkQQnhYust7kceHLN9xiXkD+3v/i+XOLl77ejavx9i+zQBpYhI1kwQIITzc8fFa9+vLTGv5u+U76D6GlsMfLpY3rm4o9wxqTZ8W9RnZvUllFlOISiH1YlGr7TmZxewNh+nRrC63f7zOnd5KHeYV6zR2mNrQ/opXMXtMeFv5xMUs2Hqci9s3QCklI5NEjSUBQtRqt320loOp2V7PXw4mnzetb5CPhadCJvK11bWY3qRruxARYqFxnVDG9UsIUImFqDwSIEStVjCr2eqxgN4Tls/pYDrI+PzHOE60O/2m85tVevmECCTpgxC1WsFI1YLH015iWs94ywI+sA9jibOHe56DELWR1CBErZUwcZ7XdkNSedn6LlucCfzPfhOABAhRq0mAELVS0cd8/rrtGDOs7xGMjQdsD3DnoPas2pvCPy5tF6ASChF4EiBErZRfZMntUeYlXGTexNO2cezTjXl8aPsSjhSi9pA+CFEr5Xksud2YFP5t+YyVjo586rg0gKUSomqRACFqpcJnMmgmWd/HjJPH7X9DY6JLXJ2Alk2IqkKamEStVFCDuMG8lIvMm3jKNp5DuiFPDGvPNT2LP1daiNrIbzUIpVSIUmqNUmqjUmqrUuo5I72FUmq1Umq3UupLpVSQkR5sbCca+xP8VTYhRk1bST0yedLyOWuc7fjMcQkAd13UigaRIQEunRBVgz+bmPKAi7XW3YDuwFClVF/gJeA1rXUbIA24w8h/B5CmtW4NvGbkE+KcZWTb6PPfhfx5MM2ddjg9h4mWWUSQw79tt6OltVWIYvz2X6FdsoxNq/GlgYuBb4z0GcDVxuuRxjbG/iFKFX3SrxBn7q+kdI5n5nHt2ytImDgPrTW91E5utCxhumM4u3TTQBdRiCrJrx+blFJmpdRfwAngV2APkK61thtZkoCCBt844BCAsT8DPNY5KDznBKXUOqXUupMnTxbdLUQxRR/mk5eXx3+sH5KkY5hivwaAK7s14ffHBweieEJUWX4NEFprh9a6OxAPnA908JXN+O6rtlBsGqvW+j2tdW+tde/Y2Fgfhwjhrehk6Pw/3qa96RDP2saRQwid46J44+YeNK0fFpgCClFFVUrDq9Y6HVgC9AXqKqUKRk/FA0eM10lAUwBjfx0gtTLKJ2q2PLvD/TqWdMJXTuZXR08WOnsBMOnaroEqmhBVmj9HMcUqpeoar0OBS4DtwGLgeiPbOOAH4/WPxjbG/t900fUQhDgLebbCSXH/sHyNtufxon2sOy3IIh3UQvjiz3kQjYEZSikzrkD0ldZ6rlJqGzBLKfUisAGYbuSfDnyqlErEVXO4yY9lE7VA6ul86ocHuec8dFT7udG8hE/0cA7owudJB5klQAjhi98ChNZ6E9DDR/peXP0RRdNzgRv8VR5Ru2w5nMGIN5Yz5abuxoqsmn9bPiOdcF7Nu9orb3y90MAUUogqTj46iRrpUGo2AHM3HSXP7uRS03r6mbfxmv16Mgl351v1xBAsUoMQwidZakPUSCHG0NbkrDy+WJHIFMvCias+AAAgAElEQVRMdjvj+NwxxCtfozoya1qIkkiAEDVSbr5r5NKGg+ncZv6ZFtbjjM9/HAeuwDH5hm5c2a1xIIsoRJUndWtRI+XYXAEinBzut8xmuaMTS5zd3fujQiwEW8wlHS6EQAKEqKGyjRrEHeafiVaneNl+o9f+sCCpPAtRFgkQokbKtTmoRyZ/s8zjF8d5bNStvfaHBkntQYiySIAQNdLmwxncbZlDGLlMthcfPS1zH4Qom/yXiBrneGYuq/7awjjzAr53Xkiijsds8l7qyyR/+UKUSf5NRI2z6/gp/m75DhNOXrdfB0DTeqEseHggY/o0AyA2MjiQRRSiWpCeOlFjnDiVywOfb2BwbCZ3mpfwvWUYSXmuFX810LZhJM+P7Mzfh7SRp8YJUQ5SgxA1xmcrD7B6XyqxG6ZiVxac/R927ytY9tFsUjSIkuAgRHlIgBA1hs2paa6OcbXpD77ickYN7s20sb0CXSwhqi0JEKLGcDg195tnY8PCDHUlSim6xNcBoGGU9DkIcaakD0LUGDnHd3ONeTkzHJeTbKoHQFzdUF6+viuD2zcIcOmEqH4kQIhq7dNVB/hjdzL/uKwtnfZMx2E2M80+Am0pfNbUDb2bBrCEQlRfEiBEtfbU7C0AjG6nuc78O585LuEk9YgMcLmEqAmkD0LUCEErX8eJYpr9SqBw1JIQ4uxJDUJUe01IpmfqT8xyXMwH91/F9OV7GdO3eaCLJUS1JwFCVHv3WH4EYH/7CdwaX4fXbyr2pFshxFmQJiZRrTUgjVHmJXzjuIiWrdsFujhC1CgSIES1drvlZyw4eMdxJZEhUiEWoiLJf5SodpLSsknPttG5vmaseRFznRdwSDckXB4CJESFkv8oUa2knc5nwEuLAfhr8GbqqhzetY8AIEJqEEJUKPmPEtXKte+sACCYfGwr3mKJsxvbdAIAEcHy5yxERZI+CFGt7Es+DcAN5qXEqkzesV/l3hcuAUKICiUBQlQ7ZhxMMM9lg7M1q3V7d3rjOrKMtxAVqcyPXEqphsB/gSZa62FKqY7ABVrr6X4vnRA+DDetppnpJC/mjwUUj17WlhCrmRCrOdBFE6JGKU8N4mNgPtDE2N4FPOSvAglRElfzkuYeyxwSnU341el61sP9F7fhzgtbBrZwQtRA5QkQMVrrrwAngNbaDjj8WiohfBg8eQkXmTbR0XSAdx0j0Jj46q4LAl0sIWqs8gSI00qpaFyP9UUp1RfI8GuphCjiVK4NgLvNcziq6zPbMQCAsCBpVhLCX8oz7OMfwI9AK6XUH0AscL1fSyVEEYfTc+ihdnOBeRsv2MZiM/50gy0yzkIIfykzQGit/1RKXQS0AxSwU2tt83vJhPBwLCOXuy1zSNfhfOG42J0eJAFCCL8pzyimW4sk9VRKobX+xE9lEsLLwZRsjiZu5GbzOqbYryEbGc4qRGUoTxPTeR6vQ4AhwJ+ABAhRKQa+vJj/Wd4lxxzEDPvlXvtiIoIDVCohar7yNDE94LmtlKoDfFrWcUqppriCSCNcI6De01pPUUrVB74EEoD9wCitdZpSSgFTgOFANjBea/3nGV2NqFFybQ6CzCYakcLV5uV86byEuROvJjYyGLNSmEwq0EUUokY7m7UJsoE25chnBx4x+jAigfVKqV+B8cAirfUkpdREYCLwT2CYcd42QB/gHeO7qGX+9f1mZq4+CMDfLmzBHZafMaF51z6cW+qGBrh0QtQe5emDmIMxxBXXsNiOwFdlHae1PgocNV6fUkptB+KAkcAgI9sMYAmuADES+ERrrYFVSqm6SqnGxnlELVIQHAC++n0zK4IX8aOzH+3bdwpgqYSofcpTg5js8doOHNBaJ53JD1FKJQA9gNVAw4Kbvtb6qFKqgZEtDjjkcViSkSYBoha71byAcJXHu/YR/DC6Z6CLI0StUp4+iKXn8gOUUhHAt8BDWutMV1eD76y+fryP800AJgA0a9bsXIomqiBXBdIlhDzGW+azyNGDjMi2staSEJWsxEHkSqlTSqlMH1+nlFKZ5Tm5UsqKKzjM1Fp/ZyQfV0o1NvY3Bk4Y6UlAU4/D44EjRc+ptX5Pa91ba907Nja2PMUQ1chfh9Ldr0eZlxCtTjHNfiVvju4RwFIJUTuVGCC01pFa6ygfX5Fa66iyTmyMSpoObNdav+qx60dgnPF6HPCDR/qtyqUvkCH9D7WXBTsTLPNY52zLWt2e3gn1A10kIWqdco9iMvoK3DOUtNYHS8kO0B+4BdislPrLSHsSmAR8pZS6AzgI3GDs+wnXENdEXCOlbitv2UTNkWd3AnCFaRXxKplnbOPKOEII4S/lGcV0FfAKruW+TwDNge1AqUNKtNbL8d2vAK7JdkXza+C+ssojajZXgHAt6b3LGcdvTmlaEiJQyrOQzQtAX2CX1roFrpv7H34tlaiVEk9k8eeBNAab/qK96RDT7FeiMXF9r/hAF02IWqk8TUw2rXWKUsqklDJprRcrpV7ye8lErXPJq64Bc18GzeEYMfzo7AfA5Bu6BbJYQtRa5QkQ6cZQ1d+BmUqpE7jmQwhR4XqqXfQx7eCNoDuw557NRH8hREUpbZjrm0qp/rhmOGfjeszoL8Ae4MrKKZ6obe6xzCFNR7AwdCgAH40/r4wjhBD+UtpHtN24ZlE3xrW43hda6xmVUipR62w5nEEblcSl5vW8br8WgsKBdJy62FxJIUQlKW0exBSt9QXARUAq8JFSartS6imlVNtKK6GoFUa8sZy7LXPI1sHMsF9GnxaueQ8No+TZD0IESpmjmLTWB7TWL2mtewCjgWtxDXMVosI0IZmrTCuY5RhMGlE8cllb5j4wgM5xdQJdNCFqrfLMg7ACQ4GbcA1xXQo85+dyiVrgeGYudqdmRWIyd1p+AuAD+3Cmje1JsMUswUGIACsxQCilLgVuBq4A1gCzgAla69OVVDZRw/X57yIA6nKKFcGL+dHZjyPEEBspT4kToioorQbxJPA58KjWOrWSyiNqofGW+YSpPKbZXYPjEqLDA1wiIQSUEiC01oMrsyCidnE4XaOTIsjmNvMvLHD0YreOJ8RqIlqeMy1ElVCepTaEqHA5NgfgeiBQHZXNVPs1ADidgSyVEMKTBAgRELd/vJYwcrnT8hOLHd3YolsC4JB5D0JUGRIgRECs2ZfKGPNC6qss3jBqDyFWE89e2THAJRNCFJDFbkRABJPPBMs8ljs68ad2zbvc8cKwAJdKCOFJahDC71Ky8pi3qfDhgPM2HeVm82/EqgzesF8bwJIJIUojNQjhdxM+Xc/6A2n0bXkJGTk2Hv58NUuD57La2Z7VukOgiyeEKIHUIITfHUhxza10aE12voMbzEtprFLdfQ9CiKpJahDC74wpD6SezufRWWv5wPIjG5ytWe7sHNiCCSFKJTUI4XcFS3YPff13eqb+RLxKZor9Gkp+ZLkQoiqQACH8zmlUIYLJ5wHL96xztmWJs3uASyWEKIsECOF3BU1MY8yLaKxSecV+AwW1hzF9mgWuYEKIUkkfhPCbU7k2Zm84TFaenTByudfyA8sdnVjp7ATA/klXABBsMdM8OiyQRRVC+CABQvjNde+sYNfxLADGm38hRmXyin1UsXxPy+xpIaokaWISflMQHKI4zV2WuSx09GCDbgPAm6N7BLJoQohykBqE8IucfIf79Z2WedRR2bxqv4GdLw4l2GIOYMmEEOUlNQjhFz9vcS2tEU0Gt5t/Ya6jD9t0ggQHIaoRCRDCLyxm15/WQ5ZvCSGfV+03YJJpD0JUKxIgxFn7ddtxVu5J8bkvyGyilTrMzebfmOkYwl7dhPAgadEUojqR/1hx1v72yTqgcLiqJ6fWTLTMIodgptivA6BOmLVSyyeEODdSgxDnLCPbxolTuazZl4rd4XpmaPiRlVxqXs/b9pGkEgVAXQkQQlQrUoMQ5+yCSYvINkYtXdKhIdf3bEKbtf/lsI7mQ8dQruzWhDkbj3Bhm9gAl1QIcSYkQIhzlu0xpHXh9uNE7PyWoUG7edh2D2ufHUFUiJW7BrakQ+OoAJZSCHGmpIlJVKgIsnnS+jl/OVsy29mfEGNYa+e4OphlGJMQ1YrUIESFetDyPTFkcKftEdo1qkOQRT6DCFFd+e2/Vyn1oVLqhFJqi0dafaXUr0qp3cb3eka6UkpNVUolKqU2KaV6+qtcwn9aqyRuM//Cl45BbNKtaNswMtBFEkKcA39+vPsYGFokbSKwSGvdBlhkbAMMA9oYXxOAd/xYLuEXmucsM8gmmJftNwIQapVZ00JUZ34LEFrrZUBqkeSRwAzj9Qzgao/0T7TLKqCuUqqxv8omKt5w02r6m7fysv1G97BWIUT1VtkNxA211kcBjO8NjPQ44JBHviQjrRil1ASl1Dql1LqTJ0/6tbCifCLJ5mnrp2x1NudzxxB3emauLYClEkKcq6rSg+hreIv2lVFr/Z7WurfWundsrIyrr2xpp/PZcjjDK+0Jy+fEks5E299wevxJhQZJE5MQ1Vllj2I6rpRqrLU+ajQhnTDSk4CmHvnigSOVXDZRDle//QcHUrKJiQgCoI/azmjLb7xnv4LNuiUA+/5vON+sT+LSjg0DWVQhxDmq7BrEj8A44/U44AeP9FuN0Ux9gYyCpihRtRxIyQYgOSufYPKZZH2PA84GvGq/3p1HKcUNvZtSNywoUMUUQlQAfw5z/QJYCbRTSiUppe4AJgGXKqV2A5ca2wA/AXuBROB94F5/lUucmT0ns3A4fbb28ZDlW1qYjvOE/U5yCQZg0SMXVWbxhBB+5LcmJq31zSXsGlI0QWutgfv8VRZxdpLSshnyylLG90vg4Uvaeq3G2kvtZIJ5LrPsg1jh7AzARW1jSYgOD1RxhRAVTGZSixIVPDb04xX7+XjFfu4f3BpwLafxmvVtDusYXrDfArjmPMy4/fyAlVUIUfEkQIgSObR309KbixMBeNryKXEqmVH5T3OaUADy7I5ixwshqreqMsxVVEF5NmextMtNaxhlWcrbjpGs1+3c6SV0UwghqjEJEKJE+Q7vABGvTvKS9X02OVswxX6t176WMdL3IERNI01MwqcVicnMXH3QvR2EjbesUzDh5AHbA9iNP50L28Rw76DWtG4QEaiiCiH8RAKE8Gn0B6u9tp+yfEo3014m5D/MAd3InZ5vd3JBq+jKLp4QohJIgBBeft12nFcW7PRKG2lazi2WhUyzj2CB8zyvfXF1QyuzeEKISiQBQniZ+O0mUk7nu7e7qL1Msn7Aamd7JttHudNvPr8Zg9vF0r91TCCKKYSoBNJJLbyEBxd+ZmhMCtODJpNCFFv6TXH3OwBYzYrLOjXyyi+EqFkkQAjSs/MZ/f4qjqTnuG/4YeQyPWgyoeRxe/5jNG2a4HWMxSR/OkLUdPJfLvjuz8Os2JPCu0v3YFJgxsFU6xu0Uwe53/Yg11x+KQPbxhLhUVuwWnyt0C6EqEkkQAj3gzdSs21sO5LO/6zvcol5A8/Yx7PU2Y3b+icQYjUzfVxv92NEezStF7gCCyEqhTQg1yLvLt3Dst0nmXlnX690bSypkZqVy3OWGVxnXs7LtlF85rgUgGCL63NEn5bRbH9hKIdSs2laP6xyCy+EqHRSg6hFdhw7xeq9qWTk2Ljv8z9JSsvmm/VJvDhvO6Dpd+BtbrX8yjT7CN5yjHQfp5R3c5IEByFqB6lB1CJ5dgd2p+a9ZXuYt+koUSEWvlhzCNA8YfmcuyzzmGkfwiT7zfh+CqwQojaRGkQtcTwzl4XbXE943ZTkeqb0F2sOoXDyouVD7rLM42P7ZfzbfhsXtJS5DUIICRC1xnXvrHAvvrfxUDoAVuxMtk5jrGURb9uv4ln7OIZ1aUKPZnUDWVQhRBUhAaIGe2HuNjo/Mx+ApLQcd3pmrp0osphhneTukM7o/ySgeGJYB3o0kxFKQgjpg6gxNidl8NGKfUy+vhsmk6v/YPryfQC8ZTzop0BzdYzPQibTwHmCh/Pv4bfgi9k4rAMPX9KWEKuZpvXD+P3xwSgFzuKPhBBC1BISIGqIO2as5cSpPB6/vD2N6oSQkpXn3vfy/MLF9waZNvCa9R2CzWbG5j7JWt2eaCOghBhzHEBGKgkhJEBUe28tTqRlTDgO45FuBY8JHf/RWq98Zhw8avmKeyxz2O5sxpHB77N2TgoAFrOMWBJCFCd9ENXMtiOZJEycx1+H0snOt/Py/J3cM/NP9wqsOfmuZ0NvPpzhPiZeneCLoBe5xzKHz+0Xc3X+80TFtXXvnz7OewlvIYQAqUFUGVprrwlp329IolOTOrRtGOmV77cdxwH4ddsxLKbGxc6Ta3O4XyucjDUvZKLlC5yY+Hv+vfzgHABAvTCrO1/nuDoVei1CiJpBAoSf2RxOzEq5O45LctN7q2hcJ4TXb+qB06l5+MuNBJlN7PrPMK98BY+JNinFiDeWFztPjs3Bv2dvppU6zIuWj7jAvI1lzm5MzL+DI8Tw4/39ScnKp3WDSCKCLWTl2SvsWoUQNYsECD969setfLxiP9f2iOPVG7uXmnf1vlQAZv91hNjIYAD3vAVPBX0MJuU74MxatplOu97m2aAFZBPC47a/sSJyGEfycwHoGl84x2H5PweT41HjEEIIT9IH4Ucfr9gPwHcbDgPwyFcbSZg4r8zjTp5yjUAKC3KNKvp581HunLGO/cmncRqd0eYiNZJg8hlnns+Te8Yw3jyfLx2DGZT3Kl85BhNWwkN96oYF0biOPDJUCOGb1CAq0bd/Jp1R/naNXP0P98z8E4CF248zuk8zoDBABGHjRvNi7rX8SGOVyipnB16w3cJWneA+T6jH8FUhhCgvCRClmDx/J3XDrNx5YcsS8zz69UZaN4jg7otaeaUXDDsty+k8O9e9s8LnvujwoGJpn68+CECoLY37zLO5xfIrjVQaa5zt+IftHlY6O1J0ob0Qq5mf/36h1wN/hBCiLHLH8EFrTZ7dyZvGDOQ7L2zJo19vZFjnRgzp0NAr7zfrXbWCogHidH7pnb+H03P4aPk+PjBmO/uycPsJjmbkeKRouqk9jDEv4tqVK7FY81nm6MKjjrtZ7uxMSSuwRgRb6NA4qtTyCCFEURIgPBxOzyGubihPfr+FL9YcdKfvPHaKb9Yn8c36JPZPuqLM82Tk2Jhh9D/44nBqRr65nOSs/DLPtW5/GvHqJFeblnONeTmtTEfJ1sGsqTecp472Z4+OK/McDeuElJlHCCGKqrWd1Ccycxnxxu/sOJaJw6mZs/EI/Sf9xorEZK/gADDijd+LHW93OLH5GGUE8OR3m3n1111eaduPZrpf59udZQQHTXt1kPvMs+m76AaWB/+dR61fc5K6PGabwKXqXR7MvMUdHNY8OQSASzo0dI+A8tSjqazOKoQ4c7WyBvHxH/t4ds42AIa+/jtj+jTDanbFyh3HThXLb3MU9ieknc7n4a/+Yv3+NK9mJLvDicU4x4HU08XOMWxKYZBJ9lgnqUAs6fQxbaevaRuDzBuJV8kAbEhvzUeOG/nB0Y/DxALQICyY5FOF54iOCGbVE0OIiQgi3+Fk7f40xn24hg6No/j3FR3o1yq6/L8cIYQw1MoAUXTS2szVB7mlb3MAVu1NKfXYYVN+51hmbrH0i15ewuH0HK7tEceWw5k+jix08f9+pZNKorNpH91UIn1MO2hlOgrAKR3KKmdHpjqvYbGjByepy03nNaWvQ7tHQZ32mNw26doumE2KRkYzksVsolVsOAA9m9Wlf2t5+I8Q4uzUygARFWItlnY0w3XTX7DteKnH+goO4Oq/gMI5DwAmnMSpk7RSR2iljtJKHaazaT/t1CGClesmn6nDWONsxyzbYFY5O7JNN8eBmQkDW3Jy2V4AJl3XFSgcJns6v3ByW+sGEcXKEl8vjLkPDKBNw+L7hBCivKpUgFBKDQWmAGbgA631JH/8nDqhxQPEwu2lBwZfFE5iyaCJSqGJSqaxSiFOpdBEpZCgjtFCHSNY2dz503QEW53N+cg5jC3OBOq3Pp9Pdym0j66gWy9oTr7dydU9CjuhPxp/HiaT4s3fdrN2fxoAYUG+30JZX0kIca6qTIBQSpmBt4BLgSRgrVLqR631tor+WVFFAoTCSRh5RJBDhMoh0vgeQQ711Snqk0m0yvR4fYpolUF9TmFV3ktVnNbBHNXR7NONWOrsSmpoAuuyYtirG5OG91DT8TEJ6F373dt/u7AF3284THJWPsEWM89e1ckr/+D2DQDoHl+Xvv+3iBybo9iMaiGEqChVJkAA5wOJWuu9AEqpWcBIoMIDRPz+b1kQ9Ko7CESQi0mVPrEtU4eSqqNIJZJUSwM257cghSiO6GjjK4Yjuj6ZhOM5H6FjRBTbThXvkxjQOoaJw9pjdzpZuuskh1Jz6NA4ivaNonji+80+azkF6oRZeWJ4e57+YSsNfIxaEkKIilCVAkQccMhjOwno448fFFonhnU6jixnGFmEkkUop7Tre5Yu3D5NKGk6gjQiycd1w27TIIKhnRvxxm+JZfwUl8s7NWKbxxDXvi3rs2pvKjf0jifEaubFq7vwyoKdvPFbIodSc/j7JW24rld8mee9pW9zbunb3GuJcCGEqEhVKUD4utMV+1ivlJoATABo1qzZWf2gsC5Xcd+skj+hl8ah9RnNSh7SoQGvLXTNiejdvJ57OK3nshd3DmhJ4oksbj6/abnPK4FBCOFvVWmiXBLgeYeMB44UzaS1fk9r3Vtr3Ts2NvasflDBfAWA+wa3olOTkm/457eo77VtczhpVuR5zZ3jSj6+UZ0QXjeW+jYphdNYrjvYUriAXp0wK++M7UWDKJnxLISoOqpSgFgLtFFKtVBKBQE3AT/6+4c+dnl7r3kFAN/f28/9+qu7LvDa1zgqlLi63ktk3zmgJR/c2tu9/dfTl2I1nvMcZDHRq3k9wNXJbDcm3Zmq0m9eCCF8qDJNTFpru1LqfmA+rmGuH2qtt/rr5/368EAijfkQGTk2r31BFt9379du7MaA1rHUDbNy18CWLNh2nH3JpzGbFJd0bEjPZnX582A6dUKtXNQ2loXbTxBsMdG0fhhr/jWE2Ihgd3+ErKwqhKjqqtRdSmv9E/BTZfysNh7Pem5UJ5S07MIgkVvkKWu/PHQhO4+dYmT3wjkJTwzvwIiuTXhrcaK7hvDx7edzJD0HpRRv3NyTQ2nZ7qakBpGu5qP/XNOZC1vH0EXmKQghqjildfmeW1AV9e7dW69bt+6cz3M8M5dlu07y2DebANjxwlAOpWZzNCOXgW3Prp9DCCGqKqXUeq1177LySUs40DAqhOFdGru3Q6xm2jSMlOAghKjVJEAYZEayEEJ4kwBhsEiAEEIILxIgDFKDEEIIbxIgDDIzWQghvFWpYa6B9tSIjlzQUp6+JoQQIAHCyx0DWgS6CEIIUWVIE5MQQgifJEAIIYTwSQKEEEIInyRACCGE8EkChBBCCJ8kQAghhPBJAoQQQgifJEAIIYTwqVo/D0IpdRI4cJaHxwDJFVic6kCuuXaQa64dzuWam2uty3yeQbUOEOdCKbWuPA/MqEnkmmsHuebaoTKuWZqYhBBC+CQBQgghhE+1OUC8F+gCBIBcc+0g11w7+P2aa20fhBBCiNLV5hqEEEKIUtTKAKGUGqqU2qmUSlRKTQx0eSqCUqqpUmqxUmq7UmqrUurvRnp9pdSvSqndxvd6RrpSSk01fgeblFI9A3sFZ08pZVZKbVBKzTW2WyilVhvX/KVSKshIDza2E439CYEs99lSStVVSn2jlNphvN8X1PT3WSn1sPF3vUUp9YVSKqSmvc9KqQ+VUieUUls80s74fVVKjTPy71ZKjTuXMtW6AKGUMgNvAcOAjsDNSqmOgS1VhbADj2itOwB9gfuM65oILNJatwEWGdvguv42xtcE4J3KL3KF+Tuw3WP7JeA145rTgDuM9DuANK11a+A1I191NAX4RWvdHuiG69pr7PuslIoDHgR6a607A2bgJmre+/wxMLRI2hm9r0qp+sAzQB/gfOCZgqByVrTWteoLuACY77H9BPBEoMvlh+v8AbgU2Ak0NtIaAzuN1+8CN3vkd+erTl9AvPGPczEwF1C4Jg9Zir7fwHzgAuO1xcinAn0NZ3i9UcC+ouWuye8zEAccAuob79tc4PKa+D4DCcCWs31fgZuBdz3SvfKd6Vetq0FQ+MdWIMlIqzGMKnUPYDXQUGt9FMD43sDIVlN+D68DjwNOYzsaSNda241tz+tyX7OxP8PIX520BE4CHxnNah8opcKpwe+z1vowMBk4CBzF9b6tp2a/zwXO9H2t0Pe7NgYI5SOtxgzlUkpFAN8CD2mtM0vL6iOtWv0elFIjgBNa6/WeyT6y6nLsqy4sQE/gHa11D+A0hc0OvlT7azaaSEYCLYAmQDiuJpaiatL7XJaSrrFCr702BogkoKnHdjxwJEBlqVBKKSuu4DBTa/2dkXxcKdXY2N8YOGGk14TfQ3/gKqXUfmAWrmam14G6SimLkcfzutzXbOyvA6RWZoErQBKQpLVebWx/gytg1OT3+RJgn9b6pNbaBnwH9KNmv88FzvR9rdD3uzYGiLVAG2MERBCuzq4fA1ymc6aUUsB0YLvW+lWPXT8CBSMZxuHqmyhIv9UYDdEXyCioylYXWusntNbxWusEXO/jb1rrMcBi4HojW9FrLvhdXG/kr1afLLXWx4BDSql2RtIQYBs1+H3G1bTUVykVZvydF1xzjX2fPZzp+zofuEwpVc+oeV1mpJ2dQHfKBKgjaDiwC9gD/CvQ5amgaxqAqyq5CfjL+BqOq+11EbDb+F7fyK9wjebaA2zGNUIk4NdxDtc/CJhrvG4JrAESga+BYCM9xNhONPa3DHS5z/JauwPrjPd6NlCvpr/PwHPADmAL8CkQXNPeZ+ALXH0sNlw1gTvO5n0FbjeuPRG47VzKJDOphRBC+FQbm5iEEEKUgwQIIUqp+SwAAAFKSURBVIQQPkmAEEII4ZMECCGEED5JgBBCCOGTpewsQgilVMFwQ4BGgAPXkhcA2VrrfgEpmBB+JMNchThDSqlngSyt9eRAl0UIf5ImJiHOkVIqy/g+SCm1VCn1lVJql1JqklJqjFJqjVJqs1KqlZEvVin1rVJqrfHVP7BXIIRvEiCEqFjdcD2fogtwC9BWa30+8AHwgJFnCq7nGJwHXGfsE6LKkT4IISrWWm2sdaSU2gMsMNI3A4ON15cAHV3LCgEQpZSK1FqfqtSSClEGCRBCVKw8j9dOj20nhf9vJlwPtMmpzIIJcaakiUmIyrcAuL9gQynVPYBlEaJEEiCEqHwPAr2Nh81vA+4OdIGE8EWGuQohhPBJahBCCCF8kgAhhBDCJwkQQgghfJIAIYQQwicJEEIIIXySACGEEMInCRBCCCF8kgAhhBDCp/8HR6gDSWowLJ0AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Create some toy data\n", "real_parameters = [0.015, 6]\n", "times = np.linspace(0, 1000, 1000)\n", "values = transformed_model.simulate(real_parameters, times)\n", "\n", "# Add noise\n", "values += np.random.normal(0, 10, values.shape)\n", "\n", "# Create an object with links to the model and time series\n", "problem = pints.SingleOutputProblem(transformed_model, times, values)\n", "\n", "# Select a score function\n", "score = pints.SumOfSquaresError(problem)\n", "\n", "# Select some boundaries\n", "boundaries = pints.RectangularBoundaries([0, -6.0], [0.03, 20.0])\n", "\n", "# Perform an optimization with boundaries and hints\n", "x0 = 0.01,5.0\n", "sigma0 = [0.01, 2.0]\n", "found_parameters, found_value = pints.optimise(\n", " score,\n", " x0,\n", " sigma0,\n", " boundaries,\n", " method=pints.CMAES,\n", ")\n", "\n", "# Show score of true solution\n", "print('Score at true solution: ')\n", "print(score(real_parameters))\n", "\n", "# Compare parameters with original\n", "print('Found solution: True parameters:' )\n", "for k, x in enumerate(found_parameters):\n", " print(pints.strfloat(x) + ' ' + pints.strfloat(real_parameters[k]))\n", "\n", "# Show quality of fit\n", "pl.figure()\n", "pl.xlabel('Time')\n", "pl.ylabel('Value')\n", "pl.plot(times, values, label='Nosiy data')\n", "pl.plot(times, problem.evaluate(found_parameters), label='Fit')\n", "pl.legend()\n", "pl.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 }