{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" }, "colab": { "name": "Gaussian_Mixture_Model.ipynb", "provenance": [] } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "nTKnq86xzpLl" }, "source": [ "# Gaussian Mixture Model and the EM algorithm\n", "\n", "Code implementing the process of clustering via the Gaussian mixture model, using the expectation maximization (EM) algorithm (see chapter 10 of [Bishop's PRML](https://www.microsoft.com/en-us/research/people/cmbishop/prml-book/)). " ] }, { "cell_type": "code", "metadata": { "collapsed": true, "id": "kMhxOfOrzpLv" }, "source": [ "import numpy as np\n", "from scipy.optimize import minimize\n", "from matplotlib import pyplot as plt\n", "\n", "\n", "# Defining epsilon for limits\n", "eps = 1e-6\n", " \n", "# Configuring matplotlib\n", "import matplotlib.pyplot as plt\n", "plt.rcParams[\"figure.figsize\"] = (8,8)\n", "plt.rcParams['axes.labelsize'] = 14\n", "plt.rcParams['xtick.labelsize'] = 12\n", "plt.rcParams['ytick.labelsize'] = 12\n", "plt.rcParams['lines.linewidth'] = 3.0\n", "plt.style.use('dark_background')\n", "plt.rcParams[\"image.cmap\"] = 'spring'\n", "\n", "# Fix np.random seed for replicability\n", "np.random.seed(0)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "K3MOPa3SzpLw" }, "source": [ "# Setting and Model\n", "\n", "## Input\n", "* $N \\in \\mathbb{N}$: number of data points, and\n", "* $d \\in \\mathbb{N}$: dimension of each data point\n", "* Input data: $x_0, x_1, \\ldots , x_{N-1} \\in \\mathbb{R}^d$. \n", "* Input data represented as matrix $X$, where $X_{n,i}$ is $i$-th component of $x_n$.\n", "\n", "## Gaussian Mixture Model\n", "\n", "The pdf of $N$ data points drawn from a Gaussian mixture model with $K \\in \\mathbb{N}$ components is given by\n", "\n", "\\begin{align}\n", " p\\left(X,Z \\middle| \\mu, \\Sigma, \\pi \\right) = \n", " \\prod_{n=0}^{N-1} \\prod_{k=0}^{K-1} \n", " \\left[\\pi_k \\mathcal{N}\\left(x_n \\middle| \\mu_k , \\Sigma_k \\right) \\right]^{z_{n,k}}\n", "\\end{align}\n", "\n", "where \n", "* $z_{n,k}\\in \\{0,1 \\}, \\sum_{k=0}^{K-1}z_{n,k}=1$ denotes the latent cluster vector for data point $n$ \n", "* $\\pi_k \\geq 0, \\sum_{k=0}^{K-1} \\pi_k = 1$ are the cluster probabilities\n", "* $\\mathcal{N}\\left(x \\middle| \\mu_k , \\Sigma_k \\right)$ is a gaussian pdf with mean $\\mu_k \\in \\mathbb{R}^d$ and covariance matrix $\\Sigma_k$\n", "\n", "Marginalizing over $Z$, we get\n", "\n", "\n", "\\begin{align}\n", " p \\left( X \\middle| \\mu, \\Sigma, \\pi \\right) =\n", " \\prod_{n=0}^{N-1} \\left[ \\sum_{k=0}^{K-1} \\pi_k \\mathcal{N}\\left(x_n \\middle| \\mu_k, \\Sigma_k \\right) \\right]\n", "\\end{align}\n", "\n" ] }, { "cell_type": "code", "metadata": { "collapsed": true, "id": "MGzEyvvzzpL6" }, "source": [ "class GaussianMixtureModel:\n", " def __init__(self, K):\n", " self.K = K\n", " self.Pi = None\n", " self.Mu = None\n", " self.Sigma = None\n", " \n", " def _init_params(self, X, random_state=None):\n", " '''\n", " Method for initializing model parameterse based on the size and variance of the input data array. \n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " '''\n", " n_samples, n_features = np.shape(X)\n", " rnd = np.random.RandomState(seed=random_state)\n", " \n", " self.Pi = np.ones(self.K)/self.K\n", " self.Mu = X[rnd.choice(n_samples, size=self.K, replace=False)]\n", " self.Sigma = np.tile(np.diag(np.var(X, axis=0)), (self.K, 1, 1))\n", "\n", " \n", " def _calc_nmat(self, X):\n", " '''\n", " Method for calculating array corresponding $\\mathcal{N}(x_n | \\mu_k)$\n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " \n", " Returns\n", " ----------\n", " Nmat : 2D numpy array\n", " 2-D numpy array representing probability density for each sample and each component, \n", " where Nmat[n, k] = $\\mathcal{N}(x_n | \\mu_k)$.\n", " \n", " '''\n", " n_samples, n_features = np.shape(X)\n", " \n", " Diff = np.reshape(X, (n_samples, 1, n_features) ) - np.reshape(self.Mu, (1, self.K, n_features) )\n", " L = np.linalg.inv(self.Sigma)\n", " exponent = np.einsum(\"nkj,nkj->nk\", np.einsum(\"nki,kij->nkj\", Diff, L), Diff)\n", " Nmat = np.exp(-0.5*exponent)/np.sqrt(np.linalg.det(self.Sigma)) / (2*np.pi)**(n_features/2)\n", " return Nmat\n", " \n", " def _Estep(self, X):\n", " '''\n", " Method for calculating the array corresponding to responsibility.\n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " \n", " Returns\n", " ----------\n", " Gam : 2D numpy array\n", " 2-D numpy array representing responsibility of each component for each sample in X, \n", " where Gamt[n, k] = $\\gamma_{n, k}$.\n", " \n", " '''\n", " n_samples, n_features = np.shape(X)\n", " Nmat = self._calc_nmat(X)\n", " tmp = Nmat * self.Pi\n", " Gam = tmp/np.reshape(np.sum(tmp, axis=1), (n_samples, 1) )\n", " return Gam\n", " \n", " def _Mstep(self, X, Gam):\n", " '''\n", " Method for calculating the model parameters based on the responsibility gamma.\n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " Gam : 2D numpy array\n", " 2-D numpy array representing responsibility of each component for each sample in X, \n", " where Gamt[n, k] = $\\gamma_{n, k}$.\n", " '''\n", " n_samples, n_features = np.shape(X)\n", " Diff = np.reshape(X, (n_samples, 1, n_features) ) - np.reshape(self.Mu, (1, self.K, n_features) )\n", " Nk = np.sum(Gam, axis=0)\n", " self.Pi = Nk/n_samples\n", " self.Mu = Gam.T @ X / np.reshape(Nk, (self.K, 1))\n", " self.Sigma = np.einsum(\"nki,nkj->kij\", np.einsum(\"nk,nki->nki\", Gam, Diff), Diff) / np.reshape(Nk, (self.K, 1, 1))\n", " \n", " def calc_prob_density(self, X):\n", " '''\n", " Method for calculating the probablity density $\\sum_k \\pi_k \\mathcal{N}(x_n | \\mu_k)$\n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " \n", " Returns\n", " ----------\n", " prob_density : 2D numpy array\n", "\n", " '''\n", " prob_density = self._calc_nmat(X) @ self.Pi\n", " return prob_density\n", " \n", " \n", " def calc_log_likelihood(self, X):\n", " '''\n", " Method for calculating the log-likelihood for the input X and current model parameters.\n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " Returns\n", " ----------\n", " loglikelihood : float\n", " The log-likelihood of the input data X with respect to current parameter set.\n", " \n", " '''\n", " log_likelihood = np.sum(np.log(self.calc_prob_density(X)))\n", " return log_likelihood\n", " \n", " \n", " def fit(self, X, max_iter, tol, disp_message, random_state=None):\n", " '''\n", " Method for performing learning. \n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " max_iter : int\n", " Maximum number of iteration\n", " tol : float, positive\n", " Precision. If the change of parameter is below this value, the iteration is stopped\n", " disp_message : Boolean\n", " Whether or not to show the message about the number of iteration\n", " '''\n", " self._init_params(X, random_state=random_state)\n", " log_likelihood = - np.float(\"inf\")\n", " \n", " for i in range(max_iter):\n", " Gam = self._Estep(X)\n", " self._Mstep(X, Gam)\n", " log_likelihood_old = log_likelihood\n", " log_likelihood = self.calc_log_likelihood(X)\n", " \n", " if log_likelihood - log_likelihood_old < tol:\n", " break\n", " if disp_message:\n", " print(f\"n_iter : {i+1}\")\n", " print(f\"log_likelihood : {log_likelihood}\")\n", " \n", " def predict_proba(self, X):\n", " '''\n", " Method for calculating the array corresponding to responsibility. Just a different name for _Estep\n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " \n", " Returns\n", " ----------\n", " Gam : 2D numpy array\n", " 2-D numpy array representing responsibility of each component for each sample in X, \n", " where Gamt[n, k] = $\\gamma_{n, k}$.\n", " \n", " '''\n", " Gam = self._Estep(X)\n", " return Gam\n", " \n", " def predict(self, X):\n", " '''\n", " Method for make prediction about which cluster input points are assigned to.\n", " \n", " Parameters\n", " ----------\n", " X : 2D numpy array\n", " 2-D numpy array representing input data, where X[n, i] represents the i-th element of n-th point in X.\n", " \n", " Returns\n", " ----------\n", " pred : 1D numpy array\n", " 1D numpy array, with dtype=int, representing which class input points are assigned to.\n", " '''\n", " pred = np.argmax(self.predict_proba(X), axis=1)\n", " return pred\n", " " ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "collapsed": true, "id": "8Txq0agtzpMB" }, "source": [ "def get_meshgrid(x, y, nx, ny, margin=0.1):\n", " x_min, x_max = (1 + margin) * x.min() - margin * x.max(), (1 + margin) * x.max() - margin * x.min()\n", " y_min, y_max = (1 + margin) * y.min() - margin * y.max(), (1 + margin) * y.max() - margin * y.min()\n", " xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),\n", " np.linspace(y_min, y_max, ny))\n", " return xx, yy" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "scrolled": false, "colab": { "base_uri": "https://localhost:8080/", "height": 486 }, "id": "620KSOcQzpMB", "outputId": "bf4a3578-f084-400b-fdd7-e778c137d48a" }, "source": [ "from sklearn.datasets import make_blobs\n", "\n", "X, y = make_blobs(centers = [[0,0],[3,1],[1,3]], n_features=2, n_samples=100,cluster_std=[0.5, 0.5, 0.5], random_state=1)\n", "plt.scatter(X[:,0], X[:,1], c=y)\n", "\n", "xx, yy = get_meshgrid(X[:, 0], X[:, 1], nx=100, ny=100, margin=0.1)" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAHVCAYAAADPSuPPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3hUZfr/8ffUlEno0ntVpEgVQSniCqsisoiKK3ZcXcXuuha+K67uruVnWUVUsGEXFBHFguCCiEAEBAVEkN6kJSQzKZOZOb8/DgRCJpRkZs7M5PPKdV+7PDPznDtRufOc8xQbYCAiIiKWs1udgIiIiJhUlEVEROKEirKIiEicUFEWERGJEyrKIiIicUJFWUREJE44rU5g165dbNq0yeo0REREYqJZs2bUrVs37GuWF+VNmzbRo0cPq9MQERGJiaysrHJf0+1rERGROKGiLCIiEidUlEVEROKEirKIiEicUFEWERGJEyrKIiIicUJFWUREJE6oKIuIiMQJFWUREZE4oaIsIiISJ1SURURE4oSKsoiISJxQURYREYkTKsoiIiJxQkVZJE65XHDrrbB8Ofz0E9x9N6SkWJ2ViEST5ecpi0h4M2bAmWeCx2P++aGHYOhQ6NsXDMPS1EQkSjRSFolDffqYcbAgg/n/O3eGc8+1Li8RiS4VZZE4dMYZ4HaXbc/IMIu1iCQnFWWROLRjBxQWlm3Pz4ft22Ofj4jEhoqySBz66CMIBCAUKt0eCMB771mTk4hEn4qySBwqKIB+/WDtWvD5zFi/HgYOhJwcq7MTkWiJSlFu3bo1BQUFvPnmm9HoXqRK+PlnOPlk6NQJTjsNWrWCJUuszkpEoikqS6LGjx9PVlZWNLoWqXLWr7c6AxGJlYiPlC+99FJycnKYPXt2pLsWERFJahEtypmZmTz88MPceeedkexWRESkSojo7et//vOfvPLKK2zbtu2o7xs9ejQ33HADAHXq1IlkCiIiIgkrYkW5c+fOnHPOOXTp0uWY7504cSITJ04E0LNnERGRAyJWlPv370/z5s3ZvHkzABkZGTgcDtq3b0+3bt0idRkREZGkFbGi/PLLL/PeYbsa3H333TRv3pybbropUpcQERFJahErygUFBRQUFJT82ev1UlhYyJ49eyJ1CRERkaQWtaMbx40bF62uRUREkpK22RQREYkTKsoiIiJxQkVZREQkTqgoi4iIxAkVZRERkTihoiwiIhInVJRFRETihIqyiIhInFBRFhERiRMqyiIiInFCRVlERCROqCiLiIjECRVlERGROKGiLCIiEidUlEVEROKEirKIiEicUFEWqSSbDRo2BI/H6kxEJNGpKItUwkUXwfbtsHYt7NkDb70F6elWZyUiicppdQIiiapXL7MIHz5C/tOfzD8PG2ZdXiKSuDRSFqmge++FtLTSbWlpMGgQNGhgTU4ikthUlEUqqFUrsIf5L8jvh0aNYp+PiCQ+FWWRCpo3zyzAR3K5YM2a2OcjIolPRVmkgh5/HPLzIRA41Ob1mu15edblJSKJS0VZpII2b4bu3WHKFNi5E376CW68EcaNszqz6KtWDc46C1q3tjoTkeSi2dcilfDbb3D55VZnEVv33w8PPghFReat+h9/hKFDYe9eqzMTSXwaKYvIcRs6FO67z5xlXqOGufyre3f44AOrMxNJDirKIhYYNgy+/BLmz4dbboGUFKszOj533QUZGaXbUlLgjDO0DEwkEnT7WiTGnn0Wrr32UHE77TQYNQr69Ck9aSwenXRS+PbiYqhVC3bsiG0+IslGI2WRGGreHEaPLj3a9HjglFNg+HDL0jpuM2eaz5KPFAxqGZhIJKgoi8TQWWeFHw1nZpo7gcW7xx6DffugsND8cygEPp95Cz7eR/kiiUC3r0ViaO9es5Adye83l1XFu127oGNHuO02OPdcc1nY//t/sGiR1ZmJJAcbYFiZQFZWFj169LAyBbFIs2Zw+unmc8j588Gw9N/E2HA6YetW89ns4Vt0+nzQubO5xEpEktvR6p5GyhJzNhtMmABXXmlOEAJzBHb22bBli7W5RVsgAAMGwKefmoU5FDJ/GbnqKhVkEVFRFgtceSX8+c/mWteDpyylp8PUqebIOdmtXm0eZtGpk/l9L1ly6JcTEanaVJQl5m65pexaV6fTfFbZuLF5e7cqWLHC6gxEJN5o9rXE3JEF+aBgsPzXRESqAhVlibmpU6GgoGy716u1riJStakoS8w9+aQ5ocvrNf/s95uzj6+6qmrMwD5ctWpwxx0wY4b5c2nZ0uqMRMRKeqYsMbd/v7m15J//bK513bgRXnqp6s0+rlfPnOR18GCHc881j348/3yYO9fq7ETECirKYomCApg0yYx4VLOm+Yw7Nzd613joIXNZlNtt/tntNuO11zRiFqmqdPta5DCnnmqOXnfsMNdOz50LTZpE51pDhhwqyIerXx8aNozONUUkvmmkLHJA9erw7bfm/x7cbat3b3O3sVatIr+3c15e+Ha7HfLzI3stEUkMGilLldSjB7zxBnz+OdxwA6Smms+43e7S2186nWaR/uMfI5/D88+bE9wO5/fD//4HOTmRv56IxD+NlKXKGT0ann7aLMQOh3ly0803m7eqPZ6y709JMY9cjLQJE6BrVxg50izGdjusWwdXXBH5a4lIYtBIWaqUjAyzIHs8ZkEG8/+3amW+Fu6WciAAS5dGPpdQCK67zjxL+ZprzD2xu3aFPXsify0RSQwaKUuV0rNn+H2mPR5o0QK2bTP/NyXFbM/PNwvyd99FL6dNm8wQEdFIWaqUnJxDI+TDhULmCPWMM+DFF83Z11u2mGcFDxoU+zxFpGrSSFmqlKVL4fffzdOZDi/OBQUwfrxZtG+/3QwRkVjTSFmqnMGDYfNmc2OQnByzID/0kDnrWUTEShopS5Xz22/mjlm9ekGtWrBgQeSXILVvb47Ef/656u3nLSIVp6IsVdbChZHv87TTYNo0qFPHLMZ5eXDJJdGdKCYiyUNFWSRC0tNhzhxz3+yDMjPNDUpatIC9eyN3rdNOgwYN4IcfYPfuyPUrItbSM2WRCBk2LPzMbofD3CAkEurXhx9/NLcDfecdcynVY49Fpm8RsZ6KspRSowZ06mSe8ysnpl69Q+ubD5eWFrkDJqZMMZ9XZ2SY/6zS0uCvf4VLL41M/yJiLRVlAczR3IsvwvbtMG8e7NwJTz4JNpvVmSWOb78NvzGJ1xuZmd2NGkG3buBylW7PyNASLpFkoaIsADz8sLnnclqaeQBDWhrceCPccYfVmSWOrCyYNav0IRM+n3kU5KxZle+/evXyT6o6/Dm2iCQuFWUBYMyYsocxeDxw993W5JOIUlLg/fdhxgz46SdzEtbf/w7nnhuZZVFr1pgHVxypsBA+/rjy/YuI9TT7WnA6w5+OBOY6Xjm2U04xT5lKSTGPfwwE4JtvzEcCFTmHuX176NvXnFn96adQVATBIFx/Pbz1lnkdp9Pcm3vPHnj88ch/TyJiDcPKyMrKsvT6CjNWrsQwjLLx7bfW55YIsWoVRjBY+mfn9WLcfPOJ9WOzYbz2GobPZ35+/36MPXswOnY89J6OHTEmTsT46iuMu+/GqFbN+u9foVAcfxyj7sV1cooYxYABZhEIBMyCEgiYf+7e3frc4j1atDCLaLhfapYtO7G+Ro7EyMsr3UcwiLFhg/Xfp0KhiEwcre7pmbIA5q3Wvn1h+nTz2eUHH8Dpp5vPReXonE6wl/Nf0pEzpY/lL38xZ1Mfzm6H2rWhY8eK5SciiUPPlKXE0qUwfLjVWSSecePM58hHys+HN988sb7C9QNgGOHXQItIctFIWaQSzjsPLrig7EjZMGD1anj22RPr7+23Sy+pOsjvh2XLKp6niCQGFWWRSrj0UnN/6yPl55vbXxYWhv+c3Q5NmpT97MSJZvHNyzP/XFhoFunLLzdnX4tIctPta5FKOLhU6cg9rwOB8gvyxRfD+PHmMjSHwzxVavRos/j6/dCvnzn6Pucc2LEDJk+Gbdui/72ISHyI2IyyN99809i+fbuxf/9+Y82aNcZ1111XqVloCkW8R+/e5iz1I2dd79+PkZp6fO/Pz8eYNs3670WhUMQmYrYkqn379obb7TYAo127dsaOHTuMrl27ViY5hSLu46GHzMJ6cF1xXh7GwIHh3zt9etn1zAcLc/361n8vCoUi+hGzJVGrVq3Cf2AfwIN/27Rq1SqSlxCJOw89ZO7odeedcNNN5jnHs2eHf2/LluGXT/n9kTtJSkQSV8Qneo0fPx6fz8eaNWvYsWMHM2fOjPQlRGKid29zm8yJE83nu0ezaRO8/LJ5xrHXW/77vv02/P7VLpe5PlxEJOJDc7vdbvTp08d44IEHDKfTWeb10aNHG1lZWUZWVpaxYcMGy28lKBI/7HaM887DeOwxjNtuw6hTp3L9PfrooR3OgkHzlvSkSYde79kT49lnMcaPx+jb9/j7bdoUIzsbo7j40K3rvDyMBx+0/meoUChiE5ZtszlhwgRjzJgxlUlOoThmuN0Y8+Zh5OaaRc7vNwvpjh0Y992H4XKdWH9t2oTfNjMvD+P00zHGjStbsJ9//vj7b9kS4623MLZtM7fhvOwy63+GCoUidnG0uhfVJVFOp1PPlCXqbroJunY9dNLVwa0t69eHBx80b0MPGXL8/Q0eDDZb2fa0NLjySrj6akhPP9SekWG2vfaaeXbysaxfb55dLSJypIg9Uz7ppJO49NJL8Xg82O12zj33XEaOHMns8ma8iETIlVeWf/RkejoMGACdOh1/f/n54TfqCATMDT/CTdRKT4eXXoLmzY//OiIiR4pYUTYMg5tuuomtW7eSnZ3Nk08+ye23386MGTMidQmRsI51XrFhmCPp4zVtWvjCGwzC99+HL9g2G5x2Gvz4I7RuffzXEhE5UtzeW1cojieuvz78Bh6Hb+QxYID53qZNMV59FWPzZowffsC45JLwfZ53nvmsOCfH/LzPh/HnP5sTyMo7pvHgkZfvvmv9z0ShUMRv6DxlRVKH3Y4xdapZmEOh0kXS78dYswbDZsNo2BBjzx6z7fDJWw88EL5fjwfjT3/CGDECo3r1Q+1Dh5qF+chrHYwdO6z/mSRLnHoqxv/+Z/6yk5eH8dxz4XdKUygSKVSUFVUiunXDeOIJcxRcWGjGl19i1Ktnvv7002bbkUXU58PIyCjdV1oaRseOGCedFP5aDRpgFBWFL8o//mjdz8Bmwzj3XHNp2F13JfYuYQ0amHcqDt8BLT8f44svrM9NoahMqCgrqlzUro2RmVm6bdmy8EU0J8dcd3zwfXfeeejWdUEBxkcfmaPmI68xebJZJI5cNjVyZOn3dexojq6bN4/u9+x0YsyaZS4NCwYPbf1Z3paf8R6PPGL+/MP9EtWunfX5KRQVDRVlhQLz0Ify9p1u3Nh8z0UXmYX1yNfff79sf2lpGFOmmIUjJ8csFvfff+j16tUx5s83C2N2ttnP229jOBzR+f7Ke7a+e3f0rhnNmDkz/C9R2dnmLzlW56dQVDRUlBUKMM44o2zRKijA+PzzQ+9ZuDB8ISgoKP1c+fCoWxfjtNPKjqbfe6/sSM/rxfjb36Lz/c2bV/6dgNNPt/7nf6Lxz3+WP1Ju29b6/BSKikbMDqQQiWfffw/XXAO7d5v7UxcWwsyZcMklh97ToEH4zwYCULt2+Nd27TKXQvl8h9pSUuCiiyA1tfR7PR645ZbKfR/lKS4O326zlf9aPBs/3vxnFAodasvPh2++gV9/tS4vkWhSUZYqZcoUc6evzp3NAjx8OOTlHXr9m2/Cr3v2+81DJ45XSkr4XcHA3AEsGiZODH8YRm4uLFsWnWtG086d5m5sc+aYv1Ts329u0DJ8uNWZiURX3A7jFYpYR4sW5jPLw5dNeb0Yo0adeF8rV5a99VpcbD5XjkbuNhvGG2+Y+ebnm+urs7PNWelW/1wVCsWh0DNlheIEomlTjBdfxFi92lx+079/xfrp1cucNHZwGZbPh7FrF0aTJtHNv1MnjDFjzFng6enW/zwVCkXpsOxACpFEtHkz3Hhj5ftZuBA6doSbb4b27WHBAvN85r17K9/30axYYYaIJB4VZZEo2rgR7rnH6ixEJFFoopeIiEicUFEWERGJEyrKkrSqV4caNazOQkTk+KkoS9Jp0QLmzzc39di5ExYvhnbtrM7qkBEjzHXDO3bA1KnxlZuIWC9up4YrFCcabjfG9u3mUX8H1wYHg+b+z+EOlYh13HVX6b21AwFzPXHr1tbnplAoYhPaZlOqjKFDza0sHY5DbXa7ud3lpZdalxeYOfzjH6V39HI4ID0dxo61Li8RiR8qypJUmjcvu980mIWwefNYZ1NaixZgGGXbnU7o0yf2+YhI/FFRlqSyZAkUFZVtz8szX7PSzp3gdod/bcOG2OYiIvFJRVmSyjffwMqVUFBwqK2w0NzE49NPLUsLgOxsc2JXfn7pdp8P/vUva3ISkfiioixJxTDg7LPhqadg61bYvh2ef968PRwMWp0dXH89vPee+UtDfr45Q3z0aPOXCRERG+aML8tkZWXRo0cPK1MQibn0dKhZ01wWdfh5wSKS/I5W97T3tYgF8vPL3sYWEdHtaxERkTihoiwiIhInVJRFRETihIqyiIhInFBRFhERiRMqyhJXBgyAL76AVavghRegSROrMxIRiR0tiZK4cdVVMH68eaAEQKtW5iESXbrA5s3W5iYiEgsaKUtccDrh6acPFWQw94nOzDRPVhIRqQpUlCUutGhhFuYjuVzmLW0RkapARVniwt69ZgEOZ+fOyF3HZoNatcL/AiAiYjUVZYkL+/bBzJmlT3cC8HrhP/+JzDWuvtos8Nu2mdd7+GGzSEdC/frQuzecdFJk+hORqsuwMrKysiy9viJ+IiMD4+OPMfLzMXJyMPLyMO64IzJ9DxuG4fViGMah8HoxHn64cv26XBhvvmnmnJ1t/u+kSRgOh/U/T4VCEZ9xjLoX18kpqmDUq4fRuTNGamrk+ly+vHRBPhh5eRhOZ8X7ffxxDJ+vbLEfO9b6n6NCoYjPOFrd0+1riTu//w7Ll0NhYeT6bNo0fLvTac7wrqgbbzSPYTycxwNjxlS8TxGpulSUpUpYsSJ8e14e5ORUrE+brfQSrsNVr16xPkWkalNRlirh738Hn690m88H994LhlGxPg3DHNGHs3BhxfoUkapNRVmqhO+/h3POgblzzZHxihUwahS89lrl+v3rX83iXlxs/rm42Bx933Zb5XMWkapHqzWlyli4EPr3j3yf3bvDPfdA586wZAk88QSsWxfZ64hI1aCiLFJJv/wC111ndRYikgx0+1pERCROqCiLiIjECRVlERGROKGiLCIiEidUlEVEROKEirKIiEicUFEWERGJEyrKIiIicUJFWUREJE6oKItEmctldQYikihUlEWi5LrrYPt281zoHTvg+uutzkhE4p32vhaJgquugmefPXTecv368MwzEAjA669bmpqIxDGNlEWi4OGHDxXkgzwes11EpDwqyiJR0KjRibWLiICKskhUbNx4Yu0iIqCiLBIVf/sb+Hyl23w+s11EpDwqyiJR8NFH8Oc/w+rVUFQEv/wCV1wBH35odWYiEs80+1okSqZPN0NE5HhppCwiIhInVJRFRETihIqyiESV0wn33mvOPP/9d3jlFahXz+qsROKTnimLSFS99x4MHnxoM5VRo2DQIDjlFMjLszY3kXijkbKIRE27dvDHP5be3czlgho14OqrLUtLJG6pKItI1HTtau73fSSPB848M/b5iMS7iBVlt9vNpEmT2LhxI7m5uSxbtozBgwdHqnsRSUAbNoDNVra9sBDWrIl9PiLxLmLPlJ1OJ1u2bKFfv35s3ryZ8847jw8++ICOHTuyadOmSF1GRBLIwoWwfr35/NjtPtReXAwvvWRdXhIf+jVvxh29etEwM5OZa9fx30WL2FdQYHValorYSDk/P59x48axadMmDMPgs88+Y8OGDXTr1i1SlxCRBDRwIHz1lbmzWVERrFxptm3bZnVmYqUbunXls8svZ0i7dvRo1Ii/9enNjzf+hVppaVanZqmoPVOuW7cubdu2ZeXKldG6hIgkgL17YcgQqF0bGjaEDh0gK8vqrMRKaS4n/2/QIDxuN/YDzzfSXC7qpKdze6/TLc7OWlEpyk6nk7fffps33niDNWEeHI0ePZqsrCyysrKoU6dONFIQkTjj88G+fVZnIfGgU716BEOhMu1pLhfnt21rQUbxI+JF2Waz8eabb+L3+7nlllvCvmfixIn06NGDHj16sGfPnkinICIicWy3Lx+XwxH2tZ1eb4yziS8RL8qvvPIK9erVY/jw4QTCrYUQEZEqbX12Nst37sQfDJZq9/r9PPX99xZlFR8iWpQnTJjAKaecwpAhQygsLIxk1yIikkQueu99lu3Ygc9fTE5hIT5/MffPns3s9RusTs1yRiSiadOmhmEYRkFBgZGXl1cSl19++VE/l5WVFZHrKxQKhSLxok3tWkbvJk0Mj9tleS6xiqPVvYitU968eTO2cLsEiIiIlGPt3n2s3asZgAdpm00REZE4oaIsIiISJ1SURURE4oSKsoiISJxQURYREYkTKsoiIiJxImJLokRipWFDuPZaaNEC5s6F9983Tx8SEUl0KsqSUPr0gS++AKcTUlNhxAi47z44/XTIzbU6OxGRytHta0kob70FGRlmQQbIzITmzeHeey1NS0QkIlSUJWE0bw5165ZtT02FSy+NeToiIhGnoiwJo7AQytvJtaAgtrmIiESDirIkjJ074aef4MgTQX0+ePFFa3ISEYkkFWVJKCNGwNat5qQurxfy8+HTT2HCBKszExGpPM2+loSyeTO0agUDB0LjxrBoEaxaZXVWIiKRoaIsCScUglmzrM5CRCTydPtaREQkTqgoi4iIxAkVZRERkTihoiwiUdW2dm1OrlPH6jREEoImeolIVHSuX48PL7mE+hkZGMDe/HxGTJlC1rbtVqcWVV0a1OfePn1oW7s28zdv4YkF37FlvzZml+OjoiwiEedxu/jmqquonpqK/cA2bBluN7NGjaLZM8+wvzA5j/X6Y5vWTBkxglSnE4fdzql16zKqcyd6vDyRdfv2WZ2eJADdvhaRiBt+SnucdntJQT7IYbdzWYcOFmUVfS9ecAEetxuH3fyr1e1wkOl28++BAy3OTBKFirKIRFz9jAxSnWVvxHlcLhpkZFqQUfSd5EmnrsdTpt1ht9O/RfOIX69RtUza1ald7n7wkph0+1pEIm7+5s0UBYO4HI5S7V6/n/mbN1uUVXR5/f5yX9uXH7kTUxpVy2TqiEvoXL8eQcPA6/dz1bSP+eq33yJ2DbGORsoiEnELtmzh202b8B1WqHx+P0t27GD2hvUWZhY9BcUBPli5koLi4lLtXr+fJxZ8F7HrzLnqKro1bECay0WG2039jAw+uvQSWteqFbFriHVUlEUkKi589z3u/fprlu7YwY87dvLAnDkMevMtDMPqzKLnxk8/5avffqOguJicwkIKiot5ISuLSUuXRaT/Pk2bUD8jo8wdCJfDwY3du0XkGmIt3b4WkagIhEKMX5zF+MVZVqcSMwXFAS56730aVcukSbXq/LJnDzmFhRHrv0FGJkaY32rcDgctatSM2HXEOirKIiIRti03j225eRHvd/G2bbiPGCWDeYt81vrkfCxQ1ej2tYhIgti8fz9vLF9ealJZYSDATq+XycuXW5iZRIpGyiIiCeSmTz/j+y1bufX0nmSmpDBl5SqeXLCA/CMmmEliUlEWEUkwk5cv18g4Sen2tYiISJzQSFlEkorDbqNp9erszS8gtyg599gG8/scfkp7hp7cjj2+fCYtW8pPv++yOi2pJBVlEUkaV3XuzFODB+F2OHDa7UxdtYobZsygoDhgdWoR5bTb+frKK+naoD6ZKSkEQiGu69qFMTM/57Uff7Q6PakE3b4WkaRwTsuWjD//PGqlpZHhdpPqdDL8lFN4behQq1MDoGO9ukwZMYK1t47h48suo3vDhhXua2THDnRr2IDMlBTALNIet5vnzvsjGW53pFIWC6goi0hSuP+sM/EcUZDSXC6GnnwyNdNSLcrKdHrjRnx/3XUMO+VkWteqxZC2bfjf1VdxdosWFervsg4dwhbf4lCIM5s2rWy6YiEVZRFJCs1q1Ajb7g8GqefJiHE2pT0zeHCpIx3th41sKyK3qIhQuJ297Hba1amNw66joxKVirKIJIVvN20mEAqFfW19dnaMsymtS/36YdtPrlOnQgX0pR+WhF2XnOpyMa5/f7bdeSddGoS/psQ3FWURSQoPz52Lz+8vVZh9fj8PzpmDPxi0MDPYVxD+6Ma8oiKCoRM/oeN/Gzfy2Pz5FBQX4/X7S/bDtttsVE9NpV5GBl+NGoXTrr/iE43+iYlIUlifnU33lyfy/s8/s2X/fhZt3cqfP/qI5xYtjuh16mdk8PpFF7Hnb/ew9c47+L9+/XA5jv5X6ZMLFpQ6xhLMXxj+u2hRhfN4ZN63tHj2Wb7bvCXsrWyX3c7AlhV7Zi3W0ZIoEUka6/bt44qPpkWt/8wUN0v+cgMnpaeXHJ94b5/e9GjUkCHvvFvu555euJB6GRnc0rMnxcEgboeD13/8kXFz51Yqn9+9PvYVFJQ8qz5StQOzsyVxqCiLiBynKzt3pnpKSqnzjNPdbgY0b06HunX5eVf4zTsMA+6d9TX/nDuPZjWqs2V/bsQ2Nvn4l18Y0q5tmdnYKU4n/9u4MSLXkNjR7WsRkePUu0mTMsuuAEKGQef69Y75ea/fz8pduyO609iHq1eRtW1byclRoVAIn9/PP775ht2+/IhdR2JDI2URkeO0evduCoqLSXO5yry2ITvHgowgGDI49823GHFqey499VT2FxXx0g9LWLBlS8l7RnftyoP9+tIgI4Nf9uzhrq++YtZvOn85HtmAE5/6F0FZWVn06NHDyhREJEn0adqERpnVyNq+LSpFsl6Gh1/HjCHD7cZuM5cyFQUCrNm7l84TXoz49SLhjl69ePjsAaVub+cXF3P+2+/o9rZFjlb3NFIWkYTXMDOTOVddSYPMTMCcefzuzz9z/b1K2D4AACAASURBVCefEGZicoX97vXR//XXeeXCoZxa9yQMw+CLdeu47pNPIneRCHLYbfxf/35lnjenu1w8evbZ9Hn1VYsyk/KoKItIwvtgxMW0qlWr1LrcS049le+3bGXS0qURvdayHTvp+tJLZKa4KQ6GKAzE72EXNVPTSHWG/2v+5JPqxDgbOR6a6CUiCa1ehoduDRuW2Sgjw+3mlp7RezSWV+SP64IMkF1YUO7GKev27YtxNnI8VJRFJKF5XG6C5WyvWdVPTAqGDP7z7fySmdkHHdzpTOKPirKIJLQNOdlkFxaWaS8MBPhw9WoLMoq9FKcD92Frpw/37/nzGTtnDrt8PgzD4Ld9+7j8w49iPvvapjMyjoueKYtIQjMMuHLaNGaMHInTbifF6cTn9/O7z8d/5s+3Or2oalGzBq8OHcqZTZtiGAZfr9/A9Z98wva8vFLve2bhIp5ZuAibjYhOfDseg1q34tnBg2lbuzbZhYU8Pv87Hl/wXczzSBRaEiUiSaFFzRrc2L07LWrUZM6GDUxevjzsSUoV1bJmTf7Rvx99mzVjW24u/54/n89+XRux/k9UusvF+ttupU56esk2m8XBINvy8mjz3+fKPTErls5s2pQvrvhzqQ1XvAf2/H5gdtW9fa4lUSKS9DZk53DvrK+j0neLmjVY+pcb8LjdOO12mteowfsXX8w9X81iwg8/ROWax3Jph1PxuFyl9r12ORzUSk3l/LZtmP7LGkvyOty4Af3L7ICW4XZz2+mn8/DcuRQFKnd6V7Ma1enesCHbcvNYuHVrpfqKFyrKIpJQ6mV4GH5Ke9JdLj5b+yurd++J+jX/r1+/koJ8kMft5j9/OIdXli2z5GjI1rVqkRHmwIkUp5PWtWrFPJ9wTq4TftmVgXna1qac/RXq12aDly64gCs6dcIfDGK32diSm8s5kyezI89biYytp4leIpIwhp1yMutvu40n/nAOj5w9gB9uuIHH/nBO1K/bt1mzsGcT2zBH0VZYvvN38sLsoe0PBlnx++8WZFTWyl27wh4rCVSqeF7bpQsjO3YkzeWiemoqmSkptKlVi/cuvrjCfcYLFWURSQjVUlJ4c9gw0l0u0t1uUpxO0l0u/tqjB72bNInqtbfsDz+iczkc7PL5onrt8kz7ZTU7vV6KDlsrXRgIsG7fPr5eHx/7Wo/95hsKikuv5fb6/Tz+3XeVurswpmfPMsvdXA4HPRs1oq7HU+F+44GKsogkhMGtWxMMM+pKczq5olOnqF773/Pn4ztirW9BcTGf/LKG7IKyy7FioTgYotekV5i8fDnZBYXszc/n5R+W0O/11+NmZvOirdu44J13WLpjB8XBINtz87jva/MIy8rILOec6GAolPBr0/VMWUQSQnnrXG2APcprYL9c9xu3ffEFT557Lg6bDZfDwfQ1a7h2+vToXvgY9hUUcMOMT7lhxqeW5nE0/9u4kW4vvRzRPqetXs0tPXuScsQWotmFhWzIyY7otWJNRVlEEsIX69aFfa6bHwjw9k8/Rf36ryxdxuTly2lWvQZ78vPJCbNhiRxb94YNuazDqdhsNt7/eSWLt2074T7+9e18hrdvz0np6XjcbooCAQKhEFd//HHc3CWoKBVlEUkI+wuLuHb6dF4dOhS7zYbTbqcoEODVZcv4dtPmsJ9xOew8PGAAf+nWjQy3mwVbtjDm88/56fddFcqhOBjSntGV8MjZA7i9V6+SQzL+0q0bzy/O4u9fn9hStn0FBXR8YQJXndaZgS1a8Nu+bF5c8gO/7UvsUTJo8xARSTCNqmVyyamnku5y8emvv7J8Z/kzjd+/+GIuaNuG9APPGUOGgc/vp8MLE9hczuStROSw26iRmkp2QWG5s52t1q5ObZb+5S+ku1yl2vOLi+k5cSIrd+22KLPY0+YhSSSVVAYzmHTS+Zqv2UXFfuMXSVTbcvN4+vuFx3xfk+rVGNKuLWmHFQG7zUaKw8EdvXpxx5dfRjPNmLnvzDO598wzSXE6KCguZtzcuTy7cJHVaZUxpG07HGEmBrjsdi5s165KFeWjUVFOIGdyJp9iTuiwYcOFiwd5kKd4yuLMROJPu9p1KAoGSxVlALfTSbeGDSzKKrJu73U69/c9q2TGcarTySNnn43X7+eVpctKvbdWWhqXdehAXY+HbzdvYvb6DTHNtSgYCDt7PhgKVXpnr2SiopwgUkhhBjOoTvVS7Q/zMHOZyxKWWJSZSHz6de9eUsKcnOQPBlm2c2dMc3E57HRt0ICC4kBEN/a4/6yzyiwBynC7Gdu3b6mifEaTxnx5xRXYbTbSXC58fj+Ltm3jvLffpjgYmz2yP1y1msfOKbvRSwiYsmplTHJIBCrKCeJczsVG2Vs/KaRwDdeoKIscYfP+/cxcu5Y/tmlT6jlmUSBwXLe/I+WCtm2ZPOwi7DYbDrsdr9/PLq+XzJQUZq1fzz/nzmNrbu4J92uzQZ309LCvNcjMLPW+qZdcUmptb2ZKCmc0bsx1XbryYoz27t6el8d1n3zCpAsvLDn/2mm3c9Onn7Jl/4l//8kqopuH3HzzzWRlZVFYWMhrr70Wya6rvHTSwxZlJ04yyLAgI5H4d/mHHzEhK4u8oiKCoRALtmyh3+uvszEnJybXb1mzJu9dPJyaaWlUT00lw+2mnsdDp/r1aVGzJtecdhrLbvwLDTJP/L9hwzAP4Qjnlz2H9gPvWLcemWE21PC43VzT5bQTvm5lvPvTzzR56mlunjmTW2Z+TpOnn2by8hUxzSHeRXSkvH37dh555BEGDRpEWlpaJLuu8mYzGxeuMu155DGVqRZkJBL//MEgd381i7u/mmXJ9a/r2gXXEWurbYdNdnI5HGS63dzduzd3ffnVCfd/51df8s6f/lQyuxzM2cx3f3Wor5BhlLrm4ayYqb2voIA3VYjLFdGR8rRp05g+fTp79+6NZLcC7GEP93APPnwEMPeSzSOP2czmMz6zODsRCadBRiZu59HHPilOJ2e3aFGh/qf/soZh739A1rZtZBcU8P2WLZz/9jvM+u3Q3tc/79rFvoKCMp/1+v1MWrq0QteV6NEz5QQynvHMZz5XczXVqMZHfMRMZmJYu9RcRMrxxbp1XNz+lHL3agZztLqxnNvQx+Or337jq99+O+p7/vT++3x95ZU4bDZSnE78wSCz16/n9R9/rPB1JTosKcqjR4/mhhtuAKBOOedtSnjLWc4d3GF1GiJyHD5avZq7zjiDU+uehOfALWbjiNvJBcXFPLFgQVTzWLJ9B02eepo/nXIK9TI8zNu0iUVbT3x7S4k+S4ryxIkTmThxImDubCIikowCoRB9X3+N0V27cnnHjhQGgmS4XXSsV49AMIg/GOKWz2eyYMuWqOfi9fuZvHx51K8jlaPb1yIiUVQUCPL84iyeX3xoAFI7PY1aaWmsz84mGNLjJzkkokXZ4XDgdDpxOBw4HA5SUlIIBAIEK3GYtYhIstmbX8De/AJsNnOtbiAUmw08JP5FdPb1gw8+SGFhIffddx+jRo2isLCQBx98MJKXEBFJeOkuFxOHDCH/gQcofPABFo2+ni4N6ludlsQBnRIlIhJjX185it5NmpTalzuvqIhTX3hBu1tVAUerexEdKYuIyNGdclIdejVuXOagDJfDwS09e1qUlcQLFWURkRhqU6s2xWGeIac6nXSsW8+CjCSeqCgfYMfOWMayj30ECbKCFfSjn9VpiUiSWbl7V9jTqwqKi8narrXDVZ2K8gFP8zT3ci81qYkdOx3pyGd8Rle6Wp2aiCSR3/Zl88W6deQXF5e0BUMhCgMBxi/Wvg1VnYoykEkmoxmNB0+p9jTSGMtYi7ISkWR16dSpPPX99+z2+cgvLubzdevoOXESu3w+q1MTi2nzEKApTSmmmDRKn2xlx04HOliUlYgkq+JgiLFzvmHsnG+sTkXijEbKwCY2hT0WMUiQ5WhbOhERiQ0VZcCLlxd4AR+lbx0VUsg/+adFWYmISFWj29cH3MM9/M7v3M3d1KQmy1nObdymkXIFNKMZwxmOCxfTmc4v/GJ1SiIiCcOwMrKysiy9viKycT3XGz58RgEFRhFFhg+f8Q/+YXleisSO7g0bGk8NOtd4ZvAg44wmjS3PR6GoTByt7mmkLBHTgAb8l/+WmjDnxs093MM0prGCFRZmJ4lq3ID+3HXGGaQ6zb+uruvalUlLlnLHl19anJlI5OmZskTMEIYQouxORSmkMIIRFmQkia5N7Vrc3bs3Hrcbh92Ow24nw+1mdLeudK6v3a8k+agoS8QYGBV6TaQ857dpi91mK9Oe4nRyYbt2FmQkEl0qyhIxM5iBPcy/Un78TGGKBRlJoisKBgiF2Sf64A5YiS7N5WRsv778OuYWfrnlZu49sw/uMFtwStWhoiwRs5Od3MzNFBz4KqKIfPL5N//mJ36yOj1JQB+uWg1hRsohw+CDlSstyOjEuB0OLuvQgUfOPptRnTuVPBcHsNts/O+qq/n7mWfSpnZt2tWpw//17cuXV1xhYcZiNU30koh6jdeYxaxSS6LWstbqtCRB7fL5uHLaNCYPG0bwwIjZabdz06efsSlnv8XZHV1dj4dFo6+ndloamSkp5BUV8Z9zzuH0iZPYmpvLoNatOOWkOqQfdoRjuttNt4YN6NusGfM2bbIwe7GKirJE3Fa28izPWp2GJIkPV61m9voNnN+2DXabjZlr17I3v8DqtI7pmcGDaJSZievA7ejMlBTSXC4mXHA+Q955l9MbNSbD7S7zuVSnk16NG6soV1EqyiIS93IKC3l7RWI9Ahl68sklBfkgp93O4NatsdlgW14uvuLiMoW5MBBga25uLFOVOKJnyoepRz1GH/iqT32r0xGRBBYywq84MA60v//zSoqDoVLvCx2YwPbR6tUxyVHij4ryAddwDRvYwNMHvtaznmu51uq0RCRBTVm5kqIjZoj7g0E+WbMGw4DcoiL6v/46a/bsoaC4mILiYn7evZu+r71eMrO8SfVqjOrcifPbtsHl0F/XVYENrF1AmpWVRY8ePaxMgSY0YQ1ryhzdmE8+J3MyW9hiUWYikqhqpKby3bXX0rh6NVIdDgqDQXb5fPR55dUy5yY3rlYNA4NtuXklbY//4Rxu6dmTQMgcTfuDQQZOnsxPv++K9bciEXa0uqdnysBwhmOj7LILGzYu5mKe5mkLshKRRJZTWEjHCRMY1LoVp55Ul1/27OHzdWsJhsqOg458hnx+2zbc1KMHaYfNzA4ZBp9d/meaPfM05dwZlySQdEW5OtX5A3+gmGK+4isKOPYsTReusJte2LGHPWdZROR4hAyDz9eu4/O1607oc3/p1q3MBDC7zUaN1BS6N2xI1rbtkUzzhKW5nLSrXYedXi87vV5Lc0k2SVWUr+IqXuAFiikGzKI6nOHMYtZRPzed6TzEQ7gp/R9BkCDTmR61fEVEwgm3VArMIn/4umYr3HXGGYwb0J9AKITb4eDr9Ru4/MMP8fr9luaVLJJm5kArWvECL5BOOtUPfGWSyUd8RDWqHfWzv/Irj/EYPnwEDnz58PE4j7OGNTH6DkRETO/89DO+MEXOZrOxcOtWCzIyXXTyyTw0oD8et5vqqamkuVz8oWUL3hh2kWU5JZukKcpXcAXOMAN/A4OhDD3m5x/mYc7gDP5z4Ks3vRnHuGikKiJyVJOXL2fZzp3kFRUB5qzt/OJirp0+naJA0LK87u3Tp8woPtXl4rw2baiZlmpRVsklaW5fe/CELcp27HjwHFcfPx34EhGxkj8YpP/rrzPs5FM4r00bdvl8vLJsKWv37ot5Lv2bN+eaLqeR4nDQslbNsO8JhELUTksnu6Awxtkln6Qpyp/yKTdxExlklGq3Y+cLvrAoKxGRigmGDKauWsXUVassy+FfA89mzOmnk+5yYbfZKA4GCYVC2O2lb7L6g0E25uRYlGVySZrb1/OYx8d8TB7mOr8gQbx4eYIn2MhGa5MTEUkwLWvW5PZevchwu0vOtHY5HNgOFGcwJ575/H5u//wLAmGO2JQTlzQjZYBRjGIQg7iMy/Dj5w3eYAELrE5LRCTh/KFVy7BbhRrAqt27cdrtbN6/n8e++465G3V4RqQkVVEG+PLAl0ROV7rSkIb8wA/sZKfV6YhIDOQV+QmGKcrFwSDTfvmFcf+ba0FWyS9pbl9L5NWjHj/yI3OZy5u8yQY28BRPWZ2WiMTAJ2vWhNnnEIKGwZvLV8Q8n6pCRTnBNaUpYxnLeMYzlKFhdyarqClMoT3tySCDGtQglVRGM5rLuTxi1xCR+OT1+xny7rvkFBay/0D4/H6u/+QT1mdnW51eUjOsjKysLEuvn8jxR/5oePEahRQaBoaRS64xj3mGG3el+25AAyOffMPAKBOLWWz5965QKGITbofDGNS6lTGkXVsjw135v1sUR697GiknKCdO3uZtPHhIIQWATDLpSleu4ZpK91+NagQIhH2tOtUr3b+IJAZ/MMiX635jxppftZVmDKgoJ6judA97q9qDhyu4otL9r2Vt2MM8CinUfuAiIlGiopygiigq9/lxIZXfVSdEiOu4Dh++kgM+8slnF7t4jMcq3b+IiJSVdEuiqoplLCObbDx4ShVnL15e4qWIXONTPqUnPRnDGJrTnK/4ild4hVxyj/3hKGpJSxw4WMtaS/MQEYk0FeUENoQhzGY2btw4cGDDxtu8zVSmRuwaq1jFTdwUsf4qowMdmMpUGtMYgF3s4lIuJYssizMTEYkMFeUEtoIVNKIR53EedajDPObxK79anVZUpJHG//gfNalZcmegBS2YxSya05wctO+uSDh10tNJcTrYlptndSpyHFSUE5wfPx/zsdVpRN0whuHCVeY5uhMnIxnJBCZYlJlIfGpcrRrvXTycbg0bYhgG2/LyuOKjj1i0dZvVqclRaKKXJIQGNChZ+nU4Dx4a0ciCjETil91mY941V3N6o0akOp2kuVy0rlWLWaNGUT8j49gdiGVUlBNAIxoxhjHcwR20opXV6VhiAQtKZoEfLo885jPfgoxE4tfZLVpQOz0dp8NRqt1pt3Ntly4WZSXHQ7evY6ge9ehFL3axi+/5/qjvTSGF1rRmAAN4nMcBsGHjUR7lYR7mP/wnFinHje/5nm/5lr70xYMHMJdorWCFDiAROULT6tVLjls8XJrLRataNS3IKPE0yMzg6UGDuKBtW4KGwTs//cTfZs0iryi6G6ioKMfIozzKHdyBHz927OxiFwMZyCbKHnl2K7fyCI8AkEEGtiO2hR/LWGYwg5WsjEnu8eJCLuRGbuQ6rsOBg8lM5jmew8CwOjWRuLJ427awRTmvqIh5m3TM4rGkuZwsHj2aeh4PrgN3G64+7TR6NGxI95cnRvXaun0dAzdwA7dzO2mkUZ3qZJJJM5oxgxll3nshF/Iv/kXmga8jCzKACxcjGBGL1ONKgADP8zxd6EInOvEkT1JEkdVpicSdn3ft4st16/Adti1mYSDADq+X937+2cLMEsNlHTpQPSWlpCADpDqdtK1dm77NmkX12irKUdSEJixhCS/wAumkl3rNiZOWtKQtbUu138d9Jbdny2PDFtHToE5EIxrRlKaWXFtEjt+IKVN4cM43rNmzh005OTy/aDGnT5xEUSBodWpxr0v9+mSmlJ1Y6rDb6VivblSvrdvXUfQ1X5fsPhVOgADVqFaqrQENjtmvH39ENwg5Hm1owwd8QDvaYWCwla2MZCRLWRrTPETk+ARDBs8sXMgzCxdanUrC+XnXbrx+Pxlud6n2QCjEmj17o3ptjZSjpBe9aEADnEf5vcfAYDnLS7XNZW7Y05kMDIopJp98nuAJVhC7Q8ZTSOFbvqUjHUkjjXTSaUtb5jCHmmjSiIgkl3d++on84mKCoVBJW1EgwLbcXGZvWB/Va6soR0l96hMiFPa1ECF8+BjN6DLLfB7iIbx4S7Xnk88UpvAP/kEPevAQD0Uz9TKGMIQ00sqM+J04uZzLY5qLiEi0ef1+ek2axJwNGwmEQviDQT5Zs4azXnsNI8rzSnX7OkoWsSjsZhcBAixgATdzMz9TdsLFBjZwGqfxAA/Qn/5sYhP/5t/MYU4s0g6rEY1w4y7T7sGj58sikpQ2ZOdw7ptvYrfZMDCiXowPUlGOkh3s4Hme50ZuJANzB50CCtjCFv7IH8knv9zPbmITN3BDrFI9psUsDntLPY88FrDAgoxEkl+/5s34a/ce1ExLY8qqlUxevlyTtCwQilU1PkBFOYru4R4Ws5gxjKEGNZjKVJ7hmaMW5Hj0Pd/zHd9xJmeWzAwvoIC1rOVTPrU4O5Hk87c+vfm/fv1Ic7mw22yc0aQxN3TtRp9XX8UfVGFOZirKUTblwFeiG8IQxjCG67keJ07e4i2e5EmC6C8IkUiqnZ7GQ/37k+ZylbRluN2cfFIdRnbswBs/Lj/KpyXRaaKXHJdiinmKp2hPe9rSlod5OOFG/CKJ4MymTcOOhjPcboadfLIFGUksqSiLiMSR7IJCbGG2yAyGQuz26RfhZJd0RdmBg1RSrU5DRKRC5m/eTG5REaFQ6SWVhYEAE374waKsJFaSpiink85EJuLFSx55rGAFZ3CG1WmJiJyQkGHwh8lvsiU3l9yiInIKC/H5i7n9iy9YumOH1elJlCXNRK+pTKU//UtGyR3pyFd8RRe6sI51FmcnInL8ftmzhxbPPkvPRo2olpLC91u24vVH98hAiQ9JMVJuSUv60Y800kq1p5DC7dxuUVYiIhVnGLBo6zZm/bZeBbkKSYqi3IpW+Cn7L60LFx3oYEFGIiIiJy4pivIqVoXd0rKIIhaxyIKMRERETlxSFOVtbGMKU/DhK2kLEqSAAp7lWQszExEROX5JUZQBruVaHuVRtrMdL15mMpOe9GQ7261OTURE5LgkzezrIEH+feBLREQkESXNSFlERCTRRbQo16xZk48++giv18vGjRsZOXJkJLuPqQY04AmeYBGLeIu36Exnq1MSETluHevV5crOnTmzqc48TyQRvX09fvx4/H4/9erV47TTTuOzzz5j+fLlrFq1KpKXibrmNGcJS/DgIYUUutGNi7iIEYzgcz63Oj0RkXK5HQ6mXXYp/Zo1KzkLeGNODgPeeIO9+QUWZyfHErGRcnp6OsOHD2fs2LH4fD6+++47PvnkE0aNGhWpS8TMIzxCNaqVLLNy4MCDh5d52eLMKiaNNIYylBGMoAY1rE5HRKLo/rPOpH/z5njcbjJTUshMSaFd7dq8OnSo1anJcYhYUW7bti2BQIC1a9eWtC1fvpxTTz21zHtHjx5NVlYWWVlZ1KlTJ1IpRMxABuIMcxOhNrWpT30LMqq4gQxkJzt5gzeYxCS2s51ruMbqtEQkSkZ360b6YWcxA7idTga3bk2qM2nm9iatiBXljIwMcnNzS7Xt37+fzMzMMu+dOHEiPXr0oEePHuzZsydSKURMNtlh223YyCMvxtlUXCaZfMzHVKMa1alONaqRRhrP8zxtaWt1eiISBSkOR7mvuRya2xvvIvZPyOv1Uq1atVJt1apVIy8vcYrYQU/yJF68pdoKKeQTPim1QUm8G8IQQoTKtDtxMorEe6wgIsc249dfKQ4GS7WFDIOfd+0ir0h7aMe7iBXlX3/9FafTSevWrUvaOnfuzMqVKyN1iZh5lVeZwAQKKCCHHAooYC5zuY7rrE7thHjw4KDsb81OnGSQYUFGIhJt9309m10+H74Dh1gUFBfj9fu5dvp0izOT42VEKt59913jnXfeMdLT043evXsbOTk5Rvv27Y/6maysrIhdP9JRi1pGX/oazWlueS4ViSY0MfLJNwyMUpFHntGf/pbnp1AoohMZbrfx1x7djTeHDTPuO+tMo67HY3lOikNxjLoXuQvVrFnTmDZtmuH1eo1NmzYZI0eOrGxyikrGWMYaXrxGgEBJQX6HdyzPS6FQKKpqxKwoRyE5xXFEYxobL/KisZa1xjzmGedxXqnXe9HLeJEXjdd4zRjMYMvzVSgUiqocR6t7mh+f4BrTmB/5kWpUw4WL1rSmC124n/t5jucAWHjgS0RE4pvmxye4v/N3MsnExaF1iRlk8CiPkkoqYG4e8l/+y372U0ghn/EZLWlpVcoiIlKOpCrKDhy4cVudRkydzdlhv+cQoZK1yJ/xGddzfckuZYMYxGIWU4tasU5XRESOIimKck1q8h7vkU8+Pnx8x3ecStmdxJLRbnaHbXfj5nd+pxOd6ElP0kgrec2BgzTSEm6Jl4hIskuKojyb2VzERbhx48RJL3oxn/nUoeJbeGaSSV/6cjInV+izf+NvfM/3fMInnMM5Fc7jaKpRjQ50wMAo1R4kyGxm8zu/cyqnEiRY5rPppNOd7lHJS0REKibhi/IZnEFrWpccHgFgx44bN9dybYX6vIu72MlOpjOdH/iBLLKoR73j+qwHDz/wA//gH/SiF0MYwjSmcQ/3VCiXo7mKq0ghBRu2Uu0GBv/iXwD8wi/Yw/xjLqCAZSyLeE4iIlJxCV+U29AmbHs66XSk4wn3N4hBPMRDpJNODWrgwUNnOvMxHx/X56/nehrTmHTSS9oyyGAc46hO9RPO52hO53Q8eMq0F1BAa8yd1ZaxjOUsp5DCktdDhCikkElMimg+IiJSOQlflFewIuxI0IuXRSw64f7u4I4yW1C6cNGJTsc1Y3kIQ0oV5IOKKKInPU84n6P5mZ8pIPz5qGs5dFrXYAbzFm9RQAFBgsxlLr3pzR7i7zAQEZGqLOGL8o/8yEIWlipOAQJ48TKZySfcX13qhm0vppja1D7m53eyM+wzXCfOiBfBV3iFIopKHTpRRBHrWFdqXbIXL6MZTTrpuHBxNmfzC79ENBcREam8hC/KABdwAc/xHHvYQx55TGUq3elOLrnH/vARZjCj1K3eg2zYWMGKY37+OZ4r8/kAAbawJeLPcHezm7M4i8UsJkAAP34+5mMGMrDczxw5KUxEROJLA3wWzgAADYpJREFU3G43ZkXUpKaxiU0lBzkECRpevMa1XHvcfVzLtUYeeUYOOUYeecYKVhhNaRrVvFNIMZw4Lf/5KRQKheLooW02T0A22XSmM3/lr1zABWxlK8/wDAtYcNx9vMqrvMu7dKUr2WSzilVRzNhURFHUryEiItFlw6zOlsnKyqJHjx5WpiAJIp10GtGIrWwtd4KbiEi8O1rdS4pnysmqHe1oT3ur07CcDRuP8zi72c0SlrCHPTzKo1anJSIScSrKcagTnVjHOn7gBxaxiM1sphe9rE7LMn/n79zETaSTTiaZpJPOrdzKndxpdWoiIhGl29dxJp10trCFGtQotf46l1xa0IJ97LMst3rU40IuxIaNGcxgBztict297A17eMZOdtKABjHJQUQkUnT7OoEMYxguXGU2RHHgYCQjLcoKruZqNrCBpw58/cZvjGZ0TK5dgxph2yuzt7mISDxSUY4z9akf9ijGdNItGxU2ohHjGU8aaWSQgQcPaaTxLM/SnOZRv/5qVodtP5514yIiiURFOc58y7cUU1ym3YuXecyzICP4E38qc+gFmAd/XMzFUb/+rdyKD1/JzmUhQvjwcQd3MJCB/JN/cjM363xoEUl4KspxZjGLmcMcfPhK2nz4WMpSZjHLkpxcuMotyk6iv9R9DnM4m7P5ki/ZzGZmMpOBDGQsY5nGNB7gAR7ncTaykT70iXo+IiLRFLc7m1TVcOAwbuAGYzGLjSUsMW7lVsON27J8WtPa8OEzDIxS4cNnnMIpluR0EzcZXrxlctrOdsOGzfJ/hgqFQlFeHKPuxXVylkUNahjDGW4MYYiRSqrl+Vgd93O/4cNn+PEbxRQbPnzGOMZZls9CFpYpyAaGkUuu0ZnOlv+8FAqForzQNpsn6Cqu4gVeKHm2a8PGRVzEN3wTkf49eKhJTbazvdQJT/HsX/yLT/iES7gEGzamMpXlLLcsn/J+bjZsYU/pEhFJFHH7G4MV0YY2YW/V5pJrZJBRqb670MVYxzojSNDw4zf2ste4nMst/54TMQ4e+nHkP6dNbLI8N4VCoThaHK3uaaLXEa7kyrCTlwwMhjCkwv32pjeLWUxLWmLHjgsXtajFJCYxgAH/v737D627uv84/kxvSMxPKMR0w65xa21dhAlqGKuO2c2pRIrryloi6/6pZUg6BmX/jE3utyqO1TH/0ha0izZYhUKyKpXVKZVKJuxCXTuaLbbDptJ22krzs0mTNOf7h7OSGqtJbnpObp+P8obwCTn3lUPoi5z7uTcziXxVeo7neI3XGGCA85ynn3566GEVq2JHk6Rp8/j6ElVUTVrKGTJUUjntdbexbdJ1yyhjN7tZwxr+wl+mvf7VZpxxVrOaBhr4Ht/jAz6gjbYJd61L0lyU7K/xMWYFKyY9Fj3HuWn/TeR5zAsXuDDpjUmfzAADoZHG6N+/4ziOM7vj8fUU7GMfe9hDP/0AXOACgwzye37PcY5Pa81P3uziciqoYAtbprW+JKkweHw9iSaaaKSRtaxlmGGe4zn+xt9mtOY2tvFLfjnpW2h+4gZumNFjSJLmNkt5EoHAnv/9y5ff8Bu+wldYwxpKKJn0HbJOcCJvjydJmns8vp6iSir5Nb/mAAfYxz5Ws/pLfd0oo/yMn1FHHX/kj585zh5kkN/y29mILEmaQ5J9wju1KaMsdNIZznHu4g1a/fSHLWyZ8lq/4BfhNKfDKKPhv/w3rGd99O8vX/MjfhRe5/Xwd/4efsWvQjnl0TM5juOkMr7NZp5mAxs+987sr/LVaa1ZRln07yufs4UtE/ZokMFwiEOhlNLo2RzHcVIY777Ok0YaJ32t8ggjfIfvTGvNIYZmGisZ13EdG9k4YY/KKefrfJ0mmiImk6S5wVKegvd5f9K/dVxEER/wQYREaVnO8kn3p5JKGmmMkEiS5hZLeQq2sY0RRiZcG2OMM5yZ8UumCsFpTk96fZRRTnLyCqeRpLnHUp6CTjpZxzrOcpY++hhkkE46+QE/IBBix4tuP/vpoeczf6VphBG2sS1SKkmaO3yd8hS1084rvMK3+BZ99HGUo7EjJWOccb7P93mFV/gaX+MCFwgE1rOef/Pv2PEkKXmW8jSMMcYBDsSOkaT/8B/qqeebfJNKKvkH/5j0eWZJ0mdZypoV/+JfsSNI0pzjc8qSJCXCUpYkKRGWsiRJibCUJUlKhKUsSVIiLGVJkhJhKUuSlAhLWZKkRFjKkiQlwlKWJCkRlrIkSYmwlCVJSoSlLElSIixlSZISYSlLkpQIS1mSpERYypIkJcJSliQpEZayJEmJsJQlSUqEpSxJUiIsZUmSEmEpS5KUCEtZkqREWMqSJCUiL6Xc3NxMLpdjeHiYlpaWfCwpSdJVpzgfi5w8eZLHHnuMe+65h7KysnwsKUnSVScvpdze3g7AbbfdxsKFC/OxpK5C1VTzU37KjdxIjhy72MUww7FjSdIVk5dSlmZqGcvooINruIYKKuinn81s5tt8m9Ocjh1Pkq6IKDd6bdiwgVwuRy6Xo6amJkYEJeZP/In5zKeCCgCqqOI6ruN3/C5yMkm6cr6wlPft20cIYdJ56623pvWgzzzzDA0NDTQ0NHDmzJlpraHCUUYZDTQw75IfxxJK+DE/jpRKkq68Lzy+XrFixZXIoavYOOMEwqSfG2X0CqeRpHjycnydyWQoLS0lk8lM+Fj6Ms5zntd5/TMFPMQQO9gRKZUkxRFmOtlsNlwqm81+qa/N5XIzfnxn7s8CFoQjHAm99IZBBkMffaGDjlBGWfRsjuM4+ZzL9V7R/z6IJpfL0dDQEDOCEjGPefyQH7KYxRzkIB10xI4kSXl3ud7zJVFKxjjj7GVv7BiSFI2lrC+UIUMjjSxlKYc5zGu8xjjjsWNJUsGxlHVZtdTSQQe11HIN1zDMMMc5znf5Lj30xI4nSQXFvxKly9rKVhaxiGqqKaGEaqq5gRv4A3+IHU2SCo6lrM9VRBErWUkJJROul1LKGtZESiVJhctS1mUVUTTp9UvffUuSNHP+z6rPFQj8lb8yxtiE66OM8mf+HCmVJBUuS1mX9XN+zmlO008/AP30c4ITbGJT5GSSVHi8+1qX9T7v8w2+wU/4CTdyI//kn7TRxggjsaNJUsGxlPWFhhmmldbYMSSp4Hl8LUlSIixlSZISYSlLkpQIS1mSpERYypIkJcJSliQpEZayJEmJsJQlSUqEpSxJUiIsZUmSEmEpS5KUCEtZkqREWMqSJCXCUpYkKRGWsiRJiSgCQswAH374Id3d3TEjJKmmpoYzZ87EjlGQ3NvZ497OHvd29lzpva2rq6O2tvZzPx+c9CaXy0XPUKjj3rq3c3Hc26tjbz2+liQpEZayJEmJyAD/FzuEJnfgwIHYEQqWezt73NvZ497OnlT2NvqNXpIk6WMeX0uSlAhLWZKkRFjKiZk/fz5tbW0MDAxw7NgxmpqaYkcqCM3NzeRyOYaHh2lpaYkdp6CUlJTw7LPPcuzYMfr6+njnnXe49957Y8cqGK2trZw8eZLe3l66urpYv3597EgFZ8mSJQwNDdHa2ho7CpDA67KcT2fnzp3hpZdeChUVFeH2228PPT09ob6+PnquuT6rVq0K999/f3j66adDS0tL9DyFNOXl5SGbzYa6urpQVFQU7rvvvtDX1xfq6uqiZyuEqa+vDyUlJQEIy5YtC6dOnQq33HJL9FyFNHv37g379+8Pra2t0bP4m3JCysvLWb16NQ8//DCDg4N0dHTw8ssvs27dutjR5rz29nZ2797NRx99FDtKwTl37hybN2+mu7ubEAJ79uzhvffe49Zbb40drSB0dnYyMjICQAiBEAKLFy+OnKpwrF27lp6eHt54443YUQCPr5OydOlSxsbGOHLkyMVrBw8e5KabboqYSpqa2tpali5dyuHDh2NHKRhPPfUUg4ODdHV1cerUKV599dXYkQpCVVUVjzzyCJs2bYod5SJLOSGVlZX09fVNuNbb20tVVVWkRNLUFBcX88ILL/D888/T1dUVO07BaG5upqqqijvuuIO2tjbOnz8fO1JBePTRR9m+fTsnTpyIHeUiSzkhAwMDVFdXT7hWXV1Nf39/pETSl1dUVERraysjIyNs3LgxdpyCMz4+TkdHBwsXLuShhx6KHWfOu/nmm7nrrrt48sknY0eZoDh2AH3q3Xffpbi4mCVLlnD06FHg4x8cjwE1F2zfvp0FCxbQ2NjI2NhY7DgFq7i42OeU8+DOO+/k+uuv5/jx48DHJ5WZTIb6+vro90NEv9vM+XRefPHFsHPnzlBeXh6WL1/u3dd5mkwmE0pLS8Pjjz8eduzYEUpLS0Mmk4meq1Bm69at4e233w4VFRXRsxTSXHvttWHt2rWhoqIizJs3L9x9991hYGAgrFy5Mnq2uT5lZWVhwYIFF+eJJ54Iu3btCjU1NbGzxd8c59OZP39+aG9vDwMDA6G7uzs0NTVFz1QIk81mw6Wy2Wz0XIUwixYtCiGEMDQ0FPr7+y/OAw88ED3bXJ+amprw5ptvhrNnz4be3t5w6NCh8OCDD0bPVYiTzWaTeEmU730tSVIivNFLkqREWMqSJCXCUpYkKRGWsiRJibCUJUlKhKUsSVIiLGVJkhJhKUuSlAhLWZKkRPw/VvV+ZlgJX9QAAAAASUVORK5CYII=\n", "text/plain": [ "