{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "name": "GPClassification.ipynb", "provenance": [], "collapsed_sections": [] }, "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "SVOCZvFgNhn8" }, "source": [ "# Gaussian Process Classification\n", "\n", "Code implementing GP classification via the Laplace approximation, based on chapter 6 of [Bishop's PRML](https://www.microsoft.com/en-us/research/people/cmbishop/prml-book/). " ] }, { "cell_type": "code", "metadata": { "id": "NxM449iDcXjJ" }, "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\"] = 'Set3'\n", "\n", "# Fix np.random seed for replicability\n", "np.random.seed(0)" ], "execution_count": 71, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "2VlVOyoMkW4g" }, "source": [ "## Choice of kernel\n", "\n", "In class, we discussed about the different kernels that we can use to define a GP. A nice example kernel is given in Eq. 6.63 in Bishop\n", "\n", "\\begin{align}\n", " k_{\\theta}(x, y) = \\theta_0 \\exp\\left( -\\frac{\\theta_1}{2} \\| x - y \\|^2 \\right) + \\theta_2 + \\theta_3 x^T y\n", "\\end{align}\n", "\n", "with non-negative hyper parameters $\\theta_0, \\theta_1, \\theta_2, \\theta_3$\n", "(At least one of $\\theta_0, \\theta_2, \\theta_3$ has to be positive for $K$ to be positive definite).\n", "Below we will implement this kernel." ] }, { "cell_type": "code", "metadata": { "id": "glDYa-7XvpSV" }, "source": [ "class MyKernel:\n", " \n", " def __init__(self, theta, bounds=None):\n", " self.theta = theta\n", " self.bounds = bounds\n", " \n", " def __call__(self, X, Y, eval_gradient=False):\n", " '''\n", " Compute kernel matrix for two data vectors\n", " \n", " Parameters\n", " ----------\n", " X, Y : 2-D numpy array\n", " Input points. X[n, i] (resp. Y[n, i]) = i-th element of n-th point in X (resp Y).\n", " eval_gradient : bool\n", " If True, return gradient of kernel matrix w.r.t. to hyperparams\n", " \n", " Returns\n", " ----------\n", " K : 2-D numpy array, shape = (len(X), len(Y))\n", " Kernel matrix K[i, j]= k(X[i], Y[j])\n", " gradK : 3-D numpy array, shape = (len(self.theta), len(X), len(Y)), optional\n", " Gradient of kernel matrix. gradK[i, m, n] = derivative of K[m, n] w.r.t. self.theta[i]\n", " Returned only if return_std is True.\n", " '''\n", " \n", " tmp = np.reshape(np.sum(X**2,axis=1), (len(X), 1)) + np.sum(Y**2, axis=1) -2 * (X @ Y.T)\n", " K = self.theta[0]*np.exp(-self.theta[1]/2*tmp) + self.theta[2] + self.theta[3]*(X @ Y.T)\n", " \n", " if not(eval_gradient):\n", " return K\n", " else:\n", " gradK = np.zeros((len(self.theta), len(X), len(Y)))\n", " gradK[0] = np.exp(-self.theta[1]/2*tmp)\n", " gradK[1] = -self.theta[0]/2*tmp*np.exp(-self.theta[1]/2*tmp)\n", " gradK[2] = np.ones((len(X), len(Y)))\n", " gradK[3] = X @ Y.T\n", " return K, gradK" ], "execution_count": 4, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "QACicVMVw0xj" }, "source": [ "from scipy.special import expit as sigmoid\n", "\n", "def W_mat(a):\n", " \"\"\"Helper to compute matrix W\"\"\"\n", " r = sigmoid(a) * (1 - sigmoid(a))\n", " return np.diag(r.ravel())\n", "\n", "def posterior_mode(X,t,K,max_iter=100, tol=1e-9):\n", " \"\"\"\n", " Computes the mode of posterior p(a|t) (via iterated Newton-Raphson)\n", " Parameters\n", " ----------\n", " X : 2-D numpy array\n", " Array representing training input data, with X[n, i] being the i-th element of n-th point in X.\n", " t : 1-D numpy array\n", " Array representing training label data.\n", " optimize_hparams : bool\n", " If True, optimization of kernel hyperparameters is performed.\n", " \"\"\"\n", " a_h = np.zeros(np.shape(t))\n", " I = np.eye(X.shape[0])\n", " for i in range(max_iter):\n", " W = W_mat(a_h)\n", " Q_inv = np.linalg.pinv(I + W @ K)\n", " a_h_new = (K @ Q_inv).dot(t - sigmoid(a_h) + W.dot(a_h))\n", " a_h_diff = np.abs(a_h_new - a_h)\n", " a_h = a_h_new\n", " if not np.any(a_h_diff > tol):\n", " break\n", " return a_h\n", " \n", "class GPClassification:\n", " \n", " def __init__(self, kernel):\n", " self.kernel = kernel\n", " \n", " def fit(self, X, t, optimize_hparams=False):\n", " '''\n", " Parameters\n", " ----------\n", " X : 2-D numpy array\n", " Array representing training input data, with X[n, i] being the i-th element of n-th point in X.\n", " t : 1-D numpy array\n", " Array representing training label data.\n", " optimize_hparams : bool\n", " If True, optimization of kernel hyperparameters is performed.\n", " '''\n", " self.X_train = X\n", " self.t_train = t\n", " if optimize_hparams:\n", " bounds_full = self.kernel.bounds\n", " result = minimize(x0=self.kernel.theta,\n", " fun=lambda x : neglog_evidence(theta=x, kernel=self.kernel, X=self.X_train, t=self.t_train, return_grad=True), \n", " jac=True,\n", " bounds=bounds_full,\n", " method = 'TNC',\n", " tol=eps)\n", " print(result.message)\n", " self.kernel.theta = result.x\n", " \n", " def predict_y(self, X_test,nu = 1e-5):\n", " ''' \n", " Parameters\n", " ----------\n", " X_test : 2-D numpy array\n", " Array representing test input data, with X_test[n, i] being the i-th element of n-th point in X.\n", " \n", " Returns\n", " ----------\n", " mean : 1-D numpy array\n", " Array representing predictive mean.\n", " '''\n", " K = self.kernel(self.X_train, self.X_train) + nu * np.eye(self.X_train.shape[0])\n", " k = self.kernel(self.X_train,X_test)\n", " c = self.kernel(X_test,X_test)\n", "\n", " y_lap = posterior_mode(self.X_train,self.t_train,K)\n", "\n", " pred_mean = k.T @ (self.t_train - sigmoid(y_lap))\n", "\n", " Winv = np.linalg.pinv(W_mat(y_lap))\n", " pred_cov = c - k.T @ np.linalg.pinv(K + Winv ) @ k\n", "\n", " return pred_mean, pred_cov\n", "\n", "\n", " def predict_prob(self, X_test):\n", " K = self.kernel(self.X_train, self.X_train)\n", " y_lap = posterior_mode(self.X_train,self.t_train,K)\n", "\n", " k = self.kernel(self.X_train,X_test)\n", " pred_mean = k.T @ (self.t_train - sigmoid(y_lap))\n", "\n", " Winv = np.linalg.pinv(W_mat(y_lap))\n", "\n", " # computing only the diagonal elements of Ker(x_test,x_test) for scalability\n", " c = np.array([kernel(np.asarray([x]),np.asarray([x])) for x in X_test]).ravel()\n", " \n", " pred_var = c - np.diagonal(k.T @ np.linalg.pinv(K + Winv ) @ k)\n", " kappa = 1.0 / np.sqrt(1.0 + np.pi * pred_var / 8)\n", "\n", " return sigmoid(kappa * pred_mean) \n", "\n", "\n", "def plot_probs(pred_probs, Xtest, ax=None):\n", " ''' \n", " Parameters\n", " ----------\n", " pred_probs : 1-D numpy array\n", " Array containing the predicted probabilities\n", " ax : Axes for plotting\n", " '''\n", " if ax is None:\n", " fig = plt.figure(figsize=(16,10))\n", " ax = fig.add_subplot(111)\n", "\n", " ax.hlines(0.5,np.min(Xtest),np.max(Xtest),'b')\n", " ax.plot(X, t,'og',markersize=10,label='training data')\n", " ax.plot(Xtest, pred_probs, label='predicted probs')\n", " ax.plot(Xtest, sigmoid(truef(Xtest)), ':r', label='ground truth probs')\n", "\n", " ax.set_xlabel(r'$x$')\n", " ax.set_ylabel(r'$t$')\n", " ax.legend()\n", "\n", "def plot_GP(pred_mean,pred_cov,Xtest,ax=None):\n", " ''' \n", " Parameters\n", " ----------\n", " pred_mean : 1-D numpy array\n", " Array containing the predicted mean\n", " pred_cov : 2-D numpy array\n", " Array containing the predicted covariance\n", " ax : Axes for plotting\n", " '''\n", " if ax is None:\n", " fig = plt.figure(figsize=(16,10))\n", " ax = fig.add_subplot(111)\n", " \n", " ax.plot(Xtest, pred_mean, label='predictive mean')\n", " ax.plot(Xtest, truef(Xtest), ':r', label='ground truth')\n", " pred_stddev = np.sqrt(np.diagonal(pred_cov))\n", " upper_cb = pred_mean + 2*pred_stddev\n", " lower_cb = pred_mean - 2*pred_stddev\n", " Xtest = Xtest.ravel()\n", " ax.fill_between(Xtest, upper_cb, lower_cb, alpha=0.2)\n", " ax.set_xlabel(r'$x$')\n", " ax.set_ylabel(r'$t$')\n", " ax.legend()" ], "execution_count": 90, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Y9AuObZvno3t" }, "source": [ "Empirical Bayes (and any general model comparison) also needs the marginal likelihood. Indeed, one way to think about the `cost' of a given vector of hyperparameters is in terms of the *negative log of the marginal likelihood*. Given dataset $D=\\{(x_i,t_i)\\}_{i\\in[N]}$, and with $y_{L}$ denoting the posterior mode, under the Laplace approximation we have (all logarithms in base $e$):\n", "\n", "\\begin{align}\n", "-\\log p(y|D) = \\frac{1}{2}\\log|I + K_DW|+ \\frac{1}{2}y_L^T(K_D)^{-1}y_L - y_L^Tt + \\sum_{i=1}^N\\log\\Big[1+e^{t[i]y_L[i]}\\Big]\n", "\\end{align}\n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "G6Vg4OXCwwPi" }, "source": [ "def neglog_evidence(theta, kernel, X, t, return_grad = False):\n", " '''\n", " 'cost' of hyperparameters is the (negative of the log of) the marginal likelihood\n", " \n", " Parameters\n", " ----------\n", " theta : 1-D numpy array\n", " kernel parameters\n", " kernel : kernel object\n", " the kernel function\n", " X : 2-D numpy array\n", " input data, X[n, i] = i-th element of n-th point in X.\n", " t : 1-D numpy array\n", " label data\n", " \n", " Returns\n", " ----------\n", " val : float\n", " value of the cost function\n", " '''\n", "\n", " kernel.theta = theta\n", " K, gradK = kernel(X, X, eval_gradient=True)\n", " Kinv = np.linalg.pinv(K)\n", " y_lap = posterior_mode(X, t, K)\n", " W = W_mat(y_lap)\n", " C = np.eye(len(y_lap)) + K @ W\n", "\n", " val = 0.5*(y_lap.T @ Kinv @ y_lap + np.linalg.slogdet(C)[1]) - t @ y_lap + np.sum(np.log(1.0+np.exp(y_lap)))\n", " \n", " if not(return_grad):\n", " return val\n", " else:\n", " grad = np.zeros(len(theta))\n", " Cinv = np.linalg.pinv(C)\n", " sig_lap = sigmoid(y_lap)\n", " for cnt in range(len(theta)):\n", " # Using the formula for the gradient from Bishop 6.4.6\n", " z = Cinv @ gradK[cnt] @ (t - sig_lap)\n", " grad[cnt] = 0.5 * np.trace(Cinv @ W @ gradK[cnt]) - 0.5 * y_lap @ Kinv @ gradK[cnt] @ Kinv @ y_lap \\\n", " + 0.5*np.sum(np.diagonal(Cinv @ K) * sig_lap * (1.0-sig_lap) * (1.0-2.0*sig_lap) * z)\n", " return val, grad" ], "execution_count": 104, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "AMfiUPCYyWRp" }, "source": [ "## GP regression with fixed hyperparameters " ] }, { "cell_type": "markdown", "metadata": { "id": "JRqWIs4Tzou1" }, "source": [ "Generating the ground truth and data" ] }, { "cell_type": "code", "metadata": { "id": "0uQCN5O7yrwq", "colab": { "base_uri": "https://localhost:8080/", "height": 395 }, "outputId": "ca781aea-adf8-4922-de7d-d68dc5874ad9" }, "source": [ "def truef(x):\n", " return 5*np.sin(x**2 *0.5) \n", " #return 5*(x-2)**2 - 5\n", "\n", "\n", "x_min = 0\n", "x_max = 5\n", "\n", "\n", "N = 200\n", "\n", "X = np.random.uniform(x_min,x_max, N)\n", "\n", "t = np.random.binomial(1,sigmoid(truef(X)))\n", "\n", "\n", "# Generate grid for plotting\n", "Xcont = np.linspace(x_min,x_max,200) \n", "\n", "plt.figure(figsize=(8,6))\n", "plt.plot(X, t,'og', label='training data')\n", "plt.plot(Xcont, sigmoid(truef(Xcont)), '-r',label='ground truth')\n", "plt.xlabel(r'$x$')\n", "plt.ylabel(r'$t$')\n", "plt.legend()\n", "plt.show()\n", "\n", "X = np.reshape(X,(len(X),1))\n", "Xtest = np.reshape(Xcont,(len(Xcont),1))" ], "execution_count": 42, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAF6CAYAAAD1UEqsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9b3/8VdICBBC2CIoslURLVSDC1UQi4pbpS6U29uiRe1VtGgX3Kr1Xi5KN6311pZae6WUKrdQawsVLy7tz9qKuSpBaVhSAS1BRRbZsrCFJOf3x5lkvjOZSWYyM+d75sz7+Xh8H/mek2/O+cz6yfd7zvmePMBBREREAqWL7QBEREQk/ZTgRUREAkgJXkREJICU4EVERAJICV5ERCSAlOBFREQCqMB2AOm0a9cutm7dajsMERERzwwbNowBAwa0WR+oBL9161bGjh1rOwwRERHPVFRUxFyvIXoREZEAUoIXEREJICV4ERGRAFKCFxERCSAleBERkQBSghcREQkgJXgREZEAUoIXEREJICV4ERGRAPJ0JrvbbruNG264gVNPPZUlS5bwla98JW7bWbNmcc8991BUVMTvf/97Zs6cSUNDg4fRRjkVmAT0BmqAl4F1adrWRuDk0LID5EW1Pwi8CJQBJxrrjwLPGXGY223G/fftCFAY2mYzsAU4Digytr2+g/07oX0VEn7sg4GzQvtoBlYDLxh/81nj9w7QAHQz4mp53J9qJ5aWtgdDvy8i8rmP95rEeh4OAvmhGMzH9c+o5yP6cecZ2wa4LE5botpF7z/6dTbXx9t2y+tuvs9iPbZY27qSyE/3e6Gf5vsn+nHGEisGk/k6t7wPPox6LO1t49QO2qbzc+c3iTy29t7jsd6L5mfN/OzH+n30697R58Fs1/K5iP4MHwn9rjC0HP2ZrgF2A58g/J7ZApQS+/sw217z9l5TC+/lPNyXyxNTpkyhubmZSy+9lB49esRN8JdccglPPfUUF154IR999BHLli3jjTfe4Nvf/na726+oqMjMVLWnAlcQftOC+yExk2sq22rvC7ZFc6hNdLtmYFmoHr3dWOIl8I72b2rC/XBGfzmswk3ynwU+ncA2OxtLA7AGOJ22r0ms9cnGEEtjqF1+Au0g9r/O8fbV3rYbgWcJf6l39Bo34r420WNzLZ/yZF7nWDGYYr3OTqhE7z/WNk4FrqLtc9XSFtL3ufObRL5T4rVZA5xJZrpnyX4XJNI+1j8SySxny2ve3mtKO79Lw+OKl/s8HaJftmwZzz77LHv27Gm33fXXX8+CBQuoqqpi//79fOc73+GGG27wJshYJtH2S7UwtD4d20rkAxWdUM31k+JsN5ZY20j2Sz8/xt/k4fbkCP1MZJudjaUwtI9Yr0ms9cnGEEsBHSf3lnbxvnjj7au9bRcQfp8l8hoXEPtTHeufw0SZMZhivc55cfYfaxuTiP1ctbRN5+fObxJ5bPHanEXmxl6TfY905nOe7HK2vObtvaaW3su+PAY/evRoKisrW5crKys59thj6devX5u2M2bMoKKigoqKCkpLSzMTUO8k13dmW6nonaHtJqtL1E8v9pXoei840KUZujZCt6NQ1OCWwkbIa05hu72jftoQa9/JPtfR22jv8bT3nvbDez1ViTy2eG08fI/nNUOPBiho8m6fMWXDa97ea2rpvezLu8kVFxdTU1PTutxS79WrF3v37o1oO3/+fObPnw/Ev6NOymqAPnHWp2tbqWiJI93bTVaz8TOR3m6q+4q1jxT23eswDKmFITUw4AD0O+SW/ofC9ZIjUHQUehx1f7aU7gl8ATbmwdF8ONoFjhRAbTfY3x1qukFNd7e+vzt81As+LHHLByXwgRN6ajPx3klUrPd6ss919Dbaezztvac787nzm0S+U+K1SePnq9dhmPA+nL4DTt0Jw2pgcC0UN7jvcfN93ZgHe4pgRzFs6QNVx8Bbg+Cvw2Fve8fs0yEbXvOOXlML72VfJvj6+npKSkpal1vqdXV1dgJ6mdjHT16O3TzpbaV6DL4lDtvH4FeH6qvx5zF4x/3y+uTH8Mnd7s9hNW5CH1ILvY90sN8UFThQ0Ag9QrEec7Cjv3AdzIMqYP2vYd2ZsGooVAyCI11jNM7UMfhY7/VYr3N7x+Cjt/Ey8Y/Bx3tPd/Zz5zeJfKfEa5PiMfhuR+GadfDltXDe+9A1wdGlAgcGHnBL2U64eqO7vhl4fQjMPwN+NxoORX/20nEMPhte845eUwvvZV8m+A0bNlBWVsYzzzwDQFlZGTt27GjTe/dMy0kQ6TgDMta20nUWvbldm2fRt/zM9Fn0HxL7NQmtHwx8+kP49Ecw9kM4c2fmk3hTKOymfGgOPY8FzVCYwhB9keM+lWftp/UL4XA+rBwKy0+GZSfDtj54fxZ99Ouc7Fn0LcsdtQ3iWfSJfKe01yb6OW7Rzln0+c0w63W4pxyOOZR4qIcKoGuTm+Bj6QKc+4FbfvhnuGMS/Ob00H5z6Sz6RF7TIJ9Fn5+fT0FBAXPmzGHw4MHMmDGDxsZGmpoixzcvvfRSfv3rX7eeRb906VJWrVpl7yx6yQolwAXAJcDFwElJ/v0h4INQ2QHsAfaGSku9FjgQanvQKA2Ej1DEk4/7fdsV9/u35dBcH+NnP+B43H9OhgDDgYEdbLcJ+BPw38ByPPxAS9Y4DVgInBHjd2uA14C/A5tx3/81uO9r8//hQtw8PAgYiXvS+ERgLG17in8GvoT7mZHMay/3OV6VOXPmONHmzJnjDBkyxKmrq3OGDBnS2vb22293duzY4dTU1Di/+tWvnMLCwg63X1FR4dljUfFH6Q3O9eC8CM5RcJwOyh5wVoLzBDi3g3MlOKeDU+qDxxKv9AdnIji3gfNLcN5p5/GtBedffBCzin/KheDUE/k+qQbnLnCOS8P2S0Pbej9qH5XgHOODx58LpZ3cZz84Dx6kSoBKHjiTwVkKzmHiJ7sD4PwFnB+AMwWcwT6IPV3leHBuAefPcR77S+Cc6IM4VeyWy8E5RPh9cRCcb4FTkIF99QTnR+A0GfvbgPsPqu3nIehFCV4l60sRODNpvwe7Gpzvg3M+OIU+iNmLMgych8Cpo+0/ONf6ID4VO+WTuAm95f3wQWhdpvd7LTiNxn5/64PnIuhFCV4la0tPcOaAs5f4Sf1OcIb4IFabpR84Pybyy9UJrevig/hUvCtdwXmb8HvgPXCGe7j/LxH5HrzKB89JkIsSvErWlXxwZoCznbZJfT84D4Mz0gdx+q2cAU5V1PP1a9xDG7ZjU/GmfN947Q+BM9pCDL8yYvgInD4+eF6CWuLlPl/OZCdyKbAWeAI41lj/LvB13LPM7wY2eR+a770NnA0sNdZdj/tcdna2WskeI4FvGcv3AhssxHEHsD1UPy60LN5Sghdf6Q88hXsp9Chj/Ye4Sepk4GdAvfehZZU64F+A+ca6m4D2LzSVILib8ER3fwV+aimO/UQm9Rm4l4iKd5TgxTcuxp0TY7qxrha4D7dX8hQdX2suYQ5wC+71zy3mAudbiUa8MAi4zliejfs+sOUZYFuofixwtcVYcpESvFiXDzyEO1mLORz/P8AI4Ae4E8tI8hzcntPfQsv5wBI6njxHstPthGdDLcedwMamJiJHkWbaCiRHKcGLVaW4id08ZrgD+BxuT/5jG0EFTBMwDdgZWj4WeMReOJIhJbgjNi0etBVIlPm4txcAd6bJT1qMJdcowYs1pwAVwIXGuhdwp9VcYSWi4NpO5KGPa4FzLcUimTEZ6BWqV+Gfz9BHuFMot5ger6GknRK8WDEB+D/cudbBPbY+G/dLSr32zPgz8DtjeR76AgiSKUb9f7B77D3aU0Z9krUoco8+3+K5ybjJpm9ouR73rqHfxV9fSkF0F+Gb850O3GAvFEmjbsBnjeU/2gokjr/iHioC9063ve2FklOU4MVTU3Cvz+4eWt4BfAb4X2sR5ZYPcE9abHE3ujY+CCYBxaH6JuAfFmOJpQZ3fgZwT/T8jMVYcokSvHjmStwh4pazfN8DxuHerlK88yjuNcrgngfxOYuxSHqYw/PLrEXRvr8Y9QvjtpJ0UoIXT0zCTe4t943eiHsv6WpbAeWwetx7x7e4y1YgkhZdcP95buG34fkWZoLXcXhvKMFLxp2J+6XTLbS8GXeylW3x/kAy7qdAQ6j+GeDTFmOR1JwFDAjVtwNvWoylPa8BR0P1UwnHLJmjBC8ZNRT3+HrL8cEPgItwj72LPR8Bi41lTUCSvc4x6v8P/56oehB4w1g+31IcuUQJXjKmBHie8Ox0e3Gno33fWkRi+rlR/zzhEx8lu5ijL37tvbcwh+nPsxZF7lCCl4zIw70Wd3RouQF3HuqN1iKSaBW4h0vA/WdsssVYpPPONuqrrEWRGDO+0XFbSboowUtG/CdwhbH8b8BKS7FIfOYw/TXWopDO6od7vwaAI0ClxVgSUWXUR8VtJemiBC9p91ngfmP5YeA3dkKRDpgJfjKagCTbmMPzfyd84qRfbSU80dJA3NtDS+YowUtaDSJyWsr/h+5B7mebgNWhejdgqsVYJHnZdPwd3BMAzUl4dOOZzFKCl7TpgttTLw0tb8O9i1lT3L8QPzB78VdZi0I6I5uOv7fQML13lOAlbe4lfOlLE+4x3d3WopFEPWfULwC62gpEkpZtPXhQgveSErykxelEHnefC7xqJxRJ0rvAllC9F5HXVYt/nUB4tGwv7uuYDZTgvaMELynrBiwi3PP7P+B79sKRTviTUb/EWhSSjDKjvjpuK/9RgveOEryk7D8JX9N6ALgeHXfPNi8Z9UutRSHJGGnU/Xb3uPZsAQ6H6sejKzcySQleUjIa95ajLe4me4YKJewvRN6vW5cv+d9JRn2TtSiS10TkhFc6kz5zlOCl0/KAXxAemn8ttCzZp4bwSVpd0N2+soHZg8+mBA8apveKErx02o3AhFD9KHAL/r3RhXTMHKa/2FoUkqigJHhNWZs5SvDSKccADxnLDxP5oZXs84pR15n0/tYbdyY4cI9nf2Axls4w/yH5hLUogk8JXjrlEdx5sAHeA75rMRZJj7eAxlB9FO4lc+JP5vH3d8m+kbMPjfrx1qIIPiV4SdqFwHRj+VbgkKVYJH0OAutD9S7AWRZjkfZl8/A8RCb4wdaiCD4leElKF+BRY3kJkddQS3YzZ0PTML1/ZesZ9C22A82h+rFAgcVYgkwJXpLyb8CpoXo9cIfFWCT93jDqZ8dtJbZlew/+KLArVO+Cm+Ql/ZTgJWHFwHeM5YeAHZZikcwwe/BK8P6V7QkeNEzvBSV4SdjdhP/T/hD3RDsJlndwr4kH97UeajEWic9M8JutRZGabUZdJ9plhhK8JOR44C5j+d/RiXVB5AAVxrKOw/vPQKAkVK8hPNSdbdSDzzwleEnId4GiUP1t3JvLSDDpOLy/BWF4HnSpnBeU4KVDY4DrjOW7yL7rbiVxq4z66daikHiGG/X3bAWRBuYQvXrwmaEELx16kPAb5TkiZzyT4Flr1E+N20psMXu7H8Zt5X8aos88JXhp17mEbx/aBHzLYizijfeBulC9FBhgMRZpy0zw2+K28j+dZJd5SvDSrrlGfRHuWdYSbA6wwVj+lK1AJCazt5vNPXgl+MxTgpe4zsedlhbcOcq/E7+pBMx6o64E7y9B6cEfAPaH6t1wR4skvZTgJS6z974Q+KetQMRzZoLXcXh/CcoxeNBx+ExTgpeYLgbOC9Ub0N3ico168P5UQHiyqWbcOd2zmRJ8ZinBS0xm7/2XuCdeSe4wE/xoIM9WIBLhWMJf2rsI3943W+k4fGYpwUsblxOeweww8H2LsYgdO4HdoXovNGWtXwRpeB7Ug880JXhp436j/t9k94k80nnrjLqG6f3BTIJB+FxqspvMUoKXCBcBY0P1w7iT3Ehu0nF4/wnKGfQtzHMINN9C+nma4Pv27cvSpUupr6+nurqaadOmxWxXWFjI448/zo4dO9izZw/Lly9n0KBBXoaas+4z6r9Et4PNZTqT3n+CNkS/x6j3txZFcHma4B977DEaGhoYOHAg1157LY8//jijRo1q0+6b3/wm48aN47TTTmPQoEHs27ePefPmeRlqThoHXBCqHwUethiL2Fdl1EfGbSVeCtoQvZngdR18+nmW4IuKipg6dSqzZ8/mwIEDlJeXs3z5cqZPn96m7Sc+8Qleeukldu3axZEjR3j66acZPXq0V6HmLLP3/ht05nyue9eon2QtCjEFrQe/26irB59+niX4kSNH0tjYyObNm1vXVVZWxkzcCxYs4Nxzz+W4446jR48eXHvttbzwwgtehZqTyoDPherN6Ni7uIdnDoTqfYB+FmMRV9COwe/H/b4B9z2WbzGWIPIswRcXF1NbWxuxrqamhl69erVpu3nzZj744AM++ugjamtr+eQnP8ncuXPbtAOYMWMGFRUVVFRUUFqqQZ7O+rZR/wOw0VYg4itmL36EtSikRdCG6JuBfcay/olML88SfH19PSUlJRHrSkpKqKura9P2scceo1u3bvTr14+ePXuydOnSuD34+fPnM3bsWMaOHcvu3btjtpH2nQR8wVjWde/SQgneP/oB3UP1GqDeYizpZH5rq4uWXp4l+E2bNlFQUMCIEeGvibKyMjZs2NCm7ZgxY/j1r3/Nvn37aGhoYN68eZx99tn076+jNJlwD+E3wvPA3y3GIv6iBO8fQeu9t9CZ9JnjWYI/ePAgS5cuZe7cuRQVFTF+/HiuuuoqFi1a1KZtRUUF1113HSUlJRQUFHDrrbeybds29uzZE2PLkorjAPM0R/XexWQm+BOtRSEQvBPsWijBZ46nl8ndeuut9OjRg127drFkyRJmzpxJVVUVEyZMiBiqv+uuuzh8+DCbN2/m448/5vLLL2fKlClehpozvgEUhurloSLS4j2jrh68Xcca9Wy/yYxJl8plToGXO9u3b1/MRP3aa69FnGy3d+9evvzlL3sZWk4qBr5qLOu6d4mmIXr/MGd622UtivTTpXKZo6lqc9hNuJemAGwClluMRfzpQ9wpi8FNMCXttJXMOsaof2wtivTTEH3mKMHnqAJglrH8COBYikX8ywH+aSzrOLw9uZDgNUSfXkrwOeoLwLBQfRfwlMVYxN80TO8PQU3wGqLPHCX4HHWXUf8Z4WFYkWhK8P4Q1ASvIfrMUYLPQRcCZ4TqB4GfW4xF/E8J3h+U4CVZSvA5yOy9LyTyAyYSzbxUTsfg7cmFBK9j8OmlBJ9jPgV8NlRvBn5sMRbJDluM+lBrUeS2olAB93BaUKaphcgE3w/IsxVIACnB55g7jfpSIntnIrF8YNQHoy8NG4Laewc4CrTchiwf6G0xlqDRZzWHDAKuMZZ/ZCsQySqHCU+s0pXIGdXEG2aCD+IttTRMnxlK8DnEnJZ2JfCmxVgku7xv1DVM770g9+BBl8plihJ8juhF5LS06r1LMswEPyxuK8mUoCd4nUmfGUrwOeJGwse2NgLPWYxFso968HblUoLXEH36KMHngHwip6X9LzQtrSRHCd6uXErw6sGnjxJ8DvgXwsOqH6NpaSV5SvB2BT3B6xh8ZijB54C7jfpjaFpaSd5Wo64E772gJ3j14DNDCT7gJgJnhuqH0LS00jnqwduVSwlex+DTRwk+4MyJbZ4kmF8OknkfEx756QcUW4wlFwU9we816n2tRRE8SvABdgpwhbGsaWmlsxwiZ7QbYiuQHBX0BF9j1DWTXfoowQfY7Ub9WWCTrUAkEDRMb0chUBKqNwL7LcaSKUrwmaEEH1ADgOuM5UdsBSKBoQRvR/Q0tUG8xNVM8CVxW0mylOAD6jage6i+CndqWpFUKMHbEfTheVAPPlOU4AOoB3Crsazeu6SDErwduZDgD+EefgC3Y1LYTltJnBJ8AF1P+FKTauAP9kKRADGvhdd89N7JhQQP6sVnghJ8wHQh8uS6HwNNlmKRYPnQqB9vLYrc08+oB/FWsS2U4NNPCT5grgBGhur7gV9ZjEWC5SOjPshaFLnHvC58n7UoMk8JPv2U4APGnNjmF0C9rUAkcOoIv5+K0JewV/oY9SBeItdCCT79lOAD5GzgvFD9KDDPYiwSTNuMunrx3lAPXjpLCT5AzN77YiKHVEXSwXxP6Ti8N3KxB69r4dNDCT4gTgA+byz/l61AJNB0HN576sFLZynBB8S3gPxQ/U/AWouxSHApwXsvF3vwSvDpoQQfAMcCNxjLP7AUhwSfErz3cqUHX2vUleDTQwk+AG4HuoXqbwJ/tReKBJwSvPfUg5fOUoLPcn2Amcayeu+SSUrw3upCONk1E5kEg0YJPv2U4LPcbUCvUH0DsNxiLBJ8SvDeMnvvtQTzTnItlODTTwk+i/UAvmksP0SwvwDEvu1GfRCQZyuQHGEm+CAffwcl+ExQgs9iNxG+EUU1sMReKJIjDhFONF0J39RIMsM8wS7Ix99B18FnghJ8lioA7jKWf0T4dosimaRheu+oBy+pUILPUtcTvif3LnRTGfGOpqv1Tq5cIgdK8JmgBJ+FugL/biz/F+7QqYgX1IP3Tq5cIgfujYyaQ/WeuKOUkhol+Cx0HfCJUH038DOLsUjuUYL3Ti714CFyshsdh0+dEnyWie69PwwcsBSL5CYleO/kUg8eNEyfbkrwWSa69/6YxVgkNynBeyfXevBK8OmlBJ9F1HsXP9AtY72jHrykQgk+i6j3Ln5gTnZzrLUockMu9+B1DD51SvBZogD13sUfdhr1AWg2u0xSD15SoQSfJWag3rv4wxHCX8RdiexlSnrlcg9eCT51SvBZoBiYYyw/hHrvYpfZix9oLYrgUw9eUqEEnwXuJPwl+gG67l3sU4L3hnrwkgoleJ8bSOSc87OBw5ZiEWmhBJ95RUBhqH4Y99BI0JkT3SjBp04J3udm4w7RA6wFFlmMRaSFmeB1Jn1m5NKNZlqoB59enib4vn37snTpUurr66murmbatGlx255++un87W9/o66ujh07dvCNb3zDw0j9YQRws7F8L+G5mkVsUg8+83LpVrEtlODTy9P5/B977DEaGhoYOHAgY8aMYcWKFVRWVlJVVRXRrn///rz44ovcfvvt/P73v6ewsJDBgwd7GaovfB/3LGWAV4AXLMYiYlKCz7xc78HrOvjUedaDLyoqYurUqcyePZsDBw5QXl7O8uXLmT59epu2d9xxBy+99BKLFy+moaGB+vp63nnnHa9C9YVzgC8Yy/fYCkQkBiX4zMu1E+wA6ox6L2tRBIdnCX7kyJE0NjayefPm1nWVlZWMHj26TdtzzjmHvXv3Ul5ezs6dO1m+fDlDhgyJud0ZM2ZQUVFBRUUFpaWlGYvfS12IvM79d0CFpVhEYlGCz7xcu0QOlODTzbMEX1xcTG1tbcS6mpoaevVq+zIOHjyY66+/nm9+85sMHTqULVu2sGTJkpjbnT9/PmPHjmXs2LHs3r07I7F77RbgjFD9EPAti7GIxKIEn3m52IPX7WLTy7Nj8PX19ZSURL5kJSUl1NXVtWl76NAhli1bxurVqwF44IEH2LNnDyUlJW3+SQiaY4DvGcvfA7ZaikUknujpaiX91IOXVHnWg9+0aRMFBQWMGDGidV1ZWRkbNmxo03bt2rU4jtO6bNaD7kHC/7lvBn5kMRaReA4S/jLuRmQykvTIxR78EeBoqF5IeB4A6RzPEvzBgwdZunQpc+fOpaioiPHjx3PVVVexaFHbK7sXLlzIlClTKCsro6CggNmzZ7Ny5crA997PAf7NWP46uTG5hWQnDdNnVi724EG9+HTy9Dr4W2+9lR49erBr1y6WLFnCzJkzqaqqYsKECRFD9a+88gr33XcfK1asYNeuXYwYMYJrrrnGy1A91wX4ubG8FHjJUiwiiVCCz6xc7MGDjsOnmxOUUlFRYT2GzpZZ4DihcgCcoT6ISUWlvfIHwu/ZL/ggnqCVV4zn9wIfxONVWWs87tN8EE82lHi5T1PV+sCJtD2x7n1LsYgkSj34zMrVHryG6NNHCd6yPGAB7o0lACqBh+2FI5IwJfjMytVj8BqiTx8leMu+AUwM1RuBrxA+i1TEz5TgM0s9ePXgU6UEb9GncC+La/EgsMZSLCLJUoLPnHzCvddmInu1QacEnz5K8JZ0AxYD3UPLbwNz7YUjkjTdMjZzzDup1eCeMZUrlODTRwnekkeAU0P1g8C1aGhesot68JmTq8PzoGPw6eTp7WLFNQ24zVi+E8ite+VJECjBZ45fTrDr27cvs2bNYvjw4eTl5XmyzxLg1VB9LPCUJ3v1P8dxqK6u5tFHH2XfvsT+7VOC99goYL6x/DTwC0uxiKSiHnf0qQj3UFMJuXWsOJP80oOfNWsWq1evZu7cuTQ1NXmyz2OAoaH6x+iS4Rb5+flMnjyZWbNmMWfOnIT+RkP0HuoPPAf0DC2/A9xkLxyRlKkXnxl+6cEPHz6c559/3rPkDmDuKd+zvfpfU1MTK1asYPjw4Qn/jRK8Rwpxp589IbRcD/xL6KdIttph1JXg08cvPfi8vDxPkzsowbenqakpqUMlSvAeyAMWAp8JLTcD1wBt76Mnkl3Ug88Mv/TgbTATvBJUavT8eeCnuAm9xb24Q/Ui2U4JPjP80oO3odmo9+ndm5kzZ3ZqOytWrKB3797ttnnggQeYNGlSp7bfnuuvv5558+a122bixImMGzcu7fs26SS7DPsB8DVj+edoKloJDiX4zMjaHvypwCTcC/lrgJeBdcltwuzB9+nTh1tvvZXHH3+8Tbv8/Px2Dx9Mnjy5w30lerJaJpx//vnU19fz+uuvZ2wf6sFn0H/h9tZbLCEy2YtkOyX4zMjKHvypwBW4/53khX5eQXjCjwSZKfv2Bx/kxBNPZM2aNfzwhz9k4sSJvPrqqzz77LNUVVUBsGzZMlavXs369euZMWNG699u2bKF/v37M2zYMKqqqnjiiSdYv349L730Et27u1OMLVy4kKlTp7a2v//++3nrrbdYu3YtJ598MgClpaX86U9/Yv369cyfP5/q6mr69+/fJu4bbriBjRs38uabb3Luuee2rv/c5z7HG2+8wdtvv82f//xnBgwYwLBhw/jqV7/K7bffzpo1a5gwYULMdgUc70cAACAASURBVOlg/VZ36Sp+uV1sV3AWEL7loQPOH0PrbcemopLOMpXwe3yZD+IJSnnBeF4/azGOp556KvH2s3C4P0aZldw+u4BzZqhcPmyYs27dutbfTZw40amvr3eGDx/euq5v374O4HTv3t1Zt26d069fPwdwtmzZ4vTv398ZNmyYc/ToUaesrMwBnKefftq59tprHcBZuHChM3Xq1Nb2X/va1xzAmTlzpjN//nwHcObNm+fce++9DuBceumljuM4Tv/+/SNiPvbYY52tW7c6paWlTteuXZ3XXnvNmTdvngM4ffr0aW134403Oj/60Y8cwJkzZ45z5513tv4uXrtEXpN4uU9D9GnWB/g97ihVi2dwj8E3WolIJHPUg8+MrOzBxzvc3f5h8DbMY/CxhphXrVpFdXV16/I3vvENpkyZAsCQIUM46aSTePPNNyP+ZsuWLVRWVgLw1ltvxb3UbOnSpa1tPv/5zwMwYcKE1u2/9NJL7N27t83fnX322fz1r39l9+7dADz99NOMHDkSgMGDB/P0009z3HHHUVhYyJYtW2LuO9F2ydAQfRp9GvdmMWZy/xXuzHVK7hJESvCZkZXH4GuSXN+O9i7MO3DgQGt94sSJXHTRRYwbN44xY8awZs2a1uF305EjR8LbbmqioCB237alXXttkjVv3jx+9rOfcdppp3HLLbfEjC+ZdslQgk+DQmA28Bow3Fh/H3Aj7b9ZRbKZEnxmZGUP/mWgIWpdQ2h9klq+Mw/W1dGrV/xbzvTu3Zt9+/Zx6NAhTj75ZM4555zkd9aB8vJy/vVf/xWAiy++mH79+rVp8+abbzJx4kT69etHQUEBX/jCFyJi3LZtG+CeXd+iLuqxxWuXCiX4FF0K/B33TnBdQ+v2AVNwz6AXCbJa4HCo3pPwLI2SGjPBZ00Pfh3u9b/7cY8A7w8tJ3kWPYSH6Wv27uWN8nLWrVvHD3/4wzbtXnzxRQoKCqiqquLBBx/kjTfe6Gz0cT3wwANccsklrFu3ji984Qts376durq6iDY7duzg/vvv5/XXX6e8vJx//OMfrb+7//77eeaZZ1i9enXrED7Ac889x5QpU1pPsovXLlXWTyhJV/HqJLt8cK4EZyWRJ9I54JSDM8wHz4WKilelmvD7/0QfxJPtpch4Pg9ajiWpk+zSWE4hfKJdT8vPQWFhoZOfn+8AzjnnnOOsWbPGd6+JTrJL0QDgLOCzuL3z46N+Xwd8G3icyJNERIJuJzAsVB8IvGcxliDIyuH5NPPTbHZDhw7ld7/7HV26dKGhoSHiUjy/SyjBr1u3jnPPPZfa2ty4V1Qv4A+4d8fqjXt3o7ZXPboagCeA7wPbPYlOxF90HD69svIEuzQzO0m256N/9913OeOMMyxH0TkJ/XM0atQounXr1mZ9SUkJP/vZz9IelG1HgIuBs4FTiJ3cd+Am9RHA11Fyl9ylG86kl3rwuuFMurTbg3/++edZtWoVjuMwZMgQPv7444jfFxUVccstt/C1rwVrfrYG3CRv/ktTj3tzmJXA86GfuvRNRD34dFMPXgk+XdpN8OvWrWPixInk5eWxatUq6urqqKysZM2aNaxdu5ZTTjmF7duD2Xe9HDiIewlnDW4P3bEakYg/KcGnl3rwHU92I4lpN8Hfc889ABw+fJhx48YxaNAgxowZw5gxY5g8eTIFBQV861vf8iRQr/3FdgAiWUIJPr3Ug1cPPl0SOsmuuLiYxsZG1qxZw4oVKzIdk4hkESX49FIP3l8Jfs6cOdTX1/PII49ErL/qqqvYtGlTxDXviRg2bBjjx49nyZIlgDupzVlnncXXv/71tMXcIqHRj8ZGHW0WkdiU4NNLPfjkh+jz873/N+Dqq69m1KhRMX/XXjzDhw/nmmuuyVRYEXR4Q0RSogSfXurBR/bgv/Ef/8E777zDypUrWbx4MXfeeScAr7zyCj/+8Y+pqKjgm9/8JhdeeCFvv/02a9euZcGCBRQWFgLh28YCnHnmmbzyyiuA2zNfsGABr7zyCu+9915ED/q+++5j48aNrFy5svW2saZx48Zx5ZVX8vDDD7NmzRpOOOGENvGYt6IFWme/e/DBBznvvPNYs2YNs2bNAmDQoEG88MILbNq0iYceeihtz6MmuhGRlOwnfNVJL6AI9wRV6Ry/9uAzeZJxXtRyS4IfddZZXDZ1KmVlZXTt2pW3336bt956q7VdYWEhY8eOpVu3bmzevJlJkyaxefNmnnzySWbOnMlPfvKTdvd7yimncMEFF9CrVy82btzI448/zmmnncaXvvQlxowZQ0FBQZt9Arz++ussX76c//3f/+UPf/hDm3jAvdd8LPfeey933XUXV1xxBeAO0Y8ZM4bTTz+dI0eOsHHjRubNm8eHH36YwDPXPvXgRSRlu4y6evGpUQ8+nODLzj2Xvzz7LEeOHKG+vp7nnnsuot3TTz8NwMknn8yWLVvYvHkzAE8++SSf+cxnOtzPihUraGhoYM+ePezatYuBAwdy3nnnsWzZMg4dOkRdXR3Lly9POO6WeJL18ssvU1tby5EjR6iqqmLYsGEd/1EClOBFJGWa7CZ9zB58rib4RI/Bm7eOjaexsZEuXdytRN+CNdHbyCbKjMfcb15eXushg1jSHUcLJXgRSZmOw6ePX+8kl5fBEq2lB19ZXs75V1xBt27d6NmzJ5/73OdixrZx40aGDx/OiSeeCMD06dP529/+BkB1dTVnnnkmQMQx8XheffVVrr76arp3705xcXHrUHq06Nu9RjP3e+WVV7Ym+I7+Lp2U4EUkZUrw6aMh+nCCr1q9mleXL2ft2rW88MILrFu3jpqamjbtjxw5wle+8hWeeeYZ1q5dS3NzM7/4xS8A93avP/nJT6ioqKCpqanN30Zbs2YNTz/9NJWVlbzwwgtUVFTEbPfb3/6Wu+++m7fffpsTTjihze/nz5/PxIkT+fvf/864ceOor68HYO3atTQ1NfH3v/+99SS7TLJ667t0Fq9uF6uiohJZvkf4FqezfRBPtpZ843lsAifPcjy2bhcL4dvFTujZ0wGcHj16OBUVFc7pp59u/XXy22ui28WKSMaYPfhjrUWR/aLPoHdsBeIDzbhDzP/+xBMcO2oU3bt358knn2TNmjW2Q8saSvAikjKdZJcefr1EzoYm3AT/H9dey1rgqOV4spGOwYtIynQMPj38dvzdcRwrs8RB5GQ3SlSu/Px8HCfxcR09byKSMiX49PBbD766uprJkydbSfLmpXK256P3g/z8fCZPnkx1dXXCf6MhehFJmRJ8evitB//oo48ya9Yspk6dSl5erAvaMmcg0HLV+k7gsKd79x/HcaiurubRRx9N+G+U4EUkZfuABqAQKAF6AIesRpSd/NaD37dvH3PmzLGy72eBS0L1q0PLkhwN0YtIWqgXnzq/9eBtqjPq3kwLEzxK8CKSFkrwqfNbD94mJfjUKcGLSFoowadOPfgwJfjUKcGLSFpospvUqQcfVmvUS6xFkd10kp2IpIV68KlTDz4saD343sAtQD3uxFBLPdinEryIpIVms0udEnxY0BL8IOChUP0feJPgNUQvImmhHnzqNEQfZg7RByHBm4+hLm6r9FKCF5G0UIJPnXrwYWYSDMIxeDPB13u0T08TfN++fVm6dCn19fVUV1czbdq0dtt37dqVqqoqPvjgA48iFJHO0kl2qVMPPixoQ/TFRt2rHrynx+Afe+wxGhoaGDhwIGPGjGHFihVUVlZSVVUVs/3dd9/Nxx9/TK9eQXh5RYJNPfjU9AS6huoHcWcGzGVBS/CB7sEXFRUxdepUZs+ezYEDBygvL2f58uVMnz49Zvvhw4fz5S9/mR/84AdehSgiKWiZrhbcIdXu7bSVttR7jxS0y+Rs9OA9S/AjR46ksbGRzZs3t66rrKxk9OjRMdvPmzeP++67j0OH2p/ResaMGVRUVFBRUUFpaWlaYxaRxDnALmNZvfjk6Ph7pCD34AOX4IuLi6mtrY1YV1NTE3P4/eqrryY/P58//vGPHW53/vz5jB07lrFjx7J79+60xSsiydMwfeepBx8pyAneqyF6z47B19fXU1ISOdBSUlJCXV3k/zJFRUX88Ic/5PLLL/cqNBFJE51o13nqwUc6inuL2O64iao72X3L2ECfZLdp0yYKCgoYMWIE7777LgBlZWVs2LAhot1JJ53E8OHDWblyJQCFhYX07t2b7du3c84557B161avQhaRJKkH33nqwbdVR/hcjhKyO8EHeoj+4MGDLF26lLlz51JUVMT48eO56qqrWLRoUUS79evXM2TIEMaMGcOYMWO46aab2LlzJ2PGjNHlciI+p9nsOk89+LaCNEwf6LPoAW699VZ69OjBrl27WLJkCTNnzqSqqooJEya0DtU3NTWxc+fO1rJ3716am5vZuXMnzc3NXoYrIklSD77zzASvHrwrSAk+0EP0APv27WPKlClt1r/22mtxr3X/29/+xpAhQzIdmoikgY7Bd56Z4Pdai8JfgnSpXKCH6EUk+NSD77z+Rn2PtSj8JUg9+MAP0YtIsCnBd14/o64evCtICT7QE92ISPDpJLvOU4JvK0gJXkP0IpLV9uFevwzQG01Xmwwl+LaCdAze7MFriF5Eso6mq+08HYNvKyg9+DwiE/wBj/arBC8iaaXj8MnLQ9fBxxKUBB/de3c82q8SvIiklRJ88noT/jKuAZosxuInQRmit3GCHSjBi0ia6US75JnH3zU8HxaUHryNE+xACV5E0kw9+OSZx991gl1YEBO8VyfYgRK8iKSZZrNLns6gjy0oCV5D9CISCOrBJ08JPragHIPXEL2IBIISfPJ0iVxsNUa9t7UoUmfjGnhQgheRNNNJdslTDz428656faxFkTr14EUkENSDT54SfGzRQ/R5tgJJkRK8iATCXsLT1fYBeliMJVvoMrnYmggn+S5k73F4DdGLSCA4wHZjeZCtQLKILpOLz5zVL1uH6dWDF5HA2GbUleA7piH6+IJwHF4JXkQC4yOjfry1KLKHEnx8QUjwGqIXkcBQDz45OgYfn5ng+8Zt5W/qwYtIYJgJXj349ulOcu1TD77zlOBFJO3MIXr14NvXG8gP1XUnubaCkODVgxeRwFAPPnEanm+fEnznKcGLSNqpB584XSLXviAkeA3Ri0hgqAefOJ1B374gJHj14EUkMOoJz0DWncgkJpGU4NsXhASvHryIBIqG6ROjY/Dty/YE34PwSZSH8PYkSiV4EckIDdMnRsfg25ftCd7W8DwowYtIhqgHnxgl+PZl+1z05vC8EryIBIJ68Ikxb6m7M26r3BWkHryXx99BCV5EMkTT1SZmgFHfZS0K/6oFmkP13mRf0jJvcasevIgEgm44kxgzwasH35ZD+IoMcJN8NjFHHbyehlgJXkQyQj34xJhD9OrBx5bNw/Q27zOgBC8iGaEefMfygdJQvRldJhdPUBL8/ritMkMJXkQyYjvhY6cDgQKLsfjVMUZ9N7rRTDxBSfDqwYtIIDQSPqbcBfXiY9EJdolRgu8cJXgRyZitRn2YtSj8S5fIJUYJvnOU4EUkY5Tg26cefGKU4DtHCV5EMkYJvn1K8IlRgu8cJXgRyRgl+PZpiD4xSvCdowQvIhmjBN8+9eATk83z0SvBi0ggKcG3Tz34xKgH3zlK8CKSMWaCHwrk2QrEp9SDT0y2JvjuoQLQgHs/eC8pwYtIxtQR7rV0JzKhiRJ8orI1wdvsvYMSvIhkmIbp49M89Ikxk2Np3Fb+owQvIoGmBB9bCdAtVK8HDlqMxe8+NupK8IlTgheRjFKCj00n2CXuMOF7qReSPbeMVYIXkUBTgo9Nx9+TYz5H2XIuh817wYMSvIhkmBJ8bDr+npxsTPDqwYtIoCnBx2YmKQ3Rd8w8Dq8EnxhPE3zfvn1ZunQp9fX1VFdXM23atJjt7rrrLtatW0dtbS3//Oc/ueuuu7wMU0TSSAk+Ng3RJ8d8jo6xFkVybCf4Ai939thjj9HQ0MDAgQMZM2YMK1asoLKykqqqqoh2eXl5XHfddaxdu5YTTzyRP/3pT3zwwQc8/fTTXoYrImnwMe4Z4kW4J0f1IfK65lx1rFFXD75j2T5Eb+M971kPvqioiKlTpzJ79mwOHDhAeXk5y5cvZ/r06W3aPvzww6xZs4ampiY2bdrEs88+y7nnnutVqCKSZv806iOsReEvQ4z6NmtRZI9sT/CBHqIfOXIkjY2NbN68uXVdZWUlo0eP7vBvzzvvPDZs2JDJ8EQkgzYZ9ZOsReEvZoJ/31oU2UMJPnmeDdEXFxdTW1sbsa6mpoZevXq1+3f3338/Xbp0YeHChTF/P2PGDG6++WYASkuzaQoEkdxhJviR1qLwl6FG/QNrUWQP8yQ7HYNPjGc9+Pr6ekpKSiLWlZSUUFdXF+cv4LbbbuO6665j8uTJNDQ0xGwzf/58xo4dy9ixY9m9e3daYxaR9Nhs1JXgoZjwl/9hIpOXxKYefPI8S/CbNm2ioKCAESPCR+DKysriDr1/5Stf4d5772XSpEls26YjVCLZTEP0kczh+Q8Ax1YgWUQJvnMcr8qSJUucxYsXO0VFRc748eOd/fv3O6NGjWrT7pprrnG2b9/unHLKKUltv6KiwrPHoqKikngZCI4TKvt9EI/tconxfLzsg3iyoRQYz1kTOF18EFN7pdCI92iG9xUv93l6Hfytt95Kjx492LVrF0uWLGHmzJlUVVUxYcKEiKH67373u/Tv35+Kigrq6uqoq6vj8ccf9zJUEUmjnUDLGTi9yZ4eWKaYx991gl1iGoG9oXoXoJ/FWBLhh967p9fB79u3jylTprRZ/9prr0WcbHfCCSd4GZaIeGATcFaofhK5PbmLTrDrnF2EE/sAwM9nXfkhwWuqWhHxhE60C9Mlcp2TTcfhbU9yA0rwIuIRnWgXph5852RTgu9v1PfGbZVZSvAi4gn14MPUg++cbLoWfpBR/8hSDErwIuIJTXYTFn2ZnCQmm3rwSvAikjPMHvwIIM9WIJYdA3QP1fcB9RZjyTbZlOCPM+pK8CISaPsJD7H2IHdvHavj752XTQlePXgRySnrjfpp1qKwS8ffO888Bp9NCX67pRiU4EXEM5VGvcxaFHapB995Zg9eJ9l1TAleRDyjBK8efCp2GPXjrUXRsXwiRxh2xGuYYUrwIuIZJXg40agrwSdnH9AyqXkxkdea+8lAwsl1J+40uzYowYuIZ6oIf9mNwP2SzjWfNOrvWIsie20x6sNtBdEBPwzPgxK8iHjoCJFJ7VRbgVhSQGQPXgk+edVGfbilGDqiBC8iOSmXh+lHAF1D9a3AQYuxZKtqo/4JW0F0wA9n0IMSvIh4LJcTvIbnU6ch+sQpwYuIp5TgXf+wFkV2qzbqwy3F0BE/zGIHSvAi4jEzwZ9Kbk1Ze4pRV4LvnGqjPtxSDB1RD15EctJOwtcFF5NbN55RDz511UZ9uKUYOqIELyI5a5VRP9daFN7KI7IHr2PwnbM/VAB64s8Z7XSSnYjkrNeM+nnWovDWYMLX/e8hcl51SU61UR9uKYZ4uhKexa4Zd8TKFiV4EfGcmeAnWIvCWxqeT59qoz7cUgzxDDTqO4EmW4GgBC8iFrwFHArVRwDHWozFKzrBLn38fKmcX46/gxK8iFjQQORx+Fzoxesa+PSpNup+m+zGvFugEryI5KSVRj0XjsOfZdTXW4siGKqN+nBLMcRj/iO3yVoULiV4EbEil47DFwFjQvVm4E2LsQRBtVEfbimGeEYb9Q3WonApwYuIFf9H+ASkMqCXxVgybSzujWbA/dKvsRhLEJjH4E8gPL+/H4wy6lXWonApwYuIFXWEZ7XLBy6yGEumjTfq/2ctiuCoI5zkuxGZVG0qAE42lpXgRSRnPW/Ur7IWReYpwaffW0b9TGtRRDoRKAzVP8D9R8QmJXgRseZZo/453J580OQB44xlJfj0MBP8GdaiiOSn4++gBC8iFr0FbAvV+xPMk+1G4j42cGeve9diLEHytlH3Sw/eT8ffQQleRCxygOXGchCH6TU8nxlmgi/DH6M/6sGLiBj+aNSDmODNUQkl+PTZDbwfqvcg8vpzW9SDFxEx/BWoDdVPAE63F0radQEmG8sr4zWUTvHTiXb5+OsMelCCFxHLGogcpr/ZViAZMIHwzUe2owlu0s0cprd9ot2JuJfsAXxI+J9Wm5TgRcS6J4z6lwnOpDf/YtSX4s5iJ+njpx58mVH3w/F3UIIXER9YSXh+9mLcJJ/t8oCpxvLvbQUSYGYP/nSgu61AgIuNeoW1KCIpwYuIL/zCqM+0FkX6jCN869Bd6Ph7JuwENobqRUQmWa9dZtRfsBZFJCV4EfGFRcCBUP1U7H5Zp8O/GvVlhOfdl/RaatSnWIphFDAkVN+Hf861UIIXEV+oBZ40lh/CHebORr2BG4xlDc9njpngryR8Ux8vfdao/xn//DOnBC8ivvFd4GCofjpwrcVYUnEbbpIH+AfwssVYgm417rzv4M4Y+BkLMZjD8y9a2H88SvAi4hvbgUeM5e9h98SpzigCZhnLD+LO2CeZs8yoez1M3xM4z1hWghcRieNh3JPSAIYC/2Uxls6YARwTqlcDi+2FkjPMYfqpeHt/+MsIX/9eiftPql8owYuIr9QB9xrLM4m83MzPhgAPGMs/BBotxZJLXsM9ox7gOOBGD/d9j1FfHreVHUrwIuI7C4HfGcu/JHKebz/Kwz1JsOXY+7u4j0Myrwn4kbH8H3hzaOcSYGyofhh4zIN9JkMJXkR86WZgS6jeB/gL/rihSDz3AheE6k3Adbhf+uKNxwgPjx8P3OrBPmcb9fmERxH8QgleRHypBneq15Y5vQfiJvnxcf/CnlnA943lHwCvW4olVx3CvQqjxX/gzg+fKVcSvlNgA+7hGL9RghcR33ob9ySmutDyscDfcI972rjeOVo+7jH3HxvrVgJz7YST8+YD/wzV+wL/izv6k26jcCdmavFr3BvM+I0SvIj42uu4E4nsCS0X4F56tha3F2VrMpzTgHLgP411K4HLgaNWIpKjwJdwe/MApwDPEb6jXzoMBZ4FSkLLW3FHC/xICV5EfK8cGAP8n7Huk7hftBuBO4HhHsSRD1yEe911JXC28bu/4Cb3eg/ikPgqgOuN5Qm4/wx+ntT/Gfwi7us+IrRcj/tP5scpbjeTnKCUiooK6zGoqKhkrhSA8y1wasBxYpT14DwOzg3gjAOnNMX9lYJzDjgzwHkKnB0x9nkEnHvA6eKD50clXGaB0xT1WlWBczs4pySxnRPAmQnOuqhtHQXnKh88Toif+/JClUCoqKhg7NixHTcUkax2DO5x+Bvp+Bjrftyz8ffi3ghkH+4x/WbC34T5uPegL8G9zO0Y3F5aR9teCswhfKtb8ZcLgN/gXhsfbS9QBbwP7MYd1m8GeuC+7sOBkYTvCGj6J3AN/rmpTLzc5+l5Kn379mXBggVccskl7N69m29/+9ssWbIkZtsHH3yQm266CYBf/vKX3HvvvTHbZdypwCTcT33LJNlFUfUa3Mmm12U4jstC+yO0/5Y5EVviay+OeI8j0b+J166jNubvm3EPCnnxfEny4r2Wibx3zDYtmnEnCo++d2b09gqAwtDvHNxx1JZtEyMm4ONJcFdv+M/D8MV/wOf/ARf/E7rFuMtHH9x57dPlo57wu0/BL8+ADQOMmE0NwN+Bk6NiH4x74XR0+5bPc65+JmJ9v60n/PwdxP0vrJvx+1jPV9R2XjkIp70F36qDWyqhpCHctB/hs+ATVVcIC06GOZ+F2h6hlQ24Z/O199p9FjgL9/sv3uciAzztwS9evJguXbpw4403MmbMGFasWMH48eOpqqqKaHfzzTdzxx13MGnSJBzH4c9//jM//elP+e///u92t5/2HvypwBWEv3za04B7NkcmPqCnAlfR9t+xJtxXz1wfK46OHkeifxPdrqM27e03k8+XJC/ea7kGNzu2995pr40DrCL8ZZbMZ6plCriCDtaF9DwCZ2+Dc9+HT+2Ck/bCSXugOIUz3uq7wrv9YHN/eKM7vDoG3h4MzYmcvRSd+Jtwv+DjHQhuxD2pINc+E/G+32L942SKfr462E7JYfjierjsXbhwC/Q5klh4dYXw6jB4cQQ8VQa13WLE1QT8kdiv3WeBT0f9TfTnIkXWe/BFRUVMnTqVT33qUxw4cIDy8nKWL1/O9OnT+fa3vx3R9vrrr+eRRx5h27ZtADzyyCPMmDGjwwSfdpNI7IuIULtJZObDOYnYr1R+gnF09DgS/Zvodh21aW+/mXy+JHnxXsuziP0+S7RNXuj3LV9kyXymYr3n2/nGOtAN/nKCW1o5cGw9HF8LfQ9D30Puz+IGyHOgi+OG2JznfnHXdoPaQtjfA/7ZF3YUE/5i7ijhRItu297z2PLYcvEzEe/7raPnOvr56mA7td1h/lluwYHBtXDKbvf9UXoQCpsgvxkOdYWabrC9F2zuB9V9oKmj1y6f+K/dWTEeS/TnIkM8S/AjR46ksbGRzZs3t66rrKxk4sSJbdqOHj2aysrKiHajR4+Oud0ZM2Zw8803A1BaWpreoHt33CSl9pnabnT7RP4+0b/pnUSbjvabqedLkhfvtUikp9pRG/P3Xr/mebCjl1sSlmwiT6dc/Eyk8pgT+T6KJQ8+7O2WtEn2M+TBNWyeXSZXXFxMbW1txLqamhp69Wr7ySsuLqampqbDdgDz589n7NixjB07lt27d6c36JqOm6TUPlPbjW6fyN8n+jc1SbTpaL+Zer4kefFei+YE/rajNubv9Zq3Lxefn1QecyLfR15J9jOUyGcrRZ4l+Pr6ekpKSiLWlZSUUFdX12HbeO0y7mXcY4yJaCB8UlAm4oh1S6qmGOtjxdHR40j0b6LbddSmvf1m8vmS5MV7LVfHWJ9MGyf0+/b2E08jbd/fsdYlItEzjRqB92K0d0Lrk9l3uExZNAAAB/xJREFU9DaaYqyL3ncufibifb919JpFP1/JbCfW6xvjBM2YYm2vifiv3eo4+1sdo22aeZbgN23aREFBASNGjGhdV1ZWxoYNG9q03bBhA2VlZR22y7h1uCeC7cd9QQ6ESnR9P5k9YWwd7skkLftr2f8fQ+tb4osXR3uPI9G/idWuozbRv29qZ1tiV7zX8gU6fu9Et2kpTbQ9kSjWe/GI8TfNxrafpe37O3pdU9RP89o3s+xMoM2B0Lb/JxS3ue1VofXRn8N42zsS+hsz9j+G1sVq37LvXPxMxPt+M5+/A7h37mnv+WpvO9Hr3qPt6/tHYr/Po1/rnUS+hkeIf4IduO//WO8nD86iB48uxAecJUuWOIsXL3aKioqc8ePHO/v373dGjRrVpt0tt9ziVFVVOYMGDXKOO+44Z/369c4tt9zS4fY10Y2KioqKSq6VdnKfd0H07dvXWbZsmVNfX+9s3brVmTZtmgM4EyZMcOrq6iLaPvTQQ86ePXucPXv2OA899FCqD1JFRUVFRSWQRTPZiYiIBFC83KebzYiIiASQEryIiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgABWqim127drF169a0ba+0tDT9d6jLQXoeU6fnMHV6DlOn5zB1mXgOhw0bxoABA2L+zvo0e34tmvpWz6Nfip5DPYd+KHoOs+s51BC9iIhIACnBi4iIBFA+cL/tIPzs7bffth1CIOh5TJ2ew9TpOUydnsPUefUcBuokOxEREXFpiF5ERCSAlOBFREQCSAk+hr59+7J06VLq6+uprq5m2rRptkPKOrfddhsVFRUcPnyYhQsX2g4nKxUWFvLLX/6S6upqamtrWbNmDZdddpntsLLOokWL+Oijj6ipqWHjxo3ceOONtkPKWiNGjODQoUMsWrTIdihZ55VXXuHQoUPU1dVRV1fHO++848l+rV8X6LeyePFi57e//a3Ts2dP59xzz3X279/vjBo1ynpc2VSmTJniXHXVVc7Pf/5zZ+HChdbjycZSVFTkzJkzxxk2bJiTl5fnTJ482amtrXWGDRtmPbZsKqNGjXIKCwsdwDn55JOd7du3O2eccYb1uLKxvPTSS86rr77qLFq0yHos2VZeeeUV58Ybb/R0n+rBRykqKmLq1KnMnj2bAwcOUF5ezvLly5k+fbrt0LLKsmXLePbZZ9mzZ4/tULLWwYMHeeCBB9i6dSuO47BixQq2bNnCmWeeaTu0rFJVVUVDQwMAjuPgOA4nnnii5aiyzxe/+EX279/Pyy+/bDsUSZASfJSRI0fS2NjI5s2bW9dVVlYyevRoi1GJwIABAxg5ciQbNmywHUrWeeyxxzhw4AAbN25k+/btPP/887ZDyiq9evVi7ty53HHHHbZDyWo/+MEP+Pjjj3nttdeYOHFixvenBB+luLiY2traiHU1NTX06tXLUkQiUFBQwG9+8xuefPJJNm7caDucrHPbbbfRq1cvJkyYwNKlSzly5IjtkLLKd77zHRYsWMC2bdtsh5K17rnnHk444QSOP/54nnjiCZ577jlOOOGEjO5TCT5KfX09JSUlEetKSkqoq6uzFJHkury8PBYtWkRDQwNf+9rXbIeTtZqbmykvL2fw4MHMnDnTdjhZo6ysjIsuuogf//jHtkPJaqtWraK+vp6GhgaeeuopysvLufzyyzO6z4KMbj0Lbdq0iYKCAkaMGMG7774LuG9wDYuKLQsWLGDgwIFcfvnlNDY22g4n6xUUFOgYfBLOP/98hg8fzvvvvw+4o5z5+fmMGjVK54OkwHEc8vLyMr8flciyZMkSZ/HixU5RUZEzfvx4nUXfiZKfn+9069bN+f73v+889dRTTrdu3Zz8/HzrcWVbefzxx53XX3/d6dmzp/VYsrEcc8wxzhe/+EWnZ8+eTpcuXZxLLrnEqa+vd6644grrsWVL6dGjhzNw4MDW8vDDDzvPPPOMU1paaj22bCm9e/d2LrnkktbvwWuuucapr693TjrppEzv2/6D91vp27evs2zZMqe+vt7ZunWrM23aNOsxZVuZM2eOE23OnDnW48qmMnToUMdxHOfQoUNOXV1da7nmmmusx5YtpbS01PnrX//q7Nu3z6mpqXHWrl3r3HTTTdbjyuYyZ84cXSaXZCktLXVWrVrl1NbWOvv27XNef/1156KLLsr4fjUXvYiISADpJDsREZEAUoIXEREJICV4ERGRAFKCFxERCSAleBERkQBSghcREQkgJXgREZEAUoIXEREJICV4ERGRAFKCF5FOmTp1KocPH2bo0KGt6x599FHeffddBgwYYDEyEWlhfZ5eFRWV7CwVFRXOE0884QDOnXfe6ezcudMZMWKE9bhUVFRw8EEAKioqWVouvvhip6Ghwbnnnnucmpoa56yzzrIek4qKilt0sxkRSUl5eTmf/vSnueKKK3jxxRdthyMiIToGLyKddsEFF1BWVkZeXh47d+60HY6IGNSDF5FOOe2003j11Ve5/fbbmTx5MsXFxVx22WW2wxIRg/XjBCoqKtlVhg4d6mzbts2ZPXu2AzijR492mpqanIkTJ1qPTUVFxS3qwYtIUvr27Ut5eTmvvvoqX/3qV1vX//a3v2Xo0KGMHz/eYnQi0kIJXkREJIB0kp2IiEgAKcGLiIgEkBK8iIhIACnBi4iIBJASvIiISAApwYuIiASQEryIiEgAKcGLiIgEkBK8iIhIAP1/2cBybbBTMO4AAAAASUVORK5CYII=\n", "text/plain": [ "