diff --git a/examples/howto/howto_debugging.ipynb b/examples/howto/howto_debugging.ipynb index 406bdb356..4ad78c8f5 100644 --- a/examples/howto/howto_debugging.ipynb +++ b/examples/howto/howto_debugging.ipynb @@ -4,332 +4,553 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "(howto_debugging)=\n", "# How to debug a model\n", "\n", + ":::{post} August 2, 2022\n", + ":tags: debugging, aesara\n", + ":category: beginner\n", + ":author: Thomas Wiecki, Igor Kuvychko\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Introduction\n", "There are various levels on which to debug a model. One of the simplest is to just print out the values that different variables are taking on.\n", "\n", - "Because `PyMC3` uses `Theano` expressions to build the model, and not functions, there is no way to place a `print` statement into a likelihood function. Instead, you can use the `Theano` `Print` operatator. For more information, see: theano Print operator for this before: http://deeplearning.net/software/theano/tutorial/debug_faq.html#how-do-i-print-an-intermediate-value-in-a-function.\n", - "\n", - "Let's build a simple model with just two parameters:" + "Because `PyMC` uses `Aesara` expressions to build the model, and not functions, there is no way to place a `print` statement into a likelihood function. Instead, you can use the `aesara.printing.Print` class to print intermediate values." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, + "outputs": [], + "source": [ + "import arviz as az\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pymc as pm" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "%config InlineBackend.figure_format = \"retina\"\n", + "\n", + "RANDOM_SEED = 8927" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### How to print intermediate values of `Aesara` functions\n", + "Since `Aesara` functions are compiled to C, you have to use `aesara.printing.Print` class to print intermediate values (imported below as `Print`). Python `print` function will not work. Below is a simple example of using `Print`. For more information, see {ref}`Debugging Aesara `." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import aesara.tensor as at\n", + "\n", + "from aesara import function\n", + "from aesara.printing import Print" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, "outputs": [ { - "name": "stderr", + "data": { + "text/plain": [ + "array([ inf, 0.5 , 0.25])" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = at.dvector(\"x\")\n", + "y = at.dvector(\"y\")\n", + "func = function([x, y], 1 / (x - y))\n", + "func([1, 2, 3], [1, 0, -1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To see what causes the `inf` value in the output, we can print intermediate values of $(x-y)$ using `Print`. `Print` class simply passes along its caller but prints out its value along a user-define message:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", "output_type": "stream", "text": [ - "Multiprocess sampling (4 chains in 4 jobs)\n", - "CompoundStep\n", - ">Metropolis: [sd]\n", - ">Metropolis: [mu]\n" + "x - y = __str__ = [0. 2. 4.]\n" ] }, { "data": { - "text/html": [ - "\n", - "
\n", - " \n", - " \n", - " 100.00% [24000/24000 00:14<00:00 Sampling 4 chains, 0 divergences]\n", - "
\n", - " " - ], "text/plain": [ - "" + "array([ inf, 0.5 , 0.25])" ] }, + "execution_count": 5, "metadata": {}, - "output_type": "display_data" - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 15 seconds.\n", - "/env/miniconda3/lib/python3.7/site-packages/xarray/core/nputils.py:215: RuntimeWarning: All-NaN slice encountered\n", - " result = getattr(npmodule, name)(values, axis=axis, **kwargs)\n" - ] + "output_type": "execute_result" } ], "source": [ - "%matplotlib inline\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import pandas as pd\n", - "import pymc3 as pm\n", - "import seaborn as sns\n", - "import theano.tensor as tt\n", - "\n", - "x = np.random.randn(100)\n", - "\n", - "with pm.Model() as model:\n", - " mu = pm.Normal(\"mu\", mu=0, sigma=1)\n", - " sd = pm.Normal(\"sd\", mu=0, sigma=1)\n", + "z_with_print = Print(\"x - y = \")(x - y)\n", + "func_with_print = function([x, y], 1 / z_with_print)\n", + "func_with_print([1, 2, 3], [1, 0, -1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`Print` reveals the root cause: $(x-y)$ takes a zero value when $x=1, y=1$, causing the `inf` output." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### How to capture `Print` output for further analysis\n", "\n", - " obs = pm.Normal(\"obs\", mu=mu, sigma=sd, observed=x)\n", - " step = pm.Metropolis()\n", - " trace = pm.sample(5000, step)" + "When we expect many rows of output from `Print`, it can be desirable to redirect the output to a string buffer and access the values later on (thanks to **Lindley Lentati** for inspiring this example). Here is a toy example using Python `print` function:" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([0., 0., 0., ..., 0., 0., 0.])" + "['Test values: 0',\n", + " 'Test values: 1',\n", + " 'Test values: 2',\n", + " 'Test values: 3',\n", + " 'Test values: 4',\n", + " '']" ] }, - "execution_count": 2, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "trace[\"mu\"]" + "import sys\n", + "\n", + "from io import StringIO\n", + "\n", + "old_stdout = sys.stdout\n", + "mystdout = sys.stdout = StringIO()\n", + "\n", + "for i in range(5):\n", + " print(f\"Test values: {i}\")\n", + "\n", + "output = mystdout.getvalue().split(\"\\n\")\n", + "sys.stdout = old_stdout # setting sys.stdout back\n", + "output" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Hm, looks like something has gone wrong, but what? Let's look at the values getting proposed using the `Print` operator:" + "### Troubleshooting a toy PyMC model" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "rng = np.random.default_rng(RANDOM_SEED)\n", + "x = rng.normal(size=100)\n", + "\n", + "with pm.Model() as model:\n", + " # priors\n", + " mu = pm.Normal(\"mu\", mu=0, sigma=1)\n", + " sd = pm.Normal(\"sd\", mu=0, sigma=1)\n", + "\n", + " # setting out printing for mu and sd\n", + " mu_print = Print(\"mu\")(mu)\n", + " sd_print = Print(\"sd\")(sd)\n", + "\n", + " # likelihood\n", + " obs = pm.Normal(\"obs\", mu=mu_print, sigma=sd_print, observed=x)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "mu __str__ = 0.0\n", - "sd __str__ = 0.0\n" - ] - }, + "data": { + "image/svg+xml": [ + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "%3\r\n", + "\r\n", + "cluster100\r\n", + "\r\n", + "100\r\n", + "\r\n", + "\r\n", + "sd\r\n", + "\r\n", + "sd\r\n", + "~\r\n", + "Normal\r\n", + "\r\n", + "\r\n", + "obs\r\n", + "\r\n", + "obs\r\n", + "~\r\n", + "Normal\r\n", + "\r\n", + "\r\n", + "sd->obs\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "mu\r\n", + "\r\n", + "mu\r\n", + "~\r\n", + "Normal\r\n", + "\r\n", + "\r\n", + "mu->obs\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n", + "\r\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "pm.model_to_graphviz(model)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "Only 3 samples in chain.\n", - "Sequential sampling (1 chains in 1 job)\n", - "CompoundStep\n", - ">Metropolis: [sd]\n", - ">Metropolis: [mu]\n", - "Sampling 1 chain for 0 tune and 3 draw iterations (0 + 3 draws total) took 0 seconds.\n", - "/dependencies/arviz/arviz/data/base.py:146: UserWarning: More chains (3) than draws (2). Passed array should have shape (chains, draws, *shape)\n", - " UserWarning,\n" + "Only 5 samples in chain.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "sd __str__ = 1.753782660690967\n", - "mu __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "mu __str__ = -1.179953715933347\n", - "mu __str__ = 0.0\n", - "sd __str__ = -0.7921743241812005\n", - "mu __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "mu __str__ = -1.6164899771368246\n", - "mu __str__ = 0.0\n", - "sd __str__ = 1.7649365130019956\n", - "mu __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "mu __str__ = 0.5633904409221188\n", - "mu __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "mu __str__ = 0.0\n", - "sd __str__ = 0.0\n", - "mu __str__ = 0.0\n", "sd __str__ = 0.0\n", "mu __str__ = 0.0\n" ] }, { - "name": "stderr", - "output_type": "stream", - "text": [ - "/dependencies/pymc3/pymc3/sampling.py:603: UserWarning: The number of samples is too small to check convergence reliably.\n", - " warnings.warn(\"The number of samples is too small to check convergence reliably.\")\n" + "ename": "SamplingError", + "evalue": "Initial evaluation of model at starting point failed!\nStarting values:\n{'mu': array(0.), 'sd': array(0.)}\n\nInitial evaluation results:\n{'mu': -0.92, 'sd': -0.92, 'obs': -inf}", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mSamplingError\u001b[0m Traceback (most recent call last)", + "Input \u001b[1;32mIn [9]\u001b[0m, in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m model:\n\u001b[0;32m 2\u001b[0m step \u001b[38;5;241m=\u001b[39m pm\u001b[38;5;241m.\u001b[39mMetropolis()\n\u001b[1;32m----> 3\u001b[0m trace \u001b[38;5;241m=\u001b[39m \u001b[43mpm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msample\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m5\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstep\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtune\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mchains\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mprogressbar\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrandom_seed\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mRANDOM_SEED\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32mc:\\users\\igork\\pycharmprojects\\pymc\\pymc\\sampling.py:558\u001b[0m, in \u001b[0;36msample\u001b[1;34m(draws, step, init, n_init, initvals, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, **kwargs)\u001b[0m\n\u001b[0;32m 556\u001b[0m \u001b[38;5;66;03m# One final check that shapes and logps at the starting points are okay.\u001b[39;00m\n\u001b[0;32m 557\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m ip \u001b[38;5;129;01min\u001b[39;00m initial_points:\n\u001b[1;32m--> 558\u001b[0m \u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcheck_start_vals\u001b[49m\u001b[43m(\u001b[49m\u001b[43mip\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m 559\u001b[0m _check_start_shape(model, ip)\n\u001b[0;32m 561\u001b[0m sample_args \u001b[38;5;241m=\u001b[39m {\n\u001b[0;32m 562\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdraws\u001b[39m\u001b[38;5;124m\"\u001b[39m: draws,\n\u001b[0;32m 563\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mstep\u001b[39m\u001b[38;5;124m\"\u001b[39m: step,\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 573\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdiscard_tuned_samples\u001b[39m\u001b[38;5;124m\"\u001b[39m: discard_tuned_samples,\n\u001b[0;32m 574\u001b[0m }\n", + "File \u001b[1;32mc:\\users\\igork\\pycharmprojects\\pymc\\pymc\\model.py:1794\u001b[0m, in \u001b[0;36mModel.check_start_vals\u001b[1;34m(self, start)\u001b[0m\n\u001b[0;32m 1791\u001b[0m initial_eval \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpoint_logps(point\u001b[38;5;241m=\u001b[39melem)\n\u001b[0;32m 1793\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mall\u001b[39m(np\u001b[38;5;241m.\u001b[39misfinite(v) \u001b[38;5;28;01mfor\u001b[39;00m v \u001b[38;5;129;01min\u001b[39;00m initial_eval\u001b[38;5;241m.\u001b[39mvalues()):\n\u001b[1;32m-> 1794\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m SamplingError(\n\u001b[0;32m 1795\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mInitial evaluation of model at starting point failed!\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 1796\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mStarting values:\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;132;01m{\u001b[39;00melem\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 1797\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mInitial evaluation results:\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;132;01m{\u001b[39;00minitial_eval\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 1798\u001b[0m )\n", + "\u001b[1;31mSamplingError\u001b[0m: Initial evaluation of model at starting point failed!\nStarting values:\n{'mu': array(0.), 'sd': array(0.)}\n\nInitial evaluation results:\n{'mu': -0.92, 'sd': -0.92, 'obs': -inf}" ] } ], "source": [ - "with pm.Model() as model:\n", - " mu = pm.Normal(\"mu\", mu=0, sigma=1)\n", - " sd = pm.Normal(\"sd\", mu=0, sigma=1)\n", - "\n", - " mu_print = tt.printing.Print(\"mu\")(mu)\n", - " sd_print = tt.printing.Print(\"sd\")(sd)\n", - "\n", - " obs = pm.Normal(\"obs\", mu=mu_print, sigma=sd_print, observed=x)\n", + "with model:\n", " step = pm.Metropolis()\n", - " trace = pm.sample(\n", - " 3, step, tune=0, chains=1, progressbar=False\n", - " ) # Make sure not to draw too many samples" + " trace = pm.sample(5, step, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In the code above, we set the `tune=0, chains=1, progressbar=False` in the `pm.sample`, this is done so that the output is cleaner.\n", - "\n", - "Looks like `sd` is always `0` which will cause the logp to go to `-inf`. Of course, we should not have used a prior that has negative mass for `sd` but instead something like a `HalfNormal`.\n", - "\n", - "We can also redirect the output to a string buffer and access the proposed values later on (thanks to [Lindley Lentati](https://github.com/LindleyLentati) for providing this example):" + "Exception handling of PyMC v4 has improved, so now SamplingError exception prints out the intermediate values of `mu` and `sd` which led to likelihood of `-inf`. However, this technique of printing intermediate values with `aeasara.printing.Print` can be valuable in more complicated cases." ] }, { - "cell_type": "code", - "execution_count": 4, + "cell_type": "markdown", "metadata": {}, + "source": [ + "### Bringing it all together" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "scrolled": true + }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "Only 5 samples in chain.\n", + "Only 10 samples in chain.\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", "Sequential sampling (1 chains in 1 job)\n", - "CompoundStep\n", - ">Metropolis: [sd]\n", - ">Metropolis: [mu]\n", - "Sampling 1 chain for 0 tune and 5 draw iterations (0 + 5 draws total) took 0 seconds.\n", - "/dependencies/arviz/arviz/data/base.py:146: UserWarning: More chains (5) than draws (2). Passed array should have shape (chains, draws, *shape)\n", - " UserWarning,\n" + "NUTS: [mu, a, b]\n", + "Sampling 1 chain for 0 tune and 10 draw iterations (0 + 10 draws total) took 1 seconds.\n", + "The chain contains only diverging samples. The model is probably misspecified.\n", + "The acceptance probability does not match the target. It is 0, but should be close to 0.8. Try to increase the number of tuning steps.\n", + "C:\\Users\\igork\\AppData\\Local\\Temp\\ipykernel_14804\\1992602661.py:15: UserWarning: The number of samples is too small to check convergence reliably.\n", + " trace = pm.sample(draws=10, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED)\n" ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAEGCAYAAAC0IuZwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9e3xdZZXw/12535ukTShtAuVSLi0ilwriyEVQB5lRQJSLzoi+Og4qo+NlHPw5L4Oo7wjeZ3SG11FeL8MICOKAIOooiqOgtFBKQ6G0XEyaQtLm5J6cXM76/bGfJ9k9PUlOkrNP9mnX9/PJJ/s8+9nPXnv3ss66PGuJqmIYhmEYxtJStNQCGIZhGIZhCtkwDMMwYoEpZMMwDMOIAaaQDcMwDCMGmEI2DMMwjBhQstQCGIXLihUrdM2aNUsthmEYRsGwadOmParalOmcKWRjwaxZs4aNGzcutRiGYRgFg4i8MNM5c1kbhmEYRgwwhWwYhmEYMcAUsmEYhmHEAFPIhmEYhhEDTCEbhmEYRgzIm0IWkfNF5GkR2SEi12Q4/04R6RaRze7nPW78JBF5SETaRGSLiFwWukZE5LMisl1EtonIB93434XW2SoikyLS6M7dLCJdIrI17f6fdutvFpGficgqN35haHyjiLw6dM2NTq5tIvLPIiJu/DJ3TZuI3BCa/+WQXNtFpNeNHy4ij7rxNhG5KnTNFSLyhFvvfhFZ4cbf6uamRGRDaP7rRGSTu2aTiJwbOlcmIt9w935KRC5x44eJyAMi8pi7zwUL+1M2DMMwFoyqRv4DFAM7gSOBMuBxYF3anHcCX8tw7THAWne8CtgN1LvP7wK+CxS5z80Zrn8j8MvQ57OAU4CtafPqQscfBG5yxzWAuOMTgafc8auA37pnKwYeAs4BlgN/BJrcvO8A52WQ62+Am91xGVAeut/z7llLgC5ghTt3I3CdOz4eOBb4FbAhtO7JwCp3fAKwK3TuU8Bn3HFRaN1vAO9zx+uA57P5cz311FPVMAzDyB5go87wf2q+LOTTgB2q+qyqjgG3Ahdmc6GqblfVZ9xxJ4GC8puq3wdcr6opd74rwxJXAN8Prfcg0JPhPv2hj9WAuvFB9xL3GXe/K3DKFCgFXiL40vGMqna7ef8NXDKbXKo6pqpJN17OtOdC3E+1s77rgE53zTZVfTrDczzm3hNAG1ApIuXu8/8C/snNS6nqntCz1LnjZf4ehmEUBr95ppvn9gwttRjGIsmXQl4NtIc+d7ixdC5xLtM7RKQ1/aSInEagAHe6oaOAy5wr+ScisjZtfhVwPnBnNkI693c78Hbg2tD4xSLyFHAvgVJDVR8CHiCw2HcDP1XVbcAO4FgRWSMiJcBFQGvafQ4HjgB+GRprFZEtBO/pBlXtVNVxgi8dTxAoyXXAt7J5FsclwKOqmhSRejf2aece/4GIHOLGrgP+QkQ6gPsIrPeZ3tF73fve2N3dPdM0wzDyyIdu3cxNv9o590Qj1sQpqeseYI2qngj8nMDVO4WIHAp8D3iXt4gJrMlRVd0A/Dtwc9qabwR+q6r7WcSZUNVPqmorcAtwdWj8LlU9jkC5ftrJczSB27iF4MvFuSJypqomCJTobcBvCNzPk2m3uhy4Q1WnxlW13T370cCVInKIiJS6tU4mcGFvAT6RzbOIyHrgBuCv3VCJk/V3qnoKgYv9C+7cFcC3VbUFuAD4nohk/Luhqt9Q1Q2quqGpKWP1N8Mw8shkSkkMj7F3KDn3ZCPW5Esh72JfK7HFjU2hqntDbttvAqf6cyJSR2CdflJVHw5d1gH80B3fRRDjDXM5IXf1PLiFDG5m5+4+0iVWXQw87Fzag8BPgDPcvHtU9XRVPQN4GtierVzO3bwVOBM4yY3tdG7z2wli17MiIi0E7+Mdquq/Nu8Fhpl+Xz8giKUDvNut7S3/CmDFXPcxDGPp6RsZRxUSw+NLLYqxSPKlkB8B1orIESJSRqCQ7g5PcBaw503ANjdeRqBcvquqd6St+yPgNe74bEKKT0SWubH/ykbANHf3hcBTbvzoUPb0KQRW+V6CxK2zRaTEWbJnh2Rudr8bgPcTfMHw9zkOaCCwUP1Yi4hUhq55NYEi3wWsExFvir7O32OW56gn+PJyjar+1o87hX4PQeIZwHnAk+74j+4zInI8gUI2f7RhFAA9Q2MAJNxvo3DJS3MJVZ0QkauBnxJkJN+sqm0icj1BxtndwAdF5E3ABEHS1Tvd5ZcSZEYvFxE/9k5V3Qx8DrhFRD4MDALvCd32YuBnqrpPpoOIfJ9AKa1wMdN/VNVvAZ8TkWOBFPAC4LceXQK8Q0TGgRHgMlVVEbkDOJcgvqvA/ap6j7vmqyLycnd8vaqGLeTLgVtDiWIQuL6/KCJKkMT1BVV9wsn7KeBBd/8X/HsRkYuBfyFIcLtXRDar6p8SuNqPBq4VER8Hf71LePt7Anf0VwgU7rvc+Y8C/+7eo7r3G5bPMIyYkhgOFHHPsCnkQkfs/11joWzYsEGt25NhLC0/a3uR935vEyKw47MXUFwkSy2SMQsissnlPe1HnJK6DMMwjHniLWTVIJ5sFC6mkA3DMAqYnqHx0LG5rQsZU8iGYRgFTCIUO05YHLmgMYVsGIZRwISzqy3TurAxhWwYhlHAJIbHqK8qnTo2ChdTyIZhGAVMz9AYRzXVuGNL6ipkTCEbhmEUML3D4xy6rILykiJ6zUIuaEwhG4ZhFDA9w2M0VpfRWF1mWdYFjilkwzCMAmUypfSNjNNQVUZDVZnFkAscU8iGYRgFim8s0VBVSkN1qTWYKHBMIRuGYRQo3kXdUO0sZHNZFzSmkA3DMAoU76KeiiGby7qgMYVsGIZRoExZyC6G3DcyzmTKGgYVKqaQDcMwChS/zSlwWZdag4kCxxSyYRhGgeILgTRWldFQXebGzG1dqJhCNgzDKFASw2NUlBZRWVZMo1PItvWpcDGFbBiGUaAkhsZoqAoUsf9tmdaFiylkwzCMAiUxHFLIZiEXPKaQDcMwCpSeobEpV3VjlY8hW1JXoWIK2TAMo0BJDI9PtV6sLCumvKTILOQCxhSyYRhGgZIYnraQISgQYjHkwsUUsmEYRgEyMZmaaizhsQYThY0pZMMwjALEN5ZIt5BtH3LhYgrZMAyjAPGdnXwM2R/3WsengsUUsmEYRgESbizhsQYThY0pZMMwjAIk3FjC4xtMTEymlkosYxGYQjYMwyhAEqFeyB5rMFHYmEI2DMMoQHwMuTFsIU9V6zKFXIiYQjYMwyhAwo0lPNZgorAxhWwYhlGA9AyN7WMdw3Q82bY+FSamkA3DMAqQ3uEx6tMVsrOQe81CLkjyopBF5DgReUhEkiLysVnmHSEivxeRHSJym4iUufEvi8hm97NdRHrTrqsTkQ4R+Zr7XBuav1lE9ojIV9KuuUREVEQ2uM+lIvIdEXlCRLaJyCdCc28WkS4R2Zq2xudF5CkR2SIid4lIvRtfIyIjofvfFLrmMje/TURuCI2Xu2fe4d7BmtC5E937a3PyVbjxz4pIu4gMpsl1uIj8wt3nVyLSMtv7cmNXuLW3iMj9IrJipj8nwzCWnnBjCY81mChs8mUh9wAfBL4wx7wbgC+r6tFAAng3gKp+WFVPUtWTgH8Bfph23aeBB/0HVR3w8901L4SvEZFa4EPA70NrvBUoV9WXAacCfx1Sit8Gzs8g78+BE1T1RGA78InQuZ0hGa5y910OfB44T1XXAytF5Dw3/91Awj37l927QERKgP8ArnLXnAP4f233AKdlkOsLwHedXNcD/zTb+3L3+CrwGnfNFuDqDOsahhETEsPj+2RYQ9BgoqLUGkwUKnlRyKrapaqPMK1I9kNEBDgXuMMNfQe4KMPUK4Dvh647FTgE+NkM6x4DNAO/CQ1/mkDhjYbFBKqdcqoExoB+J/+DBF8q0p/rZ6o64T4+DLSkz0njSOAZVe12n/8buMQdX0jwzBC8g/PcO3k9sEVVH3f33Kuqk+74YVXdneE+64BfuuMH3NrAjO9L3E+1u2cd0DnHsxjGAcftG9t5x81/WGoxsqJnaIyGUJUuT0OVlc8sVOIUQ14O9IYUXAewOjxBRA4HjsApGxEpAr4IzOgGBy4HblNVddecArSq6r1p8+4AhoDdwB+BL6jqfkp4Fv4X8JPQ5yNE5DER+bWInOnGdgDHOpd2CcEXjlZ3bjXQDuDeQR/BOzkGUBH5qYg8KiIfz0KWx4E3u+OLgVoRWT7T+1LVceB9wBMEingd8K1MC4vIe0Vko4hs7O7uzjTFMAqWh3fu5cHt3YxNxLuwxsRkiv7RfRtLeBqqyiyGXKDESSFnw+XAHd5CBN4P3KeqHXNc832YUuBfAj6aYd5pwCSwikDpf1REjsxGKBH5JDAB3OKGdgOHqerJwEeA/xSROlVNECi+2wgs9ufdPWejBHg18Hb3++KQm3smPgacLSKPAWcDu9x9Mr4vESl1cp1M8Pxb2Nf9PoWqfkNVN6jqhqampjnEMIzComsgCcCeweQSSzI7mRpLeKzBROFSEtXCIvIB4K/cxwtUdS4X6F6gXkRKnIXYQqBIwlwOfCD0+QzgTBF5P1ADlInIoKpe42R4OVCiqpvc/FrgBOBXgWeWlcDdIvIm4G3A/c5a7BKR3wIbgGfneM53An9OEBdWAFVNAkl3vElEdhJYuhtV9R6C2C8i8l6mFfIuAmu5w1nPy9w76QAeVNU97pr7gFOAX8wkk3vXb3bza4BLVLVXRDK+L+BOd91Od83twDWzPbdhHIh0DYy630lW1VcusTQz42PE6TFkP7ardyTfIhk5IDILWVW/HkpqmjMe6ZTZA8Bb3NCVwH/58yJyHNAAPBS65u2qepiqriGwCr/rlbFjn3izqvap6gpVXeOueRh4k6puJHBTn+vuVQ28EnhqNplF5Hzg426N4dB4k4gUu+MjgbU4xS4ize53A4HF+k132d3umXHv4JfunfwUeJmIVDlFfTbw5BxyrXDeAAgs3ZvneF+7gHUi4k3e1wHbZruHYRyIdDsL2f+OK74SV+YYcqkldRUo+dr2tFJEOgjct//gttzUuXP3icgqN/XvgY+IyA6C+Gk4jnk5cKu3QrPkUkIKeQ6+DtSISBvwCPD/VHWLk/H7BF8EjnWyv9td8zUCq/vnadubzgK2iMhmgtj0VaF49FdF5Engt8DnVHW7G/8WsNw9+0dwFqpzc3/JybQZeNTHv0XkRvdeq5xc17m1zgGeFpHtBAlcn53twd0Xpk8BD4rIFuAk4P9k+d4M44BgbCI1pei8pRxXMjWW8FiDicIlMpd1GFV9kRkykFX1gtDxs2TexoOqXjfHPb5NsD0pPDZrDFhVzwkdDxJsfco074oZxo+eYfxOnBt4HmuNznL//yDY+pQ+/nECCz19/A6ms9Uzkv6+VPUm4KaZ5hvGgU53KG7c1R9zC3lo/9aLnsbqsqkGE8tryvMtmrEICi2pyzAMIxLCburumCd1+Z7HmSzkeufGNrd14WEK2TAMA+jqD9zURRJ/C7l3eHy/xhKeRuv4VLCYQjYMw2DaKj66uSb+FnKGxhIeazBRuJhCNgzDILCKReC4lXV098c7qSsxtH9jCc9UT2RTyAWHKWTDMAyCvceNVWUcWl9B92CS+W3oyC+J4f0bS3i85Wwu68LDFLJhGAZBUldTbTnNtRWMTyq9MVZomRpLeKzBROFiCtkwDAPoHhilua6C5tpgq1BXjIuDBDHk/YuCeBqtwURBYgrZMAyDQAE31ZTTNKWQ4xlHnphM0TcyPmMMGaC+qsxiyAWIKWTDMA56Uillz2CS5rryKQs5ruUz+0YCV/pMMWR/zlzWhYcpZMMwDnp6R8YZn1Saa8tprqsA4uuynq2xhKehusySugoQU8iGYRz0ePd0U2051WXFVJYWx7Y4SM/QzI0lPA1VpRZDLkBMIRuGcdDj3dPNtRWICM115bEtDpKYpWymp6GqjP5RazBRaJhCNgzjoMdbwz6hq7m2fKqUZtyYrbGEJ9xgwigcTCEbhnHQ461hn9DVVBtfC3m2xhKeqWpdlthVUJhCNgzjoKerP0l1WTHV5UFH2ubaCrpjGkNODI3N2FjC4+PLPt5sFAamkA3DOOjpGhidcldDYCEPJCcYGZtcQqkykxgen7GxhKehyizkQsQUsmEYBz3dA0maayumPjfFeC9yYmhs1i1PEGrBaJnWBYUpZMMwDnq6B5I01U1byM0xrtbVMzw2a/wYQi0YzUIuKEwhG4Zx0OPLZnqaYlzPuneWxhIe32Aizg0yjP0xhWwYxkHN8NgEg8kJmvexkAP3dRxd1nM1lvBYg4nCwxSyYRgHNV7phi3kxuoyioskdi7rbBpLeKzBROFhCtkwjIOaqSpdddNJXcVFwvLqsthZyL1ZNJbwNFaXWQy5wDCFbBjGQU3XwL5FQTzNdeWxiyH3ZtFYwtNQXWYx5ALDFLJhGAc1vkRmU7pCrq2IXYMJX+hjrn3IwRxrMFFomEI2DOOgpnswSUmR7KfkmmriVz7TK9j6LJK66qvK6BuxBhOFhClkwzAOarr6k6yoKaeoSPYZb64rZ+9gksmULpFk++Nd1tnGkMEaTBQSppANwzio6RpI7ueuhsCFnVLYGyMrOZvGEh5rMFF4mEI2DOOgJiibub9Cbo5hcZBsGkt4rMFE4WEK2TCMg5qugeQ+RUE8TTEsDtIzNHdjCc9U+UxL7CoYTCEbhnHQMjGZYu/QvmUzPXGsZ907PHdjCY+PIfeay7pgyJtCFpHzReRpEdkhItdkOP9lEdnsfraLSG/o3GTo3N2h8d+ExjtF5Edu/BwR6Qudu3YuOSTgs+7e20Tkg6Hxf3bzt4jIKW78cBF51K3fJiJXhdb6rIi0i8jgDO/iEhFREdngPr89JOtmEUmJyEnu3BUi8oS79/0issKN3xaa/7yIbHbjrxORTe6aTSJybui+M631VvcMKS+TYRwM9AyNoQpNoaIgnjh2fOoZHssqoQuswUQhUpKPm4hIMfB14HVAB/CIiNytqk/6Oar64dD8vwFODi0xoqonpa+rqmeGrrkT+K/Q6d+o6p/PQ453Aq3AcaqaEpFmd9kbgLXu53Tg39zv3cAZqpoUkRpgq1urE7gH+BrwTIZ3UQt8CPh96DluAW5x518G/EhVN4tICfBVYJ2q7hGRG4GrgetU9bLQml8E+tzHPcAbVbVTRE4Afgqsnm0tYCvwZuD/pstrGAcyXRnKZnoqSoupqyiJXQy5paEqq7m+wYSVzywc5rSQReQYEfmFiGx1n08UkX+Y531OA3ao6rOqOgbcClw4y/wrgO9nu7iI1AHnAj9ahBzvA65X1RSAqna58QuB72rAw0C9iByqqmOq6v+llhN6l6r6sKrunkGGTwM3ADP5wa5wcgGI+6kWEQHqgM60ZxfgUtz7UtXH3JcCgDagUkTKZ1tLVbep6tMzyGMY80ZVuX/ri7HfAztdNnN/hQyBlRwnCzkxPJ5VYwlPY1UZCavWVTBk47L+d+ATwDiAqm4BLp/nfVYD7aHPHW5sP0TkcOAI4Jeh4QoR2SgiD4vIRRkuuwj4har2h8bOEJHHReQnIrI+CzmOAi5z9/mJiKyd6xoRaRWRLe78DSFFmBHn7m5V1XtnmXYZ08p1nOCLwhMEynMd8K20+WcCL6nqftY4cAnwqKoms1xrTkTkve4dbezu7p7v5cZBwqYXElz1H5t44Ol4/x3x8eFMWdbBeEVsLGTfWCLbGDIEW5/MQi4cslHIVar6h7SxiSiEcVwO3KGqk6Gxw1V1A/A24CsiclTaNekW9aPumpcD/8LcljMEVu6ou8+/AzfPdYGqtqvqicDRwJUicshMc0WkCPgS8NFZ5pwODKuq90aUEijRk4FVwBaCL0dhMnoT3JeQG4C/nsdac6Kq31DVDaq6oampab6XGwcJz+0ZAuCFvUNLLMns+NKYKzK4rMHXs45HUpdvLJHNHmRPQ5U1mCgkslHIe5wCVAAReQtB/HQ+7CKIz3pa3FgmLidNwajqLvf7WeBXhOLLLjHpNODe0Px+VR10x/cBpW7ebHJ0AD90x3cBJ2Yru7OMtxJYqzNRC5wA/EpEngdeCdydlkSV/uwnufV3qqoCtwOvCj17CUHs97bwjUSkxT3DO1R1ZzZrGUYuaU+MBL97hpdYktnpHkyyrLKUitLM+3qbagKXdfBPZmnxlq5ZyAcu2SjkDxAk+xwnIruAvyWwtObDI8BaETlCRMoIFM/d6ZNE5DigAXgoNNbgYqBe+f4J8GTosrcAP1bV0dA1K12cFBE5jeA5984hx4+A17jjs4Ht7vhu4B0u2/qVQJ+q7haRFhGp9DICrwZmjMOqap+qrlDVNaq6BngYeJOqbnRrFBHEgm8NXbYLWCci3hR9HbAtdP61wFOq2hF69nqCLyfXqOpv57GWYeSMDqeIvWKOK139mYuCeJrryhkdTzGQjNIpmB0+FpztPuRgbqnFkAuIObOsnVX6WhGpBopUdWC+N1HVCRG5miDjtxi4WVXbROR6YKOqeqV4OXCr7vt19Hjg/4pIikCxfi6cne2u+VzaLd8CvE9EJoAR4HK3ZkY53DWfA24RkQ8Dg8B73Ph9wAXADmAYeFdIri+KiBIkS31BVZ8AcBnMbwOqRKQD+KaqXjfHazoLaHfv27+3ThH5FPCgiIwDLxBkg4efPd1dfTWBC/1amd7u9frZ1hKRiwlc+03AvSKyWVX/dA55DWNG2hNOIcfcQu4aGM1YNtPjz3X1J6mryD6ZKgrm01jCE24wUVJsZSfizpwKOfSfuv8MgKpeP58bOdfxfWlj16Z9vi7Ddb8DXjbLuudkGPsawbajrORw473An2UYVwIvQfr4z5l2a6ef+zjw8ZlkziS3qv6KwI2dPu8m4KYZ1nhnhrHPAJ+ZYX7GtVT1LgIXt2HkhA5nGXckRlDVqf834kb3YJJTD2uY8XxzqFrX0c01+RIrI4l5NJbwTBUHGRmfMU5uxIdsvjINhX4mCfblrolQJsMwCpjkxCQv9o/SWF3GyPgke2Maw1RVuvozN5bwxKlaV2IejSU8DVatq6DIxmX9xfBnEfkCgcvXMAxjPzp7R1GFVx7ZyH1PvEh7z3AsrbP+0QmSE6kpKzgTcarWlRgao7K0OKvGEp7GqXrWFkcuBBYSVKgiyDQ2DMPYDx83PuPI5cHnmCZ2zVUUBGBZZSllJUWxUMg9Q+NTHZyypX6q45NZyIVANjHkJ3BbnggSoZqAecWPDcM4ePAJXWcc5RRyTBO7vBs6U9lMj4jQVFMei+IgiXk0lvA0Wk/kgiKbWtbhetATBFWhln4PgGEYsaS9Z4TSYuGIFTUsry6jIxFPhZyNhQzxKZ+ZmEdjCY+PN5tCLgxmVMgi0ugO07c51YkIqtoTnViGYRQqHYlhVtVXUlwktDRW0d4Tb5d10ywxZAgSu56PQcWx+TSW8FiDicJiNgt5E4GrOtN+BQWOjEQiwzAKmvbECK1OcbQ2VLJ1V98cVywNXQNJykqKqKuY3VHYXFfOH55fevujZ2hsXo0lPI1VZZbUVSDM+DdRVY/IpyCGYRwYdPQM8/r1QVn31sYqftr2IpMppbgoXnuRuweCKl1z7ZFuqqmgd3ic5MQk5SXZZzjnkonJFP2jE/OOIUOw9cm2PRUGWfVDdqUh1wJTvh1VfTAqoQzDKEyGkhPsDblWWxoqGZ9UXuofZVV95RJLty9dA6Ozls30+BjznsExVi/RM/jGEvONIftrrMFEYZBNP+T3AA8S7D3+lPt9XbRiGYZRiPgKXS0NgeLyrus4ZlrPVRTE47Owu/qXrjhIYqps5vwVcn2VNZgoFLLZh/wh4BXAC6r6GoJOS72RSmUYRkHiFW9rY9U+v+O4F7l7MDlrURCPt5CXMtPa7yOeT2MJT2NVqe1DLhCyUcijvpOSiJSr6lPAsdGKZRhGIeL3IHvLeFV9BSLxs5CTE5P0Do9nZSF7pb2Ue5F9x6aG6vkndTVUl9E/OsHEZCrXYhk5JpsYcodr6fcj4OcikiDoFGQYhrEPHYkRKkuLWVETWHLlJcWsrKuYUtRxYWoPchYKeXlNGSJLrZDnX8fa46+xBhPxJ5ta1he7w+tE5AFgGXB/pFIZhlGQtPcM09JQuU/mcmtD1VRsOS5kWxQEoLS4iMaqsli4rBekkH21rqExU8gxJ5ukrn8WkVcBqOqvVfVuVbWAhGEY+9GeGJmKG3taGivpiJnL2lu7TTVzx5DBV+tauqSu3uH5N5bwNE5V67K9yHEnmxjyJuAfRGSniHxBRDZELZRhGIWHqtLhLOQwLQ1V7O4fZWwiPjHM+VjIsPTlMxfSWMJjDSYKhzkVsqp+R1UvIMi0fhq4QUSeiVwywzAKir6RcQaSE1MJXZ7WhkpUobM3Pm7rroEkIrA8y329zbUVSx5DXkhRELAGE4XEfNovHg0cBxwOPBWNOIZhFCq+ZnVr474W8vTWp/i4rbsHRlleXUZJcXb/BTbXBRZyKqVzT46AhTSW8FiDicIhmxjyjc4ivh54Atigqm+MXDLDMAoK39UpvQHClEKOUZOJ7oHknE0lwjTVlDOR0qmKWfkmMTS2oIQuCBpMVJYWW3GQAiCbbU87gTNUdU/UwhiGUbhM7UFOS+paWVdBabHEqg1j10B2Vbo8PtbcNTC6YEt1MfQMjS04hgzQUFVqDSYKgGxiyP/XlLFhGHPR3jNCXUUJyyr3VRzFRcKq+spYVevq6k9mtQfZM10+M/9x5MU0lvA0VJeZy7oAmE8M2TAMY0baE8P7WceelobK2FTrSqWUPYPzU8jNdYF7eykyrRfTWMLTaAq5IDCFbBhGTmjPsOXJExQHiYdCTgyPMZHS+bmsa73LOv8KeTGNJTzWYKIwmDOGLCKHZRpX1T/mXhzDMAoRVaUjMcJrjm3OeL61sYo9g2MMj01QVZZV19fI6Joqm5l9Uld1eQlVZcV0LUFxkMU0lvBYg4nCIJt/GfcCCghBP+QjCPYjr49QLsMwCojuwSTJidSsLmsIal0fc0htPkXbj/kWBfE0L1FxkMU0lvCEG0xku9XLyD/ZJHW9TFVPdL/XAqcBD0UvmmEYhcJMe5A901uflusx0/gAACAASURBVN5tPV02c74KeWmKg/jY72JjyMCSbdsysmPeX5VU9VHg9AhkMQyjQOlIa7uYjh+PQ5OJhVrITbXl7FkChbyYxhIeH3+2OHK8ySaG/JHQxyLgFKAzMokMwyg4vOWbXhTEs6KmjIrSophYyKPUlJfMO5bdVFvOr7cvTVJXZWkxFaXzbyzh8fFniyPHm2z+RoYDPhMEMeU7oxHHMIxCpL1nhBU15TN2IxIRWhqqYlE+c75FQTzNdeUMJifynpiWGB5fdDESH3+2jk/xJpt+yJ/KhyCGYRQu7YmZtzx5WhsqY1E+s3uBCtnHnLsHkhy+PJ8KeWyqY9NCsXrWhUE2taw3iMhdIvKoiGzxP/kQzjCMwmC2oiCe1sZ4WMgLVci+OEi+E7t6hhbeWMLTYC7rgiCbpK5bgP8HXAK8MfSzaETkOBF5SESSIvKxWeZ9W0SeE5HN7uckN/529wXhCRH5nYi8PHTN+SLytIjsEJFrQuNHiMjv3fhtIlLmxq9y62wWkf8RkXVu/HUissmd2yQi54bWOtWN7xCRfxYRCZ37GxF5SkTaRORGN7ZGREZCz3FTaP5l7lnaROSGtOe/VESedOf+MzR+oxvbFr6/iFzh5NoiIveLyAo3/lY3P5Wpr7WIHCYig7P9WRhGOhOTKXb3jtI6p4VcxcDoBH1L7Dbt6h+dV5Uuz1KVz+wdXnhjCY9vMNFrFnKsyUYhd6vq3ar6nKq+4H9ydP8e4IPAF7KY+3eqepL72ezGngPOVtWXAZ8GvgEgIsXA14E3AOuAK7yCBW4AvqyqRwMJ4N1u/D/d1q6TgBuBL7nxPcAb3T2uBL4XkunfgL8C1rqf8939XwNcCLxcVdenPd/O0HNc5eYvBz4PnOfmrxSR89y5tcAngD9x5/7Wjb8K+BPgROAEgn7VZ4tICfBV4DWqeiKwBbja3Xsr8GbgwRne8ZeAn8xwzjAy8mL/KBMpzcJCDhT2UlrJQ8kJhsYm51UUxOOzsrvzXBwkFxYyBFufrMFEvMlGIf+jiHzTWV1v9j+5uLmqdqnqI8CC/pao6u9UNeE+Pgy0uOPTgB2q+qyqjgG3Ahc6C/Jc4A437zvARW6t/tDS1QTFUFDVx1TVZ5W3AZUiUi4ihwJ1qvqwqirwXb8W8D7gc6qa9M85x6McCTyjqt3u838TeCQgUPhf988ZWksJCrWUAeVAKfASQQEXAard89bhsuJVdZuqPp1JABG5iOALTtscshp55po7t/C9h55fajFmZGoP8gwZ1p6Wqa1PS6eQ/ZanhbisG6vKKC6SvLqsx11jicXGkAHqq0othhxzslHI7wJOIrD+vLv6z6MUagY+61ywXxaRTP+a3s20dbcaaA+d63Bjy4FeVZ1IGwdARD4gIjsJLOQPZrjHJcCjTtGudten3wPgGOBM5xr/tYi8IjTvCBF5zI2f6cZ2AMc6l3YJgWJvDa11jIj8VkQeFpHzAVT1IeABYLf7+alTuOMEXwieIFDE64BvZXiWKUSkBvh7YM4EPhF5r4hsFJGN3d3dc003FklyYpI7NnVw12O7llqUGZluuzi3yxqWti9y96Avmzl/hVxUJKyoKctrta7e4cU3lvAEFrIp5DiTTargK1T12MglmZ1PAC8SWIPfIFAe1/uTzkX8buDVi7mJqn4d+LqIvA34BwIXtb/HegJ39+uzWKoEaAReSeBKvl1EjiRQnIep6l4RORX4kYisV9WEiLwPuA1IAb8DjgqttRY4h8AD8KCIvAxYARzPtFfg507BP0ygkE8GngX+heD9fWYWea8jcOMPhsLgGVHVb+BCAxs2bNAs3oWxCLa/OMhEStm2e4DJlFJcNPufz1LQ0TNMkcChy2ZXyMuqSqmtKFlSl7WP/863KIgn39W6fMx3sTFkv0Yc9oEbM5ONhfy7UPx10Tgr1Cc1rcrmGlXdrQFJggSz00LrnQh8E7hQVfe64V1MW5gQKK1dwF6g3lmh4fF0bmXa/YyItAB3Ae9Q1Z2he7SErgmv1QH80Mn8BwIlu0JVk15GVd0E7CSwgFHVe1T1dFU9g6BW+PbQWner6riqPufG1wIXAw+r6qCqDhJ4B84g8GagqjudK/124FUZX+w0pwM3isjzBDHq/09Erp79EiMftHX2ATAyPslze4aWWJrMtCdGWFlXQVnJ3P+dtDZULalS8M0h5ls209NcW55XhZyLKl2eBmswEXuyUcivBDa7jGWf0bzgbU+q+vVQUlNWFb9cvBYXE72IIDnJd6L6IfCXqro9dMkjwFoJMqrLgMsJlJoSuHnf4uZdCfyXW2tt6Po/A55x4/UExVCuUdXfhp5jN9AvIq90cr3DrwX8CHiNu/4YAst+j4g0uYQznMW8lsCKRUSa3e8G4P0EXzL8Wue4cysIFPizwB9xSVwiUgqcDWwj+FKwTkSa3PWvc+MzoqpnquoaVV0DfAX4P6r6tdmuMfLDVqeQYVo5x42OxDAtcyR0eVobK2lfwvKZ3QNJSopkwQquKc8NJnzMdzGNJTzhBhNGPMnGZX1+VDcXkZXARoLEo5SI/C2wTlX7ReQ+4D1Oad/iFIwAm4Gr3BLXEsSF/9W5WidUdYOqTjgL76dAMXCzqvpkpb8HbhWRzwCPMR1fvVpEXkuQYJZg2l19NXA0cK2IXOvGXu+Sq94PfBuoJLBQfQz7ZuBmEdkKjAFXqqqKyFnA9SIyTmA1X6WqPe6ar8r0tq3rQ18wfgq8XkSeBCYJss33isgdBAlqTxAkeN2vqve49/opAtf2OPAC8E43fjGBC7sJuFdENqvqn87xx2QsIW2d/Ww4vIEtu/po6+znwpNWz31RnmnvGeFPjl6R1dzWhip+vb0bVWWu8EgUdA0kWVFTTtECXf/NteXsHUrmrWtSIscxZAgaTKxYoIfAiJZsKnXlaotTprVfZF+3b/jcBaHjc2eY8x7gPTOcuw+4L8P4s4Rc3qHxD82wzmeYIf6qqhsJthylj48Bf5Fh/E5mKDuqqlfMMK7AR9xPeHwS+OsZrrkJuCnD+F0ErvcZUdXrZjtv5I/JlLJtdz9vO+1wkhOpWFrIyYlJXhoYnTOhy9PaWMXoeIruweSCth4tlq6B5ILjxwBNdRWowt6hMQ6pi17+XLqsww0mTCHHE2uMaRgx5dnuQUbHU6xfVcf6VXVs3dVP8P0sPuxKjKA695YnT7gv8lLQPZBcUIa1J1w+Mx/korGExxpMxB9TyIYRU9o6g63xJ6xexvrVy+gbGWdX79LXgg7j48FzFQXxLHVf5O6B0QXtQfZ467orT8VBeoZzUxQEwg0mTCHHlfxVSDcMY160dfZRXlLEUU3VDI1NuLH+GVscLgVesWbrsl5KC3liMsXeoTGaFuEq99Z1vspn9g6P5yShC6ZjyNbxKb6YhWwYMaWts5/jVtZSUlzE8SvrKJJpqzkutCeGKS2WrOPBVWUlrKhZmv2we4fGUF1YURDPijy7rHuGFl/H2mMNJuKPKWTDiCGqytZdfaxbtQwImgMc1VRD2654JXZ1JEZYXV85r4IlS9UX2Vu1i3FZV5QWs6yyNG97kRM5aCzhqSgNGkwkTCHHFlPIhhFDOhIj9I9OcMLquqmx9avqYmchd/TM3XYxndbGqiUpn9k9GMR9F2MhQ373Iidy1FjC01hdZi7rGGMK2TBiiFe8652F7I9f7B9l72B+2//NRntiZN4x7daGSjp7R5hM5TdjfLps5uK2KwXVuqJP6splYwmPNZiIN6aQDSOGPNnZR3GRcNzK2qmx9asCazkuVvJQcoKeobGsE7o8LQ1VTKSUF/vz28bQu5lX1CzO4sxX+cxcNpbwWIOJeGMK2TBiyNbOfo5qqt5n/6m3lrfGpEDIVJen+VrIvi9ynhO7ugeS1FeVUl6yuD293mUd9Z7wRA4bS3gaqsrMQo4xppANI4a0dfZxQshdDUG3pJaGythYyFN9kOcbQ25Ymr3IXQOjC24qEaa5toLkROBOjhKffJXzGLJZyLHFFLJhxIzugSQv9SdZt6puv3PrV9XxZEwUcoezkP3e4mxZVV+JCHlvMrHYspkev0Z3xHFkb8nmOobcPzrBuDWYiCWmkA0jZvia1evTLGQ/9tyeIQZGlz5Ttr1nhMrSYpbP04IrKyni0LoKOpbAZZ2L+tneyo46jtwzFE0MGabj00a8MIVsGDHDu6QzWch+G9S23QN5lSkT7YlhWhsrF9S1qaUxv3uRVZWugeSi9iB7pi3kaBVyVDFkgF6LI8cSU8iGETOe7OznsMYqllXu76r0VnMcOj+19wzPO6HL09JQmde9yP0jE4xNpBa9BxmYKr0ZdfnMXDaW8Fi1rnhjCtkwYkZbZ9/UFqd0mmvLWVFTtuSJXapKR2Jk3gldntaGKl4aGCU5MZljyTLji4LkwkKuqyihrKSI7oj3g+eysYTHGkzEG1PIhhEj+kfHeX7v8IwKWURYv2oZW5e4hGbv8DiDyYl5J3R5WhurUIXO3vzsRc5F2UyPiAR7kSPeR50YGstZYwmPV/A+Pm3EC1PIhhEjtvkKXav3T+jyrF9Vx46uwbxZl5lon8qwXqiFnN+9yN6azUVSF7i9yBFbyInh8ZzGj2HaZW0WcjwxhWwYMWK6ZGZmCzk4t4yJlLL9xcF8ibUfHVN9kBduIQN5S+zKpYUMrlpX1DHkHDaW8FiDiXhjCtkwYkRbZz9NteWzWnLTJTSXzm093Qd5YRbyIXUVlBZL3hK7ugZGKS8poq4iNy3gm2sr8rDtKfcxZHDlM81CjiWmkA0jRsyW0OU5rLGK2vKSJS2h2Z4YZlllKXUVC4txFhcJq+sr82Yhd7uiIAvZopWJptpy+kbGIwsbjE+mGBidyLmFDEFil+1DjiemkA0jJoyOT/JM1+B+JTPTKSoSjl/iVoztPSMLdld7WhqqplzfUdM1kMxJ2UyP3z4V1V5krzBzndQFQRzZtj3FE1PIhhETtr80wGRK57SQIXBbb9vdn/cWhp72xML3IHtaGyvzVq2rK0dVujy+OEhUbusoioJ4rMFEfDGFbBgxYeuu/Xsgz8QJq5YxOp7i2e78J3alUovbg+xpaahi79AYQ8lomzTAtMs6VzTVVEytGwU9ETSW8FgLxvhiCtkwYkJbZx+1FSVZuYLXr1663sh7BpOMTaQWvAfZ4xV61G7r0fFJ+kbGc+uyjthC7o2gsYSnvqqUAWswEUtMIRtGTGjr7Gf9qrqsEo+OaqqhrKRoSTKtF9oHOZ187UX2VmwuLeTl1WWIQHdExUGiaCzhsQYT8cUUsmHEgInJFE+92J+VuxqgtLiI41bWLomFPN0HOTcWctSZ1r6AR672IAOUFBexvLossuIgUceQw/cw4oMpZMOIAc/uGWJ0PJVVQpfHl9BUzW9il7doF1qly7O8uozK0uLI9yL7Ah65TOqCoMlEVMVBEkNjVJXltrGEx1vIVhwkfphCNowY4F3PJ8xSMjOd9avq6B+dyNvWIU97Ypim2vJFKwsRoaWhko48Wci56PQUJsrymT0RVOny+Li0WcjxwxSyYcSAtl39lJcUceSK6qyvma7YlV+3dXvPyFT8d7G0NlbRHvEXiu7+UURyH4+NsnxmFI0lPNZgIr6YQjaMGNDW2c9xh9ZRUpz9P8njD62juEh4Ms+JXe2J4UVvefK0NgR7kaN0u3cNJFleXT6vd5sNzbXl7BlMkopgL3hPBI0lPBZDji95UcgicrOIdInI1hnOv11EtojIEyLyOxF5eejc8258s4hsDI1/2l2zWUR+JiKrQufOceNtIvLr0PiH3dhWEfm+iFS48XNF5FE3/h0RKclCrg+5+W0i8reh8etEZJe7/2YRucCNLxeRB0RkUES+lvb8l7n7tInIDaHxw9w1j7nzfq01IjISusdNoWs+KyLtIjKYdo/DReQXbp1fiUiLGz9JRB5y994iIpfN/Sdq5BJVzapkZjoVpcUc1VTN1jxayBOTKXb3jS56y5OntbGKgeQEfSPRWWvdA8mcu6shcFlPpDQSxdYbocvaGkzEl3xZyN8Gzp/l/HPA2ar6MuDTwDfSzr9GVU9S1Q2hsc+r6omqehLwY+BaABGpB/4VeJOqrgfe6sZXAx8ENqjqCUAxcLmIFAHfAS534y8AV84ml4icAPwVcBrwcuDPReTokGxfdvKepKr3ubFR4H8DHws/mIgsBz4PnOfkXSki57nT/wDcrqonA5e75/LsDN3jqtD4PU6udL4AfFdVTwSuB/7JjQ8D73D3Ph/4inuHRp7oSIzQPzoxZ8nMTKxftSyvW592940ymdJFb3ny+MSwKBO7ugaSOc2w9vgksSj2IkfVWMJjDSbiSV4Usqo+CPTMcv53qppwHx8GWrJYM2wWVAPeb/Q24Ieq+kc3rys0rwSodBZwFdAJLAfGVHW7m/Nz4JI55Doe+L2qDqvqBPBr4M1zyDukqv9DoJjDHAk8o6rd7vN/+/u7Z/Jm0zIn76yo6sOqujvDqXXAL93xA8CFbv52VX3GHXcCXUDTXPcpRPqGx9nc3rvUYuyHV6jztZD9NS/1JyOrGJXO1B7kHLmsvaUd5danroHRSCzkqIqDRNlYwtNQXWoWcgyJYwz53cBPQp8V+JmIbBKR94Ynevcs8HachQwcAzQ4t+wmEXkHgKruIrAS/wjsBvpU9WfAHqBERLz1/RagdQ65tgJnOjd0FXBB2jVXO/fvzSLSMMfz7gCOdW7oEuCi0FrXAX8hIh3AfcDfhK47wrmyfy0iZ85xD4DHmf7ScDFQ66zzKUTkNKAM2DnTIiLyXhHZKCIbu7u7Z5oWSz5z75O85d9+R1dExRwWytZd/RQXCceurJ33tX7fcr6s5A6/BzlHFvLUXuSIioOkUsqewbGcFgXx+Mpfuf4yFGVjCU9Qz9qSuuJGrBSyiLyGQPH9fWj41ap6CvAG4AMicpY/oaqfVNVW4BbgajdcApwK/Bnwp8D/FpFjnGK8EDgCWAVUi8hfaJBNcjnwZRH5AzAA7NNTLV0uVd0G3AD8DLgf2By65t+Ao4CTCBT/F2d7ZmeBvw+4DfgN8HxorSuAb6tqC4HS/55zse8GDnOu7I8A/ykic5lXHwPOFpHHgLOBXeHnFJFDge8B71LVGWvqqeo3VHWDqm5oaiocQ3owOcGPt+xmIqX88LFdSy3OPrR19nF0U82CthGty3OmdXtimCKBQ+tzs6c3aOFYEtnWrZ7hMSZTmtOymZ5pCzm3X/CiLArisQYT8SQ2CllETgS+CVyoqnv9uLNsvev5LjLHR29h2s3bAfzUuYj3AA8SxHlfCzynqt2qOg78EHiVW/shVT1TVU9z8737eja5vqWqp6rqWUDCX6OqL6nqpFNq/z6DvPugqveo6umqegbwdOj+7wZu9zICFcAKVU16WVR1E4FFe8wc9+hU1Tc7Jf5JN9brnrEOuBf4pKo+PJe8hciPH+9kZHySQ+rKuf2R9rwX05iNts7+qdrU82VZZSmtjZU8mS+F3DPMocsqKc1hxnKw9SkaC3mqKEhdbouCAFSVlVBTXpLzrU9RNpbwWIOJeBILhSwihxEoyL8MxXIRkWoRqfXHwOsJ3MWIyNrQEhcCT7nj/wJeLSIlzp18OrCNwFX9ShGpEhEBznPjiEiz+11OYAXfNJtcadccRuAK/k/3+dDQtIu9vHM8v1+rAXg/wRcAnMznuXPHEyjkbhFpEpFiN34ksBZ4do57rHDWNcAngJvdeBnBF53vquodc8laqNy2sZ21zTV89PXH8uyeITa9kJj7ojzQNTBK10Ay65KZmTghj4ldQZen3GRYe1obqiJzWUdRNjNMFMVBfGw3agvZGkzEj3xte/o+8BBBrLRDRN4tIleJiM8OvpYguepfZd/tTYcA/yMijwN/AO5V1fvduc+5bUdbCBT1h2DKnXw/sMVd801V3aqqvwfuAB4FniB4dp/N/Xciss1dc4+q+uSnmeQCuFNEniTIav6AtzaBGyXYJrUFeA3w4dB7eB74EvBO9x7WuVNfdWv9FvhcSPl/FPgr9/zfB97pXOxnAVtEZLN7pqtUtcfd40YXc65y97jOrXUO8LSIbHfv9bNu/FK33jtD26hO4gDimZcGeOyPvVy6oZU/e9mhVJcVc9sj7UstFjDtal5IQpdn/ao6nt87TP9o9DHB9sTwoktmptPaWElHYiQSr4XPF4giqQucQs6xhZzIRwzZrW0NJuJFST5uoqpXzHH+PcB7Mow/S+BuznTNJZnG3bnPE2wlSh//R+AfM4z/HfB32crlzmVMpFLVv5xFrjUzjGd8P6r6JPAnGcbvBO6c4ZqPAx/PMH4HgfJOH/8P4D9mkvlA4PaN7ZQUCRefsprq8hLe+PJV3P14J//4pvXUlOfln8CMeFfzukUp5MC63tbZz+lHLp9j9sIZHZ/kpf5kzhK6PK2NVSQnUq5ncW5dy/mwkHMdLshXDNnfK6p3Y8yfWLisDSMqxidT/PDRXZx3fDMrXGLPWze0Mjw2yb1b5txFFjltnX0cvryKuoqFW0P5KqG5qzc3XZ7SiXLrU1d/kpryEqrKovniFZTPzG1SV0+EjSU80+UzLY4cJ0whGwc0v9jWxd6hMS57xfSutFMOq+fo5ppYuK19D+TF0FxXQVNteeQK2cd5c7UH2dMaYXGQqKp0eZprKxgam2QoOZGzNRMRVuny+AYTvZZpHStMIRsHND/Y2M4hdeWctXZ6i5aIcOmGFh79Yy87ugaWTLb+0XFe2Du8qIQuz/pVdZEndvkmELl2WfuYdBRdn7ojqtLl8Wvnci9ylI0lPNZgIp6YQjYOWF7qH+WBp7u45JSW/RoLXHxyCyVFwu0bO5ZIuun48WItZL/GM12DjI5Pzj15gXQkhikrKcq5xVlZVsyKmvJILOSugdFIFbJ/F7ms1hVlYwmPNZiIJ6aQjQOWOx/tIKVw6Yb9C6811ZZz7nHN/PDRjiXb+jGdYZ0LC3kZkyll+0vRWfwdPSO01FdSVCQ5X7u1sTKaGPJAcqrmdBREURwkEXEdawgaTFSVFVsMOWaYQjYOSFSVH2zs4LQjGlkzQ4/hy17Ryp7BMX75VFfG81HTtquP5trynFhwvjHF1l3RxZHbE8OszlGXp3RaG3JfHGQoOcHw2GS0LusIymfmI4YMVq0rjphCNg5I/vBcD8/tGeKyDNax5+xjmmiuLecHG5cmuSsXCV2e1sZKaitKIo0jt/fkrg9yOi0NlXT2jjKRQ2+FdyNHmdTVUFVGSZHkzGWdj8YSHmswET9MIRsHJLdv7KCmvIQ3vGzljHNKiou45NQWHni6O+8NJ0bHJ9nRPcgJqxfvroYgUW3doXWRZVoPJidIDI/nPKHL09pYxWRK2d2Xuz8Hb7VG0VjCU1QkQXGQHClkb7E2RpzUBcGXiR4rDBIrTCEbBxwDo+Pc98Ru3vjyVXPuP710QyuTKeXOR/PbcOLpFweYTGnOLGSAE1Yv46kX+3NqZXqmtzxF57KG3O5F9nHdqAtfNNWW58xC9pWz6vPksrZtT/HCFLJxwPHjLbsZGZ/k0g1zttXmiBXVnLamkR9szG/Dia1TPZBzYyEHa9UxOp7i2T1DOVvTM6WQI7OQA0Wfy65PU40lIkzqCtbPXXGQfDSW8FiDifhhCtk44LjtkXaOOaSGk1rrs5p/6StaeXbPEBvz2HCirbOfuoqSqSpVuSDK3sheUUYVQ15VX0mRQEcOm0x0DyYpLRbqK6N1/zbVlrMnRw0m8tFYwmMNJuKHKWTjgGL7SwNsbg8aSQRNvebmgpetpKa8JK+Vu4KErmVZy5gNRzVVU15SRFsEmdbtiWGqyoppqIpGuZUWF3Hossqp4iO5oKs/yYqa8ki2aYVpqq1g79BYTkIFPb6OdT5iyO4elmkdH0whGwcUtz/STmmxcPHJq7O+pqqshDe+/FDu3bKbwRyWQJyJickUT+3OXYa1p6S4iOMiSuxq7xmhtaEqp18g0mlpqMxpG8augdFIM6w9zbXlqMKewcUrNh9DzpeFHL6nsfSYQjYOGMYmUtz12C5ee/whLK+Z33/Eb93Qysj4JD9+PPqGEzu7h0hOpFi/OrcKGaZLaOY6Ht6RGI4socvTkuO9yEHZzGjjx5Db8pn5aCzhsQYT8cMUsnHA8MunXmLv0FjGylxzcXJrPWuba7g9D3uSfYz3hBwmdHnWr6qjf3Qip8lRqkp7T+77IKfT2ljJS/3JnJX/jLqOtWe6fObiE7sSQ/kpCgKh8pmmkGODKWTjgOG2R9pZWVfBWcc0zT05jaDhRGteGk60dfZTUVrEkU01OV/7hAgSuxLD4wyNTUaW0OXxGdydvYv/MjE+mWLv0Fh+XNauh3Mutj4lhqNvLOGZjiGbyzoumEI2Dghe7Bvl19u7ueTU1RQvMInn4lNW56XhxNZdfRy3sm7Bcs7GsStrKS6SnJbQ9F2YWiMqm+nxCj8XiV17XTw3HxbyiprA0syJyzoPjSU81mAifphCNg4IfCOJt546f3e1Z0VNOecdH23DCVXlyQgSujwVpcUc3VSTUwvZd2GK3EJ2MepcJHZ593E+LOTykmLqq0pz5rLOxx5ksAYTccQUslHwqCq3b2zn9FkaSWSLbzjxi23RNJxo7xlhYHQiZyUzMxEkduXOQvaJVrncM52JQ2orKC2WnCR2TZfNjD6pCwLFnwsLOV+NJTzWYCJemEI2Cp7fP9fDC3uHuewVC7eOPWetjbbhRNtUha5oLGSA9auX0TWQzFlLwPaeYeqrSqmtiDa2WVQkrK6vpCMHfZF9PDcfLmt/n8XGkPPZWMJjDSbihSlko+C5fWM7teUlvOGEQxe9VklxEW85tYUHnu7ipQgaTmzt7KO4SDjmkNqcr+3xyj5XVnJ7YiSykpnptDbmZuuTL5vZNM/tbwulubZi6p4LJZ+NJTzWYCJemEI2Cpp+30jipFVUluVm7+ZbN7SS0iAunWvaqQ4RzwAAD2FJREFUOvtZ21wT6T7TdU4hP5kjhdzRE/0eZE9LQ1VOYsjdg6M0VJVSVpKf/+KaasvpHkwuav93YsgVBclTDBmCvchmIccHU8hGQfPjx3czOp5a0N7jmThiRTWnHdHIDzZ25LzAhi+ZGSV1FaUcvrwqJ4ldqZTSkVcLuZLE8PiiK6Z19ednD7KnubacsYkU/SMLl9tbyBZDPngxhWwUNLdtbOfYQ2p5eUtuldylG1p5bs8Qjzyfu4YTXf2jdA8kI40fe3KV2NU9mGRsMkVLxBnWHq/4Oxbptu4aSEbe5SlMUw6Kg+SzsYTHGkzEC1PIRsHy9IsDPN7ey6WvyL6RRLZE0XDCK8j8KORlvLB3mP7RxcUHp9su5sdlPbUXeZGJXfmq0uXJRfnMnuH8tV70NFqDiVhhCtkoWG7fOP9GEtniG07c98RuBhap1DzehbwuDwo5V3Hk6S1P+bKQF78XWVXpHkjmZQ+yx1vji8m09k0e6iPqqJWJemswEStMIRsFiW8k8bp1h0RmUVzqG05s2Z2T9do6+1mzvCry7UMQLqG5SIXsLNWo9yB7GqvLqCwtXlSmdd/IOGOTqfzGkOsW77LOZ2MJjzWYiBemkI2C5BfbXqJnaIy35jCZK52TctxwYmtnX+QJXZ6m2nKaa8tp27W4xK72nmGaa8vzpiREhNbGykW5rLvzvAcZoLa8hPKSokW5rPPZWMJjDSbihSlkoyC5baNrJLF2/o0kskVEuOwVrTz2x16eeWlxDSf6RsZp7xnJi7vak4vErvbEcOQlM9NpbahaVFKXdxvnM6lLRGiuW1xxkJ7h/JXN9ExZyBZDjgWmkI2CY3ffCA9u7+Ytp7ZE0qAhzEUn+4YTi7OSfSw3ypKZ6axftYwd3YOLamcYbHnKj7va09pYRUdiZMFbzqbLZubPQobgC8CiLOTh8bzGj2E6Xm0x5HhgCtkoOO7c5BpJbGiJ/F4rasp57fGH8MNHdzE2sfCtIfkomZnOCavrmEwpT7+4MOt+YjLF7r7RvFvILQ2VDCYnFqwkfBw3ny5rCKqCLcZCzmdjCY81mIgXppAPAkTkfBF5WkR2iMg1Gc6Xi8ht7vzvRWRN/qXMjlRKuX1jB688spHDly+ukUS2XPqKFvYOjfHLpxbecKKts59D6spZkadSjsBUvHrrAguE7O4bZTKleSsK4pluw7gwt3VXf5KK0iJqy0tyKdacNNeV07WIcqtLEUMGVxzEFHIsMIV8gCMixcDXgTcA64ArRGRd2rR3AwlVPRr4MnBDfqXMnt8/18Mfe3LTSCJbzlrbxCF15YtyW7flMaHL09JQSV1FyYLjyH7rUb4yrD0tU1ufFpbY1T0YFAXJ9d70uWiuLad/dGJBIYLxyRQDyfw2lvA0VJfaPuSYkN+vkMZScBqwQ1WfBRCRW4ELgSdDcy4ErnPHdwBfExHRXNeNdLzxX/5nwXHNnqExastLOH/94htJZItvOPFvv9rJ67706wWtsaN7kPPXr8yxZLMjIqxftYwfPbaLR57rmff1A6NBGci8J3W5+113Txtf+e/t875+V+8Ixx+av9CAx7vIL/jqb+ad2zDp/qk15LGxhKehqoyHnt274L/bByMNVWXcftUZOV/XFPKBz2ogbNp1AKfPNEdVJ0SkD1gO7ElfTETeC7wX4LDDDluQQEc1VTO2iFJ95x53SM4aSWTLlWesobN3lOTEwr5IHH9oHRdFUMBkLt579pGLaiW5sq6S1fX5tZDrKkr5m3OPZmf34IKuX3tIDX9+4qocSzU3Zx/TzMUnr17w35ETVy/j3OOacyzV3Fx5xhpqK0wVzIe6iGoJSERGkBETROQtwPmq+h73+S+B01X16tCcrW5Oh/u8083ZTyGH2bBhg27cuDE64Q3DMA4wRGSTqm7IdM5iyAc+u4BwwLXFjWWcIyIlwDJgb16kMwzDMABTyAcDjwBrReQIESkDLgfuTptzN3ClO34L8Muo4seGYRhGZixwcIDjYsJXAz8FioGbVbVNRK4HNqrq3cC3gO+JyA6gh0BpG4ZhGHnEFPJBgKreB9yXNnZt6HgUeGu+5TIMwzCmMZe1YRiGYcQAU8iGYRiGEQNMIRuGYRhGDDCFbBiGYRgxwAqDGAtGRLqBFxZ4+QoyVAKLASbX/DC55ofJNT8ORLkOV9WMjdxNIRtLgohsnKlazVJics0Pk2t+mFzz42CTy1zWhmEYhhEDTCEbhmEYRgwwhWwsFd9YagFmwOSaHybX/DC55sdBJZfFkA3DMAwjBpiFbBiGYRgxwBSyYRiGYcQAU8hGpIjI+SLytIjsEJFrMpwvF5Hb3Pnfi8iaPMjUKiIPiMiTItImIh/KMOccEekTkc3u59pMa0Ug2/Mi8oS758YM50VE/tm9ry0ickoeZDo29B42i0i/iPxt2py8vC8RuVlEukRka2isUUR+LiLPuN8NM1x7pZvzjIhcmWlOjuX6vIg85f6c7hKR+hmunfXPPAK5rhORXaE/qwtmuHbWf7sRyHVbSKbnRWTzDNdG+b4y/t+Qt79jqmo/9hPJD0G7x53AkUAZ8DiwLm3O+4Gb3PHlwG15kOtQ4BR3XAtszyDXOcCPl+CdPQ+smOX8BcBPAAFeCfx+Cf5MXyQobpD39wWcBZwCbA2N3Qhc446vAW7IcF0j8Kz73eCOGyKW6/VAiTu+IZNc2fyZRyDXdcDHsvhznvXfbq7lSjv/ReDaJXhfGf9vyNffMbOQjSg5Ddihqs+q6hhwK3Bh2pwLge+44zuA80REohRKVXer6qPueADYBqyO8p455ELguxrwMFAvIofm8f7nATtVdaEV2haFqj5I0LM7TPjv0HeAizJc+qfAz1W1R1UT8P+3d/cxdlR1GMe/T1+MsYW1AkHlJYBBFBRabRteutiYQsAXkvpKS7CIiWkiNvYfaWwiWBMFNDQgIAkvEdsVjIK1VoxVaGsFKy+tXaqt0hgSIHWrRFsL2sjy849zhp3e3rt72+7MvdDnk9zcmTNnZs49d+6cPWdm58evgIuqLFdErI6Il/PsBuD40drfoZSrTe38dispV/79fwq4d7T2165hzg21HGNukK1KxwHPluafY/+G79U8+eS1CziqltIBeYh8CvD7JovPkbRZ0i8knVFTkQJYLelJSZ9vsrydOq3SpbQ+UXaivgCOjYgdefpvwLFN8nS63q4kjWw0M9J3XoWr8lD63S2GXztZX73AQEQ83WJ5LfXVcG6o5Rhzg2yHLUkTgfuBL0XE7obFG0nDsmcB3wFW1FSsGRHxPuBi4AuSzq9pvyOS9AbgEuBHTRZ3qr72EWnssKv+l1PSYuBloK9Flrq/8+8C7wAmAztIw8PdZA7D944rr6/hzg1VHmNukK1KzwMnlOaPz2lN80gaB/QAL1RdMEnjST+4voh4oHF5ROyOiD15+kFgvKSjqy5XRDyf33cCPyENHZa1U6dVuRjYGBEDjQs6VV/ZQDFsn993NsnTkXqTdAXwEeCyfCLfTxvf+aiKiIGIGIyIV4A7WuyvU/U1DvgY8MNWeaqurxbnhlqOMTfIVqXHgVMlnZx7V5cCKxvyrASKuxE/ATzc6sQ1WvI1qruArRFxY4s8by2uZUuaTvqtVPqHgqQJko4opkk3BW1pyLYS+IySs4FdpaG0qrXsuXSivkrKx9A84KdN8vwSuFDSpDxEe2FOq4yki4AvA5dExEst8rTznY92ucr3HMxusb92frtVmAVsi4jnmi2sur6GOTfUc4xVcaeaX34VL9JdwX8h3bG5OKctIZ2kAN5IGgLdDjwGnFJDmWaQhpz6gT/k14eA+cD8nOcq4I+ku0s3AOfWUK5T8v42530X9VUul4Bbc30+BUyt6XucQGpge0pptdcX6Q+CHcD/SNfoPke65+Ah4Gng18Bbct6pwJ2lda/Mx9l24LM1lGs76ZpicYwV/03wduDB4b7zisu1LB87/aSG5m2N5crz+/12qyxXTv9ecUyV8tZZX63ODbUcY350ppmZWRfwkLWZmVkXcINsZmbWBdwgm5mZdQE3yGZmZl3ADbKZmVkXcINsZq95StGmVjVJnyrp5lKec0dxnydJmttsX2YHY1ynC2Bmrx+SxkbEYKfLUYiIJ4AiRN9MYA/waLvrSxoXQwEiGp0EzAV+0GRfZgfMPWQzG1HuDW6T1Cdpq6QfS3pTXvaMpOslbQQ+KWlOjle7RdL1pW3skbQ0x5l9SNIxOX2ypA0aihs8KacvyHFp+yXdl9OmS/qdpE2SHpV02gjlnilpVQ4UMB9YqBRHt1fSMZLul/R4fp2X17lW0jJJjwDL8mdfL2ljfhW97OuA3ry9heVeulL83BW57BsknVna9t2S1kr6q6QFo/Ud2WufG2Qza9dpwG0R8W5gNymWdeGFSA/8/w0p9u8HScELpkkqQtVNAJ6IiDOAdcA1Of37wNURcSbpCVJF+iJgSk6fn9O2Ab0RMQX4KvCNdgoeEc8AtwNLI2JyRKwHbsrz04CPA3eWVjkdmBURc0jPLb4gf75PA8Ww9CJgfd7e0oZdfg3YlMv+lfwZC+8iheqbDlyTn51s5iFrM2vbsxHxSJ5eDiwAvp3ni2AA04C1EfF3AEl9pGD0K4BXSvmWAw9I6gHeHBHrcvo9DEWT6gf6JK1gKHpUD3CPpFNJjzg8lMZsFnC6hsJvH6kU5QdgZUT8J0+PB26RNBkYBN7ZxrZnkBp5IuJhSUdJOjIv+3lE7AX2StpJCuXX9NnNdnhxg2xm7Wp8zm55/sVR2F6jD5Ma848CiyW9F/g6sCYiZudh6LUHsd/CGODsiPhvOTE30OXPsxAYAM7K6+yT/yDsLU0P4vOwZR6yNrN2nSjpnDw9F/htkzyPAR+QdLSksaQIUUXvdwwpoter60fELuCfknpz+uXAOkljgBMiYg1wNalnPDG/FyHtrjjA8v8bOKI0vxr4YjGTe8DN9AA7IoUrvBwY22J7ZeuBy/J2ZwL/iP1jbpvtww2ymbXrz6SA8FuBSaRA9/uIFApyEbCGFJHnyYgoQtW9CEyXtIV0jXlJTp8HfEtSP+m68xJSo7dc0lPAJuDmiPgXcAPwTUmbOPCe5c+A2cVNXaQh96n5xqs/MXSdutFtwDxJm0nXf4vecz8wKGmzpIUN61wLvD9/pusYCt1n1pKjPZnZiPLw8KqIeM8hbGNPREwcOafZ4ck9ZDMzsy7gHrKZmVkXcA/ZzMysC7hBNjMz6wJukM3MzLqAG2QzM7Mu4AbZzMysC/wfKOxWDwWsb6AAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" } ], "source": [ - "import sys\n", - "\n", - "from io import StringIO\n", - "\n", - "x = np.random.randn(100)\n", + "rng = np.random.default_rng(RANDOM_SEED)\n", + "y = rng.normal(loc=5, size=20)\n", "\n", "old_stdout = sys.stdout\n", "mystdout = sys.stdout = StringIO()\n", "\n", "with pm.Model() as model:\n", - " mu = pm.Normal(\"mu\", mu=0, sigma=1)\n", - " sd = pm.Normal(\"sd\", mu=0, sigma=1)\n", - "\n", - " mu_print = tt.printing.Print(\"mu\")(mu)\n", - " sd_print = tt.printing.Print(\"sd\")(sd)\n", + " mu = pm.Normal(\"mu\", mu=0, sigma=10)\n", + " a = pm.Normal(\"a\", mu=0, sigma=10, initval=0.1)\n", + " b = pm.Normal(\"b\", mu=0, sigma=10, initval=0.1)\n", + " sd_print = Print(\"Delta\")(a / b)\n", + " obs = pm.Normal(\"obs\", mu=mu, sigma=sd_print, observed=y)\n", "\n", - " obs = pm.Normal(\"obs\", mu=mu_print, sigma=sd_print, observed=x)\n", - " step = pm.Metropolis()\n", - " trace = pm.sample(\n", - " 5, step, tune=0, chains=1, progressbar=False\n", - " ) # Make sure not to draw too many samples\n", - "\n", - "sys.stdout = old_stdout\n", + " # limiting number of samples and chains to simplify output\n", + " trace = pm.sample(draws=10, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED)\n", "\n", - "output = mystdout.getvalue().split(\"\\n\")\n", - "mulines = [s for s in output if \"mu\" in s]\n", + "output = mystdout.getvalue()\n", + "sys.stdout = old_stdout # setting sys.stdout back" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "'Delta __str__ = -85.74093608165128\\nDelta __str__ = -9.182002291671038\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.315734173890055\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.312485782438435\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.314669656412736\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.31581619157038\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.315114719133609\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.31511040479387\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.314077394936474\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.313673830463395\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.31561025339713\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = -9.31526569370057\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\nDelta __str__ = 0.10737295473067673\\n'" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "output" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Raw output is a bit messy and requires some cleanup and formatting to convert to `numpy.ndarray`. In the example below regex is used to clean up the output, and then it is evaluated with `eval` to give a list of floats. Code below also works with higher-dimensional outputs (in case you want to experiment with different models)." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", "\n", - "muvals = [line.split()[-1] for line in mulines]\n", - "plt.plot(np.arange(0, len(muvals)), muvals)\n", - "plt.xlabel(\"proposal iteration\")\n", - "plt.ylabel(\"mu value\");" + "# output cleanup and conversion to numpy array\n", + "# this is code accepts more complicated inputs\n", + "pattern = re.compile(\"Delta __str__ = \")\n", + "output = re.sub(pattern, \" \", output)\n", + "pattern = re.compile(\"\\\\s+\")\n", + "output = re.sub(pattern, \",\", output)\n", + "pattern = re.compile(r\"\\[,\")\n", + "output = re.sub(pattern, \"[\", output)\n", + "output += \"]\"\n", + "output = \"[\" + output[1:]\n", + "output = eval(output)\n", + "output = np.array(output)" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([0., 0., 0., 0., 0.])" + "array([-85.74093608, -9.18200229, 0.10737295, 0.10737295,\n", + " 0.10737295, -9.31573417, 0.10737295, -9.31248578,\n", + " 0.10737295, -9.31466966, 0.10737295, -9.31581619,\n", + " 0.10737295, -9.31511472, 0.10737295, -9.3151104 ,\n", + " 0.10737295, -9.31407739, 0.10737295, -9.31367383,\n", + " 0.10737295, -9.31561025, 0.10737295, -9.31526569,\n", + " 0.10737295, 0.10737295, 0.10737295, 0.10737295,\n", + " 0.10737295, 0.10737295, 0.10737295, 0.10737295,\n", + " 0.10737295, 0.10737295])" ] }, - "execution_count": 5, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "trace[\"mu\"]" + "output" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Notice that for each iteration, 3 values were printed and recorded. The printed values are the original value (last sample), the proposed value and the accepted value. Plus the starting value in the very beginning, we recorded in total `1+3*5=16` value above." + "Notice that we requested 5 draws, but got 34 sets of $a/b$ values. The reason is that for each iteration, all proposed values are printed (not just the accepted values). Negative values are clearly problematic." ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(34,)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "output.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Authors\n", + "\n", + "* Authored by Thomas Wiecki in July, 2016\n", + "* Updated by Igor Kuvychko in August, 2022 ([pymc#406] (https://github.com/pymc-devs/pymc-examples/pull/406))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Watermark" + ] + }, + { + "cell_type": "code", + "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "seaborn 0.10.1\n", - "pandas 1.0.4\n", - "pymc3 3.9.0\n", - "numpy 1.18.5\n", - "last updated: Fri Jun 12 2020 \n", + "Last updated: Tue Aug 02 2022\n", + "\n", + "Python implementation: CPython\n", + "Python version : 3.10.5\n", + "IPython version : 8.4.0\n", "\n", - "CPython 3.7.7\n", - "IPython 7.15.0\n", - "watermark 2.0.2\n" + "aesara: 2.7.5\n", + "xarray: 2022.3.0\n", + "\n", + "matplotlib: 3.5.2\n", + "aesara : 2.7.5\n", + "numpy : 1.23.0\n", + "arviz : 0.12.1\n", + "sys : 3.10.5 | packaged by conda-forge | (main, Jun 14 2022, 06:57:19) [MSC v.1929 64 bit (AMD64)]\n", + "re : 2.2.1\n", + "pandas : 1.4.3\n", + "pymc : 4.1.2\n", + "\n", + "Watermark: 2.3.1\n", + "\n" ] } ], "source": [ "%load_ext watermark\n", - "%watermark -n -u -v -iv -w" + "%watermark -n -u -v -iv -w -p aesara,xarray" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{include} ../page_footer.md\n", + ":::" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -343,7 +564,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.10" + "version": "3.10.5" } }, "nbformat": 4, diff --git a/myst_nbs/howto/howto_debugging.myst.md b/myst_nbs/howto/howto_debugging.myst.md index 62bae83e2..02ae6be03 100644 --- a/myst_nbs/howto/howto_debugging.myst.md +++ b/myst_nbs/howto/howto_debugging.myst.md @@ -6,108 +6,193 @@ jupytext: format_version: 0.13 jupytext_version: 1.13.7 kernelspec: - display_name: Python 3 + display_name: Python 3 (ipykernel) language: python name: python3 --- +(howto_debugging)= # How to debug a model -There are various levels on which to debug a model. One of the simplest is to just print out the values that different variables are taking on. +:::{post} August 2, 2022 +:tags: debugging, aesara +:category: beginner +:author: Thomas Wiecki, Igor Kuvychko +::: -Because `PyMC3` uses `Theano` expressions to build the model, and not functions, there is no way to place a `print` statement into a likelihood function. Instead, you can use the `Theano` `Print` operatator. For more information, see: theano Print operator for this before: http://deeplearning.net/software/theano/tutorial/debug_faq.html#how-do-i-print-an-intermediate-value-in-a-function. ++++ -Let's build a simple model with just two parameters: +## Introduction +There are various levels on which to debug a model. One of the simplest is to just print out the values that different variables are taking on. -```{code-cell} ipython3 -%matplotlib inline +Because `PyMC` uses `Aesara` expressions to build the model, and not functions, there is no way to place a `print` statement into a likelihood function. Instead, you can use the `aesara.printing.Print` class to print intermediate values. +```{code-cell} ipython3 +import arviz as az import matplotlib.pyplot as plt import numpy as np import pandas as pd -import pymc3 as pm -import seaborn as sns -import theano.tensor as tt - -x = np.random.randn(100) +import pymc as pm +``` -with pm.Model() as model: - mu = pm.Normal("mu", mu=0, sigma=1) - sd = pm.Normal("sd", mu=0, sigma=1) +```{code-cell} ipython3 +%matplotlib inline +%config InlineBackend.figure_format = "retina" - obs = pm.Normal("obs", mu=mu, sigma=sd, observed=x) - step = pm.Metropolis() - trace = pm.sample(5000, step) +RANDOM_SEED = 8927 ``` +### How to print intermediate values of `Aesara` functions +Since `Aesara` functions are compiled to C, you have to use `aesara.printing.Print` class to print intermediate values (imported below as `Print`). Python `print` function will not work. Below is a simple example of using `Print`. For more information, see {ref}`Debugging Aesara `. + ```{code-cell} ipython3 -trace["mu"] -``` +import aesara.tensor as at -Hm, looks like something has gone wrong, but what? Let's look at the values getting proposed using the `Print` operator: +from aesara import function +from aesara.printing import Print +``` ```{code-cell} ipython3 -with pm.Model() as model: - mu = pm.Normal("mu", mu=0, sigma=1) - sd = pm.Normal("sd", mu=0, sigma=1) +x = at.dvector("x") +y = at.dvector("y") +func = function([x, y], 1 / (x - y)) +func([1, 2, 3], [1, 0, -1]) +``` - mu_print = tt.printing.Print("mu")(mu) - sd_print = tt.printing.Print("sd")(sd) +To see what causes the `inf` value in the output, we can print intermediate values of $(x-y)$ using `Print`. `Print` class simply passes along its caller but prints out its value along a user-define message: - obs = pm.Normal("obs", mu=mu_print, sigma=sd_print, observed=x) - step = pm.Metropolis() - trace = pm.sample( - 3, step, tune=0, chains=1, progressbar=False - ) # Make sure not to draw too many samples +```{code-cell} ipython3 +z_with_print = Print("x - y = ")(x - y) +func_with_print = function([x, y], 1 / z_with_print) +func_with_print([1, 2, 3], [1, 0, -1]) ``` -In the code above, we set the `tune=0, chains=1, progressbar=False` in the `pm.sample`, this is done so that the output is cleaner. +`Print` reveals the root cause: $(x-y)$ takes a zero value when $x=1, y=1$, causing the `inf` output. + ++++ -Looks like `sd` is always `0` which will cause the logp to go to `-inf`. Of course, we should not have used a prior that has negative mass for `sd` but instead something like a `HalfNormal`. +### How to capture `Print` output for further analysis -We can also redirect the output to a string buffer and access the proposed values later on (thanks to [Lindley Lentati](https://github.com/LindleyLentati) for providing this example): +When we expect many rows of output from `Print`, it can be desirable to redirect the output to a string buffer and access the values later on (thanks to **Lindley Lentati** for inspiring this example). Here is a toy example using Python `print` function: ```{code-cell} ipython3 import sys from io import StringIO -x = np.random.randn(100) - old_stdout = sys.stdout mystdout = sys.stdout = StringIO() +for i in range(5): + print(f"Test values: {i}") + +output = mystdout.getvalue().split("\n") +sys.stdout = old_stdout # setting sys.stdout back +output +``` + +### Troubleshooting a toy PyMC model + +```{code-cell} ipython3 +rng = np.random.default_rng(RANDOM_SEED) +x = rng.normal(size=100) + with pm.Model() as model: + # priors mu = pm.Normal("mu", mu=0, sigma=1) sd = pm.Normal("sd", mu=0, sigma=1) - mu_print = tt.printing.Print("mu")(mu) - sd_print = tt.printing.Print("sd")(sd) + # setting out printing for mu and sd + mu_print = Print("mu")(mu) + sd_print = Print("sd")(sd) + # likelihood obs = pm.Normal("obs", mu=mu_print, sigma=sd_print, observed=x) +``` + +```{code-cell} ipython3 +pm.model_to_graphviz(model) +``` + +```{code-cell} ipython3 +with model: step = pm.Metropolis() - trace = pm.sample( - 5, step, tune=0, chains=1, progressbar=False - ) # Make sure not to draw too many samples + trace = pm.sample(5, step, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED) +``` -sys.stdout = old_stdout +Exception handling of PyMC v4 has improved, so now SamplingError exception prints out the intermediate values of `mu` and `sd` which led to likelihood of `-inf`. However, this technique of printing intermediate values with `aeasara.printing.Print` can be valuable in more complicated cases. -output = mystdout.getvalue().split("\n") -mulines = [s for s in output if "mu" in s] ++++ + +### Bringing it all together + +```{code-cell} ipython3 +rng = np.random.default_rng(RANDOM_SEED) +y = rng.normal(loc=5, size=20) + +old_stdout = sys.stdout +mystdout = sys.stdout = StringIO() + +with pm.Model() as model: + mu = pm.Normal("mu", mu=0, sigma=10) + a = pm.Normal("a", mu=0, sigma=10, initval=0.1) + b = pm.Normal("b", mu=0, sigma=10, initval=0.1) + sd_print = Print("Delta")(a / b) + obs = pm.Normal("obs", mu=mu, sigma=sd_print, observed=y) -muvals = [line.split()[-1] for line in mulines] -plt.plot(np.arange(0, len(muvals)), muvals) -plt.xlabel("proposal iteration") -plt.ylabel("mu value"); + # limiting number of samples and chains to simplify output + trace = pm.sample(draws=10, tune=0, chains=1, progressbar=False, random_seed=RANDOM_SEED) + +output = mystdout.getvalue() +sys.stdout = old_stdout # setting sys.stdout back +``` + +```{code-cell} ipython3 +output +``` + +Raw output is a bit messy and requires some cleanup and formatting to convert to `numpy.ndarray`. In the example below regex is used to clean up the output, and then it is evaluated with `eval` to give a list of floats. Code below also works with higher-dimensional outputs (in case you want to experiment with different models). + +```{code-cell} ipython3 +import re + +# output cleanup and conversion to numpy array +# this is code accepts more complicated inputs +pattern = re.compile("Delta __str__ = ") +output = re.sub(pattern, " ", output) +pattern = re.compile("\\s+") +output = re.sub(pattern, ",", output) +pattern = re.compile(r"\[,") +output = re.sub(pattern, "[", output) +output += "]" +output = "[" + output[1:] +output = eval(output) +output = np.array(output) ``` ```{code-cell} ipython3 -trace["mu"] +output ``` -Notice that for each iteration, 3 values were printed and recorded. The printed values are the original value (last sample), the proposed value and the accepted value. Plus the starting value in the very beginning, we recorded in total `1+3*5=16` value above. +Notice that we requested 5 draws, but got 34 sets of $a/b$ values. The reason is that for each iteration, all proposed values are printed (not just the accepted values). Negative values are clearly problematic. + +```{code-cell} ipython3 +output.shape +``` + +## Authors + +* Authored by Thomas Wiecki in July, 2016 +* Updated by Igor Kuvychko in August, 2022 ([pymc#406] (https://github.com/pymc-devs/pymc-examples/pull/406)) + ++++ + +## Watermark ```{code-cell} ipython3 %load_ext watermark -%watermark -n -u -v -iv -w +%watermark -n -u -v -iv -w -p aesara,xarray ``` + +:::{include} ../page_footer.md +:::