{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "# Monod Example\n", "Author(s): Paul Miles | Date Created: June 18, 2018\n", "\n", "This example is from P.M. Berthouex and L.C. Brown: \"Statistics for Environmental Engineers\", CRC Press, 2002.\n", "\n", "We fit the Monod model\n", "$$y(t;q) = q_1 \\frac{t}{q_2 + t} + \\epsilon, \\quad \\epsilon\\sim N(0,\\sigma^2), \\quad q = [\\mu_{max}, K_x]$$\n", "to observations:\n", "\n", "| x (mg/L COD) | y (1/h) |\n", "|----------|:---------|\n", "| 28 | 0.053|\n", "| 55 | 0.060|\n", "| 83 | 0.112|\n", "| 110| 0.105|\n", "| 138| 0.099|\n", "| 225| 0.122|\n", "| 375| 0.125|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Import required packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9.0\n" ] } ], "source": [ "import numpy as np\n", "import scipy.optimize\n", "from pymcmcstat.MCMC import MCMC\n", "import matplotlib.pyplot as plt\n", "import pymcmcstat\n", "print(pymcmcstat.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Initialize MCMC object, and create data structure" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Initialize MCMC object\n", "mcstat = MCMC()\n", "# Next, create a data structure for the observations and control\n", "# variables. Typically one could make a structure |data| that\n", "# contains fields |xdata| and |ydata|.\n", "ndp = 7\n", "x = np.array([28, 55, 83, 110, 138, 225, 375]) # (mg / L COD)\n", "x = x.reshape(ndp, 1) # enforce column vector\n", "y = np.array([0.053, 0.060, 0.112, 0.105, 0.099, 0.122, 0.125]) # (1 / h)\n", "y = y.reshape(ndp, 1) # enforce column vector\n", "# data structure \n", "mcstat.data.add_data_set(x,y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define model and sum-of-squares function\n", "For the MCMC run we need the sum of squares function. For the plots we shall also need a function that returns the model. Both the model and the sum of squares functions are easy to write as follows" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def modelfun(x, theta):\n", " return theta[0]*x/(theta[1] + x)\n", "\n", "def ssfun(theta,data):\n", " res = data.ydata[0] - modelfun(data.xdata[0], theta)\n", " return (res ** 2).sum(axis=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Perform initial optimization study for estimating the error variance and computing the covariance matrix\n", "We define a residuals function and minimize it using scipy's optimize least-squares." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def residuals(p, x, y):\n", " return y - modelfun(x, p)\n", "\n", "theta0, ssmin = scipy.optimize.leastsq(\n", " residuals,\n", " x0=[0.15, 100],\n", " args=(mcstat.data.xdata[0].reshape(ndp,),\n", " mcstat.data.ydata[0].reshape(ndp,)))\n", "n = mcstat.data.n[0] # number of data points in model\n", "p = len(theta0); # number of model parameters (dof)\n", "ssmin = ssfun(theta0, mcstat.data) # calculate the sum-of-squares error\n", "mse = ssmin/(n-p) # estimate for the error variance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Jacobian matrix of the model function is easy to calculate so we use it to produce estimate of the covariance of theta. This can be used as the initial proposal covariance for the MCMC simulation by assigning it to `qcov` in the `simulation_options` (see below)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "tcov = [[2.44707212e-04 2.50115314e-01]\n", " [2.50115314e-01 3.20837579e+02]]\n" ] } ], "source": [ "J = np.array([[x/(theta0[1]+x)], [-theta0[0]*x/(theta0[1]+x)**2]])\n", "J = J.transpose()\n", "J = J.reshape(n,p)\n", "tcov = np.linalg.inv(np.dot(J.transpose(), J))*mse\n", "print('tcov = {}'.format(tcov))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Setup parameters, model settings, and simulation options\n", "- Parameters: \n", " - $\\mu_{max} \\in [0,\\infty]$\n", " - $K_x \\in [0,\\infty]$\n", "- Model Settings: \n", " - $\\sigma^2_0 = 0.01^2$ (initial error variance)\n", " - sum-of-squares function defined above\n", "- Simulation Options:\n", " - Number of simulation: 5000\n", " - Update error variance\n", " - Initial covariance matrix from least squares optimization" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "mcstat.parameters.add_model_parameter(\n", " name='$\\mu_{max}$',\n", " theta0=theta0[0],\n", " minimum=0)\n", "mcstat.parameters.add_model_parameter(\n", " name='$K_x$',\n", " theta0=theta0[1],\n", " minimum=0)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "mcstat.simulation_options.define_simulation_options(\n", " nsimu=5.0e3,\n", " updatesigma=1,\n", " qcov=tcov)\n", "mcstat.model_settings.define_model_settings(\n", " sos_function=ssfun,\n", " sigma2=0.01**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Run Simulation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Sampling these parameters:\n", " name start [ min, max] N( mu, sigma^2)\n", "$\\mu_{max}$: 0.15 [ 0.00e+00, inf] N( 0.00e+00, inf)\n", " $K_x$: 49.05 [ 0.00e+00, inf] N( 0.00e+00, inf)\n", " [-----------------100%-----------------] 5000 of 5000 complete in 1.2 sec" ] } ], "source": [ "mcstat.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Extract results:\n", "- Display chain statistics" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "------------------------------\n", "name : mean std MC_err tau geweke\n", "$\\mu_{max}$: 0.1565 0.0360 0.0031 34.2737 0.9768\n", "$K_x$ : 65.6914 47.4967 3.9469 32.1899 0.8879\n", "------------------------------\n", "==============================\n", "Acceptance rate information\n", "---------------\n", "Results dictionary:\n", "Stage 1: 21.40%\n", "Stage 2: 51.36%\n", "Net : 72.76% -> 3638/5000\n", "---------------\n", "Chain provided:\n", "Net : 72.76% -> 3638/5000\n", "---------------\n", "Note, the net acceptance rate from the results dictionary\n", "may be different if you only provided a subset of the chain,\n", "e.g., removed the first part for burnin-in.\n", "------------------------------\n" ] } ], "source": [ "results = mcstat.simulation_results.results\n", "names = results['names']\n", "chain = results['chain']\n", "s2chain = results['s2chain']\n", "names = results['names'] # parameter names\n", "mcstat.chainstats(chain, results)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot mean model response and compare with data" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEJCAYAAABVFBp5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5xW897/8den6UzCdFDNqNEUSoddk0hqkwgph9iFrWjTZrPdbPf+seW4HW4bOdxSUilEJRvdKjrbhDQldFBGx7HpMHSi0zSf3x/XavY0ZqZpmmvWNde8n4/Hesy6vuu7rut9LZrPrNN3mbsjIiJSEpXCDiAiIuWXioiIiJSYioiIiJSYioiIiJSYioiIiJSYioiIiJRYVIuImfUwsxVmlmFmdxawvIuZLTKzbDPrU8Dyo8ws08yei2ZOEREpmagVETNLAIYC5wMtgH5m1iJft3XAAOC1Qt7m78C/opVRREQOT+UovvepQIa7rwIws/FAb2DZ/g7uviZYlpN/ZTNrD9QH3gPSDvZhderU8SZNmpRGbhGRCmPhwoWb3b1uSdePZhFpBKzP8zoT6FicFc2sEvAkcDVwTnHWadKkCenp6YeaUUSkQjOztYezfqyeWL8JmOrumUV1MrMbzCzdzNI3bdpURtFERGS/aO6JfAck53mdFLQVx+nAmWZ2E3AkUNXMdrj7ASfn3X0EMAIgLS1Ng4CJiJSxaBaRBUAzM0shUjz6AlcWZ0V3v2r/vJkNANLyFxAREQlf1IqIu2eb2c3A+0ACMNrdl5rZg0C6u082sw7AW8AxwEVm9oC7tyytDHv37iUzM5Ndu3aV1luWS9WrVycpKYkqVaqEHUVE4ozFy1DwaWlpnv/E+urVq6lVqxaJiYmYWUjJwuXuZGVlsX37dlJSUsKOIyIxxswWuvtBr4AtTKyeWC8Vu3btqtAFBMDMSExMrPB7YyISHXFdRIAKXUD20zYQkWiJ5ol1ERGJEbt372bTpk1s3LiRjRs35s4fLhWRMnb//fdz5JFHcscddxS4/O2336Z58+a0aJF/hBgRkf9wd7Zv386GDRsOmDZu3HjAz/3z27Zti0oOFZFA2kMz2Lxjz6/a6xxZlfTB3cssx9tvv03Pnj1VREQqqN27d7Nhwwa+//57fvjhB3744YcD5n/44Yfc4rBz584C3yMxMZH69etTv3592rdvT7169ahfvz5169alXr16uVPdunU5+uijDyuvikigoAJSVPuhePjhhxk7diz16tUjOTmZ9u3b8+KLLzJixAj27NlDamoqr7zyCosXL2by5Ml88MEHPPTQQ7z55pvMnj37V/1q1qx52JlEpGzt2rWL77//nn//+9+50/7X33//fe70448/Frh+3bp1qV+/Pscddxypqam58/uLxf6pbt26VK5cdr/aVUSibOHChYwfP57FixeTnZ1Nu3btaN++PZdeeinXX389AIMHD2bUqFHccsst9OrVi549e9KnT2Rk/KOPPrrAfiISG3Jycti8eTOZmZl89913uVNmZuYBBaOg4lClShUaNGhAgwYNaNasGV27ds19fdxxx+VO9erVi9n7vFREouzDDz/kkksuyd176NWrFwBLlixh8ODBbNmyhR07dnDeeecVuH5x+4lI6cvJyWHjxo2sW7eOzMzMA6b169fnFo69e/cesF6lSpU47rjjaNSoEampqXTp0oWGDRseMDVo0CAubkFQEQnJgAEDePvtt2nTpg1jxoxh7ty5h9VPRA7djh07WLduHWvXrmXdunW50/r163MLR/4CUa1aNZKSkkhKSuKMM86gUaNGJCUlHfCzfv36ZXpIKUwV41uGqEuXLgwYMIC77rqL7Oxs/u///o9Bgwaxfft2GjRowN69exk3bhyNGjUCoFatWmzfvj13/cL6icjBbdmyhTVr1hwwrV27NnfKf4gpISGBRo0akZyczGmnnUZycjLHH388ycnJJCcnk5SURJ06dcr93kNpUhEJ1DmyaqFXZx2Odu3a8bvf/Y42bdpQr149OnToAMDf//53OnbsSN26denYsWNu4ejbty/XX389zz77LJMmTSq0n4jAzp07WbNmDatWrWL16tW5P1evXs2aNWvYunXrAf2POOIImjRpQuPGjenYsSONGzfm+OOPz/3ZsGFDEhISQvo25VNcj521fPlyTj755JASxRZtCymP3J3Nmzfz7bff/mpatWoV33///QH9a9SoQUpKSu7UpEmTA6Zjjz1WexH5HO7YWdoTEZFQuTsbN24kIyODb775hm+++SZ3PiMj41d7340aNaJp06acd955NG3alJSUFE444QRSUlKoX7++ikQZUxERkTKxfft2Vq5cmTutWLEidz5voUhISCAlJYXU1FTOOOMMUlNTadq0aW7BqF69eojfQvJTERGRUuPuZGZm8vXXXx8wLV++/IBDT2ZG48aNad68Of3796dZs2a5U+PGjWP2nogwxMpoGoVRERGRQ5aTk8Pq1atZtmxZ7rR06VK+/vprfv7559x+tWvX5qSTTuLcc8/lxBNP5MQTT6R58+Y0bdqUGjVqhPgNyo9ojqZRGlRERKRQ7s769etZsmQJX331FUuWLMktFnnHbWrUqBEtWrRg4MCBnHzyyZx00kmcdNJJOkdRAaiIiAgAW7du5csvv8ydlixZwpIlSw4Y/TUpKYmWLVty1lln0aJFi9ypdu3aISaXMKmIlCNNmjQhPT2dOnXqHFYfKb9K4/i4u7N69WoWL17MF198kTutWbMmt88xxxxDq1at+P3vf88pp5zCKaecQsuWLTnmmGNK66tInFARESlHDvX4eHZ2Nl9//TWff/75AdP+m/AqVarEiSeeSMeOHbnhhhto06YNrVu3plGjRjoMJcWiIhJla9asoUePHpx22ml8/PHHdOjQgWuvvZb77ruPjRs3Mm7cOFJTU7nuuutYtWoVNWvWZMSIEbRu3ZqsrCz69evHd999x+mnn07eG0NfffVVnn32Wfbs2UPHjh15/vnndadtBbdv3z6+/vprFi5cSHp6Ounp6SxevDj33EWNGjVo3bo1/fr1o127drRt25ZTTjlFJ7hjXLRG0ygtUS0iZtYDeAZIAEa6+//kW94FeBpoDfR190lBe1tgGHAUsA942N0nHE6W//qv/2Lx4sWH8xa/0rZtW55++umD9svIyOCNN95g9OjRdOjQgddee42PPvqIyZMn88gjj5CcnMxvfvMb3n77bWbPns0111zD4sWLeeCBB+jcuTP33nsvU6ZMYdSoUUDk7vMJEyYwb948qlSpwk033cS4ceO45pprSvX7SflSu3bt3CujjjjiCNq3b88f//hH2rVrR7t27WjevHmFGRQwnsTCZbxFidr/UWaWAAwFugOZwAIzm+zuy/J0WwcMAPI/K/YX4Bp3/8bMGgILzex9d98SrbzRlJKSQqtWrQBo2bIl3bp1w8xo1apV7oBwb775JgBnn302WVlZbNu2jX/961/885//BODCCy/MPR49a9YsFi5cmDsO186dO6lXr14I30zKytatW/nss8+K7DNw4EDS0tJIS0ujefPm2jOVMhHNP0tOBTLcfRWAmY0HegO5RcTd1wTLcvKu6O4r88z/28w2AnWBEheR4uwxREu1atVy5ytVqpT7ulKlSmRnZx/yjVXuTv/+/Xn00UdLNafEBnfn22+/5eOPP2bevHl8/PHHLF26FHen8f97t9D1nnnmmTJMKRJRKYrv3QhYn+d1ZtB2SMzsVKAq8G0By24ws3QzS9+0aVOJg4btzDPPZNy4cQDMnTuXOnXqcNRRR9GlSxdee+01AKZNm8ZPP/0EQLdu3Zg0aRIbN24E4Mcff2Tt2rXhhJfDtnfvXubPn88TTzzBxRdfTP369WnWrBn9+/dnwoQJJCUl8cADDzB9+nQSjyj4D45YOT4uFU9MHyA1swbAK0B/d8/Jv9zdRwAjIDKKbxnHKzX3338/1113Ha1bt6ZmzZqMHTsWgPvuu49+/frRsmVLOnXqxPHHHw9AixYteOihhzj33HPJycmhSpUqDB06lMaNG4f5NaSYdu7cyfz58/nXv/7Fhx9+yMcff8wvv/wCQLNmzbjwwgvp1KkTnTp14uSTT6ZSpf/8rdc9tg+PSwUUtaHgzex04H53Py94fReAu//qGIyZjQHe3X9iPWg7CpgLPJK3vTAaCr5o2hbh2bVrF59++ilz5sxhzpw5fPrpp+zduxczo02bNnTp0oUzzzyTM888k/r164cdVyqYWB4KfgHQzMxSgO+AvsCVxVnRzKoCbwEvF6eAiMSS7OxsPvvsM+bMmcPs2bP5+OOP2bVrF5UqVaJ9+/bcdtttdO3alU6dOnH00UeHHVfksEStiLh7tpndDLxP5BLf0e6+1MweBNLdfbKZdSBSLI4BLjKzB9y9JXAF0AVINLMBwVsOcPfSvUZXpBTsPxE+ffp0ZsyYwezZs3OHCmnTpg033ngjZ511Fl26dNHwIBJ3onpOxN2nAlPztd2bZ34BkFTAeq8Cr5ZShgp/5228PL0ylmzdupVZs2Yxffp0pk+fzurVqwFo3Lgxffv2pXv37vz2t7/V8DMS92L6xPrhql69OllZWSQmJlbYQuLuZGVl6UE+h8ndWb58OVOnTmXKlCl89NFHZGdnU6tWLc4++2zuuOMOunfvTmpqaoX9f00qprguIklJSWRmZlKeL/8tDdWrVycp6Vc7fHIQu3btYvbs2UyZMoWpU6fmDlDYqlUr/vKXv3DBBRdw+umn6wFKUqHFdRGpUqUKKSkpYceQcuSnn35iypQpvPPOO0ybNo2ff/6ZmjVrcs4553DnnXdywQUXkJycHHZMkZgR10Wkoor1x2nGmvXr1/P222/z9ttv88EHH7Bv3z4aNGjA1VdfTe/evTnrrLN0OFCkECoicSjWH6cZC9atW8ekSZN44403+PTTT4HITZx//etfufjii0lLSzvgJj8RKZiKiFQYa9euzS0c8+fPB+A3v/kNjz76KJdeeinNmzcPOaFI+aMiInFt06ZNTJw4kXHjxvHJJ58A0K5dOx599FH69OlDampqyAlFyjcVEYk7v/zyC++88w7jxo3j/fffJzs7m1atWvHII49wxRVX0LRp07AjisQNFRGJCzk5OcyZM4exY8fy1ltvsWPHDpKSkrj99tu56qqraN26ddgRReKSikgcivXHaZamdevWMWbMGF566SXWrFlD7dq16du3L1dddRVdunTRyXGRKFMRiUPxfhnv7t27eeeddxg1ahQzZszA3enWrRuPPPIIF198sZ4ZLlKGVESk3Pj2228ZPnw4L730EllZWSQnJ3PPPfcwYMAA3VQqEhIVEYlp+/btY9q0aQwdOpT33nuPypUrc/HFF3P99dfTrVs3PUdcJGQqInJIyupu+E2bNjFq1CiGDx/O2rVradiwIffffz/XX389DRs2LLXPEZHDoyIihyTad8MvW7aMIUOG8Oqrr7J7927OOussnnzySXr16qWBDkVikIqIhM7dmT17Nk8++STTpk2jRo0aXHfdddxyyy16pK9IjFMRkdDs2bOHCRMmMGTIEBYvXky9evV48MEHufHGG/UwJ5FyQkVEytzOnTsZOXIk//jHP8jMzKRFixaMHDmSq666SqPlipQzKiJSZn7++WdeeOEFHn/8cX744Qc6d+7MCy+8QI8ePXRToEg5pSIih6Qkd8Nv376d559/nieffJJNmzZx9tln8/rrr9O1a1c9SlaknItqETGzHsAzQAIw0t3/J9/yLsDTQGugr7tPyrOsPzA4ePmQu4+NZlYpnkO5jPfnn3/m6aefZsiQIfz444+cd9553HPPPZxxxhlRTCgiZSlqRcTMEoChQHcgE1hgZpPdfVmebuuAAcAd+dY9FrgPSAMcWBis+1O08krp2bNnDyNHjuTBBx9kw4YN9OzZk3vuuYdTTz017GgiUsqiuSdyKpDh7qsAzGw80BvILSLuviZYlpNv3fOAGe7+Y7B8BtADeD2KeeUw5eTkMGHCBAYPHsyqVavo0qULb731Fqeffnqx1tdjfUXKn2iezWwErM/zOjNoi/a6Usbcnffee4/27dtz5ZVXUqtWLaZOncrcuXOLXUBAj/UVKY/K9SUxZnaDmaWbWfqmTZvCjlMhLVmyhO7du3P++eezbds2xo0bx6JFizj//PN10lykAohmEfkOSM7zOiloK7V13X2Eu6e5e1rdunVLHFQO3U8//cSf//xn2rZty6JFi3j22WdZvnw5V155pS7XFalAovmvfQHQzMxSzKwq0BeYXMx13wfONbNjzOwY4NygTUK2b98+RowYQbNmzRg6dCiDBg3im2++4ZZbbqFq1fh76JWIFC1qRcTds4GbifzyXw5MdPelZvagmfUCMLMOZpYJXA68YGZLg3V/BP5OpBAtAB7cf5JdwvPRRx/RoUMHBg0aRMuWLVm0aBFDhw4lMTEx7GgiEpKo3ifi7lOBqfna7s0zv4DIoaqC1h0NjI5mPimerKws/vKXvzB27FiSk5OZMGECl19+eamf86hIj/UViRe6Y10K5e5MnDiRP//5z/z444/cdddd3H333RxxxBFR+TxdxitS/qiISIEyMzO58cYbeffdd+nQoQMzZsygdevWYccSkRijy2jkADk5OTz//PO0aNGC2bNnM2TIED755BMVEBEpkPZEJFdGRgYDBgxg3rx5dO/enRdeeIGUlJSwY4lIDNOeiODujB49mrZt27J06VLGjBnD+++/rwIiIgelPZEKLisri0GDBvHmm29y1lln8fLLL5OUVOAFcyIiv6I9kQps1qxZtG7dmsmTJ/OPf/yDmTNnqoCIyCFREamAdu/ezR133ME555zDUUcdxfz58/nv//5vDVciIodMh7MqmIyMDC6//HIWL17MTTfdxOOPP07NmjXDjiUi5ZSKSAUyZcoUrrrqKhISEpg8eTIXXXRR2JFEpJzT8YsKICcnhwceeICePXtywgknsHDhQhUQESkV2hOJc1u2bOHqq69mypQp9O/fn2HDhlGjRo2wY4lInFARiWNfffUVl1xyCWvXrmXo0KHceOONelCUiJQqFZE4NX78eAYOHEjt2rX54IMP6NSpU9iRRCQO6ZxInHF3Hn74Yfr160e7du1YtGiRCoiIRI32ROJIdnY2N910Ey+++CJXX301o0aN0tMGRSSqtCcSJ3bs2EHv3r158cUX+dvf/sbLL7+sAiIiUac9kTiwYcMGLrzwQj7//HOGDx/OoEGDwo4kIhWEikg5t2LFCs4//3w2bNjAO++8Q8+ePcOOJCIViIpIOTZv3jx69epFQkICc+fOpUOHDmFHEpEKJqrnRMysh5mtMLMMM7uzgOXVzGxCsHy+mTUJ2quY2Vgz+8rMlpvZXdHMWR7NmTOHc889l8TERD755BMVEBEJRdSKiJklAEOB84EWQD8za5Gv20DgJ3dPBZ4CHgvaLwequXsroD0waH+BkcgQ7hdeeCEpKSl8+OGHNG3aNOxIIlJBRXNP5FQgw91XufseYDzQO1+f3sDYYH4S0M0it1Q7cISZVQZqAHuAbVHMWm7MnDmTnj170rRpU2bPnk39+vXDjiQiFVg0i0gjYH2e15lBW4F93D0b2AokEikoPwPfA+uAJ9z9xyhmLRemT5/ORRddRPPmzZk9ezb16tULO5KIVHCxep/IqcA+oCGQAvzFzE7I38nMbjCzdDNL37RpU1lnLFPvvfcevXr14sQTT2TWrFnUrVs37EgiIlEtIt8ByXleJwVtBfYJDl3VBrKAK4H33H2vu28E5gFp+T/A3Ue4e5q7p8XzL9WpU6fSu3dvTj75ZGbNmkWdOnXCjiQiAkS3iCwAmplZiplVBfoCk/P1mQz0D+b7ALPd3YkcwjobwMyOAE4Dvo5i1pg1bdo0LrnkEk455RRmzZpFYmJi2JFERHJFrYgE5zhuBt4HlgMT3X2pmT1oZr2CbqOARDPLAG4H9l8GPBQ40syWEilGL7n7l9HKGqvmz5/PZZddximnnMLMmTM59thjw44kInIAi/zhX4yOkT2CXe6+L7qRSiYtLc3T09PDjlFqVq5cSadOnTj66KP5+OOPdRJdRKLCzBa6+69OFxRXoXsiZlbJzK40sylmtpHI4aTvzWyZmT1uZqkl/VAp2g8//ECPHj2oVKkS7733ngqIiMSsog5nzQGaAncBx7l7srvXAzoDnwKPmdnVZZCxQtm+fTsXXnghGzZs4N133yU1VbVaRGJXUWNnnePue/M3BvdrvAm8aWZVopasAtq7dy99+vThiy++YPLkyZx66qlhRxIRKVKhRSRvAQmGMKmft7+7ryuoyEjJuDt/+MMfmD59OqNGjeKCCy4IO5KIyEEddBRfM7sFuA/YAOQEzQ60jmKuCmfw4MG8/PLLPPjgg1x33XVhxxERKZbiDAV/K3Ciu2dFO0xFNXLkSB555BFuuOEGBg8eHHYcEZFiK859IuuJjGklUbBgwQL+9Kc/ce655zJ06FAi40+KiJQPhe6JmNntwewqYK6ZTQF271/u7kOinC3ubd68mcsuu4wGDRrw2muvUbmynhEmIuVLUb+1agU/1wVT1WCSUrBv3z769evHxo0bmTdvnoYzEZFyqagisheY5u6fl1WYiuS+++5j5syZjBw5kvbt24cdR0SkRIoqIt8Ct5pZG+ALYBow3d1/KpNkcWzy5Mk8/PDDDBw4kIEDB4YdR0SkxIo1dpaZ/QboAZwLJAAziQzV/ll04xVfeRk7KyMjg7S0NFJTU/noo4+oXr162JFEpAI73LGzinUmNzik9TnwqJkdBXQH/gDETBEpD3755RcuvfRSEhISmDRpkgqIiJR7h3w5kLtvM7Nt7n5DNALFK3dn0KBBLFmyhGnTptGkSZOwI4mIHLaSPk9kVKmmqADGjBnDq6++ygMPPMB5550XdhwRkVJR1H0i+Z9CmLsI0PWoh2DdunXceuutdO3albvvvjvsOCIipaaow1lnAlcDO/K1G6DhZYvJ3Rk4cCA5OTm89NJLVKoUzScSi4iUraKKyKfAL+7+Qf4FZrYiepHiy/Dhw5k5cybDhw8nJSUl7DgiIqWq2I/HjXWxeInvt99+S+vWrencuTPvvfeexsUSkZgTzcfjHvQ3XnH6VFQ5OTlce+21VK5cmZEjR6qAiEhcKvLxuGZ2i5kdn7fRzKqa2dlmNhboH9145dczzzzDhx9+yDPPPENycnLYcUREoqKoItID2Ae8bmb/NrNlZrYK+AboBzzt7mOKenMz62FmK8wsw8zuLGB5NTObECyfb2ZN8ixrbWafmNlSM/vKzMrNnXkrVqzgb3/7Gz179qR/f9VZEYlfRT0edxfwPPB88Cz1OsBOd99SnDcOHqk7lMjd7ZnAAjOb7O7L8nQbCPzk7qlm1hd4DPidmVUGXgV+7+5fmFkikQEhY152djb9+/enRo0ajBgxosDDWGkPzWDzjj2/aq9zZFXSB3cvi5giIqWiWNebuvted/++uAUkcCqQ4e6r3H0PMB7ona9Pb2BsMD8J6BacZzkX+NLdvwg+P8vd9x3CZ4fmiSeeYP78+QwdOpQGDRoU2KegAlJUu4hIrIrmTQuNiDwVcb/MoK3APu6eTeQJiolAc8DN7H0zW2Rmf41izlKzYsUK7rvvPvr06UPfvn3DjiMiEnWx+ii9ykBnoAPwCzAruAxtVt5OZnYDcAPA8ccf/6s3KWu333471atX57nnntPVWCJSIRx0TyS4QuuYErz3d0Dey5KSgrYC+wTnQWoDWUT2Wv7l7pvd/RdgKtAu/we4+wh3T3P3tLp165YgYumZOnUqU6dO5d5776V+/fqhZhERKSvFOZxVn8hJ8YnB1VbF/RN7AdDMzFLMrCrQF8g/Htdk/nOZcB9gtkfufnwfaGVmNYPi0hVYRozas2cPt912G82bN+eWW24JO46ISJk5aBFx98FAMyIj9w4AvjGzR8ys6UHWywZuJlIQlgMT3X2pmT1oZr2CbqOARDPLAG4H7gzW/QkYQqQQLQYWufuUEny/MvG///u/rFy5kqeeeoqqVQ/+GPo6Rxbcp7B2EZFYVexhT4LH5F5L5P6ROcBpwAx3j4mT3mENe7JhwwaaN2/OGWecwdSpU8v880VEDkfUn2xoZrcC1wCbgZHAf7v7XjOrROTGw5goImG5++67+eWXX3jqqafCjiIiUuaKc3XWscCl7r42b6O755hZz+jEKh8WLlzI6NGjue222zjxxBPDjiMiUuY0im9R71nEneUL7j6HM888k5UrV/LNN99Qu3btUv1sEZGyEPXDWRVZUXeWjx8/nnnz5vHiiy+qgIhIhaXH7JXQX//6V9q1a8e1114bdhQRkdBoT6SEMjMzef3110lISAg7iohIaLQnUkJ9+/alc+fOYccQEQmVikgJ/f3vfw87gohI6FREilDYHeRV9u0kNTW1jNOIiMQenRMpQv4HRN1111089thjLFsWs8N4iYiUKe2JFFNWVhbPPfccv/vd7zjppJPCjiMiEhNURIrp6aefZseOHQwePDjsKCIiMUNFpBi2bNnCs88+y2WXXUbLli3DjiMiEjNURIrh2WefZdu2bdoLERHJR0XkILZt28ZTTz1F7969adu2bdhxRERiiorIQTz33HNs2bKFe+65J+woIiIxR0WkCDt27GDIkCFccMEFtG/fPuw4IiIxR0WkCMOGDSMrK0t7ISIihVARKcQvv/zC448/Tvfu3TnttNPCjiMiEpNURArxwgsvsGnTJu69996wo4iIxKyoFhEz62FmK8wsw8zuLGB5NTObECyfb2ZN8i0/3sx2mNkd0cyZX3Z2NkOGDKFr164aqVdEpAhRKyJmlgAMBc4HWgD9zKxFvm4DgZ/cPRV4Cngs3/IhwLRoZSzMlClTyMzM5NZbby3rjxYRKVeiuSdyKpDh7qvcfQ8wHuidr09vYGwwPwnoZmYGYGYXA6uBpVHMWKBhw4bRqFEjLrroorL+aBGRciWaRaQRsD7P68ygrcA+7p4NbAUSzexI4P8BD0QxX4G+/fZb3n//fa6//noqV9YgxyIiRYnVE+v3A0+5+46iOpnZDWaWbmbpm/hkQ1oAAA9/SURBVDZtKpUPfuGFF0hISOAPf/hDqbyfiEg8i+af2t8ByXleJwVtBfXJNLPKQG0gC+gI9DGzfwBHAzlmtsvdn8u7sruPAEYApKWl+eEG3rVrF6NHj6ZXr140apR/p0lERPKLZhFZADQzsxQixaIvcGW+PpOB/sAnQB9gtrs7cOb+DmZ2P7AjfwGJhjfffJOsrCxuvPHGaH+UiEhciFoRcfdsM7sZeB9IAEa7+1IzexBId/fJwCjgFTPLAH4kUmhCM2zYMFJTU+nWrVuYMUREyo2onjl296nA1Hxt9+aZ3wVcfpD3uD8q4fL56quvmDdvHk888QSVKsXqqSIRkdii35aBYcOGUa1aNQYMGBB2FBGRckNFBNi+fTuvvPIKV1xxBYmJiWHHEREpN1REgHHjxrFjxw6dUBcROUQVvoi4O8OGDaNNmzYarVdE5BBV+FuyP/30U7788kuGDx9OMOKKiIgUU4XfExk2bBi1atXiyivz38IiIiIHU6GLSFZWFhMnTuTqq6+mVq1aYccRESl3KnQRGTNmDLt379YJdRGREqrQReSVV16hY8eOtGrVKuwoIiLlUoUtIitWrOCLL76gb99QR1oRESnXKmwReeONNwDo06dPyElERMqvCl1EzjjjDJKSksKOIiJSblXIIvL111/z5ZdfcvnlRY79KCIiB1Ehi4gOZYmIlI4KWUQmTpxI586d9fRCEZHDVOGKyLJly1iyZAlXXHFF2FFERMq9CldE3njjDcyMyy67LOwoIiLlXoUrIhMnTuTMM8+kYcOGYUcRESn3KlQRWbp0KcuWLdOhLBGRUlKhiogOZYmIlK4KU0TcnYkTJ9K1a1eOO+64sOOIiMSFqBYRM+thZivMLMPM7ixgeTUzmxAsn29mTYL27ma20My+Cn6efbhZli5dyvLly3WDoYhIKYpaETGzBGAocD7QAuhnZi3ydRsI/OTuqcBTwGNB+2bgIndvBfQHXjncPBMnTqRSpUpceumlh/tWIiISiOaeyKlAhruvcvc9wHigd74+vYGxwfwkoJuZmbt/7u7/DtqXAjXMrFpJg+hQlohIdESziDQC1ud5nRm0FdjH3bOBrUBivj6XAYvcfXdJg3z11VesWLFCV2WJiJSyymEHKIqZtSRyiOvcQpbfANwAcPzxxxf6Pm+88YYOZYmIREE090S+A5LzvE4K2grsY2aVgdpAVvA6CXgLuMbdvy3oA9x9hLunuXta3bp1Cwyx/1DWWWedRb169Q7n+4iISD7RLCILgGZmlmJmVYG+wOR8fSYTOXEO0AeY7e5uZkcDU4A73X3e4YT48ssvWblypQ5liYhEQdSKSHCO42bgfWA5MNHdl5rZg2bWK+g2Ckg0swzgdmD/ZcA3A6nAvWa2OJhKtBvx1ltvUalSJS655JLD+j4iIvJr5u5hZygVaWlpnp6e/qv2Tp06sW/fPubPnx9CKhGR2GZmC909raTrx/Ud61u3buWzzz6je/fuYUcREYlLcV1E5s6dy759+1RERESiJK6LyIwZM6hZsyann3562FFEROJS3BeRrl27UrVq1bCjiIjEpbgtIuvWrWPlypU6lCUiEkVxW0RmzpwJoCIiIhJFcVtEZsyYwXHHHUfLli3DjiIiErfisojk5OQwc+ZMzjnnHMws7DgiInErLovIF198webNm3UoS0QkyuKyiOw/H3LOOeeEnEREJL7FZRGZMWMGLVu2pGHDhmFHERGJa3FXRHbt2sWHH36ovRARkTIQd0Xko48+YteuXTofIiJSBuKuiMycOZMqVarQtWvXsKOIiMS9uCsiM2bM4PTTT+fII48MO4qISNyLqyKyefNmPv/8c50PEREpI3FVRGbNmoW763yIiEgZiasiMnPmTGrXrk1aWokf0iUiIocgrorIjBkzOPvss6lcuXLYUUREKoS4KSK7d+9m7dq1Oh8iIlKGolpEzKyHma0wswwzu7OA5dXMbEKwfL6ZNcmz7K6gfYWZnXewz9q2bRugod9FRMpS1IqImSUAQ4HzgRZAPzNrka/bQOAnd08FngIeC9ZtAfQFWgI9gOeD9yvUtm3baNy4MampqaX7RUREpFDR3BM5Fchw91XuvgcYD/TO16c3MDaYnwR0s8jY7b2B8e6+291XAxnB+xVq+/btdO/eXUO/i4iUoWgWkUbA+jyvM4O2Avu4ezawFUgs5roH2Ldvn86HiIiUsXJ9GZOZ3QDcAFC9enW6desWciIRkYolmnsi3wHJeV4nBW0F9jGzykBtIKuY6+LuI9w9zd3TWrZsSZ06dUoxvoiIHEw0i8gCoJmZpZhZVSInyifn6zMZ6B/M9wFmu7sH7X2Dq7dSgGbAZ1HMKiIiJRC1w1nunm1mNwPvAwnAaHdfamYPAunuPhkYBbxiZhnAj0QKDUG/icAyIBv4k7vvi1ZWEREpGYv84V/+paWleXp6etgxRETKFTNb6O4lHisqbu5YFxGRsqciIiIiJaYiIiIiJaYiIiIiJaYiIiIiJRY3V2eZ2XZgRdg5iqEOsDnsEMWgnKVLOUtXechZHjICnOjutUq6crke9iSfFYdzmVpZMbN05Sw9ylm6lLP0lIeMEMl5OOvrcJaIiJSYioiIiJRYPBWREWEHKCblLF3KWbqUs/SUh4xwmDnj5sS6iIiUvXjaExERkTIWF0XEzHqY2QozyzCzO8POk5eZrTGzr8xs8f6rIMzsWDObYWbfBD+PCSHXaDPbaGZL8rQVmMsing2275dm1i7knPeb2XfBNl1sZhfkWXZXkHOFmZ1XRhmTzWyOmS0zs6VmdmvQHlPbs4icsbY9q5vZZ2b2RZDzgaA9xczmB3kmBI+YIHhkxISgfb6ZNQk55xgzW51ne7YN2sP8d5RgZp+b2bvB69Lblu5ericiw8x/C5wAVAW+AFqEnStPvjVAnXxt/wDuDObvBB4LIVcXoB2w5GC5gAuAaYABpwHzQ855P3BHAX1bBP/9qwEpwf8XCWWQsQHQLpivBawMssTU9iwiZ6xtTwOODOarAPOD7TQR6Bu0DwduDOZvAoYH832BCWW0PQvLOQboU0D/MP8d3Q68BrwbvC61bRkPeyKnAhnuvsrd9wDjgd4hZzqY3sDYYH4scHFZB3D3fxF5hkteheXqDbzsEZ8CR5tZgxBzFqY3MN7dd7v7aiCDyP8fUeXu37v7omB+O7AcaESMbc8ichYmrO3p7r4jeFklmBw4G5gUtOffnvu38ySgm5lZiDkLE8p/dzNLAi4ERgavjVLclvFQRBoB6/O8zqTofxhlzYHpZrbQIs+EB6jv7t8H8z8A9cOJ9iuF5YrFbXxzcEhgdJ7DgaHnDHb/f0Pkr9KY3Z75ckKMbc/g8MtiYCMwg8he0BZ3zy4gS27OYPlWIDGMnO6+f3s+HGzPp8ysWv6cgbLank8DfwVygteJlOK2jIciEus6u3s74HzgT2bWJe9Cj+w3xtwlcrGaKzAMaAq0Bb4Hngw3ToSZHQm8CfyXu2/LuyyWtmcBOWNue7r7PndvCyQR2fs5KeRIBcqf08xOAe4ikrcDcCzw/8LKZ2Y9gY3uvjBanxEPReQ7IDnP66SgLSa4+3fBz43AW0T+QWzYvxsb/NwYXsIDFJYrpraxu28I/vHmAC/yn0MsoeU0sypEfjGPc/d/Bs0xtz0LyhmL23M/d98CzAFOJ3L4Z/9QTXmz5OYMltcGskLK2SM4bOjuvht4iXC35xlALzNbQ+RQ/9nAM5TitoyHIrIAaBZcbVCVyMmgySFnAsDMjjCzWvvngXOBJUTy9Q+69QfeCSfhrxSWazJwTXB1yWnA1jyHacpcvuPIlxDZphDJ2Te4wiQFaAZ8VgZ5DBgFLHf3IXkWxdT2LCxnDG7PumZ2dDBfA+hO5PzNHKBP0C3/9ty/nfsAs4M9vzByfp3nDwcjcq4h7/Ys0//u7n6Xuye5exMivxtnu/tVlOa2jPZVAWUxEbnqYSWR46Z3h50nT64TiFzd8gWwdH82IscYZwHfADOBY0PI9jqRQxd7iRwTHVhYLiJXkwwNtu9XQFrIOV8JcnwZ/E/fIE//u4OcK4DzyyhjZyKHqr4EFgfTBbG2PYvIGWvbszXweZBnCXBv0H4CkSKWAbwBVAvaqwevM4LlJ4Scc3awPZcAr/KfK7hC+3cUfP5v+c/VWaW2LXXHuoiIlFg8HM4SEZGQqIiIiEiJqYiIiEiJqYiIiEiJqYiIiEiJqYiIELnO38w+MLOEUnzPaWaWZGZzzazIZ22bWXMzm2qREX8XmdlEM6sfLOtskdFivw6mG/Ksl3cE3m/M7J9m1iLP8vFm1qy0vpNIfioiIhHXAf90932l8WbBzWeJ7p5ZjL7VgSnAMHdv5pFhcp4H6prZcURGX/2ju59E5F6PQWZ2YZ63eMrd27p7M2ACMNvM6gbLhhEZN0kkKlREJK6ZWYdgILzqwQgCS4PxjfK7iuCuXTP7bbBX8o6ZrTKz/zGzq4K9ga/MrGnQr6mZfRq0PWRmO/K832+BucWMeSXwibv/3/4Gd5/r7kuAPwFj/D+j724mUhQKfG6Ou08ApgfvCfAhcE6eIS5ESpWKiMQ1d19A5C7sh4g83+PV4JdzrmC4nBPcfU2e5jbAH4GTgd8Dzd39VCLDad8S9HkGeMbdWxG5mz6v84H3ihnzFKCwAfJaFrAsPWgvzCKCAQs9Mh5WBpHvI1LqVESkIniQyLhGaUQKSX51gC352hZ4ZCC93USGqZgetH8FNAnmTycyRAREDjnldQbw0eHFLrH8z3/YCDQMI4jEPxURqQgSgSOJPM2vegHLdxbQvjvPfE6e1zlAkYeGzOwEYL1HHpJWHEuB9oUsW1bAsvbBOoX5DZEBC/erTuQ7ipQ6FRGpCF4A7gHGAY/lX+juPwEJwQnuQ/EpcFkw3zdP+6EcyoLIXkynvCfLzaxLcO5mKDDA/vOc7kQi36GgPSrM7DIio0W/nqe5Of8ZSVakVKmISFwzs2uAve7+GvA/QAczO7uArtOJXPl0KP4LuN3MvgRSiTwFDqAHvy4iU8wsM5jeyLvA3XcCPYFbgst0lxF51vUmjwwVfjXwopl9DXwMjM57Eh64bf8lvkHfs919U/D96wM73f2HQ/xuIsWiUXxFADNrB9zm7r8/hHVqEvkF7WbWF+gHXAHMc/ci7wspK2Z2G7DN3UeFnUXiky77EwHcfZGZzTGzhEO4V6Q98Fzw8KEtwHXBifiYKCCBLUSeFyISFdoTERGREtM5ERERKTEVERERKTEVERERKTEVERERKTEVERERKTEVERERKbH/D+syYcTSM8wtAAAAAElFTkSuQmCC\n", "text/plain": [ "