{ "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": [], "source": [ "import numpy as np\n", "import scipy.optimize\n", "from pymcmcstat.MCMC import MCMC\n", "import matplotlib.pyplot as plt" ] }, { "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.1 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", "name : mean std MC_err tau geweke\n", "$\\mu_{max}$: 0.1573 0.0341 0.0039 73.9652 0.9094\n", "$K_x$ : 68.1782 50.6119 6.0702 88.0618 0.7321\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": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEKCAYAAADTgGjXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3X98z/X+//HbwxojIj9jG5vQB6EY6hQV5VdChxNSiKIcfXR+fCudPqcU1TnnUzk+kfxM6DOFnPlRhoXyQTYpv8KMGDI/0tBhtj2+f7xfW+9mm9n23uu99x7Xy+V92fv9fD1f7/f9/ap57PV6vl7Pl6gqxhhjTGGUczuAMcaY0suKiDHGmEKzImKMMabQrIgYY4wpNCsixhhjCs2KiDHGmEKzImKMMabQfFpERKSbiOwRkUQReT6X5R1FZKuIpItIv1yWXyciySLyji9zGmOMKRyfFRERCQImA92BZsBAEWmWo9shYCjwYR5v8yqw3lcZjTHGFM01PnzvdkCiqiYBiEg00BvYldVBVQ86yzJzriwibYA6wGdA1JU+rGbNmhoREVEcuY0xpsxISEg4qaq1Cru+L4tIKHDY63Uy0L4gK4pIOeBN4BHg3oKsExERQXx8/NVmNMaYMk1Evi/K+v46sD4KWKGqyfl1EpERIhIvIvEnTpwooWjGGGOy+HJP5AgQ7vU6zGkriNuBDiIyCqgMlBeRc6r6q8F5VZ0GTAOIioqymSSNMaaE+bKIbAEai0gknuIxAHi4ICuq6qCs5yIyFIjKWUCMMca4z2dFRFXTRWQ0sBIIAmap6k4ReQWIV9UYEWkLfAJcDzwgIuNUtXlxZbh06RLJyclcuHChuN6yVAoJCSEsLIzg4GC3oxhjAowEyv1EoqKiNOfA+oEDB6hSpQo1atRARFxK5i5V5dSpU5w9e5bIyEi34xhj/IyIJKjqFc+AzYu/DqwXiwsXLpTpAgIgItSoUaPM740ZY3wjoIsIUKYLSBbbBsYYX/HlwLoxxhg/kZ6ezokTJ0hJSeH48eOkpKSQkpJS5Pe1IlLCXn75ZSpXrsyf//znXJcvWbKEJk2a0KxZzhlijDHm17IKww8//JD9OH78OMePH//V8+PHj3Pq1CmfZLAi4ogav4qT59Iua69ZuTzxL95XYjmWLFlCz549rYgYU4adP3+eY8eO5frwLhgnTpwgt5OjKleuTJ06dahTpw433XQTHTt2pE6dOtSuXfuyn9dff32RsloRceRWQPJrvxoTJkxgzpw51K5dm/DwcNq0acP06dOZNm0aaWlpNGrUiLlz57Jt2zZiYmJYt24d48ePZ9GiRcTFxV3Wr1KlSkXOZIwpeRcvXuTo0aMcOXKEo0eP5vk4e/bsZesGBwdzww03ULduXSIiIrjtttu44YYbstvq1KnDDTfcQJ06dUr03wgrIj6WkJBAdHQ027ZtIz09ndatW9OmTRt++9vf8sQTTwDw4osvMnPmTJ5++ml69epFz5496dfPMzN+tWrVcu1njPEvqampJCcnZz8OHz5McnIyR44cyX7kdkipQoUK1KtXj3r16tGyZUu6du1KvXr1qFu3bnaBqFu3LtWrV6dcOf87F8qKiI998cUXPPjgg9l/GfTq1QuAHTt28OKLL3LmzBnOnTtH165dc12/oP2MMb5z/vx5Dh8+nOcjOTk5172H2rVrExoaSv369bn99tsJDQ3NftSrV4/Q0FCuv/76Un0GpRURlwwdOpQlS5bQqlUr3n//fdauXVukfsaYwlFVUlJS+P777/n+++85dOgQhw4d+tXznHsQIsINN9xAeHg4TZs2pUuXLoSFhf3qUa9ePcqXL+/Styo5VkR8rGPHjgwdOpSxY8eSnp7O0qVLGTlyJGfPnqVu3bpcunSJ+fPnExoaCkCVKlV+9RdNXv2MMQWTVSQOHDjAgQMHOHjwIAcPHuT777/P/pnzYtwqVarQoEED6tevz2233UZ4eDj169cnPDyc8PBwQkNDy0SBKAgrIo6alcvneXZWUbRu3Zr+/fvTqlUrateuTdu2bQF49dVXad++PbVq1aJ9+/bZhWPAgAE88cQTTJo0iYULF+bZzxjzi/Pnz3PgwAGSkpLYv38/SUlJvyoaP//886/616xZk4iICG6++WZ69uxJREQEDRo0yC4c1apVc+mblD4BPXfW7t27adq0qUuJ/IttC1OaqSqnT58mMTHxV4+kpCSSkpL44YcfftW/SpUqREZG0rBhQyIjI3/1aNCgAZUrV3bpm/ifos6dZXsixhi/8eOPP7Jv3z727t3L3r172bdvX3bBOHPmTHY/ESEsLIwbb7yRHj160LBhQ2688cbsn9WrVy/Vg9WliRURY0yJSktLY//+/Xz33Xfs2bOHPXv2ZBeNkydPZvcrV64c9evXp3HjxgwcOJBGjRrRqFEjGjduTGRkJCEhIS5+C5PFiogxxifOnDnD7t27sx9ZRSMpKYmMjIzsfnXr1qVJkyY8+OCDNGnShMaNG9OkSRMaNmxIhQoVXPwG/sFfZtPIixURY0yRnDp1ip07d7Jjxw527dqVXTSOHTuW3adChQo0adKEVq1a0b9/f2666absx3XXXediev/ny9k0ioMVEWNMgZw7d44dO3awffv27KKxc+fOXw1qV6lSJfu6iWbNmtG0aVOaNm1KZGQkQUFBLqY3vmJFxBjzKxkZGSQmJrJ9+3a+/fZbvv32W7Zv305SUlJ2n2uvvZZmzZrRrVs3br75Zpo3b07z5s0JCwuzAe0yxopIKRIREUF8fDw1a9YsUh9TehX38fHz58+zfft2tm3blv3Yvn179nUV5cqVo0mTJkRFRTFs2DBatGhBixYtaNCggV/O42RKnhURY0qRohwf/+mnn9i6dStbt24lISGBrVu3snfv3uypxKtVq8Ytt9zCE088QatWrWjVqhVNmzalYsWKxfodTGCxIuJjBw8epFu3btx222383//9H23btuWxxx7jpZdeIiUlhfnz59OoUSOGDRtGUlISlSpVYtq0abRs2ZJTp04xcOBAjhw5wu233/6r+wbMmzePSZMmkZaWRvv27ZkyZYodczbZUlNTSUhIYMuWLSQkJJCQkMD+/fuzl4eHh9O6dWsGDhzILbfcwi233EL9+vXtUJQf8tVsGsXFp0VERLoB/wSCgBmq+kaO5R2BiUBLYICqLnTabwHeBa4DMoAJqrqgKFmeeeYZtm3bVpS3uMwtt9zCxIkTr9gvMTGRjz/+mFmzZtG2bVs+/PBDvvzyS2JiYnjttdcIDw/n1ltvZcmSJcTFxTF48GC2bdvGuHHjuPPOO/nrX//K8uXLmTlzJuC5+nzBggVs2LCB4OBgRo0axfz58xk8eHCxfj9TuvzP//wPW7ZsYcuWLezZsyf7j46IiAjatGnD8OHDad26Na1bt6ZWrVoupzUF5Q+n8ebHZ0VERIKAycB9QDKwRURiVHWXV7dDwFAg571ifwYGq+o+EakHJIjISlU9QykUGRlJixYtAGjevDmdO3dGRGjRokX2BHCLFi0CoFOnTpw6dYrU1FTWr1/P4sWLAbj//vuz70C2Zs0aEhISsufh+ve//03t2rVd+GbGn/znf/4nderUoV27djz88MO0a9eOqKgoatSo4XY0E8B8uSfSDkhU1SQAEYkGegPZRURVDzrLMr1XVNW9Xs+PikgKUAsodBEpyB6Dr3hfMFWuXLns1+XKlSM9PZ3g4OCrej9VZciQIbz++uvFmtP4n4sXL5KQkMCXX37Jxo0b4abH8+x76NAhOzvKlDhfnl4RChz2ep3stF0VEWkHlAf257JshIjEi0j8iRMnCh3UbR06dGD+/PkArF27lpo1a3LdddfRsWNHPvzwQwA+/fRTfvzxRwA6d+7MwoULSUlJAeD06dN8//337oQ3xer06dMsW7aMsWPH0qFDB6pWrcodd9zBc889x44dOwjO+Heu69WsXJ7w8HArIKbE+fXAuojUBeYCQ1Q1M+dyVZ0GTAPPLL4lHK/YvPzyywwbNoyWLVtSqVIl5syZA8BLL73EwIEDad68Ob/5zW+oX78+AM2aNWP8+PF06dKFzMxMgoODmTx5Mg0aNHDza5hCOH78OOvXr2fdunWsW7eOHTt2AJ77abdp04ann36aO+64g9/85jd2yNL4JZ9NBS8itwMvq2pX5/VYAFW97BiMiLwPLMsaWHfargPWAq95t+fFpoLPn20L//DDDz+wdu1a1q5dy7p16/juu+8Az8V7d9xxBx07dqRDhw60bdvWTq01JcKfp4LfAjQWkUjgCDAAeLggK4pIeeAT4IOCFBBj/NVPP/3EunXrWLNmDWvWrGHnzp2AZ3qQDh068Nhjj3HXXXfRunXrqx4bM8Yf+KyIqGq6iIwGVuI5xXeWqu4UkVeAeFWNEZG2eIrF9cADIjJOVZsDDwEdgRoiMtR5y6GqWrzn6BpTzC5dusTGjRuJjY1l9erVxMfHk5GRQcWKFbnzzjt59NFH6dSpE7feeivXXOPXR5ONKRCf/l+sqiuAFTna/ur1fAsQlst684B5xZShzA82BsrdK/3VgQMHWLlyJStXrmTNmjWcPXuWoKAg2rVrx9ixY+ncuTO33367TWtuAlJA/ykUEhLCqVOnqFGjRpktJKrKqVOn7AY+xejixYusX7+eZcuW8dlnn7F3r+eM9AYNGjBw4EC6detGp06dqFq1qstJjfG9gC4iYWFhJCcnU5pP/y0OISEhhIVdtsNnrsLx48dZsWIFy5YtIzY2lnPnzhESEsI999zDqFGj6NatG02aNCmzf6yYsiugi0hwcDCRkZFuxzClkKqyc+dOlixZwtKlS/nqq68Azx8mjzzyCD179uSee+6hUqVKLic1xl0BXUTKKn+/naa/yszM5KuvvmLx4sV88sknJCYmAtC+fXteffVVHnjgAVq2bGl7G8Z4sSISgPz9dpr+JD09nXXr1rFo0SKWLFnCsWPHuOaaa+jUqRN/+tOf6N27N3Xr1nU7pjF+y4qIKXMyMzP54osvWLBgAYsWLSIlJYVKlSrRvXt3HnzwQe6//36qVavmdkxjSgUrIqZMUFU2bdrEggUL+Pjjjzl69CgVK1akZ8+e9O/fn+7du9v4hjGFYEXEBLS9e/cyd+5c5s2bx8GDB6lQoQLdu3enf//+9OzZk8qVK7sd0ZhSzYqICTgnT55kwYIFzJ07l82bN1OuXDnuvfdexo0bR58+fbjuuuvcjmhMwLAiEoD8/XaavnDp0iVWrFjB7NmzWb58Oenp6bRs2ZJ//OMfPPzww9SrV8/tiMYEJCsiAagsnca7f/9+Zs6cyezZs/nhhx+44YYbGDNmDI8++iitWrVyO54xAc+KiCl1Lly4wCeffMKMGTOIi4ujXLly3H///Tz++OP06NHDJjY0pgTZb5spNQ4cOMCUKVOYNWsWp0+fJiIigvHjxzN06FBCQ6/6ppnGmGJgRcRclZK+Gj4zM5PVq1fzzjvvsGzZMsqVK0efPn0YOXIknTt3plw5X97h2RhzJVZEzFUpqavhU1NTmTNnDpMnT2bPnj3Url2bv/zlL4wcOdImkzTGj1gRMX7l8OHDTJw4kenTp3P27Fnat2/PvHnz6Nevn92Pwxg/ZEXE+IVvv/2Wf/zjH0RHR6Oq9O/fn2eeeYa2bdu6Hc0Ykw8rIsY1qkpcXBx///vfiY2N5dprr2X06NE888wzNGjQwO14xpgCsCJiSpyqEhMTw6uvvkpCQgJ16tRhwoQJPPXUU1x//fVuxzPGXAU7tcVclbyuei/I1fCqypIlS2jdujV9+vThzJkzTJs2jYMHD/LCCy9YATGmFPLpnoiIdAP+CQQBM1T1jRzLOwITgZbAAFVd6LVsCPCi83K8qs7xZVZTMIU5jTczM5MlS5bwyiuv8M0339CoUSPef/99Bg0aZBcGGlPK+WxPRESCgMlAd6AZMFBEmuXodggYCnyYY93qwEtAe6Ad8JKI2J+ppYyq8sknn3DrrbfSt29ffv75Zz744AN2797NkCFDrIAYEwB8+VvcDkhU1SQAEYkGegO7sjqo6kFnWWaOdbsCq1T1tLN8FdAN+F8f5jXF6IsvvuDZZ59l06ZNNGnShLlz5zJgwIB8C4fd1teY0seXYyKhwGGv18lOm6/XNS7atWsXvXr1omPHjhw6dIgZM2awc+dOHnnkkSvuedhtfY0pfUr1wLqIjBCReBGJP3HihNtxyrTk5GSGDx9OixYtWLduHa+99hr79u1j+PDhdtjKmADmyyJyBAj3eh3mtBXbuqo6TVWjVDWqVq1ahQ5qCu/8+fP813/9F40bN2bevHmMGTOG/fv3M3bsWLvdrDFlgC+LyBagsYhEikh5YAAQU8B1VwJdROR6Z0C9i9Nm/ISqsnDhQpo2bcr48eP57W9/y549e3jrrbeoWbOm2/GMMSXEZ0VEVdOB0Xj+8d8NfKSqO0XkFRHpBSAibUUkGfgd8J6I7HTWPQ28iqcQbQFeyRpkN+7btWsX9913H7/73e+oXr06X3zxBfPnzyciIsLtaMaYEubTg9WqugJYkaPtr17Pt+A5VJXburOAWb7MZ65Oamoq48aNY9KkSVSuXJnJkyczYsSIYhvzKIu39TWmtLMRT3NFqspHH33EmDFjSElJ4fHHH2fChAkU9ziUncZrTOljRcTk68iRI4waNYqYmBiioqJYunSpzaxrjMlWqk/xNb6jqsyYMYPmzZuzatUq/vu//5tNmzZZATHG/IrtiZjLJCUl8cQTTxAXF8fdd9/N9OnTadSokduxjDF+yPZETLaMjAwmTpxIixYt2LJlC1OnTmXNmjVWQIwxebI9EQN4rjgfPHgwn3/+OT169GDq1KmEh4dfeUVjTJlmRcSwePFiHn/8cdLS0pg5cyaPPfYYIuJ2LGNMKWCHs8qw8+fPM3LkSPr27cuNN97I119/zbBhw6yAGGMKzIpIGfX111/Tpk0bpk+fznPPPceGDRto3Lix27GMMaWMFZEyRlV56623aN++PWfPnmX16tW88cYblC9vV4UbY66ejYmUIWfPnuWxxx5j0aJF9OnThxkzZlCjRg23YxljSjErImXEvn376NOnD9999x1vvvkmf/jDH2zswxhTZFZEyoBly5bxyCOPEBwczKpVq+jUqZPbkYwxAcLGRAJYZmYmr7zyCg888AA33ngj8fHxVkCMMcXK9kQC1E8//cTgwYOJiYlh8ODBTJ06lYoVK7odyxgTYKyIBKADBw7Qo0cP9u3bx6RJkxg9erSNfxhjfMKKSIDZunUrPXr0IC0tjdWrV3P33Xe7HckYE8BsTCSAxMbGctddd1GhQgU2bNhgBcQY43NWRALEBx98wP3330/Dhg3ZuHEjTZs2dTuSMaYMsCJSyqkqr732GkOGDOGuu+5i/fr11KtXz+1YxpgywsZESrGMjAxGjx7N1KlTGTRoELNmzbLpS4wxJcqneyIi0k1E9ohIoog8n8vyCiKywFm+WUQinPZgEZkjIttFZLeIjPVlztIoLS2Nhx56iKlTp/Lcc8/xwQcfWAExxpQ4nxUREQkCJgPdgWbAQBFplqPbcOBHVW0EvA38zWn/HVBBVVsAbYCRWQXG/FJAFi9ezMSJE3njjTcoV86OTBpjSp4v/+VpBySqapKqpgHRQO8cfXoDc5znC4HO4rmgQYFrReQaoCKQBqT6MGupkZaWxu9+9zv+9a9/8c477zBmzBi3IxljyjBfFpFQ4LDX62SnLdc+qpoO/ATUwFNQzgPHgEPAf6vq6ZwfICIjRCReROJPnDhR/N/Az1y8eJF+/foRExPD5MmT+f3vf+92JGNMGeevx0DaARlAPSAS+JOINMzZSVWnqWqUqkbVqlWrpDOWqKwCsnTpUqZMmcKoUaPcjmSMMT4tIkeAcK/XYU5brn2cQ1dVgVPAw8BnqnpJVVOADUCUD7P6tawCsmzZMt59912eeuoptyMZYwzg2yKyBWgsIpEiUh4YAMTk6BMDDHGe9wPiVFXxHMLqBCAi1wK3Ad/5MKvfunjxIn379s0uIE8++aTbkYwxJluBi4iIXOuccVUgzhjHaGAlsBv4SFV3isgrItLL6TYTqCEiicAfgazTgCcDlUVkJ55iNFtVvy3oZweKjIwMBgwYwPLly5k6daoVEGOM38nzYkMRKYdn72EQ0Ba4CFQQkZPAcuA9VU3M781VdQWwIkfbX72eX8BzOm/O9c7l1l6WqCpPP/00S5YsYdKkSYwcOdLtSMYYc5n89kQ+B24ExgI3qGq4qtYG7gQ2AX8TkUdKIGOZ9MYbb/Duu+/y7LPP8vTTT7sdxxhjcpXftCf3quqlnI3OqbaLgEUiEuyzZGXY3LlzeeGFFxg0aBCvv/6623GMMSZPeRYR7wLijIXU8e6vqodyKzKmaFatWsWwYcPo1KkTs2bNsivRjTF+7YoTMIrI08BLwHEg02lWoKUPc5VJ27Zto2/fvjRr1ozFixfbXFjGGL9XkFl8xwA3qeopX4cpyw4ePEj37t2pVq0an376KVWrVnU7kjHGXFFBishhPNORGB85ffo03bt358KFC6xZs8buB2KMKTXyO8X3j87TJGCtiCzHc5ovAKr6lo+zlQlZ14IkJSWxatUqmjXLOdGxMcb4r/z2RKo4Pw85j/LOwxSjl156iVWrVjFjxgw6duzodhxjjLkq+RWRS8Cnqvp1SYUpa5YuXcqECRMYPnw4w4cPdzuOMcZctfyKyH5gjIi0Ar4BPgViVfXHEkkW4BITE3n00Udp06YN77zzjttxjDGmUPK7TmQBsABARG4FugGLnWtGVuOZZferEkkZYH7++Wf69u1LUFAQCxcuJCQkxO1IxhhTKAU5OwvnkNbXwOsich1wH/A4YEXkKqkqTz75JNu3b2fFihVERES4HckYYwrtqi+HVtVUIFVVR/ggT8CbOnUqc+fO5eWXX6Zbt25uxzHGmCIp7JwaM4s1RRmxefNmxowZQ48ePXjxxRfdjmOMMUWW33UiOW8glb0Iz33QzVU4efIk/fr1IzQ0lLlz59qcWMaYgJDfmEgH4BHgXI52wXMPdFNAWeMgKSkpbNy4kerVq7sdyRhjikV+RWQT8LOqrsu5QET2+C5S4ImOjmbRokW88cYbtG7d2u04xhhTbMRzS/PSLyoqSuPj492OcZmjR49y8803c9NNN/Hll18SFFTgOwwbY4zPiUiCqkYVdv08D8yLiBTgw6/YpyxTVUaMGMGFCxeYM2eOFRBjTMDJ9/a4IvK0iNT3bhSR8iLSSUTmAEN8G690mz17NsuXL+f111+nSZMmbscxxphil18R6QZkAP8rIkdFZJeIJAH7gIHARFV9P783F5FuIrJHRBJF5PlcllcQkQXO8s0iEuG1rKWIbBSRnSKyXURK1WXd33//Pc888wx33XWX3SPdGBOw8pv25AIwBZji3Eu9JvBvVT1TkDd2pkeZjOfq9mRgi4jEqOour27DgR9VtZGIDAD+BvQXkWuAecCjqvqNiNTAMyFkqZCZmcmwYcPIzMxk9uzZl53OGzV+FSfPpV22Xs3K5Yl/8b6SimmMMUVWoIsVVPWSqh4raAFxtAMSVTVJVdOAaKB3jj69gTnO84VAZ2ecpQvwrap+43z+KVXNuIrPdtW7775LXFwcb775JpGRkZctz62A5NdujDH+ypdXvIXiuStilmSnLdc+qpqO5w6KNYAmgIrIShHZKiLP+jBnsUpMTOTZZ5+lS5cujBhhM8MYYwJbgSZgdME1wJ1AW+BnYI1zGtoa704iMgIYAVC/fv3L3qSkZWRkMHToUIKDg5k5cyZ28poxJtBdcU/EOUPr+kK89xEg3Ot1mNOWax9nHKQqcArPXst6VT2pqj8DK4DLrtJT1WmqGqWqUbVq1SpExOI1e/ZsNmzYwMSJEwkLC3M7jjHG+FxBDmfVwTMo/pFztlVB/7zeAjQWkUgRKQ8MAHLOxxXDL6cJ9wPi1HP140qghYhUcorLXcAu/NiZM2d44YUXuOOOOxgyxM58NsaUDVcsIqr6ItAYz8y9Q4F9IvKaiNx4hfXSgdF4CsJu4CNV3Skir4hIL6fbTKCGiCQCfwSed9b9EXgLTyHaBmxV1eWF+H4lZty4cZw8eZJJkyZd8TBWzcq536o+r3ZjjPFXBZ72xLlN7mN4rh/5HLgNWKWqfjHo7ea0J7t376Zly5Y89thjTJs2zZUMxhhTGEWd9uSKA+siMgYYDJwEZgD/T1UviUg5PBce+kURcYuq8swzz3DttdcyYcIEt+MYY0yJKsjZWdWB36rq996NqpopIj19E6v0WLp0KbGxsUycOBF/GNw3xpiSZLP45veeV7iy/MKFCzRv3pyQkBC2bdtGcHBwsX6+Mcb4ms8PZ5VlV7qy/O233yYpKYnY2FgrIMaYMsnu0VpIR44cYcKECfTp04f77rP5rowxZZMVkUJ6/vnnSU9P580333Q7ijHGuMaKSCHNmzePP/3pTzRs2NDtKMYY4xorIoUUGhrK2LFj3Y5hjDGusiKSj7yuIM849yMTJkygcuXKJZzIGGP8i52dlY+cN4hSVdq0acPZs2cZNHG3S6mMMcZ/WBG5CsuWLePrr79m9uzZXHONbTpjjLHDWQWkqowbN46GDRvyyCOPuB3HGGP8gv05XUArVqwgISGBmTNn2l6IMcY4bE+kALL2QiIjI3n00UfdjmOMMX7D/qQugM8++4wtW7Ywffp0m97EGGO82J7IFWTthTRo0IDBgwe7HccYY/yK7YlcQWxsLJs3b+a9996jfHm786AxxnizPZF8ZO2FhIeHM3ToULfjGGOM37E9kXysXr2ajRs3MmXKFNsLMcaYXNieSB6y9kLCwsIYNmyY23GMMcYv+bSIiEg3EdkjIoki8nwuyyuIyAJn+WYRicixvL6InBORP/syZ27i4uLYsGEDzz//PBUqVCjpjzfGmFLBZ0VERIKAyUB3oBkwUESa5eg2HPhRVRsBbwN/y7H8LeBTX2XMz/jx46lXrx7Dhw934+ONMaZU8OWeSDsgUVWTVDUNiAZ65+jTG5jjPF8IdBYRARCRPsABYKcPM+Zqx44drF27ljFjxhASElLSH2+MMaWGL4tcyInKAAAQEUlEQVRIKHDY63Wy05ZrH1VNB34CaohIZeA5YJwP8+VpypQpVKhQwcZCjDHmCvx1YP1l4G1VPZdfJxEZISLxIhJ/4sSJYvng1NRU5s6dy4ABA6hZs2axvKcxxgQqX57iewQI93od5rTl1idZRK4BqgKngPZAPxH5O1ANyBSRC6r6jvfKqjoNmAYQFRWlxRF63rx5nDt3jlGjRhXH2xljTEDzZRHZAjQWkUg8xWIA8HCOPjHAEGAj0A+IU1UFOmR1EJGXgXM5C4gvqCqTJ0+mTZs2tG3b1tcfZ4wxpZ7PioiqpovIaGAlEATMUtWdIvIKEK+qMcBMYK6IJAKn8RQa16xfv55du3Yxa9YsnPF9Y4wx+RDPH/6lX1RUlMbHxxfpPfr378+qVatITk6mUqVKxZTMGGP8l4gkqGpUYdf314H1Enfs2DEWL17MsGHDrIAYY0wBWRFxTJ8+nfT0dJ588km3oxhjTKlhRQS4dOkS7733Hl27dqVRo0ZuxzHGmFLDZvEFYmJiOHr0KFOnTnU7ijHGlCq2J4LnCvUGDRrQo0cPt6MYY0ypUuaLyO7du4mLi+PJJ58kKCjI7TjGGFOqlPki8u6771K+fHmbrdcYYwqhTBeRc+fOMWfOHB566CFq1arldhxjjCl1ynQRmT9/PqmpqTZPljHGFFKZLiLvv/8+LVq04LbbbnM7ijHGlEpltogcPHiQTZs2MWjQIJsnyxhjCqnMFpGPPvoIgIceesjlJMYYU3qV2SISHR1N+/btiYyMdDuKMcaUWmWyiOzdu5evv/6a/v37ux3FGGNKtTJZRBYsWICI2KEsY4wpojJZRKKjo+nQoQOhoaFuRzHGmFKtzBWRHTt2sGvXLjuUZYwxxaDMFZHo6GjKlStHv3793I5ijDGlXpkqIqpKdHQ0nTp1onbt2m7HMcaYUq9MFZGtW7eyf/9+BgwY4HYUY4wJCGWqiERHRxMcHMyDDz7odhRjjAkIPi0iItJNRPaISKKIPJ/L8goissBZvllEIpz2+0QkQUS2Oz87FTVLZmYmCxYsoEuXLlSvXr2ob2eMMQYfFhERCQImA92BZsBAEWmWo9tw4EdVbQS8DfzNaT8JPKCqLYAhwNyi5tm0aROHDx+2Q1nGGFOMfLkn0g5IVNUkVU0DooHeOfr0BuY4zxcCnUVEVPVrVT3qtO8EKopIhaKEiY6OJiQkhF69ehXlbYwxxnjxZREJBQ57vU522nLto6rpwE9AjRx9+gJbVfVizg8QkREiEi8i8SdOnMgzSEZGBh9//DE9evTguuuuu/pvYowxJld+PbAuIs3xHOIamdtyVZ2mqlGqGpXfnQnXr1/PDz/8YIeyjDGmmPmyiBwBwr1ehzltufYRkWuAqsAp53UY8AkwWFX3FyVIdHQ01157Lffff39R3sYYY0wOviwiW4DGIhIpIuWBAUBMjj4xeAbOAfoBcaqqIlINWA48r6obihLi0qVLLFq0iF69elGpUqWivJUxxpgcfFZEnDGO0cBKYDfwkaruFJFXRCRrdHsmUENEEoE/AlmnAY8GGgF/FZFtzqNQl5ivWbOGU6dO2aEsY4zxAVFVtzMUi6ioKI2Pj7+s/amnnmLu3LmcOnWKChWKdIKXMcYEHBFJUNWowq7v1wPrxSE2NpZOnTpZATHGGB8I6CKSmJhIUlISXbt2dTuKMcYEpIAuIitXrgSgS5cuLicxxpjAFNBFJDY2lsjISBo1auR2FGOMCUgBW0TS0tKIi4uja9euiIjbcYwxJiAFbBHZuHEj586ds/EQY4zxoYAtIrGxsQQFBXHPPfe4HcUYYwJWwBaRlStXcvvtt1O1alW3oxhjTMAKyCJy4sQJtm7daoeyjDHGxwKyiKxevRpVtVN7jTHGxwKyiKxcuZLq1avTpk0bt6MYY0xAC7gioqrExsZy3333ERQU5HYcY4wJaAFXRLZv386xY8dsPMQYY0pAwBWR2NhYAO677z6XkxhjTOALuCKycuVKmjdvTlhYmNtRjDEm4AVUEfn555/54osv7FCWMcaUkIAqIuvXr+fixYt2aq8xxpSQgCoiK1euJCQkhI4dO7odxRhjyoSAKyIdO3akYsWKbkcxxpgyIWCKSFpaGrt377bxEGOMKUE+LSIi0k1E9ohIoog8n8vyCiKywFm+WUQivJaNddr3iMgVK0NqaipgdzE0xpiS5LMiIiJBwGSgO9AMGCgizXJ0Gw78qKqNgLeBvznrNgMGAM2BbsAU5/3ylJqaSmhoKM2bNy/eL2KMMSZPvtwTaQckqmqSqqYB0UDvHH16A3Oc5wuBzuK5DWFvIFpVL6rqASDReb88paam0qVLF7uLoTHGlCBfFpFQ4LDX62SnLdc+qpoO/ATUKOC6v5KRkWGHsowxpoRd43aAohCREcAIgJCQEO69916XExljTNniyz2RI0C41+swpy3XPiJyDVAVOFXAdVHVaaoapapRzZs3p2bNmsUY3xhjzJX4sohsARqLSKSIlMczUB6To08MMMR53g+IU1V12gc4Z29FAo2Br3yY1RhjTCH47HCWqqaLyGhgJRAEzFLVnSLyChCvqjHATGCuiCQCp/EUGpx+HwG7gHTg96qa4ausxhhjCkc8f/iXflFRURofH+92DGOMKVVEJEFVowq7fsBcsW6MMabkWRExxhhTaFZEjDHGFJoVEWOMMYVmRcQYY0yhBczZWSJyFtjjdo4CqAmcdDtEAVjO4mU5i1dpyFkaMgLcpKpVCrtyqZ72JIc9RTlNraSISLzlLD6Ws3hZzuJTGjKCJ2dR1rfDWcYYYwrNiogxxphCC6QiMs3tAAVkOYuX5SxelrP4lIaMUMScATOwbowxpuQF0p6IMcaYEhYQRUREuonIHhFJFJHn3c7jTUQOish2EdmWdRaEiFQXkVUiss/5eb0LuWaJSIqI7PBqyzWXeExytu+3ItLa5Zwvi8gRZ5tuE5EeXsvGOjn3iEjXEsoYLiKfi8guEdkpImOcdr/anvnk9LftGSIiX4nIN07OcU57pIhsdvIscG4xgXPLiAVO+2YRiXA55/sicsBre97itLv5exQkIl+LyDLndfFtS1Ut1Q8808zvBxoC5YFvgGZu5/LKdxComaPt78DzzvPngb+5kKsj0BrYcaVcQA/gU0CA24DNLud8GfhzLn2bOf/9KwCRzv8XQSWQsS7Q2nleBdjrZPGr7ZlPTn/bngJUdp4HA5ud7fQRMMBpnwo85TwfBUx1ng8AFpTQ9swr5/tAv1z6u/l79EfgQ2CZ87rYtmUg7Im0AxJVNUlV04BooLfLma6kNzDHeT4H6FPSAVR1PZ57uHjLK1dv4AP12ARUE5G6LubMS28gWlUvquoBIBHP/x8+parHVHWr8/wssBsIxc+2Zz458+LW9lRVPee8DHYeCnQCFjrtObdn1nZeCHQWEXExZ15c+e8uImHA/cAM57VQjNsyEIpIKHDY63Uy+f9ilDQFYkUkQTz3hAeoo6rHnOc/AHXciXaZvHL54zYe7RwSmOV1OND1nM7u/614/ir12+2ZIyf42fZ0Dr9sA1KAVXj2gs6oanouWbJzOst/Amq4kVNVs7bnBGd7vi0iFXLmdJTU9pwIPAtkOq9rUIzbMhCKiL+7U1VbA92B34tIR++F6tlv9LtT5Pw1l+Nd4EbgFuAY8Ka7cTxEpDKwCHhGVVO9l/nT9swlp99tT1XNUNVbgDA8ez//4XKkXOXMKSI3A2Px5G0LVAeecyufiPQEUlQ1wVefEQhF5AgQ7vU6zGnzC6p6xPmZAnyC5xfieNZurPMzxb2Ev5JXLr/axqp63PnlzQSm88shFtdyikgwnn+Y56vqYqfZ77Znbjn9cXtmUdUzwOfA7XgO/2RN1eSdJTuns7wqcMqlnN2cw4aqqheB2bi7Pe8AeonIQTyH+jsB/6QYt2UgFJEtQGPnbIPyeAaDYlzOBICIXCsiVbKeA12AHXjyDXG6DQH+5U7Cy+SVKwYY7Jxdchvwk9dhmhKX4zjyg3i2KXhyDnDOMIkEGgNflUAeAWYCu1X1La9FfrU988rph9uzlohUc55XBO7DM37zOdDP6ZZze2Zt535AnLPn50bO77z+cBA8Yw3e27NE/7ur6lhVDVPVCDz/Nsap6iCKc1v6+qyAknjgOethL57jpn9xO49XroZ4zm75BtiZlQ3PMcY1wD5gNVDdhWz/i+fQxSU8x0SH55ULz9kkk53tux2IcjnnXCfHt87/9HW9+v/FybkH6F5CGe/Ec6jqW2Cb8+jhb9szn5z+tj1bAl87eXYAf3XaG+IpYonAx0AFpz3EeZ3oLG/ocs44Z3vuAObxyxlcrv0eOZ9/N7+cnVVs29KuWDfGGFNogXA4yxhjjEusiBhjjCk0KyLGGGMKzYqIMcaYQrMiYowxptCsiBiD5zx/EVknIkHF+J6fikiYiKwVkXzvtS0iTURkhXhm/N0qIh+JSB1n2Z3imS32O+cxwms97xl494nIYhFp5rU8WkQaF9d3MiYnKyLGeAwDFqtqRnG8mXPxWQ1VTS5A3xBgOfCuqjZWzzQ5U4BaInIDntlXn1TV/8BzrcdIEbnf6y3eVtVbVLUxsACIE5FazrJ38cybZIxPWBExAU1E2joT4YU4MwjsdOY3ymkQzlW7InK3s1fyLxFJEpE3RGSQszewXURudPrdKCKbnLbxInLO6/3uBtYWMObDwEZVXZrVoKprVXUH8Hvgff1l9t2TeIpCrvfNUdUFQKzzngBfAPd6TXFhTLGyImICmqpuwXMV9ng89/eY5/zjnM2ZLqehqh70am4FPAk0BR4FmqhqOzzTaT/t9Pkn8E9VbYHnanpv3YHPChjzZiCvCfKa57Is3mnPy1acCQvVMx9WIp7vY0yxsyJiyoJX8MxrFIWnkORUEziTo22LeibSu4hnmopYp307EOE8vx3PFBHgOeTk7Q7gy6LFLrSc939IAeq5EcQEPisipiyoAVTGcze/kFyW/zuX9otezzO9XmcC+R4aEpGGwGH13CStIHYCbfJYtiuXZW2cdfJyK54JC7OE4PmOxhQ7KyKmLHgP+C9gPvC3nAtV9UcgyBngvhqbgL7O8wFe7VdzKAs8ezG/8R4sF5GOztjNZGCo/HKf7hp4vkNue1SISF88s0X/r1dzE36ZSdaYYmVFxAQ0ERkMXFLVD4E3gLYi0imXrrF4zny6Gs8AfxSRb4FGeO4CB9CNy4vIchFJdh4fey9Q1X8DPYGnndN0d+G51/UJ9UwV/ggwXUS+A/4PmOU9CA/8IesUX6dvJ1U94Xz/OsC/VfWHq/xuxhSIzeJrDCAirYE/qOqjV7FOJTz/QKuIDAAGAg8BG1Q13+tCSoqI/AFIVdWZbmcxgclO+zMGUNWtIvK5iARdxbUibYB3nJsPnQGGOQPxflFAHGfw3C/EGJ+wPRFjjDGFZmMixhhjCs2KiDHGmEKzImKMMabQrIgYY4wpNCsixhhjCs2KiDHGmEL7/7NWbj+2+Mg5AAAAAElFTkSuQmCC\n", "text/plain": [ "