diff --git a/docs/source/mcmc/differential_evolution_mcmc.rst b/docs/source/mcmc/differential_evolution_mcmc.rst new file mode 100644 index 000000000..18ec0d688 --- /dev/null +++ b/docs/source/mcmc/differential_evolution_mcmc.rst @@ -0,0 +1,9 @@ +*************************** +Differential Evolution MCMC +*************************** + +.. module:: pints + +.. autoclass:: DifferentialEvolutionMCMC + +.. autoclass:: DreamMCMC diff --git a/docs/source/mcmc/index.rst b/docs/source/mcmc/index.rst index 8a688392c..7cd7060b4 100644 --- a/docs/source/mcmc/index.rst +++ b/docs/source/mcmc/index.rst @@ -11,6 +11,6 @@ parameters that maximise a :class:`LogLikelihood`. .. toctree:: adaptive_covariance_mcmc + differential_evolution_mcmc mcmc - diff --git a/examples/inference-differential-evolution-mcmc.ipynb b/examples/inference-differential-evolution-mcmc.ipynb new file mode 100644 index 000000000..ad0b75dd3 --- /dev/null +++ b/examples/inference-differential-evolution-mcmc.ipynb @@ -0,0 +1,155 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Inference: Differential Evolution (DE) MCMC\n", + "\n", + "This example shows you how to perform Bayesian inference on a time series, using [Differential Evolution MCMC](http://pints.readthedocs.io/en/latest/mcmc/differential_evolution_mcmc.html).\n", + "\n", + "It follows on from the [basic MCMC example](./inference-first-example.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running...\n", + "Done!\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAI0CAYAAAAZXPP3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xu8nWV54P3flXMICQSyQSAJoYqI\noqJsgco4VRgRQQWrtjhSKKPmtUNnnNcONXTasdVa49tardXqRKUED3WsFsGChxSPnfFAgJSzEDFg\nDJJEhAA5J9f7x3q2bsJ97+yV7LX23mv/vp/P/qy1rvUc7rX2s8PF89zXc0VmIkmS1MsmjfYAJEmS\nOs2ER5Ik9TwTHkmS1PNMeCRJUs8z4ZEkST3PhEeSJPU8Ex5JktTzTHgkSVLPM+GRJEk9b8poD6Db\n5s2bl4sWLRrtYagH3HjjjRszs2+0xyFJ2rsJl/AsWrSIlStXjvYw1AMi4r7RHoMkaXi8pCVJknqe\nCY8kSep5JjySJKnnmfBIkqSeZ8IjSZJ63oSr0tK+W7Tk2mJ8zdJzujwSSZLa4xkeSZLU80x4JElS\nzzPhkSRJPc+ER5Ik9TwTHkmS1PNMeCRJUs+zLH2Cs9RckjQReIZHkiT1PBMeSZLU87ykNUHULl2N\n1PKSJI1lnuGRJEk9z4RHkiT1vI4lPBGxICK+ERF3RsTtEfHWJn5IRKyIiHuax7lNPCLigxGxOiJu\niYjnD9rWRc3y90TERYPiJ0XErc06H4yI6NTnkSRJ41cnz/DsBP4gM48HTgUuiYhnAkuA6zPzWOD6\n5jXAy4Fjm5/FwEeglSAB7wBOAU4G3jGQJDXLLB603lkd/DySJGmc6ljCk5kPZOZNzfNHgTuBo4Bz\ngeXNYsuB85rn5wJXZsv3gIMj4gjgZcCKzHwoM38BrADOat6bk5nfzcwErhy0LUmSpF/qyhyeiFgE\nPA/4PnB4Zj4AraQIOKxZ7CjgJ4NWW9vEhoqvLcRL+18cESsjYuWGDRv29+NIkqRxpuMJT0QcCHwB\n+G+ZuWmoRQux3If4k4OZyzKzPzP7+/r69jZkSZLUYzqa8ETEVFrJzqcz85+a8IPN5Siax/VNfC2w\nYNDq84F1e4nPL8QlSZKeoJNVWgF8ArgzM/960FvXAAOVVhcBVw+KX9hUa50KPNJc8voqcGZEzG0m\nK58JfLV579GIOLXZ14WDtiVJkvRLnbzT8mnA7wC3RsSqJvZHwFLgcxHxRuB+4HXNe9cBZwOrgc3A\nxQCZ+VBEvAu4oVnunZn5UPP894ArgJnAl5sfSZKkJ4hWgdPE0d/fnytXrhztYXRdt1tFTIRu6xFx\nY2b2j/Y4JEl7552WJUlSzzPhkSRJPc+ER5Ik9TwTHkmS1PNMeCRJUs8z4ZEkST3PhEeSJPU8Ex5J\nktTzTHgkSVLPM+GRJEk9r5O9tDRKut1GQpKksc4zPJIkqeeZ8EiSpJ5nwiNJknqeCY8kSep5JjyS\nJKnnmfBIkqSeZ8IjSZJ6ngmPJEnqecNKeCLihE4PRJIkqVOGe4bnoxHxg4j4zxFxcEdHJEmSNMKG\nlfBk5r8D3gAsAFZGxGci4qVDrRMRl0fE+oi4bVDsTyPipxGxqvk5e9B7l0XE6oj4YUS8bFD8rCa2\nOiKWDIofExHfj4h7IuJ/R8S0Nj63JEmaQIY9hycz7wH+GHg78BvAByPiroj4zcoqVwBnFeLvz8wT\nm5/rACLimcD5wLOadf4uIiZHxGTgw8DLgWcCr2+WBXhvs61jgV8AbxzuZ5EkSRPLcOfwPCci3g/c\nCZwOvDIzj2+ev7+0TmZ+G3homOM4F/hsZm7LzB8Dq4GTm5/VmXlvZm4HPgucGxHR7PvzzfrLgfOG\nuS9JkjTBDPcMz4eAm4DnZuYlmXkTQGauo3XWpx2/HxG3NJe85jaxo4CfDFpmbROrxQ8FHs7MnXvE\nJUmSnmS4Cc/ZwGcycwtAREyKiAMAMvOTbezvI8BTgROBB4D3NfEoLJv7EC+KiMURsTIiVm7YsKGN\n4UqSpF4w3ITnX4CZg14f0MTakpkPZuauzNwNfIzWJStonaFZMGjR+cC6IeIbgYMjYsoe8dp+l2Vm\nf2b29/X1tTtsSZI0zg034ZmRmY8NvGieH9DuziLiiEEvXw0MVHBdA5wfEdMj4hjgWOAHwA3AsU1F\n1jRaE5uvycwEvgG8tln/IuDqdscjSZImhil7XwSAxyPi+QNzdyLiJGDLUCtExD8ALwbmRcRa4B3A\niyPiRFqXn9YA/w9AZt4eEZ8D7gB2Apdk5q5mO78PfBWYDFyembc3u3g78NmI+HPgZuATw/wskiRp\nghluwvPfgH+MiIHLRkcAvz3UCpn5+kK4mpRk5ruBdxfi1wHXFeL38qtLYpIkSVXDSngy84aIeAZw\nHK0Jw3dl5o6OjkySJGmEDPcMD8ALgEXNOs+LCDLzyo6MSuPeoiXXFuNrlp7T5ZFIkjTMhCciPkmr\nnHwVsKsJJ2DCI0mSxrzhnuHpB57ZVEdJkiSNK8MtS78NeEonByJJktQpwz3DMw+4IyJ+AGwbCGbm\nqzoyKkmSpBE03ITnTzs5CEmSpE4abln6tyLiaODYzPyXpo/W5M4OTXtTq4SSJElPNKw5PBHxZuDz\nwP9qQkcBX+zUoCRJkkbScCctXwKcBmwCyMx7gMM6NShJkqSRNNyEZ1tmbh940XQpt0RdkiSNC8NN\neL4VEX8EzIyIlwL/CHypc8OSJEkaOcNNeJYAG4BbaXU4vw74404NSpIkaSQNt0prN/Cx5keSJGlc\nGW4vrR9TmLOTmb824iOSJEkaYe300howA3gdcMjID0eSJGnkDWsOT2b+fNDPTzPzA8DpHR6bJEnS\niBjuJa3nD3o5idYZn9kdGZEkSdIIG+4lrfcNer4TWAP81oiPRpIkqQOGW6X1kk4PRJIkqVOGe0nr\nbUO9n5l/PTLDkSRJGnnDvfFgP/B7tJqGHgW8BXgmrXk8xbk8EXF5RKyPiNsGxQ6JiBURcU/zOLeJ\nR0R8MCJWR8Qtg+cMRcRFzfL3RMRFg+InRcStzTofjIho98NLkqSJYbgJzzzg+Zn5B5n5B8BJwPzM\n/LPM/LPKOlcAZ+0RWwJcn5nHAtc3rwFeDhzb/CwGPgKtBAl4B3AKcDLwjoEkqVlm8aD19tyXJEkS\nMPyEZyGwfdDr7cCioVbIzG8DD+0RPhdY3jxfDpw3KH5ltnwPODgijgBeBqzIzIcy8xfACuCs5r05\nmfndzEzgykHbkiRJeoLhVml9EvhBRFxF647Lr6aVZLTr8Mx8ACAzH4iIw5r4UcBPBi23ll9dPqvF\n1xbikiRJTzLcKq13R8SXgRc1oYsz8+YRHEdp/k3uQ7y88YjFtC5/sXDhwn0ZnyRJGseGe0kL4ABg\nU2b+DbA2Io7Zh/092FyOonlc38TXAgsGLTcfWLeX+PxCvCgzl2Vmf2b29/X17cOwJUnSeDashCci\n3gG8HbisCU0FPrUP+7sGGKi0ugi4elD8wqZa61TgkebS11eBMyNibjNZ+Uzgq817j0bEqU111oWD\ntiVJkvQEw53D82rgecBNAJm5LiKGbC0REf8AvBiYFxFraVVbLQU+FxFvBO6n1YQU4DrgbGA1sBm4\nuNnPQxHxLuCGZrl3ZubAROjfo1UJNhP4cvMjSZL0JMNNeLZnZkZEAkTErL2tkJmvr7x1RmHZBC6p\nbOdy4PJCfCVwwt7GIUmSNNyE53MR8b9olYu/GfhPwMc6Nyz1qkVLri3G1yw9p8sjkSRNJMOt0vqr\niHgpsAk4DvifmbmioyOTJEkaIXtNeCJiMq2Jwv+B1o3/JEmSxpW9Vmll5i5gc0Qc1IXxSJIkjbjh\nzuHZCtwaESuAxweCmflfOzIqSZKkETTchOfa5kejpDbZV5Ik7d2QCU9ELMzM+zNz+VDLSZIkjWV7\nm8PzxYEnEfGFDo9FkiSpI/aW8Axu0vlrnRyIJElSp+wt4cnKc0mSpHFjb5OWnxsRm2id6ZnZPKd5\nnZk5p6OjkyRJGgFDJjyZOblbA5EkSeqUvd54UJIkabwz4ZEkST3PhEeSJPU8Ex5JktTzTHgkSVLP\nG24vLXXJRO2ZVfrca5aeMwojkST1Is/wSJKknmfCI0mSep4JjyRJ6nmjkvBExJqIuDUiVkXEyiZ2\nSESsiIh7mse5TTwi4oMRsToibomI5w/azkXN8vdExEWj8VkkSdLYN5pneF6SmSdmZn/zeglwfWYe\nC1zfvAZ4OXBs87MY+Ai0EiTgHcApwMnAOwaSJEmSpMHG0iWtc4HlzfPlwHmD4ldmy/eAgyPiCOBl\nwIrMfCgzfwGsAM7q9qAlSdLYN1oJTwJfi4gbI2JxEzs8Mx8AaB4Pa+JHAT8ZtO7aJlaLS5IkPcFo\n3YfntMxcFxGHASsi4q4hlo1CLIeIP3kDraRqMcDChQvbHaskSRrnRuUMT2auax7XA1fRmoPzYHOp\niuZxfbP4WmDBoNXnA+uGiJf2tywz+zOzv6+vbyQ/iiRJGge6nvBExKyImD3wHDgTuA24BhiotLoI\nuLp5fg1wYVOtdSrwSHPJ66vAmRExt5msfGYTkyRJeoLRuKR1OHBVRAzs/zOZ+ZWIuAH4XES8Ebgf\neF2z/HXA2cBqYDNwMUBmPhQR7wJuaJZ7Z2Y+1L2PIUmSxouuJzyZeS/w3EL858AZhXgCl1S2dTlw\n+UiPUZIk9ZaxVJYuSZLUEXZL15hV6xxvF3VJUrs8wyNJknqeCY8kSep5JjySJKnnmfBIkqSe56Tl\nUVKbkCtJkkaeZ3gkSVLPM+GRJEk9z4RHkiT1PBMeSZLU85y0rHGnNOHbuy9LkobiGR5JktTzPMPT\nBZagS5I0ujzDI0mSep4JjyRJ6nle0lJPqF02dDKzJAk8wyNJkiYAEx5JktTzTHgkSVLPcw7PCLL8\nXJKksckzPJIkqeeN+4QnIs6KiB9GxOqIWDLa45EkSWPPuL6kFRGTgQ8DLwXWAjdExDWZeUen9+3l\nq/HBcnVJEozzhAc4GVidmfcCRMRngXOBjic8Gt9sQCpJE8t4T3iOAn4y6PVa4JQ9F4qIxcDi5uVj\nEfHDLoytZB6wcZT23Uk98bnivcXwUJ/t6I4NRpI0osZ7whOFWD4pkLkMWNb54QwtIlZmZv9oj2Ok\n9erngt7+bJI0kYz3SctrgQWDXs8H1o3SWCRJ0hg13hOeG4BjI+KYiJgGnA9cM8pjkiRJY8y4vqSV\nmTsj4veBrwKTgcsz8/ZRHtZQRv2yWof06ueC3v5skjRhROaTprxIkiT1lPF+SUuSJGmvTHgkSVLP\nM+GRJEk9z4RHkiT1PBMeSZLU80x4JElSzzPhkSRJPc+ER5Ik9TwTHkmS1PNMeCRJUs8z4ZEkST3P\nhEeSJPU8Ex5JktTzTHgkSVLPM+GRJEk9z4RHkiT1PBMeSZLU80x4JElSzzPhkSRJPc+ER5Ik9TwT\nHkmS1PNMeCRJUs8z4ZEkST1vymgPoNvmzZuXRx+9qKP7yEo8Rmn5oYzktiaa++5bw8aNG7v+Vc2b\nNy8XLVrU7d2qR914440bM7NvtMchddqES3iOPnoR/+f7Kzu6j8xyGhFR/m9jp5cfyu7d5W1NmmTK\nszenndI/KvtdtGgRK1d29hjWxBER9432GKRu8JKWJEnqeSY8kiSp55nwSJKknmfCI0mSet6Em7Tc\nDTt3lScCT51Sngjc7mTjfZmcXOPkZEnSRGDCI2ncWLTk2ifF1iw9ZxRGImm88ZKWJEnqeSY8kiSp\n55nwSJKknmfCI0mSep6TljtgyuT2Kp92Vdo71NRaS0yZXM9fa+tUwm1Xb41kuwtJkkaaCY+kca1U\nuQVWb0l6oq5d0oqINRFxa0SsioiVTeyQiFgREfc0j3ObeETEByNidUTcEhHPr2zzmxHxw2abqyLi\nsG59HkmSNH50ew7PSzLzxMwcaDO9BLg+M48Frm9eA7wcOLb5WQx8ZIhtvqHZ5omZub5TA5ckSePX\naE9aPhdY3jxfDpw3KH5ltnwPODgijhiNAUqSpPGvmwlPAl+LiBsjYnETOzwzHwBoHgcuSR0F/GTQ\numubWMnfN5ez/iQqM2QjYnFErIyIlRs2btj/TyJJksaVbk5aPi0z1zXzbFZExF1DLFtKXEplQG/I\nzJ9GxGzgC8DvAFc+acXMZcAygJNO6m+vJGof1CqTtu3YVYzXqqsmVyqlHt9W3k6t4grKXx5AtRhr\nd2X5ygpWY0mSxrKuneHJzHXN43rgKuBk4MGBS1XN48AcnLXAgkGrzwfWFbb50+bxUeAzzTYlSZKe\noCsJT0TMas7CEBGzgDOB24BrgIuaxS4Crm6eXwNc2FRrnQo8MnDpa9A2p0TEvOb5VOAVzTYlSZKe\noFuXtA4Hrmoue0wBPpOZX4mIG4DPRcQbgfuB1zXLXwecDawGNgMXD2woIlZl5onAdOCrTbIzGfgX\n4GNd+jySJGkc6UrCk5n3As8txH8OnFGIJ3BJZVsnNo+PAyeN7EglSVIvGu2ydEmSpI6ztUSj3V5Q\nteWh3htrZyW+O8slUdt3luMzp00uxh/evKM6pmlTyrnt9KnleK16q7znekWZRkZzK4fFAAsXLhzl\n0UjS+OMZHmkcyMxlmdmfmf19fX2jPRxJGndMeCRJUs8z4ZEkST3PhEeSJPU8Ex5JktTzJlyVVlKu\nsGq3F1Sl4GpIM6aWa5xqBU5bK723Ht+2s63tQL1CbFulsmvOAVPL29lRrhzbsr081rmzplXHVOvL\ntXNXeR+TKr+j2nYkSRrgGR5JktTzTHgkSVLPM+GRJEk9b8LN4ZE0MSxacu2TYmuWnjMKI5E0FnTt\nDE9ErImIWyNiVUSsbGKHRMSKiLineZzbxCMiPhgRqyPiloh4fmWbJzXbXN0s7+xVSZL0JN0+w/OS\nzNw46PUS4PrMXBoRS5rXbwdeDhzb/JwCfKR53NNHaPUX+h5wHXAW8OWhBhC0V5FV65lV65cFsLu2\nTqX6qFZB9fi2cuXTgdPL1V6PbClXbwFMnVz+zLXKsa2VqqtaRdTmyvIzp5XjADt2lT/3rMrnq32v\nkzDPlSQNbbTn8JwLLG+eLwfOGxS/Mlu+BxwcEUcMXrF5PSczv5utrOTKQetLkiT9UjcTngS+FhE3\nNp2fAQ7PzAcAmsfDmvhRwE8Grbu2iQ12VBMfahlJkqSuXtI6LTPXRcRhwIqIuGuIZUvXKPa8njGc\nZVoLthKsxQALFi4czlglSVIP6doZnsxc1zyuB64CTgYeHLhU1TyubxZfCywYtPp8YN0em1zbxIda\nZmDfyzKzPzP7++b17e9HkSRJ40xXEp6ImBURsweeA2cCtwHXABc1i10EXN08vwa4sKnWOhV4ZODS\n14Dm9aMRcWpTnXXhoPUlSZJ+qVuXtA4Hrmqqo6YAn8nMr0TEDcDnIuKNwP3A65rlrwPOBlYDm4GL\nBzYUEasy88Tm5e8BVwAzaVVnDVmhNZKGat9UKSbi/p9vLsZrlVKbtpT7XG3ZXv61bd9ZrgIDWLdl\nezH+8LZy/ITDDyrGa5VV06e0nztPqVSO1arobJklSdpXw0p4ImIOcBmty0ZfzszPDHrv7zLzPw+1\nfmbeCzy3EP85cEYhnsAllW2dOOj5SuCE4XwGSZI0cQ33f8v/ntYk4S8A50fEFyJievPeqR0ZmSRJ\n0ggZbsLz1MxckplfzMxXATcBX4+IQzs4NkmSpBEx3Dk80yNiUmbuBsjMd0fEWuDbwIEdG52kCanU\nB6uT27XHltT7hpvwfAk4HfiXgUBmLo+IB4G/7cTAJKlbTISk3jeshCcz/7AS/wqtflcTTq3PFdR7\nYy049IBifOOj24rxpz2lfPLs7gceK8Y/tvInxTjAW04u33Dx4BlTi/FaNVatt9hBB5S3s3VHvXJs\ncqXsaivl77ZWCVar9pIkacBo99KSJEnqOBMeSZLU84ad8ETEpIh4YScHI0mS1AnDTniaCq33dXAs\nkiRJHdHuJa2vRcRronbvf0mSpDGo3V5abwNmAbsiYgutuy9nZs4Z8ZGNEbW+WDOnlftfAeyo9LTa\nXdlYrVppxd0PFuNfvvPnxfj6R7ZUx/Suf7m7GD/+yPKv7tBZ5UNj1U82VfdR8gcv+rXqe7UqtEc2\nl3uIzZw9vRivVY5ViuWq3/dYFhGLgcUACxeWK+4kSXVtneHJzNmZOSkzp2bmnOZ1zyY70liRmcsy\nsz8z+/v6+kZ7OJI07rSV8ETLBRHxJ83rBRFxchvrT46ImyPin5vXp0fETRFxW0Qsj4gpTXxuRFwV\nEbdExA8iotggNCKuiIgfR8Sq5ufE0nKSJGlia3cOz98Bvw78x+b1Y8CH21j/rcCd0Kr6ApYD52fm\nCcB9wEXNcn8ErMrM5wAXAn8zxDYvzcwTm59VbYxFkiRNEO0mPKdk5iXAVoDM/AUwbTgrRsR84Bzg\n403oUGBbZg5MLlkBvKZ5/kzg+mYfdwGLIuLwNscqSZIEtJ/w7IiIyUACREQfUO8d8EQfAP5w0PIb\ngakR0d+8fi2woHn+b8BvNvs4GTgamF/Z7rubS1/vj4jyrFZJkjShtVul9UHgKuCwiHg3rSTlT/a2\nUkS8AlifmTdGxIuhVdoVEecDA4nK14CdzSpLgb+JiFXArcDNg94b7DLgZ7TOMi0D3g68s7D/X1a4\nLGizwqVS6MOUISp9tre5re2Vqq6TjpxbjN/7ULka62vfuqc6pi2Pl9f59YtPLcYPmFrOhQ85sJxT\n1iqfalVuAD97eGsxXus59ujW0iEAB84oH8b7Uo1Vq/jyTgySNL61lfBk5qcj4kbgDFol6edl5p3D\nWPU04FURcTYwA5gTEZ/KzAuAFwFExJnA05v9bAIubuIB/Lj52XM8DzRPt0XE3wP/vTLuZbQSIk46\nqX+I/wRLkqRe1G6V1icz867M/HBmfigz74yIT+5tvcy8LDPnZ+Yi4Hzg65l5QUQc1mx3Oq2zMx9t\nXh8cEQNzg94EfLtJgvYczxHNYwDnAbe183kkSdLE0O4cnmcNftHM5zlpP/Z/aUTcCdwCfCkzv97E\njwduj4i7gJfTqu4a2Od1EXFk8/LTEXErrcte84A/34+xSJKkHjWsS1oRcRmtUvGZEbGJ1uUsaE1X\nWdbODjPzm8A3m+eXApcWlvkucGxl/bMHPT+9nX1LGj2Lllz7pNiapeeMwkgkTUTDSngy8z3AeyLi\nPZl5WYfHJGmCKCVBktQJ7VZp/Y+IuAA4JjPfFRELgCMy8wcdGNuYUKva2VFr1ATs2FWuujpo5tRi\n/Ge7y9VKqx54uBj/9aPK1Vuv/JOXVcd07l9/uxi/6ntri/Fjjip3DFmzrtxL67HHyrVpLzy63nlk\n98PtzR8/9MDyLZ9q9VO7K7+jSUNUb1mNJUm9qd05PB9m/+60LEmS1HXtnuE5JTOfHxE3Q+tOy4Oq\nqSRJksakbt5pWZIkaVS0m/DseaflfwX+YsRHJUmSNIK6dadlSZKkUdPuHB6AB4HvNOvOjIjnZ+ZN\nIzussaNWtTOj0msKYPqU8ns/r1Qy/fzxcnzG5MnF+Oad5Z5Sl16zujqmbVu3FeO3ff66Yvz+k36j\nGN/0ox8W49+5/L8U4zc/WK40A5g1rfz5Ht+2qxifO6u8nc3by8tPmVz+3U0b4sTmUBVckqTxq62E\nJyLeBfwu8CN+1QszAW8AKEmSxqx2z/D8FvDUzKw1BJckSRpz2p20fBtwcCcGIkmS1CntnuF5D3Bz\nRNwG/HJSSGa+akRHJUmSNILaTXiWA++l1Z3c++9IkqRxod2EZ2NmfnBfd9bctHAl8NPMfEVEnA78\nFTANuBF4Y2bujIi5wOXAU4GtwH/KzNsK2zsG+CxwCHAT8DsjPb9ocqVqp9ZjC2D7znIuuKVSTXTY\n7OnF+IHbyr+ez9z2QDH+31/y1OqYLrnvF8X4nLNeWYyvufmOYnzW0U8rxl+y5IvF+HHPXlgd01/8\n5gnV99pR+17nzir3LhuqXdauSv+t2nEgSRof2p3Dc2NEvCcifj0inj/w08b6bwXuBIiISbTOGJ2f\nmScA9wEXNcv9EbAqM58DXAj8TWV77wXen5nHAr8A3tjm55EkSRNAuwnP84BTad1d+X3Nz18NZ8WI\nmA+cA3y8CR0KbMvMu5vXK4DXNM+fCVwPkJl3AYsi4vA9the0yuE/34SWA+e1+XkkSdIE0O6dll+y\nH/v6APCHwOzm9UZgakT0Z+ZK4LXAgua9fwN+E/jXiDgZOBqYT+umhwMOBR7OzIG78K0FjtqP8UmS\npB7V9p2WI+Ic4FnAjIFYZr5zL+u8AlifmTdGxIubdTIizgfeHxHTga8BA8nLUuBvImIVrQnSNw96\n75ebLeyqOAEjIhYDiwEWLKzPKZEkSb2p3TstfxQ4AHgJrUtTrwV+MIxVTwNeFRFn00qU5kTEpzLz\nAuBFzbbPBJ4OkJmbgIubeAA/bn4G2wgcHBFTmrM884F1pZ1n5jJgGcBJJ/XXZxtLkqSe1O4Znhdm\n5nMi4pbM/LOIeB/wT3tbKTMvAy4DaM7w/PfMvCAiDsvM9c0ZnrcD726WORjY3FRcvQn4dpMEDd5m\nRsQ3aCVdn6U14fnqNj/PPhuiSKv63mFzytVYDzy8tRg/8pCZxfhLf+2QYnzH7vqdAt589tOL8Vt/\n+lgxfuLxhxXj9/3s0WL85OP6ivHjDy9/BoAHHt9SjB93+Oxi/MFHyt/TgkMPKMYf2byjGD9wRv2w\nnz613N9rtA0+S7nQs5SS1LZ2Jy0P/Bdnc0QcCewAjtmP/V8aEXcCtwBfysyvN/Hjgdsj4i7g5bSq\nuwCIiOuafUMrSXpbRKymNafnE/sxFmnMysxlmdmfmf19feXkUpJU1+4Zni81Z1/+ktZ9bxL4WDsb\nyMxvAt9snl8KXFpY5rvAsZX1zx70/F7g5Hb2L0mSJp5hJzzNfXOuz8yHgS9ExD8DMzLzkY6NTpIk\naQQM+5JWZu6mdd+dgdfbTHYkSdJ40O4cnq9FxGuayilJkqRxod05PG8DZgE7I2IrrXvhZGbOGfGR\njRE7d5Urn6ZMrueK06eW39u5q1y+Na/SS2tbpSfXgdPKv7YZU+oVRi88stxX6jl95Yqov/7mvcX4\n8UfPLcbPf9YRxfg/3vGz6pg19NgyAAAgAElEQVRe/+zyOmsfKldv9VW+p1rvsmlTyr+Hofpi7a70\n0ppkLy1JGtfavdNy+b+OkiRJY9i+3Gl5Lq0KqsF3Wv72SA5KkiRpJLV7p+U30bonznxgFa1Got+l\n1cRTkiRpTGp30vJbgRcA9zWNRJ8HbBjxUUmSJI2gtu+0nJlbASJiembeBRw38sOSJEkaOe3O4Vnb\n3Gn5i8CKiPgFlYadvWKoaqyaWi+tWqXPlEo8KlVGCw4p94667YH6bZFmViq4rrjxp8X48xYdXIw/\nZfa0Yvxv/u+evV1bXnD0QdUxTa18t+1WV+2sVFbVvtdJQ9xVYXfllzcJq7QmokVLrn1SbM3Sc0Zh\nJJL2V7tVWq9unv5p07jzIOArIz4qSRqjSkkQmAhJY92wEp6ImAG8BXgacCvwicz8VicHJkmSNFKG\ne71mOdBPK9l5OYNaTLQjIiZHxM1NHy4i4vSIuCkibouI5RExpYkfFBFfioh/i4jbI+Liyva+GRE/\njIhVzc9h+zIuSZLU24Z7SeuZmflsgIj4BPCDfdzfW4E7gTlNM9LlwBmZeXdEvBO4CPgEcAlwR2a+\nMiL6gB9GxKczc3thm2/IzJX7OB5JkjQBDPcMz46BJ5m5c192FBHzgXOAjzehQ4FtmXl383oF8JqB\n3QCzm55dBwIPAfu0X0mSpOGe4XluRGxqngcws3ndTi+tDwB/CAy0p9gITI2I/uYMzWuBBc17HwKu\noVUBNhv47aZbe8nfR8Qu4AvAn2fWaqS6p1YEVKvzmVqpSirXQ9U39Kyn1H8NtcqnPz+rfFeBbTvK\nX/dPHyn3uXr6IbPK8UqvLoAdlT5lOyo9x6rVW5Xv44Bp5cq0ofpiWY0lSb1pWGd4MnNyZs5pfmZn\n5pRBz/ea7ETEK4D1mXnjoG0mcD7w/oj4AfAovzqL8zJad3I+EjgR+FBElPbzhuZS24uan9+p7H9x\nRKyMiJUbNnqfREmSJpr2bzKzb04DXhURa4DPAqdHxKcy87uZ+aLMPBn4NnBPs/zFwD9ly2rgx8Az\n9txoZv60eXwU+AxwcmnnmbksM/szs79vXt9IfzZJkjTGdSXhyczLMnN+Zi6idVbn65l5wUBVVURM\nB94OfLRZ5X7gjOa9w2ndzfnewduMiCkRMa95PhV4BXBbFz6OJEkaZ9rulj7CLm0ud00CPpKZX2/i\n7wKuiIhbac1YeXtmbgSIiFWZeSIwHfhqk+xMBv4F+FjXP4GkJ6jdmE+SRlPXE57M/Cbwzeb5pcCl\nhWXWAWdW1j+xeXwcOKlT45QkSb1jtM/w9KQYoldTSe4uVytVipWYObVSfTSrvt/Ht5Wr+ittqNhV\neePoSh+vLdt3tbUdgNkzyodf7XPPmFq+Alvrd9bmr2FIpeK/US8HlCQNW7cmLUuSJI0aEx5JktTz\nTHgkSVLPM+GRJEk9z4RHkiT1PKu0OqDWzqtWvVWrMmp3+0NVRB18wNRifHKlr9RDj+8oxqdVGlcd\nML18KO2s9Msaap3plZ5ZtUqw2j4mRa2qqzqkqtLvrptdtyJiMbAYYOHChV3csyT1Bs/wSOPAE9qj\n9NkeRZLaZcIjSZJ6npe0JO0z20hIGi9MeCRpBJSSvzVLzxmFkUgq8ZKWJEnqeZ7h6YB2e2mN1Pan\nVCquhrJtZ7nC6ZBZ5aquhx7bXozPqlRcRdRLoqa22QNrxrT2yqtq38ZQ1Wy1qjVJ0vjW1TM8ETE5\nIm6OiH9uXp8eETdFxG0RsTwipjTxgyLiSxHxbxFxe0RcXNneSRFxa0SsjogPRqczDUmSNC51+5LW\nW4E7ASJiErAcOD8zTwDuAy5qlrsEuCMznwu8GHhfREwrbO8jtO5Ncmzzc1ZHRy9JksalriU8ETEf\nOAf4eBM6FNiWmXc3r1cAr2meJzC7OWNzIPAQsHOP7R0BzMnM72brTnxXAud19lNIkqTxqJtneD4A\n/CEwMGlkIzA1Ivqb168FFjTPPwQcD6wDbgXempl7TjY5Clg76PXaJvYkEbE4IlZGxMoNGzfs9weR\nJEnjS1cSnoh4BbA+M28ciDVnZc4H3h8RPwAe5VdncV4GrAKOBE4EPhQRc/bcbGFXxdmoT7hL7Tzv\nUitJ0kTTrSqt04BXRcTZwAxgTkR8KjMvAF4EEBFnAk9vlr8YWNokRasj4sfAM4AfDNrmWmD+oNfz\naZ0R6nm7K1VGQ1Ux1fpv1WzeVu5bdVClJ9ekynzxnUNURNXUKsdmTC1/vlovrdqup1Z6dUmSeldX\n/uXPzMsyc35mLqJ1VufrmXlBRBwGEBHTgbcDH21WuR84o3nvcOA44N49tvkA8GhEnNrM9bkQuLob\nn0eSJI0vo/2/updGxJ3ALcCXMvPrTfxdwAsj4lbgeuDtmbkRICJWDVr/92hNgl4N/Aj4ctdGLkmS\nxo2u33gwM78JfLN5filwaWGZdcCZlfVPHPR8JXBCJ8YpSfur1mvMlhNS93mnZUl7ZZNQSePdaF/S\nkiRJ6jjP8IxD+9JAo9p/a3I5Xqu6mtRmr6lp+9Cbasak9qqxav2vap+5VuUG7X8+SdL44BkeSZLU\n8zzDI0ld5mRmqfs8wyNJknqeCY8kSep5XtKS9EuWn0vqVSY841Ct+mgktVutVOvVNVQLr3b3MWXy\nyJyQtBJLkiYeEx5pgvJsjqSJxIRHksaIUhJq5ZY0Mpy0LEmSep5neCRpDPOePdLIMOGRepxzdXpT\nO79XkyMJolZd06siYgNwXxd2NQ/Y2IX9jBUT8fPOysy+buwsIhYDi5uXxwE/7PAux/rv0/Htn8Hj\nO7pbx7E0miZcwtMtEbEyM/tHexzd4uftLWP98zm+/TPWxyd1gpOWJUlSzzPhkSRJPc+Ep3OWjfYA\nuszP21vG+udzfPtnrI9PGnHO4ZEkST3PMzySJKnnmfBIkqSeZ8IjSZJ6ngmPJEnqeSY8kiSp55nw\nSJKknmfCI0mSep4JjyRJ6nkmPJIkqeeZ8EiSpJ5nwiNJknqeCY8kSep5JjySJKnnmfBIkqSeZ8Ij\nSZJ6ngmPJEnqeSY8kiSp55nwSJKknmfCI0mSep4JjyRJ6nkmPJIkqeeZ8EiSpJ5nwiNJknrelNEe\nQLfNmzcvjz560WgPQz3gvvvWsHHjxuj2fj2G909W4l3/RQ5DN8Z60003bszMvhHc5F7NmzcvFy1a\n1M1dqofdeOPwjuEJl/AcffQi/s/3V472MNQDTjulf1T26zG8fzLLaUTE2Et5ujHWmVPjvhHb2DAt\nWrSIlSs9hjUyIoZ3DHtJS5Ik9TwTHkmS1PO6kvBExOURsT4ibhsUOyQiVkTEPc3j3Mq6uyJiVfNz\nzaD4MRHx/Wb9/x0R07rxWSRJ0vjTrTM8VwBn7RFbAlyfmccC1zevS7Zk5onNz6sGxd8LvL9Z/xfA\nG0d4zJIkqUd0JeHJzG8DD+0RPhdY3jxfDpw33O1Fa8be6cDn92V9SRNXRBR/xqLxNFZprBvNKq3D\nM/MBgMx8ICIOqyw3IyJWAjuBpZn5ReBQ4OHM3NkssxY4quMjliRphC1acm0xvmbpOV0eSW8bD2Xp\nCzNzXUT8GvD1iLgV2FRYrnbLCiJiMbAYYMHChZ0ZpdRBHsOStH9Gs0rrwYg4AqB5XF9aKDPXNY/3\nAt8EngdsBA6OiIGEbT6wrrajzFyWmf2Z2d83r6v315JGhMewJO2f0Ux4rgEuap5fBFy95wIRMTci\npjfP5wGnAXdk625c3wBeO9T6kiRJ0L2y9H8AvgscFxFrI+KNwFLgpRFxD/DS5jUR0R8RH29WPR5Y\nGRH/RivBWZqZdzTvvR14W0SspjWn5xPd+CySJGn86cocnsx8feWtMwrLrgTe1Dz/v8CzK9u8Fzh5\npMbYC3bu2l19b8rkcm5bW6e2fLv7Hmo74+kW/+qO2jEBI3dcbNm+qxifOrm8/Xb/FiSNTf4lS5Kk\nnmfCI0mSep4JjySp4yJicUSsjIiVGzZsGO3haAIy4ZEkddwTbq3Q560V1H0mPJIkqeeNhzstT1jt\nVj5t3VGv0jqwsk671VvbKvuYNaP9Q6ndqpvdu8sVPJMmWdW1L2rfJ4zcd7ptR7kiavrUycX4UMfE\njp3lY++RLTuK8Xmzp5f3PaV8zI/mcWTFotR5nuGRJEk9z4RHkiT1PBMeSZLU80x4JElSzzPhkSRJ\nPc8qrS7avG1nMX7A9PZ+DbXKl2mV6hOAXZWKnIc3lytcZk4tb+uA6eXqmpFUG+vkShVNrXpn6hDf\nR81EqgQbyc+0tXJMUi8EK3psa/lvBGBKpdfVQTOnFuPtVojVbK8cX1Dvy3XAtPI+asek1VhS53mG\nR5Ik9TwTHkmS1PNMeCRJUs8z4ZEkST2vKwlPRFweEesj4rZBsUMiYkVE3NM8zi2sd2JEfDcibo+I\nWyLitwe9d0VE/DgiVjU/J3bjs0iSpPGnW1VaVwAfAq4cFFsCXJ+ZSyNiSfP67Xustxm4MDPviYgj\ngRsj4quZ+XDz/qWZ+flODrzdiiGo98Wpabdn1rad5UqWHTsrlTLA7Q9sKsaP7TuwGD/gwGnF+NZK\nVUrte9pR+WxQryprt2qtVvnS7vcKvVmNtS/fQ7vbmlbZ1vbK8usf2VqMHzhET7bVP3usGH/KwTOK\n8Vql1I8eLG+ndjwuOPSA6phqh8u+fLfSnhYtubYYX7P0nC6PpDd05a8yM78NPLRH+FxgefN8OXBe\nYb27M/Oe5vk6YD3Q18GhSpKkHjSa/xtyeGY+ANA8HjbUwhFxMjAN+NGg8LubS13vj4hya2RJkjTh\njYvzrhFxBPBJ4OLMHDhHfhnwDOAFwCE8+XLY4PUXR8TKiFi5YeOGjo9XGmkew5K0f0Yz4XmwSWQG\nEpr1pYUiYg5wLfDHmfm9gXhmPpAt24C/B06u7Sgzl2Vmf2b2983zipjGH49hSdo/o5nwXANc1Dy/\nCLh6zwUiYhpwFXBlZv7jHu8NJEtBa/7PbXuuL0mSBF2q0oqIfwBeDMyLiLXAO4ClwOci4o3A/cDr\nmmX7gbdk5puA3wL+PXBoRPxus7nfzcxVwKcjog8IYBXwlk6MvVaNVatWGUqtdqu2j3b7Ew3Vj+f4\np8wpxh+t9C764bpHi/HZM8uHTO0zzBiib9GkynhrPcdq+6j1RqoUjk04Q1UU1tR6iu2sxGttyzZv\nKx/DW3eU/3627dy+98HtoVYUed/GzcV435zydL99+Vuv/c3VjuFaFWDt76RW8WnvrfGrVnmlzutK\nwpOZr6+8dUZh2ZXAm5rnnwI+Vdnm6SM2QEmS1NPGxaRlSZKk/WHCI0mSep4JjyRJ6nkmPJIkqed1\nq5dWz6lVGAE8XqnQqPWI2lLpT7VpS3k723eWq0Zq+wX4+eZy9cvhs8t9iGpjqvWt+sfb1hXjv3/a\nMdUx/XjD4+UxHVQe0/TKvidPqvR3qizfbq8zGN9VMbWx1yqxoP773/DotmK89vewacuOYvzQ2eVK\nqVp1E8DDW8vbqlVd1T73VbeXj9WTj3hS/2Jg6L5YBx0wtRifW4lPrWyr9n3PmNr+MTyej1Wpk4Z1\nhicinhIRH4mID0fEoRHxpxFxa0R8buB+OJIkSWPVcC9pXQHcAfwE+AawBTgH+A7w0Y6MTJIkaYQM\nN+E5PDP/NjOXAgdn5nsz8/7M/Fvg6A6OT5Ikab8NN+EZvNyV+7gNSZKkUTHcSctXR8SBmflYZv7x\nQDAingbc3ZmhSZJ6RUQsBhYDLFy4cJRH03m2kBh7hpXwZOb/rMRXA68d0RGNMbU+OkNVbsycVu6L\ns7VSibG7UnFRq9DYURnTUL2jnnlEuZfW266+vRj/rROfUozfu+mxYvylx8wrxlfc/WB1TE+fO7sY\nr/UQq1Vd1fo7Qfl7GqrCbjwXuLRbfVY77qD+PdR6PtWq9z5fqYg666mHFeNf+dH66phevPDQYvy7\n9/+8GP/e/ZuK8f/6wkXF+E0//UUxfurR5f1CvVfY9srf6ORd5S+29m9GL/XSysxlwDKA/v5+O92N\nsFqCtWbpOV0eydjl5ShJktTzvA+PJEnjiJfL9s2wz/BExKSIeGEnByNJktQJw054MnM38L4OjkWS\nJKkj2p3D87WIeE2MxxlzkiRpwmp3Ds/bgFnArojYAgSQmVkuAWpExOXAK4D1mXlCEzsE+N/AImAN\n8FuZ+aQyiYi4CBgohf/zzFzexE+idQfomcB1wFtzX5ok7UWtGqtWvQX1XleTJ5XzxAcfLvcnqn2c\nm3/2cDG+pVLdBHDXLVuK8bdWel0tu+H+Yryv0gNpd6XByNzp5Z5CUK9MeWxruZ/S7Bnlw3XHzvL3\ntKvyfc+sVBrB+Kh+SYbug7WnWjXW1h31Y3hXZfuTKt/prsrfQ/8RBxXjDz6+tRg/dFb9n6TLb/5p\nMb7q7g3F+HGLyr2x7l5frjS89s6NxfgLFhxSHVOtYrL2L9HjlaquxyvH/MGzppW3v7v+u6sZqrJU\nmgja+gvIzNmZOSkzp2bmnOb1kMlO4wrgrD1iS4DrM/NY4Prm9RM0SdE7gFOAk4F3RMTAv2IfoXVP\nh2Obnz23L0mSBLSZ8ETLBRHxJ83rBRFx8t7Wy8xvAw/tET4XWN48Xw6cV1j1ZcCKzHyoOfuzAjir\naVg6JzO/25zVubKyviRJUttzeP4O+HXgPzavHwM+vI/7PjwzHwBoHkt3IjuKVsPSAWub2FHN8z3j\nkiRJT9JuwnNKZl4CbAVozrqULzKPjNKEgRwiXt5IxOKIWBkRKzdsLF/vl8aywcfwRo9hSWpbuwnP\njoiYTJNcREQftfv3792DzaUpmsfSPeXXAgsGvZ4PrGvi8wvxosxclpn9mdnfN69vH4crjZ7Bx/A8\nj2FJalu7VVofBK4CDouId9Pqo/Un+7jva4CLgKXN49WFZb4K/MWgicpnApdl5kMR8WhEnAp8H7gQ\n+Nt9HMeQ2qmGGVCrhqhVXU2dXK582bBpezF+6IzySbXNU+pVWtt2bS7G11R6Y51+bLnCpfZtrFz7\naDG+8OByVRfA8YeX57sfWKnGqlW/Ta/1d6r8HsZBIdaQgnq1VMnuXeXf2pTKcQdDVGlVVvnRhvLx\n9cj2HcX4P99Rroh6xTPLPdkAHn6sXM34utPKjSgPO7BcIbhtV/nv5KxnlHtmPfx4+TMAPLatXF1V\nq8jsm1P+e5g1vXzMV3tpTbLiSmpXWwlPZn46Im4EzqD17+55mXnn3taLiH8AXgzMi4i1tCqvlgKf\ni4g3AvcDr2uW7QfekplvahKbdwE3NJt6Z2YOTH7+PX5Vlv7l5keSJOlJ2kp4IuKTmfk7wF2FWFVm\nvr7y1hmFZVcCbxr0+nLg8spyJwxz6JIkaQJr97zoswa/aObznDRyw5EkSRp5w0p4IuKyiHgUeE5E\nbGrmzzxKa6Jxae6NJEnSmDGshCcz35OZs4G/HHSH5dmZeWhmXtbhMUqSJO2Xdqu0/kdEXAAck5nv\niogFwBGZ+YMOjG1MmzREqU+td1Gtj06tCqRW7XXQtHL1ySe+/5NiHOBVzynd17FeLfPmFywoxj/w\nnR8X47Nqfa4qFUIAGzaVq24OObBchbazsq0jK/uuVdjVa25g2pSxX8KVlKt3av2bdla+h9rxCHDQ\nzPJ3unVH+fs5vFJ9tOLH5eOrVpX0mZXVu0vQd9CMYvz6O8r3JTr/BeUGb3/8j7cV4xecXu4rt3D2\nrOqYamqVg5u2lKu6alVxc2bWe9FJak+7c3g+zMjdaVmSJKkr2j3Dc0pmPj8ibobWnZYjopN3WpYk\nSdpvo3mnZUmSpK5oN+HZ807L/wr8xYiPSpIkaQR15U7LkiRJo6ndOTwADwLfadadGRHPz8ybRnZY\nY0etGGuo6qOtO8rVL49X+u4cOrs8DapWxfTzreUeWy8+7pDqmN5/7T3FeK1Xz8vf9/Fi/KnnnFuM\n/7sTjyzGz3l6vTfSg5u3FuO1ypRZ08uVL1t3lK+q1pYfqsJuPAggCp+h9rvcVTlWD5lVrwCqVRPV\njuEt28vH/IUnzi/Gl9+8thi/7js/q47pwbvvrb5X8p2P312Mn/mfLyrG1z9art+bMbV+Ivwnm8o9\nxH5tarmyq93efO32NIPysSGp/dYS7wJ+F/gRv+ojmcDpIzssSZLGh0VLrh3tIWgY2j3D81vAUzOz\nfIpBkiSNGaVkbM3Sc0ZhJKOv3UnLtwEHd2IgkiRJndLuGZ73ADdHxG3ALyeYZOarRnRUkiRJI6jd\nhGc58F7gVrz/jiRJGifaTXg2ZuYHOzKSMapW8RDUqy1qNRJzZ5WrsX72cLla6bEd5YqYHbvLuebd\nG8rbAZh/xOxi/IIXlKurpv/m/yjGL/nI94rxFz+1fKVz6656v6ZjDy2PqaZWGTd3Vrkaa3KllKVX\nq1hqH6v2PdQqgKDeD67Wr6u2j2+tKfe5evGiucX4Ya95VnVMS/6qXNn1slf2F+NvPvW3i/H7Nm0p\nxv/1Rw8X40NV9T28rVzZVataq1Ug1r6/Wnznrvr/b06Z3P460kTQbsJzY0S8B7iGJ17S2uey9Ih4\nK/BmWnnCxzLzA3u8fynwhkHjPR7oy8yHImIN8CiwC9iZmeV/+SRJoyoiFgOLARYuXDjKo9FE1G7C\n87zm8dRBsX0uS4+IE2glOycD24GvRMS1mfnLm8Zk5l8Cf9ks/0rg/83MhwZt5iWZWW7JLEkaEzJz\nGbAMoL+/v70bEkkjoN07Lb9khPd/PPC9zNwMEBHfAl4N/H+V5V8P/MMIj0GSJPW4tu+0HBHnAM8C\nZgzEMvOd+7j/24B3R8ShwBbgbGBlZb8HAGcBvz8onMDXIiKB/9X8H0Rp3V+eSl3gqVSNQx7DkrR/\n2roPT0R8FPht4L/QmnPzOuDofd1504frvcAK4CvAvwHlmbrwSuD/7HE567TMfD7wcuCSiPj3lf0s\ny8z+zOzvm9e3r8OVRo3HsCTtn3bP8LwwM58TEbdk5p9FxPuAf9qfAWTmJ4BPAETEXwDlUgw4nz0u\nZ2XmuuZxfURcRWsu0Lf3cRxtLT91Sj1XrG2pVhUzt9LTaHqlh8+PNj5ejP+Hp5YrXwDef98vivE/\n+PgNxfj73vSCYvyVZzy9GN9eqQA5ena5p9BQat9HrcKlVhFTCTNrev2wrx0H46GyqzbGGdPKVWxD\nVWnVvqPan0mtquuejeXKwSMPnFGMH3fIgdUxfeCys4rx9/3TXcX4NQfNLMbPfHq559ztq39ejG87\n7ZjqmE6Zf2gxPqVSXTW1UkFVqwQbycqq8d5DTtpf7d5peeBfr80RcSSwA6j/azAMEXFY87gQ+E0K\nc3Qi4iDgN4CrB8VmRcTsgefAmbQukUmSJD1Bu2d4vhQRB9OqmrqJ1smMj+3nGL7QzOHZAVySmb+I\niLcAZOZHm2VeDXwtMwef2jgcuKr5v9opwGcy8yv7ORZJktSDhp3wRMQk4PrMfJhWkvLPwIzMfGR/\nBpCZLyrEPrrH6yuAK/aI3Qs8d3/2LUmSJoZhX9LKzN3A+wa93ra/yY4kSVI3tHtJ62sR8Rrgn7Ld\nWb6SJGnULVpybTG+Zuk5XR5Jd7Wb8LwNmAXsjIittErTMzPnjPjIuqxW4VLL67bWSoCo97LZUamK\n2VzZVq2KZuHcA4rxGZWqLoDfPrncM2vdM+YV4ycdVa5kmT21XCn1yPbtxXitegegUsgyRJeysp2V\n72lfekjV1ulFQ33SWuVb7bubMrl87F30vKOK8Vpfua/e87PqmP51dfmE8hv+w1OL8TU/L/fMesoB\n5QqxT7/5lGK8VuUGMGdme9VsOyv94HZnuRrrgEq13L4cpZMm0LEtlbR7p+X2uj1KktQjamdGND7s\ny52W5wLH8sQ7Le/TvW8kSZK6oa2EJyLeBLwVmA+sotVE9LvsY/NQSZKkbmj3DM9bgRfQavj5koh4\nBvBnIz8sSZJGh5euelPbd1rOzK0AETE9M+8Cjhv5YUmSJI2cds/wrG3utPxFYEVE/AJYN/LDGjtq\nVSlD9dKqFUPUKoCmVbY1s1Idcu+D5V5a164t9wICeOVxTynGN20p92r9xr3ri/GXPu3wYrzWpmeo\nHlQHVD7fjkr/oG07y/FZ08vbqfUOmmjVKrVjeKg+TdMrx2RtnVq10rRaxWJlO888pF7w+awXHlSM\nf+v+jcX4kpc8rRh/+PFyRWGtdu+OB+u3GzvtmHKV49bK55s9o/w91X5Hm7eV/z5nV/rKSaprt0rr\n1c3TP42IbwAH0epyLkmSNGYNK+GJiBnAW4CnAbcCn8jMb3VyYJIkSSNluHN4lgP9tJKdlzOoxYQk\nSdJYN9xLWs/MzGcDRMQngB90bkiSJEkja7hneHYMPMnM8iw6SZKkMWq4Z3ieGxGbmucBzGxe90wv\nrZpaj6ChbK9UE9WqsWZMba/KaFFfuZdW35zp1THVqmJqY3pOX7ki5u71jxbjD20rV74cd8gQ3UgO\nLPdTmlXpH7RtR61/Wfn7mzmt/d/d7kq1zFiq7ErKPd5qFXG16sAY4v93aj3QDqxUGT22tfz/QQ8+\nsq26j5JaZSLAe76xuhj/8GueXYxvfLR8TP54U7nK8dmHlY/55x55cHVMkytVaFH538KtlWO41jNr\nqGpQSe0ZVsKTmfV/hSRJksa4Uf/fh4h4a0TcFhG3R8R/K7z/4oh4JCJWNT//c9B7Z0XEDyNidUQs\n6e7IJUnSeNF289CRFBEnAG8GTga2A1+JiGsz8549Fv1OZr5ij3UnAx8GXgqsBW6IiGsy844uDF2S\nJI0jo32G53hafbk2N5OhvwX/f3t3HidXVeZ//PMlidkgyBJAlqAOEFQGkEQwwyIIOAM4KBFHFFHZ\nIhiGn+MoorjwE1GCrx+CMhijjKiEIAxLEDDA8BsERZSAQNhhUNlBQhwM2ZPv/HFOkaKpW+lOqvvW\nvXner1de6b51u+s5VdW3Tp1znvNwyCp+pmFX4FHbj9leAlwMvK+f4gwhhBBChZXd4bkX2EvSRpJG\nAAcCW7U4b4KkuyX9Qv0JLU0AABm0SURBVNLb8rEtgCeaznkyH3sNSZMkzZY0+88v/LmT8YcwIJpf\nwy/EaziEEPqs1Ckt2w9ImgLcAMwH7gZ65jfcCWxte76kA0l1vLYlZYi95lcW3M80YBrAuHHji0rm\ndMzggqyYokyfolpQywoyhua93DoFpCiDBmD+ota/a7vN1m15/Kl5C1seL8qi2XGL1hkuw9pk3Sxb\n3jqmoqyioUNat68ok65VJlM6XhhSV2VjNev5Gm6VkdXX9rZr64qC56YoE2xoQabhxgWZg0U1pYqe\ne4DP7vU3LY/PX9w686no72rnzVpnXW1UkDXYrh7c8oLHqehvsSj7rfD3F1wD2j1OIYTWyh7hwfb5\ntnexvRfwIvBIj9tfsj0/f30tMETSxqQRnebRoC2peSHTEEIIIaye0js8kjbJ/48BJgIzety+mfJH\nLEm7kmKeC9wObCvpTZJeBxwGXDWQsYcQQgihGkqd0souk7QRaTfnybbnSToOwPZU4FDgeEnLgIXA\nYU5j98sknQBcR9p57t9t31dOE0IIIYTQzUrv8Njes8WxqU1fnwucW/Cz1wLX9l90IYQQOkHSJGAS\nwJgxY0qOJqyNSp/SCiGEUH+2p9keb3v86NGjyw4nrIVKH+GpqmUFtanaKcp8KaqZpdZJZ4W1t9rV\nISrKiimq+/Xi/NZ1iF4/snUmS1GWztA2tYAGr9O3TKCiLKQiRb+nTdJNpfW1ve0ez6LXZNHPDBvS\n+nkuqh311IutswDb1YP7m01HtjxeVK9rg5FDCo63fg0X1a0qypSC4ozMopdY0d9ukcjGCqFzYoQn\nhBBCCLUXIzwhhBBC4I0nX9Py+B/POGiAI+kfMcITQgghhNqLDk8IIYQQai86PCGEEEKovVjDs5qK\n6je1U5ThsqCoFlBBZtXogkyPdtkkfc1w2mLD4S2Pv35E68yXxUv7nrVWlAm0YEnrx6OoPlFR24qy\nlla0eZy6tZbWmihqb7u2Fj2mRT9T9NorytLbpqCG29KCrEEozqIqypgcUZC1WPSnUFTnauTQ4stk\n0X23q78VukvRupVQPzHCE0IIIYTaiw5PCCGEEGovOjwhhBBCqL1YwxNCCCGEQq3WOVVxb57o8IQQ\nQqi9WJwcosMzgIoyN0YM7VvWVWEGTZvMkGVFWTQFNZBWuHVMSwvqgfW1DVCc6VZUl6koI6bo9xRl\n/BRl+9TV6mSeFf3M4oLaWEW11IruemFBJl675Kali5e1PL5lQUZhkaLnv6iuXLu6eUW1rvqaIVjU\n7sj2CqFz1q4rfwghhBDWStHhCSGEEELtld7hkfR/JN0r6T5Jn25x++GS7sn/bpW0U9Ntf5Q0R9Jd\nkmYPbOQhhBBCqIpS1/BI2gE4FtgVWALMknSN7UeaTvsD8C7b8yQdAEwDdmu6fR/bLwxY0CGEEEKo\nnLJHeN4C3GZ7ge1lwC+BQ5pPsH2r7Xn529uALQc4xhBCCCFUXNlZWvcCp0vaCFgIHAi0m5o6GvhF\n0/cGrpdk4Pu2p7X6IUmTgEkAW40Z04m4O6ooE2NQYYJG6xteXtQ6iwVgUMEvK84c68xLY3BxI9r8\nTGf64XXKxuqW13BRNlaRoudycN9+TVvzC173RbEuKsgQG1ZQe6uT6lirLYSqKLXDY/sBSVOAG4D5\nwN1Ay6uXpH1IHZ49mg7vbvtpSZsAN0h60PbNLe5nGmkqjHHjxvetimYIXSBewyGEsGbKHuHB9vnA\n+QCSvgE82fMcSTsCPwQOsD236Wefzv8/L+kK0lqg13R4QgghrD1ik8HQSulj/nl0BkljgInAjB63\njwEuB46w/XDT8ZGS1mt8DbyHNEUWQgghhPAqpY/wAJflNTxLgck5G+s4ANtTga8AGwHn5fUmy2yP\nBzYFrsjHBgMX2Z5VRgNCCCGEtUnRKFo319gqvcNje88Wx6Y2fX0McEyLcx4Ddup5PIQQwtohpq5C\nX5Te4VmbFNXA6lS9nKK6WND3zKei+kFFtbSKagoNbpOVEhkr9VH0eimqpVb4elmNDL11h7W+jBVl\nbxWd399/nyGEcpW+hieEEEIIob/FCE8IIYQQOqLVNGO3rOuJDk8IIYSuUcXFsKEaosMTQgih68UC\n5erqlk5sdHhCCCGE0DX6a1pMRZkJdSXpz8CfBuCuNgbWpirua2N7R9oePdB33A+v4ao9dxFvZ209\nEK/j5npwwFjgof6+zzXU7c9bkarGDasfe69ew2tdh2egSJqdN0hcK0R7q6tqbYl4w0Co6vNW1bih\n/2OPtPQQQggh1F50eEIIIYRQe9Hh6T/Tyg5ggEV7q6tqbYl4w0Co6vNW1bihn2OPNTwhhBBCqL0Y\n4QkhhBBC7UWHJ4QQQgi1Fx2eLqEoyVxr3f78Slqn6euujjWEgSBpZNkxhM6KNTwlkrQ5MMj2E2XH\nMhAkvQvYzPbPyo5lIEh6EzDK9t1lx9KOpPcAewBLbZ9Wdjyrkl9HGwBDbF9adjyrImlvYBNgsO2L\nSg4n9IKk/YC9gdNtLyw5nNUiaTTpb+TppmNyl7/pS9oHsO2bOv27Y4SnJJImAjcB50v6D0lvlzS0\n5LD6jaS/B74NrC2du4nA9cBZ+fk9VNL6ZcfVU76wnwM8AHxS0udKDqmtHO800o6sx0s6X9LGJYdV\nKF+8ZwBjgM9IOi9/0AldStIBwBTghgp3dg4FrgFmSjpN0p6QehHdPIKbP3ydDSzpl9/f5Z29WpK0\nEeki+AXbd0g6ExgFzARutN0vT3ZZ8ifcq4Ddbc+RtC6wwvaCciPrH5JGAD8Fptj+naTJwDbAI8B0\n2/9TaoCZpCHARcAs2+fnC/0OwMPAVd32STDH+xPgGtsX5tfRw8C1wEm2X+ymT7D5jWUK8Iztb0sa\nBpxP2jr/G7af66Z4A0gaC9wNHG17uqRNgBHAurbvLTe63snvLzOBycBzwAnAMOBW25eXGVs7+frz\nA2A/2w9KGk7qo3TsfSJGeMqxCBgKjAawfRLpzfB9wHZQu3UUy4F5wKZ5FOti4EJJV0jattzQ+sVy\n4PXAWwBs/xvwG1Kn511Q7vPbuG/bS4E7gO0kvRf4MfBm4MvAufmCU7oe8T6Sjw23PR+4BHgb8LV8\nTtd0HnIsdwJjJW1qexFwLLAp8NWmc0L3+CtwLrCbpL8jfSD4EnCjpONLjaz3BpHeXxbZfpY0sv44\nMEHSO0uNrL31gM2B5yUNIn24uVTSFElv78QdRIenBLZfJv0hjZP05nzs/wELgVPz97W5ENq+BTgC\n+C7pE8cvgKOBZ4CzSgytX9heTGrXrpLekY9dQprOOzJ/X+bz21xk75ekDtpxwMW2jyet59kpH+sG\nzfH+HvgQ8GVJPyJ1LN8DbC9ppzKC60nSVpKG5g7jb0gX8h1zJ20B6TWwm6SDSw00vEZe73IOMJ+0\n5GCm7WOAg4Cvd3mHAQDbzwOXAUdL2tz2XNL7DcCB5UXWXr5GHkV6X3gI+E/gZGBd4KOduI/Bnfgl\nYdUkHUT6o3kW+DlwNfAV4B8kXWf7v23/i6RrJI22/ecy411Tub3vJb14r7N9s6RPATvlEQ+AT0ma\nJWmz/EmksiTtACyz/WA+9CAwDjg4T1v8Lk9rHCRpW9uPlBTnAcBnJd0BvAxMtf1FSfsD4yWNtP2y\npJmkjlCpesQ7nzS//zSwPfAi8G+2F0r6A+kDQ6ny634KcCupo/MZ0vT1p9PNmmP7GUk30gWPb1ip\nMb1o+wlJ5wE3256Vj8+WNIPqPGf/BXwAOEzSxbaflnQOMEPSxra7qpp602N/gaSlwFjb38+3nQ5c\n1In3xejwDID8qeAs4Juk7JLrgEOBC4BPAJvlC/pwYGv6acHWQOnR3tcDV0s6xvZMSTc3nfeRfHvp\nb1RrQtKBpA7sOZIusH237f+WdDWp03dsHn14GdiC9EZdRpxvI42yHUkauv8wcImkjwOPAf8MzM2L\ngD9KumCWpkW8HwGuBD5u+7eS1rG9QtJRwM5AaWvC8rTblsAZpDUTDwAfB34LTAC+R/6UKukpUlt+\nWEqw4RV5Sn190sgh5A6N7SclPZu/dr5W7QmcWUqgvSRpkO3l+e9jY9IU+uclTQPeCoi0pKLbSBK5\n0zO9x217ACvoQNzR4RkYW5AWjF0AkD+Nfg/4JGm9xIGk6YOlwBHdsqh1DbRq7zclLbN9TV58ehhp\nuPKfqtzevEB5Amm0bhTwwfxp5S7bd0qaC2xL6kz8FfhIHmIuw2DgV3mKEUljSBeTH5LWj10MjCWt\nI/uA7QdKirOhVby7Az+RdER+U9qV9Nh+wvaTZQWapyifkPQb0kLq522fmT+t3gq8k7Se5x2k6cJ9\nbT9cVrwBJL0f+L/Ao8CTwEOSfpxHONexvUzS60h/G18EPmT78RJDfg1Ju5EWJC+wfbvt5ZKG2F6a\nr7XPAO8mvd8sBU7Ma99K1SLuFXndzvIe5x0NnAgcbvuva3y/NVoq0rUk7QJ8Cviy7WfysUOA7wPv\nsX2XUgbHoLy+p9IK2vs+0gr8A4A5wMeAX3fBm+oayZ/s32T7MUlbkBY4zgUut31n03nDSJlpAz56\nJ2kPYCtSptwc4Lt5eu0bpCnW0cDvGxkckgbbXjbQcfYh3g2B+21fImkU6e9mXonx/iNpQfq5pOy8\ne2x/o+n2L5A6vcfn9V2hZEqZTBcC/2r7/jxKeDwpu+k7tl9qOncCKdPuj6UEWyBP936HNH21CTDX\n9tH5tqHNr7U82rOwG95fVhH3K9cepSzMY4Hrbd/XkfuODk//yE/qUNtXKmUm/QR4Cvgs6QOhJZ0E\nLLZ9TpmxdkIv2/s5YIntc8p+U11Tub3DbF/R4/hWwCmk1OPvAPsB97mEzQeVdk8eQZpWGUL6lDqH\n1JH4PbARqQN6LLC17ZMHOsZmfYx3jO0vlBTqK5T2DTkT+Lzt6yS9EbiZtLZoSj7njaS2fLLkxeoh\nU9oT62rgq7b/fz72H6Rr1m22ZyhlaY2w/Z8lhtpSHg2ZTtqi4ae5438t8KztQ5vO2wO4vVs62n2I\newIw2ykzs2MiS6sf5IvgWcBf4JWsnUnA20kpglvnU4eSNiSrtD60d1jj64p3dhrtndfjuJx2zT6d\nNB1zESnjo6N/tL1le0Uevv4xaTTx/cAutseS3nz3z8/DCmBQvhiVpo/xDpZU6pR8fkP8KTApd3Y2\nJk2NvB/4tKTPSNqOtGPvLqT1aqEL5Gn06cCRko7IC2MXAfeTsv4gXZs7MrLQabaXs3LdEbZfsr0H\naeuPxmLfkaTXXtdszNmHuPfj1dmZHRFreDosXwTPIw1f3yRpPdKmVc/kxa3fBU7Nx7cnpdhWVrRX\n65E6rosbc8xOWR4LgB2BvW3fX17EACwjXbx/BExSSpVfLOkU0mjJ54D354tRN+htvGV3mueSOrNv\nyFMkl5Jiv4+0LmocaSprPHBkmdNuoaUZpHV1+wF/sd1YVH5w7kz/rNtG5CRt17T26yngZEm/bFpb\ndAgwVdJbSJmiZ5Yxjd7TasQ9pT/ijg5P520IvAQ8nTMAzgaQ9DxpaP5Y0kVwG+AB238oK9AOifYm\nT0n6je0f5U8oo4C/79Tc8xqaCXzQ9o2SdgZOA/49LxQcChxk+6FyQ3yVSsRr+yGlNPQrgNeRFsCe\nDxxDWph8cu78bhCdne7TGOWRNMP2CgBJHyNl0g7rhsW9zZQ2B71E0lW2D3PabXws8GtJu9t+3PYL\nkpYB6+fOWjd0drom7ljD0w8kfZiUmro+aUpjFrAXsD/wOVd8j52eor2vtPfdwBedSgYM6fT88+pS\nqt10Oilb6CTSNMyupI0GLywztlYqGO9bgX28cn8pJF1HKh1zZ57qjAttl8sLlz9LysaaU3Y8zfKH\nqMuAy4G/I62X/HC+7TTgYNLIc2NLiQO74cNlt8UdHZ4OkLQj6bG8u+nYRGBz2+fm79cjbYP/qW54\nIa6JaG/b9k62/Vg5kRaT9DXSnk+Tbf9cqajlo3nNUdepWrzNJH2AtEj5QNvPlR1P6B1JW5Oqiz9a\ndiyt5A8CL5HWQk4FljZ1Hg4BNiNNo57tLqr71U1xR4dnDWllDaJfAN+z/eum2wY11kUo7flwEmnt\nwfOlBNsB0d5qtjdnj21i+478/TqNYfxuVLV44ZUtCo4kjRJ8sEumM0MN5TVj00hZrx9W2qRzvu0/\nlRxaW2XHHR2eNaC0KdXZpGydx0jD7j9pflPM551Iqh11eDf1vPsq2lv99lZteqVK8eYOz7tIKbYP\nrur8ENZEzgr8FmmqaBApQaK0zTd7q8y4Y9HyGrC9RNKXgMWkYbkNgI/lT6O3NJ36HGmH3Up/4ov2\nVr+9Vek8NFQp3hzrTWXHEdYOeaHvPaS9qfavQmcHyo07RnhWQ84cWQzgpp2Cc9bO+0hZSVOANwOP\nu+JbyEd7Xzley/aGEKpH0gakdYP/avuesuPprTLjjo0H+0hph92fA5OBSyUd2bjNqQL2laSdYi8h\npddWukcZ7a13e0MI1ZS3OvjHKnV2oNy4Y0qrl/L8/EhSocLJtq9Sqgp+oVLdkqkAth/N6Y0bA7vm\nN8nKifbWu70hhOqz3Y2Vz1eprLijw9NLeX5+vqTZwKi8z8ptkg4jjQQssn2B0vb82wMTq7Cmo0i0\nt97tDSGEtU1MafXds8C+wHAA27OBI4ATJG1je7ntiW6qlF1x0d56t7ffSVou6S5J90q6VNKIsmMC\nkPTFDvyOb0l6UNI9kq6QFPWyQuhS0eHppTzlge3zSBWdp0paP48E/Aq4h1RHpxaivfVu7wBbaHtn\n2zuQtow/rrc/qP4taNrnDk+LeG4AdrC9I/AwUHoF9xBCa9HhaUPSWEkTJA2h6bGy/aH8/dnAUZIm\nk/bfqPQbYrQ3qWt7u8QtpLpqSLpS0h2S7pM0qXGCpPmSvibpt8AESV+RdHseIZrW6JxKuknStyXd\nLOkBSe+QdLmkRyR9ven3fVTS7/Io0/clDZJ0BjA8H5tedF6reJobY/t6ryxiehuwZf89dCGENREd\nngJKpQNmAl8nFQScLGlU43bbh5Eu3qOBvYGDq7IPQivR3nq3txsoVaA+AGjUKTrK9jhSNfETlXZh\nhbR4/F7bu+XRtXNtvyOPEA0H3tv0a5fY3ou0Zf1MUnbdDsAnJG2kVH35Q8DutncGlpM2iDyZlSNP\nhxedVxBPkaNIO3KHmqn5tOxpeUr2LknXK5WCqKXYh6eF/In/QuA7tn+tVBvnnaS9Wb7lVGW3+fyh\ntheXEGpHRHvr3d6ySVrOyk7OLaT9N5ZIOhU4JB9/I6m6/G1KVZOHemXZjg+QynaMIFWr/67tMyTd\nBJySn8N3k4p17p9/5mbgRGAP0tRVo9zHcGCG7VMlzbe9bj7/hDbnvSqegjaeQuq4TazSZomhd3q8\nVqYDd9g+q5c/O6jda6dTcfXhZ14Vj6RRtl/KX58IvNV2r6edqyRGeIqNIm0wB3AFcDXwOqBR9GxX\nSbvk2/ullP0Ai/bWu71laoyk7Gz7n3NnZ29gP2CC7Z2A35OKCwIsaursDCNVUz7U9t8CP2g6D/IG\nkcCKpq8b3w8GBPy46f7H2j61RYztzlu0is7Ox0mjTodHZ2etULdp2Zeavh1JjfcWiw5PC7aXAmcB\nEyXt6VS08FfAXcBekoYDuwNP5/Mr/QKJ9ta7vV1qfWCe7QWStieNsLXS6Ny8IGld4NA+3s+NwKGS\nNgGQtKFSVWyApXm0b1XnFZL0D8DnSVOeC/oYW6iYuk7LSjpd0hP5/K904KHqStHhKXYLcD1whKS9\ncjryRcDmwOa2v2372XJD7Khob73b221mAYOVauqcRlrw+xq2/0Ia1ZlD2uX69r7cie37gS8B1+f7\nugF4Q755GnCPpOmrOK+dc4H1gBvyp+qpfYkvVMZwSXcBs4HHSev+IHVy7ia9frdi5ajxcuCypp/f\nR9JvJc0B3g28rem2q/L/c4D7bD+Tp9Afy79zX2AccHuOYV9SWZue2p3XM55XsX2K7a2A6cAJbR+J\nCouNBwvYXpSHCQ18IX8KXUxaxDq/1OD6QbS33u0tU6s1BvmCfkBvzrf9JVJnpOd5ezd9fRNNhTt7\n3PYz4Gctfv7zpNGZVZ1XuEbC9jZFt4VaWZhHTV7RY1p2QV5T1m5adrztJ/LatdWZll3Vlgftzms7\nLdvkIuAa4Ku9OLdyYoSnDaeaHz8AziT1yvcBPmr7uVID6yfR3nq3N4TQUXWZlt226duDgQf7GF9l\nxAjPKtheAvxXzvpwXu9RW9Heerc3hNAxs4Dj8hToQ7SZlpXUmJb9I6sxLSupMd26DrCUtM7nT6yc\nlr0zr+MpOq+dMySNJY0o/Yk+bAxaNZGWHkIIIYTaiymtEEIIIdRedHhCCCGEUHvR4QkhhBBC7UWH\nJ4QQQgi1Fx2eLqd6F637YN6SfYWk8Z2IK4QQQmglOjzdr7F1+A6kmk69Thls1FHpJ33u8LSI515g\nInBzRyIKIYQQCkSHp1rqVrTuAdsP9feDFkIIIUSHpyLqWrQuhBBCGAix03L3axStgzTC01y07pD8\ndaNo3VxaF607CRgBbAjcB/w83/aaonUAkhpF6/ZgZTE6SB2m51vEuG+b89oWrQshhBAGQnR4ut/a\nUrQuhBBC6DcxpVVNtShaF0IIIQyU6PBU0yxgcC5adxptitaRqoHPAa5kNYrWAY1idPcANwBvyDc3\nitZNX8V5hSQdIulJ0mLmayRd15f4QgghhN6K4qEhhBBCqL0Y4QkhhBBC7UWHJ4QQQgi1Fx2eEEII\nIdRedHhCCCGEUHvR4QkhhBBC7UWHJ4QQQgi1Fx2eEEIIIdRedHhCCCGEUHv/Cx3NXMu9uLgrAAAA\nAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import sys\n", + "import os\n", + "from __future__ import print_function\n", + "import pints\n", + "import pints.toy as toy\n", + "import pints._diagnostics as diagnostics\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# Load a forward model\n", + "model = toy.LogisticModel()\n", + "\n", + "# Create some toy data\n", + "real_parameters = [0.015, 500]\n", + "times = np.linspace(0, 1000, 1000)\n", + "org_values = model.simulate(real_parameters, times)\n", + "\n", + "# Add noise\n", + "noise = 10\n", + "values = org_values + np.random.normal(0, noise, org_values.shape)\n", + "real_parameters = np.array(real_parameters + [noise])\n", + "\n", + "# Get properties of the noise sample\n", + "noise_sample_mean = np.mean(values - org_values)\n", + "noise_sample_std = np.std(values - org_values)\n", + "\n", + "# Create an object with links to the model and time series\n", + "problem = pints.SingleSeriesProblem(model, times, values)\n", + "\n", + "# Create a log-likelihood function (adds an extra parameter!)\n", + "log_likelihood = pints.UnknownNoiseLogLikelihood(problem)\n", + "\n", + "# Create a uniform prior over both the parameters and the new noise variable\n", + "prior = pints.UniformPrior(\n", + " [0.01, 400, noise*0.1],\n", + " [0.02, 600, noise*100]\n", + " )\n", + "\n", + "# Create a Bayesian log-posterior (prior * likelihood)\n", + "log_posterior = pints.LogPosterior(prior, log_likelihood)\n", + "\n", + "# Create differential evolution MCMC routine\n", + "x0 = real_parameters * 1.1\n", + "mcmc = pints.DifferentialEvolutionMCMC(log_posterior, x0)\n", + "\n", + "# Use 4000 iterations in total\n", + "mcmc.set_iterations(4000)\n", + "\n", + "# Discard the first 2000 iterations as burn in\n", + "mcmc.set_burn_in(2000)\n", + "\n", + "# Store only every 4th sample\n", + "mcmc.set_thinning_rate(4)\n", + "\n", + "# Disable verbose mode\n", + "mcmc.set_verbose(False)\n", + "\n", + "# Set gamma\n", + "mcmc.set_gamma(0.1)\n", + "\n", + "# Run!\n", + "print('Running...')\n", + "chains = mcmc.run()\n", + "print('Done!')\n", + "\n", + "# Plot output\n", + "import pints.plot # Use Pints' diagnostic tool\n", + "\n", + "stacked = np.vstack(chains)\n", + "pints.plot.pairwise(stacked, kde=False)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[1.0235890729350972, 1.0247660109632153, 1.2770439560500759]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "diagnostics.rhat_all_params(chains)" + ] + } + ], + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/inference-dream-mcmc.ipynb b/examples/inference-dream-mcmc.ipynb new file mode 100644 index 000000000..012c16dee --- /dev/null +++ b/examples/inference-dream-mcmc.ipynb @@ -0,0 +1,155 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Inference: DREAM MCMC\n", + "\n", + "This example shows you how to perform Bayesian inference on a time series, using [DREAM MCMC](http://pints.readthedocs.io/en/latest/mcmc/dream_mcmc.html).\n", + "\n", + "It follows on from the [basic MCMC example](./inference-first-example.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running...\n", + "Done!\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAI0CAYAAAAZXPP3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xu8XXV54P/Pk3sIxCQk4ZoQtNR6\nqRc4RS0zrZeRorSCVq12WtO+oPnZ0Rn7c6ZjnGmro2ONv47j1KnVSSsD2Kp1qggtKGaoSqf1QkAq\nKFgYDBBDSQIICUnIhef3x17RY/h+d84+2ZdzVj7v12u/9t7PXpfv2medkydrrWc9kZlIkiS12YxR\nD0CSJGnQTHgkSVLrmfBIkqTWM+GRJEmtZ8IjSZJaz4RHkiS1ngmPJElqPRMeSZLUeiY8kiSp9WaN\negDDtnTp0ly1atWoh6EWuPHGG7dn5rJRj0OSdHhHXcKzatUqNm7cOOphqAUi4u5Rj0GSNDGe0pIk\nSa1nwiNJklrPhEeSJLWeCY8kSWo9Ex5JktR6R12VliZv1dqri/FN684f8kgkSeqNR3gkSVLrmfBI\nkqTWM+GRJEmtZ8IjSZJaz4RHkiS1ngmPJElqPRMeSZLUeiY8kiSp9Ux4JElS65nwSJKk1jPhkSRJ\nrWfCI0mSWs+ER5IktZ7d0lVU64wuSdJ05BEeSZLUeiY8kiSp9QaW8ETEioj4YkTcFhHfioi3NPEl\nEbEhIu5onhc38YiID0bEnRHxzYg4c9yyVjfT3xERq8fFz4qIW5p5PhgRMajtkSRJ09cgj/DsB/5t\nZj4NeD7wpoh4OrAWuC4zzwCua94DvAw4o3msAT4MnQQJeAfwPOBs4B0Hk6RmmjXj5jtvgNsjSZKm\nqYElPJl5X2be1LzeAdwGnAJcAFzWTHYZcGHz+gLg8uz4KrAoIk4Cfg7YkJkPZuZDwAbgvOazhZn5\nlcxM4PJxy5IkSfqBoVRpRcQq4LnA14ATMvM+6CRFEbG8mewU4N5xs21uYt3imwtx9cBqLEnS0WDg\nFy1HxLHAp4HfysxHuk1aiOUk4qUxrImIjRGxcdu2bYcbsiRJapmBJjwRMZtOsvPnmfmZJnx/czqK\n5nlrE98MrBg3+6nAlsPETy3EnyAz12fmWGaOLVu27Mg2SpIkTTuDrNIK4KPAbZn5X8d9dBVwsNJq\nNXDluPgbmmqt5wMPN6e+rgXOjYjFzcXK5wLXNp/tiIjnN+t6w7hlSZIk/cAgr+E5B/hV4JaIuLmJ\n/QdgHfCpiLgIuAd4TfPZNcDLgTuBXcCvA2TmgxHxbuCGZrp3ZeaDzevfBC4F5gOfax4astJ1QJvW\nnT+CkUiSVDawhCcz/w/l62wAXlKYPoE3VZZ1CXBJIb4ReOYRDFOSJB0FvNOyJElqPRMeSZLUenZL\nP0p4vx1J0tHMIzySJKn1THgkSVLrmfBIkqTWM+GRJEmtZ8IjSZJaz4RHkiS1ngmPJElqPRMeSZLU\neiY8kiSp9Ux4JElS65nwSJKk1jPhkSRJrTehhCcinjnogUiSJA3KRI/wfCQivh4R/yoiFg10RJIk\nSX02oYQnM/8Z8C+BFcDGiPh4RLx0oCOTJEnqkwlfw5OZdwC/A7wN+FnggxFxe0S8alCDkyRJ6oeJ\nXsPzrIj4AHAb8GLgFzLzac3rDwxwfJIkSUdsokd4/gi4CXh2Zr4pM28CyMwtdI76PEFEXBIRWyPi\n1nGxd0bE9yLi5ubx8nGfvT0i7oyI70TEz42Ln9fE7oyItePip0fE1yLijoj4i4iY09umS5Kko8VE\nE56XAx/PzN0AETEjIo4ByMyPVea5FDivEP9AZj6neVzTLO/pwOuAZzTz/HFEzIyImcCHgJcBTwde\n30wL8L5mWWcADwEXTXBbJEnSUWaiCc//BuaPe39ME6vKzOuBBye4/AuAT2bmY5n5XeBO4OzmcWdm\n3pWZe4FPAhdERNA5nfaXzfyXARdOcF2SJOkoM9GEZ15m7jz4pnl9zCTX+eaI+GZzymtxEzsFuHfc\nNJubWC1+PPD9zNx/SFySJOkJJprwPBoRZx58ExFnAbsnsb4PA08BngPcB7z/4CIL0+Yk4kURsSYi\nNkbExm3btvU2YkmSNO3NmuB0vwX8r4jY0rw/CfilXleWmfcffB0RfwL8dfN2M517/Bx0KnBwXaX4\ndmBRRMxqjvKMn7603vXAeoCxsbFqYiRJktppQglPZt4QET8BPJXO0ZXbM3NfryuLiJMy877m7SuB\ngxVcVwEfj4j/CpwMnAF8vVnXGRFxOvA9Ohc2/3JmZkR8EXg1net6VgNX9joeSZJ0dJjoER6AnwJW\nNfM8NyLIzMtrE0fEJ4AXAksjYjPwDuCFEfEcOqefNgH/D0BmfisiPgV8G9gPvCkzDzTLeTNwLTAT\nuCQzv9Ws4m3AJyPiPwPfAD7aw7ZowFatvboY37Tu/CGPRJKkCSY8EfExOtfe3AwcaMIJVBOezHx9\nIVxNSjLzPcB7CvFrgGsK8bvoVHFJkiR1NdEjPGPA0zPT618kSdK0M9EqrVuBEwc5EEmSpEGZ6BGe\npcC3I+LrwGMHg5n5ioGMSpIkqY8mmvC8c5CDkCRJGqSJlqV/OSJOA87IzP/d9NGaOdihSZIk9ceE\nruGJiN+g07fqfzShU4DPDmpQkiRJ/TTRi5bfBJwDPAKQmXcAywc1KEmSpH6aaMLzWNOtHICImEWX\n3lWSJElTyUQvWv5yRPwHYH5EvBT4V8BfDW5YOhK1uxxLknS0mugRnrXANuAWOu0grgF+Z1CDkiRJ\n6qeJVmk9DvxJ85AkSZpWJtpL67sUrtnJzCf3fUSSJEl91ksvrYPmAa8BlvR/OJIkSf03oWt4MvOB\ncY/vZeZ/A1484LFJkiT1xURPaZ057u0MOkd8jhvIiCRJkvpsoqe03j/u9X5gE/Davo9GkiRpACZa\npfWiQQ9EkiRpUCZ6Suut3T7PzP/an+FIkiT1Xy9VWj8FXNW8/wXgeuDeQQxKkiSpnyaa8CwFzszM\nHQAR8U7gf2XmxYMamCRJUr9MtLXESmDvuPd7gVV9H40kSdIATDTh+Rjw9Yh4Z0S8A/gacHm3GSLi\nkojYGhG3jostiYgNEXFH87y4iUdEfDAi7oyIb44vg4+I1c30d0TE6nHxsyLilmaeD0ZE9LLhkiTp\n6DHRGw++B/h14CHg+8CvZ+bvH2a2S4HzDomtBa7LzDOA65r3AC8Dzmgea4APQydBAt4BPA84G3jH\nwSSpmWbNuPkOXZckSRIw8SM8AMcAj2TmHwKbI+L0bhNn5vXAg4eELwAua15fBlw4Ln55dnwVWBQR\nJwE/B2zIzAcz8yFgA3Be89nCzPxKZiado00XIkmSVDChhKc5jfU24O1NaDbwZ5NY3wmZeR9A87y8\niZ/Cj1Z8bW5i3eKbC/Ha+NdExMaI2Lht27ZJDFuSJE1nEz3C80rgFcCjAJm5hf62lihdf5OTiBdl\n5vrMHMvMsWXLlk1yiJIkabqaaMKztzl1lAARsWCS67u/OR1F87y1iW8GVoyb7lRgy2HipxbikiRJ\nTzDRhOdTEfE/6Fxb8xvA/wb+ZBLruwo4WGm1GrhyXPwNTbXW84GHm1Ne1wLnRsTi5mLlc4Frm892\nRMTzm+qsN4xbliRJ0o+YaC+t/xIRLwUeAZ4K/F5mbug2T0R8AnghsDQiNtOptlpHJ3m6CLgHeE0z\n+TXAy4E7gV10KsLIzAcj4t3ADc1078rMgxdC/yadSrD5wOeahyRJ0hMcNuGJiJl0jqr8CzpVUhOS\nma+vfPSSwrQJvKmynEuASwrxjcAzJzoeSZJ09DrsKa3MPADsiognDWE8kiRJfTfRXlp7gFsiYgNN\npRZAZv6bgYxKkiSpjyaa8FzdPCRJkqadrglPRKzMzHsy87Ju00mSJE1lhzvC81ngTICI+HRm/uLg\nh6SJWrXWg26SJE3E4RKe8Xc0fvIgB6KjQy1J27Tu/CGPRJJ0NDlclVZWXkuSJE0bhzvC8+yIeITO\nkZ75zWua95mZCwc6OkmSpD7omvBk5sxhDUSSJGlQJtpLS5Ikadoy4ZEkSa1nwiNJklrPhEeSJLWe\nCY8kSWo9Ex5JktR6JjySJKn1THgkSVLrmfBIkqTWM+GRJEmtZ8IjSZJabyQJT0RsiohbIuLmiNjY\nxJZExIaIuKN5XtzEIyI+GBF3RsQ3I+LMcctZ3Ux/R0SsHsW2SJKkqe9w3dIH6UWZuX3c+7XAdZm5\nLiLWNu/fBrwMOKN5PA/4MPC8iFgCvAMYAxK4MSKuysyHhrkR6o9Va69+QmzTuvNHMBJJUhtNpVNa\nFwCXNa8vAy4cF788O74KLIqIk4CfAzZk5oNNkrMBOG/Yg5YkSVPfqBKeBL4QETdGxJomdkJm3gfQ\nPC9v4qcA946bd3MTq8WfICLWRMTGiNi4bdu2Pm6GJEmaDkZ1SuuczNwSEcuBDRFxe5dpoxDLLvEn\nBjPXA+sBxsbGitNIkqT2GskRnszc0jxvBa4Azgbub05V0TxvbSbfDKwYN/upwJYucUmSpB8x9CM8\nEbEAmJGZO5rX5wLvAq4CVgPrmucrm1muAt4cEZ+kc9Hyw5l5X0RcC/z+wWquZjlvH+KmDFXpol5J\nkjQxozildQJwRUQcXP/HM/PzEXED8KmIuAi4B3hNM/01wMuBO4FdwK8DZOaDEfFu4IZmundl5oPD\n2wxJkjRdDD3hycy7gGcX4g8ALynEE3hTZVmXAJf0e4ySJKldRnkfHqmr2mk8788jSerVVLoPjyRJ\n0kCY8EiSpNYz4ZEkSa1nwiNJklrPhEeSJLWeVVpTjDcYlCSp/zzCI0mSWs+ER5IktZ4JjyRJaj0T\nHkmS1HomPJIkqfVMeCRJUuuZ8EiSpNbzPjyadkr3KrKDuiSpG4/wSJKk1jPhkSRJrecprRGxhYQk\nScPjER5JktR6HuFRK9SOmHkxsyQJWpDwRMR5wB8CM4E/zcx1Ix7SE3j6SpKk0ZrWCU9EzAQ+BLwU\n2AzcEBFXZea3RzsyTRUe+ZEkwTRPeICzgTsz8y6AiPgkcAEwkoTHIznTh/fykaSjy3RPeE4B7h33\nfjPwvEMniog1wJrm7c6I+M4QxlazFNg+wvUfiVaPPd7X8zJPm+xgJEnDNd0TnijE8gmBzPXA+sEP\n5/AiYmNmjo16HJPh2CVJ09V0L0vfDKwY9/5UYMuIxiJJkqao6Z7w3ACcERGnR8Qc4HXAVSMekyRJ\nmmKm9SmtzNwfEW8GrqVTln5JZn5rxMM6nClxam2SHLskaVqKzCdc8iJJktQq0/2UliRJ0mGZ8EiS\npNYz4ZEkSa1nwiNJklrPhEeSJLWeCY8kSWo9Ex5JktR6JjySJKn1THgkSVLrmfBIkqTWM+GRJEmt\nZ8IjSZJaz4RHkiS1ngmPJElqPRMeSZLUeiY8kiSp9Ux4JElS65nwSJKk1jPhkSRJrWfCI0mSWs+E\nR5IktZ4JjyRJaj0THkmS1HqzRj2AYVu6dGmedtqqUQ9DLXD33ZvYvn17DHu9S5cuzVWrVg17tWqp\nG2+8cXtmLhv1OKRBO+oSntNOW8XffW3jqIehFjjneWMjWe+qVavYuNF9WP0REXePegzSMHhKS5Ik\ntZ4JjyRJaj0THkmS1HomPJIkqfWOuouWNRiZWYxHDL2ISZKkJzDhkTStrVp7dTG+ad35Qx6JpKnM\nU1qSJKn1plzCExGbIuKWiLg5IjY2sSURsSEi7mieFzfxfxkR32wefx8Rzx7t6CVJ0lQ05RKexosy\n8zmZefDObmuB6zLzDOC65j3Ad4GfzcxnAe8G1g9/qJIkaaqbqgnPoS4ALmteXwZcCJCZf5+ZDzXx\nrwKnjmBskiRpipuKCU8CX4iIGyNiTRM7ITPvA2ielxfmuwj4XGmBEbEmIjZGxMZt27cNZNDjZWbx\nMRUdeDyLj15FRPEhSdJUMBWrtM7JzC0RsRzYEBG3H26GiHgRnYTnn5U+z8z1NKe7zjprbGpmHlIX\nTfK/BmDlypUjHo0kTT9T7ghPZm5pnrcCVwBnA/dHxEkAzfPWg9NHxLOAPwUuyMwHhj9iafAyc31m\njmXm2LJlNraWpF5NqYQnIhZExHEHXwPnArcCVwGrm8lWA1c206wEPgP8amb+4/BHLEmSpoOpdkrr\nBOCK5tqPWcDHM/PzEXED8KmIuAi4B3hNM/3vAccDf9zMs39cZZeko1jphoTejFA6ek2phCcz7wKe\ncC+d5lTVSwrxi4GLhzA0SZI0jU2phKcteq1OerxSFTVjRnk5+/Y/XozPnlU/Q1mrvJpZWUev00uS\nNJVNqWt4JEmSBsGER5IktZ4JjyRJaj0THkmS1HomPJIkqfWs0pqk/QfKlVJQr2SqVW89XumztXdf\neR2zKsvfvfdAdUxzKhVcte2obUOtoqy2DbNmmlOrf0r31pGkifBfI0mS1HomPJIkqfVMeCRJUuuZ\n8EiSpNYz4ZEkSa1nldYkdas+ykrFUi2+70BvlU+1Sqm9lR5bAI/tK1dw1fpvzahUlNUqxGrtw2o9\nuQBqXblqPcQkSZosj/BIkqTW8wiPpKNG7T4+m9adP+SRSBo2j/BIkqTWM+GRJEmtZ8IjSZJaz2t4\nGrUKqmr/qy7VR72q9bOqraFWw7S/y5ge3rWvGD9uXnkXqFVK1Xps1ZbTrUqrVow1o7qFkiRNjkd4\nJElS65nwSJKk1jPhkSRJrec1PJKmnNr9ciRpsjzCI0mSWs8jPI1aNVZNrc8V1PtQ7a71s6r05Xpg\n595ifE9lObv3luNQ7/215aE9xfhJi+cV47Wiq+07ymNdcuyc6phqX2FtO+ZU+n7VKsckSTqor0d4\nImJhRLw3Ij4WEb98yGd/3M91SZIkTVS/T2n9Tzq3ifk08LqI+HREzG0+e/5EFhARmyLiloi4OSI2\nNrElEbEhIu5onhc38Z+IiK9ExGMR8e/6vC2SJKkl+p3wPCUz12bmZzPzFcBNwN9ExPE9LudFmfmc\nzBxr3q8FrsvMM4DrmvcADwL/Bvgv/Ri8JElqp34nPHMj4gfLzMz3AOuB64Fek57xLgAua15fBlzY\nLH9rZt4AlG8jLLVERKyJiI0RsXHbtm2jHo4kTTv9Tnj+Cnjx+EBmXgb8W6B8VesTJfCFiLgxItY0\nsRMy875mefcBy/s0XmlayMz1mTmWmWPLli0b9XAkadrpa5VWZv77SvzzwBkTXMw5mbklIpYDGyLi\n9iMdV5M4rQFYsXLlkS4OqFc9dTOnMs+BSrlSrfqoFt+8Y1d13VHpT3XPI7uL8RN3zi3Gz1hyXDE+\ne2Z5+d2K32r9t2pVa9ZiSZIma8rdhycztzTPW4ErgLOB+yPiJIDmeWuPy/zh/46X+r9jSZKONlMq\n4YmIBRFx3MHXwLnArcBVwOpmstXAlaMZoSRJmo76fuPB5qLl52fm309i9hOAK5qbAM4CPp6Zn4+I\nG4BPRcRFwD3Aa5p1nQhsBBYCj0fEbwFPz8xH+rApkiSpJfqe8GTm4xHxfuAFk5j3LuDZhfgDwEsK\n8X8CTp3MOCVJ0tFjUKe0vhARvxi99muQJEkagEH10norsAA4EBG76RTYZGYuHND6BubxSvOoGV36\nN9XmqXXf2rF7f3k5leqtv793ezH+7fvLFVcA37lvRzF+9umLi/Hb7v9+Mf7F/1uOP3X5/GL8Xzy5\nfpH4vNkzK588XoweeLz8nc+qVIjV+n7VenJJktprIAlPZpZrlyVJkkZgIP/VjY5fiYjfbd6viIiz\nB7EuSZKkwxnUsf0/pnPR8sGO6TuBDw1oXZIkSV0N6hqe52XmmRHxDYDMfCgi5gxoXZIkSV0NKuHZ\nFxEzaa7TjYhl1K5ElaQRW7X26mJ807rzhzwSSYMyqITng3TaQiyPiPcArwZ+d0DrGqhu1Vg1e/eX\nc7taldb+SjlRbTnL5s+rLKlepXX35oeL8Zu+sbkYnzu/3Evr5JPL16NvffiYYvyU48rLAXj2rEXF\n+JJjywcDH6t8H7Nmlqu9av29JElHn0FVaf15RNxI52aBAVyYmbcNYl2SJEmHM5CEJyI+lpm/Ctxe\niEmSJA3VoKq0njH+TXM9z1kDWpckSVJXfU14IuLtEbEDeFZEPBIRO5r3W7HDuSRJGpG+JjyZ+d7m\nLst/kJkLM/O45nF8Zr69n+uSJEmaqEFVaf3HiPgV4PTMfHdErABOysyvD2h9A1Pri9WtLWqtt9P+\nA5UeW5WeWbWeT7NnlOOVxQPw0IOPFuPzF5Qrvk488dhi/Marv1yMP/mnyzfS/nSXvlWrnrSgGN+9\n90Axfsricr+umgOVn93MLmm+/W4lqZ0GdQ3Ph/BOy5IkaYrwTsuSJKn1BnWExzstS5KkKWNQCc+h\nd1r+P8DvD2hdkiRJXXmnZUkjU+thJUn9NqhreADuB/62Wcf8iDgzM28a4PoGotZLq1ZZ1fmst3U8\naf7sYvyBnXuL8Uf37y/Gzzy5XPUEcNXc8jpqRUk3fuaa8ge7dxTD2+8v9+p67MnHV8e0ZWe599fT\nli8sxh+tVG/VepEdO6+8e1uJJUlHn0G1lng38GvA/+WHPTMTePEg1idJktTNoI7wvBZ4SmaWD1FI\nkiQN0aAuWr4VWDSgZUuSJPVkUEd43gt8IyJuBR47GMzMVwxofZIkSVWDSnguA94H3IL335EkSSM2\nqIRne2Z+cEDLnhK6VfrMmlmbpxzfva9cfTSzUiF22sJjivFvbvt+dUzPfcYJxfjXv/G9YvzFb3hl\nMf43l19RjJ94ypJi/NTjy2MFWDinXDlW66VVs2RBeTn79pdz7az0OgOY1a3R1ghFxBpgDcDKlStH\nPBpJmn4G9df9xoh4b0S8ICLOPPgY0Lqk1svM9Zk5lpljy5YtG/VwJGnaGdQRnuc2z88fF5tQWXpE\nbAJ2AAeA/Zk5FhFLgL8AVgGbgNc2/bkC+EPg5cAu4Nem471+JEnSYA3qTssvOsJFvCgzt497vxa4\nLjPXRcTa5v3bgJcBZzSP5wEfbp4lSZJ+YGB3Wo6I84FnAPMOxjLzXZNc3AXAC5vXlwFfopPwXABc\nnp3bHn81IhZFxEmZed9kxy1JktpnINfwRMRHgF8C/jWdXlqvAU6b4OwJfCEibmwu1AQ44WAS0zwv\nb+KnAPeOm3dzE5MkSfqBQR3h+enMfFZEfDMz/1NEvB/4zATnPSczt0TEcmBDRNzeZdpSuc0TGiuN\nr3BZMcIKl1pl19xZ5bxzRmX6WZUqo5/YX+5BBfDoj5crllYsLfffuv175d5Yp57zM8X4gQPl5W97\nZE91TLsq1Wl7d5TnWbZgbjH+WKUaa/7scrlcj63OJEktMKgqrYP/Yu2KiJOBfcDpE5kxM7c0z1uB\nK4Czgfsj4iSA5nlrM/lmYMW42U8FthSW+cMKl6VWuEiSdLQZVMLzVxGxCPgD4CY6lVWfONxMEbEg\nIo47+Bo4l06biquA1c1kq4Erm9dXAW+IjucDD3v9jiRJOlTfT2lFxAw6FVXfBz4dEX8NzMvM8jmS\nH3UCcEVz6mcW8PHM/HxE3AB8KiIuAu6hc00QwDV0StLvpFOW/uv93RpJktQGfU94MvPx5pqdFzTv\nH2NcP63DzHsX8OxC/AHgJYV4Am86ogFLkqTWG9QprS9ExC9Gt/4LkiRJQzKoKq23AguA/RGxh041\nVWZmvYzoKFDL/uZUqrTmzu6tqisX19ddW8fsmQ8V4w/sKB+U++lzzyjGf/KEcrXXzfftrI7pkX37\nivGnLjmuGJ9XqbraW6nSmlPpizWj0qOsm87BxENiPS9FkjQqg7rTcvlfLEmSpBEY5J2WF9Np+TD+\nTsvXD2p9kiRJNQNJeCLiYuAtdO6LczOdJqJfYQLNQyVJkvptUBctvwX4KeDuppHoc4FtA1qXJElS\nV4M6pbUnM/dEBBExNzNvj4inDmhdkjQQq9Ze/YTYpnXnj2Akko7UoBKezc2dlj9Lpx/WQxRaPrRV\nrRq/VqSf5SKjLssvx+fPKVcxASxaMKcYP/fHTijG/9mKpcX4nQ/uKMYfr5QsPevEY6tjetqyctHe\ngrnl7aito1aNVfueutVolaqxOst64lzec0GSpo9BVWm9snn5zoj4IvAk4PODWJckSdLh9DXhiYh5\nwBuBHwNuAT6amV/u5zokSZJ61e+Lli8DxugkOy8D3t/n5UuSJPWs36e0np6ZPwkQER8Fvt7n5UuS\nJPWs3wnPD3oFZOZ+W2lJOqhU8SRJw9LvhOfZEfFI8zqA+c17e2lRrwCaWent9HilLKk2/eyZ9QTz\nuHnlH/WByjoe2b2/GF+1qFx1tWfvgWJ8QWW9AMdUqrGOmVuep9Yza1Zlu2tVXd1aadW+j9o6JEnT\nQ18Tnsys10VLkiSNyKDutCxJkjRlmPBIkqTWM+GRJEmtZ8IjSZJab1C9tFTQrzL9GZXlzJtdv2Z8\nRvTWsGvFkvnF+K5KNdacWeXcuVtFVK2KqlYpdWyXiq+SWbXqt0q1HNQr4CRJ05tHeCRJUuuZ8EiS\npNbzlJY0DUTEGmANwMqVK0c8mqNb7Y7Rm9adP+SRSOqFR3ikaSAz12fmWGaOLVu2bNTDkaRpxyM8\nkvrKnlmSpiITnilsRh8rhmZXqqhqstK3av6cciVYrbppMn2rahVfNbWqter316VgzYa3ktROU+6U\nVkTMjIhvRMRfN+9fHBE3RcStEXFZRMxq4osj4oqI+GZEfD0injnakUuSpKlqyiU8wFuA2wAiYgZw\nGfC6zHwmcDewupnuPwA3Z+azgDcAfziCsUqSpGlgSiU8EXEqcD7wp03oeOCxzPzH5v0G4Beb108H\nrgPIzNuBVRFxwhCHK0mSpokplfAA/w349/zwKovtwOyIGGvevxpY0bz+B+BVABFxNnAacOrwhipJ\nkqaLKZPwRMTPA1sz88aDscxM4HXAByLi68AOYH/z8TpgcUTcDPxr4BvjPjt02WsiYmNEbNy2fdsg\nN0OSJE1BU6lK6xzgFRHxcmAesDAi/iwzfwX45wARcS7w4wCZ+Qjw6008gO82jyfIzPXAeoCzzhqr\nN1JqsVoVVVb6StUqpWrT16qbJtObqtd11Kav6Vb91uu6JUnTw5Q5wpOZb8/MUzNzFZ2jOn+Tmb8S\nEcsBImIu8DbgI837RRExp5mpz9+xAAAgAElEQVT9YuD6JgmSJEn6EVPpCE/Nbzenu2YAH87Mv2ni\nTwMuj4gDwLeBi0Y1QEmSNLVNyYQnM78EfKl5/dvAbxem+QpwxlAHJkmSpqUpc0pLkiRpUEx4JElS\n603JU1oanmp1VbUoafDVSrUxPV7pvdXPnmNWY2mySk1TN607fwQjkVTiER5JktR6JjySJKn1THgk\nSVLreQ2PpEkrXbciSVORR3gkSVLreYRHPan1mqq1s+pnBVU/lyVJOrp4hEeSJLWeR3gkaUBq1zh5\nfx5p+DzCI0mSWs+ER5IktZ6ntCQdluXnkqY7Ex71pNZryhZUkqSpzIRHkobMi5ml4fMaHkmS1Hom\nPJIkqfU8pSXpB7w4WVJbmfBI0hTRS8Lp9T5Sb6LWG6mtImIbcPcQV7kU2D7E9Q2a2/NDp2Xmsn4O\npiYi1gBrmrdPBb4zjPVOQW3b/yarn9/D0PZjaZSOuoRn2CJiY2aOjXoc/eL2aJT8eXX4PUi986Jl\nSZLUeiY8kiSp9Ux4Bm/9qAfQZ26PRsmfV4ffg9Qjr+GRJEmt5xEeSZLUeiY8kiSp9Ux4JElS65nw\nSJKk1jPhkSRJrWfCI0mSWs+ER5IktZ4JjyRJaj0THkmS1HomPJIkqfVMeCRJUuuZ8EiSpNYz4ZEk\nSa1nwiNJklrPhEeSJLWeCY8kSWo9Ex5JktR6JjySJKn1THgkSVLrmfBIkqTWM+GRJEmtZ8IjSZJa\nz4RHkiS13qxRD2DYli5dmqedtmrUw5iQ7HH6GMgojky3bZiK4+3F3XdvYvv27UPfjDbsw5P50npd\nVj/X3WY33XTj9sxcNsx1Ll26NFetWjXMVarFbrxxYvvwUBKeiLgE+Hlga2Y+s4m9Bngn8DTg7Mzc\nWJn3POAPgZnAn2bmuiZ+OvBJYAlwE/Crmbn3cGM57bRV/N3XiquacjJ7S3kipt6f8scfr2/DjBlT\nb7y9OOd5YyNZbxv24cnsq70u60Bl35s5zfe7fps/O+4e9jpXrVrFxo3TYx/W1BcxsX14WKe0LgXO\nOyR2K/Aq4PraTBExE/gQ8DLg6cDrI+LpzcfvAz6QmWcADwEX9XnMkiSpJYaS8GTm9cCDh8Ruy8zv\nHGbWs4E7M/Ou5ujNJ4ELovNfuhcDf9lMdxlwYZ+HLUmSWmKqX7R8CnDvuPebm9jxwPczc/8h8aKI\nWBMRGyNi47bt2wY2WGlQ3Icl6chM9YSndLI9u8SLMnN9Zo5l5tiypUO9Nk/qC/dhSToyU71KazOw\nYtz7U4EtwHZgUUTMao7yHIy3Su1izNrFm8O4QLjXC0en+4XJOjL9vJC+12XVdr3a78ko99V+Xtwt\nqWyqJzw3AGc0FVnfA14H/HJmZkR8EXg1net6VgNXjm6YkiRNzqq1Vxfjm9adP+SRtNtQTmlFxCeA\nrwBPjYjNEXFRRLwyIjYDLwCujohrm2lPjohrAJqjN28GrgVuAz6Vmd9qFvs24K0RcSeda3o+Ooxt\nkSRJ089QjvBk5usrH11RmHYL8PJx768BrilMdxedKi5JkqSupvpFy5IkSUdsql/DI0lSa9Su19Hg\nmfBMAb1WaNTij3dpRVFbR+32+7Nmlg/+1acffDWJlSzt163ScN+Bx4vx2ZV9tbZb1OL93L96rQRz\nH5YGz1NakiSp9Ux4JElS65nwSJKk1jPhkSRJrWfCI0mSWs8qrUmqVXRA7xUXtelrFVGTafmz/0Bv\nVVpQrojZvfdAMT5/TnkpXb6masXKTCtZRq7Xqr4Z1YrC8vLr+12XMVXie/eX99Va1eKcWeX/59W2\nodvvek1t364tq7YKe9FJ/eMRHkmS1HomPJIkqfVMeCRJUut5DY8kSX1mC4mpZyhHeCLikojYGhG3\njostiYgNEXFH87y4MN+LIuLmcY89EXFh89mlEfHdcZ89ZxjbIkmSpp9hHeG5FPgj4PJxsbXAdZm5\nLiLWNu/fNn6mzPwi8BzoJEjAncAXxk3y25n5lwMcd9VQ+utUltPPCpdee/vUKqhqvbf2VSpouqlt\nX+37m12putHk1X/+vU1f+5nVKqtq+xfUfx+iMs+Myhy1iqgDlQ+69Ynr9Vex3h+vt+VI6t1QEp7M\nvD4iVh0SvgB4YfP6MuBLHJLwHOLVwOcyc1efhydJ0pRTOy22ad35Qx5JO4zyv8YnZOZ9AM3z8sNM\n/zrgE4fE3hMR34yID0TE3NqMEbEmIjZGxMZt27cd2ailEXAflqQjMy3OBUTEScBPAteOC78d+Ang\np4AldDk6lJnrM3MsM8eWLV020LFKg+A+LElHZpQJz/1NInMwodnaZdrXAldk5r6Dgcy8LzseA/4n\ncPZARytJmrQfOUq5zaOUGr5RJjxXAaub16uBK7tM+3oOOZ01LlkK4ELg1sJ8kqQp4EeOUi7zKKWG\nbygXLUfEJ+hcoLw0IjYD7wDWAZ+KiIuAe4DXNNOOAW/MzIub96uAFcCXD1nsn0fEMjrFGzcDbxzE\n2Gu9byZTpdVzf50el9+t5U+tKmZupcLpsX3lnlm799Wqa8rT1/oTQb2Cp7YZs7tUy2i09h/orRqv\nVo21v0vZU23/3n+gvO/NrpSU1XpszZ8zsxjvVhVZ+zvQz78bkvpjWFVar6989JLCtBuBi8e93wSc\nUpjuxf0anyRJardpcdGyJEnSkTDhkSRJrWfCI0mSWs+ER5IktZ7d0g+jn1UVtcqNWhHIvkrlS60/\n1Y49+6vr3n+gvJK5s8s5b236WhXN7r3lSpmF8+u7WG07Fswtz7O38v3NnV2urtHk1SqTurS6Knqs\nUtW3u1IF+ODOvdVl1fq1VXtsVT6oVW/VKhnnVaq3AGbV+nhVVj6j8l/M2t+GWjWb1V5S7zzCI0mS\nWs+ER5IktZ4JjyRJaj2v4ZEk6QisWnv1qIegCfAIjyRJaj2P8AxRrRdQrRqrplahUasYgXr1y7z9\n5QqUWvegBx59rBhffty8YvzhXfuKcYBjKtUvtXXPrGx3rxVFVrh01CqDoP7d1SqZalV6tf5x+ypV\ngLXld8ZUXtbOSnVirQLxkQPl6U980txivFuPulrV4pxK37daZWJtrLUKTtvKSb2b0BGeiDgxIj4c\nER+KiOMj4p0RcUtEfOpg13JJkqSpaqKntC4Fvg3cC3wR2A2cD/wt8JGBjEySJKlPJprwnJCZ/z0z\n1wGLMvN9mXlPZv534LTDzRwRl0TE1oi4dVxsSURsiIg7mufFlXkPRMTNzeOqcfHTI+Jrzfx/ERFz\nJrgtkiTpKDPRhGf8dJdPYhmXAucdElsLXJeZZwDXNe9Ldmfmc5rHK8bF3wd8oJn/IeCiCYxDkiQd\nhSZ60fKVEXFsZu7MzN85GIyIHwP+8XAzZ+b1EbHqkPAFwAub15cBXwLeNpHBROeq0xcDvzxu/ncC\nH57I/JIkTVelMvhN684fwUimlwklPJn5e5X4ncCrJ7nuEzLzvmY590XE8sp08yJiI7AfWJeZnwWO\nB76fmQfLLTYDp0xyHH33eKW04vFKuUetMqW2nJ2PlStiatUqANt2laur9lYqxPYcKK/j0X3ldRyo\nbNvyY8vVW1DvaVSr+Kn1FTp2Xnk3rlXXdCvSqn3ntWqj6axbtdr+yn5RK1jaW6m6ysq+/f1Hy1WD\ndzy0ozqmE4+ZX4zf9mB5noWVnmwrjzumGN/0wK5i/ISF9X34uMq+t6fSQ2zB3HJlYq0aq7bb1fZT\naOe+KvXDdChLX5mZWyLiycDfRMQtwCOF6ap/ASJiDbAGYMXKlYMZpTRA7sOSdGRGeePB+w+WtDfP\nW0sTZeaW5vkuOqe9ngtsBxZFxMGE7VRgS21Fmbk+M8cyc2zZ0mX92wJpSNyHJenITDjhiYgZEfHT\nfVz3VcDq5vVq4MrCOhdHxNzm9VLgHODb2blj2hf54em04vySJEnQQ8KTmY8D75/MSiLiE8BXgKdG\nxOaIuAhYB7w0Iu4AXtq8JyLGIuJPm1mfBmyMiH+gk+Csy8xvN5+9DXhrRNxJ55qej05mbJIkqf16\nvYbnCxHxi8Bnstt96Q+Rma+vfPSSwrQbgYub138P/GRlmXcBZ090DJIk6ejVa8LzVmABcCAidgMB\nZGYu7PvIpohaXtct3at91GsPrD2VypdHK9VYteomgPse3V2M3/NQuVrmru3lipVHdpd7Y/18MS2F\nRfPq94N89LHydsybXZ6nVuVW66VV022nn84VLrV9tbbfdav0qfW6qqlVddUqB3dV9tWHulQa3v39\n7xfjX/7O9mL82SsXFeN3bi//Ljzl+HI11tJjyz22oF5hedz88l5W6wdW/TvD9N0fpammp4QnM48b\n1EAkSZIGpacqrej4lYj43eb9iojwtJIkSZrSei1L/2PgBfzwDsc7gQ/1dUSSJEl91us1PM/LzDMj\n4hsAmfmQTTslSdJU12vCsy8iZtJclxsRy4DyVXuSJLVIqYeVpo9eE54PAlcAyyPiPXRu/Pe7fR/V\nCPRQZd+ZvstntbqKWt+qWqXH/koVTa3C5cE95X5ZAH/33VI3Dtj+yJ5i/PTlx5bXsbO8jk98vXyj\n6396WrmqC+CVTzuxPKYd5cqxeXPKfYjmV+JzKr26atVeANHDbtDbHjN4vfYO6/Y9zJ5ZnunRSh+3\n2j5cW8VdD+8sxv/m9geqY6pV0G2tVBR+tVIFuG9feRsWPe/UYnxPl+rHWlHf3Nm9XS0wd1Zv00/n\nakJpVHqt0vrziLiRzv1zArgwM28byMgkSZL6pKeEJyI+lpm/CtxeiEmSJE1JvVZpPWP8m+Z6nrP6\nNxxJkqT+m1DCExFvj4gdwLMi4pGI2NG834pNOyVJ0hQ3oYQnM9/b3GX5DzJzYWYe1zyOz8y3D3iM\nkiRJR6TXKq3/GBG/Apyeme+OiBXASZn59QGMbahq/YZqZnSp0an1dqr1Lqr1LXpgR7kiqtYX66v3\nliuxALY8WK5k2fJP5WqZjTdtLsaXn/ikYvzkE8pVXfu79GTavrNcjbX4mNnF+MxKZUqtEqhWOTSr\nUoEEUGsvVVr3VKuTqVXu1CoQZ1Wq2KB7X7aSBys/y9r3ecs/lffH736vvg8vWFDeL75zy93F+Kof\nP6UY31/ZL/7+rnKvrlMWlntsAfzY4vJ+f8zc8p/W2eWCwnpvvknsZL32VJOOFr1ew/MhJnGn5Yi4\nJCK2RsSt42JLImJDRNzRPC8uzPeciPhKRHwrIr4ZEb807rNLI+K7EXFz83hOj9siSRqSiFgTERsj\nYuO2bdtGPRwdhXpNeJ6XmW8C9kDnTsvARO60fClw3iGxtcB1mXkGcF3z/lC7gDdk5jOa+f9bRIxv\ngfzbmfmc5nFzb5siSRqWzFyfmWOZObZs2bJRD0dHoaHcaTkzr4+IVYeELwBe2Ly+DPgS8LZD5vvH\nca+3RMRWYBlQPvYsSdJRqHYX6E3rzh/ySKauXo/wHHqn5f8D/P4k131CZt4H0Dwv7zZx05V9DvB/\nx4Xf05zq+kBEzJ3kOCRJUstNizstR8RJwMeA1Zl58IjS24F/opMEradzdOhdlfnXAGsAVqxcOejh\nSn3nPixJR6bXU1oA9wN/28w7PyLOzMybJrOciDgpM+9rEpqtpYkiYiFwNfA7mfnVg/GDR4eAxyLi\nfwL/rraizFxPJynirLPGemqB1GuPLahXJtV6Y9XitWqvY2aVf2wL5tQP2N3wte8W44uXlauudu0s\nV9FsurMc37lzSTF+0pJjqmPa8mh5WbNnLCjHK1VFuyt9iI6dV/6eulWrTNUWRUeyD9e2t1Y1CPBY\npd/Unn3lM9hLji1fyveP23cU4/dWqga7uemrdxbjB+7YWIzfN7980Hff3nJ/t+OPn1+MP/xYvR9c\n7Xe0thvtr/TTq1YOVi4Y6NZLy2osqazX1hLvBn6Nzmmlg7/pCbx4Euu+ClgNrGuen3ADw4iYQ+cU\n2uWZ+b8O+exgshTAhcCth84vSZIEvR/heS3wlMws33SjIiI+QecC5aURsRl4B51E51MRcRFwD/Ca\nZtox4I2ZeXGzvp8Bjo+IX2sW92tNRdafNxdNB3Az8MYet0WSJB0lek14bgUWUTn9VJOZr6989JLC\ntBuBi5vXfwb8WWWZkzmqJEmSjkK9JjzvBb7R3EDwB7cBzsxX9HVUkiRJfdRrwnMZ8D7gFiZw/x1J\nkqSpoNeEZ3tmfnAgI5miahUP3aq3ahUXM2eUG+nseqxcEbOg0o/nG/eX77t41/Zyjy2Apz/r1GL8\n4YfL/boWLllYjG+9/tpifO5pFxTjDzyypzqmVQvLfYhq1VX7KhUutR5bk6m4amOFy2QqDWu9oGB/\nMbpjTzm+a3953659zzt3lvdHgAMPP1D+4JhKpeH27cX4SU99SjG+fHG5ovCUY8vVW93UdqNu/ctK\nZlcqECX1rteE58aIeC+dCqvxp7QmU5YuSZI0FL0mPM9tnp8/LjbZsnRJkqSh6PVOyy8a1EAkSZIG\npec7LUfE+cAzgHkHY5lZbOkgSZI0FfR0RVxEfAT4JeBf07nh32uA0wYwLkmSpL7p9QjPT2fmsyLi\nm5n5nyLi/cBnBjGwqa5b4Uutv86MSulGrS9OLX7awnI1yfynlqvAAD72ULl30d0PPVSeoVJlMv8n\nf7oYX7as3P/quPmzq2O6f1e5quy4eccV47XKoVqVVv1nVP/h1eaZztVbk6k0rFW4VdrEVfty/dii\nciXejy8v749/t6fetyqOW1yML37Kk4vxfZUeWKedtqgYP3ZeeV+tVZoB7Kr0HNu7v1xRWKu8rP1t\nqOn2s5vO+6o0SL3WPB6sMd4VEScD+4DT+zskSZKk/ur1CM9fRcQi4A+Am+j8V/lP+j4qSZKkPppw\nwhMRM4DrMvP7wKcj4q+BeZn58MBGJ0mS1AcTPqWVmY8D7x/3/jGTHUmSNB30ekrrCxHxi8Bnsof7\n1UfEJcDPA1sz85lNbAnwF8AqYBPw2sx8whW0EbEa+J3m7X/OzMua+FnApcB84BrgLb2MSZKkklVr\nrx71EPqmtC2b1p0/gpGMXq8Jz1uBBcD+iNhDp5YnM7PceOmHLgX+CLh8XGwtnVNk6yJibfP+beNn\napKidwBjdK4XujEirmoSow8Da4Cv0kl4zgM+1+P2TFq3Qohav5z9lV5Qx8wpV1ft2VuuADn+mLnF\n+PY99T5EK5aVq2X2V6prHnqoXEF1yinlvkUrK8s/eVF5rAAnLSj3KHq8krfOq3xPtaq4yRSrHE0V\nLrXvGepVRrNqFXGV5Ty0Z28xfuy88s/y/J8p97kC+Mf7lhbjmzaVe8s9aVG5cvCkJeUqx2ecXN6H\nTzpmXjEO8KRj6lWIJbVqtur3WvkZHU37qdQvPVVpZeZxmTkjM+dk5sLm/eGSHTLzeuDBQ8IX0Om+\nTvN8YWHWnwM2ZOaDTZKzATgvIk4CFmbmV5qjOpdX5pckSZrUnZYXA2fwo3davn4S6z4hM+9r5r8v\nIpYXpjkFuHfc+81N7JTm9aFxSZKkJ+gp4YmIi4G3AKcCN9NpIvoVBtc8tHTcNrvEywuJWEPn9Bcr\nVq7sz8ikIXIflqQj0+uNB98C/BRwd9NI9LnAtkmu+/7m1BTN89bCNJuBFePenwpsaeKnFuJFmbk+\nM8cyc2zZ0mWTHK40Ou7DknRker7TcmbuAYiIuZl5O/DUSa77KmB183o1cGVhmmuBcyNicXMq7Vzg\n2uZU2I6IeH50rt57Q2V+SZKknq/h2dzcafmzwIaIeIguR1YOiohPAC8ElkbEZjqVV+uAT0XERcA9\ndBqREhFjwBsz8+LMfDAi3g3c0CzqXZl58OLn3+SHZemfY4gVWs046x/2WFkxI8rTL1pQrgDZXane\n2lepAgNYWllWLi9Xppx1+pJifPvOctVNrQLtx5eVK7EA5leqrpYcO6cYr/V3OnZeeTeuVcTMntVr\nnt9O3fo3za5UGkL557y4sn8dV/nZ1Mzq8mt1z/adxXitN9bCSh+3uZWf/+JjyvvjnC77S61CcN7s\n3pY1a2Z5w63Gkvqnp79GmfnK5uU7I+KLwJOAz09gvtdXPnpJYdqNwMXj3l8CXFKZ7pkTGLYkSTrK\nTSjhiYh5wBuBHwNuAT6amV8e5MAkSZL6ZaLH9i+jc/O/W4CXMa7FhCRJ0lQ30VNaT8/MnwSIiI8C\nXx/ckCRJkvprokd49h18kZn7BzQWSZKkgZjoEZ5nR8QjzesA5jfvJ9pL66gys1JOVKsy2vVYufKl\nVilz4qJyb5/nxuLqmE49rtw/aOuuPcX4rBnlde/eX64Q27WvHH/6kvqusWhBuRqrVslSq66pscCl\nuxm1HZJ6j7Xavl3rHxeVqq6llX5wzz+53ptq8fzyn6u9B8pjfWBX+f9mL1x5fDFeq5Q6Zm79z+Si\nSi+t2ve0q1JhubCybZL6Z0K/ZZlZrrGUJEmaBrwhiSRJaj0THknSwEXEmojYGBEbt22bbEciafI8\ncSxJGrjMXA+sBxgbG6s2ex62VWuvHvUQNCQe4ZEkSa3nEZ4hqvXFqfXdqVXR7NlTrj6p9aYCWDG7\nXKV10sJyr6s9lWqSHY+V1/14pX9Yrc8VwLzZtcqesn2Vapyo9CKbVfn+sjLWzrKOntKubt9DrYBr\nf6Vd2/zKz3Jnpdqr1ufq8Xn1MT1j5pOK8T37yoM6sLi8rMcqfd+OP65c/ditOrBW2VWrNKxVXh5N\n+500Kh7hkSRJrecRHkmSjiK165Y2rTt/yCMZrpEf4YmIt0TErRHxrYj4rcLnvx0RNzePWyPiQEQs\naT7bFBG3NJ9tHP7oJUnSdDDSIzwR8UzgN4Czgb3A5yPi6sy84+A0mfkHwB800/8C8P9m5oPjFvOi\nzNw+xGFLkqRpZtRHeJ4GfDUzdzU9ur4MvLLL9K8HPjGUkUmSpNYY9TU8twLviYjjgd3Ay4HiqamI\nOAY4D3jzuHACX4hOmc7/aO7zMGXVqmJqFR21yqfjKpVPXVojsbtSdbVgbrmya8mCchXNY7UynYpa\nVQrUq6hq89Sq1mp9i9Rdt8qgWvURld5Ytdq6xZWf5d5KpdS+LvtXrQqxtm9XCsSqv2/HVJZf/y7q\n+2q1QrCynMcrg+3W70xSb0aa8GTmbRHxPmADsBP4B6DWjf0XgL875HTWOZm5JSKWAxsi4vbMvP7Q\nGSNiDbAGYMXKlX3dBmkY3Icl6ciM+pQWmfnRzDwzM38GeBC4ozLp6zjkdFZmbmmetwJX0LkWqLSO\n9Zk5lpljy5Yu69/gpSFxH5akIzPyhKc5OkNErAReReEanYh4EvCzwJXjYgsi4riDr4Fz6ZwikyRJ\n+hGjvoYH4NPNNTz7gDdl5kMR8UaAzPxIM80rgS9k5qPj5jsBuKK5DmEW8PHM/PwQxy1JkqaJkSc8\nmfnPC7GPHPL+UuDSQ2J3Ac8e5NgkSTpatP2GhCNPeI4mtaqYWvXWjMr0tcKnWk+uzjzlZfXat6pa\nldKlkqWmVi1T2+5axUrt+7M/Uf/VfjY1u/eVK6hqP/tuZtYqmSpj2lepBFswt/xnr7ZpM7tsc207\nDlTGakWhNDojv4ZHkiRp0Ex4JElS65nwSJKk1jPhkSRJrWfCI0mSWs8qrSmg1345taqk2V17/pQr\nuCqLYm65lVa1X1etKqUW76b2fdT6DVmMNTy97qu9Vu/N7VK9VduV5lX21T2VCrFaz6xa77pu1X61\nfdJqrKmnVnKto4dHeCRJUuuZ8EiSpNYz4ZEkSa1nwiNJklrPhEeSJLWeVVpDNOieT7NqTbaA3XvL\nFSvzZvfWC6i2jhlRqxyrDqnnip9ep9fo1XpNTWafr1VE1czvsRqrtvgZ1Nc7s8vvnEbHiiyVjPy3\nNSLeEhG3RsS3IuK3Cp+/MCIejoibm8fvjfvsvIj4TkTcGRFrhztySZI0XYz0CE9EPBP4DeBsYC/w\n+Yi4OjPvOGTSv83Mnz9k3pnAh4CXApuBGyLiqsz89hCGLkmSppFRn9J6GvDVzNwFEBFfBl4J/H8T\nmPds4M7MvKuZ95PABYAJjyRJfVI6Rbhp3fkjGMmRGfUprVuBn4mI4yPiGODlwIrCdC+IiH+IiM9F\nxDOa2CnAveOm2dzEJEmSfsRIj/Bk5m0R8T5gA7AT+Adg/yGT3QSclpk7I+LlwGeBM4DSVY/Fqwsj\nYg2wBmDFypV9Gr00PO7DknRkRn1Ki8z8KPBRgIj4fTpHasZ//si419dExB9HxNJmuvFHg04FtlTW\nsR5YD3DWWWO9N3fqk35VY01mObVqrNqyasUnvVbKWFnVH1NlH+610rBf+zzU96XamGZW1l3bhef0\n2PdLOpr1Ugk3VU5/jTzhiYjlmbk1IlYCrwJecMjnJwL3Z2ZGxNl0TsM9AHwfOCMiTge+B7wO+OXh\njl6SNCqWn6sXI094gE9HxPHAPuBNmflQRLwRIDM/Arwa+M2I2A/sBl6Xnf/S7Y+INwPXAjOBSzLz\nW6PZBEmSNJWNPOHJzH9eiH1k3Os/Av6oMu81wDWDG50kaSrwaI6O1MgTHkmS1F61ZHXY1/aY8EiS\npgyP5Bw9hn3hc9QqHNoqIrYBdw9xlUuB7UNc36C5PT90WmYu6+dgJmIE+/BU0bZ9bxAm8x0NZT8e\nf2sF4KnAd3pcxHT4+U+HMcL0GGcvY5zQPnzUJTzDFhEbM3Ns1OPoF7dHo+LP6vDa/B1Nh22bDmOE\n6THOQYxx1HdaliRJGjgTHkmS1HomPIO3ftQD6DO3R6Piz+rw2vwdTYdtmw5jhOkxzr6P0Wt4JElS\n63mER5IktZ4JjyRJaj0THg1M9LNVtqSetel3MCIWjHoMmt5MeEYsIn42In5p1OPol4g4OSJWAGQL\nLhBr289nOouIF0bEayPil0c9lqkqIs6JiFdFxOugHb+DABHxL4C3R8T8UY9lIiJiWUScfEhsSiWf\nEfGiiHjhqMfRTUScFhE/fkhs0t+jCc8IRcTPAR8A7h31WPohIl4FfAn4aET8ZUQ8NyLmjnhYk9a2\nn890FvH/t3fm0XJVVfWYqnkAABIqSURBVB7+fiSBBELCJAityEKQUQYZQpgJIgI2SkTGoBDQRkAc\numWQQSDYMjSDzYxtNyABBxCCYANpEBlUmkFImKEdgGZqbAQihISXn3/sU6R4vKleXt6tqre/td56\ndeueqrv3ubfO3XefffbWdsCVwMrANySd3/mGMtSRtDOxsmVd4DhJp1Us0oAgaSfgVGC67Terlqc3\nJO0O3ABMkzRF0lYQxmezGD2SPgGcDcypWpbuKP14LXCZpDMlTYIF68dcpVURxbK+DtjC9kxJo4F5\ntt+oVrL+IWlZ4oZ0tO37ymA7BpgG3GK7aX9YXdFu56eVKYPbqcDzts+SNBL4AZF2/p9tvyhJ7eLN\n6A+SVid+f1+zfaekVYDvAQcBL7dq30haA3gQOND2VEnLA4sDo20/VK1076WMg9OAQ4EXgcOAkcCv\nbf+sStlqFAPy+8DHbT9WvGZqprGtTF9eBxwJPALsDWwIPGX77P5+b3p4qqMDeAVYoXhBfgRcLuma\nMni1GrOBxYD3Adg+AngS+DTwEWg+l24vtNv5aVnKzfp+YA1JK9ieDXwRWAH4dl2boc5pxdgZBrwO\nLA8sX983LfYbhNDjXGCcpM2BK4BjgVskfblSybpmGDEOzrb9AuEhfhoYL2mzSiWbz5LASsBL5Vq5\nDPippFMlbVitaO8gYAQwrBhiPwFuAlaVtEd/vzQNnoqwfQcwCTiHeBL4T+BA4HngzApF6xe2/0oM\nRhtJWrW8dwbwJnBC2W6Zm1I5P/vRJuenFZH0QUmLlSfQ3xAD9XqSRpVB8ADiRrhrpYJWiKSVJY0A\n/mj7J+Xtebb/DPwP8EZptz601m8QwPZzhKdqFjFdPs32QcAuwMlNZEQAYPsl4GrgQEkrlfNwRdm9\nc3WSzadcJ5OJsexx4L+Ao4DRxD2pcmzPIh4yvynpw7ZfB34FPAaM6+/3psEziEiaIOkESZMkfbTc\nVA8FTrR9nu1XbB8CjJD0/orF7RVJu5RYiuPLk8H1RIzFJyV9GMD214FRkga9qnijFH0uKPqMs307\ncAhwQiuen1ZG0i6EkXkO8O9ErMGVwNeArSStWOI5biG8cUOO0ke/AM4nvI9rll3Dy/9lgCVK7MNV\nrfAbrKfmjbL9DKHjrrbPKdOX9xLXQzOe+1+W/3sVo+dlwmjbTtJyFcpV36eXEEbPFbYvsj0T+A7x\nwFrpdVLnhbwGmAF8tRg9rwE/BDaT9KH+fPfw3pskA0EJEjsDuJCIbbla0oG2b5V0e127fYClCM9I\n01KerM4EvgssTbgbdwcuAfYH3i/pPmAU8CGaODgO3qPPUsD1kg6yPa0Vz0+rUga7DwCnEPEPjwJf\nAO4GxgMXUJ5CJf0vsA/wb5UIWxHd9NEk4FZJO9h+uDR9gfCuvh/4jO3/q0DchijTxWOB35W3OgBs\nPyvphfLa5Xe4FdA0gdmShtnusH13MWy2AY6UdDGwNjFNM7tSIcvl42Bqp31bAvOoWMaaF9L285Km\nEWERZ0k6lXigHkFMdTZMGjyDx2bA6bYvg3dWnVwraVfbd0kaTgRmHQXsYfvVCmXtC39HBOJdAiDp\nD8TN6B+A4wj37cHAXGC/FtXnu5Letn1DmTbYi9Y5Py1JGeyekfQb4AngJdunSZoL/Jr4Hd0PbAKs\nD2xv+4nKBK6AbvrojNJHN0uaYPtxIgZta+BTth+rUOQ+IekzwInAU8CzwOOSLrX9V0mL2H5b0qLE\nDfBbwJ62n65Q3nFEQPIbtu+x3SFphO25Zcx4HphAjItzgcPLVE2VMs4rcTsdndodCBwO7FumjwZT\nxjHFe/MebD8o6TlgInA0YTQebPv/+3Uw2/m3kP/KSToduKjuvaOIZYF3E09rY4gVFWtWLW8fdfoY\n8WS9Yt17uwEvARuU7ZHAElXLugD6fLrosxGwaDk/a1Uta7v+AX8PfJ14gvsR8K1O+48mprcWq1rW\nJu6jI4BLieDZnYBVq5a5j3otS0xhrl22JwP3EAHKYzq1HQ+sUrG8OxGLMi4mlk7/oG7fYp3aLlfF\nONiLjMPrXo8u19Q6Fcg4kViFNw5YpNO+ztujgUUX6HhVXjTt/ke4CLcrr5chnlyuAC4npoCGAecB\n65c2w6uStY/67ES4xiFWIvyYmAZahPkpDo4Avlq1rAOozzdr+jT7+WnlP+ATwAPAjmV7FWJ1y5F1\nbVYpg7eqlrfJ++j7VcvaD93GAncAE+reu4qIfdm7bG9OLKWuWtZhhLG5X9keA9wJXNWp3ZadjZ8m\nlHE8MKIiGVcpMk0vsm7c1W8b2BEYORDHzKDlhYCC5YlcBz+UtLPDBbcRYWn/DNjFdgdhta4MYPvt\nqmTujRKDdCbwFwDbbwFfInIjnEXE6UAYDitXIWMjNKDPyNrrZj4/rUxZbvxD4Eu2byrxD88CnwG+\nJukbimyr2xKeuKUqE7YiGuyjDSUtU520jeOYIp4KHCBpP0nfIWJJHiEMPYhx5eFuvmLQKOP27+q2\nX7O9JZHC4iJ4J4/MtoR3p5ll/DgllUgFzAOOsb0DcZ6PJ4Kmhxf5asHLmxNxaAtMJh5ciEg6BliC\nWEI5xfZVnfYfQExt7eAK56J7owy2lwFftj1d0pJE4q/ny5Lhc4h4sCWBNYm59aZLClaj3fRpdRTJ\n5W4hVizeSTzZv03c3F4HVgVeI54AJztWlAwphkIfSRoLfIq4Cf/FscITSTcQ08sdrvCGJekjLvFi\nZeXbUcDOtbG7GKEXEjGMjxGek0FdrNGCMo4txi6SjgM2BU6yfY9iJfOAXscZtLwQkLQIYMI78ARw\nFzClLNXusP0vkjYlXHV7NLOxU1iGGEyfK6sozgaQ9BKRDfOLwOrAasCjtv9QlaB9pN30aWlsP16W\nWF9DxEqdSGRSPogITD7K9jOSlrb9SoWiVsZQ6KOal0fSlbbnAUj6PLEKdKQHOeC3HkmfAn4i6Trb\ne9m+vBihd0nawvbTtl+W9DYwthhmg21ItJKM02zvbftVSYvanmN7SjF6vi7paWAXSds7chsNzPHT\nwzPwSJHmXpGbZkfbp0g6iYhvOc328WXVz+JukdU+kvYmlsCOJebVbyRWgOwAfNMtsOS1nnbTpx2Q\ntDYR83Ze3Xs3EeVK7q/9rqqTsHqGUh9Jmgz8E+FhrcxjVaZ+riZCETYn4nL2LvumALsSeYKWI9ID\n7DzYD0ktKuNw27UUE4uVsAIk3UZk599xoM97GjwDhKT1iP58sO69DQkX9M+JfBHXAnsAh9r+RSWC\n9pFu9JkIrGT73LK9JJHy+5Bm94K0mz5DAUmfJZYf72z7xarlaUbauY8UyeVG2H6qCWRZifAKjySm\nhObWGRS7ETEmGwFnVzX93aIyzq4ZPWX/R4jFI/vXj9UDdvw0eBac4qa7lFhWeYHtu+r2fY+wqA+y\nfY2kTwMzbf++Gml7pxd9hpWAuFrejCOIlU4D5nYcaNpNn3anBCseQDzdf87zE+klheyj6lAUCL0Y\nmGN7b0nrALNs/6li0d6hxWR80/YkSRsQq8kecWSnHvhjpsGzYCgSYZ1NJPn6PRF0dVntpippAvCq\no4J407ube9Onrt3hRG2pfZs5oLfd9BkKlJv5NsALboGEeVWQfVQtJfj3dGJqZhiwre1nq5Xq3bSY\njOMJGbdx1E9bKGTQ8gJie46kY4G3CJfh0sDnJQ23/Svbt1YrYWP0oM8ijtpfNV4E9mn2J8t202co\nUB4KbqtajmYm+6haSvDvDCKX1w7NZkhAy8q40IwdSA9Pvynut7cAbD9a9/7qxBLK1YFTieWizzhS\nvTctDerztJs8nX+76ZMkSfMgaWki3u8fbc+oWp6uSBm7OF4aPI0jaSdi7nEakVzqDNv/Ubd/NWLF\nzyFEHpcNbT9Zgah9okF91iJKR6Q+SZIMWSSNtF11MdAeSRnfTU5pNUCZN18C+Aqx0uo6RZXty8uy\nugsBbD9VllQuB4xr1ptpP/XZNPVJkmSo0+yGBKSMnUmDpwHKvPksSfcCYxSVcX8raS/gp5Jm275E\nUY12TWBiM8eEpD7NrU+SJEkycGQtrf7xArA9MArA9r3AfsBhklaz3WF7ou37qxSyAVKfZKEhqUPS\nA5IekvRTSYtXLROApG8NwHecLukxSTMkXSNpyNX5SpJWIQ2eBihTJtg+H1gcuFDS2OJJuBOYQdS3\naQlSn2SQeNP2BrbXJVLZH9zXDxZv3MKiYYOnC3mmA+vaXo8oI3P0QAiWJMnAkwZPL0haQ9J4RSmI\nd/rL9p5l+2xgsqRDibwYTX1DTX2aW58hwB1EjTIkXSvpPkkPS/pSrYGkWZJOknQ3MF7S8ZLuKR6i\ni2uGraTbJJ0l6XZJj0raRNLPJD0p6eS675sk6b+Ll+kiScMknQKMKu9N7a5dV/LUK2P7Ztu1a+q3\nwAcWXtclSbIgpMHTA4rSA9OAk4lCfYdKGlPbb3svYgB/H7EaaNdmzHVQI/Vpbn3aHUnDiXwbtfo4\nk21vRFT4PlyReRUi8Pwh2+OKZ+5c25sUD9EooqJ2jTm2tybS1E8jSrmsC+wvaVlJawF7AlvY3gDo\nIJJLHsV8z9O+3bXrRp7umExk807ajDaflp1SpmQfkHSzovxDW5LL0ruheAwuB/7V9l2KmjWbEbld\nTnenop+qK37WjKQ+za1POyOpg/lGzh1Ezo05kk4Adivvr0IUC/ytoprzYp5f8uOzRMmPxYlK9+c4\nCvLeBhxTzv8EooDmDuUztwOHA1sSU1e1UiGjgCttnyBplu3Rpf1hPbR7lzzd6HgMYbhNdA6qbUen\na2UqcJ/tM/v42WE9XTsDJVcDn3mXPJLG2H6tvD4cWNt2n6edW4n08PTMGCJBHcA1wPXAokCtINum\nkj5W9s8ZfPEaJvVJqqDmSdnA9leKsbMt8HFgvO31gd8RBQUhCgrWjJ2RRJXn3W1/FPh+XTsoySWB\neXWva9vDAQGX1h1/DdsndCFjT+1m92LsfIHwOu2bxs6QoN2mZV+r21wCaNtrOA2ebrA9FzgTmChp\nK9vzgDuBB4CtJY0CtgCeK+2b+iJJfZpbnyHIWOAV229IWpPwznVFzbh5WdJoYPcGj3MLsLuk5QEk\nLaOowg0wt3gKe2vXLZI+CRxJTJe+0aBsSYvRrtOykr4j6ZnS/vgB6KqmJA2enrkDuBnYT9LWZTnz\nFcBKwEq2z7L9QrUiNkTqkzQLNwLDFXV0phABv+/B9l8Ir85M4FrgnkYOYvsR4Fjg5nKs6cCKZffF\nwAxJU3tp1xPnAksC08tT9YWNyJe0DKMkPQDcCzxNxAxCGDkPEtfvB5nvce4Arq77/HaS7pY0E5gA\nrFO377ryfybwsO3ny/T778t3bg9sBNxTZNieKInTmZ7adZbnXdg+xvYHganAYT32RAuTiQd7wPbs\n4io0cHR5En2LCIKdValw/SD1SaqgqxiDMqDv1Jf2to8ljJHO7bate30bdcU0O+37MfDjLj5/JOGd\n6a1dtzEStlfrbl/SVrxZvCbv0Gla9o0SU9bTtOzGtp8psWv9mZbtLeVBT+16nJat4wrgBuDbfWjb\ncqSHpxdsv0I8YZ5GWObbAZNsv1ipYP0k9UmSJBkQ2mVadvW6zV2BxxqUr2VID08fsD0H+GVZ+eES\nL9KypD5JkiQLzI3AwWUK9HF6mJaVVJuW/SP9mJaVVJtuXQSYS8T5/In507L3lzie7tr1xCmS1iA8\nSn+igcSgrUYuS0+SJEmSpO3JKa0kSZIkSdqeNHiSJEmSJGl70uBJkiRJkqTtSYMnSZIkSZK2Jw2e\nJkftXbTucyUl+zxJGw+EXEmSJEnSFWnwND+11OHrEvWg+rxksFZHZSHRsMHThTwPAROB2wdEoiRJ\nkiTphjR4Wot2K1r3qO3HF3anJUmSJEkaPC1CuxatS5IkSZLBIDMtNz+1onUQHp76onW7lde1onV/\npuuidUcAiwPLAA8DPy/73lO0DkBSrWjdlswvRgdhML3UhYzb99Cux6J1SZIkSTIYpMHT/AyVonVJ\nkiRJstDIKa3WpC2K1iVJkiTJYJEGT2tyIzC8FK2bQg9F64hK4jOBa+lH0TqgVoxuBjAdWLHsrhWt\nm9pLu26RtJukZ4lg5hsk3dSIfEmSJEnSV7J4aJIkSZIkbU96eJIkSZIkaXvS4EmSJEmSpO1JgydJ\nkiRJkrYnDZ4kSZIkSdqeNHiSJEmSJGl70uBJkiRJkqTtSYMnSZIkSZK2Jw2eJEmSJEnanr8BXFmV\nm8py84EAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import sys\n", + "import os\n", + "from __future__ import print_function\n", + "import pints\n", + "import pints.toy as toy\n", + "import pints._diagnostics as diagnostics\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# Load a forward model\n", + "model = toy.LogisticModel()\n", + "\n", + "# Create some toy data\n", + "real_parameters = [0.015, 500]\n", + "times = np.linspace(0, 1000, 1000)\n", + "org_values = model.simulate(real_parameters, times)\n", + "\n", + "# Add noise\n", + "noise = 10\n", + "values = org_values + np.random.normal(0, noise, org_values.shape)\n", + "real_parameters = np.array(real_parameters + [noise])\n", + "\n", + "# Get properties of the noise sample\n", + "noise_sample_mean = np.mean(values - org_values)\n", + "noise_sample_std = np.std(values - org_values)\n", + "\n", + "# Create an object with links to the model and time series\n", + "problem = pints.SingleSeriesProblem(model, times, values)\n", + "\n", + "# Create a log-likelihood function (adds an extra parameter!)\n", + "log_likelihood = pints.UnknownNoiseLogLikelihood(problem)\n", + "\n", + "# Create a uniform prior over both the parameters and the new noise variable\n", + "prior = pints.UniformPrior(\n", + " [0.01, 400, noise*0.1],\n", + " [0.02, 600, noise*100]\n", + " )\n", + "\n", + "# Create a Bayesian log-posterior (prior * likelihood)\n", + "log_posterior = pints.LogPosterior(prior, log_likelihood)\n", + "\n", + "# Create differential evolution MCMC routine\n", + "x0 = real_parameters * 1.1\n", + "mcmc = pints.DreamMCMC(log_posterior, x0)\n", + "\n", + "# Use 4000 iterations in total\n", + "mcmc.set_iterations(4000)\n", + "\n", + "# Discard the first 2000 iterations as burn in\n", + "mcmc.set_burn_in(2000)\n", + "\n", + "# Store only every 4th sample\n", + "mcmc.set_thinning_rate(4)\n", + "\n", + "# Disable verbose mode\n", + "mcmc.set_verbose(False)\n", + "\n", + "# Set gamma\n", + "mcmc.set_gamma(0.1)\n", + "\n", + "# Run!\n", + "print('Running...')\n", + "chains = mcmc.run()\n", + "print('Done!')\n", + "\n", + "# Plot output\n", + "import pints.plot # Use Pints' diagnostic tool\n", + "\n", + "stacked = np.vstack(chains)\n", + "pints.plot.pairwise(stacked, kde=False)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[1.0227840549752976, 1.0210350415772438, 1.0086711799402364]" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "diagnostics.rhat_all_params(chains)" + ] + } + ], + "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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/pints/__init__.py b/pints/__init__.py index ee06aa414..1c89e3492 100644 --- a/pints/__init__.py +++ b/pints/__init__.py @@ -90,6 +90,8 @@ def version(formatted=False): # from ._mcmc import MCMC from ._mcmc._adaptive import AdaptiveCovarianceMCMC, adaptive_covariance_mcmc +from ._mcmc._differentialEvolution import DifferentialEvolutionMCMC, differential_evolution_mcmc +from ._mcmc._differentialEvolution import DreamMCMC from ._mcmc._result import McmcResultObject # diff --git a/pints/_mcmc/_differentialEvolution.py b/pints/_mcmc/_differentialEvolution.py new file mode 100644 index 000000000..d56678251 --- /dev/null +++ b/pints/_mcmc/_differentialEvolution.py @@ -0,0 +1,578 @@ +# +# Differential evolution MCMCs +# +# This file is part of PINTS. +# Copyright (c) 2017, University of Oxford. +# For licensing information, see the LICENSE file distributed with the PINTS +# software package. +# +# Some code in this file was adapted from Myokit (see http://myokit.org) +# +from __future__ import absolute_import, division +from __future__ import print_function, unicode_literals +import pints +import numpy as np + + +class DifferentialEvolutionMCMC(pints.MCMC): + """ + *Extends:* :class:`MCMC` + + Uses differential evolution MCMC as described in [1] + to do posterior sampling from the posterior. + + In each step of the algorithm N chains are evolved + using the evolution equation, + + x_proposed = x[i,r] + gamma * (X[i,r1] - x[i,r2]) + epsilon + + where r1 and r2 are random chain indices chosen (without + replacement) from the N available chains, which must not + equal i or each other, where i indicates the current time + step. epsilon ~ N(0,b) in d dimensions (where d is the + dimensionality of the parameter vector). + + If x_proposed / x[i,r] > u ~ U(0,1), then + x[i+1,r] = x_proposed; otherwise, x[i+1,r] = x[i]. + + [1] "A Markov Chain Monte Carlo version of the genetic + algorithm Differential Evolution: easy Bayesian computing + for real parameter spaces", 2006, Cajo J. F. Ter Braak, + Statistical Computing. + """ + def __init__(self, log_likelihood, x0, sigma0=None): + super(DifferentialEvolutionMCMC, self).__init__( + log_likelihood, x0, sigma0) + + # Total number of iterations + self._iterations = self._dimension * 2000 + + # Gamma + self._gamma = 2.38 / np.sqrt(2 * self._dimension) + + # Normal proposal std. + self._b = 0.01 + + # Number of chains to evolve + self._num_chains = 100 + + # Number of iterations to discard as burn-in + self._burn_in = int(0.5 * self._iterations) + + # Thinning: Store only one sample per X + self._thinning_rate = 1 + + def burn_in(self): + """ + Returns the number of iterations that will be discarded as burn-in in + the next run. + """ + return self._burn_in + + def iterations(self): + """ + Returns the total number of iterations that will be performed in the + next run, including the non-adaptive and burn-in iterations. + """ + return self._iterations + + def run(self): + """See: :meth:`pints.MCMC.run()`.""" + + # Report the current settings + if self._verbose: + print('Running differential evolution MCMC') + print('gamma = ' + str(self._gamma)) + print('normal proposal std. = ' + str(self._b)) + print('Total number of iterations: ' + str(self._iterations)) + print( + 'Number of iterations to discard as burn-in: ' + + str(self._burn_in)) + print('Storing one sample per ' + str(self._thinning_rate)) + + # Initial starting parameters + mu = self._x0 + current = self._x0 + current_log_likelihood = self._log_likelihood(current) + if not np.isfinite(current_log_likelihood): + raise ValueError( + 'Suggested starting position has a non-finite log-likelihood.') + + # chains of stored samples + chains = np.zeros((self._iterations, self._num_chains, + self._dimension)) + current_log_likelihood = np.zeros(self._num_chains) + + # Set initial values + for j in range(self._num_chains): + chains[0, j, :] = np.random.normal(loc=mu, scale=mu / 100.0, + size=len(mu)) + current_log_likelihood[j] = self._log_likelihood(chains[0, j, :]) + + # Go! + for i in range(1, self._iterations): + for j in range(self._num_chains): + r1, r2 = R_draw(j, self._num_chains) + proposed = chains[i - 1, j, :] \ + + self._gamma * (chains[i - 1, r1, :] # NOQA + - chains[i - 1, r2, :]) \ + + np.random.normal(loc=0, scale=self._b * mu, # NOQA + size=len(mu)) + u = np.log(np.random.rand()) + proposed_log_likelihood = self._log_likelihood(proposed) + + if u < proposed_log_likelihood - current_log_likelihood[j]: + chains[i, j, :] = proposed + current_log_likelihood[j] = proposed_log_likelihood + else: + chains[i, j, :] = chains[i - 1, j, :] + + # Report + if self._verbose and i % 50 == 0: + print('Iteration ' + str(i) + ' of ' + str(self._iterations)) + print(' In burn-in: ' + str(i < self._burn_in)) + + non_burn_in = self._iterations - self._burn_in + chains = chains[non_burn_in:, :, :] + chains = chains[::self._thinning_rate, :, :] + + # Convert 3d array to list of lists + samples = [chains[:, i, :] for i in range(0, self._num_chains)] + + # Return generated chain + return samples + + def set_burn_in(self, burn_in): + """ + Sets the number of iterations to discard as burn-in in the next run. + """ + burn_in = int(burn_in) + if burn_in < 0: + raise ValueError('Burn-in rate cannot be negative.') + self._burn_in = burn_in + + def set_iterations(self, iterations): + """ + Sets the total number of iterations to be performed in the next run + (including burn-in and non-adaptive iterations). + """ + iterations = int(iterations) + if iterations < 0: + raise ValueError('Number of iterations cannot be negative.') + self._iterations = iterations + + def set_thinning_rate(self, thinning): + """ + Sets the thinning rate. With a thinning rate of *n*, only every *n-th* + sample will be stored. + """ + thinning = int(thinning) + if thinning < 1: + raise ValueError('Thinning rate must be greater than zero.') + self._thinning_rate = thinning + + def thinning_rate(self): + """ + Returns the thinning rate that will be used in the next run. A thinning + rate of *n* indicates that only every *n-th* sample will be stored. + """ + return self._thinning_rate + + def set_gamma(self, gamma): + """ + Sets the gamma coefficient used in updating the position of each + chain. + """ + if gamma < 0: + raise ValueError('Gamma must be non-negative.') + self._gamma = gamma + + def set_b(self, b): + """ + Sets the normal scale coefficient used in updating the position of each + chain. + """ + if b < 0: + raise ValueError('normal scale coefficient must be non-negative.') + self._b = b + + def set_num_chains(self, num_chains): + """ + Sets the number of chains to evolve + """ + if num_chains < 10: + raise ValueError('This method works best with many chains (>>10,' + + 'typically).') + self._num_chains = num_chains + + +def differential_evolution_mcmc(log_likelihood, x0, sigma0=None): + """ + Runs an differential evolution MCMC routine with the default parameters. + """ + return DifferentialEvolutionMCMC(log_likelihood, x0, sigma0).run() + + +class DreamMCMC(pints.MCMC): + """ + *Extends:* :class:`MCMC` + + Uses differential evolution adaptive Metropolis (DREAM) + MCMC as described in [1] to do posterior sampling + from the posterior. + + In each step of the algorithm N chains are evolved + using the following steps: + + 1. Select proposal: + + x_proposed = x[i,r] + (1 + e) * gamma(delta, d, p_g) * + sum_j=1^delta (X[i,r1[j]] - x[i,r2[j]]) + + epsilon + + where [r1[j], r2[j]] are random chain indices chosen (without + replacement) from the N available chains, which must not + equal each other or i, where i indicates + the current time step; + delta ~ uniform_discrete(1,D) determines + the number of terms to include in the summation; + e ~ U(-b*, b*) in d dimensions; + gamma(delta, d, p_g) =: + if p_g < u1 ~ U(0,1): + 2.38 / sqrt(2 * delta * d) + else: + 1 + + epsilon ~ N(0,b) in d dimensions (where + d is the dimensionality of the parameter vector). + + 2. Modify random subsets of the proposal according to + a crossover probability CR: + + for j in 1:N: + if 1 - CR > u2 ~ U(0,1): + x_proposed[j] = x[j], + else: + x_proposed[j] = x_proposed[j] from 1. + + If x_proposed / x[i,r] > u ~ U(0,1), then + x[i+1,r] = x_proposed; otherwise, x[i+1,r] = x[i]. + + [1] "Accelerating Markov Chain Monte Carlo Simulation by + Differential Evolution with Self-Adaptive Randomized Subspace + Sampling ", 2009, Vrugt et al., International Journal of + Nonlinear Sciences and Numerical Simulation. + """ + def __init__(self, log_likelihood, x0, sigma0=None): + super(DreamMCMC, self).__init__( + log_likelihood, x0, sigma0) + + # Total number of iterations + self._iterations = self._dimension * 2000 + + # Normal proposal std. + self._b = 0.01 + + # b* distribution for e ~ U(-b*, b*) + self._b_star = 0.01 + + # Probability of longer gamma versus regular + self._p_g = 0.2 + + # Determines maximum delta to choose in sums + self._D = 3 + + # Constant crossover probability boolean + self._constant_crossover = False + + # Crossover probability for variable CR case + self._nCR = 3 + + # Constant CR probability for constant CR probability + self._CR = 0.5 + + # Number of chains to evolve + self._num_chains = 100 + + # Number of iterations to discard as burn-in + self._burn_in = int(0.5 * self._iterations) + + # Thinning: Store only one sample per X + self._thinning_rate = 1 + + def burn_in(self): + """ + Returns the number of iterations that will be discarded as burn-in in + the next run. + """ + return self._burn_in + + def iterations(self): + """ + Returns the total number of iterations that will be performed in the + next run, including the non-adaptive and burn-in iterations. + """ + return self._iterations + + def run(self): + """See: :meth:`pints.MCMC.run()`.""" + + # Report the current settings + if self._verbose: + print('Running differential evolution MCMC') + print('gamma = ' + str(self._gamma)) + print('normal proposal std. = ' + str(self._b)) + print('Total number of iterations: ' + str(self._iterations)) + print( + 'Number of iterations to discard as burn-in: ' + + str(self._burn_in)) + print('Storing one sample per ' + str(self._thinning_rate)) + + # Initial starting parameters + mu = self._x0 + current = self._x0 + current_log_likelihood = self._log_likelihood(current) + if not np.isfinite(current_log_likelihood): + raise ValueError( + 'Suggested starting position has a non-finite log-likelihood.') + + # chains of stored samples + chains = np.zeros((self._iterations, self._num_chains, + self._dimension)) + current_log_likelihood = np.zeros(self._num_chains) + + # Set initial values + for j in range(self._num_chains): + chains[0, j, :] = np.random.normal(loc=mu, scale=mu / 100.0, + size=len(mu)) + current_log_likelihood[j] = self._log_likelihood(chains[0, j, :]) + + # Go! + p = np.repeat(1.0 / self._nCR, self._nCR) + L = np.zeros(self._nCR) + Delta = np.zeros(self._nCR) + after_burn_in_indicator = 0 + if self._constant_crossover is False: + for i in range(1, self._iterations): + # Burn-in + if i < self._burn_in: + for j in range(self._num_chains): + # Step 1. Select (initial) proposal + delta = int(np.random.choice(self._D, 1)[0] + 1) + dX = 0 + u1 = np.random.rand() + if self._p_g < u1: + gamma = 2.38 / np.sqrt(2 * delta * self._dimension) + else: + gamma = 1.0 + e = np.random.uniform(low=-self._b_star * mu, + high=self._b_star * mu) + for k in range(0, delta): + r1, r2 = R_draw(j, self._num_chains) + dX += (1 + e) * gamma * (chains[i - 1, r1, :] - + chains[i - 1, r2, :]) + proposed = chains[i - 1, j, :] + dX \ + + np.random.normal(loc=0, scale=self._b * mu, size=len(mu)) # NOQA + + # Select CR from multinomial distribution + m = np.nonzero(np.random.multinomial(self._nCR, p))[0][0] # NOQA + CR = float(m + 1) / float(self._nCR) + L[m] += 1 + + # Step 2. Randomly set elements of proposal to original + for d in range(0, self._dimension): + u2 = np.random.rand() + if 1.0 - CR > u2: + proposed[d] = chains[i - 1, j, d] + + # Accept/reject + u = np.log(np.random.rand()) + proposed_log_likelihood = self._log_likelihood(proposed) # NOQA + + if u < proposed_log_likelihood - current_log_likelihood[j]: # NOQA + chains[i, j, :] = proposed + current_log_likelihood[j] = proposed_log_likelihood + else: + chains[i, j, :] = chains[i - 1, j, :] + + # Update CR distribution + for d in range(0, self._dimension): + Delta[m] += (chains[i, j, d] - chains[i - 1, j, d])**2 / np.var(chains[:, j, d]) # NOQA + for k in range(0, self._nCR): + p[k] = i * self._num_chains * (Delta[k] / float(L[k])) / np.sum(Delta) # NOQA + p = p / np.sum(p) + + # After burn-in + else: + if after_burn_in_indicator == 0: + if self._verbose: + print('Finished warm-up...starting sampling') + print('Crossover probabilities = ' + str(p)) + after_burn_in_indicator = 1 + for j in range(self._num_chains): + # Step 1. Select (initial) proposal + delta = int(np.random.choice(self._D, 1)[0] + 1) + dX = 0 + u1 = np.random.rand() + if self._p_g < u1: + gamma = 2.38 / np.sqrt(2 * delta * self._dimension) + else: + gamma = 1.0 + e = np.random.uniform(low=-self._b_star * mu, + high=self._b_star * mu) + for k in range(0, delta): + r1, r2 = R_draw(j, self._num_chains) + dX += (1 + e) * gamma * (chains[i - 1, r1, :] - + chains[i - 1, r2, :]) + proposed = chains[i - 1, j, :] + dX \ + + np.random.normal(loc=0, scale=self._b * mu, size=len(mu)) # NOQA + + # Step 2. Randomly set elements of proposal to original + # Select CR from multinomial distribution using tuned p + m = np.nonzero(np.random.multinomial(self._nCR, p))[0][0] # NOQA + CR = float(m + 1) / float(self._nCR) + for d in range(0, self._dimension): + u2 = np.random.rand() + if 1.0 - CR > u2: + proposed[d] = chains[i - 1, j, d] + + # Accept/reject + u = np.log(np.random.rand()) + proposed_log_likelihood = self._log_likelihood(proposed) # NOQA + + if u < proposed_log_likelihood - current_log_likelihood[j]: # NOQA + chains[i, j, :] = proposed + current_log_likelihood[j] = proposed_log_likelihood + else: + chains[i, j, :] = chains[i - 1, j, :] + # Report + if self._verbose and i % 50 == 0: + print('Iteration ' + str(i) + ' of ' + + str(self._iterations)) + print(' In burn-in: ' + str(i < self._burn_in)) + + # Constant crossover probability + else: + CR = self._CR + for j in range(self._num_chains): + # Step 1. Select (initial) proposal + delta = int(np.random.choice(self._D, 1)[0] + 1) + dX = 0 + u1 = np.random.rand() + if self._p_g < u1: + gamma = 2.38 / np.sqrt(2 * delta * self._dimension) + else: + gamma = 1.0 + e = np.random.uniform(low=-self._b_star * mu, + high=self._b_star * mu) + for k in range(0, delta): + r1, r2 = R_draw(j, self._num_chains) + dX += (1 + e) * gamma * (chains[i - 1, r1, :] - + chains[i - 1, r2, :]) + proposed = chains[i - 1, j, :] + dX \ + + np.random.normal(loc=0, scale=self._b * mu, size=len(mu)) # NOQA + + # Step 2. Randomly set elements of proposal to original + for d in range(0, self._dimension): + u2 = np.random.rand() + if 1.0 - CR > u2: + proposed[d] = chains[i - 1, j, d] + + # Accept/reject + u = np.log(np.random.rand()) + proposed_log_likelihood = self._log_likelihood(proposed) + + if u < proposed_log_likelihood - current_log_likelihood[j]: + chains[i, j, :] = proposed + current_log_likelihood[j] = proposed_log_likelihood + else: + chains[i, j, :] = chains[i - 1, j, :] + + # Report + if self._verbose and i % 50 == 0: + print('Iteration ' + str(i) + ' of ' + str(self._iterations)) # NOQA + print(' In burn-in: ' + str(i < self._burn_in)) + + non_burn_in = self._iterations - self._burn_in + chains = chains[non_burn_in:, :, :] + chains = chains[::self._thinning_rate, :, :] + + # Convert 3d array to list of lists + samples = [chains[:, i, :] for i in range(0, self._num_chains)] + + # Return generated chain + return samples + + def set_burn_in(self, burn_in): + """ + Sets the number of iterations to discard as burn-in in the next run. + """ + burn_in = int(burn_in) + if burn_in < 0: + raise ValueError('Burn-in rate cannot be negative.') + self._burn_in = burn_in + + def set_iterations(self, iterations): + """ + Sets the total number of iterations to be performed in the next run + (including burn-in and non-adaptive iterations). + """ + iterations = int(iterations) + if iterations < 0: + raise ValueError('Number of iterations cannot be negative.') + self._iterations = iterations + + def set_thinning_rate(self, thinning): + """ + Sets the thinning rate. With a thinning rate of *n*, only every *n-th* + sample will be stored. + """ + thinning = int(thinning) + if thinning < 1: + raise ValueError('Thinning rate must be greater than zero.') + self._thinning_rate = thinning + + def thinning_rate(self): + """ + Returns the thinning rate that will be used in the next run. A thinning + rate of *n* indicates that only every *n-th* sample will be stored. + """ + return self._thinning_rate + + def set_gamma(self, gamma): + """ + Sets the gamma coefficient used in updating the position of each + chain. + """ + if gamma < 0: + raise ValueError('Gamma must be non-negative.') + self._gamma = gamma + + def set_b(self, b): + """ + Sets the normal scale coefficient used in updating the position of each + chain. + """ + if b < 0: + raise ValueError('normal scale coefficient must be non-negative.') + self._b = b + + def set_num_chains(self, num_chains): + """ + Sets the number of chains to evolve + """ + if num_chains < 10: + raise ValueError('This method works best with many chains (>>10,' + + 'typically).') + self._num_chains = num_chains + + +def R_draw(i, num_chains): + r_both = np.random.choice(num_chains, 2, replace=False) + r1 = r_both[0] + r2 = r_both[1] + while(r1 == i or r2 == i or r1 == r2): + r_both = np.random.choice(num_chains, 2, replace=False) + r1 = r_both[0] + r2 = r_both[1] + return r1, r2