{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Beeler-Reuter AP Model\n", "\n", "In this example we use the [1977 Beeler-Reuter model](http://scholarpedia.org/article/Models_of_cardiac_cell#Beeler-Reuter_model_.281977.29_.5B8_variables.5D) of the mammalian ventricular action potential [action potential](https://en.wikipedia.org/wiki/Action_potential), and try to estimate its maximum conductance parameters from its output.\n", "\n", "The model describes several ion currents, each with a maximum conductance parameter, that together give rise to the cardiac AP and calcium transient. In this (non-trivial) \"toy\" model, we use the maximum conductances as parameters, and the AP and calcium transient as observable outputs. To make the problem tractable, all other model parameters are assumed to be known.\n", "\n", "The parameters are _scaled_: instead of passing in the conductances directly, users provide the natural log of the maximum conductances. This makes the parameters easier to find for our optimisation algorithm.\n", "\n", "For the science behind the model, see the [original 1977 paper](https://doi.org/10.1113/jphysiol.1977.sp011853).\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pints\n", "import pints.toy\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "model = pints.toy.ActionPotentialModel()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can get an example set of parameters using the `suggested_parameters()` method:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "x_true = model.suggested_parameters()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And it can also provide a suggested sequence of sampling times:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "times = model.suggested_times()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the suggested parameters and times, we can run a simulation:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "values = model.simulate(x_true, times)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gives us all we need to create a plot of current versus time:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEKCAYAAABdWiGrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4HNX18PHv2V31ajVXyXLv3WB6MRAMOJgWMBAgCQlpJJD8SIDwppAKhBIIhAQICSEE00JiesAYgzHuvfciN0m2ZctW1573jxnZslBZ29qdlXQ+zzPP7s7Ozj26lvfo3rlzr6gqxhhjjNd8XgdgjDHGgCUkY4wxUcISkjHGmKhgCckYY0xUsIRkjDEmKlhCMsYYExUsIRljjIkKlpCMMcZEBUtIxhhjokLA6wDakqysLM3Pz/c6DGOMaVMWLFhQrKrZLR1nCekY5OfnM3/+fK/DMMaYNkVEtoRynHXZGWOMiQqWkIwxxkQFS0gRpKrcNmURT3+80etQjDEm6lhCiqCNxYf47+Id/ObtVRSWVngdjjHGRBVLSBG0bnfp4edvLNnpYSTGGBN9LCFF0KHKWgDSEmL4z6LtHkdjjDHRxRJSBJVV1QBw3bg8lm3fz5pdpS18whhjOg5LSBF0qMppIV13ch4xfuGleds8jsgYY6KHJaQIKquqRQS6pyfwhSFd+PeiAiprar0OyxhjooIlpAgqq6whIcaPzydce1IeJWXVvL7QriUZYwxYQoqosupaEmOd2ZpO75vJiB5pPD59PVU1QY8jM8YY71lCiqCyyhoSY/0AiAi3n9+fgn3lPD87pGmejDGmXbOEFEFlVbUkxPgPvz5nQDbnDMjmwffWsGXPIQ8jM8YY71lCiqCgKn6fHH4tIvzuimEE/MIt/1jA/vLqiMShqhEpxxhjjoUtPxFBQQVfgz8BuqYl8Ocvj+Erf5vL9c/M5i83jKV7esIxn7u0opqNRYfYUHSQrXvL2H2gkqLSCgpLK9lfXk1ZVS0VVbWUV9dSE3QSUsAn+HyCX4TEWD+pCTGkxAecLS6GzORYclLiyU6JIycljpzUuMOv6ydWY4xpDZaQIiioil8+/0V+et8s/nLDGG57cTEXPvIxN5/RiytGdycvIxGpd7yqsutABRsKncRzeCs8xK4DR8+Nl5kUS05qPDkpcfTKSiIx1k9CTICEWB8xfh/BoFKrSm0QaoNByqpqOVBRQ2lFNaUVNRQeOMicTZXsK/t8qy3GL3RLTyC3UyK5GQn06JRIbkYiuZ0SyMtIJCMp9qi4jTEmFJaQIiioNPlFPX5gZ978/hn85q1VPDptHY9OW0daQgydU+Pw+3zsL6tiz6EqKuuNyEuJC9A7J5nT+2bRJyeJPtnJ9MlOJi8jkdhA6/TGVtUEKT5YSWFpJYUHnBbX9pJytu0tY9u+cv63Yjd7DlUd9Zm0hBh6ZyfROyuZPjnuY3YSPTOTWi0uY0z743lCEpEDLR0C7FTV/pGIJ5xUleZ6unpmJvHUjWPZtreMj9YUsmpXKXsPVlETDDK4ayoZSTHkZSbRJzuJvtnJZKfEhb0lEhvw0S09gW7NdCMeqqyhYJ+TpLbsLWNTsdNqm7m+iNcWFhw+zieQl5FI7+xk+uUk079zCgO6pNA3J5n4eoM9jDEdk+cJCdigqqOaO0BEFkUqmHAKquILIYHkZiRyw6n54Q+olSTFBRjQxUkuDZVWVLOp+NDh61t1jzPXFVNV67T2fAL5mUn075xC/y4pDHATVX5mIgG/taiM6SiiISFd2UrHRL1gkJASUnuSEh/D8B7pDO+RftT+mtogm/eUsXZ3KWt2Odva3aX8b+Uu3DEXxPp99MlJZkDnZAZ2TWVw11QGd0slKznOg5/EGBNu0ZCQfigiL6rqp00doKrtYonVoCodLB81KeD30Tcnmb45yVw8rOvh/RXVtawvPHg4Qa3ZXcrcTXv5z+Idh4/pnBrH4K6pDOmWxuBuTqLKy0jEZyP/jGnToiEhrQMeFJGuwEvAi6q62OOYwkIVGy7dgvgYP0O7pzG0e9pR+0vKqli58wArdzjbih0H+HhdMbVucyo5LsCgrilOknJbUv06JxMXsGtTxrQVnickVX0UeFREegKTgb+JSDzwIjBFVdd6GmArCqoSYwnpuKQnxnJanyxO65N1eF9FdS1rd5c6SWqnk6Renr+NMneZj4BP6Nc5heHd0xiem8aIHun075xiI/2MiVKeJ6Q6qroFuB+4X0RGAc8CPwfazZ+4oQ5qMKGJj/F/7vpUMKhs3nPocGtq2fb9vLdyFy/Nd9aeig34GNQ1lRE90hjWPY0Ruen0yU62lqsxUSBqEpKIxAATcFpJ5wEzgHsjUO4E4FGcxPeMqt4XrrKauw/JtA6fT+idnUzv7GQmDu8GOMPtt+0tZ0lBCcu272fJthJeW1DAPz5zJrVNivUzpHua25JKZ3j3NHpmJtq/lTER5nlCEpELgGuBS4C5wBTgFlUN+2yjIuIHngAuAAqAeSIyVVVXhqO8lu5DMuEhIuRlJpKXmcgXRzhJKhhUNhYfZMm2/U6SKijh+dlbqJy5CYCMpFhG53VidM90Rud1YkSPdBJi201j3Zio5HlCAn4C/Au4Q1X3Rrjsk4H1daP4RGQKMAkIS0IKascb9h2tfD6hb04KfXNSuHJMDwCqa4Os3V3Kkm37Wbh1Hwu37uODVbsB53rUoK6pjOnZiVF56Yzp2Ynu6QnWijKmFXmekFT13LrnItIJyKVeXKq6MIzFdwe21XtdAIwLV2FBayFFtRi/jyHd0hjSLY3rxuUBsO9QFYu27WPBln0s3FLCS/O28fdZmwHISYljTM9OjOnZiZN7ZTC4a6rdyGvMCfA8IdURkV8BXwE2AnUTtikwPpzFNrLvqLUZROQW4BaAvLy8EyrMriG1PZ2SYhk/sDPjB3YGnBt6V+8qZeFWN0lt3cc7y3cBztDzk/I7Ma53JuN6ZTC0exoxlqCMCVnUJCTgaqCPqla1eGTrKcBpkdXpAeyof4CqPgU8BTB27NgTWkjIriG1fQG/7/B9Uje60zvtPlDB7I17mLNpL3M27mH6miIAEmP9jM3PYFyvDE7pncGw7uk25NyYZkRTQloOpAOFESxzHtBPRHoB23FG+F0XrsJs2Hf71Dk1nkkjuzNpZHcACksrmLtpL3M27mXOpj38/r01gJOgxvXK4Mx+2ZzVP4s+2cnWYjamnmhKSL8DFonIcqCybqeqXhquAlW1RkRuBd7DGfb9rKquCFd5NqihY8hJiWfi8G6Hh53vOVjJ3E17mbVhDzPXFzN9jTNmpmtaPGf2y+KMftmc0TeLjKRYL8M2xnPRlJCew7kxdhlHriGFnaq+DbwdibJsLruOKTM5jouGdeUid86+bXvLmLm+mE/WFfHu8l28PL8AERjaLY2z+2czflAOI3uk29x8psOJpoRUrKqPeR1EOKm1kAzO8iLXnpzHtSfnURtUlhaUMHNdMR+vK+LJGRt4fPp6spJjOXdADucP7syZ/bJIjI2m/6rGhEc0/ZYvEJHfAVM5ussunMO+I8qGfZuG/D5hVF4nRuV14nvn9aOkrIoZa4v4YFUh767YxSsLCogN+DitTybnDerMhYM7k5Ma73XYxoRFNCWkukX6Tqm3L9zDviPKBjWYlqQnxh4eIFFdG2Tepr18sKqQaat389P/LOdn/13OSfkZTBzelQlDu5CTYsnJtB9Rk5Dq3yDbXgWDdh+SCV2M38dpfbM4rW8WP504iHWFB3l72U7eXraTn/13BT+fuoKT8zO4xJKTaSc8T0giMlFV3zzRY9oC67Izx0tEnCXeO6dw+/n9Wbe7lLeW7eTNpU5y+sXUFZzeN4urxvTgC4O72Lx7pk3yPCEBvxeR7TQ+a0Kd3wLtJCFZRjInrl/nFG53k9Pa3aW8sWQH/164ndumLCY5LsAlw7py5ZgenJTfyVrlps2IhoS0G3i4hWPWRSKQcAsqNpTXtLr+nVP4vy8M4Afn92f2pj28tmA7byzdwUvzt5GbkcCXxuQy+aRcGwxhop7nCUlVz/E6hkixqYNMOPl8cnhV3V9dNoR3l+/i1QUFPPz+Wh6bto4Lh3Th+lPyOLV3prWaTFTyPCF1JDZTg4mUxNgAV4zuwRWje7Cp+BAvzN7CKwsKeGvZTvrmJHP9uDyuGN2DtIQYr0M15jCb6TGCbFCD8UKvrCT+38TBzPnJefz+quEkxQW4942VnPa7afzyjZUU7CvzOkRjAGshRVQwqNZVYjwTH+PnS2Nz+dLYXJYV7OfZTzfxj88289xnm7lkWFduOas3Q7uneR2m6cCipoUkIoki8lMRedp93U9EJnodV2uyqYNMtBjWI41HrhnJxz8+l6+dns+HqwuZ+MeZXPvUbD5eW4TqCa20YsxxiZqEBPwNZ8qgU93XBcCvvQun9VmXnYk23dITuOeSwcy6ezw/uXggm4oPceOzc7niyVl8tKbQEpOJqGhKSH1U9QGgGkBVy2n+3qQ2x4Z9m2iVGh/DLWf1YcaPz+E3lw+l8EAlX/nbPC770yw+XL3bEpOJiGhKSFUikoC7hLiI9KHeJKvtgS0/YaJdXMDP9eN6Mv2Oc/jdFcPYc7CSr/19PpOe+JTpq63FZMIrmhLSz4F3gVwReQGYBvzY25Bal11DMm1FbMDHtSfnMf2Oc3jgyuGUlFXz1b/P45qnZrNw6z6vwzPtVNSMslPV90VkIc5s3wLcpqrFHofVquwakmlrYvw+rj4pl8tHd2fK3K08Om09V/xpFhcO6cyPLhxI35xkr0M07UjUJCQRGe0+3ek+5olIGrBFVWs8CqtV2Vx2pq2K8fu44dR8rhjdg2dnbuIvH2/k/ZUz+NKYXH5wQX+6pNm0RObERU1CAv4EjAaW4rSQhrrPM0XkW6r6Py+Daw1BteUnTNuWFBfge+f147pxeTwxfQP/nL2FqUt28O1z+nDLWb2Jj7FZxs3xi6ZrSJuBUao6VlXH4CzYtxw4H3jAy8BaQ93FYOuyM+1BZnIcP/viYKb939mMH5jDw++v5byHZvDW0p028MEct2hKSANVdUXdC1VdiZOgNnoYU6sJuv9HrcvOtCe5GYk8cf1optxyCqkJMXz3XwuZ/NRsVuzY73Vopg2KpoS0RkSeFJGz3e1PwFoRicO9N6ktC1oLybRjp/TO5M3vncFvLx/GusKDTPzjTH723+WUVrT5/7omgqIpIX0FWA/cDvwA2Ojuqwba/PLmdQnJriGZ9srvE64b5wwVv+nUfJ6fvYXzH57Bu8utG8+EJmoSkqqWq+pDqnq5ql6mqg+qapmqBlX1oNfxnSi1LjvTQaQlxPCLS4fwn++cTkZSHN/650K+8Y8F7Cgp9zo0E+WiJiG5k6m+KiIrRWRj3eZ1XK3FuuxMRzMiN503bj2dn1w8kJnri7jg4Rm8PG+btZZMk6ImIeFMrvokUIPTRfcP4HlPI2pFNqjBdEQBv49bzurD+z84m6Hd0/jxa0u55fkFFB9sV7OCmVYSTQkpQVWnAaKqW1T1F8B4j2NqNUeuIXkciDEeyM1I5MVvnMI9Fw9ixpoiJvzhY6at2u11WCbKRFNCqhARH7BORG4VkcuBHK+Dai01tU5CClifnemgfD7hG2f1Zur3Tic7JZ6bn5vP795ZRU1t0OvQTJSIpoR0O5AIfB8YA3wZuNHTiFpRVY3zny42YHeym45tYJdUXv/OaVw3Lo+/zNjIl/86h8LSCq/DMlEgmhJSvqoeVNUCVf2qql4J5HkdVGs5kpCiqcqN8UZ8jJ/fXj6Mh68eweJtJUx8bCaLbBbxDi+avh3vDnFfyETkFyKyXUQWu9vF9d67W0TWi8gaEbnwRMoJRVVtLQAxfuuyM6bOFaN78J/vnk58jJ/JT83m7WU7W/6Qabc8n1xVRC4CLga6i8hj9d5KxRlxd6IeUdUHG5Q5GJgMDAG6AR+ISH9VrW2F8hpVVeNcQ4qzFpIxR6nrwrvl+QV854WF3DlhIN86u7fdRN4BRcO34w5gAVDhPtZtU4FwtVwmAVNUtVJVN+HMEHFymMoCoKrWuuyMaUpmchwvfH0cE4d35f53V3PPf5ZTG7T7lToaz1tIqroEWCIi/wzTuke3isiNwHzg/1R1H9AdmF3vmAJ33+eIyC3ALQB5ecd/SevwNSS/DWowpjHxMX4emzyKHp0S+fOMDRyqrOHBL40gxm9/xHUUnickEVkGqPv8c++r6vAWPv8B0KWRt+7BudH2V+75fwU8BHwNZ72lzxXV2PlV9SngKYCxY8ce959sNqjBmJb5fMJdFw0kJT7A799bQ1lVLY9fN4o4G53aIXiekICJJ/JhVT0/lONE5GngTfdlAZBb7+0eOF2HYVM3qMESkjEt++65fUmOC/DzqSv4+nPzeeqGsSTEWlJq7zz/dnRnZdiiqltwriMNc7dyd99xE5Gu9V5ejrPgHzjXpyaLSJyI9AL6AXNPpKyW1LWQbJSdMaG56bR8HrhqODPXF3PL8/OpqA7bmCMTJTxPSHVE5GqcpPAl4GpgjohcdYKnfUBElonIUpz58X4A4C4E+DKwEngX+G44R9gBVNXaKDtjjtXVY3O574phfLKumFv/tfDwH3amfYqGLrs69wAnqWohgIhkAx8Arx7vCVX1hmbe+w3wm+M997GyQQ3GHJ9rTsqjqibIT/+7gttfWsRjk0cRsIEO7VI0JSRfXTJy7SGKWnAnygY1GHP8bjg1n8qaIL9+axUx/iU8fPVI/DYvZLsTTQnpXRF5D3jRfX0N8LaH8bSqqhob1GDMifj6mb2pqg3ywLtriAv4uO+K4fgsKbUrUZOQVPVHInIFcAbOsOynVPV1j8NqNXZjrDEn7jvn9KWiOshj09YRG/Dxq0lDbUaHdsTzhCQijwP/UtVZqvpv4N9exxQOZVVOCykhxq4hGXMifnB+Pyqra/nLxxuJD/i555JBlpTaCc8TErAOeMgdov0S8KKqLvY4plZXXl1LbMBn/d7GnCAR5+bZypogz8zcRHyMnzsuHOB1WKYVeJ6QVPVR4FER6Ykz4enfRCQe51rSFFVd62mAraS8qpZEu7HPmFYhIvxs4mAqa2p5fPp64mN83Dq+n9dhmRMUNRc03Jtj71fVUcB1ODeyrvI4rFZTVlVr3XXGtCKfT/jNZcO4YlR3HvzfWp75ZKPXIZkT5HkLqY6IxAATcFpJ5wEzgHs9DaoVlVfX2tQnxrQyn0944Krh9YaE+7jptHyvwzLHyfOEJCIXANcCl+DM1DAFuEVVD3kaWCuzLjtjwiPg9/GHySOpqg3y86krOFBeza3j+9pAhzYoGrrsfgJ8BgxS1S+q6gvtLRkBlFXVWJedMWES4/fxp+tHc8Wo7jz0/lrufWMlQVtPqc3xvIWkqud6HUMklFcHSUuI8ToMY9qtGL+PB780gk5Jsfx15iZ2lJTz8DUjSY7z/GvOhCgaWkgdQnlVDYnWQjImrHw+4f9dMoifThzMtNWFXP7Ep2wsOuh1WCZElpAipKzKBjUYEwkiws1n9OIfXzuZ4oOVTPzjTF6cuxVV68KLdpaQIqTCRtkZE1Gn983ize+fycjcdO7+9zK+9vd5bNnT7i5PtyuWkCKkrKrWuuyMibDu6Qn88+Zx/PyLg5mzaS8XPPwx972zmv1l1V6HZhphCSkCVNXuQzLGIz6f8NXTezH9jnP44ohu/HnGBk67bxq/fXsVO0rKvQ7P1GPDTyKgsiaIKpaQjPFQ59R4Hrp6BF8/sxdPfrSBZz7ZyNOfbOSMvllcOboH4wflkBpvI2G9ZAkpAupm+rYuO2O8N6hrKo9dO4ofXTiAVxcU8NrCAm5/aTEBnzCmZyfO6p/NqLx0hvdItyHjEWa1HQFlVTWAtZCMiSa5GYn84IL+3HZePxZu3cf0NYVMX13E799bA4BPoE92Mn2yk+mdnUSf7GTyMhPpnBJPTmoc8fYHZquzhBQBFdXuWkixVt3GRBufTxibn8HY/Ax+dOFASsqqWLythMXbSli+/QBrC0v5YNVuahrM/JAaHyAnNZ70hBhS4gOkxMeQmuA8JscFiAv4iIvxE+f3ERfjIy7gIzbgIy7gJ8bvw+8Dnwh+n+ATOfz8c/t9gl8En+AsXQoIR6ZFqj9DkhzeJ43sq3+cHP1mCOcJ+CTsSdi+ISPAuuyMaTvSE2M5Z0AO5wzIObyvujbI1r1lFOwrp/BABYWllRQeqGD3gUoOVFRTfLCKTcWHOFBRQ2lFNdW17e+ep4nDu/L4daPDWoYlpAjokhrPPRcPYkCXFK9DMcYchxi/73D3XUtUlaraIJU1QSqrg87z6lrndU2QKnerVSUYVGqDSlCdrTbI4f3O6yP7g+6NvXp0YZ/bV//+X23m/aP3fT6BHnUelN5ZLf/sJ8oSUgTkpMbzjbN6ex2GMSYCRIS4gJ+4gB/ivY6mbbH7kIwxxkQFS0jGGGOigtiEg6ETkSJgywmcIgsobqVwWpPFdWwsrmNjcR2b9hhXT1XNbukgS0gRJCLzVXWs13E0ZHEdG4vr2Fhcx6Yjx2VddsYYY6KCJSRjjDFRwRJSZD3ldQBNsLiOjcV1bCyuY9Nh47JrSMYYY6KCtZCMMcZEBUtIxhhjooIlpAgQkQkiskZE1ovIXR7HsllElonIYhGZ7+7LEJH3RWSd+9gpAnE8KyKFIrK83r5G4xDHY279LRWRsM3w2ERcvxCR7W6dLRaRi+u9d7cb1xoRuTCMceWKyHQRWSUiK0TkNne/p3XWTFye1pmIxIvIXBFZ4sZ1r7u/l4jMcevrJRGJdffHua/Xu+/nRziuv4vIpnr1NdLdH7Hffbc8v4gsEpE33deRrS9VtS2MG+AHNgC9gVhgCTDYw3g2A1kN9j0A3OU+vwu4PwJxnAWMBpa3FAdwMfAOzqz4pwBzIhzXL4A7Gjl2sPvvGQf0cv+d/WGKqysw2n2eAqx1y/e0zpqJy9M6c3/uZPd5DDDHrYeXgcnu/j8D33affwf4s/t8MvBSmOqrqbj+DlzVyPER+913y/sh8C/gTfd1ROvLWkjhdzKwXlU3qmoVMAWY5HFMDU0CnnOfPwdcFu4CVfVjYG+IcUwC/qGO2UC6iHSNYFxNmQRMUdVKVd0ErMf59w5HXDtVdaH7vBRYBXTH4zprJq6mRKTO3J/7oPsyxt0UGA+86u5vWF919fgqcJ5I/RWCwh5XUyL2uy8iPYBLgGfc10KE68sSUvh1B7bVe11A8/9hw02B/4nIAhG5xd3XWVV3gvMFA+Q0+enwaiqOaKjDW90uk2frdWl6EpfbPTIK56/rqKmzBnGBx3Xmdj8tBgqB93FaYyWqWtNI2Yfjct/fD2RGIi5Vrauv37j19YiIxDWMq5GYW9sfgB8DQfd1JhGuL0tI4dfYXw1ejrU/XVVHAxcB3xWRszyMJVRe1+GTQB9gJLATeMjdH/G4RCQZeA24XVUPNHdoI/vCFlsjcXleZ6paq6ojgR44rbBBzZTtWVwiMhS4GxgInARkAHdGMi4RmQgUquqC+rubKTsscVlCCr8CILfe6x7ADo9iQVV3uI+FwOs4/1F313UDuI+FHoXXVBye1qGq7na/RILA0xzpYopoXCISg/Ol/4Kq/tvd7XmdNRZXtNSZG0sJ8BHONZh0EalbB65+2Yfjct9PI/Su2xONa4Lb9amqWgn8jcjX1+nApSKyGeeywnicFlNE68sSUvjNA/q5o1VicS4ATvUiEBFJEpGUuufAF4Dlbjw3uYfdBPzXi/iaiWMqcKM74ugUYH9dN1UkNOizvxynzurimuyOOOoF9APmhikGAf4KrFLVh+u95WmdNRWX13UmItkiku4+TwDOx7m+NR24yj2sYX3V1eNVwIfqXrGPQFyr6/1RITjXaerXV9j/HVX1blXtoar5ON9RH6rq9US6vlprdIZtzY5cuRhn9NEG4B4P4+iNM8JpCbCiLhacvt9pwDr3MSMCsbyI05VTjfPX1s1NxYHTPfCEW3/LgLERjut5t9yl7n/ErvWOv8eNaw1wURjjOgOnS2QpsNjdLva6zpqJy9M6A4YDi9zylwM/q/d/YC7OYIpXgDh3f7z7er37fu8Ix/WhW1/LgX9yZCRexH7368V4DkdG2UW0vmzqIGOMMVHBuuyMMcZEBUtIxhhjooIlJGOMMVEh0PIhpk5WVpbm5+d7HYYxxrQpCxYsKFbV7JaOs4R0DPLz85k/f77XYRhjTJsiIltCOc667IwxxkSFsLaQRGQC8CjOjNfPqOp9Dd6PA/4BjAH2ANeo6mb3vbtx7gGpBb6vqu81d073JrspONNuLARuUNWqFsoYDvwFSMWZv+kkVa0IS2WYw4JBpaKmlorqIBXVtZRX11JR7bwOqhIMKkF17pELKs4+VVTB5xMCdZtf8Pt8h587+30kxPqJj/GTGOsnxm9/cxnTVoQtIYmIH+eGrgtwbjCcJyJTVXVlvcNuBvapal8RmQzcD1wjIoNx7hYeAnQDPhCR/u5nmjrn/cAjqjpFRP7snvvJZsoI4NyAdoOqLhGRTJybIc0xUlWKDlaya3+Fsx2oYOf+CopLK9lfXk1JeTUHyqspKaumpLyKiupgyydtJQGfkBDrJyHGf9RjanwMqQkxpCUESI2PIS2h7nXM4dedkmLISo4jPsYfsXiN6cjC2UI6vOwCgIjULbtQPyFNwlk3BZwpzB93p844PEU9sElE6k9R/7lzisgqnLmXrnOPec4975PNlPEFYKmqLgFQ1T2t9pO3Y4UHKli8rYTVu0rZUHSQDUUH2Vh0iLKq2qOOi/ELmUlxpCc6X+49MxMZ3sN5nhQXID7GT3zgSGsmLuAnPsZHwOfDJyAi+MRpEdW9FiCoUBtUamqD1ATVeR5UaoNBqmuV6togFdVByqtrKa+qoby6lrIqpwVWXuU8L6+upaSsiq17y9jvJsuaYNM3iKfEBchKiSMzKZas5DiyUmLJTIojKyWOLqnxdE2Lp1t6Ap0SYwjDigXGdBjhTEiNTZs+rqljVLVGROqmMO8OzG7w2bppzxs7Z8jTpNcroz+gIvIekI2TAB84vh+1/dpRUs6MtUXMXFfMoq372LH/SI9m9/QE+uQkc1J+BvmZSXQ98gDJAAAgAElEQVRNi6drWgJd0uLJTIrF52sbX86qSllVLQcqqt0EVUNJWRV7D1Wx51AVRaWVFB+sZM/BKjYUHWTOpkr2lX2+MR0f46NbWgLd0hOcukhPoEd6Aj0zE+mVlUR2SpwlLGOaEc6EFMr05E0d09T+xi4INHd8c2UEcObhOgkoA6aJyAJVnXZUgM6aQbcA5OXlNXKq9mfLnkP8Z9EO3lq2g7W7nbXEuqTGMza/E1/LTWdUXjqDuqaSGNs+BmmKCElxAZLiAnRNSwjpM9W1QfYeqmLX/gp27i9nR0kFO0rK2bm/gu0l5Xy8rojC0krqz8yVEOOnZ2Yi+ZlJ5GclkZ+ZSH5WEgM6p9ApKTZMP50xbUc4v1FCmTa97piCBlOYN/fZxvYX406T7raSGpsmvbEyZqhqMYCIvI2zdPVRCUlVnwKeAhg7dmy7nfivNqj8b8Uunv10E/M270METs7P4J6LB3H2gGz65STbX/f1xPh9dE6Np3NqPCNy0xs9pqomyM795WzZU8bmPYfYXFzGlj2HWFdYyoerC6mqPXItLSs5jgFdkumXk8KALin075xCv87JpMbHROpHMsZz4UxIh5ddALbjDFK4rsExdVOYf0a9KcxFZCrwLxF5GGdQQ90U9dLYOd3P1E2TPoXGp0lvWMZ7wI9FJBGoAs4GHglDPUS12qDyyvxtPPHRerbtLScvI5G7LhrIpSO60S09tNaCaVxswEfPzCR6ZiZxFkffE1gbVHaUlLOx+BDrdpeyZlcpa3eX8tK8bZRXH7kel5uRwPDu6Qzrkcaw7mkM7Z5GWoIlKdM+hS0huddrbgXewxmi/ayqrhCRXwLzVXUqzjoqz7uDFvbiJBjc417GGQBRA3xXVWsBGjunW+SdwBQR+TXO9O5/dfc3VcY+N+HNw+nCe1tV3wpXfUSjT9YV8es3V7FmdykjctO55+JBXDC4C/42cu2nLfP7hNyMRHIzEjm7/5FkFQwq20vKWbOrlDW7S1m54wBLt5fw1rIjS+DkZyYyrEc6I3qkcVJ+BkO6pRKw4e2mHbDlJ47B2LFjtT3M1FBaUc2v31zFS/O30TMzkbsmDGTC0C7WJRfF9h2qYtn2/Szbvp+lBSUsK9h/eIBJYqyfUXnpjO2Zwcm9MhiZm05SXPu4vmfaB/f6/NiWjrPf2g5maUEJ3/7nQnbuL+dbZ/fhBxf0Iy5g99lEu05JsZzVP5uz6rWmdu2vYN7mvczfvJe5m/fx2IfrUHVaX8N7pHFm3yzO7J/NyNx0u0HYtAnWQjoGbb2F9J9F27nztaVkJcfx2LWjGNOzk9chmVZ0oKKaBVv2MW/TXmZt2MPSghKCCslxAU7pncGZ/bI5s18WvbOTvQ7VdDDWQjJHeeaTjfz6rVWM65XBn64fTWZynNchmVaWGh/DuQNyOHdADgD7y6qZtaGYT9YXM3NdMR+sKgSgd1YS5w/uzPmDOjM6L92uP5moYS2kY9BWW0iPf7iOB/+3louHdeEP14wiNmBfQB3Rlj2HmLG2iA9WFfLZhmKqa5VOiTGcOzCHCwZ15sz+2STbtScTBqG2kJpNSO7w65bsVdWvHENsbVZbTEjPf7aZn/53BVeM6s4DVw23v4YN4Axs+WRdMR+s3M2HawopKasmNuDj7P7ZTBzelfMGdbbkZFpNa3XZDQK+3lw5OJOdmij03opd/GzqCs4flGPJyBwlJT6Gi4d15eJhXampDbJgyz7eXbGLd5bt4v2Vu4kN+Dh3QDaXDO/GeQNzbNSeiYiWWkhXq+rLzZ4ghGPai7bUQlq3u5RLH/+UAV1SePEbp5AQayPpTMuCQWXh1n28uXQnby/bSWFpJXEBH+cOyOGS4V25YHBnm/3cHLNW6bIzR2srCamsqoZJj3/K3kNVvH3bmXROjfc6JNMGBYPK/C37eGvpDt5evoui0kqS4wJMGNqFK0Z1Z1zvTLuJ2oQkIteQVPXS44itzWorCenHry7hlQUFPP+1cZzRL8vrcEw7UBtU5mzcw+uLtvPO8l0crKyhS2o8k0Z14/JR3RnYJdXrEE0Ua62EVISzdMOLwBwazJytqjNOMM42pS0kpBlri7jp2bl855w+/HjCQK/DMe1QeVUtH6zazeuLtjNjbRG1QWVQ11QuH9WNy0Z2J8da5KaB1kpIfpzVWa8FhgNvAS/Wmz+uQ4n2hHSosoYvPPIx8TE+3r7tTJuBwYTdnoOVvLl0J/9etJ0l20rw+4TzB+Vw/bienNE3q82siWXCq1VG2bkTmr4LvCsicTiJ6SMR+aWq/rF1QjWt5ZH317K9pJxXvnWqJSMTEZnJcdx0Wj43nZbPxqKDvDR/G6/ML+C9FbvJy0jk2pPzuHpsD7sR24SkxUENbiK6BCcZ5eMs5/Csqm4Pe3RRJppbSJuKD3HBwzO4akwP7rtyuNfhmA6ssqaWd5fv4oU5W5m7aS/xMT6uGZvL18/sTW5GotfhGQ+0SgtJRJ4DhgLvAPeq6vJWis+0sgfeXU1swMcPv9Df61BMBxcX8DNpZHcmjezOut2lPP3JRv41dyv/nLOVSSO68YML+ltiMo1q6RpSEDjkvqx/oACqqh1qaE20tpAWbNnLlU9+xg/O789t5/fzOhxjPmdHSTl/nbmJF+ZsIRiEm07rya3n9iMt0RYb7AjsPqQwiNaEdN3Ts1lXeJAZPzqHxFi7o95Er137K3j4/TW8sqCAtIQY7rl4EFeN6WFrcbVzoSakkOeSEZFOIjJcREbXbScWomkNC7fuY9aGPXzzrN6WjEzU65IWzwNXjeCt751Jv5xkfvTqUm762zwK9pV5HZqJAiElJBH5FbAU+CPwkLs9GMa4TIj+NH0D6YkxXHtyntehGBOywd1SeemWU/nlpCHM37yXi/7wCe8u39nyB027FmoL6Wqgj6qerarnutv4cAZmWrZ61wE+WLWbr57Wyya/NG2OzyfceGo+791+Fr2zk/jWPxfy6zdXUlMb9Do045FQE9JyID2cgZhj99dPNpEY6+em03p6HYoxxy03I5GXv3UqN53ak2dmbuIb/5jPocoar8MyHgg1If0OWCQi74nI1LotnIGZ5pWUVTF1yQ4uH9Wd9MRYr8Mx5oTEBfzcO2kov75sKDPWFjH5qdkUlVZ6HZaJsFAT0nPA/cB9HLmG9FBLHxKRCSKyRkTWi8hdjbwfJyIvue/PEZH8eu/d7e5fIyIXtnROEenlnmOde87Ylspw388TkYMickeIdREVXl1QQGVNkC+fYq0j0358+ZSePHXDWNYVlnLd07MpPmhJqSMJNSEVq+pjqjpdVWfUbc19wJ0H7wngImAwcK2IDG5w2M3APlXtCzyCk/Rwj5sMDAEmAH8SEX8L57wfeERV+wH73HM3WUY9j+Dc+NtmBIPKP2dvYWzPTgzq2qFuBTMdwPmDO/O3r5zMtn1lXP/0HPZYUuowQk1IC0TkdyJy6jEM+z4ZWK+qG1W1CpgCTGpwzCSc1hfAq8B54tyQMAmYoqqVqroJWO+er9Fzup8Z754D95yXtVAGInIZsBFoU5PFzlxfzOY9ZdxwqrWOTPt0ap9M/nrTSWzec4gbn51r15Q6iFAT0ijgFOC3hD7suzvO0hV1Ctx9jR6jqjXAfiCzmc82tT8TKHHP0bCsRssQkSTgTuDeFn6OqPOqe1PhhKFdvA7FmLA5vW8Wf/7yGFbvKuV7Ly6y0XcdQEgJqd5Q73OPYdh3Y7deN5wWoqljWmt/c2Xci9PFd7CR948EKHKLiMwXkflFRUXNHRoRBytr+N/KXXxxRFeb0du0e+cOzOHeS4fw4epC7n1jpdfhmDBrNiGJyMSWTtDMMQVAbr3XPYAdTR0jIgEgDdjbzGeb2l8MpLvnaFhWU2WMAx4Qkc3A7cBPROTWhj+Eqj6lqmNVdWx2dnYTP2rkvLNsJxXVQS4f1cPrUIyJiC+f0pNbzurN87O38NK8rV6HY8Kopbspfy8i22m8lVHnt8CbjeyfB/QTkV7AdpxBCtc1OGYqcBPwGXAV8KGqqjuk/F8i8jDQDegHzHXj+Nw53c9Md88xxT3nf5srAzizLggR+QVwUFUfb6E+PPf6ou3kZyYyOs9uCzMdx50TBrJyxwF++t8VDOmWxtDuaV6HZMKgpYS0G3i4hWPWNbZTVWvcFsd7gB9nDaUVIvJLYL6qTgX+CjwvIutxWi2T3c+uEJGXgZVADfBdd7FAGjunW+SdwBQR+TWwyD03TZXRFu3cX85nG/dw23n9bDJK06H4fcKjk0cy8Y8z+c4LC3nje2eQlmAzhbc3Ntv3MfB6tu9nPtnIr99axUd3nEN+VpJncRjjlQVb9nHNXz5j4vCu/GHyKK/DMSFq9dm+jffeXb6LQV1TLRmZDmtMz058b3w//rN4B28ttclY2xtLSG1E4YEKFmzdx0U21Nt0cN89tw8jctO55z/LKDxQ4XU4phVZQmoj3luxC1Xs3iPT4QX8Ph6+egQV1bX8+LWl2GWH9iOkNQvcKXsuAfLrf0ZVWxrwYFrJuyt20Ts7iX45yV6HYozn+mQnc+eEgdz7xkqmLtnBpJEN77k3bVGoLaQ3gK/gzIiQUm8zEbDvUBWzN+7loqFdbHSdMa4bT81nZG46v3xjJSVlVV6HY1pBqKu69VDV4WGNxDTpo7WF1AaVLwy27jpj6vh9wm8vH8YXH5/Jfe+s5r4r7SuqrQu1hfSOiHwhrJGYJk1fXURWcizD7GZAY44yuFsqXz+jF1PmbWPupr1eh2NOUKgJaTbwuoiUi8gBESkVkQPhDMw4aoPKjLVFnN0/B5/PuuuMaei28/vRo1MCP3l9GVU1NgFrWxZqQnoIOBVIVNVUVU1RVVuIJwIWb9vH/vJqzh3o/Tx6xkSjxNgAv5w0hPWFB3lu1mavwzEnINSEtA5Yrja+MuKmry7C7xPO7GsJyZimjB/YmXMHZPPotHUUltq9SW1VqAlpJ/CRu6z4D+u2cAZmHNPXFDImrxNpiTZvlzHN+enEwVTW1PL7d9d4HYo5TqEmpE3ANCAWG/YdMbsPVLBixwHOse46Y1rUOzuZr53ei1cWFLBkW4nX4ZjjENKwb1Vtc6uqtgcz1jgLAp47IMfjSIxpG24d35d/L9rOL95YwWvfOs0GArUxIbWQRGS6iHzYcAt3cB3dJ+uLyU6JY2AXa4waE4qU+BjunDCQRVtLeH3Rdq/DMcco1Btj76j3PB64EmedIhMmqspnG4o5o2+Wzc5gzDG4YlR3np+9hfveXc2FQ7uQHBfq15zxWkgtJFVdUG/7VFV/iLMEuAmTdYUHKT5YxWl9srwOxZg2xecT7r10CEWllfzxw0bXDzVRKtQuu4x6W5aIXAjYPDZhNGt9MQCn9sn0OBJj2p6Ruel8aUwPnp25iY1FB70Ox4Qo1FF2C4D57uNnwP8BN4crKAOzNuwhNyOB3IxEr0Mxpk368YSBxAf8/PLNlbZERRsRapddL1Xt7T72U9UvqOrMcAfXUdUGldkb93Bab+uuM+Z4ZafEcdv5/fhoTRHTVhV6HY4JQbNX+0RkvKp+KCJXNPa+qv47PGF1bCt3HOBARQ2n9bXuOmNOxE2n5TNl3jZ++eZKzuiXRXyM3+uQTDNaaiGd7T5+sZFtYhjj6tBmbXCvH/W2hGTMiYjx+/j5FwezdW8Zf525yetwTAuaTUiq+nP38auNbF9r6eQiMkFE1ojIehG5q5H340TkJff9OSKSX++9u939a9xBFM2eU0R6uedY554ztrkyROQCEVkgIsvcx/Et/TyRMmvDHvrmJJOTGu91KMa0eWf2y+bCIZ15/MP17Cgp9zoc04xQR9n9VkTS673uJCK/buEzfuAJ4CJgMHCtiAxucNjNwD5V7Qs8AtzvfnYwMBkYAkwA/iQi/hbOeT/wiKr2A/ZxZNBFo2UAxcAXVXUYcBPwfCh1EW5VNUHmbd7LaTa6zphW8/8uGUxQld+9s9rrUEwzQh1ld5GqHp4cSlX3ARe38JmTgfWqulFVq4ApwKQGx0wCnnOfvwqcJ85doJOAKapaqaqbgPXu+Ro9p/uZ8e45cM95WXNlqOoiVd3h7l8BxItIXEi1EUZLC0ooq6q1hGRMK8rNSOSbZ/fhjSU7mL1xj9fhmCaEmpD89b+sRSQBaOnLuzuwrd7rAndfo8eoag2wH8hs5rNN7c8EStxzNCyrqTLquxJYpKqVLfxMYTdrwx5EYFwvS0jGtKZvn92H7ukJ3PP6Miprar0OxzQi1IT0T2CaiNwsIl8D3udIq6Mpjc130/BmgKaOaa39LcYhIkNwuvG+2chxiMgtIjJfROYXFRU1dkirmrWhmMFdU+mUFBv2sozpSBJi/fz68qFsKDrEE9M3eB2OaUSo9yE9APwGGIRzXedX7r7mFAC59V73AHY0dYyIBIA0YG8zn21qfzGQ7p6jYVlNlYGI9ABeB25U1UZ/Q1X1KVUdq6pjs7PDuwxERXUtC7eUcHpfu//ImHA4d0AOl43sxpMfrWfNrlKvwzENhNpCQlXfUdU7VPX/VPW9ED4yD+jnjn6LxRmkMLXBMVNxBhQAXAV86K5KOxWY7I6Q6wX0A+Y2dU73M9Pdc+Ce87/NleEO0ngLuFtVPw21HsJpwZZ9VNUGbbogY8LoZ18c4swK/tpSaoM2g0M0aTYhiUipiBxoZCsVkQPNfda9XnMr8B6wCnhZVVeIyC9F5FL3sL8CmSKyHvghcJf72RXAy8BK4F3gu6pa29Q53XPdCfzQPVeme+4my3DP0xf4qYgsdjdPFx6ataGYgE84KT/DyzCMadcykmL52cTBLN5WwnOzNnsdjqlHbI6n0I0dO1bnz58ftvNf/qdP8Ynw2rdPC1sZxhhneZebn5vPp+uLeeN7Z9C/s605Fk4iskBVx7Z0XMhddu5Jc0Qkr247/vBMQ6UV1Swt2G+zMxgTASLC/VcOJzkuwPdfXGSj7qJEqDfGXioi64BNwAxgM/BOGOPqcOZu2kttUG3+OmMiJDsljt9/aTird5Xy4HtrvA7HEHoL6VfAKcBaVe0FnAdExUCA9mLWhj3EBXyMzuvkdSjGdBjjB3bmhlN68vQnm5i+2mYE91qoCalaVfcAPhHxqep0YGQY4+pwPl1fzNj8TjYbsTERds8lgxjcNZXbpixic/Ehr8Pp0EJNSCUikgx8DLwgIo8CNS18xoRoz8FKVu8qteXKjfFAfIyfv9wwBp9P+ObzCzhUaV9tXgk1IU0CyoAf4AzD3oCzBIVpBZ+5c2vZ/UfGeCM3I5HHJo9iXWEpt7+02O5P8khL9yH1FZHTVfWQqgZVtUZVnwMWA+nNfdaEbtaGPSTHBRjePc3rUIzpsM7qn83PJg7m/ZW7+dl/l9uy5x5oqYX0B6Cx+TXK3PdMK5i1vphxvTII+I9pFL4xppV95fRefPPs3rwwZyuPTVvvdTgdTrNLmAP5qrq04U5VnV9/MT1z/LaXlLN5TxlfPqWn16EYY4A7LxxIUWklj3ywFp/A987r53VIHUZLCam5JUsTWjOQjurT9c5y5TahqjHRwecTHrhyOCg89P5aqoPKD87vh7PsmgmnlhLSPBH5hqo+XX+niNwMLAhfWB3H9NWFdE6NY2AXm7rEmGgR8Pv4/ZdGEPALj01bR1FpBb+cNJQY61YPq5YS0u3A6yJyPUcS0FggFrg8nIF1BFU1QT5ZV8zE4V3try9joozfJ9x3xXByUuJ5fPp6tu4t44nrRpOeaGuVhUuz6V5Vd6vqacC9ONMFbQbuVdVTVXVX+MNr3+Zv2cvByhrOHejpJOPGmCb4fMIdFw7gwS+NYO6mvVz86CfM27zX67DarVAX6Juuqn90tw/DHVRHMX11IbF+H2fY9SNjotpVY3rw6rdOIybg45q/fMaD762hotomZG1t1iHqEVVl2upCxvXOICmupZ5TY4zXRuSm89b3z+TyUT14fPp6JvzhY2asLfI6rHbFEpJHVu48wMaiQ1w4pIvXoRhjQpQcF+Chq0fwwtfHISLc9Oxcrnt6Ngu2WDdea7CE5JGpS3YQ8AkXD+vqdSjGmGN0et8s3rntTH7+xcGs3X2QK5/8jOuens17K3ZRUxv0Orw2y/qKPBAMKm8s3sGZ/bLISLIRO8a0RfExfr56ei+uOSmX5z/bwnOzNvPN5xfQPT2BS0d245JhXRnSLdVG0B4DS0ge+HB1ITv2V/CTSwZ5HYox5gQlxgb45tl9uPmMXkxbXci/5mzlqY838uRHG8jPTOTs/tmc2ieLU3tnkpYY43W4Uc0SUoSpKn+esYFuafFMsOtHxrQbAb+PC4d04cIhXdh3qIr/rdzF28t28fL8Ap77bAsiMKBzCoO7pTK4aypDuqXRNyeZrORYa0W5LCFF2CsLCpi/ZR+/vXyYTaZqTDvVKSmWa07K45qT8qiqCbKkoIRZ6/eweNs+Zq4r5t8Ltx8+NjHWT15GIrkZiXRLiycrOY6slDgyk2IPPybHBUiODxAXaN8LeIY1IYnIBOBRwA88o6r3NXg/DvgHMAbYA1yjqpvd9+4GbgZqge+r6nvNnVNEegFTgAxgIXCDqlYdTxmtbe3uUu54ZQlBVVbuOMC4Xhlcc1JuOIoyxkSZ2ICPk/IzOCk/4/C+otJKVu48wKaig2zdW87WvWVs2XOI2Rv3UFrR9AKBMX4hKS7gJKi4AHExfuL8PmICQqzfR4zfR2zAd9Rzv0/wieD3gU8EEcEnznOf4L529/kEqf8eTstNBHpnJzF+YOew1lXYEpKI+IEngAuAApx58aaq6sp6h90M7FPVviIyGbgfuEZEBgOTgSFAN+ADEenvfqapc94PPKKqU0Tkz+65nzzWMlS11e92C/iETHfwwo2n5nPHhQPw+6yJbkxHlZ0Sx9kp2ZzdP/tz71VU17L3UBXFByspPljJnoNVHKqs4WBlDQcrazlUWcOhyhpKK2uorAlSXROkojpIaUUNVTVBqmqDVNUEqXYfa4KKKgRV3c25dBB094W67NPE4V3bbkICTgbWq+pGABGZgrPybP2ENAn4hfv8VeBxcTpTJwFTVLUS2CQi693z0dg5RWQVMB64zj3mOfe8Tx5HGZ+1VgXU6Z2dzN++enLLBxpjOrz4GD/d0hPolh6ZBRVU6yesI0mq1s1UdQsVBnzhv8QQzoTUHdhW73UBMK6pY1S1RkT2A5nu/tkNPtvdfd7YOTOBElWtaeT44ynDGGM6BBG3mw7ve23CmfIa++kaNg6bOqa19h9PGUcHKHKLiMwXkflFRTZNiDHGhEs4E1IBUP/KfQ9gR1PHiEgASAP2NvPZpvYXA+nuORqWdaxlHEVVn1LVsao6Njv78/29xhhjWkc4u+zmAf3c0W/bcQYQXNfgmKnATTjXba4CPlRVFZGpwL9E5GGcAQf9gLk4rZrPndP9zHT3HFPcc/73OMto0oIFC4pFZMtx1whk4STPaGNxHRuL69hYXMemPcbVM6SjnAta4dmAi4G1wAbgHnffL4FL3efxwCvAepxk0LveZ+9xP7cGuKi5c7r7e7vnWO+eM+54ywhjfcwPdxkWl8UVbZvFZXGFuolbkIkAEZmvqmO9jqMhi+vYWFzHxuI6Nh05LpsqwBhjTFSwhBRZT3kdQBMsrmNjcR0bi+vYdNi4rMvOGGNMVLAWkjHGmKhgCSkCRGSCiKwRkfUicpfHsWwWkWUislhE5rv7MkTkfRFZ5z52ikAcz4pIoYgsr7ev0TjE8Zhbf0tFZHSE4/qFiGx362yxiFxc77273bjWiMiFYYwrV0Smi8gqEVkhIre5+z2ts2bi8rTORCReROaKyBI3rnvd/b1EZI5bXy+JSKy7P859vd59Pz/Ccf1dRDbVq6+R7v6I/e675flFZJGIvOm+jmx9eT2UsL1vOLOSb8AZlh4LLAEGexjPZiCrwb4HgLvc53cB90cgjrOA0cDyluLAGer/Ds59aKcAcyIc1y+AOxo5drD77xkH9HL/nf1hiqsrMNp9noJz68Ngr+usmbg8rTP35052n8cAc9x6eBmY7O7/M/Bt9/l3gD+7zycDL4WpvpqK6+/AVY0cH7Hffbe8HwL/At50X0e0vqyFFH6HJ5lV1SqcG3cneRxTQ5NwJqTFfbws3AWq6sc4M2aEEsck4B/qmI0zK0fXCMbVlMMT9KrqJpx73cIyi66q7lTVhe7zUmAVztyLntZZM3E1JSJ15v7cB92XMe6mOJMwv+rub1hfdfX4KnCeSOuvmtdMXE2J2O++iPQALgGecV8LEa4vS0jh19gks15O4qrA/0RkgYjc4u7rrKo7wfmCAXI8iq2pOKKhDm91u0yerdel6UlcbvfIKJy/rqOmzhrEBR7Xmdv9tBgoBN7HaY2FNAkzUDcJc9jjUtW6+vqNW1+PiLOO21FxNRJza/sD8GMg6L4OedJqWqm+LCGFX0iTuEbQ6ao6GrgI+K6InOVhLKHyug6fBPoAI4GdwEPu/ojHJSLJwGvA7ap6oLlDG9kXttgaicvzOlPVWlUdiTNP5cnAoGbK9iwuERkK3A0MBE7CWWT0zkjGJSITgUJVXVB/dzNlhyUuS0jhF9IkrpGiqjvcx0LgdZz/qLvrugHcx0KPwmsqDk/rUFV3u18iQeBpjnQxRTQuEYnB+dJ/QVX/7e72vM4aiyta6syNpQT4COcazLFOwhyJuCa4XZ+qzvpsfyPy9XU6cKmIbMa5rDAep8UU0fqyhBR+hyeZdUeoTMaZ8DXiRCRJRFLqngNfAJZzZAJaOHpi2khrKo6pwI3uiKNTgP113VSR0KDP/nKcOquLa7I74qgXIUzQewIxCPBXYJWqPlzvLU/rrKm4vK4zEckWkXT3eQJwPs71rbpJmKHxSZih3iTMEYprdb0/KgTnOk39+gr7v6Oq3q2qPVQ1H+c76kNVvZ5I11drjc6wrdmRK41OCOtBHL1xRjgtAVZwZMLbTGAasM59zIhALC/idOVU4+7hNq0AAAKwSURBVPy1dXNTceB0Dzzh1t8yYGyE43reLXep+x+xa73jIzJBL3AGTpfIUmCxu13sdZ01E5endQYMBxa55S8Hflbv/8AxTcIcobg+dOtrOfBPjozEi9jvfr0Yz+HIKLuI1pfN1GCMMSYqWJedMcaYqGAJyRhjTFSwhGSMMSYqWEIyxhgTFSwhGWOMiQqWkIwxxkQFS0jGhJmIZNZbVmCXHL0sw6wwlPcVESkSkWda4VzXuEsMvNkasRnTnEDLhxhjToSq7sGZ0w0R+QVwUFUfDHOxL6nqrSd6ElV9SUR2A3e0QkzGNMtaSMZ4SEQOuo/niMgMEXlZRNaKyH0icr04i7ktE5E+7nHZIvKaiMxzt9NDKOMrIvIfEXlDnEXgbhWRH4qzENtsEclwj/u+iKx0Z5yeEt6f3JjPsxaSMdFjBM6M1HuBjcAzqnqyOKuwfg+4HXgUeERVZ4pIHvAejc9i3dBQnKUh4nGme7lTVUeJyCPAjTgTad4F9FLVyrr51oyJJEtIxkSPeepOnCkiG/5/e3ePEjEURmH4PWPr3wbcgILY2inYiOIaLK2sbF2Dla5A+wG7aQZcwrgQEYvBZq5FgpPGPyYjmeF9qpsQPpLi8nGTGw4wqM8/A4f1+AjYbmShrSdZK1U43neG9TVvSV6Bx0bt3Xo8Ah6S9IH+zE8j/ZENSeqO98Z40jieMJ2rPWC/lDKeQ+0Tqgj3M+A6yU6ZhrNJc+c3JGmxDIDPzQpJ9toomqQHbJVShlSpoZvAahu1pd9yhSQtlkvgNsmIav4+ARct1F0B7pNsUEUe3JQqQE76N8ZPSEsmyTlVbs7M277regfAVSnltI160ld8ZSctnzFw3NaPscAd8DLzXUk/cIUkSeoEV0iSpE6wIUmSOsGGJEnqBBuSJKkTbEiSpE74AOk2bu6juy+pAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "_, axes = plt.subplots(2, 1)\n", "axes[0].plot(times, values[:, 0])\n", "axes[1].plot(times, values[:, 1])\n", "axes[0].set_ylabel('Voltage [mV]')\n", "axes[1].set_ylabel('Calcium [mM]')\n", "axes[1].set_xlabel('Time [ms]')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will add some noise to generate some fake \"experimental\" data and try to recover the original parameters." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running...\n", "Minimising error measure\n", "using Covariance Matrix Adaptation Evolution Strategy (CMA-ES)\n", "Running in sequential mode.\n", "Population size: 8\n", "Iter. Eval. Best Time m:s\n", "0 8 12.49416 0:00.9\n", "1 16 11.57587 0:01.9\n", "2 24 11.57587 0:02.9\n", "3 32 11.45721 0:03.8\n", "10 80 11.45721 0:08.8\n", "Halting: Maximum number of iterations (10) reached.\n", "Found solution: True parameters:\n", " 4.17219904709831546e+00 4.20508438550409647e+00\n", " 3.14367650342651161e-03 3.15381328912807055e-03\n", " 9.39651794702475590e-02 9.46143986738421233e-02\n", " 3.64168052428662792e-01 3.67944883731608385e-01\n", " 8.17133456668796998e-01 8.41016877100819293e-01\n" ] } ], "source": [ "# Add noise\n", "values[:, 0] += np.random.normal(0, 1, values[:, 0].shape)\n", "values[:, 1] += np.random.normal(0, 5e-7, values[:, 1].shape)\n", "\n", "# Create an object with links to the model and time series\n", "problem = pints.MultiOutputProblem(model, times, values)\n", "\n", "# Add weighting to the outputs\n", "# (we believe both outputs should contribute roughly equally to the score function\n", "# so we scale it with roughly the values that it spans)\n", "weights = [1. / 70., 1. / 0.000006]\n", "\n", "# Select a score function\n", "score = pints.SumOfSquaresError(problem, weights=weights)\n", "\n", "# Select some tight boundaries\n", "lower = x_true - 1\n", "upper = x_true + 1\n", "boundaries = pints.RectangularBoundaries(lower, upper)\n", "\n", "# Perform an optimization\n", "# For this example, we will start quite close to the true solution to limit the number of iterations\n", "x0 = x_true + 0.25\n", "optimiser = pints.OptimisationController(score, x0, boundaries=boundaries, method=pints.CMAES)\n", "\n", "print('Running...')\n", "found_parameters, found_score = optimiser.run()\n", "\n", "# Compare parameters with original\n", "print('Found solution: True parameters:' )\n", "for k, x in enumerate(found_parameters):\n", " # We take an exponential here to show the actual conductance values\n", " print(pints.strfloat(np.exp(x)) + ' ' + pints.strfloat(np.exp(x0[k])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can create a figure to show how the obtained parameters differ from the original ones:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEACAYAAACDEBA8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF8pJREFUeJzt3X+U3XV95/Hnq6ClImiCAtmEiNZ0Oep6XBmBbldFQ0JaPUJdaXV1nXahWTmu2rp7WFzcRmGpEVS21mKbhazR2voDf5AexRiidLdbsAwKRV01sVWZkoVI0COr/LLv/eN+p1zHO5k7N3e+czPzfJwz597P536+3/v+TObw4vu9n/v9pqqQJKktP7PQBUiSlhaDR5LUKoNHktQqg0eS1CqDR5LUKoNHktSqkQqeJBuSfD3JniQX9nj9Z5N8uHn9C0lO7HrtTU3/15Oc2e8+JUntGpngSXIY8IfALwNPA16R5GnThp0L3FtVTwWuAN7ebPs04OXA04ENwJVJDutzn5KkFo1M8ACnAHuq6m+r6kHgQ8BZ08acBWxrnl8DrE2Spv9DVfVAVf0dsKfZXz/7lCS1aJSCZyVwR1d7sunrOaaqHga+DxxzgG372ackqUWHL3QBXdKjb/r1fGYaM1N/r2DteY2gJBuBjQBHHnnkySeddNLMlR7A7X///YG2G0X/bOXj5rzNYpn/Up47LO35L+W5w2Dzn3LLLbd8t6qeONu4UQqeSeCErvYq4M4ZxkwmORx4HLB/lm1n2ycAVbUF2AIwNjZWExMTA03ixAs/NdB2o2hi84vmvM1imf9Snjss7fkv5bnDYPOfkuTb/YwbpVNtNwNrkjw5yaPpLBbYPm3MdmC8ef4y4HPVucrpduDlzaq3JwNrgL/uc5+SpBaNzBFPVT2c5N8DO4DDgK1V9ZUkFwMTVbUduBr4QJI9dI50Xt5s+5UkHwG+CjwMvLaqfgzQa59tz02S9IiRCR6Aqvo08Olpfb/b9fx+4JwZtr0UuLSffUqSFs4onWqTJC0BBo8kqVUGjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVBo8kqVUGjySpVSN162sd+r61+UULXYKkEecRjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVIxE8SZYn2Zlkd/O4bIZx482Y3UnGu/pPTnJ7kj1J3p0kTf/lSb6W5G+SfCLJ49uakySpt5EIHuBCYFdVrQF2Ne2fkGQ5sAk4FTgF2NQVUO8FNgJrmp8NTf9O4BlV9UzgG8Cb5nMSkqTZjUrwnAVsa55vA87uMeZMYGdV7a+qe+mEyoYkK4Cjq+rGqirg/VPbV9Vnq+rhZvubgFXzOQlJ0uxGJXiOq6q9AM3jsT3GrATu6GpPNn0rm+fT+6f7t8B1Q6lWkjSw1q5ckOR64PgeL13U7y569NUB+rvf+yLgYeCDB6hvI53TdaxevbrPkiRJc9Va8FTVGTO9luSuJCuqam9z6uzuHsMmgdO72quAG5r+VdP67+za9zjwYmBtcypupvq2AFsAxsbGZhwnSTo4o3KqbTswtUptHLi2x5gdwPoky5pFBeuBHc2puR8kOa1Zzfbqqe2TbAD+E/CSqvrhfE9CkjS7UQmezcC6JLuBdU2bJGNJrgKoqv3AJcDNzc/FTR/A+cBVwB7gmzzyWc57gKOAnUluTfJHLc1HkjSDkbg6dVXdA6zt0T8BnNfV3gpsnWHcM3r0P3W4lUqSDtaoHPFIkpYIg0eS1CqDR5LUKoNHktQqg0eS1CqDR5LUKoNHktQqg0eS1CqDR5LUKoNHktQqg0eS1KqRuFabtBh8a/OLFroE6ZDgEY8kqVUGjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVLqceMpfUStKBecQjSWqVwSNJapXBI0lqlcEjSWqVwSNJatXIBE+S5Ul2JtndPC6bYdx4M2Z3kvGu/pOT3J5kT5J3J8m07f5jkkryhPmeiyRpZiMTPMCFwK6qWgPsato/IclyYBNwKnAKsKkroN4LbATWND8burY7AVgHfGc+JyBJmt0oBc9ZwLbm+Tbg7B5jzgR2VtX+qroX2AlsSLICOLqqbqyqAt4/bfsrgAuAmrfqJUl9GaXgOa6q9gI0j8f2GLMSuKOrPdn0rWyeT+8nyUuAv6+q2+ajaEnS3LR65YIk1wPH93jpon530aOvZupP8phm3+v7qG0jnVN1rF69us9yJElz1WrwVNUZM72W5K4kK6pqb3Pq7O4ewyaB07vaq4Abmv5V0/rvBH4eeDJwW7PWYBXwxSSnVNX/nVbbFmALwNjYmKfkJGmejNKptu3A1Cq1ceDaHmN2AOuTLGsWFawHdjSn5n6Q5LRmNdurgWur6vaqOraqTqyqE+kE1LOnh44kqT2jFDybgXVJdtNZgbYZIMlYkqsAqmo/cAlwc/NzcdMHcD5wFbAH+CZwXbvlS5L6MTJXp66qe4C1PfongPO62luBrTOMe8Ys73HiQRcqSToocz7iSXJkksPmoxhJ0uI3a/Ak+Zkk/zrJp5LcDXwN2JvkK0kuT7Jm/suUJC0W/RzxfJ7O6rA3AcdX1QlVdSzwXOAmYHOSV81jjZKkRaSfz3jOqKqHpnc2H+p/DPhYkkcNvTJJhxTvvqt+zXrE0yt0BhkjSRIMuJw6yR8meV/zfNarAkiSNGXQ7/E8CPxt8/yFQ6pFkrQEDBo8PwQe13y244XNJEl9GzR4NtG5OsCVwJ8OrxxJ0mI36JULDq+qK4daiSRpSRg0eN6Z5Ofo3JLga1V1+RBrkiQtYgMFT1W9FiDJUfS4RbUkSTMZdDn1RUnWAkcwQhcalSSNvkFDI8AK4Ew6126TJKkvg65q+2Xg8cDVzW0KJEnqy6DB81I6N1z71ambtEmS1I9BT7W9sqreBXwmyUnDLEiStLjNKXiSPB64AvinSe4HbqNzd9DfnIfaJEmL0JyCp6q+B/xmkjOB7wLPBD4+H4VJkhanWYMnyXLgiKq6c6qvqnY0T2+Zr8IkSYtTP4sL3gGMTzWS/FWSjyS5MMnK+StNkrQY9RM8JwObu9pHAVcDT6BzO2xJkvrWz2c8D1RVdbU/V1U7knwWuHGe6pIkLVL9HPHcn+RJU42qekPzWMCj5qswSdLi1E/wXAp8cvr3dZKswOu0SZLmaNbgaVaw/R7w+STXJbk8yeXAX/KTn/0MLMnyJDuT7G4el80wbrwZsztJ94KHk5PcnmRPkncnSddrr0vy9SRfSXLZMOqVJA1u1uBJ8ovANcDP01lUcB+wD3hpVf3ZkOq4ENhVVWuAXfS41UKzrHsTcCpwCrCpK6DeC2wE1jQ/G5ptXgCcBTyzqp5OZ4WeJGkB9XOqbZzO93W2Ao8F/riqLquq24ZYx1nAtub5NuDsHmPOBHZW1f6quhfYCWxoTvkdXVU3Np87vb9r+/OBzVX1AEBV3T3EmiVJA+jnVNtrqurZwFuAZcD7ktyY5PeSPC/JYUOo47iq2tu8317g2B5jVgJ3dLUnm76VzfPp/QC/ADw3yReS/EWS58xUQJKNSSaSTOzbt+8gpiJJOpC+FwdU1dfo3Hvniua21y8AzgHeBYzNtn2S64Hje7x0UZ8lpEdfHaAfOvNbBpwGPAf4SJKnTFse3tmgaguwBWBsbOynXpckDcecV6UleSVwc1V9Gvh0v9tV1RkH2OddSVZU1d7m1FmvU2KTwOld7VXADU3/qmn9d3Zt8/EmaP46yT/Q+eKrhzSStEAGuR/PPuDKJP8zyceTDGNl23YeuSzPOHBtjzE7gPVJljWLCtYDO5pTcz9Iclqzmu3VXdt/EnghQJJfAB5N5+KmkqQFMufgqarPAl+oqufRCYnHDqGOzcC6JLuBdU2bJGNTN5qrqv3AJcDNzc/FTR90FhFcRefmdN8Ermv6twJPSfJl4EPAeK/TbJKk9gz6BdCjk5wM3A4cebBFVNU9wNoe/RN07vcz1d5KJ0x6jXtGj/4HgVcdbH2SpOEZ9NbXbwR+Cfgj4DPDK0eStNgNGjwfoPMlzkfRuXq1JEl9GfRU241V9fsASY4ZYj2SpEVu0OA5q1mavKOqvjHMgiRJi9ugp9r+DZ3VY/9qatWZJEn9GPSI5565foFUkiQYPHje2Vw2p4CvVdXlQ6xJkrSIDRQ8VfVagCRH0eMWBpIkzWSgz3iS/Icka4Ej8PbXkqQ5GPRU2w7gWXSul/Z/hleOJGmxGzR4fh04CfghPS5hI0nSTAZdTv34qjoH+C3g9UOsR5K0yA0aPA8meTadVW0HfZFQSdLSMefgae5rcxGdq0lvAT487KIkSYvXIJ/x/BbwTb+7I0kaxCDBsx84P8lJwG3ArVX1peGWJUlarGYNniTLgSOq6k6Aqnpbkl3AN+gsqX4uYPBIkvrSzxHPO4DdwNsAkvwVMAl8EfhAVd0wb9VJkhadfhYXnAxs7mofBVwNPAF403wUJUlavPo54nmgqqqr/bmq2pHks8CN81SXJGmR6ueI5/4kT5pqVNUbmsfC67RJkuaon+C5FPhks4rtHyVZweCX3JEkLVGzBkdzWu1o4PNJbgW+3Lz0UuDN81mcJGnx6euIpao+muRTwK8ATwd+BLy0qm6bz+IkSYvPrKfakgSgqn5YVddU1Vur6rLu0JkaczCSLE+yM8nu5nHZDOPGmzG7k4x39Z+c5PYke5K8e6qmJM9KclOSW5NMJDnlYGuVJA2un894Pp/kdUlWd3cmeXSSFybZBozPsO1cXAjsqqo1wC563Nm0+TLrJuBU4BRgU1dAvRfYCKxpfjY0/ZcBb62qZwG/27QlSQukn+DZAPwY+LMkdyb5apK/o/Ol0lcAV1TV+4ZQy1nAtub5NuDsHmPOBHZW1f6quhfYCWxoFjocXVU3Nqvt3t+1fQFHN88fB9w5hFolSQPqZ3HB/cCVwJVJHkXni6M/qqrvDbmW46pqb/Oee5Mc22PMSuCOrvZk07eyeT69H+C3gR1J3kEnaP9FrzdPspHOEROrV6/uNUSSNARzWg5dVQ8Bewd9syTXA8f3eOmifnfRq6wD9AOcD/xOVX0sya/RuerCGT81uGoLnds8MDY2VtNflyQNR9/Bk+SNPbq/D9xSVbf2s4+q+qn/4Hft/64kK5qjnRXA3T2GTQKnd7VXATc0/aum9U+dUhsH3tA8/yhwVT+1SpLmx1xuBDcGvIZHTm1tpBMC/z3JBUOoZTuPLFIYB67tMWYHsD7JsmZRwXpgR3OK7gdJTmtWs726a/s7gec3z19I57MpSdICmcuptmOAZ1fVfQBJNgHXAM8DbuHgV4ttBj6S5FzgO8A5zfuMAa+pqvOqan+SS4Cbm20urqr9zfPzgfcBPwdc1/xA58Z1v5/kcOB+ms9xJEkLYy7Bsxp4sKv9EPCkqvpRkgcOtpCquofO7bSn908A53W1twJbZxj3jB79f0nnCtuSpBEwl+D5U+CmJNfS+TD/xXSWWB8JfHU+ipMkLT59B09VXZLk08C/pBM8r2mOMgBeOR/FSZIWn7leXfph4B/oLFV+aPjlSJIWu75XtSV5A/BBOl8gPRb4kySvm6/CJEmL01yOeM4FTq2q/weQ5O107kD6B/NRmCRpcZrL93hC55ptU35M7ysGSJI0o7kc8fwP4AtJPkEncM6mx7JmSZIOZC6r2t6V5Abgl+gEz3i/l8qRJGnKrMGT5Ac8csFN6Dq9lqSq6uif3kqSpN76uS3CUW0UIklaGuayuECSpINm8EiSWmXwSJJaZfBIklpl8EiSWmXwSJJaZfBIklpl8EiSWmXwSJJaZfBIklpl8EiSWmXwSJJaZfBIklpl8EiSWjUSwZNkeZKdSXY3j8tmGDfejNmdZLyr/9IkdyS5b9r4n03y4SR7knwhyYnzOxNJ0mxGIniAC4FdVbUG2NW0f0KS5cAm4FTgFGBTV0D9edM33bnAvVX1VOAK4O3zULskaQ5GJXjOArY1z7cBZ/cYcyaws6r2V9W9wE5gA0BV3VRVe2fZ7zXA2iTpMU6S1JJRCZ7jpoKjeTy2x5iVwB1d7cmm70D+cZuqehj4PnBMr4FJNiaZSDKxb9++OZYvSerXrLe+HpYk1wPH93jpon530aOvhrVNVW0BtgCMjY3Ntl9J0oBaC56qOmOm15LclWRFVe1NsgK4u8ewSeD0rvYq4IZZ3nYSOAGYTHI48Dhg/1zqliQN16icatsOTK1SGweu7TFmB7A+ybJmUcH6pq/f/b4M+FxVeTQjSQtoVIJnM7AuyW5gXdMmyViSqwCqaj9wCXBz83Nx00eSy5JMAo9JMpnkLc1+rwaOSbIHeCM9VstJktrV2qm2A6mqe4C1PfongPO62luBrT3GXQBc0KP/fuCcoRYrSTooo3LEI0laIgweSVKrDB5JUqsMHklSqwweSVKrDB5JUqsMHklSqwweSVKrDB5JUqsMHklSqwweSVKrDB5JUqsMHklSqwweSVKrDB5JUqsMHklSqwweSVKrDB5JUqsMHklSqwweSVKrDB5JUqsMHklSqwweSVKrDB5JUqtGIniSLE+yM8nu5nHZDOPGmzG7k4x39V+a5I4k900b/8YkX03yN0l2JXnSfM9FknRgIxE8wIXArqpaA+xq2j8hyXJgE3AqcAqwqSug/rzpm+5LwFhVPRO4BrhsHmqXJM3BqATPWcC25vk24OweY84EdlbV/qq6F9gJbACoqpuqau/0Darq81X1w6Z5E7Bq6JVLkuZkVILnuKngaB6P7TFmJXBHV3uy6evXucB1M72YZGOSiSQT+/btm8NuJUlzcXhbb5TkeuD4Hi9d1O8uevRVn+/9KmAMeP5MY6pqC7AFYGxsrK/9SpLmrrXgqaozZnotyV1JVlTV3iQrgLt7DJsETu9qrwJumO19k5xBJ9yeX1UPzKloSdLQjcqptu3A1Cq1ceDaHmN2AOuTLGsWFaxv+maU5J8Dfwy8pKp6hZkkqWWjEjybgXVJdgPrmjZJxpJcBVBV+4FLgJubn4ubPpJclmQSeEySySRvafZ7OfBY4KNJbk2yvc1JSZJ+Wmun2g6kqu4B1vbonwDO62pvBbb2GHcBcEGP/hlP70mSFsaoHPFIkpYIg0eS1CqDR5LUqpH4jEeSDmXf2vyihS7hkOIRjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVBo8kqVUGjySpVQaPJKlVBo8kqVWp8mab0yXZB3x7oeuYxROA7y50EQvEuS9dS3n+h8Lcn1RVT5xtkMFziEoyUVVjC13HQnDuS3PusLTnv5jm7qk2SVKrDB5JUqsMnkPXloUuYAE596VrKc9/0czdz3gkSa3yiEeS1CqDR5LUKoNHktQqg2eEJTkzyf9KMpHk9iTvS/KEha5rPi3FOfey1H8PS3n+S2HuBs+ISnIOcBkw3nxp7FnAbuCIBS1sHi3FOfey1H8PS3n+S2XurmobQUmOBL4JrK2qryx0PW1YinPuZan/Hpby/JfS3D3iGU2/Aty22P/4plmKc+5lqf8elvL8l8zcDZ7R9HTgy1ONJO9O8uUkNzXtO5Ksa55/NEkWqM5hGsqck/x6s+17kry9lcqHa7bfw6E+v9nMNv9PpOM/J/ntJE9JcnWSaxas4uE54Ny7HerzPnyhC1BPPwKOmmpU1euTbADOS3IC8BngxcBOOqdLK8kmYDnwvaratBBFH6SDnnOSXwROrarXAyR5dOuzOHgH+j30nN8i+LfvNtvfwX3AVcD7q+ovmmHnHqr/AZ7mQHN/Ip3Pfv4LcDHw76rqkJ23RzyjaQfw0iT/BKD5v/t1wBeBk4EbgMckeSrw7SQrgUcB3wNOW5CKD94w5vwbwH+b2mFVPdhW8UN0oN/DbzBtfovk377bbH8HLwA+2hU6i8mMc6+qfcB3gHcCr6+qhxauzIPnEc8Iqqpbk7wZ+EySHwMPARPAB4CNwAfp/E/DG4H/DVwCvAF4InDCghR9kIY05yOAh6f2meSwqvpxa5MYgll+D/+VafNjEfzbd+vj7+BlwIVJvlRVdy1gqUN3oLkneSzwFODhqrpvIescBle1HWKSfBh4BXA0MAk8h86HkkcDxwBfr6o/WLgKh6/fOSd5OvBmYB+dUxa/U1XfW5iqh6/X/IBzWcT/9t2SfBz4NeB44D10/iYeA1xK58jgqqp628JVOD+SHE7nAqFvpTP/m4HbOYTnbfBIklrlZzySpFYZPJKkVhk8kqRWGTySpFYZPJKkVhk8kqRWGTySpFYZPJKkVhk8kqRW/X/35nHWjMI9SgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "logGRatio = (found_parameters - x0) / np.log(10)\n", "x = range(len(logGRatio))\n", "plt.bar(x, logGRatio)\n", "plt.xticks(x, (r'$G_{Na}$', r'$G_{NaC}$', r'$G_{Ca}$', r'$G_{K1}$', r'$G_{x1}$'))\n", "plt.ylabel(r'$\\log(G_{found} / G_{true})$')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the differences are very small.\n", "\n", "We can now compare the true and fitted model output:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEKCAYAAACsUXomAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd4VMX6xz+zJbvplRJIIJRQggkQQBCQIhYEwYYIUhSuFwsWuF4VO3YQbFwUfnAFLxLFhoAioIhICyCh9xZCQk8IIT27e+b3x24aqUCS3cB8nmef3XPOnHO+Z7Z89515Z46QUqJQKBQKhSujc7YAhUKhUCgqQpmVQqFQKFweZVYKhUKhcHmUWSkUCoXC5VFmpVAoFAqXR5mVQqFQKFweZVYKhUKhcHmUWSkUCoXC5VFmpVAoFAqXx+BsAdcKQUFBMiwszNkyFAqFolYRFxeXLKWsU1E5ZVZVRFhYGFu2bHG2DIVCoahVCCESKlNONQMqFAqFwuVRZqVQKBQKl0eZlYvw15pVzJv2KlLTnC1FoVAoXA7VZ+UihPzxFD3FCVYsbs4d9z7ibDkKxXWFxWIhKSmJnJwcZ0u5ZjGbzYSEhGA0Gq9of2VWLkJDkQyA27YvSb9zGN7mK3tDFQrF5ZOUlIS3tzdhYWEIIZwt55pDSklKSgpJSUk0adLkio6hmgFdASlxIw+AHmI7J44fdbIgheL6Iicnh8DAQGVU1YQQgsDAwKuKXJVZuQKWLHRIfnfvh15IdDu+drYiheK6QxlV9XK19avMygWQuekA5ARGsMEWgfe+b0AlWigUCkUByqxcgOyMNAA8vX3Y4NufYO0M2YdWO1eUQqGoUYQQPPfccwXLU6dOZeLEieXuM3PmTObNm1elOsLCwkhOTi63zHvvvVel56wMyqxcgJzMiwDozD7cdNcoLkhP4n+b4WRVCoWiJjGZTCxcuLBCoyjK448/zsiRI6tRVekos7pOyXWYlcHdm/ZN67HI1o1myX+yetsBJytTKBQ1hcFgYMyYMXz88ccltiUkJNCnTx+ioqLo06cPx48fB2DixIlMnToVgGnTphEREUFUVBRDhgxB0zTCw8M5d+4cAJqm0bx58xJmmJKSwu2330779u157LHHkFIWbLvnnnvo0KEDbdq0YdasWQBMmDCB7Oxs2rVrx7Bhw8osV9Wo1HUXIDfLYVZmbzzcDJg7j8IU9xvbfpjMvB1PMOeRTk5WqFBcP7z58x72nrxYpceMaODDGwPaVFhu7NixREVF8cILLxRb/9RTTzFy5Egefvhh5syZwzPPPMOiRYuKlZk0aRLx8fGYTCYuXLiATqdj+PDhxMTEMG7cOFauXEnbtm0JCgoqtt+bb75J9+7def3111m6dGkxs5kzZw4BAQFkZ2fTqVMn7r//fiZNmsT06dPZvn17ueUCAwOvpKrKREVWLkBebjYAJncPAIYM6McqfVeeMCzh4IE95FptzpSnUChqCB8fH0aOHMm0adOKrY+NjeWhhx4CYMSIEaxbt67EvlFRUQwbNoz58+djMNjjkNGjRxf0ac2ZM4dRo0aV2G/NmjUMHz4cgP79++Pv71+wbdq0abRt25YuXbqQmJjIoUOHStVd2XJXg4qsXIBcixUAs7Hw7fC8axK2RbfxkfFzNhzoQ+82oZd9XJsmsdg0zEZ9sfWZuVZOXMjGbNCzJeE8db3NWDSN7s2DMOp1/BCXRJMgDzo0DijYJ8+qMXn5fhoFeODv6caibSd4oW9LWtX3ucKrVihck8pEQNXJuHHjiI6OLtVY8iktDXzp0qWsWbOGJUuW8Pbbb7Nnzx5CQ0OpV68eq1atYtOmTcTExFT6eKtXr2blypXExsbi4eFBr169Sh0nVdlyV4syKxfAarNHTm6Gwrejc/u2nMr8iBtXjmX1z89wKGAOWVZBsK+Zuj7mSh33lZ92seDvRH595mZyrDaiG/nz2+5TjJ+/niCRRiAX8RLZmLBgJo9fRR51vU2cumhhDQb+qutDPX8fGtSrS7LVzIr157goPUjHHYmOQ2fTebJXcxoHeNA21A9PU6H+HIuNNQfPcXub+gCsP5xM40APQvw9qrDmFIprj4CAAAYPHswXX3zB6NGjAejatSsLFixgxIgRxMTE0L1792L7aJpGYmIivXv3pnv37nz99ddkZGTg5+fHo48+yvDhwxkxYgR6vb7E+Xr06EFMTAyvvvoqy5YtIzU1FYC0tDT8/f3x8PBg//79bNy4sWAfo9GIxWLBaDSWW64qUWblCjjGVAl98VbZ4O7D+XlbHANS5rD287t4yfpPkmQdNr/cp8CwMnOt6HWCPScv0qGxP+cz8zifmQuZ5zm85XcG60+x9vMYmopTxIuT9BAp7DFbytaSA7g5Xl9wPOLtiw+YHHKlIBUvzmb4ce5nP05If7bhx3ldAOd1gdRv3JKj1iBWHMnhvXsjOX4+i5l/HcHTTc8Xj3SiS9NApJS0fHU5T9/SnKf7hFdZVSoU1wLPPfcc06dPL1ieNm0ao0ePZsqUKdSpU4e5c+cWK2+z2Rg+fDhpaWlIKRk/fjx+fn4ADBw4kFGjRpUZqb3xxhsMHTqU6OhoevbsSaNGjQDo27cvM2fOJCoqipYtW9KlS5eCfcaMGUNUVBTR0dHMmTOnzHJViSia+aG4cjp27Civ9OaL25bOpv3f/yZ+yF80adWu2LYj5zKY9clEXjfMw4iV37UO7DBE0bRFJPuSbRw9dRZ/0mkszhCmO0NjcYZm4iT+IqPgGLnSyFFZn6MymGyPEA5kmEmRPqTgS7p0p09kIx7t3Zo+0zajSUFUA09ubORNdIgX3206ytETpwnztDIk0pusiyn8vT+eQC5ST1ygjki1P3MBgyg+kDlNepAo65Io6xR7tvg15dG7ejFq3lYAdk28HS9HVCaEYOnOUzQO9OCGhr5XVJ8KxeWyb98+Wrdu7WwZ1cKWLVsYP348a9eudbaUUutZCBEnpexY0b4uHVkJIe6rRLEcKeWv1S6mGpHS/iOv05VsN25Wx4vvtN6szY1kQWQcnQ4upr/cDPlZ7W6FZU/JABJkPX61deaoDOaIbMDgO/swfUs2e89ksfGlPtT3NfPHvjNk5Fqp622mY5g/RkdE9+IQP9yNem6NqFdwzLYdupJrLd7v1RVYvvsUX21O5IP7o/DzMDJ73VG+WPE3jY2pvN3Th7gd25GpCYSKs4SLE/TWbccsHBFdFuR+a+A3t3rEy2AWfvAdu3Pr4hfaBs/gVnwSmwLA0ff6YdUkboay84CsNg2dEKXWnUJxvTNp0iRmzJhRZl9VbcKlIyshRAqwGCjvl6iHlLJZDUkqk6uJrOKWzKDD1gkkDF9P4+Y3lNh+5mIOuRaNRoEerNp3mtf+t5xgkcJT3RuQYjHSsH4wTVtEcOMHGwDY+9YdnEjNJrye91Vd0+WSkpGLXifw87A76IHT6exIuoDVJnnlpx3UIY3RERL3i/Fknz5AU3GKpuIUjcQZ3ERhxuN56cVR2YBDWkMOyRCatenI7b16885fKTzftxUpGXn8vvcM/76jJTe+u5KOYf58PqxDjV6r4triWo6sXIlrNrIClkkpR5dXQAgxv6bEVBvS/kOt05Xs/ASoVySh4pbW9Yl8ZSh1vE0lyn3yYDs8TQY83Aw1blQAgV7FNbWs703L+nYdQzqFYtE0TAY9UkreX7afhBwLPQe2YcCnf5GbfJQm4jRNxUmaiVM0053kdv0Whoo/4eBXcBDelJ4c2BfKIa0hF2UIQ/4KxaqF8OuuXJq+tJTHejYjLNCDwR1D1aSkCsU1hkublZRyeFWUcXWkZo9uhajcsLfSjArgnvYNq0xTVaPTCUwOMxZC8HK/wn9XK567BbiFJTtO8tbPe1j6zM2k51hIkXA2/SzfLV2BdnY/LUQSLXSJDNDH4iuyCvZPlj4c1EI4uC6ErbIxtsRu3Ni5G5+vP8mILo0J8jLh4abH02QokcavUChqBy5tVvkIIR4Alksp04UQrwLRwDtSyq1OllYlyILI6voeoz2wbQMGtm0AFIkm63nz0JARPP3NNnJC/JiTcJ6hnUKZszyWZiTyckfYsXUjLXVJDDWuw6T9BjtnY9sheEI2YN+uRvymNWavbEyCsSn/e/ouVu4/x7zYY4TX9WJY58b0blXXeRetUCgqRa0wK+A1KeX3QojuwB3AVGAG0Nm5sqoG4eg3FKWMgVBAeD1vlo/rUWzdoz3s3ZRSSuIaHKdhRD1MXm6cOX6A+Yt+wZC8lwiRQLTuEAP1sQX7Jf/Hh3CtMQ/JRuxNbcyUg2GsjL6Rm1vWp6GfB99tSeSNAREY9Nf3HweFwtWoLWaV3/veH5ghpVwshJjoRD1VSkFkVclmQEUhQgiGd2lcsFwvrDXjnmlFalYeu06k8dqGYwQZsmmQc4QbPU7icX4vwRf209n6GyaDPTsxd6eBAztC2ak1xSKb8HJSO7p07sahlDxyLRpz1sczY1g0vu5GXl28m1+e7o6HW2356ihqC3q9nsjIyILlRYsWkZyczLx585g2bRqrV6/Gzc2Nrl27OlGl86gt37gTQoj/A24FJgshTFTzvIZCiL7Ap4Ae+K+UclJ1nSu/z+p6bwasKvQ6QZCXid4t69K7ZX4TX89iZdYfOMU785bQQh4jQpdAG3GMAfpYhok/IBlyfzGyT4ayS2vKYH0Tpn2dwCHZECsG9p9Op76PGU+TAR+zQSVzKKoEd3f3YpPDgv3eUh072hPlVq9ejZeXlzIrF2cw0BeYKqW8IIQIBp6vrpMJIfTAZ8BtQBLwtxBiiZRyb7Wc0BFZocyqxujWMpipY4ewIzGN5nW96BTmD1Lywn+XkHVsC5G6o3Q2JTKYWEbYVgKQI43sk43ZOasJ38gm7NCaER3dmVsigrm9TX02x5+nQ2N/9GrMl6KKWL16NVOnTmX69OnMnDkTvV7P/Pnz+c9//sPNN9/sbHk1Sq0wKylllhDiTyBUCBHtWF35O5RdPjcCh6WURwGEEAuAu4FqMStRQeq6onpo08CXNg2KzJIhBJP/eTc2bSCAvd9K09i1ayux61YhTm0jSneU+/VreVj8DkD6bnd27GzKNBnOVi2cnrf0o1e7lni66dl6/AJ9b6hPQkomr/y0m88cTYkKF2fZBDi9q2qPWT8S7iy/cSb/HlEATZo04aeffirYFhYWxuOPP46Xlxf//ve/q1ZbLaFWmJUQ4m3gEeAIkD+KWQK3VNMpGwKJRZaTqMZkDqnlz2BRK96OaxohBAZ9kchIpyOybUci23Zk5l9H2GTVIMyPrFMHaJS9nw1/LaO97jBP6pZgMGiwbgpH1gSzxmFep6J78fbfAg0dK/ee4f4OIc67OIVLU1ozoKKQ2vLrOBhoJqXMq6HzldaOU2KqDyHEGGAMUDD545WQP4uIiqxcm8d7FpkopVkdoDui3RAWbz/J4vQL9PQ+wYbVdvPqpdvOIP0a2PUFD5jM7NCakfB7JCfEPSS430CXliEIAZpENRu6GhVEQArnUFvMajfgB5ytofMlAUVvIBUCnLy0kJRyFjAL7NMtXfHZZP6s6+pHq7bRtI4X429rUbDcpFNfzl7MoePn62kkztLD/RiR8gARcj8PZn+LfskC6ko9R8wtOebZjm/PhtI0ug+DukXQopRZRzJzrRw4k050I/8S2xTXF97e3ly8WLV3MK5N1Bazeh/YJoTYDeTmr5RSDqym8/0NhAshmgAngCHAQ9V0rgKzUpFV7aehnzsN/dxZ9VwvElKyuKlZIFLC15uPM+NgAt5n42ievYPo7L30yvmW2ww2bDs+YM/2ML7QWpFWrzNfnmhIvbr1cHfT427Usyn+fMEkxIrrlwEDBjBo0CAWL158XSZYuPREtvkIIfYA/wfsAgruQyGl/Ksaz9kP+AR76vocKeW75ZW/molsN817lc5H/0PWCyfw8PC6omMoahcTl+zhbEoKTzZPZXfscsIyttFeHMYkLNikYIdsxlotknW2SLbJ5jSq40vS+WxmP9yRnYkX+GePpmrqqCpETWRbM1zLE9nmkyylnFaTJ3TcdqRmbj2iIqvrjokDC2+d7h1xK88u2M5ztzTGfGYbcX/+RA/jbp6yLuJZw09kSDOxFyJYRyRvzU3kiGxAs7pe9IsMduIVKBQ1S20xqzghxPvAEoo3A14jcwM6+qzUOKvrksaBniwa282+0DqETr0GAJBzMYW1qxbRMuNvWhz8g9uM9o97kgzi3Pre/JnQh4WpYdT19+Wxnk2p662aCRXXLrXFrNo7noveL7k6U9drFk1FVoqSmH0C6XXPP4B/EHskhRH/XczNul300m2n+6nFuJ/+nk7SzBotiph93XhyzBMYvOuq7MIrREqpZiOpRq62y6lWmJWUsrezNVQr+c2Aam5ARRnc1CyQ7196CINOcCoth0eXbMWUuI5bdVu5Rb+Nfpmb0T76hAPGVqQ27ssm95t5bGBPNYdhJTGbzaSkpBAYGKgMqxqQUpKSkoLZfOXRv0t/koUQd0kpf7naMi6P1NCkujW7onzyb5sS6GXii0d70Oq1LJLq9GCVvwenD2zkVv1WbpVb6XrkY7ryMdt2hZPdfACh3Ycyc3su90U3pEPjACdfhWsSEhJCUlIS586dc7aUaxaz2UxIyJUPindpswKmCCFOUP5t7d8DardZoaEhMKh/dIpKYjbq+W18Dxr4uaMTEPH6WfTB0bTp3ZzTqfHsWPEldxBL+yMfwZGPuE8LZ+mWLuzsMZxWzcO5qVmgsy/BpTAajTRp0sTZMhTl4OpmdQb4qIIyh2pCSLUiJVq5fqxQlKToIOIV43pQz8eEn4cbUI/QZm2445M1NBan6a/bRH/9Rl43foV1Qwxr1kXxR9uh5DS9g+hm9dHrBHlWjRB/D+ddjEJRAS5tVlLKXs7WUCNIDVm9dzxRXOO0rO9dYnntC705lZZD29BHyLFo9HlrLvfp13Kvfh0Ndr9I2q43+dl2Ez/aerBNNufgO/3IzLXi7+nmpKtQKMrGpc3qekFIm4qsFFVOaIAHoQH2aMlk0PPZs0NIPH83Pef/TWd2c79+Dffr1zLc8Af7tVDeemMlP9m6s+61gcqwFC6HMitXQGrKrBTVTqv6PrSq78Oh9+7CpvUnLftfDPh0BYNMm+mZ/jPv6OYywfANCV8vJ+Pmx/hoh5H374tUM2UoXAJlVq6ApqGpZkBFDaLXCQI83Vj58gBgABbrRNLiN7Ny3rvclbQI04LvGaq15EjAWNr0GQ5qDKDCydSKX0ghhIcQ4jUhxGzHcrgQ4i5n66o6pOqzUjgVo0GPb/hN1Bv5JRNb/MAk2zDqkUqb9c9wfGJLVsx9C5mbwZ6TaeRYbM6Wq7gOqS2/kHOxT7N0k2M5CXjHeXKqGKmhqbR1hQvQPTyI94f1YsLbn7O+3woeyxvHOfy4I+FD0t5rwarPn2H09KVYbRoWm1bxARWKKqK2mFUzKeUHgAVASplN+WOvahVC9VkpXJCHujTl7qGP8/UNX3Bf7kQ2ahGM1S9m7oXRLHl/GMM+/NHZEhXXEbWlzypPCOGO4269QohmFJnQttajUtcVLkq/yGD6RQaz46bG3P1ZC8LEKcbof+EBVtDf8hsbPv2ZtkPfwrNumLOlKq5xassv5BvAciBUCBED/AG84FxJVYjUkCqyUrgwbUP92PLqrYweeCsvW/9Jr9yP+N7Wk47nf8H4WTTrPnqI5KSDpGVbnC1VcY1SKyIrKeXvQoit2GddF8CzUspkJ8uqQqTKBlS4PEFeJkbeFEbi+SzqerfmaNqN9Fy/hScMS3gwbQW62cv51taL/S2fZGD3aG5souYhVFQdtcKshBDRjpenHM+NhBC+QIKU0uokWVWG0FRkpag9vNI/ArDPpP36gAg0bQSjpi2mT/J8hupXYTm8jjkH+tL02akE1anrZLWKa4Xa8nf+c2AjMAuYDcQCC4CDQojbnSmsalAJForaR/6tNHQ6wYeP3om460NGeX3GJrcuPGVYjPvMDhxb8j6PzF7DvlMXnaxWUdupLWZ1DGgvpewopeyA/WaMu4FbgQ+cKaxKkBpS3ctKUYsJ8jIxoktj5j8/lN6vLOGrtvP5O68JYVsn8V7Sw/w4dypn07KcLVNRi6ktv5CtpJR78heklHuxm9dRJ2qqMoRKsFBcYwy/5y729pnDkLxXOSf9eDXvU5I+vJm7X/qU/204xpIdJ50tUVHLqC1mdUAIMUMI0dPx+Bx7E6AJx9ir2o2awUJxbSGE4MlezZk24Wm+iZzLfwOfJ0Qks9j0Ol7LnuKdb1bx4g87nS1TUYuoLb+QjwCHgXHAeOCoY50FqPW3vFeDghXXKnV9zEx6oB2PPv0q+watYq64l7t0sfxp+hf+2z7j+Nk0Z0tU1BJqRTagY8aKDx2PS8moYTlVjpA21WeluObpGdmMnpFfcnj/TuJjxjHBuIDE/27mu4jXMTfpzMC2DZwtUeHC1IpfSMfEtT8IIfYKIY7mP5ytq6oQSNVnpbhuaN4qCo+Hv+Pn1lMx5F5g0LZRJH8/niNJpyreWXHdUivMCvtEtjMAK/Zmv3nAV05VVJVIdYsQxfVFt+ZB3DnoHwxz+5SvbLfyiH4FptndWbdykbOlKVyU2vIL6S6l/AMQUsoEKeVE4BYna6oyhNS4hublVSgqhUGv45fn+zHgpa/5KmIWFqmn69pH+Pb90SyKu2YaThRVRK3oswJyhBA64JAQ4ingBHDNDI3XSSs2UVveCoWi6vBwM+DhBg8/OJi5qyIwrXqVh3J/ZPfiOL47OxWPhm2o72OmY5iauul6p7ZEVuMAD+AZoAMwHBjpVEVViE6zYBVGZ8tQKJzKqFtuYMhbP/Fz66kEixT6bxzGrwtmMGhmLHPXx7M98YKzJSqcSG0xqzApZYaUMklKOUpKeT/Q6GoOKISYKIQ4IYTY7nj0K7LtJSHEYSHEASHEHVetvgL0Wh42nVt1n0ahcHl0OsGAB/9J+qjVHNGF8bnbNF4yxPDOz7u457P1zpancCK1xaxequS6y+VjKWU7x+NXACFEBDAEaAP0BT4XQuir4FxlYpAWbDoVWSkU+YSFNafxv/7kYKMHecywlHnGSfjW/lEqiqvApTtKhBB3Av2AhkKIaUU2+WDPDKwO7gYWSClzgXghxGHgRuyT51YLepmHJlRkpVAUxdfbE9/Rs/jz2wi67n2bhW5v0PPlbGY8dR8RDXycLU9Rw7h6ZHUSiANyHM/5jyVAVTTPPSWE2CmEmCOE8HesawgkFimT5FhXbRikBZtemZVCURq9Bj/L56FTCRJp/GB8nV+WL8Vq05wtS1HDuLRZSSl3SCm/BJpJKf9X5LFQSpla0f5CiJVCiN2lPO7GPm6rGdAO+32y8mfHKC2HXJZx/DFCiC1CiC3nzp27omsEMEoLmuqzUihKRQjB+EdH8WnYZ+Ri5KmEZ3n0tfdJSlWzuF9PuHoz4C4cRpF/75yiSCmjyttfSnlrJc8zG/jFsZgEhBbZHII9wivt+LOw32OLjh07lmpolcEgLUgVWSkU5fL6qPt4erY7jyVOYJbxQ5YuDyJk6GPOlqWoIVzarIC7quvAQohgKWX+/C73Yr8/FtibGL8WQnwENADCgc3VpQPAiDIrhaIyfDjqDg4fb8ORmEEM2D+Bjz9JZdyzL5b6Z1ZxbeHqzYAJ+Q/s/VaRjke2Y93V8IEQYpcQYif2KZzGO865B/gO2AssB8ZKKW1Xea5ycZMWUM2ACkWFuBl0RDRtxLL2M4iTLXgmdRKLvvyQ7Lxq/YoqXACXNqt8hBCDsUc3DwCDgU1CiEFXc0wp5QgpZaSUMkpKObBIlIWU8l0pZTMpZUsp5bKrU18xRqwqslIoLoMHurZmyQ2fEqtFcPexd3j9zQkknld9WNcytcKsgFeATlLKh6WUI7Gnkr/mZE1VgtRsGIUNDCZnS1Eoag2hAR68O7gL0+u/wzrtBiYbZvNbzFRny1JUI7XFrHRSyrNFllOoPdrLJS832/7CoCIrheJy+WjYTRy+ZRbxvjcyKvkjFs2ZxJ8HziLlFec7KVyU2vKDv1wIsUII8YgQ4hFgKfCrkzVVCbm5OQAIvYqsFIrLpYGfO6N7RxD8+ELWapEMTJjEr/Om8MT8rc6WpqhiaoVZSSmfB/4PiALaArOklC86V1XVYMm2t7MLo7uTlSgUtRcPDy9yB33FWi2SyYbZeO9fwMKtSc6WpahCXNqshBDThRBdARwDgf8lpRwvpfzJ2dqqipysiwAIk6eTlSgUtZturUIYJ55nrRbJFOMsNv74CbPWHHG2LEUV4dJmBRwCPhRCHBNCTBZCtHO2oKomN9s+OafBrMxKobgaPE0G4iYOoOfrv7Pa1pYPjLM5vGIGORaV1n4t4NJmJaX8VEp5E9ATOA/MFULsE0K8LoRo4WR5VUJelsOsTF5OVqJQ1H50OgFGM9EvLi0wrN/nT+HwWTVje23Hpc0qH8fA4MlSyvbAQ9hnnNjnZFlVgiUnHQCjhzIrhaKq8PHypt4/f+BvYwcGJLzP7E8nkpyR62xZiqugVpiVEMIohBgghIgBlgEHgfudLKtKsOZkAmA0K7NSKKqS1o3qktxvDn/a2jLZOJs/Y6agaSqlvbbi0mYlhLhNCDEH++SyY7CnqzeTUj4opVzkXHVVgy3H3jxh9vB2shKF4tqjT2QjVkZ+yCpbOx44NYXFc97lkbmbSVFRVq3Dpc0KeBn7TQ9bSykHSCljpJSZzhZVldhy7ZdjUmalUFQ5bgYd7w7uRNiTP7LJ0JF7k6bQ8ch03v1lN3lWdU+s2oRLm5WUsreUcraU8ryztVQXMs9hVp7KrBSK6qJpcBBroz/ma+stPGVYzIA943niv6tUs2AtwqXN6npAc0RWHiqyUiiqlUd7tWJL5Bu8bnuUbrrdTDz5OA+9OpWPfjvgbGmKSqDMytlYsrBIPQY3s7OVKBTXNH4ebnz0YDveevtDkh9YjAU9C9zewX/Na/zfim2cvZij5hR0YZRZORmdNYscoeYFVChqkuA23ZkW/gVL3QfysP43Bm0YwLTJE2jz2i9YbKovyxVx9TsFX/MISxY5mFGNgApFzSGE4JMR3YHuJB/czNGYZ3gromSQAAAgAElEQVTHOJcn5BImvb6chfQmqlkjXu7Xmpb11bfTFVCRlZPRW7PJ1akmQIXCWQS1uJHIl9fxfcuPSJR1ec04nw2GJ7g3/g3em/YfTqekcSErT0VcTkZFVk5Gb80mTyizUiicibvJwAND/8GsNb3YlhTHjReW0evMCu7RbyBz2ids0G5gnXYDPXv3pWOXHvy0M5l3f92HXggGtm0AwORBUQAcPJPO3PXHePqW5jTwU3dTqCqUWTkZo5ZNnl59oBUKV2BMj2ZAM2AwWHNZ8M2XWA78Ri/dDm4zxsG6/5G3Vk8HGcpUEUy8DObo1vokyrp0j4ujb+e2xMSdJdtiI9di49W7IgjwdCM1M48z6Tm0qu9T7HxxCankWTVubBKAXieKbUvJyMXH3YhRX34DWHaeDZNBZ58X0UFCSibJGblk5tro2iwQg+MYh8+m0zjQs8QxUzPz8Pd0Q0rJzqQ02ob6XXklVhPKrJyMUcsmz+BTcUGFQlGzGEwMGfEYudZHufPjNeSkHOe16BwSdq4lQhyjnTjMXbqN6AxFMgi3w7N6D84KX9L2eLJ3rwdmL38OpwkycCexQTBnsiQ2nRuazo2953LJkwaSIxtR19+HPWdyCA7wxs/TzIe/H0JDR+tgX+6ODuXD3w9zxw3BGI0GvNxNnM+0EF7Pl/dXHKRdowBe6teaTfGppGTkMX31YcBuXkFeblg16NmiDou2nSA0wBOTUcehMxnc3KIOBp1g1f5zPHtrCz5ZeRCAf93Wgk5NAmngZ2bB3yf4KvYYF3OsSMDHbKR/VDD/7NGU9BwrDXzdwasuCFGiCqsSoVI1q4aOHTvKLVu2XPZ+x96K5IJHGO3+/XM1qFIoFFWBlBKrJjHqdcQeSWHJjpNM6NsKsy4Pt7QEpnz3Bz62FDoFWWjnl8vyTTvwJpsAfTZmLQsvkY03WXiKa3Sap1fOgPHKujOEEHFSyo4VlVORlZMJu/cNpNnX2TIUCkU5CCEw6u2Rw03NArmpWaBjixHMEbzwdESx8q062+f8lBJ+3nGSuj4m7o8O4envttIpxBOdZuGvvUk83TOUL1YfYP+JZNyw8kT3ENz1GmcvZnMqNYtdSan0aRXE73tOoUNDhyx4Dgtw5/SFDDQtf31h4NGqvhcS6NwkgL8OnOVUWg5WzZ4g0ryOF4fPZSCKlG/ga+ZUWo79WousF0iCfd05lZaNALxMBjJyrcXKAPw7R8PfWBU1XTYqsqoirjSyUigU1zdZeVY2x5+nY1gAXqbi8YOUEiEE5zPz2HY8lRmrj9A9PAgvk4H7okOwaRKTUccXa+Pp0Nif2WuP4utuZPpD0QXHsGkSAXy9+Th+HkbuimrA8ZQsNsWn0D8qmKPnMrmhoS+j5m7mzwPnWPN8bxZuS2LmX0d4Y0Abht7YCCklNk1i0OtYfziZ4+ezaBzgQV0fM1uOnefe6IaYDPoruv7KRlbKrKoIZVYKhaI2k55jIfZICre3qV+j562sWalxVgqFQqHA22yscaO6HJRZKRQKhcLlUWalUCgUCpdH9VlVEUKIc0DCFe4eBCRXoZyqQum6PJSuy0PpunxcVdvV6GospaxTUSFlVi6AEGJLZToYaxql6/JQui4PpevycVVtNaFLNQMqFAqFwuVRZqVQKBQKl0eZlWswy9kCykDpujyUrstD6bp8XFVbtetSfVYKhUKhcHlUZKVQKBQKl0eZlUKhUChcHmVWTkQI0VcIcUAIcVgIMcEF9BwTQuwSQmwXQmxxrAsQQvwuhDjkePavAR1zhBBnhRC7i6wrVYewM81RhzuFENFlH7ladE0UQpxw1Nl2IUS/Ittecug6IIS4oxp1hQoh/hRC7BNC7BFCPOtY79Q6K0eXU+tMCGEWQmwWQuxw6HrTsb6JEGKTo76+FUK4OdabHMuHHdvDaljXl0KI+CL11c6xvsY++47z6YUQ24QQvziWa7a+pJTq4YQHoAeOAE0BN2AHEOFkTceAoEvWfQBMcLyeAEyuAR09gGhgd0U6gH7AMux3musCbKphXROBf5dSNsLxnpqAJo73Wl9NuoKBaMdrb+Cg4/xOrbNydDm1zhzX7eV4bQQ2OerhO2CIY/1M4AnH6yeBmY7XQ4Bvq6m+ytL1JTColPI19tl3nO9fwNfAL47lGq0vFVk5jxuBw1LKo1LKPGABcLeTNZXG3cD/HK//B9xT3SeUUq4BzldSx93APGlnI+AnhAiuQV1lcTewQEqZK6WMBw5jf8+rQ9cpKeVWx+t0YB/QECfXWTm6yqJG6sxx3RmORaPjIYFbgB8c6y+tr/x6/AHoI0TV3xa3HF1lUWOffSFECNAf+K9jWVDD9aXMynk0BBKLLCdR/he5JpDAb0KIOCHEGMe6elLKU2D/8QHqOklbWTpcoR6fcjTDzCnSTOoUXY4ml/bY/5W7TJ1dogucXGeOJq3twFngd+xR3AUppbWUcxfocmxPAwKpBi7VJaXMr693HfX1sRDCdKmuUjRXNZ8ALwCaYzmQGq4vZVbOo7R/Gs4eR9BNShkN3AmMFUL0cLKeyuDsepwBNAPaAaeADx3ra1yXEMIL+BEYJ6W8WF7RUtZVm7ZSdDm9zqSUNillOyAEe/TWupxzO02XEOIG4CWgFdAJCABerEldQoi7gLNSyriiq8s5d7XoUmblPJKA0CLLIcBJJ2kBQEp50vF8FvgJ+5f4TH7TguP5rJPklaXDqfUopTzj+IHRgNkUNlvVqC4hhBG7IcRIKRc6Vju9zkrT5Sp15tByAViNvc/HTwiRf6veoucu0OXY7kvlm4OvVldfR3OqlFLmAnOp+frqBgwUQhzD3l1xC/ZIq0brS5mV8/gbCHdk1Lhh74hc4iwxQghPIYR3/mvgdmC3Q9PDjmIPA4udo7BMHUuAkY7MqC5AWn7TV01wSR/BvdjrLF/XEEdmVBMgHNhcTRoE8AWwT0r5UZFNTq2zsnQ5u86EEHWEEH6O1+7Ardj70/4EBjmKXVpf+fU4CFglHdkDNaBrf5E/HAJ7v1DR+qr291FK+ZKUMkRKGYb9d2qVlHIYNV1fVZUpoh5XlF3TD3uG1BHgFSdraYo9E2sHsCdfD/a25j+AQ47ngBrQ8g325iEL9n9p/yhLB/Ymh88cdbgL6FjDur5ynHen40saXKT8Kw5dB4A7q1FXd+zNLDuB7Y5HP2fXWTm6nFpnQBSwzXH+3cDrRb4Dm7EndnwPmBzrzY7lw47tTWtY1ypHfe0G5lOYMVhjn/0iGntRmA1Yo/WlpltSKBQKhctTrc2AooJBr+UNHhNlDA4s65iXO0BNCDFMFA6y2y6E0EThYLvVjnPkb3NWBpxCoVAoqEazEkLosYeod2If7DdUCBFxSbF/AKlSyubAx8Bkx74R2NtG2wB9gc8dKZ3lHXMy8LGUMhxIdRy7zHNIKWOklO2kPfNmBHBMSrm9iLZh+dulPeFAoVAoFE6iOiOrygx6LWvwWFmDA0s9ZhUMUBuKvT9CoVAoFC6IoeIiV0xpA9Y6l1VGSmkVQuQPHmsIbLxk3/wBZ6Uds9ID1IqcI7nIcR6kpJHOFULYsKfdviMr6NwLCgqSYWFh5RVRKBQKxSXExcUlSynrVFSuOs2qMgPDyipT1vrSIsHyyleoQwjRGciSUu4usn2YlPKEI5X7R+zNhPMuPYiwz/IwBqBRo0Zs2bKllFMpFAqFoiyEEAmVKVedzYCVGbBW1uCxsvYta30yVz5AbQiXNAFKKU84ntOxT9xY6vxkUspZUsqOUsqOdepU+MdAoVAoFFdIdZpVZQa9ljV4rKzBgaUe07HPZQ9QE0LogAew933hWGcQQgQ5XhuBuygchKdwYSw2jcNnMyouqFAoah3VZlaO/qOngBXYR4d/J6XcI4R4Swgx0FHsCyBQCHEY+/TzExz77sE+/fxeYDkwVtqnZyn1mI5jvQj8y3GsQMexyzyHgx5AkpTyaJF1JmCFECJ/IOMJ7FPCKEohLuE8FXTnFcNq0/ju70RsWtWP73t36T5u/egvTqVl8/veM2Tn2cotv+5QMmfTc8otc/BMOu8v23dZ16hQKKoeNSi4iujYsaO8lvusjiVnEuDlho/ZWLBu5d4zPDpvC2/fcwMjujQudb/ot38nKsSXyfdHEZ+cyc6kC7z3634+uD+KwZ1CS93nSrnz07XsO3WRKYOieP6HnTzQIYQ63iZC/D14qHOjYmWllDR56VcaB3rw1/O9yzxml/f+4PTFHP5+5VbqeJvKLKeo3VgsFpKSksjJKf/Pi+LKMZvNhISEYDQai60XQsRJKTtWtH91JlgoriF6TV1NeF0vfv9Xz4J1ialZABw+kw5Adp6NbYmp3BgWQGaeDV93I+cz81h94Bz9p60lOSOP0d2aAHAxx1LlGnWOVJrnf9gJQHxyJt/HJQGUMCuLzf4nLSElq9xjWh0R4O6TabQN8SPA063c8lJKElKyCAvy5Me4JDo09icsyPOyr2XLsfNEhfjhZlDTd9YESUlJeHt7ExYWRsmRLYqrRUpJSkoKSUlJNGnS5IqOob4JikpzqEh/0G97TvP73jNAYWrlyz/t4qHZmxg6eyNt3/ytWNNZckYeAHk2e9OcUV/yo7fnZBpjY7ZisWkltoH9A//aot1sT7xQ6na9rviPjFbk/N9sPl7wOs+qMf7b7VSGfJmj5v7NkFmxJGfkkpCSWWb52KMp9Jq6mq82JvDc9zsYOntjmWXL4sDpdAbNjOX9ZfsqLJtn1RgyK5a4hNTLPo+ikJycHAIDA5VRVRNCCAIDA68qclVmpbgsTqfZP2xjvopjw5EUAFKzLLy7dC+7T6QB8Pcx+w9nfvRSlKzc0s0qLdtC/2nrWLrrVJlmkG2x8dXGBO75bH2p2y/9oTmVVvjFeGnhroLX6w8ns3RX5San1hU55sEzGXSbtIqeU1aXKPdDXBKf/XmYX3baj/vOL3sBuJh9+RFkSmYuAHtOlndLqnxN6Ww8ep5XftpVYVlF+Sijql6utn6VWSkqpGiE1OX9P0ps/3nHSWavjS8WeQGlRkgXc+zjti8Jgth0NKVCHeUlZczfmMCOSyKuomZVlPOZeRWeKx/dJV+wXGvpUd+/v9/BlBUH+HrT8WLl/DwqbjacsmJ/MYMW+UMDHZc7eGYsN5VS72A3cABPk2rRr+0IIXjuuecKlqdOncrEiRPL3WfmzJnMm1diCOhVERYWRnJycrll3nvvvSo9Z2VQZqUoQZOXlvLMN9sKlkuLkCpDVinZeCv32ZsOi/7on0vPJTWr0EDmrj/G38eK36tNSsmi7cWH6a09dI6zF+2G9Oqiyo8uKHquitAV+YZ4m4sbwskL2XSbtIpfdpZ9vzsf9+Kdyd/+fZw3f95TsHwsJYvP/jzC4/O3Fqy7NOlp87HzxYx30rL9hE1YipSSzFy7+Xu46St9TUXJrz+F8zGZTCxcuLBCoyjK448/zsiRI6tRVekos1K4BFLCkh0nSXVEIJdGSFuPV65/JCm17OSFPw+cJccRFYz6cjMv/ljYjBWz6TgPzIzl5x2FJrDhSAqvFTGkXKuNEV9sZtDM2Eppycdq0/hiXXyly+uLRFZmY6EhWG0aR85lcOJCNhN+LLsJTtMkYROWMne9/Zwv/riLueuPFWzPtdrrIL1Iwkl+tCTLuBP4zL+OAJBn0wr+EKw9lIymSb79+zgnL2SXut+ek2k8MT8Oq+P93HQ0hRvf+4NfK9kkqqheDAYDY8aM4eOPPy6xLSEhgT59+hAVFUWfPn04ftwewU+cOJGpU6cCMG3aNCIiIoiKimLIkCFomkZ4eDjnzp0DQNM0mjdvXsIMU1JSuP3222nfvj2PPfZYsT9L99xzDx06dKBNmzbMmjULgAkTJpCdnU27du0YNmxYmeWqGtV2oCiTF37cyeyRHbFeElkNnVW5pIF7P99QYp0bFvxJ59TB4/y+4iQDWvvQ8NQawnU5eIocPMjBjAWDsJL4/QI43Qg0K0Enz/Oe4RxGbOiFxsp3Z/KhUSLTBSz6ng8MSUgEGsLxE29/bcFALkbyMCBX7yYuMYO+GalY9I710kD29mxsejNe3n7g5gluXo5nTwwUGnVGjrXg9YOzNhYkNWTkWimLfBN68+e9NPRzL1hvtWkY9DouZtv3tWkSq03jxIXsgqbES0eV3Pf5ehY+2a1gOTvPVix6/W5LIhMW7qJ5XS9WFsnazGfcgu0cOpvB0eRMWtTz5oAji3PDkWT6RQaXKH+98ubPe9hbif7CyyGigQ9vDGhTYbmxY8cSFRXFCy+8UGz9U089xciRI3n44YeZM2cOzzzzDIsWLSpWZtKkScTHx2Mymbhw4QI6nY7hw4cTExPDuHHjWLlyJW3btiUoKKjYfm+++Sbdu3fn9ddfZ+nSpcXMZs6cOQQEBJCdnU2nTp24//77mTRpEtOnT2f79u3llgsMDLySqioTZVaKYuQVaZ77fe8ZDp1JL9HvopUxNs+fi4SJMzQUyYSIczQUyTQQKQSKiwRwEX+Rgbco8q9/i/3xf2V06+RJPdbNRgwGE0E5Gn30eqzosUkd2EAIEELC0aN012chAIFEh3Q8axix4YYFN6yI1YvoDHQ2XnKiRaWc3MFKIMdkJBMzWdJMppuZi3iQdtKTNKMXadKTNOnJBezPaTjW4cl56U16bmE0NuaruILXN0xcQcyjXQoSMKya5KPfD/L56iOF1ZOQys6kwn64rccvFGRggj0CyyxilCmOSPjw2QyOp2TRKNCj2LXkv2/5saKbI8klr4x+OEXN4+Pjw8iRI5k2bRru7oV/bmJjY1m4cCEAI0aMKGFmAFFRUQwbNox77rmHe+6x33Ri9OjR3H333YwbN445c+YwatSoEvutWbOm4Nj9+/fH39+/YNu0adP46aefAEhMTOTQoUOlmlBly10NyqwU/HftUbYeT2Vg24bc1LT4B+yOT9aw7sVbiq3zNRtxz0qkozhIpC6eFiKRlrok6oi0YuUuSE9OyiAyjAHE5dXnvObDeenNeby5IL3IxIxF50GarYgZYCYHN7T8Fupc2P92Xzq8trxM/cf+1Z+uE5ZWcJWSTRN6MWPlXn7aEo8bVruJCStGrLiTi4fIxYMcPhjYjECjBfIy+XL1bnIy0/HAHvl5koMPmTQUKbQWx/Els7gBX4JF6kk2+ZIsfUiWvvYHviRrvqRsPILZrz7NxHlOpgfy+ercEvv/4Bgnls+PRZaz8mykFck2TC8S+f26+xSP92x2SQ0UJ38MV1lJI9crlYmAqpNx48YRHR1dqrHkU1pm3dKlS1mzZg1Llizh7bffZs+ePYSGhlKvXj1WrVrFpk2biImJqfTxVq9ezcqVK4mNjcXDw4NevXqVmnpe2XJXizIrBe8stY/n+XXXadZPKG5MmgSrTeJLBr112+mj30pn6wHqmuxNYFnSxEHZkD9t7TggQ4iXwSTJOpyQQWRi/2cYWceXXSeKG1nhCSrWl1JB9t6ILzaVuz2yof38nSf95VjjVbix6C+44/Xhul0IdJj2j1vWsSstjWBfc5nZhQas+JCFn8jAl0x8Hc+BIp0gkUYQafZnkUYrXSKBpOEmbPYJw4A/HBNjpEovTslATshATjke4RdasU9oJMk6nMaf5XtOF5w3O89WLLMx/ZKB1ufSc0nJzKVVfR/75Tmu77nvd1Dfx0x4PXs95FpKfxNOp+Uw+P9imf+PzsWitLQsCyajrlgfnqLqCAgIYPDgwXzxxReMHj0agK5du7JgwQJGjBhBTEwM3bt3L7aPpmkkJibSu3dvunfvztdff01GRgZ+fn48+uijDB8+nBEjRqDXl3zPevToQUxMDK+++irLli0jNdX+3U5LS8Pf3x8PDw/279/Pxo2Fzf9GoxGLxYLRaCy3XFVSrlkJIS6deLY0zkspH6kaOQpnc6FIppyZXPrpNhG0cDpxpo0YhMZZ6ccGLYItWku2aC05KEMKo6AyCK/rVWBWXZsF0jEsgGl/HKq0puT0khFHUdYeKj17qkmQJ/HJmfh5XNruVz4PztrIorHdeH3xbiw2jS5NA7ivfQgv/Liz1PJWDJzHh/PSbgpIeLlfK977dX8ZZ5D4ksn4m/zwsJxn/dadNBQpBDseISKZTroD+IlMiIf7HGaWK40kyLokyPock/Xw2nWEOmfdCREGTslAYjYdL3aWvp+sISUzj31v9cUmZUEz4M6kNHaSxm+OJsV1h5NZvvs0fW+oD8Cri3YR2dCX5Iw8jp/PImZzAi/d2brguG3f+o12oX4sGtuN8kjLspCSmcu7S/fxTJ9w2ob6lVteUchzzz3H9OnTC5anTZvG6NGjmTJlCnXq1GHu3LnFyttsNoYPH05aWhpSSsaPH4+fn72+Bw4cyKhRo8qM1N544w2GDh1KdHQ0PXv2pFEj+2wvffv2ZebMmURFRdGyZUu6dOlSsM+YMWOIiooiOjqaOXPmlFmuKqkosmoNPFrOdoH9NvOKa4Sz6bkEk8IYwy/cr1+Dj8jm5KlgfrQNZKUtmp2yKfIykkif7NWM/OFRdbxNzBzRAR+zkZl/HSHPqmHUiwpT4xPOF88qvLtdAxZvLztdPJ/mdb2IT86kjlfFc/oF+5r55MF2POhIHhn/7Xbik+1jn25pVRf3MlLD8w0R4KU7W/H+MrtBjelhb4Ir3bAEaXgxMdYK+ADdS5QY1CGEbYcS6RKUQ+Kxg4SIZBqL04SJMzQWZ7hZtxPzxl8ZC4w12fv3jsoGHJINOaiFsG1FKD4yhFTq0W3yqnLHlmXkWnl8fhyxL91CsK878zfaTe+ZPuEAmAwlr72sWUSKcscnazjtSI2PT8lk1XO9KtzneiYjo3CcYr169cjKKvzch4WFsWrVqhL7FB2HtW7dulKPu2PHDtq2bUurVq1K3R4YGMhvv/1WsFw0G3HZsmWl7jN58mQmT55cYbmqpCKzekVK+Vd5BYQQb1ahHkUNcPZiDrtOpBFe17tYenkdLtBo3Uv8ZfoJAfys3cQCa282y1Zceg9Lk0FXqb6OvjfULzCWMTc3LZgId8qgKF78cSeBniZOOFKtn+kTXmrEVXTMF1BsMt3SiArxZWdSGoGOefwuHetUGm8MaENdH3PBcr4BgT0Rwb2MJq8uTQMLyl4awfWPalBgVq/2b13Q3FoZujQNJC4hlZijAmhbYrtA499dfVkTu5nGutM0FacJF0m0E0cYYCxshsmVRo5YGnDAGMIeLYzdsgl7tDDS8ShxzBOp2QR6Fhp7flq9qcj8hJfOZJ+Ra8WgE8WaBFMz8zhyLqPAqIASGaWKmmHSpEnMmDGjzL6q2kS5ZiWl/K6iA1SmjMK1eCJma7G55AxYeUS/gmcNC/FIsjDfdguzrHdxgqu/oaSXyVAwrseoLzS8u9s1ZGDbBqzcd5ZXftrF6ud74eFmoFGAB//+fke5x7x0cteZwzvw+PzCTLvFY7uhSfhghd0oTJWYDNbTpMfHXPrXwdNkKDOyalykL8fX3Uj/yOCCmTYa+rkzulsT5qyPp14RI6wMHm76MudIBJDomLIhHWjNJlvrYts8yKG5OEELXRIthP3RRbePe/WF01Qd1eqzWzZhtxbGLtmU3VoT0nOsxQZM5/dlFa2/SwdU3/DGihKp8oNmbuDIueJTZuV/BjJyrXy68iD/uq0l7m56fttzmiBvE9GN/FFUPRMmTGDChAkVF6wFXFWflZRyYHnbFa5J0QyyRuIMnxo/o73uMH/Y2nOo3UtM2lT2uKF8KjvNVwM/d/Ic/6qNl5iGEILbIupxW0S9gnWDOoRwX/uGDPvvJmKPphDi705Sqj3y8jYZSM+1lphaqHt4UInj6gXU87YbhFc5UxEVPWZZEZi/h7FM4wj2LTQho17HZ8Oii21/tX9rxt8WzrbjFTebAdwX3ZCFW0/g4aYvuO7LJQszO2UzdtqKZwMGkkakLp424hiRuniidYcYqLcPqtakIH1pOGlB7blP588W2ZLUTPvYK1PRqKmU2T8uveHlpUYFYHEY+P82HGP22ngCPE080atZQTr/sUn9S+yz6WgKzet6EViJZlzFtU9FzYA3AYnYb/u+iUvbghS1igWbjzNh4S46Nrb/i71LF8sk42xs6Bib9wxLtS48aAvG/pZXDWajHk9HVOJdQfNdPjqdKIie3hjQhn/Os98nzNNhLJ2bBPD8HS2ZsuIAULYZDe/SmIs5Fv5xcxNsUvLJyuJNjLEv3cIdH6+xH9vNgFGvY9fE24mcaG+/9/cwkpplwctswL+MOf6KNn+VNpO8TifwNhuLTYf04xM3cf+Mwpk3Zo/sWHCN97Szm1VZtyK5PzqEH7cmlbqtIlLwZbXWjtW0A0drnj8XuUF3jHbiMDdePETbi0v4yM1ukhcO+dHP2Jy4n8M56TWMBq1v4ru/S/9sJKRksufkRbo0LX1sTX60mR9dz1kfTwO/sqNNKSUPztpY5gBnxfVHRWZVH7gNGAo8BCwFvilyd15FLWLWWvsNkc9n5vKMfiH/Mv7A31oLns17ijEDerD0572Vnjcvp0i687djuhQkJpTG+NtaUMfbRP/LmCUhv4/IbCw0AIPjh86gE4zt3Zx72zfkjKNfZO0Lvbn5gz+LHcPNoGPcrS0AGHdrC57tE06uVaOVY8xWsK97wfiSfDMpanwdwwL4fe8ZrDZJ21A/3hgQwZs/7y12jlb1vQuis9LMKp/8bWGBHnRoHMAfz/Vk4pI9PNmrebG+rpvDg1gwpguRDX1LHOP5O1oytnfzKzarS3XvP51OKj6s1aJYS5R9oDUa4eIEHXUH6WE+SiuxlzuMW+CHb9DcvOiZ3Rw3fQSxWgTZOYWfldJmoi+KxaZhtWkFfXjn0nN5dkHZt2nJN7dLozbF9UtFfVY27LeVXy6EMGE3rdVCiLeklP+pCYGKqiPXoiHQGJ87kwHG5fxo685Lln+Sh5GOYQEABenM+Xz4QFu+3HCsIPW8RT0vDruIdTYAACAASURBVJ4p/AGZfH8kkSElf1gBPrg/CrBHRI9dMkC1It659waa1/Wia7PCJr78psd802rg504DxxRGoQEeeJsNNK/rVeJYhfvbEwHubteg4Hrzj5nfJ1V0cGSov70/Kj/RoGeL4n142167DX9PN+r7mkk/m1GsT+5SWgV7M/KmxgUDdZvV8eKrf3QGILFItqMQokR00qNFHdYcPFcQxR2b1J+wMgZBLx7bjee+31HsRz7fmIoS7GsusQ7sfWEHZSgHbaF8ndkHsCfedNbt45mQUzQ+vp5bjHaTsX70AbON4cRqdvPaL0PLzBS12mTBPc0qw5VOnqy4dqlwULDDpPpjN6owYBqwsHplKaqD3DwL7xjmMsDyBzOsA5hsHQII3PS6YokCRWng586XozpxxydrSM7I48lezRECbmoWyFf/3955h0dVpQ38905JgYSEhFASSkJTOkhEQSkqIoiCoq6sLqKLurq6upZdwbWAiqIr+tlZsaGoyKooCgiKIFhQmtSA9KV3CEjKlPP9ce9MZiYzk4FkZgKc3/PMkzvnnvLmZDLvPee85actXN25UdDUHZueurRS+WvqpCRy/yVn+JV50nVYLcG/EJc90ieis7QXBnfyXj97dQeenbWWdJ/zqt6t6vFNwW5u6ZHH9+v3ckPXXADv2UlOejK392pGbXOrrn5aEuv2HPWLIBFIos3KYwPbBr1XUXqPOinGOMEsEp+5qr2f/1fr7Fo81L8VN769EDD8vRpn1PQzQAFo4BOnsCL2ks6X7q58tVFwugdRjwOca1lNV+dqulpWc7Hd6HuvqsX37nbMc7Vnvrs9+yh7iHG63WFX7UWlLj8jFoe7bOWulDotck1ZrVbatWvnff/ZZ5+xb98+3n33XV588UXmzp1LQkIC3bp1i6OU8aMiA4sJQFtgBjBKKRV5HgajfV/gBcAKvKGUGhNwPxF4F+gM7AeuVUptNu+NAIZh7K7fpZSaGa5PEckDJgEZwBJgiFKqNNQYIpKLEUNgrSnOAqXUbWZfnYF3gGRgOnC3CszbcBIw4tPl9GltOHre9M5Chts+5Hqbv6ICuDq/ITUTgn8UEmwWMlMSqV0jgX1HS8mtU5OOpnPnfX0MZeKbm2ryX7qy83BRVL5cPD3aApNhmVhClIejd+t69PYx8AD4z5DOHDpWSmZKIrPuKTsvSUu288HN59AmO400n627m87LZf66fWFXdeGoKL1Hn9b1+XTJds5qUuZUK2JEpPjD2Y38lJXdavGzlsxJr+GX2qR7izrMX7ePuqmG4m2WVZOv/t6Dns/MYUeICB0enOZDyW4y+Nx9Pp+7Df+wbPbRzbqK8y0r6GFZ7rU6XO1uwjx3e75zt2exq2VYX68JP22mR4ssWmcbjtUOH7eIvBHTWTe6X9ht1lOB5ORkv+CwYPhX5efnA0ZYo5SUFK2sQjAE+B1oCdzl8wUkgFLK47JfHhGxYjgMXwxsAxaKyFSllO+m/zDgoFKquYgMBp4GrhWR1sBgoA2QDXwjIi3NNqH6fBp4Xik1SUTGmX2/FmoMs68NSqmOQcR/DbgVWIChrPpiKOyTBqUUH/6ylQ9/2co5eRlca53DbbYveN95kZ+iwrwK9UXvKR41sA3fFuyhQ5AtP1/F1CUvoyp/DcAwSDha4mLUVOOoNDB9fVVjtUhIC7RuzeuUK7vwzHpBrdkiJZRpvcdhum/b+qx9oq+fc+73D1xYLu9VMHq3rsuanWXbfY8NbMuQN3/mmvxGdG9Rh2ZZKditlpCm+ZHgqpXDx4V1+NjVE8FNa9lCT8tyeliXM8w6ndtsX3BMJbJ3+tkMsTZnliuf3fh/TsbMWMOYGWvYPKY/z81ay4vfrve7X+xwnfLKKhhz587l2Wef5eWXX2bcuHFYrVYmTpzISy+9RPfu3eMtXkyp6MyqMp+OLsB6pdRGABGZBAwEfJXVQGCkef0x8LIY33wDgUlKqRJgk4isN/sjWJ8iUgBciGEEAjDB7Pe1MGMERUQaALWUUj+Z798FruAkU1a+Dru19y3mCdtbzHO141HnUAKNOgOz4ULZOYenn27N6vidH8WSzk3MLzZTzGgrq1jj+TgO6JDtVz73Hxewy1ztBEaRyAmyjbfm8b5+78/JyyDRZiXFZ2WVV6emNzCxbx81fFbWV3bKYcrS7UFlvTa/EQk2C+8t2OIty2+SwTQzJ5bCwiqVxypXHq+6BlKTIs61rKa7ZQU99y3jcfv3PG5/h6Xu5nzlOpuZ7nw2qzLDG4fLXU5RGeUx3NiYMRx2hc5RdkLUbwf9xoSt4skRBZCXl+eNYg7GCuu2224jJSWF+++/v2plO0mIOJCtiNQGGvm2UUotCd2CHPxtoLcB54Sqo5RyishhINMsXxDQNse8DtZnJnBIKeUMUj/UGAB5IrIUKAQeUkrNN+v7mlv59uWHiNyKsQLzxtOqLngSG6ZxlEcdz7Fd1eEOx904g/zJg6nuTo1rs2bXkXLZceOJBPw8lVg/ul+5h4ac9OSgSikUHgMMz9lWlrnVlxpBynvPyiot2c7TV7UPqayu6JRDy3opvLdgC+c2zSDJbuXM+qleZRXI7yQz292Z2e7OADST7VxiWUhf60JG2D9kBB+y0p3Lf109+dzVjctfCh4yKJyD9KlCsG1ATRkRfROJyOPAjcBGyuJkK4zVTMhmQcoCH49C1QlVHmylF65+uDF2Ao2VUvvNM6rPRKRNhHIbhUq9DrwOkJ+fX63OtAzTcsUz9tepwyEGOUYFDbEDwX/hhy9rxYAO2bTJDm7pF8gZ9VK5olNQnV5leL7Mg9hznPTYqnCLq2OjdJ4a1I7+7Y0VS0oEDxweBdeibkq5CCG+NKydTGZKot+257TlkWca3qByeNWVw6uuK8hhL5dYFzHIOp9R9gk8aHufWfvzeUv6sVS18GsX05xbFayANPEh0sfmPwDNlFKR254aK5JGPu8bAoHRRz11tomIDUgDDlTQNlj5PiBdRGzm6sq3ftAxTIOJEgCl1GIR2YBxNrfNbB9O7mqHw+XmyekF3HFBc+qkJDJ+/kYGWeZziXURY1zXs0I1DdnWExNv5ahLaPvoTMD48uraLPLkaTPv6VG5XyACPAuPUOneNQYiwh+7lK30Q8U19MWzs1rR2VX9tPKOvKGi2lstEtRS1MN2snjL1Y+3XP1oLZu5xvodg6zzuTxxAb+4z+B152V84z4LEK9xB0BhsYM73l9Cm+w0hvcLHpz1VCQ1NZXCwqrNYHwyEenj3ErgeOP7LwRaiEieiCRgGEwEhm+aCgw1r68GvjWVyFRgsIgkmlZ+LYBfQvVptplj9oHZ5+fhxhCRLNMIBBFpao6xUSm1EzgiIueaZ1s3+PRVbZldsJu3f9jME18aR4Kffr+Mh+wTWexuwX8c/UK2e/aaDtzaw1Bkvg6x1dFU2LuyOvV3hKoUz9/y4gCrR19WmWncPZae39zbg6/v6cGjl7dm9JVlJvfBjBxCbRW3za7F0K5NIpJxtcpllHMo3UpeYpRjCNmynzcSxjIl4VG6SAGjp61m/rq9ANz8ziLmr9vHuO82VNDrqcXll1/OlClT6NixI/Pnz4+3ODEn0pXVU8BSEVmJuRqB8LEBzfOhO4GZGGbmbymlVonIY8AipdRU4E3gPdOA4gCG8sGsNxnDGMMJ3GE6KBOsT3PIB4BJIvIEsNTsm1BjAD2Ax0TEiWEef5tS6oB573bKTNdncBIYV3iePB0uhcuteND2AakUMcJxc9iUHld3bhjyXnVj9JXteHJ6Ac3q1oy3KCcdqx+7xJvGPhh7zJxhfz4vD4DmdVMBaFEvlcJiB/+aEtprJVQUfKtFGDWwLUUOF5MXlR0DT7+rO1sPHuPRz1f5RWYH44zrbVc/3nX1YZB1PvfaPmZy4uPMWj+dfxTcyIS/X8kvmw8EDnVK4JsixEOvXr3o1asXAC1btmT58uA51U4HIlVWEzBMvlcQUW5XA6XUdAzTb9+yR3yui4FrQrQdDYyOpE+zfCNlFoO+5UHHUEp9AnwSYuxFGP5lJx0KxcL5X3GNbR6vOAfwm2pUcaOThM5NavPJ7aenj0llqRHCj87D60M6M3XZDq+jsy8VRa33XVm9fdPZ7DxUzINTVmAznbfTAgIEt86uRevsWqQl2xkcIkyXCyv/dfViqqsbf7Z+xV22T/k68Z98PXUzQifvA9iGvUdplnVi/m2ak4tItwH3KaVeVErNUUp953lFVTLNCbH9YBE1vxvFHpXOy84r4i2OJgYM7JhNl9zK+bf1aVOfl687K+i9cCsy8A9QfMEZdb3O0R7DjnsvPoORl7cu1+7cppksH9nHryzQsbqEBF5zDaBP6dP86m7GoB1jmWB/mkyM8F8Xjf2Oz3/djsutvBawmlOTSJXVYhF5SkS6ishZnldUJdMcFx7/0Lo7ZtPOXcDzzqsoovxh+K09mrL2ib7lyjUnLy8M7sTk27pGrf+Kzi8DrQfzm9Tmrgub87QZGzI5wcqN5+Ux/oZ8PrjZ33ulVpKdtjllsQU+ujV4SvStqh5DHCN40DGMcyxrmJ44gi5iJLNcs+sIz85aS9tHZ7J5X/n0JJpTg0i3AT3B1Hw/SRWZrmtiSKnTjQU3D9gmsd6dzWRXL7/7NRKsHCt10TwrJWiaco2mIvq2qR/y3thrOlAz0fhcWSzCvX3OKFcnlIGHr9FGqNQoBsIHrotY6m7OK/YXmJjwJCMct2CzNOcl05F40ZaD5NY5sTPN0yUGYbyobMS6iJSVUuqCSo2iiRrr9xxlytJtNEhLpr9lAc0tO7i99G7uvrgVdquF67o05tdth3j8y9Ws33OUpEqE1dGcvqx5vG/YcEdXVcJQx26ebY3701kRKYsC1YQrSh/nNfv/MTZhHC/M2wNcBQjLth5iQIfssL5iwUhKSmL//v1kZmZqhRUFlFLs37+fpKTjy5jtS0WBbC9TSn1Z2Tqa6DH0rV/YfqiIv/TI5a+2z1nnzuEr99l0q2FniBktvGfLLDad05ixX/9Gtwp8px7q3ypo6gjN6U1SBL5aJ4rdJn5jbB7Tn8+WbufvH4WO5lBkTeFGxwM8od7ibtunpHKMx5xDeG/BFg4XOXjxj51Ctg1Gw4YN2bZtG3v37j3xX0QTlqSkJBo2PPGHmopWVv8Wke2Ej3DzJKCVVZzwhKFptHc+rSxbuaf0dhQWjpX6HzbfeF4eQ7vlep8aE2yWoFEBbu4e2nlYo4kGHqtB312iKzrl8K8pK/jd53N8wRlZrNh+mH1HS0lLtrPvqOIB5y0cJZlhthnYcPGocyhTl+04bmVlt9vJy8urkt9HEx0qUla7gecqqLOugvuaKOIJ6tp6w3i2ksVUdzdSk2xcGiQrr+/2xrJH+pS7r9HEA0/SysD4f78HPHC9dePZPP/1b7z47XrqpCSayRyFx51/woGN22xfICgedt7Ehc/OZfZ9PfWW3ilERVHXe8VIDs0J8OOGfew8XExHWc9ZlvU84hhKcmIiK0ZeUmHbyqSE0GiqEs/KKlxk9dzMGoiUhV1qXjfFZ7taGOM0fP1vs33BPpXGC/uuYndhCSu2H2bGip2M/UMHrbhOcqpPSG3NcXPd+J8BuME2iyMqmU9cPVj46EVxlkqjOT48jsiBqV9e/GMnPlu6nW/X7PHmEfPEGmxazhHYUFiZHOYe+yfsI41znyq7O3JgG2ol2VFKsWHvUW+EDs3Jg1ZWJzG1kmwkFO/nMstPvO/qze8kVxipQKOpbjx46ZnkpCeVM20f0CGbAR2yWbf7CE0yDXP03q3r8Z95G+nTuh5WEZ7/5jefFsII581kSiGP2d5ml6rtTU1SXOqiVpKdKUu3c+/kZUz4cxd6tsyK1a+oqQJOv9SbpxCZKYkMtn5Lgrh419WHl687vkNljaY6kJpk584LW4RMqtmiXqrXFP3s3Aw2j+lP25w07u7dopyJuhMbdzjuYpXK5QX7K7QQIyZhkRnd4tethwDYtLd8HD5N9SYiZSUiVhEZICJ3ici9nle0hdOEJ0HcXGebzTxXOzapBpwXp0y+Gk28CBYKqogkbi29lyISGW8fSxpHvdaxnnOxqswfpokNkf7FvsBIvpgJpPq8NHEk372cbDnAhy4jkMjxOkJqNCc7oT7zu8jkL6X30ED287L9RYpKjKjyHovDiuIdaqofkR5wNFRKtY+qJJrj4r0FWzi3cAYHLSnMdhthGrWy0pxuBFM6bXNqsXJ7IUtUSx5y/pl/219n+89jIHcsHy82tgU9jsiak4dIv91miIh2zKkmdHpsFv/+bAF9LIv5zHUepRhRr20h9vw1mlMV34/8mfWNzZ63bjzbW/ZfVy/ec/Ymp+ANipaX5VCtZJg6TRyIVFktAKaISJGIFIrIERE5ffMrx5mDxxwMsP5Eojj42NXTW679SDSnG77p7m/r2YyCx/qSlZLoV+dx5xAOpbchcdqdNJLdgBH4ucTp4kixI6byak6cSJXVWKArUEMpVUsplaqUqlVRI03VM3PVLgCusX7HKncTVqnc+Aqk0cSR4f3O9F6LGM7ugQ9tpdj5odNYnC7Fq/YXSKSUUpeboW/9QruRs2ItsuYEiVRZrQNWqsrGeNdUiqJSF395bzHNZRsdLBv52NUj3iJpNHElMAtxKD7eaOWOY7fQzrKZh2wTKXW6WbDxAFA+zJOmehKpstoJzBWREcdjui4ifUVkrYisF5HhQe4nishH5v2fRSTX594Is3ytiFxSUZ8ikmf2sc7sMyHcGCJysYgsFpEV5s8Lffqaa47xq/mqG+E8RRVPJtTLrAtwK+FLV/QS7mk0JwO+yirco/R3v+3la3c+/3H2Z4jtG5bPeMN7b8+REu/15n2/s8z0xdJULyJVVpuA2UACEZqui4gVeAXoB7QG/igigbmthwEHlVLNgeeBp822rYHBQBugL/Cq6esVrs+ngeeVUi2Ag2bfIccA9gGXK6XaAUOB9wJku14p1dF87alogmJBsdMFKC6zLOBndyv2ks6Qc5vEWyyNJm40rxsYdik4nqOtfzuv5Rf3GTxlf4PmpsPwzkNF3PT2L+QOn0avZ+cy8JUfoiWuphJEpKyUUqOCvSpo1gVYr5TaqJQqBSYBAwPqDAQmmNcfAxeJseE8EJiklCpRSm0C1pv9Be3TbHOh2Qdmn1eEG0MptVQptcMsXwUkiYj/yWw1o6jUxRmyleaWHXzpNpI2a5sKzelMeo0ERl/ZFoCmWeUzBHdolE5Nn6DNTmz8rfRvHCOR1+wvUINixn23gTlrdR6r6k6kESzmiMi3ga8KmuUAW33ebzPLgtZRSjmBwxiOx6HahirPBA6ZfQSOFWoMX64CliqlSnzK3ja3AB+WOJvZOVxuZq7aRZHDxWXWBbiU8JWrzDy3fcM0UpN0TEDN6cl1XRrzw/ALad8wvdy9z+84j9bZ/rZgu8ngLsffaCo7GGMfz7x1WlGdDET6DXe/z3USxpe7M0RdD8G+4AN3lUPVCVUeTLmGq1+hHCLSBmNr0NeP7Hql1HYRSQU+AYYA7wZ2IiK3ArcCNG7cOMgwVcO4uRsY+/Vv3N6zKddYFvCjuw37SQMgN7Mmowa00WbrmtMWESEnPTnk/do1EsqV/eRuw1jnNfzTPplFjpa8S8VpdTTxJdJtwMU+rx+UUvcC51TQbBvQyOd9Q2BHqDoiYgPSgANh2oYq3wekm30EjhVqDESkITAFuEEptcHn991u/jwCfICx/VgOpdTrSql8pVR+Vlb0IjjvOFwMgGvHMppadjHN3AIEuKRtfa2oNJowZKYYyiqvTk0+vq3MKOk11wC+cXXiIdtEOonOIVvdiXQbMMPnVce0zqtfQbOFQAvTSi8Bw2BiakCdqRjGDQBXA9+a5vFTgcGmJV8e0AL4JVSfZps5Zh+YfX4ebgwRSQemASOUUt4TVRGxiUgd89oOXAasjGSeooUnk2r6pmk4lYWZrnwAZt/XM+wTpUajwZs2Z1CnHPJzM1j2aB96t6qLwsK9jtvZrTJ4OeFFalMW52DOmj3MWXt8dlXr9xxFe/dEj0itARcDi8yfPwH3UWZtFxTzfOhOYCZQAExWSq0SkcdEZIBZ7U0gU0TWA/cCw822q4DJwGrgK+AOpZQrVJ9mXw8A95p9ZZp9hxzD7Kc58HCAiXoiMFNElgO/AtuB8RHOU1SwWy2Aor9lAT+423IQYw/ebtGxADWaYCx7pA/LHjF29j0uH7VMM/e0ZDsN0oyHvEJSuN1xN3U4zAv2V7Bg+Fzd9M5Cbnp7Ybl+C3YWkjt8Glv2/w7AD+v3UbCzkMVbDtL7ue+Y8OPmaP9qpy0RnVkppfJOpHOl1HRgekDZIz7XxcA1IdqOBkZH0qdZvpEg23WhxlBKPQE8EUL0ziHK40KCzUI72UQTyx5edlzhLdfBODWa4KTVKPO/8mQXruFjFejZGgRYqZoy0jmUp+xv8k/1EWOcf/Te+2LZDi5t18Cba8sTCPfr1bu5uXtTrn/DyNb93B86AGX5sjRVT1hlJSIXKqW+FZFBwe4rpT6NjlgaXywC/a0LKFVW7xYgeFZcGo0mHPdc3BKHS3FZ+2xvWUZNf6OLD10X0kY2c5vtCzaoBvzX1QuAv324lO2Hirj5/DxKnO6Qjsf3Tl4GlPlzaaqeilZWPYFvgcuD3FOAVlYxoKTUxR+tC/je3Y5CypwgtbLSaCqmXq0kxporHw8XtarHlv3HePP7TWaJMNI5lMayhydtb7JV1WWB24g3sHL7Yf76/hJmrd7NTeflhh1L66roEVZZKaUeNX/eFBtxNMHIPLyChrKP511X+5XrBHIazYmRk57Mw5e19lFWhsPwnY67+CRhJOPsz3NV6Ug2qBy+XL7TW8ezsgplgasNLKJHpNaAT5rWc573tUUk1HmPpoppdeAbSrHxtdv/KM1m1WdWGk1VUkhN/uy4HwdW3kt4ihz8HYZ/LzHcSwU4fKx8epFQqsrlVvznuw3e9prjJ9JH835KKe/JoVLqIHBpdETS+HL49xLaHZ7LYttZFOIfTkYnW9RoKsdTg9qVK9uq6vF3+yOkUMzEhCfJosxo4kixoWxKnG46PBZ5epFZq3bx1Iw1PP3VmsoL7cPhIgc/bdhfpX1WVyJVVlbfuHkikoxh4q2JMsOeeJlM115+SuruV/7gpWdqZ2CNppKcnZsRtLwoozU3lv6TunKI9xKeog6HASg0kzUu/d/B4B2GWFq9t2CL0b6oapM93jJhEX8cv4Bjpaf+ii1SZTURmC0iw0Tkz8DXlAWH1USRy60/Uazs/FrzPJ4a1I7GGTUAsGhFpdFUmsCswp5/q6ZZKSxRLbnFcR+NZQ+TE0aRw17vyup/B44F7e9IiZNRX6zy+nZ5+NFc/VT1idaybcaqz3UCZohnj/7Gq0RPBiINt/QMhs9TK4y0HY+bZZooMnP5Ni61/sxsdyckMZU/dmnM9Lu7c2O3XK4/R6cG0WgqS1oNO/P+cQFrHu/L9Lu6c0v3pgA0MR8Kf3S35U+lI6gjhfw3cRTpxwyDjDW7jgTtb95ve3n7h8188PP/vGW+RhdVbX9RaiaO3F1YTMuHZkSci8vtVuw9UsLDn8U1OM9xEbE5mVJqhlLqfqXUfUqpmdEUSgPPfLWG9yZNJEsKmerqRpLd+FOlJNoYOaANyT4OjhqN5sRpnFmDJLuV1tm1aJdjBIhu36gsgvsS1ZLBpQ9hx8nLx/5JD8uyCvt0+2ilEmf0MhF7hpm7di+lTjdv/bApfAMTh/vEZVq76whvRzhOVRJWWYnIEREpDPI6IiKF4dpqTpzdhcW8OncDl1t+4ohKZq67I8l2rZw0mmhzeYds5v3jAnq29A9MvVrlcmXpY+xQdXjb/gy3Wr8g0k29otKyLcGpy3ZExSLQo7Qi3Q10uk58idf3hXmM+mJ1zM30wyorpVSqUqpWkFeqUqpWuLaa42fywq3kDp/GOU/OJgEHfa2/MMudTwkJeiWl0cSIxpk1/N6/cUM+E4edwzZVl0GlI5nh7sKD9g95w/6sn6WgL54zZafLzfZDRX73pi4rSz6hlCp3vuXh/Z+3cMf7SyKSWZmK0x1EWx0ucjBmxhoe+Xwl+4+WmHKduKLx6KgTOSerDMflVSoidUWksecVLaFOV974fqP3uqdlGWlyjC9cRkqDJL2y0mjiQu/W9ejWzMjXWkQSdzruYqTjBs63rGRm4j+51LKAwFXWz5v287/9x3h06ioue+l7v3sPf7aS3OHTAHh93kbOfPgrDgexEvzXlJVMW7GzXHkw3N6VlXHx0ux1/OW9RZQ4XTw1vYBx323g3Z+2MHpaAaVON3d/tPR4piAojkoovBMhokC2ZpT0sUA2sAdoghH1vE30RDs9mLJ0G263YdVzzGe74FrrHPaodL53Gym79TagRhM/LH4+jcI7rr7Md7fjOftrvJrwIvNc7XjcOYR1qiEAM1ftZuaq3UH7cpqapdjhYnaBkYZkwcb9JNgsPDRlJbPv6+n3cLr1wDFy0pMDZPDHo6Q8P9/8YROHjjn4aqWRYdzDxn2/8+vWQ8xdGz478gMfLyfRbuGxgW1D1nG43SQTu++lSFdWjwPnAr+ZEdgvAn4I30QTCePmbuS+/y7j3Z+2sO2gsV3QgP1cYPmVj1y9cJrPE1pZaTTViz2JTbiqdCSPOobSwbKBGQnDedI2nkYSXEkFcqTYSRNzy3H7wSJGTl3F9kNF7AjYNuz+zBzeMVOPHDpWGnT7LfDMypPnbu2uI35uLr9uPRTSJ6vU6SZ3+DTe/3kLHy3ayrs/hTdrdzjdrNlVGLOzq0iVlUMptR+wiIhFKTUH6BhFuU4b0n1SGXi41jYHAT4yIz+D3gbUaGLN9ec0pk5K6NgHgzrl4MTGBNcl9Cp5jomu3lxlnc+chPt4zv4q7WUD4YwwjhQ7vHeLHC4cptWgW8F9k/0tDlfvLKTdyJl0fOxrRk8rCigYwgAAEhVJREFUYM+RYkZPW+297zBN2D1nVh6FVup0E+iSue9oqd/7BRv3U1Tq8hp+jJ5W4L3naxwSyPSVu+j7f/P5y3uLveNHk0iV1SERSQHmAe+LyAvAqe8yHQNq1/BPVZCAg8HWOXzvbss2Vddbnq0zAms0MWX0le1Y9FBvn/f+W2K+D5AHqcVI5410L3mBt1z9uMSykKmJDzM94UFutH5FfcqHRDpS7PQqlaJSF6XmGdDK7Yf5ZMk2v7oWKQv19N/FW7l94hLGzy8zHy92eBRdmZICw2w+MIDAoWP+ymrw6wsY/ulyip2GYvI9jrjhrZ/ZeuAYTlMZHfWxZFxhOiTPWr3bz1Q/WkSqrAYCx4B7MDL3biB42hDNcRL41PMH61zqy0HGufyn94z6qTGUSqPRBBLoiB9st2MPtXnSeT3nlrzC2ITbcGJhpP1dFiT9jR8zHufF7K/pLGtJwMFf31/ijUDx8pz17DMt9YJt801eVKa8XG7F4i3+4Z48FoWe8zCPs3CJ00XgUdfBAGUFhoIMtopauPkg3Z+ZwzMz1wL4OR0XOcpWU7HIAFFR8sXmQD2llOd8yg1MEJEeQDoEeVzQHBeb95eFbUmmmL/aPmehuyU/ug3blTGD2uFSimZZNUN1odFoYsSbQ/MZNmERAJ2b1Abg5vPzeMMn1QjAEWpQkHMNLxX0oJls52LLYh6ovZ4B299mQCIUKztLf2/BMtWUAktj1qjGbFTZOLDxb1MxhOJYEKVSYq6K5q/bx5Fih3dbrsTpJsnmr1QP/B4kWrzCzxAjkNkFu3nw0lZMXrTVW/aFjwl+LOKUVmQN+H/Ag0HKj5n39Oqqklzcuh4FOwsBxQjbh2TLAe4uvRMjCQHUT0ui1xl1w/ah0Whiw0Wt6nmve7TM4sfhF5KdnuxVVuNvyOeWdw1ldkb9FL4p2M0GlcMGVw7Db+lP0aE93P3Mq5xjWUMXSwE3Wb4i0WZsrTmUle2qDtuK6rDVVpdtKoudKoMD1GK/qsUBlcp+alEcJIZ4sc8qZ+2uI+wuNFZpJQ43NQJ8NA/+Xn5l5Q7j7wXGam7igi18/uuOkHWiTUXKKlcptTywUCm1SERyK+pcRPoCLwBW4A2l1JiA+4nAu0BnjFXatUqpzea9EcAwwAXc5QnxFKpPEckDJgEZwBJgiFKqtCrHiAZ/T/mGRNtP5Mou+lt/YbzzUhaqM6M1nEajqUICz5Ivbl2mzFrULdu6//CWcwFITq9L70HDaFE3hW/W7uXK2QXkyS5ayf84w/I/GsleGsleeluWkCWHg45ZpBI4ShLHVBIuewr7HTZqbE7jYrtQRCLr3hjPYzYbDmw02JtK8pEksqzHcGDDgRXHGht/slpQCG4EhVDbkUhqwVqusmxCAW4suBEw66SXJvLD1DlcYhEsFgueCFLKfKjG3Rcs0TUCq0hZJYW5F/bEX0SswCvAxcA2YKGITFVKrfapNgw4qJRqLiKDgaeBa0WkNTAYw48rG/hGRFqabUL1+TTwvFJqkoiMM/t+rYrHqHIsSydyq3UNRSTykvMKnnP6ZwPWeUc1mupHHx+lBNC/fQMubuVfllenbOu+q+lUDPCH/EYALPnfIZzYWKcaGv5Z7m5+7T+5uSOdaxdz55tfU3xoDxlSyJUtE1mxbhM1KebaDhkcKTzEjs07yHAdobkcoYaUYMeJDRd2nCQVurAednCBvYJvklJgAYxNCHHfAYS6B+D6Z9yV1UIRuUUpNd63UESGAYsraNsFWK+U2mi2mYRhqOH7pT8QGGlefwy8LMbm50BgklKqBNgkIuvN/gjWp4gUABcC15l1Jpj9vlZVYwTIXXX89UdamN7svnxzb09enbve6zmv0WiqB5vH9C9X9sp1Z5Urq5Vc3i3Fl9Sk8F+/tsSakJnDxsQdrHY3Ir9JbXa1b8yTBYZZ+/XX9GfT/w7yp1d/NJRJGCy4seP0vqwoc02lsOAmOy2R5lkp/LhhDxaz3CYKlNtsr7xrLSueLccyBTjNGk6TVQ0VKau/A1NE5HrKlFM+ho69soK2OcBWn/fbgHNC1VFKOUXkMJBpli8IaJtjXgfrMxM4pJRyBqlfVWOUQ0RuBW4FaNz4xKNPnZOXwc+bDnjfWwSa103huT9oVzaN5mTFbhXG/eks6tYKvkEVaPVXI8HqZzxhNy3sPLmzSl1u6qb69xVptnA3FkpIoCTE8mjXYVhyGKDsfHx43zMZM6N8ZuPzm9fh+/X7/Ast0bcGrCiQ7W6lVDdgFLDZfI1SSnVVSu2qoO9gsxi4Fg1Vp6rKq3KM8oVKva6UyldK5WdlZQWrEhEvBzyV2WLwh9doNNElwWqhb9sGnNW4dtD7gXom0I8rwWZU8Pg2rdl5hLqp/sYV1giV1YkQzCH6zaH53HOxcVrSNie2scwjig1oRqyYc5x9bwMa+bxvCASaknjqbBMRG5AGHKigbbDyfUC6iNjM1ZVv/aoaI2oE+ihUtD2g0WiqP/YKfI8GndWQA787+GrVLpZtPVQupJqnfWqSzXAIFoKsrKL3YNu+oZHbq0PDNJZtM4w9GmfU4OAxY8+xhj2231PRfIRfCLQQkTwRScAwZpgaUGcqMNS8vhr4VhmBpqYCg0Uk0bTyawH8EqpPs80csw/MPj+vyjGqaE6CYrOWPR3945Iz+Ogv50ZzOI1GEwMSbOG/Xu1WC7f3aub1oUywWRgzqJ3ffYDpd3UHjKgUtZL9FURlV1b92tYnzedsbcnDF3uvG9ZO5vsHLmDCn7t4y9JrJHjHrJ+WxDNXt2fisKCnJFVO1FSjeT50JzATwwT8LaXUKhF5DFiklJoKvAm8Zxo3HMBQDJj1JmMYNTiBO5RSLoBgfZpDPgBMEpEngKVm31TxGFHB9wnsjguaR3MojUYTIypaWXl49LI25KQn06NFFjarheGfrvBrX7tm2TmTiHDBGVn0MJNDhjuzat2gFhecmcUrczaErPOXns2on5bE2z9s5sz6qWT4jJVks9Kwtn9ur/QaduqkpPP4FW25omM2qUnhjUiqkqiu45RS04HpAWWP+FwXA9eEaDsaGB1Jn2b5Rsqs+XzLq2yMaGG3Rt/7W6PRxJZI/6/Tati5r88Z5co9K7OaplOvZ0Xz9k1lX3PhVlb39WnJRa3qeZXVB7ecw3Xjfwbgxm65tGqQSoeGacxcZZgfeM6onr+2A+/8sNkvJUnbnFqs3F7oVaBDzvUPPRUL9OFINSAWoUo0Gk1s+PCWc5mxcmel/689Z9kiwthrOtDOPEPyJVz82IsC/L66Ns3krotaMKBDNs3rpnjLPWdlngwQV3ZqyJWdGvq1/ejWrhwIEvkilmhlpdFoNFVI12aZfk7AJ4rvyuyqzg2D1il1BQ+R5ImY4YuIcO/FLcuVJ9kNpZiSGFod1Ey0UTPM/VigbaQ1Go2mGtGynrHqicR4Iq9OCuc3r1Ou/HiUZYkZVzC9RvQdeyuDXllVEx6+rDUdG5Vf5ms0mtOLSbd2Zf2eoxFtI1otwvgb8mn1yFcnPN4BM2VIRs3YGUucCFpZVROGnZ8XbxE0Gk01IKNmAl3yMiKu79nGS7BaeP2GzuwuLD6u8TwWfU0yq3caIq2sNBqN5iRGRHjlurNo1SCVplkp5e6fnVubnYdDK7A7LjB8vQKD81Y3RMUgHfHpQH5+vlq0aFG8xdBoNJqTChFZrJTKr6ieNrDQaDQaTbVHKyuNRqPRVHu0stJoNBpNtUefWVURIrIX2HKCzetgRI6vbmi5jg8t1/Gh5Tp+qqtslZGriVKqwhxLWllVA0RkUSQHjLFGy3V8aLmODy3X8VNdZYuFXHobUKPRaDTVHq2sNBqNRlPt0cqqevB6vAUIgZbr+NByHR9aruOnusoWdbn0mZVGo9Foqj16ZaXRaDSaao9WVnFERPqKyFoRWS8iw6uBPJtFZIWI/Coii8yyDBH5WkTWmT9rx0COt0Rkj4is9CkLKocYvGjO4XIROSvGco0Uke3mnP0qIpf63BthyrVWRC6JolyNRGSOiBSIyCoRudssj+uchZErrnMmIkki8ouILDPlGmWW54nIz+Z8fSQiCWZ5ovl+vXk/N8ZyvSMim3zmq6NZHrPPvjmeVUSWisiX5vvYzpdSSr/i8AKswAagKZAALANax1mmzUCdgLJngOHm9XDg6RjI0QM4C1hZkRzApcAMQIBzgZ9jLNdI4P4gdVubf9NEIM/8W1ujJFcD4CzzOhX4zRw/rnMWRq64zpn5e6eY13bgZ3MeJgODzfJxwO3m9V+Bceb1YOCjKM1XKLneAa4OUj9mn31zvHuBD4AvzfcxnS+9soofXYD1SqmNSqlSYBIwMM4yBWMgMMG8ngBcEe0BlVLzgAMRyjEQeFcZLADSRaRBDOUKxUBgklKqRCm1CViP8TePhlw7lVJLzOsjQAGQQ5znLIxcoYjJnJm/91Hzrd18KeBC4GOzPHC+PPP4MXCRSCVz1h+fXKGI2WdfRBoC/YE3zPdCjOdLK6v4kQNs9Xm/jfD/yLFAAbNEZLGI3GqW1VNK7QTjyweoGyfZQslRHebxTnMb5i2fbdK4yGVuuXTCeCqvNnMWIBfEec7MLa1fgT3A1xiruENKKWeQsb1ymfcPA5XPWx+BXEopz3yNNufreRFJDJQriMxVzf8B/wTc5vtMYjxfWlnFj2BPGvE2zTxPKXUW0A+4Q0R6xFmeSIj3PL4GNAM6AjuBsWZ5zOUSkRTgE+DvSqnCcFWDlEVNtiByxX3OlFIupVRHoCHG6q1VmLHjJpeItAVGAGcCZwMZwAOxlEtELgP2KKUW+xaHGTsqcmllFT+2AY183jcEdsRJFgCUUjvMn3uAKRj/xLs9Wwvmzz1xEi+UHHGdR6XUbvMLxg2Mp2zbKqZyiYgdQyG8r5T61CyO+5wFk6u6zJkpyyFgLsaZT7qIeBLS+o7tlcu8n0bk28GVlauvuZ2qlFIlwNvEfr7OAwaIyGaM44oLMVZaMZ0vrazix0KghWlRk4BxEDk1XsKISE0RSfVcA32AlaZMQ81qQ4HP4yNhSDmmAjeYllHnAoc9W1+xIOCM4EqMOfPINdi0jMoDWgC/REkGAd4ECpRSz/nciuuchZIr3nMmIlkikm5eJwO9Mc7T5gBXm9UC58szj1cD3yrTeiAGcq3xeeAQjHMh3/mK+t9RKTVCKdVQKZWL8T31rVLqemI9X1VlKaJfJ2RdcymGhdQG4F9xlqUphiXWMmCVRx6MvebZwDrzZ0YMZPkQY3vIgfGUNiyUHBhbDq+Yc7gCyI+xXO+Z4y43/0kb+NT/lynXWqBfFOU6H2ObZTnwq/m6NN5zFkauuM4Z0B5Yao6/EnjE53/gFwzDjv8CiWZ5kvl+vXm/aYzl+tacr5XARMosBmP22feRsRdl1oAxnS8dwUKj0Wg01R69DajRaDSaao9WVhqNRqOp9mhlpdFoNJpqj1ZWGo1Go6n2aGWl0Wg0mmqPVlYajUajqfZoZaXRxBERyfRJ/bBL/FNn/BiF8W4Ukb0i8kYV9HWtmQbiy6qQTaMJh63iKhqNJloopfZjxMhDREYCR5VSz0Z52I+UUndWthOl1Ecishu4vwpk0mjColdWGk01RUSOmj97ich3IjJZRH4TkTEicr0YifpWiEgzs16WiHwiIgvN13kRjHGjiHwmIl+IkeDvThG5V4wkewtEJMOsd5eIrDYjf0+K7m+u0ZRHr6w0mpODDhiRwQ8AG4E3lFJdxMi++zfg78ALwPNKqe9FpDEwk+DRxANpi5G+IwkjRM4DSqlOIvI8cANG0NLhQJ5SqsQTv06jiSVaWWk0JwcLlRmkVEQ2ALPM8hXABeZ1b6C1T567WiKSqozEh+GYY9Y5IiKHgS98+m5vXi8H3heRz4DPKv3baDTHiVZWGs3JQYnPtdvnvZuy/2ML0FUpVRSFvvsDPYABwMMi0kaVJd7TaKKOPrPSaE4dZgFewwkR6VgVnYqIBWiklJqDkS02HUipir41mkjRKyuN5tThLuAVEVmO8b89D7itCvq1AhNFJA0jLcXzykgOqNHEDJ0iRKM5jRCRGzHyHlXadN3srxdwv1LqsqroT6MJhd4G1GhOL4qAflXlFAy8ChystFQaTQXolZVGo9Foqj16ZaXRaDSaao9WVhqNRqOp9mhlpdFoNJpqj1ZWGo1Go6n2aGWl0Wg0mmrP/wPH6cbHr2qnlgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Evaluate model at found parameters\n", "found_values = problem.evaluate(found_parameters)\n", "\n", "# Show quality of fit\n", "_, axes = plt.subplots(2, 1)\n", "axes[0].plot(times, values[:, 0], label='Noisy data')\n", "axes[0].plot(times, found_values[:, 0], label='Fit')\n", "axes[1].plot(times, values[:, 1], label='Noisy data')\n", "axes[1].plot(times, found_values[:, 1], label='Fit')\n", "axes[0].set_ylabel('Voltage [ms]')\n", "axes[1].set_ylabel('Calcium [mM]')\n", "axes[1].set_xlabel('Time [ms]')\n", "axes[0].legend()\n", "axes[1].legend()\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.7" } }, "nbformat": 4, "nbformat_minor": 2 }