diff --git a/arviz/wrappers/__init__.py b/arviz/wrappers/__init__.py index d8993d14c1..2184a56ddc 100644 --- a/arviz/wrappers/__init__.py +++ b/arviz/wrappers/__init__.py @@ -1,5 +1,5 @@ """Sampling wrappers.""" from .base import SamplingWrapper -from .wrap_pystan import PyStanSamplingWrapper +from .wrap_stan import PyStan2SamplingWrapper, PyStanSamplingWrapper -__all__ = ["SamplingWrapper", "PyStanSamplingWrapper"] +__all__ = ["SamplingWrapper", "PyStan2SamplingWrapper", "PyStanSamplingWrapper"] diff --git a/arviz/wrappers/base.py b/arviz/wrappers/base.py index edefac64d9..3877c23e54 100644 --- a/arviz/wrappers/base.py +++ b/arviz/wrappers/base.py @@ -1,3 +1,4 @@ +# pylint: disable=too-many-instance-attributes,too-many-arguments """Base class for sampling wrappers.""" from xarray import apply_ufunc @@ -230,11 +231,9 @@ def check_implemented_methods(self, methods): if method in supported_methods_1arg: if self._check_method_is_implemented(method, 1): continue - else: - not_implemented.append(method) + not_implemented.append(method) elif method in supported_methods_2args: if self._check_method_is_implemented(method, 1, 1): continue - else: - not_implemented.append(method) + not_implemented.append(method) return not_implemented diff --git a/arviz/wrappers/wrap_pystan.py b/arviz/wrappers/wrap_stan.py similarity index 54% rename from arviz/wrappers/wrap_pystan.py rename to arviz/wrappers/wrap_stan.py index 20823266a7..8f59bb26e0 100644 --- a/arviz/wrappers/wrap_pystan.py +++ b/arviz/wrappers/wrap_stan.py @@ -1,11 +1,14 @@ # pylint: disable=arguments-differ """Base class for PyStan wrappers.""" -from ..data import from_pystan +from typing import Union + +from ..data import from_cmdstanpy, from_pystan from .base import SamplingWrapper -class PyStanSamplingWrapper(SamplingWrapper): - """PyStan sampling wrapper base class. +# pylint: disable=abstract-method +class StanSamplingWrapper(SamplingWrapper): + """Stan sampling wrapper base class. See the documentation on :class:`~arviz.SamplingWrapper` for a more detailed description. An example of ``PyStanSamplingWrapper`` usage can be found @@ -47,16 +50,65 @@ def sel_observations(self, idx): """ raise NotImplementedError("sel_observations must be implemented on a model basis") - def sample(self, modified_observed_data): - """Resample the PyStan model stored in self.model on modified_observed_data.""" - fit = self.model.sampling(data=modified_observed_data, **self.sample_kwargs) - return fit - def get_inference_data(self, fit): """Convert the fit object returned by ``self.sample`` to InferenceData.""" - idata = from_pystan(posterior=fit, **self.idata_kwargs) + if fit.__class__.__name__ == "CmdStanMCMC": + idata = from_cmdstanpy(posterior=fit, **self.idata_kwargs) + else: + idata = from_pystan(posterior=fit, **self.idata_kwargs) return idata def log_likelihood__i(self, excluded_obs_log_like, idata__i): """Retrieve the log likelihood of the excluded observations from ``idata__i``.""" return idata__i.log_likelihood[excluded_obs_log_like] + + +class PyStan2SamplingWrapper(StanSamplingWrapper): + """PyStan (2.x) sampling wrapper base class. + + See the documentation on :class:`~arviz.SamplingWrapper` for a more detailed + description. An example of ``PyStanSamplingWrapper`` usage can be found + in the :ref:`pystan_refitting` notebook. For usage examples of other wrappers + see the user guide pages on :ref:`wrapper_guide`. + + Warnings + -------- + Sampling wrappers are an experimental feature in a very early stage. Please use them + with caution. + + See Also + -------- + SamplingWrapper + """ + + def sample(self, modified_observed_data): + """Resample the PyStan model stored in self.model on modified_observed_data.""" + fit = self.model.sampling(data=modified_observed_data, **self.sample_kwargs) + return fit + + +class PyStanSamplingWrapper(StanSamplingWrapper): + """PyStan (3.0+) sampling wrapper base class. + + See the documentation on :class:`~arviz.SamplingWrapper` for a more detailed + description. An example of ``PyStan3SamplingWrapper`` usage can be found + in the :ref:`pystan3_refitting` notebook. + + Warnings + -------- + Sampling wrappers are an experimental feature in a very early stage. Please use them + with caution. + """ + + def sample(self, modified_observed_data): + """Rebuild and resample the PyStan model on modified_observed_data.""" + import stan # pylint: disable=import-error,import-outside-toplevel + + self.model: Union[str, stan.Model] + if isinstance(self.model, str): + program_code = self.model + else: + program_code = self.model.program_code + self.model = stan.build(program_code, data=modified_observed_data) + fit = self.model.sample(**self.sample_kwargs) + return fit diff --git a/doc/source/api/wrappers.rst b/doc/source/api/wrappers.rst index ef16cc6987..5c66b7a096 100644 --- a/doc/source/api/wrappers.rst +++ b/doc/source/api/wrappers.rst @@ -11,3 +11,4 @@ Experimental feature SamplingWrapper PyStanSamplingWrapper + PyStan2SamplingWrapper diff --git a/doc/source/user_guide/pystan2_refitting.ipynb b/doc/source/user_guide/pystan2_refitting.ipynb new file mode 100644 index 0000000000..de2c82cd3c --- /dev/null +++ b/doc/source/user_guide/pystan2_refitting.ipynb @@ -0,0 +1,460 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "(pystan2_refitting)=\n", + "# Refitting PyStan (2.x) models with ArviZ\n", + "\n", + "ArviZ is backend agnostic and therefore does not sample directly. In order to take advantage of algorithms that require refitting models several times, ArviZ uses {class}`~arviz.SamplingWrapper` to convert the API of the sampling backend to a common set of functions. Hence, functions like Leave Future Out Cross Validation can be used in ArviZ independently of the sampling backend used." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Below there is one example of `SamplingWrapper` usage for PyStan exteding {class}`arviz.PyStan2SamplingWrapper` which already implements some default methods targetted to PyStan.\n", + "\n", + "Before starting, it is important to note that PyStan cannot call the C++ functions it uses. Therefore, the **code** of the model must be slightly modified in order to be compatible with the cross validation refitting functions." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import arviz as az\n", + "import pystan\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the example we will use a linear regression." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(26)\n", + "\n", + "xdata = np.linspace(0, 50, 100)\n", + "b0, b1, sigma = -2, 1, 3\n", + "ydata = np.random.normal(loc=b1 * xdata + b0, scale=sigma)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xc5Zno8d87oy6N+kiyiq3qblwxNrYxzQRILi2QQAokkJDsZnfZLNkQNjf3brKbTbmbsiWbhACJEwglgQRIQjGmGnDvcpVk9d41ozKamff+cWZGGhWrzaiMnu/n448055yZc84H8ejVc573eZXWGiGEEKHJNNMXIIQQIngkyAshRAiTIC+EECFMgrwQQoQwCfJCCBHCwmb6AgZLTU3Vubm5M30ZQggxpxw6dKhZa20dad+sCvK5ubkcPHhwpi9DCCHmFKVUxWj7JF0jhBAhTIK8EEKEMAnyQggRwiTICyFECJMgL4QQIUyCvBBChDAJ8kIIEcIkyAshxBTZ+5w8c6ASt3v2tW6XIC+EEFP00rFaHnruBHtKmmf6UoaRIC+EEFN0rsEGwOunG/y2t9j6uPvx/ZQ322fisgAJ8kIIMWXnG7sAeP1UA4NX2/v9oWreOdfE0weqZurSAhPklVLlSqkTSqmjSqmDnm3JSqldSqnznq9JgTiXEELMNucbbMRFhlHb0UtxbScAWmueO1wNwMsn65ippVYDOZK/Smu9Rmu9wfP6a8BurXURsNvzWgghQkpHTz/1nb3ctTEHpWDXKSNlU1zbybkGG2tyEqlo6eZUXeeMXF8w0zU3Azs93+8EbgniuYQQYkaUNBr5+E35KaxfmOQL8s8dribCbOJHH1+D2aR4+UT9jFxfoIK8Bl5TSh1SSt3v2Zauta7zfF8PpI/0RqXU/Uqpg0qpg01NTQG6HCGEmB7nG4x8fFGahR3L0zlV10lFi50Xj9Zy7fI08lJj2Zyfwl9OzEzKJlBBfqvWeh1wA/AlpdQVg3dq485GvDut9SNa6w1a6w1W64g974UQYtY632gjKtxEdlI01y43xrLfeKGYFruD29ZmA3DDqgzKmu2+KpzpFJAgr7Wu8XxtBP4AbAQalFILADxfGwNxLiGEmE3ONXRRmBaHyaQosMaRb43lnXNNpMRGsH2JMXC9bnkGJgV/OVE3xqcF3pSDvFIqVill8X4PXAecBF4E7vEcdg/wwlTPJYQQs01Jo43FaRbf6x3LjNH8TWsyCTcbIdZqiWRjXjIvn5yDQR4j175HKXUM2A/8WWv9CvBdYIdS6jxwree1EELMWXvLWvj2n0/52hd09vZT19FLYXqc75hb1mZhtURy18aFfu+9cdUCzjXY+O7LZ/jeK2f41z+doqIl+JOkprzGq9a6DFg9wvYW4Jqpfr4QQswG9R29/NUTh2jr7ufKJWlsKUz1VdYMHskvWxDPga9fO+z916/M4PuvnOVnb5cSYTbR73bT1t3PDz42LHwG1KxayFsIIWYjl1vz5WeO0ud0Y4kK46n9lWwpTPVV1ixOt4zxCZBmieLI/9mBAsLMJh76/XFeOl7Lv9yygpiI4IViaWsghBBj+NnbpXxQ1sI3b1rB7euzebW4nhZbH+cbBiprxiPcbCLMk6f/6Ppsuh0uXjkZ3Pp5CfJCCHERhyvb+OGuc9y0OpPb12dz18aF9LuMlgXnGm2+ypqJ2rAoiZzkaJ4/XBOEqx4gQV4IIS7iNx9UkBAdzrdvXYlSisXpFtYvSuLp/VWcb+iiKG3sVM1ITCbFbWuzea+0mdr2ngBf9aDzBO2ThRAiBByrbmf9oiQsUeG+bXdtXEhZs526jl6KBlXWTNRt67LQGv54NHijeQnyQggxis7efsqa7KzOTvDb/uFVC7BEGQ9LJzuSB1iUEsuluUk8d6g6aC0PJMgLIcQoTlZ3AHBJdqLf9ugIM7euzQJg8RRG8gC3rcumtMnOcc+5Ak2CvBBCjOKYL8gnDNv399cu5v/dfgmLUmKndI4bVy0gIszE857e84EmdfJCCDGK49XtLEqJITEmYti+5NgI7tiQM+VzJESH84u7NwxLCQWKBHkhhBjF8eoO1i0K/qJ22xcHrwOvpGuEEGIETV191LT3BG2EPV0kyAsh5q0z9Z0cq2ofcd/xamP70Ieuc42ka4QQ89aDzx6juLaTj1yygIdvXEZW4kB7gmPVHZgUrMyKn8ErnDoZyQsh5iWtNRUt3eRbY9l1qoFrfvAWv3rvgm//8ep2itIsQW0eNh0kyAsh5qWOnn5sfU4+sXEhux/czub8FP75pVO8fqoBrTXHqztYNcfz8SBBXggxT1W3Gf1ispNiyE6K4aefWs/KrHi+/OxRPihtodXumPMPXUGCvBBinqpu6wbwtQmOCjfz00+uRwGf//VBYO4/dAUJ8kKIEPLC0Ro6uvvHdax3JJ+TFOPblpMcw4/vXIPd4SLcrFi6YPJ9aWYLCfJCiJBQ1drNA08f5akDleM+3hIZRny0/4PVq5em842PLOfuzblEhpmDcanTam4/NhZChCxbn5Pdpxu4aXUmSo29KEdJk7Heqnfd1bFUt/WQlRQ94mfftzVvYhc7i8lIXggxKz21r5IHnj7K+XEG7dLGiQf5nOSYsQ+c4yTICyFmpSNVbQCcrusc1/FlzXYASptsw3qzP7bnAocr23yvtdZUt3WPe23WuUyCvBBiVjpaabQVOFvfNa7jyzzpmq5eJ01dfb7t9j4n//rnUzz6bplvW3t3P3aHi+wkGckLIcS0a+zspbajFxh/kC9tspOZEAUM5OcBTtV1ojUcqxpYlKNqSPlkKJMgL4SYdY54moblJEdzZhxBvqu3n6auPq5bkQEM5OcBTtYYwb2mvcc3wh+YCCVBftyUUmal1BGl1J88r/OUUvuUUiVKqWeUUsO77gshxAiOVrUTZlLcujabmvYeunovXvte1mTk4zcXpBAbYabU8xrgRM3ACN7bWXJgIpSkaybiAeD0oNffA36ktS4E2oD7AnguIUQIO1rZzrIF8b62AucaLj6aL/WkZwqscRSkxflV2BTXdHJZXjImNbCcX3VbD/FRYSREhwfpDmaPgAR5pVQ28GHgUc9rBVwN/N5zyE7glkCcSwgxu93/64N875Uzk36/y605Xt3OmpxElmQYM07HStmUNdkxmxQLk2MotMb5gn6Pw8X5xi4uy0tmcbrF1zu+uq1nXoziIXAj+R8DXwXcntcpQLvW2ul5XQ1kjfRGpdT9SqmDSqmDTU1NAbocIcRMcLrcvHW2iQ9KWyb9GSWNNuwOF2tyEslKjMYSGTbmw9eyZhsLk2OICDNRkBZHXUcvtj4np+s7cWtYkZXA6uxEjle3o7WmqnV+lE9CAIK8UuojQKPW+tBk3q+1fkRrvUFrvcFqDd46h0KI4LvQbMfhcvty3pNx1FMfv2ZhIkopFmdYxhzJlzbaKbDGAkbKBoySymJPPn5lVgKX5CTQ1t1PVWuPjOQnaAtwk1KqHHgaI03zH0CiUsrbNiEbqAnAuYQQs9gpz8SlZpuDbodzjKNHdrSqnYTocPJSjKC9JMPC2fquYROcvFxuzYUWO/me4F6YZryvpNHGiZoOkmMjyEyIYrWno+SbZxvp6XeRkywj+XHRWj+stc7WWucCdwJvaK0/CbwJ3O457B7ghameSwgxuw0ecXvLFCfqSGU7q3MSMZmMnjJLMyx09PTT0Nk34vG17T04nG7yU43gviglljCTorTJxsmaTlZkxqOUYkmGhcgwE38+XgfMj8oaCG6d/EPAPyilSjBy9I8F8VxCiABxutz87O1S7H0TH4mfruskzBOcq1onnrKx9zk519DFmpyBPu5L0r0PX0dub+CrrEkzRvLhZhMLU2I4VdvJuYYuVmUl+LavyIznQEUrMD9q5CHAQV5r/ZbW+iOe78u01hu11oVa6zu01iP/GhZCzCoHytv47stneLW4fsLvPVPXxeaCFODiQd7hdNM5Qu37iZoO3BrWDgrySzOMhbRHe/jqrYn3juQBCq1x7ClpxunWrMwaWN1pdU4i3qxPlgR5IcR85J3yf65hfN0cvdrsDuo7e9lamEpUuImqi6Rr/ukPJ9jynTc4MqhpmMPp5tF3yzApIxh7JcSEkxEfNWqQL2uykRAdTnLswHzLwrQ4+l1GNF81OMh78vIJ0eHER4V+jTxIkBdCDOHNpZc0jq9njNdpTzpl2YJ4spNiRh3JN3b28scjNdgcTu5+bD9Hq9pxON38zW8P8/rpRv75phV+ARuMh6+nRw3ydvKtsX594b0VNgnR4X5pGe8vj/ny0BUkyAshhqie5Ej+TJ0RhJctiCcnKXrUkfxT+6twujVP3ncZSbERfPrRfXz2V/t57VQD37p5BXdvzh32nqUZFkobbfS73MP2lTbZfEHdq9CTn1+ZFe8X/HNTYozAnzg/HrqCBHkhxBDVrUZwrmrrpsfhGvf7Ttd1khoXgdUSSU5yDNWt3cPKHvtdbp7cV8H2xVYuL0zl6fs3kRQbwXslLaMGeDB+cThcbv7p+RNUtBg5eK01J6o7aOzqI98a63d8vjUWk4JVWf4LcSul+O9PrOXLOxaP+77mOln+T4gQca6hi+auPi4vTJ3S51S3GWufdvU5KW2y+T24vJgz9V2+h6Q5STF09Tnp7HGSEDOQ+361uJ7Grj6++9FFAGQmRvP8X1/OhWY7l+Ymj/rZN6zK4FDFIp45WMVzh6vZVmSlpNFGTXsPJgXrFib5HW+JCueJ+y5j2YL4YZ+1rWh+TbqUkbwQIeIHr53lq88dn9JnOJxu6jt72bbY+EVxfpx5eafLzbmGLpYtMModvTnvqiEzX3/9fgULk2PYvjjNty01LvKiAR4gMszMv9yykj1fvYrPb8unpNHGisx4vvfRVex9+Bo25acMe8/lhakkxUrzWxnJCxEiatp7aLb1obUe18LXI6nr6MGtYWuhlV2nGsadly9vsdPndPtG8t6JRlWt3b6/BE7VdrK/vJWv37gMs2ly15cWH8XDNy7j4RuXTer985GM5IUIEbXtvfT2u+meQB59KG9lTW5qDHmpsZwfZ5A/PeihK+BbIHvwSP6JfRVEhZu4Y0P2pK9PTJwEeSFCQI/DRavdAUCLzTHpz/FW1uQkxVCUZhl3uuZMvTHTtcDTN8aoQw+jyvMQ1+XWvHqynh3LM0iMkRTKdJIgL0QIqO0YKFdssY9vcvmxqnZ2/PBt3y8HMEbyZpNiQUIURelxVLZ209s/9l8Gp+u6KEyLIzLM7NuWkxzjG8kfLG+lxe7gQyvSx3tLIkAkyAsRAmrbBwX5cY7kXy2u53yjjb1lA73fq1q7yYiPIsxsoijNgtb4rbI0kiOVbRy40DqskiVn0ISoV4sbiAgzceWStJE+QgSRBHkhQsDgID94ZH4xRyqNVZIOVwy0FjD6rBuVMYvTjQlFFwvyLxyt4eOP7CUpNoK/u6bIb19OcjTVbT243ZpXi+vZWphKXKTUekw3CfJChICa9l7f983jSNe43JpjnkWtD1f6B3nvQ9PcVKNl72jrq/5o1zkeePooa3MS+eOXtpCX6j8hKSc5hj6nm7fPN1HT3iOpmhkiv1aFCAG17T1kxEfR2ds/rnTNuYYuuh0uMhOiOFnTSZ/TyLs3dPX6RvLhZpNRYTPCSP6NMw38x+7z3L4+m3+7dRURYcPHizmeMsrH3r2AScG1yyTIzwQZyQsRAuo6eshMjCIlLmJc6RpvquYzW3JxuNycrOmktr0XrQeCM8DidAvnh4zku3r7+fofTrIk3TJqgIeBCVF7SprZkJtMSlzkZG9PTIEEeSFCQG17L5mJ0STHRtJsGztdc6SyjeTYCG5ZmwUYeXlv+eTgro2FacMrbL7/ylnqO3v57kdHD/DG5wz8svjQiowJ35MIDAnyQsxxWmtq2nvISowmNTZiXOmaI1XtrM1JJM0SRU5yNIcr23w17dnJ/iN5tx5YfWn/hVZ+s7eCz16ex9oh/WKGigo3Y7UYo/frlkuqZqZITl6IOa7F7sDhdJOZGE2r3UFx7cjL5Hl19PRT0mjjljWZgNHca29Zi+9Ba0Z8lO/YJRlGhc0tP3mPNEsUdoeT7KRovvKh8XVxLLTGsSAhyvcwV0w/CfJCzHHe8snMxGjqOnppsV+8f82xKiMf7x2Jr1+UxAtHa9l/oZXMxGi/vjIF1jh+cMdqSppsNHT00mx38DdXFRITMb7Q8aOPr5nKrYkAkCAvxBw3EOSjqGiJoN+l6epzjrq83ZHKdpSCS7KNxmHeNr2HKtq4vMC/m6NSio+un3yvmYyEqLEPEkElOXkh5jhvjXxWYrRv2byL5eWPVLWxOM2CxfNLYGmGhehwox3B4MoaERokyAsxx9W29xATYSYhOtxXptgySoWN1pojle2sXTiwYlKY2cTqHGNUP7iyRoQGCfJCzHG17T1kJkajlCLFO5IfpVb+QrOdjp5+vyAPAymb7Hm0wPV8IUFeiDnOG+QBUuIunq45UN4KMKz88fICYyWoojRLsC5TzJApB3mlVJRSar9S6phSqlgp9U3P9jyl1D6lVIlS6hmllDSRFiIIatp7yUo0HnB6c/KtI/Sv0Vqz8/0K8q2xFFrj/PZtLUpl94Pbx72eq5g7AjGS7wOu1lqvBtYA1yulNgHfA36ktS4E2oD7AnAuIcQgvf0umm19LEgwRvKRYWYskWE0jzCSf+d8M6fqOvniFQWYRlh+r2BI4BehYcpBXhu8HYzCPf80cDXwe8/2ncAtUz2XEMJffYdRWeNN14CRshkpJ//Tt0rIiI/i5rWZ03Z9YuYFJCevlDIrpY4CjcAuoBRo11o7PYdUA1mjvPd+pdRBpdTBpqamQFyOEPPG4Bp5r+TYiGHpmiOVbewta+Vz2/L8Vm8SoS8gQV5r7dJarwGygY3A0gm89xGt9Qat9Qar1RqIyxFi3qjxBPksv5F85LAHrz97u5SE6HDu3LhwWq9PzLyAVtdorduBN4HNQKJSyjujNhuoCeS5hBBG90nwn1maOiRdU9LYxavFDdyzeZGszDQPBaK6xqqUSvR8Hw3sAE5jBPvbPYfdA7ww1XMJIfzVdfRgtUT6pWCMdI0Dt1sD8MyBKiLMJu65PHeGrlLMpECM5BcAbyqljgMHgF1a6z8BDwH/oJQqAVKAxwJwLiHEIGVN9mGzVFNiI3G5NR09/QC8c66ZS/OSZNGOeWrKf7tprY8Da0fYXoaRnxdCBIGtz8nhyjY+f0W+33bfhCi7A4fLzdmGLm5bN+7HZCLESIJOiFlk5/vlhJkVVy1J8yuLHMn7Jc043Zrti/0LFlJiB/rXHPW0Fd5WJEUN85UEeSGmQV1HD9/4YzFf3lHEisyRZ5V2dPfzf18s9r1emmHhmzet4LL8lBGPf/tcE3GRYb6+M16DR/Lvnm8iNS6SpRnSrmC+kt41QkyDX71XzuunG7j3Vweo6+gZ8ZgqzxqrD9+wlIdvWEpHTz/ffOkUWuthx2qteftcE5cXpAxbZ9XbpKzZ1see881sK0odcYarmB8kyAsRZH1OF787VM2anETsfS4++8sDdPX2Dzuuus0I/lsKU/nC9gL+9uoiTtV1cqiibdixZc12qtt62L5keBomyRPk3z3fTIvdwbai1ADfkZhLJMgLEWSvFjfQanfw5R2L+emn1lHSaONLvz1Cv8vtd9zQiU23rM0kPiqMX71fPuwz3z5rzA6/YoRce7jZREJ0OG+dbQRga6EE+flMgrwQQfbbfRXkJEezrTCVbUVWvnXzSt4518Trpxr8jqtu6yY2wkxijLFiU0xEGB/bkMMrJ+tp6Oz1O/btc03kW2NHXSA7Jc5YBnBphoW0eFmCbz6TIC9EEJU22dhb1sqdly705cVvW5eFUnCuweZ3bHVbD9lJMX4LcH968yJcWvPkvkrftt5+F3vLWoZV1QyW6qmwueIix4j5QYK8EEH01L5KwkyKOzYMLIYdFW4mKzGaC83+Qb6mrYesIRObFqXEctWSNH67rxKH00jv7LvQSp/TfdEg7+0rL/l4IUFeiCDp7Xfx+8PVXLcinTSLf8okLzWWsma737bqtu4R11i95/Jcmm19fOfl0zx3qJqn9lUSGWZi0yillWD0sokON3NpbnJgbkbMWVInL8QU9bvc1Hf0DsuPv3Gmkfbufu4aofNjfmoszx+uQWuNUorO3n46e51+3SS9thWmsmxBPL98r9y37dplaUSFj94y+EtXFfLRddkXPUbMDxLkhZii/36jhJ++Xcp7D12N1TLQH+b1Uw0kxoSzeYQRd15qLF19TpptDqyWSGo85ZPZScMfpJpMipf+ZgvtPf1097mw9TlZmDLyA1cvqyXS71rE/CXpGiGmwOXWPHOgCofTzWun6v22v3m2kauWpBFmHv6/Wb5nqb0LnpRNtS/Ij9zKIMxsIjUukoUpMSzPjJeWwWLcJMgLMQXvnm+ivrOXMJPilZMDQf5oVRtt3f1cvTRtxPflpcYCUNZkPHyt8cx2HfrgVYipkiAvxBT87mA1ybERfHZLLu+XttDmWaxj9+lGzCY1agljZmI0EWEmv5F8VLjJ15JAiECRIC/EJLXZHew61cDNazK5eU0WLrdml2eC0+7TjVyam0RCdPiI7zWbFLkpMb4Km+q2HrISo/1q5IUIBAnyQkzSC0drcLjc3LE+hxWZ8eQkR/OXk3VUtXZztqGLa5elX/T9eamxvpF8TXvPiA9dhZgqCfJCTNKzB6tZlZXA8sx4lFLcuHIB75U088JRYznj0fLxXvnWOCpa7Dhd7lFr5IWYKgnyQkzCyZoOTtV1+s1kvX5lBv0uzU/eLCUvNdZXQTOavNRY+l2acw022rr75aGrCAqpwxJiiN5+F8CwiUT9LjdvnW3ixWO1vH6qgahwEzetzvTtX5OTSGZCFLUdvWOO4sGYEAWwp8ToKCnpGhEMEuSFGOKvnjhEn9PNbz+/yW/7/3mhmKf2V5IUE86t67L4xMaFJMYMVMMopbh+5QIef+8C1ywbO8h7yyjfPd8MMOJsVyGmSoK8CGn9Ljdf/f1xPnN5LqtzEsc83u3W7L/Qit3h4mRNByuzjKX6mrr6eO5QNXesz+bfbltF+AgTnAA+ty2P2EgzG8fRMyY5NoKE6HD2XWgFIEfSNSIIJCcvQtqBC6384UgNvx3Uqvdiylvs2B1GumbwYh1P7a/E4XLzxSsLRg3wYNS/P3jdkhFnuQ6llCIvNRaH001EmDGjVYhAkyAvQtprnrr1PSXNw9ZKffj547x4rNZvW3FtJwAbFiXx4rFaWmx9OJxunthbwfbFVgrGeJg6Ud68fFZitKzDKoJCgrwIWVprXj/dQITZRE17j68mHYyeMU/tr+KJvRV+7zlZ20GE2cS3bl6Jw+nm6QNVvHyyjsauPj6zJTfg1+jNy0v5pAiWKQd5pVSOUupNpdQppVSxUuoBz/ZkpdQupdR5z9ekqV+uEON3tqGL6rYePrs1FzBG814vn6wD4GhVu6+aBqC4ppMlGRaWZ8aztTCVJ/ZW8Ph75eSlxrJ9hPVUp8pbZikPXUWwBGIk7wQe1FovBzYBX1JKLQe+BuzWWhcBuz2vhZg23jVU79uSR3ZStK+KBeCVk/VEhJlwON0crWoHjJF/cW0HK7PiAWOxjrqOXo5VtXP35kVBSafISF4E25SDvNa6Tmt92PN9F3AayAJuBnZ6DtsJ3DLVc4nQsP9CK88dqsbpck/5s4bm2QfbdbqR1TmJpMVHsa0olb2lLb7ZpcerO7h3Sx5Kwb4yo7qltqOXtu5+lmcaFTVXL00jJzma2Agzt6/PHvU8U1GYFsf/Wp3JNWO0QBBisgKak1dK5QJrgX1Auta6zrOrHpCfYgHAv796lgd/d4yP/NcePihtmfTnPHOgks3feYPTdZ3D9jV2GiPwHZ569a2FVrr6nByrbve1BL7z0hyWZsSzv9y4hpM1HQCszDRG8maT4scfX8tPPrkOS9TIjcamKiLMxH/dtZZlC+KD8vlCBCzIK6XigOeAv9da+/1fp43h1ohDLqXU/Uqpg0qpg01NTYG6HDGLNdv6WJweR1evk7t+sZd/frF4wp/x2J4LPPTcCeo7e/nP3eeH7d99phGAHcszALi8IAWljIlHr5ysZ2mGhdzUWC7LS+ZQRRsOp5vi2k5MCpZmDATc9YuSuHLJ2BObhJitAhLklVLhGAH+Sa31857NDUqpBZ79C4DGkd6rtX5Ea71Ba73Bag38gy0x+zTb+ticn8LuB7dz+/psdn5QTlVr97jeq7Xmv3af51/+dIobVmbwhe35vHyynnMNXX7HvX6qgZzkaBanGw82k2IjWJWVwIvHajlU2cYNKxcAsCk/md5+Nydq2imu6aAwLY7oCFkXVYSOQFTXKOAx4LTW+oeDdr0I3OP5/h7ghameS8x9Dqebzl4nKXGRRIWb+fKOxSjg2YNV43r/66cb+cGuc9y2Lov/umstX7yigJgIMz95s8R3TLOtjz0lzVy7LN2vP/vWwlTKmuxoDTesMkb4G/OM9Vf3lrVysraDlZ58vBChIhAj+S3Ap4GrlVJHPf9uBL4L7FBKnQeu9bwW81yrZ+WklDij50tWYjTbF1t59mDVuB7EvnGmAUtUGN//6CWEmU0kxUbw6U2LeOlYLRea7TR29XLXI3tRCu5Yn+P33q1FqQDkW2MpSjNG+MmxESxOj+MvJ+po6OxjeabkxkVoCUR1zR6ttdJaX6K1XuP59xetdYvW+hqtdZHW+lqtdWsgLljMbc22PgBSYgem8N+1cSENnX28eXbsZzIflLZwWV6yX9uA+7blEW428Z2/nObOR/ZS097DLz+zcVjAXr8oidS4CD66LttvhH9ZXopvpqu3V40QoUJmvIpp1eIZyafGDXRvvHppGmmWSJ7af/H+MrXtPZS3dLO5INVve5olirs2LuS1Uw00dPSy896NbC5IGfb+yDAz73z1Kv5qe4Hf9svyB5qJyUhehBrpQimmVYt3JD+oGVeY2cTHNuTwP2+VUNveQ+Yosz+95ZaXjxDA//rKAmrbe/jC9gLWLxp9cnVMxPAf+Y15RpBflBJDfJBKJYWYKTKSF9Oqxeafk/f6+KU5aC7+AMNXefAAABe3SURBVPb90haSYsJZkm4Zti8tPopH7t5w0QA/mjRLFCsy48fVHliIuUZG8mJaNdv7iDCbsET6/+jlJMewrcjKbz6o4COXZFKY5t/tUWvN3rIWNhekBKW9wNP3b7poC2Eh5ir5qRbTqsXmICUuwu/Bp9c3PrwMpRQf//kHFNd2+O2rbO2mpr2HzfnDUzWBYIkKH7bcnxChQIK8mFYttr5hqRqvonQLz35hE5FhJu56ZC+HK9t8+7z5+KEPXYUQFydBXkxYb7+L77x8etyzVAdrsTv8yieHyrfG8ewXN5McG8GnHt3HvjIjuL9f2oLVEkmBNXbS1y3EfCRBXkyI1pqv/+EkP3+7jJeO1479hiG86ZqLyU6K4dkvbGZBQhSf+eUBPiht4YOyFk//GVk9SYiJkCAvJuTJfZU8d7gagLIm+xhH+9Na02LvG9dapmnxUTx1/yayk6K5+/F9NHX1BS0fL0QokyAvxu1wZRvffKmYK5dY2ZibTFmTbULv73a46O13kxJ78ZG8V5rFCPT5qUalzeWSjxdiwiTIi3Hp6Onnr584TEZCFD/++BoK0uL81kwdj4Ea+bFH8l6pcZE8+4XN/O6Lm1mYEjOh8wkhJMjPSv0uN3c/vp+D5bOn3c+T+yqM3u13riUxJoICayxt3f20edoUjEez3TvbdXwjea+EmHAulYlKQkyKBPlZqL6jl3fONfH66RFb8E/aa8X1PPz8iYsumTeSPqeLX71XztbCVNYuNGaU5nuqXMqax5+y8Y7kUy9SXSOECCwJ8rNQk6e/y4UJBNCx9Pa7+MYLJ3lqfyV7yyb2F8KLR2tp7Orj81fk+7Z58+SlE3j4OtC3ZmIjeSHE5EmQn4W8I96Rct6//qCcQxVtw7aP5dcflNPQ2Ud0uJlfvndh3O/TWvPouxdYkm7hiqKBB5/ZSdGEm9WEKmy8HSiTx/ngVQgxdRLkZyHviLe8pRu3eyC14nC6+dZLpyYUpAE6e/v5n7dK2b7Yyr1bc9l1umHcE5neOd/M2YYuPn9Fvl+NepjZxMLkmIv+tfFqcT1n6geW+2229REXGSbtA4SYRhLkZyHvwhoOp5vajh7f9rJmG0635nzDxNI4j757gfbufv7xQ0v49KZczEqx8/3ycb33F++UkR4fyU2rM4fty7fGjTqSf6+kmS8+cYjvv3LWt208E6GEEIElQX4WarYNVKwMTtmcrTcWqy5rto1rqTww/ip47N0yPrxqASuzEshIiOKGVQt45mAV9j4nTpebH752lqv//S1fCwGvF4/Vsqekmc9cnkdE2PAflXxrLBUt3bjc/g9yW+0O/uHZo2gNRyrbfA96W+x9466RF0IEhgT5Wcib1gD/IH/GE+T7XZqKcaZbHn/vAj39Lv7husW+bZ/dkktXr5P/eauEu36xl/98o4TWbgefemwfvztYhdaaH79+jr976ggbFiXx6c2LRvzsgtQ4HC431W0D16K15qHnjtNqd/DpTYto6+733YMxkpfKGiGmkwT5WajF5mBJhoWYCLNfOuRcfZdvRH2+oWtcn3W8uoNVWQkUWAf6s69bmMTqnER+8mYpp2o7+dHHV/P2V67isrwU/vH3x/nwf+7hx6+f57Z1WTz5+ct8v3CG8pVRDrrGJ/dVsutUAw9dv5RPbTJ+ORyubAeMv1BSJV0jxLSSID8Ltdj7sMZFkpsSS3mL/0j+iiIrwLjz8uUtdnJTh3du/Nr1S7lueTp/+rtt3Lo2m4SYcH752Uv51KaFnK7v5Gs3LOUHd6wmMmz0h6T5nl8cZZ6Renu3g3/7y2m2FaVy75Y8itLisESGcbiyDbdb02rvu2gHSiFE4MnKULNQs83BpbkRmM2KkzXG4hldvf3UtPfwicsWcqa+k3ONYwd5h9NNTVsPt67NHrZvc0HKsMWuw80m/vWWVTx0/VIs41jrNCkmnITocF8Pmyf3VdLtcPHwDct8qzetWZjI4Yo22nv6cWupkRdiuslIfpZxuty0dRu56/zUWKpau3E43ZzzjNyXpFsoSosbV7qmuq0bt4bcCfZ8GU+AB1BKkW+NpazJjsPpZuf7xqzY5ZnxvmPWLUziXEOX7y8SyckLMb0kyM8yrd0OtAZrXAR5qbG4NVS1dfsqa5ZkWFicbqGs2T5mhU1Fi/FAdFEQG3vlp8ZR1mzjpWPGrNjPbcvz279uURJuDW94WjSkSnWNENNKgvwsM7hTozeXfqHJzrmGLmIjzGQlRlOYFofD6aZyjAob7+h5UUrwVlPKt8bS0NnHT94qoSgtju2LrX771+QkArDrVAMgI3khpltAgrxS6nGlVKNS6uSgbclKqV1KqfOer0mBOFeo8zXx8qRrwCijPFPfyeIMCyaTYnG6BYDzY+TlK1q6iYsMC2ptesGgCpvPbcsbtnJTQnQ4RWlxnPWklyQnL8T0CtRI/lfA9UO2fQ3YrbUuAnZ7XosxNA9q4pUYE0FSTDhlzXbO1nexxBPcC9KMqpax8vLlLXYWpcQEdck8b4VNalwEN6/JGvGYdZ7OlUpBUowEeSGmU0CCvNb6HWBoa8ObgZ2e73cCtwTiXKHOG+S97XjzUmM5UN5KW3c/SzKMIB8XGUZWYvS4RvK5QUzVgJHvT4wJ5/Pb8kftSbN+kRHkk2MiMJtkjVYhplMwc/LpWus6z/f1QPpIByml7ldKHVRKHWxqagri5cwNzTYH4WZFfLRR3ZqbGktJ40BljVdRetxFa+WdLjdVrd1BfegKEBlm5oOvXcP9g9oQD7VukZGXl1SNENNvWh68aqN5yYgrVWitH9Fab9Bab7BarSMdMq+02IwJQ94US/6giUzekTxAUVocpU02XG6N1ppvvXSKn75V6ttf296L062DPpIHiI4wXzQllJ8aR3xUmEyEEmIGBHMyVINSaoHWuk4ptQAI7DJHIarF7iDVMjDizUv15rwj/SpTitIs9DmN0fqfT9Tx+HsXSIoJ5/4r8jGb1KDKmplfF9VkUnx5x2LpIy/EDAjmSP5F4B7P9/cALwTxXCGj2eY/9T/PM5JfkhHnd1xRuvH6Z2+X8u+vnWVhcgxt3f2c8MyQ9TYwG6mlwUz47Ja8UR/MCiGCJ1AllE8BHwBLlFLVSqn7gO8CO5RS54FrPa/FGFpsDlIHjdhzU2NQCpZmxPsdV+ipsHn6QBUrMxN45gubUAreOmv8wVTRbCcq3ESaRVIkQsxnAUnXaK3vGmXXNYH4/PlCa02Trc+vU2NMRBiP3r2BlVkJfsdaosLJSozG4XLzi7s3kJEQxersRN4628TfX7uYck9lTTDLJ4UQs580KJtFbH1OHE73sCqUa5aNWJjETz+1jviocDISogC4comV/9h9nla7g4oWu68VsBBi/pK2BrPI4Nmu43FJdqJfzv3KJWloDW+fa6SiNfg18kKI2U+C/CwyMNt1cnn0S7ISSI6N4On9VTic7qD2rBFCzA0S5IOkze7glZN1Yx84SLNvJD+5UkOTSbGtKJV9F4zJxxNtMSyECD0S5IPkqQOVfPGJw9S094z7Pb6WBlPo1HjlkoEJZYtmSfmkEGLmSJAPkipPnbp3Zafx8ObkpzJp6IoiK0pBhNlERnzUpD9HCBEaJMgHSXWbMYIvru0c93ta7H0kxoQTbp78f5aUuEguyU5kUUqMNAMTQkgJZbB4R/LFExjJG7Ndpz71/99vv4SefteUP0cIMfdJkA8Ct1v7cvETGck3D5ntOllFg7pVCiHmN0nXBEFjVx/9Lk1eaiz1nb2+B6pjabb1BSTICyGElwT5IKhuM1I1H1qRAYx/NN9ic0jPdSFEQEmQnwSjPf7ovA9dr1thtCMYT4WNw+mmo6dfRvJCiICSID8J33n5DNf96G3fiH0o7/blC+JZmBxDce3YQb5p0NquQggRKBLkJ+GNM42ca7Dx8Z/vpbzZPmx/VWsPqXGRRIWbWZEZP650ze8OVgEDi14LIUQgSJCfoG6Hk9ImGzeuyqCn38XHfv4BJY1dfsdUt3eTnRQNwMqsBCpauuns7R/1Mzt6+nl8zwV2LE9n2YL4UY8TQoiJkiA/QafrutAabl2bzdP3b0ID9+086Jenr27rISfZ6BuzPNMI2qcuMpr/1XvldPY6eeCaoqBeuxBi/pEgP0He/PqKzHgWp1t44JoiKlq6qWo1Hra63Jra9p6BkXymsdjHaA9fO3v7eWxPGdcuSx+2MIgQQkyVBPkJKq7pJCkmnAWehTouzU0G4EC50fmxsauXfpf2BXmrJZI0S+SoI/mdMooXQgSRBPkJOlnbwcqsBN+yekVpccRHhXGwwgjy3vLJ7KSBNr8rsxI4OajCprffxdGqdn6zt4JH91zg2mVprMqWUbwQIvCkrcEEOJxuzjV0ce/WPN82k0mxITeZA+VtwEDPGu9IHozUzltnG7nrkb1Ut3dT296Ly23k8DPio/jKh5ZM410IIeYTCfITcK6hi36X9uXZvTbkJvHGmUZa7Q7fSD4rcSDIX700jReP1dLndLE2J4lb1sSwIjOelVkJZCVGy2LbQoigmTdBvs/p4nBFO+sXJRERNrkslTevviLTv8zRm5c/VNFGdVs3aRajRt5r7cIk3v7HqyZ55UIIMXkhH+RP13XyzIEq/ni0hvbufh7csZi/HfSQU2vNw8+fIC4yjLs357LwIkvmnaztIDbCPGyB7FVZCUSYTRwsb6W6rccvVSOEEDMppIP8hWY7N/7nu4SbTVy3PJ2qth6e2FfBF68s8C3M8da5Jp4+YMw2fey9C1y7LJ3//eFlIy6CXVzbyfLMeExDFuOICjdzSXYCB8pbabY5WJOTGPybE0KIcQh6dY1S6nql1FmlVIlS6mvBPt9g5c12tIbf3LuR//7EOh64ppCGzj5eLa73HfPI22VkxEfx7lev4ktXFvJeSTPff+XssM9yuTWnajtZkTlyFcyG3GRO1HT41cgLIcRMC2qQV0qZgZ8ANwDLgbuUUsuDec7B6jp6AXyzT7cvTmNhcgw73y8H4Hh1Ox+UtXDv1lxykmP4yoeWcNWSNI7XtA/7rAvNdnr6XcPy8V4b85Lod2mcbu1XPimEEDMp2CP5jUCJ1rpMa+0AngZuDvI5feo7ejApY0ISgNmkuHvzIg6Ut1Fc28HP3ynDEhnGXRsX+t6zIiueqtYeOrr9e814Z7qONit1/cJk3/cykhdCzBbBDvJZQNWg19WebT5KqfuVUgeVUgebmpoCevL6zl6slki/hbHvWJ9DVLiJ7758hpdP1PGJTQuxRIX79nvLI4e2By6u7SQizERhWtyI50qICWeJZ9k9718OQggx02Z8xqvW+hGt9Qat9Qar1RrQz67r6CUjwX9UnRATzq1rs3j3fDNmk+LeLXl++70j9ZNDgvyxqnaWZlj8fmEMtSE3CZPC1/JACCFmWrCDfA2QM+h1tmfbtKjv6GVB/PCAe8/luQDcvCaL9CH7k2MjyEqM5mTNQK+ZboeTw5VtXJaXzMX87dVFPPLpDX418kIIMZOCXUJ5AChSSuVhBPc7gU8E+Zw+9R29bClMHbZ9aUY8O+/dyOpR+sWsyIz36xq570Ir/S7NtqKL/6WRkRBFhozihRCzSFBH8lprJ/A3wKvAaeBZrXVxoM9zvLqdh35/3O9hqa3PSVefc9Sgu32xlcSYkZfaW5mVQFmznS7PQh97zjcTEWZi4xgjeSGEmG2CnpPXWv9Fa71Ya12gtf52MM7RYnPwzMEqSpoGVmiq95RPTiY/vjLLKJM8XWd83p7zzVyamyRpGCHEnDPjD14DId9qzE4tbRxYb9Ub5DNGyMmPZfBCHw2dvZxt6GJrYWAfCgshxHQIibYG2UkxRJhNlDbZfNvqOoxukJPJkafFR5FmieRkbQcJ0UZ55bai4bl9IYSY7UIiyJtNirzUWL8g7x3JD62eGa+VWQmcrOlAa0iJjWC5LLAthJiDQiLIAxSkxXKmblBOvrOX5NiISefRV3oW+mi2OdhSmDqsKZkQQswFIZGTByiwxlHR2o3D6QaMkfxk8vFeK7IScGtotTvYNkIZphBCzAUhFeRdbk1lq/Hwta6jd0ozT1cN6lGzVfLxQog5KmSCvLfCpsRTYVPf2Uv6FIL8goQokmMjKLDGkpkoDceEEHNTyOTk861G47DSJhu9/S5a7Y4RWxqMl1KKB69bTNIoE6aEEGIuCJkgHxcZRkZ8FKVNNho7+4DJlU8O9snLFgXi0oQQYsaETLoGjAqb0ia7r0Z+QYKkWYQQ81toBXlrHGVNNt+KUBkJkTN8RUIIMbNCLsh39To54ekgObSXvBBCzDchFeS9FTbvlTRjiQwjLjJkHjkIIcSkhFSQL/BU2Jyp75K+7kIIQYgF+Yz4KGIijDYGEuSFECLEgrzJpHwpm6m0NBBCiFARUkEeBlI2spi2EEKEYJDPTzWCvFTWCCFECAb5gjRPukZq5IUQIvSC/JVL0vjc1jw25afM9KUIIcSMC7lC8rjIMP73R5bP9GUIIcSsEHIjeSGEEAMkyAshRAiTIC+EECFMgrwQQoSwKQV5pdQdSqlipZRbKbVhyL6HlVIlSqmzSqkPTe0yhRBCTMZUq2tOArcBPx+8USm1HLgTWAFkAq8rpRZrrV1TPJ8QQogJmNJIXmt9Wmt9doRdNwNPa637tNYXgBJg41TOJYQQYuKClZPPAqoGva72bBtGKXW/UuqgUupgU1NTkC5HCCHmpzHTNUqp14GMEXZ9XWv9wlQvQGv9CPCI51xNSqmKSX5UKtA81euZY+Se5we55/lhKve8aLQdYwZ5rfW1kzhhDZAz6HW2Z9tY57JO4lwAKKUOaq03jH1k6JB7nh/knueHYN1zsNI1LwJ3KqUilVJ5QBGwP0jnEkIIMYqpllDeqpSqBjYDf1ZKvQqgtS4GngVOAa8AX5LKGiGEmH5TKqHUWv8B+MMo+74NfHsqnz9Bj0zjuWYLuef5Qe55fgjKPSutdTA+VwghxCwgbQ2EECKESZAXQogQFhJBXil1vadHTolS6mszfT3BoJR6XCnVqJQ6OWhbslJql1LqvOdr0kxeY6AppXKUUm8qpU55eiQ94NkesvetlIpSSu1XSh3z3PM3PdvzlFL7PD/jzyilImb6WgNJKWVWSh1RSv3J8zrU77dcKXVCKXVUKXXQsy0oP9dzPsgrpczAT4AbgOXAXZ7eOaHmV8D1Q7Z9DdittS4CdntehxIn8KDWejmwCfiS579tKN93H3C11no1sAa4Xim1Cfge8COtdSHQBtw3g9cYDA8Apwe9DvX7BbhKa71mUG18UH6u53yQx+iJU6K1LtNaO4CnMXrnhBSt9TtA65DNNwM7Pd/vBG6Z1osKMq11ndb6sOf7LowgkEUI37c22Dwvwz3/NHA18HvP9pC6Z6VUNvBh4FHPa0UI3+9FBOXnOhSC/Lj75ISgdK11nef7eiB9Ji8mmJRSucBaYB8hft+e1MVRoBHYBZQC7Vprp+eQUPsZ/zHwVcDteZ1CaN8vGL+4X1NKHVJK3e/ZFpSf65BbyHu+0lprpVRI1sMqpeKA54C/11p3GgM9Qyjet2fi4BqlVCLGPJSlM3xJQaOU+gjQqLU+pJS6cqavZxpt1VrXKKXSgF1KqTODdwby5zoURvKT6pMTIhqUUgsAPF8bZ/h6Ak4pFY4R4J/UWj/v2Rzy9w2gtW4H3sSYUZ6olPIOykLpZ3wLcJNSqhwj1Xo18B+E7v0CoLWu8XxtxPhFvpEg/VyHQpA/ABR5nsZHYCxW8uIMX9N0eRG4x/P9PcCUu4LOJp7c7GPAaa31DwftCtn7VkpZPSN4lFLRwA6MZxFvArd7DguZe9ZaP6y1ztZa52L8v/uG1vqThOj9AiilYpVSFu/3wHUYCzAF5ec6JGa8KqVuxMjrmYHHPS0VQopS6ingSox2pA3A/wX+iNEjaCFQAXxMaz304eycpZTaCrwLnGAgX/tPGHn5kLxvpdQlGA/dzBiDsGe11t9SSuVjjHSTgSPAp7TWfTN3pYHnSdd8RWv9kVC+X8+9edvBhAG/1Vp/WymVQhB+rkMiyAshhBhZKKRrhBBCjEKCvBBChDAJ8kIIEcIkyAshRAiTIC+EECFMgrwQQoQwCfJCCBHC/j/p3iPkbBu3NgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(xdata, ydata)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we will write the Stan code, keeping in mind that it must be able to compute the pointwise log likelihood on excluded data, that is, data which is not used to fit the model. Thus, the backbone of the code must look like:\n", + "\n", + "```\n", + "data {\n", + " data_for_fitting\n", + " excluded_data\n", + " ...\n", + "}\n", + "model {\n", + " // fit against data_for_fitting\n", + " ...\n", + "}\n", + "generated quantities {\n", + " ....\n", + " log_lik for data_for_fitting\n", + " log_lik_excluded for excluded_data\n", + "}\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "refit_lr_code = \"\"\"\n", + "data {\n", + " // Define data for fitting\n", + " int N;\n", + " vector[N] x;\n", + " vector[N] y;\n", + " // Define excluded data. It will not be used when fitting.\n", + " int N_ex;\n", + " vector[N_ex] x_ex;\n", + " vector[N_ex] y_ex;\n", + "}\n", + "\n", + "parameters {\n", + " real b0;\n", + " real b1;\n", + " real sigma_e;\n", + "}\n", + "\n", + "model {\n", + " b0 ~ normal(0, 10);\n", + " b1 ~ normal(0, 10);\n", + " sigma_e ~ normal(0, 10);\n", + " for (i in 1:N) {\n", + " y[i] ~ normal(b0 + b1 * x[i], sigma_e); // use only data for fitting\n", + " }\n", + " \n", + "}\n", + "\n", + "generated quantities {\n", + " vector[N] log_lik;\n", + " vector[N_ex] log_lik_ex;\n", + " vector[N] y_hat;\n", + " \n", + " for (i in 1:N) {\n", + " // calculate log likelihood and posterior predictive, there are \n", + " // no restrictions on adding more generated quantities\n", + " log_lik[i] = normal_lpdf(y[i] | b0 + b1 * x[i], sigma_e);\n", + " y_hat[i] = normal_rng(b0 + b1 * x[i], sigma_e);\n", + " }\n", + " for (j in 1:N_ex) {\n", + " // calculate the log likelihood of the exluded data given data_for_fitting\n", + " log_lik_ex[j] = normal_lpdf(y_ex[j] | b0 + b1 * x_ex[j], sigma_e);\n", + " }\n", + "}\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_4275bea8cf61cb4b45f01fa01c73d194 NOW.\n" + ] + } + ], + "source": [ + "sm = pystan.StanModel(model_code=refit_lr_code)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "data_dict = {\n", + " \"N\": len(ydata),\n", + " \"y\": ydata,\n", + " \"x\": xdata,\n", + " # No excluded data in initial fit\n", + " \"N_ex\": 0,\n", + " \"x_ex\": [],\n", + " \"y_ex\": [],\n", + "}\n", + "sample_kwargs = {\"iter\": 1000, \"chains\": 4}\n", + "fit = sm.sampling(data=data_dict, **sample_kwargs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have defined a dictionary `sample_kwargs` that will be passed to the `SamplingWrapper` in order to make sure that all\n", + "refits use the same sampler parameters. We follow the same pattern with `az.from_pystan`." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "dims = {\"y\": [\"time\"], \"x\": [\"time\"], \"log_likelihood\": [\"time\"], \"y_hat\": [\"time\"]}\n", + "idata_kwargs = {\n", + " \"posterior_predictive\": [\"y_hat\"],\n", + " \"observed_data\": \"y\",\n", + " \"constant_data\": \"x\",\n", + " \"log_likelihood\": [\"log_lik\", \"log_lik_ex\"],\n", + " \"dims\": dims,\n", + "}\n", + "idata = az.from_pystan(posterior=fit, **idata_kwargs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will create a subclass of {class}`~arviz.PyStan2SamplingWrapper`. Therefore, instead of having to implement all functions required by {func}`~arviz.reloo` we only have to implement `sel_observations`. As explained in its docs, it takes one argument which are the indices of the data to be excluded and returns `modified_observed_data` which is passed as `data` to `sampling` function of PyStan model and `excluded_observed_data` which is used to retrieve the log likelihood of the excluded data (as passing the excluded data would make no sense)." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "class LinearRegressionWrapper(az.PyStan2SamplingWrapper):\n", + " def sel_observations(self, idx):\n", + " xdata = self.idata_orig.constant_data.x.values\n", + " ydata = self.idata_orig.observed_data.y.values\n", + " mask = np.full_like(xdata, True, dtype=bool)\n", + " mask[idx] = False\n", + " N_obs = len(mask)\n", + " N_ex = np.sum(~mask)\n", + " observations = {\n", + " \"N\": N_obs - N_ex,\n", + " \"x\": xdata[mask],\n", + " \"y\": ydata[mask],\n", + " \"N_ex\": N_ex,\n", + " \"x_ex\": xdata[~mask],\n", + " \"y_ex\": ydata[~mask],\n", + " }\n", + " return observations, \"log_lik_ex\"" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Computed from 2000 by 100 log-likelihood matrix\n", + "\n", + " Estimate SE\n", + "elpd_loo -250.66 7.17\n", + "p_loo 2.85 -\n", + "------\n", + "\n", + "Pareto k diagnostic values:\n", + " Count Pct.\n", + "(-Inf, 0.5] (good) 100 100.0%\n", + " (0.5, 0.7] (ok) 0 0.0%\n", + " (0.7, 1] (bad) 0 0.0%\n", + " (1, Inf) (very bad) 0 0.0%\n", + "\n", + "\n", + "The scale is now log by default. Use 'scale' argument or 'stats.ic_scale' rcParam if\n", + "you rely on a specific value.\n", + "A higher log-score (or a lower deviance) indicates a model with better predictive\n", + "accuracy." + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "loo_orig = az.loo(idata, pointwise=True)\n", + "loo_orig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this case, the Leave-One-Out Cross Validation (LOO-CV) approximation using Pareto Smoothed Importance Sampling (PSIS) works for all observations, so we will use modify `loo_orig` in order to make {func}`~arviz.reloo` believe that PSIS failed for some observations. This will also serve as a validation of our wrapper, as the PSIS LOO-CV already returned the correct value." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "loo_orig.pareto_k[[13, 42, 56, 73]] = np.array([0.8, 1.2, 2.6, 0.9])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We initialize our sampling wrapper" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "pystan_wrapper = LinearRegressionWrapper(\n", + " sm, idata_orig=idata, sample_kwargs=sample_kwargs, idata_kwargs=idata_kwargs\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And eventually, we can use this wrapper to call `az.reloo`, and compare the results with the PSIS LOO-CV results." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/percy/anaconda3/envs/arviz/lib/python3.6/site-packages/arviz/stats/stats_refitting.py:98: UserWarning: reloo is an experimental and untested feature\n", + " warnings.warn(\"reloo is an experimental and untested feature\", UserWarning)\n", + "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 13\n", + "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 13\n", + "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 42\n", + "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 42\n", + "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 56\n", + "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 56\n", + "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 73\n", + "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 73\n" + ] + } + ], + "source": [ + "loo_relooed = az.reloo(pystan_wrapper, loo_orig=loo_orig)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Computed from 2000 by 100 log-likelihood matrix\n", + "\n", + " Estimate SE\n", + "elpd_loo -250.67 7.17\n", + "p_loo 2.86 -\n", + "------\n", + "\n", + "Pareto k diagnostic values:\n", + " Count Pct.\n", + "(-Inf, 0.5] (good) 100 100.0%\n", + " (0.5, 0.7] (ok) 0 0.0%\n", + " (0.7, 1] (bad) 0 0.0%\n", + " (1, Inf) (very bad) 0 0.0%\n", + "\n", + "\n", + "The scale is now log by default. Use 'scale' argument or 'stats.ic_scale' rcParam if\n", + "you rely on a specific value.\n", + "A higher log-score (or a lower deviance) indicates a model with better predictive\n", + "accuracy." + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "loo_relooed" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Computed from 2000 by 100 log-likelihood matrix\n", + "\n", + " Estimate SE\n", + "elpd_loo -250.66 7.17\n", + "p_loo 2.85 -\n", + "------\n", + "\n", + "Pareto k diagnostic values:\n", + " Count Pct.\n", + "(-Inf, 0.5] (good) 96 96.0%\n", + " (0.5, 0.7] (ok) 0 0.0%\n", + " (0.7, 1] (bad) 2 2.0%\n", + " (1, Inf) (very bad) 2 2.0%\n", + "\n", + "\n", + "The scale is now log by default. Use 'scale' argument or 'stats.ic_scale' rcParam if\n", + "you rely on a specific value.\n", + "A higher log-score (or a lower deviance) indicates a model with better predictive\n", + "accuracy." + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "loo_orig" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/doc/source/user_guide/pystan_refitting_xr_lik.ipynb b/doc/source/user_guide/pystan2_refitting_xr_lik.ipynb similarity index 99% rename from doc/source/user_guide/pystan_refitting_xr_lik.ipynb rename to doc/source/user_guide/pystan2_refitting_xr_lik.ipynb index c5f2b39678..2cadee432e 100644 --- a/doc/source/user_guide/pystan_refitting_xr_lik.ipynb +++ b/doc/source/user_guide/pystan2_refitting_xr_lik.ipynb @@ -4,8 +4,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "(pystan_refitting_xr)=\n", - "# Refitting PyStan models with ArviZ (and xarray)\n", + "(pystan2_refitting_xr)=\n", + "# Refitting PyStan (2.x) models with ArviZ (and xarray)\n", "\n", "ArviZ is backend agnostic and therefore does not sample directly. In order to take advantage of algorithms that require refitting models several times, ArviZ uses {class}`~arviz.SamplingWrapper`s to convert the API of the sampling backend to a common set of functions. Hence, functions like Leave Future Out Cross Validation can be used in ArviZ independently of the sampling backend used." ] @@ -14,7 +14,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Below there is one example of `SamplingWrapper` usage for PyStan." + "Below there is one example of `SamplingWrapper` usage for PyStan (2.x)." ] }, { @@ -124,7 +124,7 @@ " vector[N] y_hat;\n", " \n", " for (i in 1:N) {\n", - " // pointwise log likelihood will be calculated outside stan, \n", + " // pointwise log likelihood will be calculated outside Stan, \n", " // posterior predictive however will be generated here, there are \n", " // no restrictions on adding more generated quantities\n", " y_hat[i] = normal_rng(b0 + b1 * x[i], sigma_e);\n", @@ -195,7 +195,7 @@ "source": [ "We are now missing the `log_likelihood` group because we have not used the `log_likelihood` argument in `idata_kwargs`. We are doing this to ease the job of the sampling wrapper. Instead of going out of our way to get Stan to calculate the pointwise log likelihood values for each refit and for the excluded observation at every refit, we will compromise and manually write a function to calculate the pointwise log likelihood.\n", "\n", - "Even though it is not ideal to lose part of the straight out of the box capabilities of PyStan-ArviZ integration, this should generally not be a problem. We are basically moving the pointwise log likelihood calculation from the Stan code to the Python code, in both cases we need to manyally write the function to calculate the pointwise log likelihood.\n", + "Even though it is not ideal to lose part of the straight out of the box capabilities of PyStan-ArviZ integration, this should generally not be a problem. We are basically moving the pointwise log likelihood calculation from the Stan code to the Python code, in both cases we need to manually write the function to calculate the pointwise log likelihood.\n", "\n", "Moreover, the Python computation could even be written to be compatible with Dask. Thus it will work even in cases where the large number of observations makes it impossible to store pointwise log likelihood values (with shape `n_samples * n_observations`) in memory." ] @@ -3181,14 +3181,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We will create a subclass of {class}`~arviz.SamplingWrapper`. Therefore, instead of having to implement all functions required by {func}`~arviz.reloo` we only have to implement `sel_observations` (we are cloning `sample` and `get_inference_data` from the `PyStanSamplingWrapper` in order to use `apply_ufunc` instead of assuming the log likelihood is calculated within Stan). \n", + "We will create a subclass of {class}`~arviz.SamplingWrapper`. Therefore, instead of having to implement all functions required by {func}`~arviz.reloo` we only have to implement `sel_observations` (we are cloning `sample` and `get_inference_data` from the `PyStan2SamplingWrapper` in order to use `apply_ufunc` instead of assuming the log likelihood is calculated within Stan). \n", "\n", "Note that of the 2 outputs of `sel_observations`, `data__i` is a dictionary because it is an argument of `sample` which will pass it as is to `model.sampling`, whereas `data_ex` is a list because it is an argument to `log_likelihood__i` which will pass it as `*data_ex` to `apply_ufunc`. More on `data_ex` and `apply_ufunc` integration below." ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -3202,7 +3202,7 @@ " return data__i, data_ex\n", " \n", " def sample(self, modified_observed_data):\n", - " #Cloned from PyStanSamplingWrapper.\n", + " #Cloned from PyStan2SamplingWrapper.\n", " fit = self.model.sampling(data=modified_observed_data, **self.sample_kwargs)\n", " return fit\n", "\n", @@ -3265,7 +3265,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We initialize our sampling wrapper. Let's stop and analize each of the arguments. \n", + "We initialize our sampling wrapper. Let's stop and analyze each of the arguments. \n", "\n", "We then use the `log_lik_fun` and `posterior_vars` argument to tell the wrapper how to call {func}`~xarray:xarray.apply_ufunc`. `log_lik_fun` is the function to be called, which is then called with the following positional arguments:\n", "\n", @@ -3419,5 +3419,5 @@ } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/doc/source/user_guide/pystan_refitting.ipynb b/doc/source/user_guide/pystan_refitting.ipynb index b6c70f1f03..e3f8b2732b 100644 --- a/doc/source/user_guide/pystan_refitting.ipynb +++ b/doc/source/user_guide/pystan_refitting.ipynb @@ -5,7 +5,7 @@ "metadata": {}, "source": [ "(pystan_refitting)=\n", - "# Refitting PyStan models with ArviZ\n", + "# Refitting PyStan (3.0+) models with ArviZ\n", "\n", "ArviZ is backend agnostic and therefore does not sample directly. In order to take advantage of algorithms that require refitting models several times, ArviZ uses {class}`~arviz.SamplingWrapper` to convert the API of the sampling backend to a common set of functions. Hence, functions like Leave Future Out Cross Validation can be used in ArviZ independently of the sampling backend used." ] @@ -26,11 +26,22 @@ "outputs": [], "source": [ "import arviz as az\n", - "import pystan\n", + "import stan\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# enable PyStan on Jupyter IDE\n", + "import nest_asyncio\n", + "nest_asyncio.apply()" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -40,7 +51,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -53,22 +64,22 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "[]" + "[]" ] }, - "execution_count": 3, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xc5Zno8d87oy6N+kiyiq3qblwxNrYxzQRILi2QQAokkJDsZnfZLNkQNjf3brKbTbmbsiWbhACJEwglgQRIQjGmGnDvcpVk9d41ozKamff+cWZGGhWrzaiMnu/n448055yZc84H8ejVc573eZXWGiGEEKHJNNMXIIQQIngkyAshRAiTIC+EECFMgrwQQoQwCfJCCBHCwmb6AgZLTU3Vubm5M30ZQggxpxw6dKhZa20dad+sCvK5ubkcPHhwpi9DCCHmFKVUxWj7JF0jhBAhTIK8EEKEMAnyQggRwiTICyFECJMgL4QQIUyCvBBChDAJ8kIIEcIkyAshxBTZ+5w8c6ASt3v2tW6XIC+EEFP00rFaHnruBHtKmmf6UoaRIC+EEFN0rsEGwOunG/y2t9j6uPvx/ZQ322fisgAJ8kIIMWXnG7sAeP1UA4NX2/v9oWreOdfE0weqZurSAhPklVLlSqkTSqmjSqmDnm3JSqldSqnznq9JgTiXEELMNucbbMRFhlHb0UtxbScAWmueO1wNwMsn65ippVYDOZK/Smu9Rmu9wfP6a8BurXURsNvzWgghQkpHTz/1nb3ctTEHpWDXKSNlU1zbybkGG2tyEqlo6eZUXeeMXF8w0zU3Azs93+8EbgniuYQQYkaUNBr5+E35KaxfmOQL8s8dribCbOJHH1+D2aR4+UT9jFxfoIK8Bl5TSh1SSt3v2Zauta7zfF8PpI/0RqXU/Uqpg0qpg01NTQG6HCGEmB7nG4x8fFGahR3L0zlV10lFi50Xj9Zy7fI08lJj2Zyfwl9OzEzKJlBBfqvWeh1wA/AlpdQVg3dq485GvDut9SNa6w1a6w1W64g974UQYtY632gjKtxEdlI01y43xrLfeKGYFruD29ZmA3DDqgzKmu2+KpzpFJAgr7Wu8XxtBP4AbAQalFILADxfGwNxLiGEmE3ONXRRmBaHyaQosMaRb43lnXNNpMRGsH2JMXC9bnkGJgV/OVE3xqcF3pSDvFIqVill8X4PXAecBF4E7vEcdg/wwlTPJYQQs01Jo43FaRbf6x3LjNH8TWsyCTcbIdZqiWRjXjIvn5yDQR4j175HKXUM2A/8WWv9CvBdYIdS6jxwree1EELMWXvLWvj2n0/52hd09vZT19FLYXqc75hb1mZhtURy18aFfu+9cdUCzjXY+O7LZ/jeK2f41z+doqIl+JOkprzGq9a6DFg9wvYW4Jqpfr4QQswG9R29/NUTh2jr7ufKJWlsKUz1VdYMHskvWxDPga9fO+z916/M4PuvnOVnb5cSYTbR73bT1t3PDz42LHwG1KxayFsIIWYjl1vz5WeO0ud0Y4kK46n9lWwpTPVV1ixOt4zxCZBmieLI/9mBAsLMJh76/XFeOl7Lv9yygpiI4IViaWsghBBj+NnbpXxQ1sI3b1rB7euzebW4nhZbH+cbBiprxiPcbCLMk6f/6Ppsuh0uXjkZ3Pp5CfJCCHERhyvb+OGuc9y0OpPb12dz18aF9LuMlgXnGm2+ypqJ2rAoiZzkaJ4/XBOEqx4gQV4IIS7iNx9UkBAdzrdvXYlSisXpFtYvSuLp/VWcb+iiKG3sVM1ITCbFbWuzea+0mdr2ngBf9aDzBO2ThRAiBByrbmf9oiQsUeG+bXdtXEhZs526jl6KBlXWTNRt67LQGv54NHijeQnyQggxis7efsqa7KzOTvDb/uFVC7BEGQ9LJzuSB1iUEsuluUk8d6g6aC0PJMgLIcQoTlZ3AHBJdqLf9ugIM7euzQJg8RRG8gC3rcumtMnOcc+5Ak2CvBBCjOKYL8gnDNv399cu5v/dfgmLUmKndI4bVy0gIszE857e84EmdfJCCDGK49XtLEqJITEmYti+5NgI7tiQM+VzJESH84u7NwxLCQWKBHkhhBjF8eoO1i0K/qJ22xcHrwOvpGuEEGIETV191LT3BG2EPV0kyAsh5q0z9Z0cq2ofcd/xamP70Ieuc42ka4QQ89aDzx6juLaTj1yygIdvXEZW4kB7gmPVHZgUrMyKn8ErnDoZyQsh5iWtNRUt3eRbY9l1qoFrfvAWv3rvgm//8ep2itIsQW0eNh0kyAsh5qWOnn5sfU4+sXEhux/czub8FP75pVO8fqoBrTXHqztYNcfz8SBBXggxT1W3Gf1ispNiyE6K4aefWs/KrHi+/OxRPihtodXumPMPXUGCvBBinqpu6wbwtQmOCjfz00+uRwGf//VBYO4/dAUJ8kKIEPLC0Ro6uvvHdax3JJ+TFOPblpMcw4/vXIPd4SLcrFi6YPJ9aWYLCfJCiJBQ1drNA08f5akDleM+3hIZRny0/4PVq5em842PLOfuzblEhpmDcanTam4/NhZChCxbn5Pdpxu4aXUmSo29KEdJk7Heqnfd1bFUt/WQlRQ94mfftzVvYhc7i8lIXggxKz21r5IHnj7K+XEG7dLGiQf5nOSYsQ+c4yTICyFmpSNVbQCcrusc1/FlzXYASptsw3qzP7bnAocr23yvtdZUt3WPe23WuUyCvBBiVjpaabQVOFvfNa7jyzzpmq5eJ01dfb7t9j4n//rnUzz6bplvW3t3P3aHi+wkGckLIcS0a+zspbajFxh/kC9tspOZEAUM5OcBTtV1ojUcqxpYlKNqSPlkKJMgL4SYdY54moblJEdzZhxBvqu3n6auPq5bkQEM5OcBTtYYwb2mvcc3wh+YCCVBftyUUmal1BGl1J88r/OUUvuUUiVKqWeUUsO77gshxAiOVrUTZlLcujabmvYeunovXvte1mTk4zcXpBAbYabU8xrgRM3ACN7bWXJgIpSkaybiAeD0oNffA36ktS4E2oD7AnguIUQIO1rZzrIF8b62AucaLj6aL/WkZwqscRSkxflV2BTXdHJZXjImNbCcX3VbD/FRYSREhwfpDmaPgAR5pVQ28GHgUc9rBVwN/N5zyE7glkCcSwgxu93/64N875Uzk36/y605Xt3OmpxElmQYM07HStmUNdkxmxQLk2MotMb5gn6Pw8X5xi4uy0tmcbrF1zu+uq1nXoziIXAj+R8DXwXcntcpQLvW2ul5XQ1kjfRGpdT9SqmDSqmDTU1NAbocIcRMcLrcvHW2iQ9KWyb9GSWNNuwOF2tyEslKjMYSGTbmw9eyZhsLk2OICDNRkBZHXUcvtj4np+s7cWtYkZXA6uxEjle3o7WmqnV+lE9CAIK8UuojQKPW+tBk3q+1fkRrvUFrvcFqDd46h0KI4LvQbMfhcvty3pNx1FMfv2ZhIkopFmdYxhzJlzbaKbDGAkbKBoySymJPPn5lVgKX5CTQ1t1PVWuPjOQnaAtwk1KqHHgaI03zH0CiUsrbNiEbqAnAuYQQs9gpz8SlZpuDbodzjKNHdrSqnYTocPJSjKC9JMPC2fquYROcvFxuzYUWO/me4F6YZryvpNHGiZoOkmMjyEyIYrWno+SbZxvp6XeRkywj+XHRWj+stc7WWucCdwJvaK0/CbwJ3O457B7ghameSwgxuw0ecXvLFCfqSGU7q3MSMZmMnjJLMyx09PTT0Nk34vG17T04nG7yU43gviglljCTorTJxsmaTlZkxqOUYkmGhcgwE38+XgfMj8oaCG6d/EPAPyilSjBy9I8F8VxCiABxutz87O1S7H0TH4mfruskzBOcq1onnrKx9zk519DFmpyBPu5L0r0PX0dub+CrrEkzRvLhZhMLU2I4VdvJuYYuVmUl+LavyIznQEUrMD9q5CHAQV5r/ZbW+iOe78u01hu11oVa6zu01iP/GhZCzCoHytv47stneLW4fsLvPVPXxeaCFODiQd7hdNM5Qu37iZoO3BrWDgrySzOMhbRHe/jqrYn3juQBCq1x7ClpxunWrMwaWN1pdU4i3qxPlgR5IcR85J3yf65hfN0cvdrsDuo7e9lamEpUuImqi6Rr/ukPJ9jynTc4MqhpmMPp5tF3yzApIxh7JcSEkxEfNWqQL2uykRAdTnLswHzLwrQ4+l1GNF81OMh78vIJ0eHER4V+jTxIkBdCDOHNpZc0jq9njNdpTzpl2YJ4spNiRh3JN3b28scjNdgcTu5+bD9Hq9pxON38zW8P8/rpRv75phV+ARuMh6+nRw3ydvKtsX594b0VNgnR4X5pGe8vj/ny0BUkyAshhqie5Ej+TJ0RhJctiCcnKXrUkfxT+6twujVP3ncZSbERfPrRfXz2V/t57VQD37p5BXdvzh32nqUZFkobbfS73MP2lTbZfEHdq9CTn1+ZFe8X/HNTYozAnzg/HrqCBHkhxBDVrUZwrmrrpsfhGvf7Ttd1khoXgdUSSU5yDNWt3cPKHvtdbp7cV8H2xVYuL0zl6fs3kRQbwXslLaMGeDB+cThcbv7p+RNUtBg5eK01J6o7aOzqI98a63d8vjUWk4JVWf4LcSul+O9PrOXLOxaP+77mOln+T4gQca6hi+auPi4vTJ3S51S3GWufdvU5KW2y+T24vJgz9V2+h6Q5STF09Tnp7HGSEDOQ+361uJ7Grj6++9FFAGQmRvP8X1/OhWY7l+Ymj/rZN6zK4FDFIp45WMVzh6vZVmSlpNFGTXsPJgXrFib5HW+JCueJ+y5j2YL4YZ+1rWh+TbqUkbwQIeIHr53lq88dn9JnOJxu6jt72bbY+EVxfpx5eafLzbmGLpYtMModvTnvqiEzX3/9fgULk2PYvjjNty01LvKiAR4gMszMv9yykj1fvYrPb8unpNHGisx4vvfRVex9+Bo25acMe8/lhakkxUrzWxnJCxEiatp7aLb1obUe18LXI6nr6MGtYWuhlV2nGsadly9vsdPndPtG8t6JRlWt3b6/BE7VdrK/vJWv37gMs2ly15cWH8XDNy7j4RuXTer985GM5IUIEbXtvfT2u+meQB59KG9lTW5qDHmpsZwfZ5A/PeihK+BbIHvwSP6JfRVEhZu4Y0P2pK9PTJwEeSFCQI/DRavdAUCLzTHpz/FW1uQkxVCUZhl3uuZMvTHTtcDTN8aoQw+jyvMQ1+XWvHqynh3LM0iMkRTKdJIgL0QIqO0YKFdssY9vcvmxqnZ2/PBt3y8HMEbyZpNiQUIURelxVLZ209s/9l8Gp+u6KEyLIzLM7NuWkxzjG8kfLG+lxe7gQyvSx3tLIkAkyAsRAmrbBwX5cY7kXy2u53yjjb1lA73fq1q7yYiPIsxsoijNgtb4rbI0kiOVbRy40DqskiVn0ISoV4sbiAgzceWStJE+QgSRBHkhQsDgID94ZH4xRyqNVZIOVwy0FjD6rBuVMYvTjQlFFwvyLxyt4eOP7CUpNoK/u6bIb19OcjTVbT243ZpXi+vZWphKXKTUekw3CfJChICa9l7f983jSNe43JpjnkWtD1f6B3nvQ9PcVKNl72jrq/5o1zkeePooa3MS+eOXtpCX6j8hKSc5hj6nm7fPN1HT3iOpmhkiv1aFCAG17T1kxEfR2ds/rnTNuYYuuh0uMhOiOFnTSZ/TyLs3dPX6RvLhZpNRYTPCSP6NMw38x+7z3L4+m3+7dRURYcPHizmeMsrH3r2AScG1yyTIzwQZyQsRAuo6eshMjCIlLmJc6RpvquYzW3JxuNycrOmktr0XrQeCM8DidAvnh4zku3r7+fofTrIk3TJqgIeBCVF7SprZkJtMSlzkZG9PTIEEeSFCQG17L5mJ0STHRtJsGztdc6SyjeTYCG5ZmwUYeXlv+eTgro2FacMrbL7/ylnqO3v57kdHD/DG5wz8svjQiowJ35MIDAnyQsxxWmtq2nvISowmNTZiXOmaI1XtrM1JJM0SRU5yNIcr23w17dnJ/iN5tx5YfWn/hVZ+s7eCz16ex9oh/WKGigo3Y7UYo/frlkuqZqZITl6IOa7F7sDhdJOZGE2r3UFx7cjL5Hl19PRT0mjjljWZgNHca29Zi+9Ba0Z8lO/YJRlGhc0tP3mPNEsUdoeT7KRovvKh8XVxLLTGsSAhyvcwV0w/CfJCzHHe8snMxGjqOnppsV+8f82xKiMf7x2Jr1+UxAtHa9l/oZXMxGi/vjIF1jh+cMdqSppsNHT00mx38DdXFRITMb7Q8aOPr5nKrYkAkCAvxBw3EOSjqGiJoN+l6epzjrq83ZHKdpSCS7KNxmHeNr2HKtq4vMC/m6NSio+un3yvmYyEqLEPEkElOXkh5jhvjXxWYrRv2byL5eWPVLWxOM2CxfNLYGmGhehwox3B4MoaERokyAsxx9W29xATYSYhOtxXptgySoWN1pojle2sXTiwYlKY2cTqHGNUP7iyRoQGCfJCzHG17T1kJkajlCLFO5IfpVb+QrOdjp5+vyAPAymb7Hm0wPV8IUFeiDnOG+QBUuIunq45UN4KMKz88fICYyWoojRLsC5TzJApB3mlVJRSar9S6phSqlgp9U3P9jyl1D6lVIlS6hmllDSRFiIIatp7yUo0HnB6c/KtI/Sv0Vqz8/0K8q2xFFrj/PZtLUpl94Pbx72eq5g7AjGS7wOu1lqvBtYA1yulNgHfA36ktS4E2oD7AnAuIcQgvf0umm19LEgwRvKRYWYskWE0jzCSf+d8M6fqOvniFQWYRlh+r2BI4BehYcpBXhu8HYzCPf80cDXwe8/2ncAtUz2XEMJffYdRWeNN14CRshkpJ//Tt0rIiI/i5rWZ03Z9YuYFJCevlDIrpY4CjcAuoBRo11o7PYdUA1mjvPd+pdRBpdTBpqamQFyOEPPG4Bp5r+TYiGHpmiOVbewta+Vz2/L8Vm8SoS8gQV5r7dJarwGygY3A0gm89xGt9Qat9Qar1RqIyxFi3qjxBPksv5F85LAHrz97u5SE6HDu3LhwWq9PzLyAVtdorduBN4HNQKJSyjujNhuoCeS5hBBG90nwn1maOiRdU9LYxavFDdyzeZGszDQPBaK6xqqUSvR8Hw3sAE5jBPvbPYfdA7ww1XMJIfzVdfRgtUT6pWCMdI0Dt1sD8MyBKiLMJu65PHeGrlLMpECM5BcAbyqljgMHgF1a6z8BDwH/oJQqAVKAxwJwLiHEIGVN9mGzVFNiI3G5NR09/QC8c66ZS/OSZNGOeWrKf7tprY8Da0fYXoaRnxdCBIGtz8nhyjY+f0W+33bfhCi7A4fLzdmGLm5bN+7HZCLESIJOiFlk5/vlhJkVVy1J8yuLHMn7Jc043Zrti/0LFlJiB/rXHPW0Fd5WJEUN85UEeSGmQV1HD9/4YzFf3lHEisyRZ5V2dPfzf18s9r1emmHhmzet4LL8lBGPf/tcE3GRYb6+M16DR/Lvnm8iNS6SpRnSrmC+kt41QkyDX71XzuunG7j3Vweo6+gZ8ZgqzxqrD9+wlIdvWEpHTz/ffOkUWuthx2qteftcE5cXpAxbZ9XbpKzZ1see881sK0odcYarmB8kyAsRZH1OF787VM2anETsfS4++8sDdPX2Dzuuus0I/lsKU/nC9gL+9uoiTtV1cqiibdixZc12qtt62L5keBomyRPk3z3fTIvdwbai1ADfkZhLJMgLEWSvFjfQanfw5R2L+emn1lHSaONLvz1Cv8vtd9zQiU23rM0kPiqMX71fPuwz3z5rzA6/YoRce7jZREJ0OG+dbQRga6EE+flMgrwQQfbbfRXkJEezrTCVbUVWvnXzSt4518Trpxr8jqtu6yY2wkxijLFiU0xEGB/bkMMrJ+tp6Oz1O/btc03kW2NHXSA7Jc5YBnBphoW0eFmCbz6TIC9EEJU22dhb1sqdly705cVvW5eFUnCuweZ3bHVbD9lJMX4LcH968yJcWvPkvkrftt5+F3vLWoZV1QyW6qmwueIix4j5QYK8EEH01L5KwkyKOzYMLIYdFW4mKzGaC83+Qb6mrYesIRObFqXEctWSNH67rxKH00jv7LvQSp/TfdEg7+0rL/l4IUFeiCDp7Xfx+8PVXLcinTSLf8okLzWWsma737bqtu4R11i95/Jcmm19fOfl0zx3qJqn9lUSGWZi0yillWD0sokON3NpbnJgbkbMWVInL8QU9bvc1Hf0DsuPv3Gmkfbufu4aofNjfmoszx+uQWuNUorO3n46e51+3SS9thWmsmxBPL98r9y37dplaUSFj94y+EtXFfLRddkXPUbMDxLkhZii/36jhJ++Xcp7D12N1TLQH+b1Uw0kxoSzeYQRd15qLF19TpptDqyWSGo85ZPZScMfpJpMipf+ZgvtPf1097mw9TlZmDLyA1cvqyXS71rE/CXpGiGmwOXWPHOgCofTzWun6v22v3m2kauWpBFmHv6/Wb5nqb0LnpRNtS/Ij9zKIMxsIjUukoUpMSzPjJeWwWLcJMgLMQXvnm+ivrOXMJPilZMDQf5oVRtt3f1cvTRtxPflpcYCUNZkPHyt8cx2HfrgVYipkiAvxBT87mA1ybERfHZLLu+XttDmWaxj9+lGzCY1agljZmI0EWEmv5F8VLjJ15JAiECRIC/EJLXZHew61cDNazK5eU0WLrdml2eC0+7TjVyam0RCdPiI7zWbFLkpMb4Km+q2HrISo/1q5IUIBAnyQkzSC0drcLjc3LE+hxWZ8eQkR/OXk3VUtXZztqGLa5elX/T9eamxvpF8TXvPiA9dhZgqCfJCTNKzB6tZlZXA8sx4lFLcuHIB75U088JRYznj0fLxXvnWOCpa7Dhd7lFr5IWYKgnyQkzCyZoOTtV1+s1kvX5lBv0uzU/eLCUvNdZXQTOavNRY+l2acw022rr75aGrCAqpwxJiiN5+F8CwiUT9LjdvnW3ixWO1vH6qgahwEzetzvTtX5OTSGZCFLUdvWOO4sGYEAWwp8ToKCnpGhEMEuSFGOKvnjhEn9PNbz+/yW/7/3mhmKf2V5IUE86t67L4xMaFJMYMVMMopbh+5QIef+8C1ywbO8h7yyjfPd8MMOJsVyGmSoK8CGn9Ljdf/f1xPnN5LqtzEsc83u3W7L/Qit3h4mRNByuzjKX6mrr6eO5QNXesz+bfbltF+AgTnAA+ty2P2EgzG8fRMyY5NoKE6HD2XWgFIEfSNSIIJCcvQtqBC6384UgNvx3Uqvdiylvs2B1GumbwYh1P7a/E4XLzxSsLRg3wYNS/P3jdkhFnuQ6llCIvNRaH001EmDGjVYhAkyAvQtprnrr1PSXNw9ZKffj547x4rNZvW3FtJwAbFiXx4rFaWmx9OJxunthbwfbFVgrGeJg6Ud68fFZitKzDKoJCgrwIWVprXj/dQITZRE17j68mHYyeMU/tr+KJvRV+7zlZ20GE2cS3bl6Jw+nm6QNVvHyyjsauPj6zJTfg1+jNy0v5pAiWKQd5pVSOUupNpdQppVSxUuoBz/ZkpdQupdR5z9ekqV+uEON3tqGL6rYePrs1FzBG814vn6wD4GhVu6+aBqC4ppMlGRaWZ8aztTCVJ/ZW8Ph75eSlxrJ9hPVUp8pbZikPXUWwBGIk7wQe1FovBzYBX1JKLQe+BuzWWhcBuz2vhZg23jVU79uSR3ZStK+KBeCVk/VEhJlwON0crWoHjJF/cW0HK7PiAWOxjrqOXo5VtXP35kVBSafISF4E25SDvNa6Tmt92PN9F3AayAJuBnZ6DtsJ3DLVc4nQsP9CK88dqsbpck/5s4bm2QfbdbqR1TmJpMVHsa0olb2lLb7ZpcerO7h3Sx5Kwb4yo7qltqOXtu5+lmcaFTVXL00jJzma2Agzt6/PHvU8U1GYFsf/Wp3JNWO0QBBisgKak1dK5QJrgX1Auta6zrOrHpCfYgHAv796lgd/d4yP/NcePihtmfTnPHOgks3feYPTdZ3D9jV2GiPwHZ569a2FVrr6nByrbve1BL7z0hyWZsSzv9y4hpM1HQCszDRG8maT4scfX8tPPrkOS9TIjcamKiLMxH/dtZZlC+KD8vlCBCzIK6XigOeAv9da+/1fp43h1ohDLqXU/Uqpg0qpg01NTYG6HDGLNdv6WJweR1evk7t+sZd/frF4wp/x2J4LPPTcCeo7e/nP3eeH7d99phGAHcszALi8IAWljIlHr5ysZ2mGhdzUWC7LS+ZQRRsOp5vi2k5MCpZmDATc9YuSuHLJ2BObhJitAhLklVLhGAH+Sa31857NDUqpBZ79C4DGkd6rtX5Ea71Ba73Bag38gy0x+zTb+ticn8LuB7dz+/psdn5QTlVr97jeq7Xmv3af51/+dIobVmbwhe35vHyynnMNXX7HvX6qgZzkaBanGw82k2IjWJWVwIvHajlU2cYNKxcAsCk/md5+Nydq2imu6aAwLY7oCFkXVYSOQFTXKOAx4LTW+oeDdr0I3OP5/h7ghameS8x9Dqebzl4nKXGRRIWb+fKOxSjg2YNV43r/66cb+cGuc9y2Lov/umstX7yigJgIMz95s8R3TLOtjz0lzVy7LN2vP/vWwlTKmuxoDTesMkb4G/OM9Vf3lrVysraDlZ58vBChIhAj+S3Ap4GrlVJHPf9uBL4L7FBKnQeu9bwW81yrZ+WklDij50tWYjTbF1t59mDVuB7EvnGmAUtUGN//6CWEmU0kxUbw6U2LeOlYLRea7TR29XLXI3tRCu5Yn+P33q1FqQDkW2MpSjNG+MmxESxOj+MvJ+po6OxjeabkxkVoCUR1zR6ttdJaX6K1XuP59xetdYvW+hqtdZHW+lqtdWsgLljMbc22PgBSYgem8N+1cSENnX28eXbsZzIflLZwWV6yX9uA+7blEW428Z2/nObOR/ZS097DLz+zcVjAXr8oidS4CD66LttvhH9ZXopvpqu3V40QoUJmvIpp1eIZyafGDXRvvHppGmmWSJ7af/H+MrXtPZS3dLO5INVve5olirs2LuS1Uw00dPSy896NbC5IGfb+yDAz73z1Kv5qe4Hf9svyB5qJyUhehBrpQimmVYt3JD+oGVeY2cTHNuTwP2+VUNveQ+Yosz+95ZaXjxDA//rKAmrbe/jC9gLWLxp9cnVMxPAf+Y15RpBflBJDfJBKJYWYKTKSF9Oqxeafk/f6+KU5aC7+AMNXefAAABe3SURBVPb90haSYsJZkm4Zti8tPopH7t5w0QA/mjRLFCsy48fVHliIuUZG8mJaNdv7iDCbsET6/+jlJMewrcjKbz6o4COXZFKY5t/tUWvN3rIWNhekBKW9wNP3b7poC2Eh5ir5qRbTqsXmICUuwu/Bp9c3PrwMpRQf//kHFNd2+O2rbO2mpr2HzfnDUzWBYIkKH7bcnxChQIK8mFYttr5hqRqvonQLz35hE5FhJu56ZC+HK9t8+7z5+KEPXYUQFydBXkxYb7+L77x8etyzVAdrsTv8yieHyrfG8ewXN5McG8GnHt3HvjIjuL9f2oLVEkmBNXbS1y3EfCRBXkyI1pqv/+EkP3+7jJeO1479hiG86ZqLyU6K4dkvbGZBQhSf+eUBPiht4YOyFk//GVk9SYiJkCAvJuTJfZU8d7gagLIm+xhH+9Na02LvG9dapmnxUTx1/yayk6K5+/F9NHX1BS0fL0QokyAvxu1wZRvffKmYK5dY2ZibTFmTbULv73a46O13kxJ78ZG8V5rFCPT5qUalzeWSjxdiwiTIi3Hp6Onnr584TEZCFD/++BoK0uL81kwdj4Ea+bFH8l6pcZE8+4XN/O6Lm1mYEjOh8wkhJMjPSv0uN3c/vp+D5bOn3c+T+yqM3u13riUxJoICayxt3f20edoUjEez3TvbdXwjea+EmHAulYlKQkyKBPlZqL6jl3fONfH66RFb8E/aa8X1PPz8iYsumTeSPqeLX71XztbCVNYuNGaU5nuqXMqax5+y8Y7kUy9SXSOECCwJ8rNQk6e/y4UJBNCx9Pa7+MYLJ3lqfyV7yyb2F8KLR2tp7Orj81fk+7Z58+SlE3j4OtC3ZmIjeSHE5EmQn4W8I96Rct6//qCcQxVtw7aP5dcflNPQ2Ud0uJlfvndh3O/TWvPouxdYkm7hiqKBB5/ZSdGEm9WEKmy8HSiTx/ngVQgxdRLkZyHviLe8pRu3eyC14nC6+dZLpyYUpAE6e/v5n7dK2b7Yyr1bc9l1umHcE5neOd/M2YYuPn9Fvl+NepjZxMLkmIv+tfFqcT1n6geW+2229REXGSbtA4SYRhLkZyHvwhoOp5vajh7f9rJmG0635nzDxNI4j757gfbufv7xQ0v49KZczEqx8/3ycb33F++UkR4fyU2rM4fty7fGjTqSf6+kmS8+cYjvv3LWt208E6GEEIElQX4WarYNVKwMTtmcrTcWqy5rto1rqTww/ip47N0yPrxqASuzEshIiOKGVQt45mAV9j4nTpebH752lqv//S1fCwGvF4/Vsqekmc9cnkdE2PAflXxrLBUt3bjc/g9yW+0O/uHZo2gNRyrbfA96W+x9466RF0IEhgT5Wcib1gD/IH/GE+T7XZqKcaZbHn/vAj39Lv7husW+bZ/dkktXr5P/eauEu36xl/98o4TWbgefemwfvztYhdaaH79+jr976ggbFiXx6c2LRvzsgtQ4HC431W0D16K15qHnjtNqd/DpTYto6+733YMxkpfKGiGmkwT5WajF5mBJhoWYCLNfOuRcfZdvRH2+oWtcn3W8uoNVWQkUWAf6s69bmMTqnER+8mYpp2o7+dHHV/P2V67isrwU/vH3x/nwf+7hx6+f57Z1WTz5+ct8v3CG8pVRDrrGJ/dVsutUAw9dv5RPbTJ+ORyubAeMv1BSJV0jxLSSID8Ltdj7sMZFkpsSS3mL/0j+iiIrwLjz8uUtdnJTh3du/Nr1S7lueTp/+rtt3Lo2m4SYcH752Uv51KaFnK7v5Gs3LOUHd6wmMmz0h6T5nl8cZZ6Renu3g3/7y2m2FaVy75Y8itLisESGcbiyDbdb02rvu2gHSiFE4MnKULNQs83BpbkRmM2KkzXG4hldvf3UtPfwicsWcqa+k3ONYwd5h9NNTVsPt67NHrZvc0HKsMWuw80m/vWWVTx0/VIs41jrNCkmnITocF8Pmyf3VdLtcPHwDct8qzetWZjI4Yo22nv6cWupkRdiuslIfpZxuty0dRu56/zUWKpau3E43ZzzjNyXpFsoSosbV7qmuq0bt4bcCfZ8GU+AB1BKkW+NpazJjsPpZuf7xqzY5ZnxvmPWLUziXEOX7y8SyckLMb0kyM8yrd0OtAZrXAR5qbG4NVS1dfsqa5ZkWFicbqGs2T5mhU1Fi/FAdFEQG3vlp8ZR1mzjpWPGrNjPbcvz279uURJuDW94WjSkSnWNENNKgvwsM7hTozeXfqHJzrmGLmIjzGQlRlOYFofD6aZyjAob7+h5UUrwVlPKt8bS0NnHT94qoSgtju2LrX771+QkArDrVAMgI3khpltAgrxS6nGlVKNS6uSgbclKqV1KqfOer0mBOFeo8zXx8qRrwCijPFPfyeIMCyaTYnG6BYDzY+TlK1q6iYsMC2ptesGgCpvPbcsbtnJTQnQ4RWlxnPWklyQnL8T0CtRI/lfA9UO2fQ3YrbUuAnZ7XosxNA9q4pUYE0FSTDhlzXbO1nexxBPcC9KMqpax8vLlLXYWpcQEdck8b4VNalwEN6/JGvGYdZ7OlUpBUowEeSGmU0CCvNb6HWBoa8ObgZ2e73cCtwTiXKHOG+S97XjzUmM5UN5KW3c/SzKMIB8XGUZWYvS4RvK5QUzVgJHvT4wJ5/Pb8kftSbN+kRHkk2MiMJtkjVYhplMwc/LpWus6z/f1QPpIByml7ldKHVRKHWxqagri5cwNzTYH4WZFfLRR3ZqbGktJ40BljVdRetxFa+WdLjdVrd1BfegKEBlm5oOvXcP9g9oQD7VukZGXl1SNENNvWh68aqN5yYgrVWitH9Fab9Bab7BarSMdMq+02IwJQ94US/6giUzekTxAUVocpU02XG6N1ppvvXSKn75V6ttf296L062DPpIHiI4wXzQllJ8aR3xUmEyEEmIGBHMyVINSaoHWuk4ptQAI7DJHIarF7iDVMjDizUv15rwj/SpTitIs9DmN0fqfT9Tx+HsXSIoJ5/4r8jGb1KDKmplfF9VkUnx5x2LpIy/EDAjmSP5F4B7P9/cALwTxXCGj2eY/9T/PM5JfkhHnd1xRuvH6Z2+X8u+vnWVhcgxt3f2c8MyQ9TYwG6mlwUz47Ja8UR/MCiGCJ1AllE8BHwBLlFLVSqn7gO8CO5RS54FrPa/FGFpsDlIHjdhzU2NQCpZmxPsdV+ipsHn6QBUrMxN45gubUAreOmv8wVTRbCcq3ESaRVIkQsxnAUnXaK3vGmXXNYH4/PlCa02Trc+vU2NMRBiP3r2BlVkJfsdaosLJSozG4XLzi7s3kJEQxersRN4628TfX7uYck9lTTDLJ4UQs580KJtFbH1OHE73sCqUa5aNWJjETz+1jviocDISogC4comV/9h9nla7g4oWu68VsBBi/pK2BrPI4Nmu43FJdqJfzv3KJWloDW+fa6SiNfg18kKI2U+C/CwyMNt1cnn0S7ISSI6N4On9VTic7qD2rBFCzA0S5IOkze7glZN1Yx84SLNvJD+5UkOTSbGtKJV9F4zJxxNtMSyECD0S5IPkqQOVfPGJw9S094z7Pb6WBlPo1HjlkoEJZYtmSfmkEGLmSJAPkipPnbp3Zafx8ObkpzJp6IoiK0pBhNlERnzUpD9HCBEaJMgHSXWbMYIvru0c93ta7H0kxoQTbp78f5aUuEguyU5kUUqMNAMTQkgJZbB4R/LFExjJG7Ndpz71/99vv4SefteUP0cIMfdJkA8Ct1v7cvETGck3D5ntOllFg7pVCiHmN0nXBEFjVx/9Lk1eaiz1nb2+B6pjabb1BSTICyGElwT5IKhuM1I1H1qRAYx/NN9ic0jPdSFEQEmQnwSjPf7ovA9dr1thtCMYT4WNw+mmo6dfRvJCiICSID8J33n5DNf96G3fiH0o7/blC+JZmBxDce3YQb5p0NquQggRKBLkJ+GNM42ca7Dx8Z/vpbzZPmx/VWsPqXGRRIWbWZEZP650ze8OVgEDi14LIUQgSJCfoG6Hk9ImGzeuyqCn38XHfv4BJY1dfsdUt3eTnRQNwMqsBCpauuns7R/1Mzt6+nl8zwV2LE9n2YL4UY8TQoiJkiA/QafrutAabl2bzdP3b0ID9+086Jenr27rISfZ6BuzPNMI2qcuMpr/1XvldPY6eeCaoqBeuxBi/pEgP0He/PqKzHgWp1t44JoiKlq6qWo1Hra63Jra9p6BkXymsdjHaA9fO3v7eWxPGdcuSx+2MIgQQkyVBPkJKq7pJCkmnAWehTouzU0G4EC50fmxsauXfpf2BXmrJZI0S+SoI/mdMooXQgSRBPkJOlnbwcqsBN+yekVpccRHhXGwwgjy3vLJ7KSBNr8rsxI4OajCprffxdGqdn6zt4JH91zg2mVprMqWUbwQIvCkrcEEOJxuzjV0ce/WPN82k0mxITeZA+VtwEDPGu9IHozUzltnG7nrkb1Ut3dT296Ly23k8DPio/jKh5ZM410IIeYTCfITcK6hi36X9uXZvTbkJvHGmUZa7Q7fSD4rcSDIX700jReP1dLndLE2J4lb1sSwIjOelVkJZCVGy2LbQoigmTdBvs/p4nBFO+sXJRERNrkslTevviLTv8zRm5c/VNFGdVs3aRajRt5r7cIk3v7HqyZ55UIIMXkhH+RP13XyzIEq/ni0hvbufh7csZi/HfSQU2vNw8+fIC4yjLs357LwIkvmnaztIDbCPGyB7FVZCUSYTRwsb6W6rccvVSOEEDMppIP8hWY7N/7nu4SbTVy3PJ2qth6e2FfBF68s8C3M8da5Jp4+YMw2fey9C1y7LJ3//eFlIy6CXVzbyfLMeExDFuOICjdzSXYCB8pbabY5WJOTGPybE0KIcQh6dY1S6nql1FmlVIlS6mvBPt9g5c12tIbf3LuR//7EOh64ppCGzj5eLa73HfPI22VkxEfx7lev4ktXFvJeSTPff+XssM9yuTWnajtZkTlyFcyG3GRO1HT41cgLIcRMC2qQV0qZgZ8ANwDLgbuUUsuDec7B6jp6AXyzT7cvTmNhcgw73y8H4Hh1Ox+UtXDv1lxykmP4yoeWcNWSNI7XtA/7rAvNdnr6XcPy8V4b85Lod2mcbu1XPimEEDMp2CP5jUCJ1rpMa+0AngZuDvI5feo7ejApY0ISgNmkuHvzIg6Ut1Fc28HP3ynDEhnGXRsX+t6zIiueqtYeOrr9e814Z7qONit1/cJk3/cykhdCzBbBDvJZQNWg19WebT5KqfuVUgeVUgebmpoCevL6zl6slki/hbHvWJ9DVLiJ7758hpdP1PGJTQuxRIX79nvLI4e2By6u7SQizERhWtyI50qICWeJZ9k9718OQggx02Z8xqvW+hGt9Qat9Qar1RrQz67r6CUjwX9UnRATzq1rs3j3fDNmk+LeLXl++70j9ZNDgvyxqnaWZlj8fmEMtSE3CZPC1/JACCFmWrCDfA2QM+h1tmfbtKjv6GVB/PCAe8/luQDcvCaL9CH7k2MjyEqM5mTNQK+ZboeTw5VtXJaXzMX87dVFPPLpDX418kIIMZOCXUJ5AChSSuVhBPc7gU8E+Zw+9R29bClMHbZ9aUY8O+/dyOpR+sWsyIz36xq570Ir/S7NtqKL/6WRkRBFhozihRCzSFBH8lprJ/A3wKvAaeBZrXVxoM9zvLqdh35/3O9hqa3PSVefc9Sgu32xlcSYkZfaW5mVQFmznS7PQh97zjcTEWZi4xgjeSGEmG2CnpPXWv9Fa71Ya12gtf52MM7RYnPwzMEqSpoGVmiq95RPTiY/vjLLKJM8XWd83p7zzVyamyRpGCHEnDPjD14DId9qzE4tbRxYb9Ub5DNGyMmPZfBCHw2dvZxt6GJrYWAfCgshxHQIibYG2UkxRJhNlDbZfNvqOoxukJPJkafFR5FmieRkbQcJ0UZ55bai4bl9IYSY7UIiyJtNirzUWL8g7x3JD62eGa+VWQmcrOlAa0iJjWC5LLAthJiDQiLIAxSkxXKmblBOvrOX5NiISefRV3oW+mi2OdhSmDqsKZkQQswFIZGTByiwxlHR2o3D6QaMkfxk8vFeK7IScGtotTvYNkIZphBCzAUhFeRdbk1lq/Hwta6jd0ozT1cN6lGzVfLxQog5KmSCvLfCpsRTYVPf2Uv6FIL8goQokmMjKLDGkpkoDceEEHNTyOTk861G47DSJhu9/S5a7Y4RWxqMl1KKB69bTNIoE6aEEGIuCJkgHxcZRkZ8FKVNNho7+4DJlU8O9snLFgXi0oQQYsaETLoGjAqb0ia7r0Z+QYKkWYQQ81toBXlrHGVNNt+KUBkJkTN8RUIIMbNCLsh39To54ekgObSXvBBCzDchFeS9FTbvlTRjiQwjLjJkHjkIIcSkhFSQL/BU2Jyp75K+7kIIQYgF+Yz4KGIijDYGEuSFECLEgrzJpHwpm6m0NBBCiFARUkEeBlI2spi2EEKEYJDPTzWCvFTWCCFECAb5gjRPukZq5IUQIvSC/JVL0vjc1jw25afM9KUIIcSMC7lC8rjIMP73R5bP9GUIIcSsEHIjeSGEEAMkyAshRAiTIC+EECFMgrwQQoSwKQV5pdQdSqlipZRbKbVhyL6HlVIlSqmzSqkPTe0yhRBCTMZUq2tOArcBPx+8USm1HLgTWAFkAq8rpRZrrV1TPJ8QQogJmNJIXmt9Wmt9doRdNwNPa637tNYXgBJg41TOJYQQYuKClZPPAqoGva72bBtGKXW/UuqgUupgU1NTkC5HCCHmpzHTNUqp14GMEXZ9XWv9wlQvQGv9CPCI51xNSqmKSX5UKtA81euZY+Se5we55/lhKve8aLQdYwZ5rfW1kzhhDZAz6HW2Z9tY57JO4lwAKKUOaq03jH1k6JB7nh/knueHYN1zsNI1LwJ3KqUilVJ5QBGwP0jnEkIIMYqpllDeqpSqBjYDf1ZKvQqgtS4GngVOAa8AX5LKGiGEmH5TKqHUWv8B+MMo+74NfHsqnz9Bj0zjuWYLuef5Qe55fgjKPSutdTA+VwghxCwgbQ2EECKESZAXQogQFhJBXil1vadHTolS6mszfT3BoJR6XCnVqJQ6OWhbslJql1LqvOdr0kxeY6AppXKUUm8qpU55eiQ94NkesvetlIpSSu1XSh3z3PM3PdvzlFL7PD/jzyilImb6WgNJKWVWSh1RSv3J8zrU77dcKXVCKXVUKXXQsy0oP9dzPsgrpczAT4AbgOXAXZ7eOaHmV8D1Q7Z9DdittS4CdntehxIn8KDWejmwCfiS579tKN93H3C11no1sAa4Xim1Cfge8COtdSHQBtw3g9cYDA8Apwe9DvX7BbhKa71mUG18UH6u53yQx+iJU6K1LtNaO4CnMXrnhBSt9TtA65DNNwM7Pd/vBG6Z1osKMq11ndb6sOf7LowgkEUI37c22Dwvwz3/NHA18HvP9pC6Z6VUNvBh4FHPa0UI3+9FBOXnOhSC/Lj75ISgdK11nef7eiB9Ji8mmJRSucBaYB8hft+e1MVRoBHYBZQC7Vprp+eQUPsZ/zHwVcDteZ1CaN8vGL+4X1NKHVJK3e/ZFpSf65BbyHu+0lprpVRI1sMqpeKA54C/11p3GgM9Qyjet2fi4BqlVCLGPJSlM3xJQaOU+gjQqLU+pJS6cqavZxpt1VrXKKXSgF1KqTODdwby5zoURvKT6pMTIhqUUgsAPF8bZ/h6Ak4pFY4R4J/UWj/v2Rzy9w2gtW4H3sSYUZ6olPIOykLpZ3wLcJNSqhwj1Xo18B+E7v0CoLWu8XxtxPhFvpEg/VyHQpA/ABR5nsZHYCxW8uIMX9N0eRG4x/P9PcCUu4LOJp7c7GPAaa31DwftCtn7VkpZPSN4lFLRwA6MZxFvArd7DguZe9ZaP6y1ztZa52L8v/uG1vqThOj9AiilYpVSFu/3wHUYCzAF5ec6JGa8KqVuxMjrmYHHPS0VQopS6ingSox2pA3A/wX+iNEjaCFQAXxMaz304eycpZTaCrwLnGAgX/tPGHn5kLxvpdQlGA/dzBiDsGe11t9SSuVjjHSTgSPAp7TWfTN3pYHnSdd8RWv9kVC+X8+9edvBhAG/1Vp/WymVQhB+rkMiyAshhBhZKKRrhBBCjEKCvBBChDAJ8kIIEcIkyAshRAiTIC+EECFMgrwQQoQwCfJCCBHC/j/p3iPkbBu3NgAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8+yak3AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3tElEQVR4nO3dd3xc5Zno8d87oy6N+kiyiq3qblwxNrYxzQRILi2QQAokkJDsZnfZLNkQNjf3brKbTbmbsiWbhACJEwglgQRIQjGmGnDvcpVk9d41ozKamff+cWZGGhWrzaiMnu/n448055yZc84H8ejVc573eZXWGiGEEKHJNNMXIIQQIngkyAshRAiTIC+EECFMgrwQQoQwCfJCCBHCwmb6AgZLTU3Vubm5M30ZQggxpxw6dKhZa20dad+sCvK5ubkcPHhwpi9DCCHmFKVUxWj7JF0jhBAhTIK8EEKEMAnyQggRwiTICyFECJMgL4QQIUyCvBBChDAJ8kIIEcIkyAshxBTZ+5w8c6ASt3v2tW6XIC+EEFP00rFaHnruBHtKmmf6UoaRIC+EEFN0rsEGwOunG/y2t9j6uPvx/ZQ322fisgAJ8kIIMWXnG7sAeP1UA4NX2/v9oWreOdfE0weqZurSAhPklVLlSqkTSqmjSqmDnm3JSqldSqnznq9JgTiXEELMNucbbMRFhlHb0UtxbScAWmueO1wNwMsn65ippVYDOZK/Smu9Rmu9wfP6a8BurXURsNvzWgghQkpHTz/1nb3ctTEHpWDXKSNlU1zbybkGG2tyEqlo6eZUXeeMXF8w0zU3Azs93+8EbgniuYQQYkaUNBr5+E35KaxfmOQL8s8dribCbOJHH1+D2aR4+UT9jFxfoIK8Bl5TSh1SSt3v2Zauta7zfF8PpI/0RqXU/Uqpg0qpg01NTQG6HCGEmB7nG4x8fFGahR3L0zlV10lFi50Xj9Zy7fI08lJj2Zyfwl9OzEzKJlBBfqvWeh1wA/AlpdQVg3dq485GvDut9SNa6w1a6w1W64g974UQYtY632gjKtxEdlI01y43xrLfeKGYFruD29ZmA3DDqgzKmu2+KpzpFJAgr7Wu8XxtBP4AbAQalFILADxfGwNxLiGEmE3ONXRRmBaHyaQosMaRb43lnXNNpMRGsH2JMXC9bnkGJgV/OVE3xqcF3pSDvFIqVill8X4PXAecBF4E7vEcdg/wwlTPJYQQs01Jo43FaRbf6x3LjNH8TWsyCTcbIdZqiWRjXjIvn5yDQR4j175HKXUM2A/8WWv9CvBdYIdS6jxwree1EELMWXvLWvj2n0/52hd09vZT19FLYXqc75hb1mZhtURy18aFfu+9cdUCzjXY+O7LZ/jeK2f41z+doqIl+JOkprzGq9a6DFg9wvYW4Jqpfr4QQswG9R29/NUTh2jr7ufKJWlsKUz1VdYMHskvWxDPga9fO+z916/M4PuvnOVnb5cSYTbR73bT1t3PDz42LHwG1KxayFsIIWYjl1vz5WeO0ud0Y4kK46n9lWwpTPVV1ixOt4zxCZBmieLI/9mBAsLMJh76/XFeOl7Lv9yygpiI4IViaWsghBBj+NnbpXxQ1sI3b1rB7euzebW4nhZbH+cbBiprxiPcbCLMk6f/6Ppsuh0uXjkZ3Pp5CfJCCHERhyvb+OGuc9y0OpPb12dz18aF9LuMlgXnGm2+ypqJ2rAoiZzkaJ4/XBOEqx4gQV4IIS7iNx9UkBAdzrdvXYlSisXpFtYvSuLp/VWcb+iiKG3sVM1ITCbFbWuzea+0mdr2ngBf9aDzBO2ThRAiBByrbmf9oiQsUeG+bXdtXEhZs526jl6KBlXWTNRt67LQGv54NHijeQnyQggxis7efsqa7KzOTvDb/uFVC7BEGQ9LJzuSB1iUEsuluUk8d6g6aC0PJMgLIcQoTlZ3AHBJdqLf9ugIM7euzQJg8RRG8gC3rcumtMnOcc+5Ak2CvBBCjOKYL8gnDNv399cu5v/dfgmLUmKndI4bVy0gIszE857e84EmdfJCCDGK49XtLEqJITEmYti+5NgI7tiQM+VzJESH84u7NwxLCQWKBHkhhBjF8eoO1i0K/qJ22xcHrwOvpGuEEGIETV191LT3BG2EPV0kyAsh5q0z9Z0cq2ofcd/xamP70Ieuc42ka4QQ89aDzx6juLaTj1yygIdvXEZW4kB7gmPVHZgUrMyKn8ErnDoZyQsh5iWtNRUt3eRbY9l1qoFrfvAWv3rvgm//8ep2itIsQW0eNh0kyAsh5qWOnn5sfU4+sXEhux/czub8FP75pVO8fqoBrTXHqztYNcfz8SBBXggxT1W3Gf1ispNiyE6K4aefWs/KrHi+/OxRPihtodXumPMPXUGCvBBinqpu6wbwtQmOCjfz00+uRwGf//VBYO4/dAUJ8kKIEPLC0Ro6uvvHdax3JJ+TFOPblpMcw4/vXIPd4SLcrFi6YPJ9aWYLCfJCiJBQ1drNA08f5akDleM+3hIZRny0/4PVq5em842PLOfuzblEhpmDcanTam4/NhZChCxbn5Pdpxu4aXUmSo29KEdJk7Heqnfd1bFUt/WQlRQ94mfftzVvYhc7i8lIXggxKz21r5IHnj7K+XEG7dLGiQf5nOSYsQ+c4yTICyFmpSNVbQCcrusc1/FlzXYASptsw3qzP7bnAocr23yvtdZUt3WPe23WuUyCvBBiVjpa2Q7A2fqucR1f5knXdPU6aerq82239zn51z+f4tF3y3zb2rv7sTtcZCfJSF4IIaZdY2cvtR29wPiDfGmTncyEKGAgPw9wqq4TreFY1cCiHFVDyidDmQR5IcSsc8TTNCwnOZoz4wjyXb39NHX1cd2KDGAgPw9wssYI7jXtPb4R/sBEKAny46aUMiuljiil/uR5naeU2qeUKlFKPaOUGt51XwghRnC0qp0wk+LWtdnUtPfQ1Xvx2veyJiMfv7kghdgIM6We1wAnagZG8N7OkgMToSRdMxEPAKcHvf4e8COtdSHQBtwXwHMJIULY0cp2li2I97UVONdw8dF8qSc9U2CNoyAtzq/Cprimk8vykjGpgeX8qtt6iI8KIyE6PEh3MHsEJMgrpbKBDwOPel4r4Grg955DdgK3BOJcQojZ7f5fH+R7r5yZ9Ptdbs3x6nbW5CSyJMOYcTpWyqasyY7ZpFiYHEOhNc4X9HscLs43dnFZXjKL0y2+3vHVbT3zYhQPgRvJ/xj4KuD2vE4B2rXWTs/raiBrpDcqpe5XSh1USh1samoK0OUIIWaC0+XmrbNNfFDaMunPKGm0YXe4WJOTSFZiNJbIsDEfvpY121iYHENEmImCtDjqOnqx9Tk5Xd+JW8OKrARWZydyvLodrTVVrfOjfBICEOSVUh8BGrXWhybzfq31I1rrDVrrDVZr8NY5FEIE34VmOw6X25fznoyjnvr4NQsTUUqxOMMy5ki+tNFOgTUWMFI2YJRUFnvy8SuzErgkJ4G27n6qWntkJD9BW4CblFLlwNMYaZr/ABKVUt62CdlATQDOJYSYxU55Ji412xx0O5xjHD2yo1XtJESHk5diBO0lGRbO1ncNm+Dk5XJrLrTYyfcE98I0430ljTZO1HSQHBtBZkIUqz0dJd8820hPv4ucZBnJj4vW+mGtdbbWOhe4E3hDa/1J4E3gds9h9wAvTPVcQojZbfCI21umOFFHKttZnZOIyWT0lFmaYaGjp5+Gzr4Rj69t78HhdJOfagT3RSmxhJkUpU02TtZ0siIzHqUUSzIsRIaZ+PPxOmB+VNZAcOvkHwL+QSlVgpGjfyyI5xJCBIjT5eZnb5di75v4SPx0XSdhnuBc1TrxlI29z8m5hi7W5CT6ti1J9z58Hbm9ga+yJs0YyYebTSxMieFUbSfnGrpYlZXg274iM54DFa3A/KiRhwAHea31W1rrj3i+L9Nab9RaF2qt79Baj/xrWAgxqxwob+O7L5/h1eL6Cb/3TF0XmwtSgIsHeYfTTecIte8najpwa1g7KMgvzTAW0h7t4au3Jt47kgcotMaxp6QZp1uzMmtgdafVOYl4sz5ZEuSFEPORd8r/uYbxdXP0arM7qO/sZWthKlHhJqoukq75pz+cYMt33uDIoKZhDqebR98tw6SMYOyVEBNORnzUqEG+rMlGQnQ4ybED8y0L0+LodxnRfNXgIO/JyydEhxMfFfo18iBBXggxhDeXXtI4vp4xXqc96ZRlC+LJTooZdSTf2NnLH4/UYHM4ufux/RytasfhdPM3vz3M66cb+eebVvgFbDAevp4eNcjbybfG+vWF91bYJESH+6VlvL885stDV5AgL4QYonqSI/kzdUYQXrYgnpyk6FFH8k/tr8Lp1jx532UkxUbw6Uf38dlf7ee1Uw186+YV3L05d9h7lmZYKG200e9yD9tX2mTzBXWvQk9+fmVWvF/wz02JMQJ/4vx46AoS5IUQQ1S3GsG5qq2bHodr3O87XddJalwEVkskOckxVLd2Dyt77He5eXJfBdsXW7m8MJWn799EUmwE75W0jBrgwfjF4XC5+afnT1DRYuTgtdacqO6gsauPfGus3/H51lhMClZlJfptV0rx359Yy5d3LB73fc11svyfECHiXEMXzV19XF6YOqXPqW4z1j7t6nNS2mTze3B5MWfqu3wPSXOSYujqc9LZ4yQhZiD3/WpxPY1dfXz3o4sAyEyM5vm/vpwLzXYuzU0e9bNvWJXBoYpFPHOwiucOV7OtyEpJo42a9h5MCtYtTPI73hIVzhP3XcayBfHDPmtb0fyadCkjeSFCxA9eO8tXnzs+pc9wON3Ud/aybbHxi+L8OPPyTpebcw1dLFtglDt6c95VQ2a+/vr9ChYmx7B9cZpvW2pc5EUDPEBkmJl/uWUle756FZ/flk9Jo40VmfF876Or2PvwNWzKTxn2nssLU0mKlea3MpIXIkTUtPfQbOtDaz2uha9HUtfRg1vD1kIru041jDsvX95ip8/p9o3kvRONqlq7fX8JnKrtZH95K1+/cRlm0+SuLy0+iodvXMbDNy6b1PvnIxnJCxEiatt76e130z2BPPpQ3sqa3NQY8lJjOT/OIH960ENXwLdA9uCR/BP7KogKN3HHhuxJX5+YOAnyQoSAHoeLVrsDgBabY9Kf462syUmKoSjNMu50zZl6Y6ZrgadvjFGHHkaV5yGuy6159WQ9O5ZnkBgjKZTpJEFeiBBQ2zFQrthiH9/k8mNV7ez44du+Xw5gjOTNJsWChCiK0uOobO2mt3/svwxO13VRmBZHZJjZty0nOcY3kj9Y3kqL3cGHVqSP95ZEgEiQFyIE1LYPCvLjHMm/WlzP+UYbe8sGer9XtXaTER9FmNlEUZoFrfFbZWkkRyrbOHChdVglS86gCVGvFjcQEWbiyiVpI32ECCIJ8kKEgMFBfvDI/GKOVLYDcLhioLWA0WfdqIxZnG5MKLpYkH/haA0ff2QvSbER/N01RX77cpKjqW7rwe3WvFpcz9bCVOIipdZjukmQFyIE1LT3+r5vHke6xuXWHPMsan240j/Iex+a5qYaLXtHW1/1R7vO8cDTR1mbk8gfv7SFvFT/CUk5yTH0Od28fb6JmvYeSdXMEPm1KkQIqG3vISM+is7e/nGla841dNHtcJGZEMXJmk76nEbevaGr1zeSDzebjAqbEUbyb5xp4D92n+f29dn8262riAgbPl7M8ZRRPvbuBUwKrl0mQX4myEheiBBQ19FDZmIUKXER40rXeFM1n9mSi8Pl5mRNJ7XtvWg9EJwBFqdbOD9kJN/V28/X/3CSJemWUQM8DEyI2lPSzIbcZFLiIid5d2IqJMgLEQJq23vJTIwmOTaSZtvY6ZojlW0kx0Zwy9oswMjLe8snB3dtLEwbXmHz/VfOUt/Zy3c/OnqANz5n4JfFh1ZkTPieRGBIkBdijtNaU9PeQ1ZiNKmxEeNK1xypamdtTiJplihykqM5XNnmq2nPTvYfybv1wOpL+y+08pu9FXz28jzWDukXM1RUuBmrxRi9X7dcUjUzRXLyQsxxLXYHDqebzMRoWu0OimtHXibPq6Onn5JGG7esyQSM5l57y1p8D1oz4qN8xy7JMCpsbvnJe6RZorA7nGQnRfOVD42vi2OhNY4FCVG+h7li+kmQF2KO85ZPZiZGU9fRS4v94v1rjlW1A/hG4usXJfHC0Vr2X2glMzHar69MgTWOH9yxmpImGw0dvTTbHfzNVYXERIwvdPzo42smf2MiICTICzHHDQT5KCpaIuh3abr6nKMub3eksh2l4JJso3GYt03voYo2Li/w7+aolOKj6yffayYjIWrsg0RQSU5eiDnOWyOflRjtWzbvYnn5I1VtLE6zYPH8EliaYSE63GhHMLiyRoQGCfJCzHG17T3ERJhJiA73lSm2jFJho7XmSGU7axcm+raFmU2szjFG9YMra0RokCAvxBxX295DZmI0SilSvCP5UWrlLzTb6ejp9wvyMJCyyZ5HC1zPFxLkhZjjvEEeICXu4umaA+WtAMPKHy8vMFaCKkqzBOsyxQyZcpBXSkUppfYrpY4ppYqVUt/0bM9TSu1TSpUopZ5RSkkTaSGCoKa9l6xE4wGnNyffOkL/Gq01O9+vIN8aS6E1zm/f1qJUdj+4fdzruYq5IxAj+T7gaq31amANcL1SahPwPeBHWutCoA24LwDnEkIM0tvvotnWx4IEYyQfGWbGEhlG8wgj+XfON3OqrpMvXlGAaYTl9wqGBH4RGqYc5LXB28Eo3PNPA1cDv/ds3wncMtVzCSH81XcYlTXedA0YKZuRcvI/fauEjPgobl6bOW3XJ2ZeQHLySimzUuoo0AjsAkqBdq2103NINZA1ynvvV0odVEodbGpqCsTlCDFvDK6R90qOjRiWrjlS2cbeslY+ty3Pb/UmEfoCEuS11i6t9RogG9gILJ3Aex/RWm/QWm+wWq2BuBwh5o0aT5DP8hvJRw578Pqzt0tJiA7nzo0Lp/X6xMwLaHWN1rodeBPYDCQqpbwzarOBmkCeSwhhdJ8E/5mlqUPSNSWNXbxa3MA9mxfJykzzUCCqa6xKqUTP99HADuA0RrC/3XPYPcALUz2XEMJfXUcPVkukXwrGSNc4cLs1AM8cqCLCbOKey3Nn6CrFTArESH4B8KZS6jhwANiltf4T8BDwD0qpEiAFeCwA5xJCDFLWZB82SzUlNhKXW9PR0w/AO+eauTQvSRbtmKem/Leb1vo4sHaE7WUY+XkhRBDY+pwcrmzj81fk+233TYiyO3C43Jxt6OK2deN+TCZCjCTohJhFdr5fTphZcdWSNL+yyJG8X9KM063Zvti/YCEldqB/zVFPW+FtRVLUMF9JkBdiGtR19PCNPxbz5R1FrMgceVZpR3c///fFYt/rpRkWvnnTCi7LTxnx+LfPNREXGebrO+M1eCT/7vkmUuMiWZoh7QrmK+ldI8Q0+NV75bx+uoF7f3WAuo6eEY+p8qyx+vANS3n4hqV09PTzzZdOobUedqzWmrfPNXF5QcqwdVa9TcqabX3sOd/MtqLUEWe4ivlBgrwQQdbndPG7Q9WsyUnE3ufis788QFdv/7DjqtuM4L+lMJUvbC/gb68u4lRdJ4cq2oYdW9Zsp7qth+1LhqdhkjxB/t3zzbTYHWwrSg3wHYm5RIK8EEH2anEDrXYHX96xmJ9+ah0ljTa+9Nsj9LvcfscNndh0y9pM4qPC+NX75cM+8+2zxuzwK0bItYebTSREh/PW2UYAthZKkJ/PJMgLEWS/3VdBTnI02wpT2VZk5Vs3r+Sdc028fqrB77jqtm5iI8wkxhgrNsVEhPGxDTm8crKehs5ev2PfPtdEvjV21AWyU+KMZQCXZlhIi5cl+OYzCfJCBFFpk429Za3ceelCX178tnVZKAXnGmx+x1a39ZCdFOO3APenNy/CpTVP7qv0bevtd7G3rGVYVc1gqZ4KmysucoyYHyTICxFET+2rJMykuGPDwGLYUeFmshKjudDsH+Rr2nrIGjKxaVFKLFctSeO3+ypxOI30zr4LrfQ53RcN8t6+8pKPFxLkhQiS3n4Xvz9czXUr0kmz+KdM8lJjKWu2+22rbusecY3Vey7PpdnWx3dePs1zh6p5al8lkWEmNo1SWglGL5vocDOX5iYH5mbEnCV18kJMUb/LTX1H77D8+BtnGmnv7ueuETo/5qfG8vzhGrTWKKXo7O2ns9fp103Sa1thKssWxPPL98p9265dlkZU+Ogtg790VSEfXZd90WPE/CBBXogp+u83Svjp26W899DVWC0D/WFeP9VAYkw4m0cYceelxtLV56TZ5sBqiaTGUz6ZnTT8QarJpHjpb7bQ3tNPd58LW5+ThSkjP3D1sloi/a5FzF+SrhFiClxuzTMHqnA43bx2qt5v+5tnG7lqSRph5uH/m+V7ltq74EnZVPuC/MitDMLMJlLjIlmYEsPyzHhpGSzGTYK8EFPw7vkm6jt7CTMpXjk5EOSPVrXR1t3P1UvTRnxfXmosAGVNxsPXGs9s16EPXoWYKgnyQkzB7w5WkxwbwWe35PJ+aQttnsU6dp9uxGxSo5YwZiZGExFm8hvJR4WbfC0JhAgUCfJCTFKb3cGuUw3cvCaTm9dk4XJrdnkmOO0+3ciluUkkRIeP+F6zSZGbEuOrsKlu6yErMdqvRl6IQJAgL8QkvXC0BofLzR3rc1iRGU9OcjR/OVlHVWs3Zxu6uHZZ+kXfn5ca6xvJ17T3jPjQVYipkiAvxCQ9e7CaVVkJLM+MRynFjSsX8F5JMy8cNZYzHi0f75VvjaOixY7T5R61Rl6IqZIgL8QknKzp4FRdp99M1utXZtDv0vzkzVLyUmN9FTSjyUuNpd+lOddgo627Xx66iqCQOiwhhujtdwEMm0jU73Lz1tkmXjxWy+unGogKN3HT6kzf/jU5iWQmRFHb0TvmKB6MCVEAe0qMjpKSrhHBIEFeiCH+6olD9Dnd/Pbzm/y2/58XinlqfyVJMeHcui6LT2xcSGLMQDWMUorrVy7g8fcucM2ysYO8t4zy3fPNACPOdhViqiTIi5DW73Lz1d8f5zOX57I6J3HM491uzf4LrdgdLk7WdLAyy1iqr6mrj+cOVXPH+mz+7bZVhI8wwQngc9vyiI00s3EcPWOSYyNIiA5n34VWAHIkXSOCQHLyIqQduNDKH47U8NtBrXovprzFjt1hpGsGL9bx1P5KHC43X7yyYNQAD0b9+4PXLRlxlutQSinyUmNxON1EhBkzWoUINAnyIqS95qlb31PSPGyt1IefP86Lx2r9thXXdgKwYVESLx6rpcXWh8Pp5om9FWxfbKVgjIepE+XNy2clRss6rCIoJMiLkKW15vXTDUSYTdS09/hq0sHoGfPU/iqe2Fvh956TtR1EmE186+aVOJxunj5Qxcsn62js6uMzW3IDfo3evLyUT4pgmXKQV0rlKKXeVEqdUkoVK6Ue8GxPVkrtUkqd93xNmvrlCjF+Zxu6qG7r4bNbcwFjNO/18sk6AI5WtfuqaQCKazpZkmFheWY8WwtTeWJvBY+/V05eaizbR1hPdaq8ZZby0FUESyBG8k7gQa31cmAT8CWl1HLga8BurXURsNvzWohp411D9b4teWQnRfuqWABeOVlPRJgJh9PN0ap2wBj5F9d2sDIrHjAW66jr6OVYVTt3b14UlHSKjORFsE05yGut67TWhz3fdwGngSzgZmCn57CdwC1TPZcIDfsvtPLcoWqcLveUP2tonn2wXacbWZ2TSFp8FNuKUtlb2uKbXXq8uoN7t+ShFOwrM6pbajt6aevuZ3mmUVFz9dI0cpKjiY0wc/v67FHPMxWFaXH8r9WZXDNGCwQhJiugOXmlVC6wFtgHpGut6zy76gH5KRYA/PurZ3nwd8f4yH/t4YPSlkl/zjMHKtn8nTc4Xdc5bF9jpzEC3+GpV99aaKWrz8mx6nZfS+A7L81haUY8+8uNazhZ0wHAykxjJG82KX788bX85JPrsESN3GhsqiLCTPzXXWtZtiA+KJ8vRMCCvFIqDngO+Huttd//ddoYbo045FJK3a+UOqiUOtjU1BSoyxGzWLOtj8XpcXT1OrnrF3v55xeLJ/wZj+25wEPPnaC+s5f/3H1+2P7dZxoB2LE8A4DLC1JQyph49MrJepZmWMhNjeWyvGQOVbThcLopru3EpGBpxkDAXb8oiSuXjD2xSYjZKiBBXikVjhHgn9RaP+/Z3KCUWuDZvwBoHOm9WutHtNYbtNYbrNbAP9gSs0+zrY/N+SnsfnA7t6/PZucH5VS1do/rvVpr/mv3ef7lT6e4YWUGX9iez8sn6znX0OV33OunGshJjmZxuvFgMyk2glVZCbx4rJZDlW3csHIBAJvyk+ntd3Oipp3img4K0+KIjpB1UUXoCER1jQIeA05rrX84aNeLwD2e7+8BXpjqucTc53C66ex1khIXSVS4mS/vWIwCnj1YNa73v366kR/sOsdt67L4r7vW8sUrCoiJMPOTN0t8xzTb+thT0sy1y9L9+rNvLUylrMmO1nDDKmOEvzHPWH91b1krJ2s7WOnJxwsRKgIxkt8CfBq4Wil11PPvRuC7wA6l1HngWs9rMc+1elZOSokzer5kJUazfbGVZw9WjetB7BtnGrBEhfH9j15CmNlEUmwEn960iJeO1XKh2U5jVy93PbIXpeCO9Tl+791alApAvjWWojRjhJ8cG8Hi9Dj+cqKOhs4+lmdKblyElkBU1+zRWiut9SVa6zWef3/RWrdora/RWhdpra/VWrcG4oLF3NZs6wMgJXZgCv9dGxfS0NnHm2fHfibzQWkLl+Ul+7UNuG9bHuFmE9/5y2nufGQvNe09/PIzG4cF7PWLkkiNi+Cj67L9RviX5aX4Zrp6e9UIESpkxquYVi2ekXxq3ED3xquXppFmieSp/RfvL1Pb3kN5SzebC1L9tqdZorhr40JeO9VAQ0cvO+/dyOaClGHvjwwz885Xr+Kvthf4bb8sf6CZmIzkRaiRLpRiWrV4R/KDmnGFmU18bEMO//NWCbXtPWSOMvvTW255+QgB/K+vLKC2vYcvbC9g/aLRJ1fHRAz/kd+YZwT5RSkxxAepVFKImSIjeTGtWmz+OXmvj1+ag+biD2DfL20hKSacJemWYfvS4qN45O4NFw3wo0mzRLEiM35c7YGFmGtkJC+mVbO9jwizCUuk/49eTnIM24qs/OaDCj5ySSaFaf7dHrXW7C1rYXNBSlDaCzx9/6aLthAWYq6Sn2oxrVpsDlLiIvwefHp948PLUErx8Z9/QHFth9++ytZuatp72Jw/PFUTCJao8GHL/QkRCiTIi2nVYusblqrxKkq38OwXNhEZZuKuR/ZyuLLNt8+bjx/60FUIcXES5MWE9fa7+M7Lp8c9S3WwFrvDr3xyqHxrHM9+cTPJsRF86tF97Cszgvv7pS1YLZEUWGMnfd1CzEcS5MWEaK35+h9O8vO3y3jpeO3YbxjCm665mOykGJ79wmYWJETxmV8e4IPSFj4oa/H0n5HVk4SYCAnyYkKe3FfJc4erAShrso9xtD+tNS32vnGtZZoWH8VT928iOymaux/fR1NXX9Dy8UKEMgnyYtwOV7bxzZeKuXKJlY25yZQ12Sb0/m6Hi95+NymxFx/Je6VZjECfn2pU2lwu+XghJkyCvBiXjp5+/vqJw2QkRPHjj6+hIC3Ob83U8RiokR97JO+VGhfJs1/YzO++uJmFKTETOp8QQoL8rNTvcnP34/s5WD572v08ua/C6N1+51oSYyIosMbS1t1Pm6dNwXg0272zXcc3kvdKiAnnUpmoJMSkSJCfheo7ennnXBOvnx6xBf+kvVZcz8PPn7joknkj6XO6+NV75WwtTGXtQmNGab6nyqWsefwpG+9IPvUi1TVCiMCSID8LNXn6u1yYQAAdS2+/i2+8cJKn9leyt2xifyG8eLSWxq4+Pn9Fvm+bN09eOoGHrwN9ayY2khdCTJ4E+VnIO+IdKef96w/KOVTRNmz7WH79QTkNnX1Eh5v55XsXxv0+rTWPvnuBJekWrigaePCZnRRNuFlNqMLG24EyeZwPXoUQUydBfhbyjnjLW7pxuwdSKw6nm2+9dGpCQRqgs7ef/3mrlO2Lrdy7NZddpxvGPZHpnfPNnG3o4vNX5PvVqIeZTSxMjrnoXxuvFtdzpn5gud9mWx9xkWHSPkCIaSRBfhbyLqzhcLqp7ejxbS9rtuF0a843TCyN8+i7F2jv7ucfP7SET2/KxawUO98vH9d7f/FOGenxkdy0OnPYvnxr3Kgj+fdKmvniE4f4/itnfdvGMxFKCBFYEuRnoWbbQMXK4JTN2XpjseqyZtu4lsoD46+Cx94t48OrFrAyK4GMhChuWLWAZw5WYe9z4nS5+eFrZ7n639/ytRDwevFYLXtKmvnM5XlEhA3/Ucm3xlLR0o3L7f8gt9Xu4B+ePYrWcKSyzfegt8XeN+4aeSFEYEiQn4W8aQ3wD/JnPEG+36WpGGe65fH3LtDT7+Ifrlvs2/bZLbl09Tr5n7dKuOsXe/nPN0po7Xbwqcf28buDVWit+fHr5/i7p46wYVESn968aMTPLkiNw+FyU902cC1aax567jitdgef3rSItu5+3z0YI3mprBFiOkmQn4VabA6WZFiIiTD7pUPO1Xf5RtTnG7rG9VnHqztYlZVAgXWgP/u6hUmszknkJ2+Wcqq2kx99fDVvf+UqLstL4R9/f5wP/+cefvz6eW5bl8WTn7/M9wtnKF8Z5aBrfHJfJbtONfDQ9Uv51Cbjl8PhynbA+AslVdI1QkwrCfKzUIu9D2tcJLkpsZS3+I/kryiyAow7L1/eYic3dXjnxq9dv5Trlqfzp7/bxq1rs0mICeeXn72UT21ayOn6Tr52w1J+cMdqIsNGf0ia7/nFUeYZqbd3O/i3v5xmW1Eq927JoygtDktkGIcr23C7Na32vot2oBRCBJ6sDDULNdscXJobgdmsOFljLJ7R1dtPTXsPn7hsIWfqOznXOHaQdzjd1LT1cOva7GH7NhekDFvsOtxs4l9vWcVD1y/FMo61TpNiwkmIDvf1sHlyXyXdDhcP37DMt3rTmoWJHK5oo72nH7eWGnkhppuM5GcZp8tNW7eRu85PjaWqtRuH0805z8h9SbqForS4caVrqtu6cWvInWDPl/EEeAClFPnWWMqa7Dicbna+b8yKXZ4Z7ztm3cIkzjV0+f4ikZy8ENNLgvws09rtQGuwxkWQlxqLW0NVW7evsmZJhoXF6RbKmu1jVthUtBgPRBcFsbFXfmocZc02XjpmzIr93LY8v/3rFiXh1vCGp0VDqlTXCDGtJMjPMoM7NXpz6Rea7Jxr6CI2wkxWYjSFaXE4nG4qx6iw8Y6eF6UEbzWlfGssDZ19/OStEorS4ti+2Oq3f01OIgC7TjUAMpIXYroFJMgrpR5XSjUqpU4O2paslNqllDrv+ZoUiHOFOl8TL0+6BowyyjP1nSzOsGAyKRanWwA4P0ZevqKlm7jIsKDWphcMqrD53La8YSs3JUSHU5QWx1lPekly8kJMr0CN5H8FXD9k29eA3VrrImC357UYQ/OgJl6JMREkxYRT1mznbH0XSzzBvSDNqGoZKy9f3mJnUUpMUJfM81bYpMZFcPOarBGPWefpXKkUJMVIkBdiOgUkyGut3wGGtja8Gdjp+X4ncEsgzhXqvEHe2443LzWWA+WttHX3syTDCPJxkWFkJUaPaySfG8RUDRj5/sSYcD6/LX/UnjTrFxlBPjkmArNJ1mgVYjoFMyefrrWu83xfD6SPdJBS6n6l1EGl1MGmpqYgXs7c0GxzEG5WxEcb1a25qbGUNA5U1ngVpcddtFbe6XJT1dod1IeuAJFhZj742jXcP6gN8VDrFiUCkqoRYiZMy4NXbTQvGXGlCq31I1rrDVrrDVardaRD5pUWmzFhyJtiyR80kck7kgcoSoujtMmGy63RWvOtl07x07dKfftr23txunXQR/IA0RHmi6aE8lPjiI8Kk4lQQsyAYE6GalBKLdBa1ymlFgCBXeYoRLXYHaRaBka8eanenHekX2VKUZqFPqcxWv/ziToef+8CSTHh3H9FPmaTGlRZM/ProppMii/vWCx95IWYAcEcyb8I3OP5/h7ghSCeK2Q02/yn/ud5RvJLMuL8jitKN17/7O1S/v21syxMjqGtu58Tnhmy3gZmI7U0mAmf3ZI36oNZIUTwBKqE8ingA2CJUqpaKXUf8F1gh1LqPHCt57UYQ4vNQeqgEXtuagxKwdKMeL/jCj0VNk8fqGJlZgLPfGETSsFbZ40/mCqa7USFm0izSIpEiPksIOkarfVdo+y6JhCfP19orWmy9fl1aoyJCOPRuzewMivB71hLVDhZidE4XG5+cfcGMhKiWJ2dyFtnm/j7axdT7qmsCWb5pBBi9pMGZbOIrc+Jw+keVoVyzbIRC5P46afWER8VTkZCFABXLrHyH7vP02p3UNFi97UCFkLMX9LWYBYZPNt1PC7JTvTLuV+5JA2t4e1zjVS0Br9GXggx+0mQn0UGZrtOLo9+SVYCybERPL2/CofTHdSeNUKIuUGCfJC02R28crJu7AMHafaN5CdXamgyKbYVpbLvgjH5eKIthoUQoUeCfJA8daCSLz5xmJr2nnG/x9fSYAqdGq9cMjChbNEsKZ8UQswcCfJBUuWpU/eu7DQe3pz8VCYNXVFkRSmIMJvIiI+a9OcIIUKDBPkgqW4zRvDFtZ3jfk+LvY/EmHDCzZP/z5ISF8kl2YksSomRZmBCCCmhDBbvSL54AiN5Y7br1Kf+//vtl9DT75ry5wgh5j4J8kHgdmtfLn4iI/nmIbNdJ6toULdKIcT8JumaIGjs6qPfpclLjaW+s9f3QHUszba+gAR5IYTwkiAfBNVtRqrmQysygPGP5ltsDum5LoQIKAnyk2C0xx+d96HrdSuMdgTjqbBxON109PTLSF4IEVAS5CfhOy+f4bofve0bsQ/l3b58QTwLk2Morh07yDcNWttVCCECRYL8JLxxppFzDTY+/vO9lDfbh+2vau0hNS6SqHAzKzLjx5Wu+d3BKmBg0WshhAgECfIT1O1wUtpk48ZVGfT0u/jYzz+gpLHL75jq9m6yk6IBWJmVQEVLN529/aN+ZkdPP4/vucCO5eksWxA/6nFCCDFREuQn6HRdF1rDrWuzefr+TWjgvp0H/fL01W095CQbfWOWZxpB+9RFRvO/eq+czl4nD1xTFNRrF0LMPxLkJ8ibX1+RGc/idAsPXFNERUs3Va3Gw1aXW1Pb3jMwks80FvsY7eFrZ28/j+0p49pl6cMWBhFCiKmSID9BxTWdJMWEs8CzUMeluckAHCg3Oj82dvXS79K+IG+1RJJmiRx1JL9TRvFCiCCSID9BJ2s7WJmV4FtWrygtjvioMA5WGEHeWz6ZnTTQ5ndlVgInB1XY9Pa7OFrVzm/2VvDongtcuyyNVdkyihdCBJ60NZgAh9PNuYYu7t2a59tmMik25CZzoLwNGOhZ4x3Jg5HaeetsI3c9spfq9m5q23txuY0cfkZ8FF/50JJpvAshxHwiQX4CzjV00e/Svjy714bcJN4400ir3eEbyWclDgT5q5em8eKxWvqcLtbmJHHLmhhWZMazMiuBrMRoWWxbCBE08ybI9zldHK5oZ/2iJCLCJpel8ubVV2T6lzl68/KHKtqobusmzWLUyHutXZjE2/941SSvXAghJi/kg/zpuk6eOVDFH4/W0N7dz4M7FvO3gx5yaq15+PkTxEWGcffmXBZeZMm8k7UdxEaYhy2QvSorgQiziYPlrVS39filaoQQYiaFdJC/0Gznxv98l3CzieuWp1PV1sMT+yr44pUFvoU53jrXxNMHjNmmj713gWuXpfO/P7xsxEWwi2s7WZ4Zj2nIYhxR4WYuyU7gQHkrzTYHa3ISg35vQggxHkGvrlFKXa+UOquUKlFKfS3Y5xusvNmO1vCbezfy359YxwPXFNLQ2cerxfW+Yx55u4yM+Cje/epVfOnKQt4raeb7r5wd9lkut+ZUbScrMkeugtmQm8yJmg6/GnkhhJhpQQ3ySikz8BPgBmA5cJdSankwzzlYXUcvgG/26fbFaSxMjmHn++UAHK9u54OyFu7dmktOcgxf+dASrlqSxvGa9mGfdaHZTk+/a1g+3mtjXhL9Lo3Trf3KJ4UQYiYFeyS/ESjRWpdprR3A08DNQT6nT31HDyZlTEgCMJsUd29exIHyNoprO/j5O2VYIsO4a+NC33tWZMVT1dpDR7d/rxnvTNfRZqWuX5js+15G8kKI2SLYQT4LqBr0utqzzUcpdb9S6qBS6mBTU1NAT17f2YvVEum3MPYd63OICjfx3ZfP8PKJOj6xaSGWqHDffm955ND2wMW1nUSEmShMixvxXAkx4SzxLLvn/ctBCCFm2ozPeNVaP6K13qC13mC1WgP62XUdvWQk+I+qE2LCuXVtFu+eb8ZsUty7Jc9vv3ekfnJIkD9W1c7SDIvfL4yhNuQmYVL4Wh4IIcRMC3aQrwFyBr3O9mybFvUdvSyIHx5w77k8F4Cb12SRPmR/cmwEWYnRnKwZ6DXT7XByuLKNy/KSuZi/vbqIRz69wa9GXgghZlKwSygPAEVKqTyM4H4n8Ikgn9OnvqOXLYWpw7YvzYhn570bWT1Kv5gVmfF+XSP3XWil36XZVnTxvzQyEqLIkFG8EGIWCepIXmvtBP4GeBU4DTyrtS4O9HmOV7fz0O+P+z0stfU56epzjhp0ty+2khgz8lJ7K7MSKGu20+VZ6GPP+WYiwkxsHGMkL4QQs03Qc/Ja679orRdrrQu01t8OxjlabA6eOVhFSdPACk31nvLJyeTHV2YZZZKn64zP23O+mUtzkyQNI4SYc2b8wWsg5FuN2amljQPrrXqDfMYIOfmxDF7oo6Gzl7MNXWwtDOxDYSGEmA4h0dYgOymGCLOJ0iabb1tdh9ENcjI58rT4KNIskZys7SAh2iiv3FY0PLcvhBCzXUgEebNJkZca6xfkvSP5odUz47UyK4GTNR1oDSmxESyXBbaFEHNQSAR5gIK0WM7UDcrJd/aSHBsx6Tz6Ss9CH802B1sKU4c1JRNCiLkgJHLyAAXWOCpau3E43YAxkp9MPt5rRVYCbg2tdgfbRijDFEKIuSCkgrzLralsNR6+1nX0Tmnm6apBPWq2Sj5eCDFHhUyQ91bYlHgqbOo7e0mfQpBfkBBFcmwEBdZYMhOl4ZgQYm4KmZx8vtVoHFbaZKO330Wr3TFiS4PxUkrx4HWLSRplwpQQQswFIRPk4yLDyIiPorTJRmNnHzC58snBPnnZokBcmhBCzJiQSdeAUWFT2mT31cgvSJA0ixBifgutIG+No6zJ5lsRKiMhcoavSAghZlbIBfmuXicnPB0kh/aSF0KI+Sakgry3wua9kmYskWHERYbMIwchhJiUkAryBZ4KmzP1XdLXXQghCLEgnxEfRUyE0cZAgrwQQoRYkDeZlC9lM5WWBkIIESpCKsjDQMpGFtMWQogQDPL5qUaQl8oaIYQIwSBfkOZJ10iNvBBChF6Qv3JJGp/bmsem/JSZvhQhhJhxIVdIHhcZxv/+yPKZvgwhhJgVQm4kL4QQYoAEeSGECGES5IUQIoRJkBdCiBA2pSCvlLpDKVWslHIrpTYM2fewUqpEKXVWKfWhqV2mEEKIyZhqdc1J4Dbg54M3KqWWA3cCK4BM4HWl1GKttWuK5xNCCDEBUxrJa61Pa63PjrDrZuBprXWf1voCUAJsnMq5hBBCTFywcvJZQNWg19WebcMope5XSh1USh1samoK0uUIIcT8NGa6Rin1OpAxwq6va61fmOoFaK0fAR7xnKtJKVUxyY9KBZqnej1zjNzz/CD3PD9M5Z4XjbZjzCCvtb52EiesAXIGvc72bBvrXNZJnAsApdRBrfWGsY8MHXLP84Pc8/wQrHsOVrrmReBOpVSkUioPKAL2B+lcQgghRjHVEspblVLVwGbgz0qpVwG01sXAs8Ap4BXgS1JZI4QQ029KJZRa6z8Afxhl37eBb0/l8yfokWk812wh9zw/yD3PD0G5Z6W1DsbnCiGEmAWkrYEQQoQwCfJCCBHCQiLIK6Wu9/TIKVFKfW2mrycYlFKPK6UalVInB21LVkrtUkqd93xNmslrDDSlVI5S6k2l1ClPj6QHPNtD9r6VUlFKqf1KqWOee/6mZ3ueUmqf52f8GaVUxExfayAppcxKqSNKqT95Xof6/ZYrpU4opY4qpQ56tgXl53rOB3mllBn4CXADsBy4y9M7J9T8Crh+yLavAbu11kXAbs/rUOIEHtRaLwc2AV/y/LcN5fvuA67WWq8G1gDXK6U2Ad8DfqS1LgTagPtm7hKD4gHg9KDXoX6/AFdprdcMqo0Pys/1nA/yGD1xSrTWZVprB/A0Ru+ckKK1fgdoHbL5ZmCn5/udwC3TeU3BprWu01of9nzfhREEsgjh+9YGm+dluOefBq4Gfu/ZHlL3rJTKBj4MPOp5rQjh+72IoPxch0KQH3efnBCUrrWu83xfD6TP5MUEk1IqF1gL7CPE79uTujgKNAK7gFKgXWvt9BwSaj/jPwa+Crg9r1MI7fsF4xf3a0qpQ0qp+z3bgvJzHXILec9XWmutlArJelilVBzwHPD3WutOY6BnCMX79kwcXKOUSsSYh7J0Zq8oeJRSHwEatdaHlFJXzvDlTKetWusapVQasEspdWbwzkD+XIfCSH5SfXJCRINSagGA52vjDF9PwCmlwjEC/JNa6+c9m0P+vgG01u3AmxgzyhOVUt5BWSj9jG8BblJKlWOkWq8G/oPQvV8AtNY1nq+NGL/INxKkn+tQCPIHgCLP0/gIjMVKXpzha5ouLwL3eL6/B5hyV9DZxJObfQw4rbX+4aBdIXvfSimrZwSPUioa2IHxLOJN4HbPYSFzz1rrh7XW2VrrXIz/d9/QWn+SEL1fAKVUrFLK4v0euA5jAaag/FyHxIxXpdSNGHk9M/C4p6VCSFFKPQVcidGOtAH4v8AfMXoELQQqgI9prYc+nJ2zlFJbgXeBEwzka/8JIy8fkvetlLoE46GbGWMQ9qzW+ltKqXyMkW4ycAT4lNa6b+auNPA86ZqvaK0/Esr367k3bzuYMOC3WutvK6VSCMLPdUgEeSGEECMLhXSNEEKIUUiQF0KIECZBXgghQpgEeSGECGES5IUQIoRJkBdCiBAmQV4IIULY/wfp3iPkf5jG9QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -109,7 +120,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -162,26 +173,37 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_4275bea8cf61cb4b45f01fa01c73d194 NOW.\n" + "Building... This may take some time.\n", + "Done.\n", + "Sampling...\n", + " 0/8000 [>---------------------------] 0% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 Messages received during sampling:\n", + " Gradient evaluation took 4.2e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 4.2e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 0.00014 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 1.4 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 3.9e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.\n", + " Adjust your expectations accordingly!\n", + "\n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + "Done.\n" ] } ], - "source": [ - "sm = pystan.StanModel(model_code=refit_lr_code)" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], "source": [ "data_dict = {\n", " \"N\": len(ydata),\n", @@ -192,8 +214,9 @@ " \"x_ex\": [],\n", " \"y_ex\": [],\n", "}\n", - "sample_kwargs = {\"iter\": 1000, \"chains\": 4}\n", - "fit = sm.sampling(data=data_dict, **sample_kwargs)" + "sm = stan.build(program_code=refit_lr_code, data=data_dict)\n", + "sample_kwargs = {\"num_samples\": 1000, \"num_chains\": 4}\n", + "fit = sm.sample(**sample_kwargs)" ] }, { @@ -218,7 +241,7 @@ " \"log_likelihood\": [\"log_lik\", \"log_lik_ex\"],\n", " \"dims\": dims,\n", "}\n", - "idata = az.from_pystan(posterior=fit, **idata_kwargs)" + "idata = az.from_pystan(posterior=fit, posterior_model=sm, **idata_kwargs)" ] }, { @@ -243,10 +266,10 @@ " N_obs = len(mask)\n", " N_ex = np.sum(~mask)\n", " observations = {\n", - " \"N\": N_obs - N_ex,\n", + " \"N\": int(N_obs - N_ex),\n", " \"x\": xdata[mask],\n", " \"y\": ydata[mask],\n", - " \"N_ex\": N_ex,\n", + " \"N_ex\": int(N_ex),\n", " \"x_ex\": xdata[~mask],\n", " \"y_ex\": ydata[~mask],\n", " }\n", @@ -261,11 +284,11 @@ { "data": { "text/plain": [ - "Computed from 2000 by 100 log-likelihood matrix\n", + "Computed from 4000 by 100 log-likelihood matrix\n", "\n", " Estimate SE\n", - "elpd_loo -250.66 7.17\n", - "p_loo 2.85 -\n", + "elpd_loo -250.85 7.20\n", + "p_loo 3.05 -\n", "------\n", "\n", "Pareto k diagnostic values:\n", @@ -273,13 +296,7 @@ "(-Inf, 0.5] (good) 100 100.0%\n", " (0.5, 0.7] (ok) 0 0.0%\n", " (0.7, 1] (bad) 0 0.0%\n", - " (1, Inf) (very bad) 0 0.0%\n", - "\n", - "\n", - "The scale is now log by default. Use 'scale' argument or 'stats.ic_scale' rcParam if\n", - "you rely on a specific value.\n", - "A higher log-score (or a lower deviance) indicates a model with better predictive\n", - "accuracy." + " (1, Inf) (very bad) 0 0.0%" ] }, "execution_count": 9, @@ -322,7 +339,7 @@ "outputs": [], "source": [ "pystan_wrapper = LinearRegressionWrapper(\n", - " sm, idata_orig=idata, sample_kwargs=sample_kwargs, idata_kwargs=idata_kwargs\n", + " refit_lr_code, idata_orig=idata, sample_kwargs=sample_kwargs, idata_kwargs=idata_kwargs\n", ")" ] }, @@ -342,16 +359,93 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/percy/anaconda3/envs/arviz/lib/python3.6/site-packages/arviz/stats/stats_refitting.py:98: UserWarning: reloo is an experimental and untested feature\n", + "Building...\n", + "/home/ahartikainen/github_ubuntu/arviz/arviz/stats/stats_refitting.py:99: UserWarning: reloo is an experimental and untested feature\n", " warnings.warn(\"reloo is an experimental and untested feature\", UserWarning)\n", - "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 13\n", - "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 13\n", - "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 42\n", - "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 42\n", - "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 56\n", - "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 56\n", - "arviz.stats.stats_refitting - INFO - Refitting model excluding observation 73\n", - "INFO:arviz.stats.stats_refitting:Refitting model excluding observation 73\n" + "Building...\n", + "Found model in cache. Done.\n", + "Sampling...\n", + " 0/8000 [>---------------------------] 0% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 Messages received during sampling:\n", + " Gradient evaluation took 4e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.4 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 5.9e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 4.4e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 3.8e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.\n", + " Adjust your expectations accordingly!\n", + "\n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + "Done.\n", + "Building...\n", + "Found model in cache. Done.\n", + "Sampling...\n", + " 0/8000 [>---------------------------] 0% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 Messages received during sampling:\n", + " Gradient evaluation took 2e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 1.9e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 2.2e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 2.4e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.\n", + " Adjust your expectations accordingly!\n", + "\n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + "Done.\n", + "Building...\n", + "Found model in cache. Done.\n", + "Sampling...\n", + " 0/8000 [>---------------------------] 0% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 Messages received during sampling:\n", + " Gradient evaluation took 1.9e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 1.6e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 1.8e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 1.3e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.\n", + " Adjust your expectations accordingly!\n", + "\n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + "Done.\n", + "Building...\n", + "Found model in cache. Done.\n", + "Sampling...\n", + " 0/8000 [>---------------------------] 0% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + " 8000/8000 [============================] 100% 1 sec/0 Messages received during sampling:\n", + " Gradient evaluation took 1.7e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 1.8e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 2.2e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.\n", + " Adjust your expectations accordingly!\n", + " Gradient evaluation took 2.6e-05 seconds\n", + " 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.\n", + " Adjust your expectations accordingly!\n", + "\n", + " 8000/8000 [============================] 100% 1 sec/0 \n", + "Done.\n" ] } ], @@ -367,11 +461,11 @@ { "data": { "text/plain": [ - "Computed from 2000 by 100 log-likelihood matrix\n", + "Computed from 4000 by 100 log-likelihood matrix\n", "\n", " Estimate SE\n", - "elpd_loo -250.67 7.17\n", - "p_loo 2.86 -\n", + "elpd_loo -250.85 7.20\n", + "p_loo 3.05 -\n", "------\n", "\n", "Pareto k diagnostic values:\n", @@ -379,13 +473,7 @@ "(-Inf, 0.5] (good) 100 100.0%\n", " (0.5, 0.7] (ok) 0 0.0%\n", " (0.7, 1] (bad) 0 0.0%\n", - " (1, Inf) (very bad) 0 0.0%\n", - "\n", - "\n", - "The scale is now log by default. Use 'scale' argument or 'stats.ic_scale' rcParam if\n", - "you rely on a specific value.\n", - "A higher log-score (or a lower deviance) indicates a model with better predictive\n", - "accuracy." + " (1, Inf) (very bad) 0 0.0%" ] }, "execution_count": 13, @@ -405,11 +493,11 @@ { "data": { "text/plain": [ - "Computed from 2000 by 100 log-likelihood matrix\n", + "Computed from 4000 by 100 log-likelihood matrix\n", "\n", " Estimate SE\n", - "elpd_loo -250.66 7.17\n", - "p_loo 2.85 -\n", + "elpd_loo -250.85 7.20\n", + "p_loo 3.05 -\n", "------\n", "\n", "Pareto k diagnostic values:\n", @@ -417,13 +505,7 @@ "(-Inf, 0.5] (good) 96 96.0%\n", " (0.5, 0.7] (ok) 0 0.0%\n", " (0.7, 1] (bad) 2 2.0%\n", - " (1, Inf) (very bad) 2 2.0%\n", - "\n", - "\n", - "The scale is now log by default. Use 'scale' argument or 'stats.ic_scale' rcParam if\n", - "you rely on a specific value.\n", - "A higher log-score (or a lower deviance) indicates a model with better predictive\n", - "accuracy." + " (1, Inf) (very bad) 2 2.0%" ] }, "execution_count": 14, @@ -434,6 +516,13 @@ "source": [ "loo_orig" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -456,5 +545,5 @@ } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/doc/source/user_guide/sampling_wrappers.md b/doc/source/user_guide/sampling_wrappers.md index 1421cf34cc..9c665e97b4 100644 --- a/doc/source/user_guide/sampling_wrappers.md +++ b/doc/source/user_guide/sampling_wrappers.md @@ -11,10 +11,11 @@ whereas the second one externalizes the computation of the pointwise log likelihood to the user who is expected to write it with xarray/numpy. ```{toctree} +pystan2_refitting pystan_refitting pymc3_refitting numpyro_refitting -pystan_refitting_xr_lik +pystan2_refitting_xr_lik pymc3_refitting_xr_lik numpyro_refitting_xr_lik ```