TD_enonce.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b874d320-25b0-44ec-b0e5-64b8cdae4de8",
   "metadata": {},
   "source": [
    "# Transfer Learning Tutorial : Unsupervised Instance-Based Approaches"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc01bdce-baf0-409f-a1a5-10693ff04e0c",
   "metadata": {},
   "source": [
    "In this tutorial, we are going to apply instance-based approaches to perform transfer learning.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a51320a9-806b-468e-88bf-d940c2fdb0e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a73a1f5-e80a-4109-ae69-09c9f13b37c8",
   "metadata": {},
   "source": [
    "# Exercise 1 : Gaussians 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d8845cfc-97bb-4e0f-bf54-4f35fe7d0c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "def show(Xs=None, Xt=None, ys=None, yt=None, model=None, weights=None):\n",
    "    \"\"\"\n",
    "    This is the plotting function \n",
    "    \"\"\"\n",
    "\n",
    "    if ys is None and yt is None:\n",
    "        fig, ax = plt.subplots(1, 1, figsize=(8, 3))\n",
    "\n",
    "        alpha = 0.1\n",
    "        for x in Xs:\n",
    "            ax.plot([x, x], [0, alpha], c=\"C0\")\n",
    "        ax.plot([x, x], [0, alpha], c=\"C0\", label=r\"$x_i$ observations\")\n",
    "        \n",
    "        alpha = 0.1\n",
    "        for x in Xt:\n",
    "            ax.plot([x, x], [0, alpha], c=\"C1\")\n",
    "        ax.plot([x, x], [0, alpha], c=\"C1\", label=r\"$x'_i$ observations\")\n",
    "\n",
    "        ax.tick_params(left = False, right = False , labelleft = True ,\n",
    "                        labelbottom = True, bottom = False, top = False)\n",
    "        ax.spines.top.set_visible(False)\n",
    "        ax.spines.right.set_visible(False)\n",
    "\n",
    "        ylim_max = 0.45\n",
    "\n",
    "        ax.plot(1, 0.0, ls=\"\", marker=\">\", ms=10, color=\"k\", clip_on=False,\n",
    "                transform=ax.get_yaxis_transform())\n",
    "        ax.plot(-4.3, ylim_max, ls=\"\", marker=\"^\", ms=10, color=\"k\", clip_on=False)\n",
    "\n",
    "        # ax.set_ylim(0., ylim_max)\n",
    "        ax.set_xlim(-4.3, 4.3)\n",
    "\n",
    "        ax.set_xlabel(\"X ~ P(X)\", fontsize=16)\n",
    "        ax.set_ylabel(r\"$\\widehat{p_s}(x)$\", fontsize=16)\n",
    "        \n",
    "        if weights is None:\n",
    "            sns.kdeplot(Xs, label=r\"$\\widehat{p_s}(x)$\", ax=ax, shade=True)\n",
    "        else:\n",
    "            np.random.seed(123)\n",
    "            bs_index = np.random.choice(len(Xs), 3 * len(Xs), p=weights/weights.sum())\n",
    "            sns.kdeplot(Xs[bs_index], label=r\"$w(x) \\widehat{p_s}(x)$\", ax=ax, shade=True)\n",
    "        sns.kdeplot(Xt, label=r\"$\\widehat{p_t}(x)$\", ax=ax, shade=True)\n",
    "\n",
    "        ax.legend(fontsize=16, loc=\"upper right\", bbox_to_anchor=(1.4, 1.05))\n",
    "    \n",
    "    else:\n",
    "        fig = plt.figure(figsize=(8, 8))\n",
    "\n",
    "        gs = fig.add_gridspec(2, 1, height_ratios=(2, 1),\n",
    "                              left=0.1, right=0.9, bottom=0.1, top=0.9,\n",
    "                              wspace=0.05, hspace=0.0)\n",
    "        # Create the Axes.\n",
    "        ax2 = fig.add_subplot(gs[1, 0])\n",
    "        ax1 = fig.add_subplot(gs[0, 0], sharex=ax2)\n",
    "        \n",
    "        if weights is None:\n",
    "            lns1 = ax1.plot(Xs, ys, '.', ms=15, alpha=0.7, c=\"C0\",\n",
    "                            markeredgecolor=\"C0\", label=r\"$(x_i, y_i)$ observations\")\n",
    "        else:\n",
    "            lns1 = ax1.scatter(Xs, ys, s=100*weights, alpha=0.7, c=\"C0\",\n",
    "                            edgecolor=\"C0\", label=r\"$(x_i, y_i)$ observations\")\n",
    "            lns1 = [lns1]\n",
    "        lns8 = ax1.plot(Xt, yt, '.', ms=15, alpha=0.7, c=\"C1\",\n",
    "                        markeredgecolor=\"C1\", label=r\"$(x'_i, y'_i)$ !not available!\")\n",
    "\n",
    "        lin = np.linspace(-4.2, 4.2, 100)\n",
    "\n",
    "        ax2.plot(lin, gaussian(lin, mu=-1), color=\"C0\")\n",
    "        lns2 = ax2.fill_between(lin, gaussian(lin, mu=-1),\n",
    "                                alpha=0.2, color=\"C0\", label=r\"$p_s(x)$\")\n",
    "        ax2.plot(lin, gaussian(lin, mu=1.), color=\"C1\")\n",
    "        lns6 = ax2.fill_between(lin, gaussian(lin, mu=1.),\n",
    "                                alpha=0.2, color=\"C1\", label=r\"$p_t(x)$\")\n",
    "\n",
    "        alpha = 0.1\n",
    "        for x in Xs:\n",
    "            ax2.plot([x, x], [0, alpha], c=\"C0\")\n",
    "        lns3 = ax2.plot([x, x], [0, alpha], c=\"C0\", label=r\"$x_i$ observations\")\n",
    "\n",
    "        alpha = 0.1\n",
    "        for x in Xt:\n",
    "            ax2.plot([x, x], [0, alpha], c=\"C1\")\n",
    "        lns7 = ax2.plot([x, x], [0, alpha], c=\"C1\", label=r\"$x'_i$ observations\")\n",
    "\n",
    "        ax2.tick_params(left = False, right = False , labelleft = False ,\n",
    "                        labelbottom = True, bottom = False, top = False)\n",
    "        ax2.spines.top.set_visible(False)\n",
    "        ax2.spines.right.set_visible(False)\n",
    "\n",
    "        ax1.tick_params(left = True, right = False , labelleft = True ,\n",
    "                        labelbottom = False, bottom = False, top = False)\n",
    "        ax1.spines.top.set_visible(False)\n",
    "        ax1.spines.right.set_visible(False)\n",
    "        ax1.spines.bottom.set_visible(False)\n",
    "\n",
    "        ylim_max_ax2 = 0.45\n",
    "        ylim_min_ax2 = 0.\n",
    "        ylim_max_ax1 = ax1.get_ylim()[1]\n",
    "        ylim_min_ax1 = ax1.get_ylim()[0]\n",
    "\n",
    "        if model is not None:   \n",
    "            if hasattr(model, \"predict\"):\n",
    "                yp = model.predict(lin.reshape(-1, 1))\n",
    "            else:\n",
    "                yp = model(lin.reshape(-1, 1))\n",
    "            lns4 = ax1.plot(lin, yp, c=\"k\", label=r\"$\\widehat{h^*}$  predictive model\")\n",
    "\n",
    "        ax2.plot(1, 0.0, ls=\"\", marker=\">\", ms=10, color=\"k\", clip_on=False,\n",
    "                transform=ax2.get_yaxis_transform())\n",
    "        ax1.plot(-4.3, ylim_max_ax1, ls=\"\", marker=\"^\", ms=10, color=\"k\", clip_on=False)\n",
    "\n",
    "        ax2.set_ylim(ylim_min_ax2, ylim_max_ax2)\n",
    "        ax1.set_ylim(ylim_min_ax1, ylim_max_ax1)\n",
    "        ax2.set_xlim(-4.3, 4.3)\n",
    "        ax2.set_xlim(-4.3, 4.3)\n",
    "\n",
    "        ax2.set_xlabel(\"X ~ P(X)\", fontsize=16)\n",
    "        ax1.set_ylabel(\"Y ~ P(Y|X)\", fontsize=16)\n",
    "\n",
    "        lns = lns1 + [lns2] + lns3 + lns8 + [lns6] + lns7\n",
    "        labs = [l.get_label() for l in lns]\n",
    "        ax1.legend(lns, labs, fontsize=16, loc=\"upper right\", bbox_to_anchor=(1.4, 1.05))\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76cc23ad-bfe9-4642-91dc-8cdbf6c85a5c",
   "metadata": {},
   "source": [
    "## 1. Weighted Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a306dde1-1e31-4d2a-a182-939560c81994",
   "metadata": {},
   "source": [
    "In this exercise, we consider two domains defined on $\\mathbb{R}$ : the source domain $P_S(X, Y)$ and the target domain $P_T(X, Y)$ where:\n",
    "\n",
    "$$P_S(X) = \\mathcal{N}(-1, 1)$$\n",
    "\n",
    "$$P_T(X) = \\mathcal{N}(1, 1)$$\n",
    "\n",
    "$$P_S(Y|X=x) = P_T(Y|X=x) = \\mathcal{N}(|x|, 0.5)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "571555c9-6791-49ff-9b53-da0fa56ba28f",
   "metadata": {},
   "source": [
    "We remind the definition of the gaussian density function of mean $\\mu$ and standard deviation $\\sigma$:\n",
    "\n",
    "$$p(x) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} \\exp \\left( \\frac{-(x-\\mu)^2}{2 \\sigma^2} \\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4877a8cf-527a-4004-8d89-372970e362a9",
   "metadata": {},
   "source": [
    "We consider the scenario where the learner has access to a source **labeled** sample `(Xs, ys)` of size $m=30$ drawn according to $P_S(X, Y)$ and a target **unlabeled** sample `Xt` of size $n=30$ drawn according to $P_T(X)$.\n",
    "\n",
    "**NB** : For the tutorial purpose, we sample also the associated labels `yt` which are not available at training time but will be used to evaluate the models on the target domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "50cb163b-6684-4915-86b7-f39565f853ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "def label_func(x, noise=0.1):\n",
    "    return np.abs(x) + noise * np.random.randn(len(x))\n",
    "\n",
    "def gaussian(x, mu=0., s=1.):\n",
    "    return 1./np.sqrt( 2. * np.pi * s**2 ) * np.exp( -(x-mu)**2 / ( 2. * s**2 ) )\n",
    "\n",
    "np.random.seed(0)\n",
    "Xs = np.random.randn(30) * 1. - 1.\n",
    "Xt = np.random.randn(30) * 1. + 1.\n",
    "ys = label_func(Xs, 0.5)\n",
    "yt = label_func(Xt, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6ca37293-9c60-4b93-a7a0-5bf289eb7a8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqwAAAIGCAYAAACRRd24AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACfOklEQVR4nOzdeVxU9foH8M+ZBYZhX2UXREAUUdHU3NPcbmZupamVmpVlLnlNzW5lprb9tLyWlZlaauVycatr3rTFxF3ccBcXUET2TRgYZs7vjxFj35lzgM/79eJFc9aHU+nD9zzf5yuIoggiIiIiIrlSSB0AEREREVFFmLASERERkawxYSUiIiIiWWPCSkRERESyxoSViIiIiGSNCSsRERERyZpK6gCIiIjMJSoqaqBKpXpHFEV3cNCGSA6MgiCki6J4uqCgYEnHjh2vlHWQwD6sRETUFERFRQ20tLT8zM/PL9/KykqnUCj4FyCRxERRhF6vV2VmZlonJCSIeXl50zt27PhLyeNk/9ulIAgegiD8KQiCu9SxEBFRw6VSqd7x8/PLt7a2zmWySiQPgiDAwsKiwMXFJcPPz69ArVbPKes42SesAN4C0OP+dyIiohoRRdHdyspKJ3UcRFQ2a2vrHFEU/cvaJ+uEVRAEDwATYYpzIkdZiYioFhQcWSWSL0EQAEAoa5+sE1aYRlULA1eAo6xERERETY5sE9Yio6uW9zdZgqOsRERERE2OnNtaFR1dBQBYWFhYTZ48+Y5E8RARUQNmYWEBAM2ljqOenZA6AKL6IMsR1jJGVwEA+fn5WLt2LRISEqQJjIiIiIjMTpYJK8oYXS1kMBjw3nvvmTkcIiIi+Zo+fTqGDBkidRilLFy40C0oKKi1wWCos2vOmjXLUxCEjnq9vs6uKQfr1693WLBgQbOy9hX+zOaOSU5kl7CWN7paiKOsREQkBzFJ2Vj5+1W8EXEGK3+/ipikbGniiInBl19+iQULFkhy/4rMmjUrKTU1Vf3ZZ585Sx2L3G3fvt3hiy++KDNhnTp1atLevXsvmjsmOZFdwooKRlcLcZSViIiktC3qFl778RR+iU5ATGI29pxLwGs/nsK2qFtmj+XTTz9Fu3bt0KlTJ7PfuzI2NjbiqFGjUlasWNGoJkzn5uZWmKfUtYCAAH2/fv3umfOeciPHhPVhlDO6Wig/Px8HDx40UzhERER/i0nKxtrIG3DQquFqawlbjRouNpZw0Kqx9uANs4605uXlYcOGDRg7dmyx7dHR0ZYqlSp85syZnkW3jxs3ztfa2rrD/v37tdW9V/fu3QPbtWvXquT2o0ePWqlUqvAvvvjCqazzxo8fnxoTE6P59ddfrSu7x9atW+3at2/fSqPRhNva2rZ/9NFHA06fPl1mTnDy5EmrLl26BFlZWXVwdXUNmzlzpmfR0oMzZ85Y9u/fP8DJyamdpaVluIeHR9vBgwe3KFlKcOjQIau+ffu2tLOza6/RaMLDw8Nb/fLLLzaF+wtfxx87dkzTo0ePQK1W22HIkCEt1qxZ4ygIQscjR45YlYytd+/eLYODg1sXfo6OjrYcNmyYv5eXV1uNRhPu7e3ddty4cb5JSUlKABg5cqRfRESEc2JioloQhI6CIHT08vJqWzKG6j6rwvPOnj1r2adPn5ZarbaDp6dn29mzZ3vU5FlJSXYJqyiKHURRFIp+dezYEaIoFvs6efKk1KESEVETtCfaVJKmVhb/K1StVEAost8cDh8+jPT0dPTs2bPY9tDQ0LzRo0cnr169ullCQoISAGbPnu2xZcsWl40bN17t1atXTnXv1bVr1+wLFy5oi44uGo1GvPzyy74dOnS49/LLL6eWdd7DDz+cY21tbfz555/tK7r+1q1b7UaPHh2o1WoN33zzTczHH38ce/nyZas+ffq0un79urrk8SNHjgzo06dP5saNG2OGDRuWunz5co/XX3/9QYL++OOPByYkJFgsW7bsZkRExOV33nnntoWFhWgwGB7Ef+DAAW3fvn1bZWRkKJcvX37z22+/jXFwcCgYOnRo0F9//VUsqR8xYkTLHj16ZP34449XX3vttcQxY8ak29jYGNauXVus3CEuLk4VGRlpP2bMmJQi29Te3t75H374Ydy2bdsuz5kzJz4yMtKuf//+gQCwcOHCO717985wdHQs2Lt378W9e/de3Lx589W6elYjRoxo2bt376zvv//+6sCBA9OXLl3qWbRMoyrPSmpybmtFREQkO3FpOdCoyx7vsVQpcCu92rlgjR0+fBiCICAsLKzUviVLltyJiIhwfueddzyCg4N1n3zyiedXX311bdiwYVk1uVfPnj2zly1bJhw8eFBb+Hr6888/dz59+rT1oUOHzpd3nlKpRHBwcM6xY8cqHGFdsGCBl7e3d96ff/55Ra025Vx9+vTJDg0NDV28eHGz1atXF6u3eOaZZ5KXLFmSAAAjRozIzMrKUn711VfN5s+ff1ev1wuxsbGWGzZsuDpu3LiMwnOmTJlSLKl+/fXXvT08PPIjIyMvazQaEQBGjhyZERQU1Obdd9/12Lt3b0zhsS+++GLiW2+9lVj0/Mceeyxt27ZtTp9//vktpVIJAFi7dq0TAEyaNOnBvQYPHpw9ePDgB0Pvjz76aHZwcHDeoEGDgiMjI626d++e6+zsXKBWq8WqvPqv7rN69dVXE2bMmJECAMOGDcuKjIy03bx5s9OMGTNS7ty5o6rKs5Ka7EZYiYiI5MzHUYu8AmOZ+/IKjPB2qPbb9hqLj4+HnZ1dYY/ZYpo3b66fPHly4rp169zmzJnTfNGiRbGTJ09OK9yfmpqq6NSpU3BVZ/A/8sgj95RKJSIjI60BIDk5WblgwQLv5557Lumhhx7SVXRNZ2fngsTExNJB3peZmak4f/689oknnkgtTMAAoFWrVvnh4eH3Dh06ZFvynPHjxxdLqJ5++unUnJwcxYkTJ6yaNWtW4O3tnff22297L1261OXs2bOlygqys7OFY8eO2Q4dOjRNqVSKer0eer0eoiiiZ8+emceOHSt2zzFjxqSXvMaECRNSEhMT1bt27Xpw7KZNm5y7du2a2bx58wfv03U6nTBv3jx3f3//NhqNJtzCwqLjoEGDggHg3LlzmvKeS109q1GjRmUU/RwcHJwbHx9vAQBVeVZywISViIioGgaGukMUAb2heNKqNxgh3t9vLjqdDpaW5ecXgYGBuvz8fKFdu3bZb7zxRlLRfU5OTsbjx49fKhwZrIy9vb0xODg45+jRozYA8Nprr3kpFArx448/vl3ZNTUajVGn05X7ejkpKUkpiiI8PDxKFU26ubnp09PTS70R9vb2Lij62dPTUw8AsbGxaoVCgT179lwOCwu7t2jRIq+wsLBQb2/vth9++KFrkXuqDAYDli9f7mFhYdGx6Nd3333nlpmZqSyaePv6+paKbcCAAdmenp753333nTMAREVFac6fP68dN25cStHjpk2b5rV06VLPJ598MnXz5s1X/vjjjwvffvttDADodLpq5WI1eVZubm7FnpWFhYWYn5+vAICqPCs5YEkAERHJX/IV4MJOIC0WcPQFQoYCLoGShBLgaoOJ3f2w9uANCDCVAeQVmJLVid38EOBqU9kl6oyzszPS09PL3Ldjxw7b2bNnN2/fvv29kydP2hw5csSqS5cuuYX7Cydkffrpp/FVvd9DDz2U/euvvzocOHBA+/3337uuWLHiupOT04PMvbxrpqenqxwdHQtKXq+Qq6urQRAEJCQklKq/TExMVDs4OJQ699atW6rWrVvnF36Oj49XA38nlq1bt87ftm3bDaPRiMOHD1stX77cbd68eb7+/v55Tz31VKazs7NBoVDgmWeeSZw0aVJKyesDpnKGQgqFQiy5X6FQYNSoUSmrV69ulpWVFbtmzRpnrVZrHD9+fHrR43bs2OE0YsSIlI8++ujBap27du2q2m8KJdTkWVWmsmdVkzjrGkdYiYhI3k5vAiJeAC78BCRfBi78bPp8epNkIQ0P98Yno9tjYBt3tGxmg4Ft3PHJ6PYYHu5t1jhatWqF/Px83LpVvJ3WgQMHtOPGjWs5evTo5EOHDl3y8PDInzt3rlfRY06dOqV96KGHqtUqqUePHtnx8fEWEydO9O/QoUP2K6+8Uuy1fHnXjIuLswgICNCVd107Oztj69atc3bu3OlYUPB3vnX58mWLkydPWnfr1q1U3e2GDRuKdSX44YcfnLRarbFTp065RbcrFAp069Yt94svvogDgLNnz1oV3rNjx47Z586d03bv3j2nV69epb6q8kyef/75lJycHMX69esdIiIinAYOHJhma2tbbPhdp9MpVCpVsYR3zZo1xSZrWVpainl5eZXmZTV5VlVV3rOSA46wUqVikrKxJzoBcWk58HHUYmCou1lHEIioCUu+Ahz5ArByApRFBpQMeuDIl4BXuKQjra880lKSexfq1asXAODo0aPw9jYly9HR0ZZPPPFEYI8ePTLWrVsXq1QqMWfOnPiZM2f67d6926Zw8k90dLR1t27dbhZe69KlSxatWrVq+9prr91ZtmxZmaOu/fr1ywaA69evayIjI0tNtCp5TcBU63rz5k3NtGnT7lb0syxYsOD26NGjA/v27Rs4ZcqUxKysLOX777/vaWNjY3jzzTdLnbt+/XoXo9GILl265Ozevdtu06ZNLrNmzYp3dnY2HDlyxGr69Ok+I0eOTAsKCtIZDAZh3bp1zkqlUhwwYMCDhG7ZsmVxAwYMCO7Zs2fghAkTkr28vPRJSUmqEydOWBsMBqxcufJ2yfuWFBYWlhcWFnbv3Xff9U5MTFQ/99xzpUZre/funREREeH8wQcf5AYFBeVt3brV4cSJE8X+Ig0JCcn94YcfXD788EPXrl273rOyshI7d+6cW/JaNXlWFanqs5IaR1ipQnJqjk1ETdCFnQCE4skqcP+zcH9/0+Xn54fOnTtj165dAIDY2FjVoEGDggICAnTbtm27XvhKe+rUqSn+/v66+fPnewNATEyMWhAE0d/f/0EdZFZWlgIA3N3dy22+aWdnZ1Sr1eKzzz6bWLS8oLxrAsCWLVvs1Wq1OHbs2DRUYNSoUZmbNm26kpmZqZw0aVLA7NmzfQMCAnL/+OOPi35+fqViioiIuPr777/bPf300y0jIiKcp0+ffqfwlbu3t7fey8srf+XKlc2efvrplpMmTWqRkJBgsWnTpqs9e/Z8MHLao0ePnAMHDlxwcnIyzJs3z3f48OFBc+fO9T137pxV7969q9xQd8yYMSmJiYlqNzc3/ZAhQ0oleatWrYrr169fxpIlS7wmTJjQIjs7W7l+/fprRY+ZMWNG8pAhQ1IXL17s1adPn5Dhw4eX+9tQdZ9VRar6rKQmiGKpkgzZ6dSpk3j8+HGpw2hyYpKy8dqPp+CgVRfrN6g3GJGeq8cno9tzpJWI6tfOGaYyAI1d6X26TMA1GHj80ypd6sKFCwgJCanb+GRg3bp1mDFjBu7cuQOtVnuiKuds2LDB4ZtvvnH5/fffH/T6/L//+z+XxYsXe8XGxp4t+Uq70AsvvOC9fft2p4sXL55zdnYu1gqgrGsCQK9evQKdnJwKtm/ffr0mPx81LadPn3Zp166dX8ntHGGlcsmpOTYRNVGOvkBBXtn7CvIAB1/zxiND48ePh6enJ1auXFnlc44fP65t3759sdGz/fv327700kt3SyarWVlZir1791r/61//arZ27dpmS5cujS2ZrJZ3zYMHD1odPnzYdtGiRVWe2EVUFtawUrnk1BybiJqokKHAhV2mmtWSNawQgZDHJQtNLlQqFdauXYuoqKgqn1NWZ4CdO3eWOQK6c+dO2/Hjx7d0c3PTv/fee7HPPvtselWvGR8fr16xYsWN0NDQcn7rIKoaJqxULh9HLc7HZ6JUB2KYvzk2ETVRLoFAl5dNE6wgACrL+yOuItBlimQTruSma9eu6Nq1a71ce9y4cRnjxo2rUqlBSaNGjZJFSyRq+JiwUrkGhrrjl+gE6A3GUjWs5m6OTURNWLvRpm4AF3YC6XGAg4+kfViJyPyYsFK55NQcm4iaOJdAoOc/pY6CiCTChJUqNDzcG2E+DtgTnYBb6TnwdmAfViIiIjIvJqxUKTk0xyYiIqKmi22tiIiIiEjWmLASERERkawxYSUiIiIiWWPCSkRERESyxoSViIiIiGSNCSsRERERyRoTViIiIiKSNfZhJSKiJi0qNg2ZuXrJ7m9npUa4r6Nk9weACRMm+Ny8edPy999/v1qd8xYuXOi2YcMGlwsXLpxXKpX1FR4RE1YiImraMnP1cLa2lOz+KffyJLs3AJw7d85y48aNrvv27btY3XNnzZqV9O9//9vjs88+c54xY0ZKfcRHBLAkgIiIqEn76KOP3IKDg3N79eqVU91zbWxsxFGjRqWsWLHCvT5iIyrEhJWIiKgB69KlC5566im8/fbb8PHxCbW0tAwPCgpqvWPHDtvKzs3NzRW2bdvm/OSTTxYbHY2OjrZUqVThM2fO9Cy6fdy4cb7W1tYd9u/fry3cNn78+NSYmBjNr7/+al13PxVRcUxYiYgISL4C/LUU2DnD9D35itQRURUUFBTgzJkz+Pnnn3HgwAF8+OGHcV9//fU1g8EgjB8/PiAhIaHCwtLffvvNOisrS/nII49kF90eGhqaN3r06OTVq1c3K7zG7NmzPbZs2eKycePGq0VHYx9++OEca2tr488//2xfPz8lEWtYiYjo9CbgyBcABEClARLOAhd2AV1eBtqNljo6qsD58+eh0+nQq1cv/Prrr1AqlRkA4OzsXDBkyJDg//3vf7bPPvtsennnHzx40EYQBHTu3Dm35L4lS5bciYiIcH7nnXc8goODdZ988onnV199dW3YsGFZRY9TKpUIDg7OOXbsGEdYqd4wYSUiasqSr5iSVSsnQKn+e7tBDxz5EvAKB1wCpYuPKhQVFQUAWLJkCYrO0g8LC9MBQHJycoV/z8fHx6utra0NGo1GLLmvefPm+smTJyeuWrWqmcFgEBYtWhQ7efLktLKu4+zsXHD9+nVNbX4WooqwJICIqCm7sBOAUDxZBe5/Fu7vJ7k6ceIEPD090b1792LbY2Nj1QDg6+ubX9H5eXl5CgsLi1LJaqHAwEBdfn6+0K5du+w33ngjqbzjNBqNUafTCdUMn6jKmLASETVlabGmMoCyqCyB9DjzxkPVEhUVBS8vr1LbN27c6KTRaIwDBgzIBoBdu3bZhoeHt2rdunWIn59f6OLFi90AwMnJqSArK6vMOtcdO3bYzp49u3n79u3vnTx50ubIkSNW5cWRnp6ucnR0LKijH4uoFCasRERNmaMvUFBOH9CCPMDB17zxUJUZjUacPn0aN2/eREHB37nijRs31N99953bhAkTEu3s7IwAMHHixBZbtmyJOX/+/IVr165Fv/jiiykA0KpVK51erxdiYmKKDbEfOHBAO27cuJajR49OPnTo0CUPD4/8uXPnls6M74uLi7MICAjQ1dOPSmT+hFUQBI0gCEcFQTgtCMI5QRDeNXcMRER0X8hQAEZTzWpRBj0AEQh5XIqoqAouXryIe/fuwdLSEhMmTMC+ffuwYsUK5549ewb7+/vrli5dGl94rLu7e/5LL73ku2rVKsesrCyFq6urAQD69++fBQAHDhx4MGEqOjra8oknngjs0aNHxrp162I1Go04Z86c+D///NN+9+7dNiXjSE5OVt68eVPTs2fP7JL7iOqKFCOseQD6iqLYDkB7AIMEQegqQRxEROQSaOoGkJsGZCcBukzT99w0oMuUJjHhys5KjZR7eZJ92VmpKw+yDIUTrv773/8iPT0djz/+OBYsWODdt2/fjP3791/SarUPalNPnjx5Yfr06Yl//fWXbcuWLUNTUlKUABAcHJzftm3be7t27XIAgNjYWNWgQYOCAgICdNu2bbteOJFr6tSpKf7+/rr58+d7l4xjy5Yt9mq1Whw7dmyZE7KI6oLZuwSIoigCKPwtTH3/q9yCbzmLScrGnugExKXlwMdRi4Gh7ghwLfXLJxGRvLUbbeoGcGGnqWbVwcc08toEklUACPd1lDqEGomKioK3tzdCQ0Px008/FW4+XfK4kydPasLCwnRDhw7Natu2re6nn35yVKvVD/7enTx5ctKbb77pk5WVpfD19S24devW2ZLXUKlUuHbt2rmy4vjhhx+cBw8enObu7m6oq5+NqCRJ2loJgqAEcAJASwCfi6J4RIo4amNb1C2sjbwBANCoFTgfn4lfohMwsbsfhoeX+gWUiEjeXAKBnv+UOgqqhhMnTqBjx46VHvfxxx83O3jwoK2VlZXR0tLS+PXXX18vrG0FgJdffjnl008/df/4449dFy5ceLc6MRw8eNDq8OHDtlFRUWUms0R1RZKEVRRFA4D2giA4ANgmCEKoKIrRRY8RBOFFAC8CgK+vvIr+Y5KysTbyBhy0aqiVpqoKWwB6gxFrD95AmI8DR1qJiKjeiKKIU6dOYfbs2ZUe+/3339+saL9arcbXX39949ixY9qKjitLfHy8esWKFTdCQ0PLmblHVDckXThAFMV0QRB+BzAIQHSJfasArAKATp06yapkYE90AgA8SFYLqZUKCPf3v/JISwkiIyKipkAQBGRkZNTZ9fr163evX79+96p73qhRozLrLAiiCkjRJcD1/sgqBEGwAtAfwEVzx1EbcWk50KjLfnSWKgVupeeUuY+IiIiIqk+KLgEeAH4XBOEMgGMAfhVF8adKzpEVH0ct8gqMZe7LKzDC26Hab1WIiIiIqBxSdAk4A6CDue9blwaGuuOX6AToDcZiZQF6gxHi/f1EREREVDe40lUNBLjaYGJ3P6Tn6pGcnYcsnel7eq4eE7v5ccIVERERUR2SdNJVQzY83BthPg7YE52AW+k58HZgH1YiIiKi+sCEtRYCXG3YDYCIiIionrEkgIiIiIhkjQkrEREREckaE1YiIiIikjUmrEREREQka0xYiYiIqEpmzZrlKQhCR71eL3UodWr9+vUOCxYsaFZye+HPK0VMVBwTViIiImrStm/f7vDFF1+USlinTp2atHfv3ga1fHxjxbZWRERE1GDk5uYKVlZWojnuFRAQoA8ICGhcw8kNFEdYiYiIGqirV69CrVbj7bffLrZ93LhxvtbW1h3279+vreq1tm7date+fftWGo0m3NbWtv2jjz4acPr0acuyjj158qRVly5dgqysrDq4urqGzZw509NgMDzYf+bMGcv+/fsHODk5tbO0tAz38PBoO3jw4BZFSwkOHTpk1bdv35Z2dnbtNRpNeHh4eKtffvml2Oo7ha/kjx07punRo0egVqvtMGTIkBZr1qxxFASh45EjR6xKxta7d++WwcHBrQEgOjractiwYf5eXl5tNRpNuLe3d9tx48b5JiUlKQuPHzlypF9ERIRzYmKiWhCEjoIgdPTy8mpb9P41eVaF5549e9ayT58+LbVabQdPT8+2s2fP9qjusyImrERERA1Wy5YtMXnyZHz66adISUkBAMyePdtjy5YtLhs3brzaq1evnKpcZ+vWrXajR48O1Gq1hm+++Sbm448/jr18+bJVnz59Wl2/fl1d8viRI0cG9OnTJ3Pjxo0xw4YNS12+fLnH66+/7lm4//HHHw9MSEiwWLZs2c2IiIjL77zzzm0LCwvRYDAIAHDgwAFt3759W2VkZCiXL19+89tvv41xcHAoGDp0aNBff/1VKskeMWJEyx49emT9+OOPV1977bXEMWPGpNvY2BjWrl3rXPS4uLg4VWRkpP2YMWNS7n9We3t753/44Ydx27Ztuzxnzpz4yMhIu/79+wcWnrNw4cI7vXv3znB0dCzYu3fvxb17917cvHnz1bp6ViNGjGjZu3fvrO+///7qwIED05cuXer52WefPYi7smdFJiwJICKiJuvdXedwPj5T0hhae9rhncfb1Pj8t99+G9999x0++OADeHh4uHzyySeeX3311bVhw4ZlVfUaCxYs8PL29s77888/r6jVppyrT58+2aGhoaGLFy9utnr16ltFj3/mmWeSlyxZkgAAI0aMyMzKylJ+9dVXzebPn39Xr9cLsbGxlhs2bLg6bty4jMJzpkyZklr4z6+//rq3h4dHfmRk5GWNRiMCwMiRIzOCgoLavPvuux579+6NKXq/F198MfGtt95KLLrtscceS9u2bZvT559/fkupNA2Yrl271gkAJk2alAoAgwcPzh48eHB24TmPPvpodnBwcN6gQYOCIyMjrbp3757bpk2bPGdn5wK1Wi3269fvXl0/q1dffTVhxowZKQAwbNiwrMjISNvNmzc7zZgxI+XOnTuqyp4VmXCElYiIqAHz8PDAzJkzsWLFCsyZM6f5okWLYidPnpxWuD81NVXRqVOn4KKvoYvKzMxUnD9/XvvEE0+kFiZgANCqVav88PDwe4cOHbItec748eOLJVRPP/10ak5OjuLEiRNWzZo1K/D29s57++23vZcuXepy9uzZYq/Ks7OzhWPHjtkOHTo0TalUinq9Hnq9HqIoomfPnpnHjh0rdb8xY8akl9w2YcKElMTERPWuXbseHL9p0ybnrl27ZjZv3lwPADqdTpg3b567v79/G41GE25hYdFx0KBBwQBw7tw5TbkPtRw1eVajRo3KKPo5ODg4Nz4+3gIAKntW9DeOsBIRUZNVm5FNOQkMDEReXh7Cw8Oz33jjjaSi+5ycnIzHjx+/VN65SUlJSlEU4eHhUapo0s3NTX/y5Enrktu9vb0Lin729PTUA0BsbKxaoVBgz549l998803PRYsWec2ePVvl5eWVP23atIS5c+cmJSUlqQwGA5YvX+6xfPlyj7JiMhgMKBw1BQBfX99SsQ0YMCDb09Mz/7vvvnMeNmxYVlRUlOb8+fPalStXXi88Ztq0aV7r1q1ze+211+706NEj297e3nDz5k2L5557LkCn01V70K4mz8rNza3Ys7KwsBDz8/MVAFDZs6pufI0ZE1YiIqIGbN++fXjppZfw8MMP4/DhwzZHjhyx6tKlS27h/pkzZ3oCwKeffhpf1vmurq4GQRCQkJBQqv4yMTFR7eDgUFBy+61bt1StW7fOL/wcHx+vBv5OLFu3bp2/bdu2G0ajEYcPH7Zavny527x583z9/f3zBg0alK1QKPDMM88kTpo0KaWsmIomqwCgUChKdQVQKBQYNWpUyurVq5tlZWXFrlmzxlmr1RrHjx+fXnjMjh07nEaMGJHy0Ucf3SnctmvXLmXJa1VVTZ5VZSp6Vk899ZS09SoywpIAIiKiBioqKgrDhw/H5MmT8ccff8DDwyN/7ty5XkWPOXXqlPahhx4qtzbTzs7O2Lp165ydO3c6FhT8nW9dvnzZ4uTJk9bdunUrVQu7YcMGp6Kff/jhByetVmvs1KlTbtHtCoUC3bp1y/3iiy/iAODs2bNWdnZ2xo4dO2afO3dO271795xevXqV+qrqz//888+n5OTkKNavX+8QERHhNHDgwDRbW1tj4X6dTqdQqVTFkt01a9Y4l7yOpaWlmJeXV2lOVJNnVVVlPauaXqsx4ggrERFRA3T16lUMHjwYAwYMwIoVK6BQKDBnzpz4mTNn+u3evdumcLJRdHS0dbdu3W5WdK0FCxbcHj16dGDfvn0Dp0yZkpiVlaV8//33PW1sbAxvvvnm3ZLHr1+/3sVoNKJLly45u3fvttu0aZPLrFmz4p2dnQ1Hjhyxmj59us/IkSPTgoKCdAaDQVi3bp2zUqkUBwwYkAUAy5YtixswYEBwz549AydMmJDs5eWlT0pKUp04ccLaYDBg5cqVt6vyDMLCwvLCwsLuvfvuu96JiYnq5557rtiIbe/evTMiIiKcP/jgg9ygoKC8rVu3Opw4ccKm5HVCQkJyf/jhB5cPP/zQtWvXrvesrKzEzp0755Y8ribPqiJVeVZkwoSViIiogUlISMCAAQMQEhKCjRs3QqEwDQ5OnTo1Zfny5e7z58/3Hjx48MWYmBi1IAiiv79/hU09R40alblp06YrixYt8pw0aVKAWq02du7cOWvZsmW3/Pz8Sp0bERFx9dVXX/X99NNPPW1sbAzTp0+/U/ja3dvbW+/l5ZW/cuXKZnfv3lVbWFiIQUFBuZs2bbras2fPHADo0aNHzoEDBy68/fbbnvPmzfPNzs5WOjo6FrRp0ybnpZdeqlbt5pgxY1Lmz5/v6+bmph8yZEixJG/VqlVxL7zwgrBkyRIvAOjTp0/G+vXrr/Xp0yek6HEzZsxIPnr0qPXixYu9srKylJ6envm3b98+WxfPqiJVeVZkIoiiWRaLqJVOnTqJx48flzoMIiJqwC5cuICQkJDKD2zYThT9sGHDBodvvvnG5ffffy+3ryiRnJw+fdqlXbt2fiW3s4aViIiokTp+/Li2ffv2HKmjBo8lAURERI1UeZ0BiBoajrASERERkawxYSUiIiIiWWPCSkRERESyxoSViIiIiGSNCSsRERERyRoTViIiIiKSNSasRERERCRrTFiJiIiISNaYsBIRERGRrDFhJSIiIiJZY8JKRERERLLGhJWIiIiIZI0JKxERERHJGhNWIiKiBm769OkYMmRItc+7ffu2atWqVY71EFKV77Nw4UK3oKCg1gaDoUrXmjVrlqeXl1fbOg2wDOvXr3dYsGBBs/q+T03NmjXLUxCEjuV9rsm19Hp9te5Z3XvU5t8bE1YiIqIGLCYmBl9++SUWLFhQ7XNXrFjh8vnnn9d7UlbRfWbNmpWUmpqq/uyzz5zrO47q2L59u8MXX3wh24R16tSpSXv37r0odRzmwoSViIioJpKvAH8tBXbOMH1PviJJGJ9++inatWuHTp06VfvcvXv32k+ZMiWxHsKq8n1sbGzEUaNGpaxYscK9vuNoTAICAvT9+vW7J3Uc5sKElYiIqLpObwIiXgAu/AQkXwYu/Gz6fHqTWcPIy8vDhg0bMHbs2GLbo6OjLVUqVfjMmTM9i24fN26cr7W1dYf9+/drU1NTFXfu3LGYNGlSWk3u3b1798B27dq1Krn96NGjViqVKvyLL75wAoCq3Gf8+PGpMTExml9//dW6JrEAf7+uPnv2rGWfPn1aarXaDp6enm1nz57tUbLcYOvWrXbt27dvpdFowm1tbds/+uijAadPn7Ys3D9y5Ei/iIgI58TERLUgCB0FQehY0evs6Ohoy2HDhvl7eXm11Wg04d7e3m3HjRvnm5SUpCw8Zs2aNY6CIHQ8cuSIVcnze/fu3TI4OLh1Va9V9Oet6JlU9VqFTp48adWlS5cgKyurDq6urmEzZ870rKxU49ChQ1Z9+/ZtaWdn116j0YSHh4e3+uWXX2wqPKkGmLASERFVR/IV4MgXgJUTYOMGaOwAG1fT5yNfmnWk9fDhw0hPT0fPnj2LbQ8NDc0bPXp08urVq5slJCQoAWD27NkeW7Zscdm4cePVXr165fz0009248ePT7K0tBRrcu+uXbtmX7hwQZubmysUbjMajXj55Zd9O3TocO/ll19OBYCq3Ofhhx/Osba2Nv7888/2NYmlqBEjRrTs3bt31vfff3914MCB6UuXLvUsWm6wdetWu9GjRwdqtVrDN998E/Pxxx/HXr582apPnz6trl+/rgaAhQsX3undu3eGo6Njwd69ey/u3bv34ubNm6+Wd8+4uDi1t7d3/ocffhi3bdu2y3PmzImPjIy069+/f2DhMWPGjEm3sbExrF271rnEuarIyEj7MWPGpFT1WlVV3WuNHDkyoE+fPpkbN26MGTZsWOry5cs9Xn/9dc+yjgWAAwcOaPv27dsqIyNDuXz58pvffvttjIODQ8HQoUOD/vrrL211462Iqi4vRkRE1Ohd2AlAAJTq4tuVatP2CzuBnv80SyiHDx+GIAgICwsrtW/JkiV3IiIinN955x2P4OBg3SeffOL51VdfXRs2bFgWABw7dsx6/vz5d2t67549e2YvW7ZMOHjwoLbw1fTnn3/ufPr0aetDhw6dLzyuKvdRKpUIDg7OOXbsWI1HWAu9+uqrCTNmzEgBgGHDhmVFRkbabt682alw24IFC7y8vb3z/vzzzytqtenfYZ8+fbJDQ0NDFy9e3Gz16tW32rRpk+fs7FygVqvFqrx2Hzx4cPbgwYOzCz8/+uij2cHBwXmDBg0KjoyMtOrevXuuVqsVH3vssbRt27Y5ff7557eUStMg59q1a50AYNKkSalVvVZVn0V1r/XMM88kL1myJAEARowYkZmVlaX86quvms2fP/+ui4tLqaHW119/3dvDwyM/MjLyskajEQFg5MiRGUFBQW3effddj71798ZUNdbKcISViIioOtJiAZWm7H0qSyA9zmyhxMfHw87ODhYWFqX2NW/eXD958uTEdevWuc2ZM6f5okWLYidPnvzgtfyKFStue3h4FNT03o888sg9pVKJyMhIawBITk5WLliwwPu5555Leuihh3TVvY+zs3NBYmJi6R+kmkaNGpVR9HNwcHBufHy8BQBkZmYqzp8/r33iiSdSC5NVAGjVqlV+eHj4vUOHDtnW5J46nU6YN2+eu7+/fxuNRhNuYWHRcdCgQcEAcO7cuQf/sUyYMCElMTFRvWvXrgf32bRpk3PXrl0zmzdvrq/OteoyrkLjx49PLfr56aefTs3JyVGcOHGiVBlDdna2cOzYMduhQ4emKZVKUa/XQ6/XQxRF9OzZM/PYsWM1epblYcJKREQmMplEJHuOvkBBXtn7CvIAB1+zhaLT6WBpaVnu/sDAQF1+fr7Qrl277DfeeCOpsuulpqYqOnXqFFyVFlP29vbG4ODgnKNHj9oAwGuvvealUCjEjz/++HZN7qHRaIw6nU4o57Qqc3NzK5YcW1hYiPn5+QoASEpKUoqiCA8Pj1L9m9zc3PTp6ek1evM8bdo0r6VLl3o++eSTqZs3b77yxx9/XPj2229jAECn0z3ItQYMGJDt6emZ/9133zkDQFRUlOb8+fPacePGpVT3WnUZVyFvb+9iz87T01MPALGxseqSxyYlJakMBgOWL1/uYWFh0bHo13fffeeWmZmprGqrsqpgSQAREZkmCx35AoBgGj1MOAtc2AV0eRloN1rq6OQlZKjp2Rj0xcsCDHoAIhDyuNlCcXZ2Rnp6epn7duzYYTt79uzm7du3v3fy5EmbI0eOWHXp0qXC18lOTk7G48ePX6rq/R966KHsX3/91eHAgQPa77//3nXFihXXnZycjDW5R3p6usrR0bHSkdhly5bFL1u2LL6qMRbl6upqEAQBCQkJpRKwxMREtYODQ41GnHfs2OE0YsSIlI8++uhO4bZdu3aVmtikUCgwatSolNWrVzfLysqKXbNmjbNWqzWOHz8+vbrXqsu4Ct26dUvVunXr/MLP8fHxagDw9fUtleA7OzsbFAoFnnnmmcRJkyallNwPmEo9CtXm3xvAEVYiIpLRJKIGwSXQlMjnpgHZSYAu0/Q9Nw3oMsW030xatWqF/Px83Lp1q9j2AwcOaMeNG9dy9OjRyYcOHbrk4eGRP3fuXK/Krjdz5kzPkp0FKtKjR4/s+Ph4i4kTJ/p36NAh+5VXXkmt7Jzy7hEXF2cREBCgK+ucumJnZ2ds3bp1zs6dOx0LCv7OTS9fvmxx8uRJ627dumUVbrO0tBTz8vKqlCfpdDqFSqUqNqlszZo1ZfaVff7551NycnIU69evd4iIiHAaOHBgmq2t7YMkvzrXqsu4AGDDhg1ORT//8MMPTlqt1tipU6dSv+jY2dkZO3bsmH3u3Dlt9+7dc3r16lXqqyYxl4cJKxFRU1eVSURUXLvRwIhVQMhjgGuw6fuIVWYfje7VqxcA4OjRow+2RUdHWz7xxBOBPXr0yFi3bl2sRqMR58yZE//nn3/a7969u8J2Q6dOndI+9NBD9wDg0qVLFoIgdJw1a1a5CWy/fv2yAeD69euazz//PLYqMRe9R6Hk5GTlzZs3NT179swu77xCTz31VPParJi0YMGC27GxsZq+ffsGfv/99/ZfffWV04ABA4JsbGwMb7755oPJYSEhIbkZGRnKDz/80PXPP//UHj16tFQdZ6HevXtnREREOH/wwQeuERERdmPHjvU9ceJEmc86LCwsLyws7N67777rfefOHYvnnnsupabXqkx1r7V+/XqXefPmuW/bts3uxRdf9N60aZPLlClTEpydnct8t79s2bK4c+fOaXv27Bm4atUqx59//tlm3bp1DtOmTfN65ZVXiv2CVNt/b0xYiYiaOhlNImpQXAJN3QAe/9T03Ywjq4X8/PzQuXNn7Nq1CwAQGxurGjRoUFBAQIBu27Zt1wtfyU6dOjXF399fN3/+fO+KrhcdHW3drVu3HADIyspSAIC7u3u563Xa2dkZ1Wq1+OyzzyZWVm5Q1j0KbdmyxV6tVotjx46ttCeswWAQCgoKalzrOmrUqMxNmzZdyczMVE6aNClg9uzZvgEBAbl//PHHRT8/vwc/64wZM5KHDBmSunjxYq8+ffqEDB8+vGV511y1alVcv379MpYsWeI1YcKEFtnZ2cr169dfK+/4MWPGpCQmJqrd3Nz0Q4YMyarNtSpS3WtFRERc/f333+2efvrplhEREc7Tp0+/U7ScoKQePXrkHDhw4IKTk5Nh3rx5vsOHDw+aO3eu77lz56x69+5d7JeP2v57E0SxRu3XzKpTp07i8ePHpQ6DiKhx+mupqfG9jWvpfdlJQMgQoOcs88dVxy5cuICQkBCpw6hz69atw4wZM3Dnzh1otdoTNb1OTEyMumvXriFJSUlnAOD//u//XBYvXuwVGxt7tugr66JeeOEF7+3btztdvHjxXHmjcBXdo1CvXr0CnZycCrZv3369pvFT43D69GmXdu3a+ZXczhFWIqKmLmQoAOP9SUNFSDCJiKpv/Pjx8PT0xMqVK2t1nUOHDlmHhoY+GPncv3+/7UsvvXS3ZLKalZWl2Lt3r/W//vWvZmvXrm22dOnS2Kokq2XdAwAOHjxodfjwYdtFixbVeEIONX7sEkBE1NQVTiI68iVMXQIs77dtEs0+iYiqT6VSYe3atYiKiqrVdY4fP65t3779g2Ry586dZY527ty503b8+PEt3dzc9O+9917ss88+m17TewCmmegrVqy4ERoaWk6vMCKWBBARUaHkK6YJVulxgIOPaeS1ESWrjbUkoIQalwQQyUF5JQEcYSUiIpPCSURERDLDGlYiIiIikjUmrEREREQka0xYiYiIiEjWmLASERERkaxx0hURUS3EJGVjT3QC4tJy4OOoxcBQdwS41mgVRSIiKgcTViKiGtoWdQtrI28AADRqBc7HZ+KX6ARM7O6H4eEVroBJRETVYPaEVRAEHwDfAWgGQASwShTF5eaOg4ioNmKSsrE28gYctGqolabqKlsAeoMRaw/eQEebZPje3QekxQKOvo2upykRkTlJMcJaAOCfoihGCYJgC+CEIAi/iqJ4XoJYiIhqZE90AgA8SFYLqZUK9Nb9Dsud/wXsNIBKAyScBS7sMq0m1W60FOESETVoZp90JYriHVEUo+7/cxaACwC8zB0HEVFtxKXlQKMu/UdoM30cntDtQLpgA9i4ARo7wMYVsHIyLX2afEWCaImIGjZJuwQIguAHoAOAI2Xse1EQhOOCIBxPSkoye2xERBXxcdQir8BYanv4vQMwArBUWxbfoVQDEExLnxIRUbVIlrAKgmAD4D8AZoqimFlyvyiKq0RR7CSKYidXV1fzB0hEVIGBoe4QRVPNalFO+jvIhwVc7SxLn6SyBNLjzBQhEVHjIUmXAEEQ1DAlqxtFUYyQIgYiotoIcLXBxO5+WHvwBgQAlioF8gqMiBNd0d4mDtYWZfzxWpAHOPiaPVaqRNxRQJch3f019oBPZ8luP2HCBJ+bN29a/v7771eres7ChQvdNmzY4HLhwoXzSqWyPsMjAiBNlwABwDcALoiiuMzc9yciqivDw70R5uOAPdEJuJWeA28HLR7xmgSnP84ABv39MoD7DHoAIhDyuGTxUjl0GYC1i3T3v5cs2a3PnTtnuXHjRtd9+/ZdrM55s2bNSvr3v//t8dlnnznPmDEjpb7iIyokxQhrdwDPADgrCMKp+9vmi6L4XwlioQaCzdlJrgJcbfDKIy2Lb8x92TTBCoKpDKAgD4AIdJnC1lYkKx999JFbcHBwbq9evXKqc56NjY04atSolBUrVrgzYSVzkKJLwAFRFAVRFMNEUWx//4vJKpVrW9QtvPbjKfwSnYCYxGzsOZeA1348hW1Rt6QOjahs7UYDI1YBIY8BrsGm7yNWsaUV1YsuXbrgqaeewttvvw0fH59QS0vL8KCgoNY7duywrei83NxcYdu2bc5PPvlksYQzOjraUqVShc+cOdOz6PZx48b5Wltbd9i/f78WAMaPH58aExOj+fXXX63r/qciKq7KCasgCBaCIHQVBGGEIAjjBEEYeH+WP1G9Kdqc3dXWErYaNVxsLOGgVWPtwRuIScqWOkSisrkEAj3/CTz+qek7R1apHhQUFODMmTP4+eefceDAAXz44YdxX3/99TWDwSCMHz8+ICEhodwC099++806KytL+cgjjxT7gzQ0NDRv9OjRyatXr25WeP7s2bM9tmzZ4rJx48arhaOxDz/8cI61tbXx559/tq/fn5KokpIAQRCUAIYDmAygNwALAEKRQ0RBEG4D+AHA16IoVrlgm6gqKmrOLtzfX+p1LJEcJF8xtbDiSldUj86fPw+dTodevXrh119/hVKpzAAAZ2fngiFDhgT/73//s3322WfTyzr34MGDNoIgoHPnzrkl9y1ZsuRORESE8zvvvOMRHBys++STTzy/+uqra8OGDcsqPEapVCI4ODjn2LFjHGGlelduwioIwigA7wPwAbAHwL8AnASQBCAXgBMAfwBdYEpqZwmCsA7Av0RRvFu/YVNTUV5zdsA0K/tWerXKrojM4/Qm4MgXMNWwaoBbx4CjXwPOLYGAR5i8Up2JiooCACxZsgRFZ+uHhYXpACA5OVkFAEajEXPnzvV4//3376hUpr/64+Pj1dbW1gaNRiOWvG7z5s31kydPTly1alUzg8EgLFq0KHby5MlpJY9zdnYuuH79uqZefjiiIioqCfg3gM8BuIui+IQoiktFUfxNFMWzoiheFUXxqCiKm0RRnCWKYhCAHgCcAbxojsCpaSivOTsA5BUY4e2gNXNERJVIvmJKVq2cTCtd5WUBGXFAXjZwOwo4+x8g4gVTUktUSydOnICnpye6d+9ebHtsbKwaAHx9ffMBICoqSrNlyxbnwmQVAPLy8hQWFhalktVCgYGBuvz8fKFdu3bZb7zxRpkr+Gg0GqNOpxPK2kdUlypKWFuIovipKIrpVbmQKIpHRFEcAeDjOomMCOU3Z9cbjBDv7yeSlQs7AQimllb52UD6TdM/qzWAQgUY87lMK9WZqKgoeHmVXt1848aNThqNxjhgwIDsEydOaP7xj38EZWRkKFu1atX6+eef9wEAJyengqysrDJrXHfs2GE7e/bs5u3bt7938uRJmyNHjliVdVx6errK0dGxoE5/KKIylJuwiqKoq8kFa3oeUUxSNlb+fhVvRJzByt+vIiYp+0Fz9vRcPZKz85ClM31Pz9VjYjc/trYi+UmLNZUBAEDWXQACINz/o1ZQAHodl2mlOmE0GnH69GncvHkTBQV/54w3btxQf/fdd24TJkxItLOzM3bs2FE3ZMiQtGnTpiVcvHjx/DfffBMHAK1atdLp9XohJiZGXfS6Bw4c0I4bN67l6NGjkw8dOnTJw8Mjf+7cuaWzYgBxcXEWAQEB/Huf6l25CasgCJOqcgFBEJwEQVhfdyFRU1RR66rh4d74ZHR7DGzjjpbNbDCwjTs+Gd0ew8O9pQ6bqDRH3/t9VwHoc/9OVgFANALq+wNVXKaVaunixYu4d+8eLC0tMWHCBOzbtw8rVqxw7tmzZ7C/v79u6dKl8YXHnj59Wtu5c+diRf/9+/fPAoADBw48mDQVHR1t+cQTTwT26NEjY926dbEajUacM2dO/J9//mm/e/fuYiMEycnJyps3b2p69uzJdi1U7yrqErBaEITRAF4QRTG2rAMEQRgHYBmAMhbNJqqaoq2rCrsB2ML02n/twRsI83Eouzk7kRyFDAUu7DKtbKW2MtWwAoAoAhABm2amz1ymVT409pKuNgVNzbpCFU64+u9//4t58+bh8ccfh5WVlfeQIUNSP/3009tarVYETK2vLl68qO3evfu9oucHBwfnt23b9t6uXbscnnvuufTY2FjVoEGDggICAnTbtm27XjiJa+rUqSnLly93nz9/vvfgwYMfrIi1ZcsWe7VaLY4dO7bUZCyiulZRwjocwEoA5wRBmC+K4orCHYIg+AD4CsAgALsAvFKvUVKjxtZV1Ki4BAJd7q90pbQAxAJAL5oaAjr6ARbWXKZVbnw6Sx1BjURFRcHb2xuhoaH46aefCjefLnnczZs31Wq1WrS3ty81g3Xy5MlJb775pk9WVpbC19e34NatW2dLHqNSqXDt2rVzJbf/8MMPzoMHD05zd3c31MXPQ1SRimpYdwBoDWALgOWCIBwQBKGNIAjTAZwHEA5gzP0OArfNEy41RmxdRY1O4UpXbUcCnuGApQ1g7wNY2ADZSUBuGpdppVo7ceIEOnbsWOlxzZs317dv3z47MDCwzaRJk3yK7nv55ZdTXF1d9R9//LFrde598OBBq8OHD9suWrQovvKjiWqvwoUDRFHMADBJEIQfAKwBcOb+rm8B/FMURb4GoFrzcdTifHwmylpDkK2rqMEqXOmq5z//XkQgPQ5w8GEfVqo1URRx6tQpzJ49u9JjVSoV/vjjjzIX9lGr1fj6669vHDt2rFp/0MbHx6tXrFhxIzQ0NK865xHVVIUJKwAIgqAA8BAAVwD5MK12lQNAX7+hUVMxMNQdv0QnQG8wFisLYOsqajQKk1eiOiIIAjIyMurkWv369bvXr1+/e5Uf+bdRo0Zl1snNiaqooj6sEAShPYDjABYC+BKAG4A5ACbAVNs6uJ7joyaArauIiIioIhUtzboYwGwAVwH0FEXx0P1dSwVB2A5gFYCf7pcLzBBFMaW+g6XGa3i4N8J8HLAnOgG30nPg7aDFwFB3JqtERERUYUnAbAAfAXhPFMX8ojtEUYwB0E8QhBfuH9MfQLN6i5KaBLauIiIiorJUlLB2EkWxVHuLokRR/FoQhJ8BfF63YRERERERmZSbsFaWrBY5Lh6mnq1ERERERHWuwklXRERERERSq2jS1du1ubAoigtrcz4REREREVBxDevEWlxXhKkVFhHVpcIG9GmxgKMvG9ATEVGTUFENq785AyGiSpzeBBz5AoAAqDRAwlngwi7TuvXtRksdHRERUb2pqCRAEEVRNGcwRI1VTFI29kQnIC4tBz6ONegxm3zFlKxaOQFK9d/bDXrgyJeAV3jjH2k15+gyR7KJiGSlopKAw4IgTBRF8bzZoiFqhLZF3cLayBsAAI1agfPxmfglOgETu/theLh31S5yYScAoXiyCtz/LJj2N+alP2syulzTpJMj2UREslNRwmoHIEoQhEUA3hdF0WCmmIgajZikbKyNvAEHrRpqpakphy0AvcGItQdvIMzHoWojrWmxpuSpLCpLID2u7oKWm5qMLtc06azoXpGfAndOAfk5HHUlIjKzitpatQPwfwDeAnBMEIT2ZomIqBHZE50AAA+S1UJqpQJCkf2VcvQFCvLK3leQBzj41iJKmavK6HJRRZNOGzdAYwfYuJo+H/nStL+697qXBKTdAM5tB5IvAxd+BiJeMCXGRA3A7du3VatWrXKs6+vOmjXLUxCEjnq9vq4vLan169c7LFiwoMwVPAt/ZnPH1NSVm7CKopgviuK/ADwEwADgiCAIiwRB0AiCoCj5ZbaIiRqQuLQcaNRl/+9hqVLgVnpO1S4UMhSA0TTSV5RBD0AEQh6vVZyyVt3R5eomuJXdK/+eKVlVWgCCsnoJMJFMrFixwuXzzz/nEupVtH37docvvviizOc1derUpL179140d0xNXaWJpiiKZwB0gWm09Q0A9wDoS3zl12OMRA1WW8tEDMn8Ec8kf4LB6T+gmf7v5CqvwAhvB23VLuQSaHqdnZsGZCcBukzT99w0oMuUxv1qurqjy7UpnyjrXtn3R8FFEVBb/b29KgkwkUzs3bvXfsqUKYlSx1FfcnNzBXPdKyAgQN+vX7975rofmVR1ZHQEgOcBpAH4AKYeq0W/3quX6IgastObMOLGQnTNOwj3/JsIzzmAF5LeR5esvdAbjBABDAx1r/r12o0GRqwCQh4DXINN30esavwTgao7ulyb8omy7pWfC0Aw3cvWrfjxjb1+mGTv6tWrUKvVePvt4mv9jBs3ztfa2rrD/v37tampqYo7d+5YTJo0Ka061966datd+/btW2k0mnBbW9v2jz76aMDp06ctyzr25MmTVl26dAmysrLq4OrqGjZz5kxPg+HvqS9nzpyx7N+/f4CTk1M7S0vLcA8Pj7aDBw9uUbKU4NChQ1Z9+/ZtaWdn116j0YSHh4e3+uWXX4oV+he+kj927JimR48egVqttoOTk1N7QRA6HjlyxAol9O7du2VwcHDrws/R0dGWw4YN8/fy8mqr0WjCvb29244bN843KSlJCQAjR470i4iIcE5MTFQLgtBREISOXl5ebUvev7rPqvC8s2fPWvbp06elVqvt4Onp2Xb27NkeNXlWTU1Fk64gCII7gJUAhgH4D4Cpoig22t/QiOrM/TpKjZ0LbBUFiEvNBQBYCAXok/EfnDT4Y3ivHtVrbQWYRlIbczeAshSOLh/5EqZJVJb3E1Kx7NHlkKGmCVYGfemJU5WVT5R1L9EIGPWAUwBgUeLfV2OvHybZa9myJSZPnoxPP/0UM2bMgLOzM2bPnu2xZcsWl82bN1/p1atXznfffecwfvz4JEtLyyq3qty6davd6NGjA7t06ZL5zTffxGRlZSmXLFni2adPn1ZRUVHn/f39i2VPI0eODBg7dmzy3LlzE3bv3m23fPlyD4VCgWXLlsUDwOOPPx5oZ2dnWLZs2U1XV9eCuLg4i927d9sbDAZBrVaLAHDgwAFt//79g1u3bp2zfPnym1qt1vjVV1+5Dh06NGjfvn0Xe/bsWayGasSIES3HjRuXPGfOnIT8/Hxh4sSJLdauXevcpUuXW4XHxMXFqSIjI+3ffPPNotvU3t7e+aNGjYpzdnYuuHLliuWyZcs8+vfvrz116tTFhQsX3klJSVGdOXPGesuWLVcBQKPRGOvqWY0YMaLl2LFjk2fOnHl3x44dDkuXLvX08fHJnzFjRkpVn1VTVFEf1skAPgKQB+BJURT/Y7aoiBq6InWUHvZq2FmpkZSZh9wCFeyMeVjU6jqcw8dIHWXD0W60qRvAhZ2mEU0Hn/Jn6Vc3wa3sXp4dgGu/mSZwFdUU6oebgt3zTF0kpOTeFhj8QY1Pf/vtt/Hdd9/hgw8+gIeHh8snn3zi+dVXX10bNmxYFgAcO3bMev78+Xerc80FCxZ4eXt75/35559X1GrTL359+vTJDg0NDV28eHGz1atX3yp6/DPPPJO8ZMmSBAAYMWJEZlZWlvKrr75qNn/+/Lt6vV6IjY213LBhw9Vx48ZlFJ4zZcqU1KLXeP311709PDzyIyMjL2s0GhEARo4cmREUFNTm3Xff9di7d29M0eNffPHFxLfeeuvBINrWrVvTtm3b5vT555/fUiqVAIC1a9c6AcCkSZMe3Gvw4MHZgwcPzi78/Oijj2YHBwfnDRo0KDgyMtKqe/fuuc7OzgVqtVqsyqv/6j6rV199NaEwOR02bFhWZGSk7ebNm51mzJiRcufOHVVVnlVTVNEI6yoAGwHMEEWxyT8oomopUUdpbaGCtcv9/910APTV+ruDgOqNLlcnwa3KvU5vqnkCXJe4oAGVwcPDAzNnzsTSpUtRUFDQfNGiRbGTJ09+8Pp/xYoVt6tzvczMTMX58+e1r7766p3CBAwAWrVqlR8eHn7v0KFDtiXPGT9+fLE84emnn07dtGmTy4kTJ6z69++f7e3tnff22297JyQkqAcMGJDVtm3bYnU72dnZwrFjx2xfffXVO0qlUiz6+rtnz56Z27dvdy55zzFjxqQX/TxhwoSUTZs2uezatcu2MFnftGmTc9euXTObN2/+4II6nU5YsGBBs02bNjnfuXPHMi8v70H967lz5zTdu3fPrc9nNWrUqIyin4ODg3Ojo6O1ANCsWbOCyp5VU1VRwvq4KIo/my0SogamwtWrHH3LH7Xha2TzqMvyidomwHWBCxrUj1qMbMpJYGAg8vLyEB4env3GG28kVXZ8amqqYsCAAYFHjhy5VDgaWSgpKUkpiiI8PDxKFU26ubnpT548aV1yu7e3d0HRz56ennoAiI2NVSsUCuzZs+fym2++6blo0SKv2bNnq7y8vPKnTZuWMHfu3KT791QZDAYsX77cY/ny5R5lxWwwGFA0Vl9f32LxDRgwINvT0zP/u+++cx42bFhWVFSU5vz589qVK1deL3rctGnTvNatW+f22muv3enRo0e2vb294ebNmxbPPfdcgE6nq1bXo5o8Kzc3t2LPysLCQszPz1cAQFWeVVNVbsLKZJWofJWuXlWbOkqSJynrh7k0L1Vg3759eOmll/Dwww/j8OHDNkeOHLHq0qVLhaOETk5OxuPHj18qa5+rq6tBEAQkJCSoS+5LTExUOzg4FJTcfuvWLVXr1q0fdAyKj49XA38nla1bt87ftm3bDaPRiMOHD1stX77cbd68eb7+/v55Tz31VKazs7NBoVDgmWeeSZw0aVJKWXGVTKwVCoVY4jNGjRqVsnr16mZZWVmxa9ascdZqtcbx48enFz1ux44dTiNGjEj56KOP7hRu27VrV/GLV1FNnlVlKntWNYmzMSj3NwlBEE4JgjBcEIQqtYoQBMFbEIR/C4Iwp+7CI5KfoqtXudpawlajhouNJRy0aqw9eAMxSdlNuw0V1b3a9JalRi0qKgrDhw/H5MmT8ccff8DDwyN/7ty5XpWdN3PmTM+ZM2d6lrXPzs7O2Lp165ydO3c6FhT8nW9dvnzZ4uTJk9bdunXLKnnOhg0bnIp+/uGHH5y0Wq2xU6dOxRJnhUKBbt265X7xxRdxAHD27Fmrwnt27Ngx+9y5c9ru3bvn9OrVq9RXVZ7H888/n5KTk6NYv369Q0REhNPAgQPTbG1ti02Y0ul0CpVKVSzZXbNmTbGSA0tLSzEvL6/S0daaPKuqKu9ZNVUVlQR8B+BrAJ8JgrAZwF8ATgNIgmkiliOAFgA6A3gcQG8A+wB8Vp8BE0mtKqtXvfJIS3m8RqbGoSkvzUvlunr1KgYPHowBAwZgxYoVUCgUmDNnTvzMmTP9du/ebVN0YlFJp06d0r7wwgvlvmJesGDB7dGjRwf27ds3cMqUKYlZWVnK999/39PGxsbw5ptvlirCX79+vYvRaESXLl1ydu/ebbdp0yaXWbNmxTs7OxuOHDliNX36dJ+RI0emBQUF6QwGg7Bu3TpnpVIpDhgw4EFCt2zZsrgBAwYE9+zZM3DChAnJXl5e+qSkJNWJEyesDQYDVq5cWWktblhYWF5YWNi9d9991zsxMVH93HPPlRqt7d27d0ZERITzBx98kBsUFJS3detWhxMnThRrARISEpL7ww8/uHz44YeuXbt2vWdlZSV27ty5zFHr6j6rilT1WTVFFZUELBME4RsAk2HqwToDQMl2CgJMyesOAP1EUfyzvgIlkotqrV7VFNtQUd1jTTSVkJCQgAEDBiAkJAQbN26EQmH6M2nq1Kkpy5cvd58/f7734MGDy12NKTo62rpbt243y9s/atSozE2bNl1ZtGiR56RJkwLUarWxc+fOWcuWLbvl5+dXql4zIiLi6quvvur76aefetrY2BimT59+p/CVu7e3t97Lyyt/5cqVze7evau2sLAQg4KCcjdt2nS1aKuqHj165Bw4cODC22+/7Tlv3jzf7OxspaOjY0GbNm1yXnrppSrXb44ZMyZl/vz5vm5ubvohQ4aUSvJWrVoV98ILLwhLlizxAoA+ffpkrF+//lqfPn1CCo+ZMWNG8tGjR60XL17slZWVpfT09My/fft2mf8TVvdZVaSqz6opEkSxai29BEFoDtOKV54ANABSAFwEcFQUxXqdwdapUyfx+PHj9XkLoipb+ftV7DmXABeb0v2zk7PzMLCNu2mElaiuJF8BIl4ou4Y1N820gARH7it14cIFhISEVH5gw3aisgNiYmLUXbt2DUlKSjpjjoCIquP06dMu7dq18yu5vcKFA4oSRfEmgHJ/GyNqKgaGuuOX6AToDcZiZQE1Wr2KqCpq21uWqIhDhw5Zh4aGNunROmp4KlvpygXAeAAtAaQD2CqK4qn6D4tIvgJcbTCxux/WHrwBAaYygLwCU7I6sZtf9VevIqoK1kRTHTl+/Li2ffv2TFipQalopatgAPsBuBbZPFcQhFGiKO6o98iIZGx4uDfCfBywJzoBt9Jz4O1Qog8rUX1gTTTVgU8//TRe6hiIqquiEdZFMK3J0wfAMZg6AqwCsAymSVZETVqAqw1rVYmIiMygoh5jXQC8LYriflEUc0VRPAfgJQB+giC4VnAeEREREVGdqWiE1QtAyVUwLsHUysoTpn6sRFRDFS7tSg1P8hVTfWlarKkNFetLiYjqTEUJqwDAUGJb4WoR1Vprl4iKq3RpV2pYTm8yLZ0KwdTgP+GsaWneLi+bJkuRbIiiiCou4EhEZna/1WqZ/VYra2v1riAIyUU+F/5f/p4gCKnF7yE+V/MQiZqOoku7FrbFsoWpLdbagzcQ5uPAkdaGJPmKKVktq0fqkS9NM/s50ioLarUaubm50Gq1UodCRGW4d++eVhCE62XtqyhhjQVQVoflmwDalNhWtdUHiKjqS7tSw3BhJwCheLIK3P8smPZzZr8suLm54fbt2/Dy8oKVlRVHWolkQBRF6PV6VWZmpk1CQoJRr9d/VNZxFS3N6ldv0RE1YdVa2pXkLy3WVAZQFpWlqWcqyYKdnR0AID4+Hnp9tVbMbDDy8/NdpI6BqJqMgiCkiaK4v6Cg4P2OHTteKeugKq90RUR1w8dRi/PxmbAtY19egRHeDnxd2aA4+ppqVstSkAc4+Jo3HqqQnZ3dg8S1kfKTOgCi+sDJU0RmNjDUHaJoqlktSvKlXZOvAH8tBXbOMH1PLvOXXCopZCgAo6lmtSiDHoAIhDwuRVRERI0KR1iJzEyWS7tylnvNuQSantORL2F6fpamkVWIQJcpnHBFRFQHhPstBGStU6dO4vHjx6UOgwhA3fVPLbyO5Eu7Jl8BIl4oe5Z7bhowYpU8ky659T0tjCc9DnDwkT4eaqo4k4waJSasRNVQsn9qXoERooiG3T/1r6XAhZ8AG7fS+7KTgJDH5DfLveSIcEEeACNHhImYsFIjxZIAoipqtP1TG9osd/Y9JSJqcjjpiqiKqtI/tUFy9L0/QlkGOc5yr0rfUyIialSqlbAKgmAhCMIKQRDa11M8RLLVaPun1tUsd3N1GWhoI8JERFRr1UpYRVHMB3AHwIr6CYdIvnwctcgrMJa5r0H3Ty2c5Z6bZqpZ1WWavuemVX2W++lNpolbF34Cki8DF342fT69qe7jbWgjwkREVGs1KQn4DkA3QRCa1XUwRHIm2/6pdaHdaFM3gJDHANdg0/cRq6o2galoTamNG6CxA2xcTZ+PfFn3I63se0pE1OTUZNJVJkyzEJ0B3K3uyYIgrAEwBECiKIqhNbg/kSRk2T+1LrkE1qwbQFVqSqt53Qpbh7HvKRFRk1OThLUrgHwAN2p4z3UAPoNppJaoQRke7o0wHwd59E81p4p6ntZxTWnJ1mHn4zPxS3RC8dZh7UabugGw7ykRUZNQk4R1JIDfRVGs0QwTURT3C4LgV5NziUqqqyb+1RHgaoNXHmlZr/eQlcpWwXL0NW0rSzVrSqvVOqymI8JERNTg1CRhHQxgYV0HQlRdVRqJo9qpSs/TkKGmBNagL31MNWtKq9I6rKa/LNT6lxu5raxFRNSE1CRh9QBwvq4DKUkQhBcBvAgAvr6c9UvFNdom/nJT1frUOqopra/WYbX+5aayUWYiIqpXNUlYrwLwA3CwbkMpThTFVQBWAaalWevzXtTw1OdIXEWkKEGQVFXrU+uoptTHUYvz8ZmwLWNfTVuH1fqXG66sRUQkuZq0tVoFYKYgCFyvmCQjRRP/bVG38NqPp/BLdAJiErOx51wCXvvxFLZF3arze8lGdXqeFtaUPv6p6XsNkrj6aB1W6xXKuLIWEZHkapKwrrj/VaM+rIIg/ADgEIBgQRBuCYLwfE2uQ02buZv4Fx2lc7W1hK1GDRcbSzho1Vh78AZikrLr9H6yYeaep4Wtw9Jz9UjOzkOWzvQ9PVdfceuwClbZqvUvN1xZi4hIctUuCRBFsQDA+preUBTFp2t6LlGhgaHu+CU6AXqDsdjIWX018ZeqBEFyEvQ8rXbrsErqS2tdZlCHXRCIiKhmalLDSiQ5czfxl6IEQTYk6Hla5dZhVagvHRjqUbtfbuqwCwIREdUME1ZqsMzZxL8+JgM1KHLteVqF+tKAnv+s3S83XFmLiEhyTFipQTNXE39zlyBQFSWcBXKSgexEwMIKsHEHLKxN+4rUl1b4y01V+qtyZS0iIkkJoij/jlGdOnUSjx8/LnUY1MRti7pV7igdFyqQwOlNwN53AF0GoLQERCMAEXD0A+w8gewkIGQI0HNWxdcoWv9akAfAyP6q1JCxgw81ShxhJaoic5Yg1LXy+sc22L6yhbWrDs2BxPOAQgEIKkAUgbQbgNoaldaXsr8qEVGDUW7CKgjCNFEUV5gzGCK5M1cJQl0qb5Wntt52OHsrs9T2BrG0bWHtqpWDKWlNv2n6LCgAY4Hp86MLKk44q7qKFxERSa6iPqyfCoKwXxCEhvW3MxE9UF7/WAuVgIioeFiohIbZV7Zob1R7L8AjDLBxAyxsAK0L4Nu18lf67K9KRNRgVJSw9gPgAeC0IAizubIVkbzEJGVj5e9X8UbEGaz8/WqZSWZ5/WPTckwtmUzf/1bl1Z+kVnIFLgsbwDkAcG8DWLsC7m2rf42i2F+ViEhWyk1YRVH8A0BbAJ8BeB/AIUEQWpspLiKqQFWXiS2vf2yu3gClQgFdvqHUvgbRV7YuVuAy8ypeRERUcxVOuhJFUQdgriAImwCsBhAlCMJ2APmlDxWfq58Qiaiooq/5C0dObWFqsbX24A2E+Tg8mDhVXv9YK7USaffyobFQlrp+g+grWxe9UdlflYiowahql4ArAE4BaA+gJ8pIWOsuJCKqSHWWiS2vf6yjVo34dB0ctcUnHDWovrJ10RuV/VWJiBqEShNWQRCGAvgcgDWAF0RR/KbeoyKiclVnmdiKlrAdEe6Js7czkZydV+9L29abuliBS66reBER0QMVtbVyhal+dRSA/wJ4SRTFeHMFRkRlq+4ysRX1jy3sw9rQ+spWaXUqIiJqNMpd6UoQhBSYXvXPEEVxo1mjKoErXRH9LSYpG6/9eKpYDStgep2fnqvHJ6PbN4yks1B1k0+uTkVUEXb0oUapopKA3wBMFUUx0VzBEDVU5lwxqqLX/A3qdT5QOvlMOAtc2FV+8snVqYiImqRyE1ZRFJ80ZyBEDVV5K0nV54pRDXmZ2AdqknxydSoioiapql0CiKgM1WkxVZf3LDqa+3yPFg0rUS1Uk+STq1MRETVJFa10RUSVqEqLqbpU1QUDGoSaJJ9cnYqIqEliwkpUC9VpMVVbRUdzXW0tYatRw8XGEg5aNdYevFHm0qyyVpPkk6tTERE1SUxYiWrBx1GLvAJjmfvqesUoc4/m1ruaJJ+Fq1PlpgHZSYAu0/Q9N42rUxERNWKsYSWqhfJWkqqPFaPMOZprFjVdGpWrUxERNTlMWIlqwZwtpqq7YECDUNPkk6tTERE1KUxYiWqpPltMFe0IoLVQITffYJbRXLNi8klERJVgwkpUBwJcbfDKIy3r9Jol+7vmFRiRqzfgevI9OFlbNOwFAypgzkUYiIioYWDCSiRD5fV3tbdS426mDp38nJCrL2iYCwZUQIpFGIiISP6YsBLJUEUdATRqJZytLfDKI62lCK3eSLEIAxERNQxsa0UkQ42uI0AVNLq2XUREVGeYsBLJkDn7u8pFU0zSiYioapiwEsnQwFB3iKLpdXhRDb4jQAWaYpJORERVwxpWIhmqqL/rY209GuUs+rpchIGdBoiIGhdBFEWpY6hUp06dxOPHj0sdBpHZFSZehf1dlQLw81lTLWdhqytRRKOZRb8t6la5izBU9ecrqx1YY3pGRJUQpA6AqD4wYSVqIGKSsvHaj6eKzaIHTCOQ6bl6fDK6faMYRSyZpFdndLSpPCOiCjBhpUaJJQFEDURVZtHX9eIFUqjNIgzmfEYsOyAiMh8mrEQNhDln0TfUZMxcz4gLHBARmRcTVqIGwsdRi/PxmbAtY19dzqJvyMmYOZ4RFzggIjI/trUiaiDM0eqqaDLmamsJW40aLjaWcNCqsfbgDcQkZdf6HvXJHM+ICxwQEZkfR1iJZK7463krXEnKhpVaWWoWfV2M6jX0OtmK2oHV1TPiAgdERObHhJVIxspq0QQA3o5WsNeqqz2LvjKNIRkbHu6NMB+HGncaqIy5SjOIiOhvTFiJZKqiWsm4tFzMGhBc57WSjSUZq02ngcrU5QIHRERUNaxhJZIpKWolm+KSsNVVWHaQnqtHcnYesnSm7+m5+jorOyAiouI4wkokgaq0jZLi9bw5akAbg/ouOyAiouKYsBKZWVXbRkn1ep7JWNXUZ9kBEREVx4SVyIyq08NTylpJJmNERCQnrGElMqPq1KU2uFrJ5CvAX0uBnTNM35OvSB0RERE1EhxhJTKj6talNpjX86c3AUe+ACAAKg2QcBa4sAvo8jLQbrTU0RERUQPHhJXIjGpSlyr71/PJV0zJqpUToFT/vd2gB458CXiFAy6B0sVHREQNHksCiMyoUbaNurATgFA8WQXufxbu7yciIqo5JqxEZtTg6lKrIi3WVAZQFpUlkB5n3niIiKjRYUkAkZk1mLrUqnL0NdWslqUgD3DwNW88RETU6DBhJZKA7OtSqyNkqGmClUFfuoYVIhDyuGShERFR48CSACKqHZdAUzeA3DQgOwnQZZq+56YBXaZwwhUREdUaR1iJqPbajTZ1A7iw01Sz6uBjGnllskpERHWACSsR1Q2XQKDnP6WOgoiIGiEmrERylXzFNGKZFmua2NTERixjL5/C7UNboMi4CaN9c3g9/CR8g9pLHRYREUlAEEVR6hgq1alTJ/H48eNSh0FkPiVXjirIA2Cs25WjZJwQH9vxBezPfAMAKFBYQCnqIYgiMsKex0NPvCxxdESyJkgdAFF94AgrkdyYY+UoGS+lGnv5FOzPfAOdyg5Gxd8/v8Koh/3ZbxAb8jBHWomImhh2CSCSm/peOapoQmzjBmjsABtX0+cjX5r2S+j2oS0AUCxZLfwsQniwn4iImg4mrERyU98rR8l8KVVFxk0UKCzK3GcQ1FBkxJo5IiIikpokCasgCIMEQbgkCMJVQRDmSREDkWw5+t6vWS1DXawcJfOlVI32zaEU9cW2GYwicvMNMOjzcMPggpikbImiIyIiKZg9YRUEQQngcwCDAbQG8LQgCK3NHQeRbIUMBWC8v1JUEXW1clR9J8S15PXwkxBEEQqj6efPKzAgS6eHoSAPEEX8z9AJr/14CtuibkkaJxERmY8UI6ydAVwVRfGaKIr5AH4E8IQEcRDJU32vHFXfCXEt+Qa1R0bY89AYMqHJS4EiPwtOyIC9cA+7rYcBLkFw0Kqx9uANjrQSETURUnQJ8AJQ9J3jLQBdSh4kCMKLAF4EAF9faUd8iMyuPleOKkyIj3wJU5cAy/sjrqJsllJ96ImXERvyMA7/tA6KzFtIs3DHeYdeyNT6AwDUSgUEAHuiE/DKIy2lDZaIiOqdbNtaiaK4CsAqwNSHVeJwiMyvPleOagBLqfoGtccXfpMQk5gNW4261H5LlQK30nMkiIyIiMxNioT1NgCfIp+9728jInNqAEup+jhqcT4+E7Zl7MsrMMLbQWv2mIiIyPykSFiPAQgUBMEfpkR1DICxEsRBRHKWfAWj87aiRdYpZOZ74JRtL9xVm37X1RuMEAEMDHWXNkYiIjILsyesoigWCILwKoA9AJQA1oiieM7ccRCRjN1ficsZAjpaC0jNvIb29w7gF5th2KfuAxHAxG5+CHC1kTpSIiIyA0lqWEVR/C+A/0pxbyKSuRJL07raAFo7R6Sk38PwvF1wCeqKLp27MlklImpCZDvpioiaqDJW4rK2UMHazR7IzsdYu9OA66PSxUdERGbHpVmJSF5kvhIXERGZHxNWIpIXma/ERURE5seElYjkReYrcRERkfkxYSUieanvpWmJiKjB4aQrIpKfBrASFxERmQ8TViKSpwawEhcREZkHE1aixij5iml0Mi3WNImJo5NERNSAMWElamzurxIFCKb2UAlngQu7THWh7UZLHR0REVG1MWElakxKrBL1gEEPHPnSVBfKkVYiImpg2CWAqDEpY5UoAPc/C/f3ExERNSxMWIkaE64SRUREjRATVqLGhKtEERFRI8SElagx4SpRRETUCDFhJWpMuEoUERE1QuwSQNTYcJUoIiJqZJiwEjVG9bRKVExSNvZEJyAuLQc+jloMDHVHgKtNnd+HiIioKCasRFQl26JuYW3kDQCARq3A+fhM/BKdgInd/TA83Fva4IiIqFFjwkpElYpJysbayBtw0KqhVppK320B6A1GrD14A2E+DhxpJSKiesNJV0RUqT3RCQDwIFktpFYqIBTZT0REVB+YsBJRpeLScqBRl/3HhaVKgVvpOWaOiIiImhImrERUKR9HLfIKjGXuyyswwttBa+aIiIioKWHCSkSVGhjqDlE01awWpTcYId7fT0REVF+YsBJRpQJcbTCxux/Sc/VIzs5Dls70PT1Xj4nd/DjhioiI6hW7BBBRlQwP90aYjwP2RCfgVnoOvB3Yh5WIiMyDCSsRVVmAqw1eeaSl1GEQEVETw5IAIiIiIpI1JqxEREREJGtMWImIiIhI1piwEhEREZGsMWElIiIiIlljwkpEREREssaElYiIiIhkjQkrEREREckaE1YiIiIikjUmrEREREQka0xYiYiIiEjWmLASERERkawxYSUiIiIiWWPCSkRERESyxoSViIiIiGRNJXUARCR/MUnZ2BOdgLi0HPg4ajEw1B0BrjZSh0VERE0EE1YiqtC2qFtYG3kDAKBRK3A+PhO/RCdgYnc/DA/3ljY4IiJqEpiwElG5YpKysTbyBhy0aqiVpgoiWwB6gxFrD95AmI8DR1qJiKjesYaViMq1JzoBAB4kq4XUSgWEIvuJiIjqExNWIipXXFoONOqy/5iwVClwKz3HzBEREVFTxISViMrl46hFXoGxzH15BUZ4O2jNHBERETVFTFiJqFwDQ90hiqaa1aL0BiPE+/uJiIjqGxNWIipXgKsNJnb3Q3quHsnZecjSmb6n5+oxsZsfJ1wREZFZsEsAEVVoeLg3wnwcsCc6AbfSc+DtwD6sRERkXkxYiahSAa42eOWRllKHQURETRRLAoiIiIhI1piwEhEREZGsMWElIiIiIlljwkpEREREssaElYiIiIhkjQkrEREREckaE1YiIiIikjUmrEREREQkaw0iYb1x44bUIRARERGRRBpEwurn5yd1CEREREQkkQaRsBIRERFR08WElYiIiIhkjQkrEREREckaE1YiIiIikjUmrEREREQka0xYiYiIiEjWmLASERERkawxYSUiIiIiWWPCSkRERESyxoSViIiIiGSNCSsRERERyRoTViIiIiKSNSasRERERCRrTFiJiIiISNaYsBLVI1EUodMbkJ6TD6NRlDocKo/RCOSkAgV5gMh/T0REcqOSOgCixqDAYER0fCYirybjUEwKriRmISffgJx8Awz3E1UrtRIt3WwQ7G6LoGY2aOvlgC7+TlAoBImjb2IMBcDNA0D8KSDpIpB4Hki6DBTkmvYrVICFNaC2BtxaAS36AP69APcwQKGUMnIioiZLEBvAaEKnTp3E48ePSx0GUTGiKOLQtRR8e/AGDsakIEtXAADwddLC10kLK7USGrUClmoF1EoFkrPycDtdh9vpucjI1QMA3O00eLKTN57s6ANfZ62UP07jl3wVOLUBOPUDkJ1g2qZ1Buy9TV9aV8CQDxToTF/6HNM5mbdMx2rsgRaPAF1fBny7SvdzEFWMvwFTo8SElaiajEYR/zufgC/+vIbTcelw0KrRztsBbTxs0drTHs3sNNCoyx+JMxhFpGTnISo2DX9eTsK5+EyIADr7O+Hl3gF4pJWb+X6Yxk4UgYs/AQc/A+IOA4IS8AoH/HsDHu0BaydAaVn+yKkoAlnxQNxRIOEscOsYkJcFeHcGes0GAgcAAvMDkhX+B0mNEhNWoioSRRE/n72DZb9exrWke3C302BAm2boE+QGFxsLCDVMXOLTc7HnfAL+upyMpOw89Ap0wYKhbdDC1aaOf4ImJvECsHsOcH2/aQTVv4/p9b6jL6C0qNk18+8B0duA89uBnGTAtRXQ7x2g1T/qLm6i2mHCSo0SE1aiKkjKysO/tp/FnnN34eesxaBQd/QIcIWdlarGiWpJ+QUG/CfqNiJO3kKBQcTE7n6Y8WgQbCxZal4tuWnAHx8AR7821aK2fRJoNQTQOtbdPQx64OJ/gbObgczbQJvhwGPLAK1T3d2DqGaYsFKjxISVqAKiKGLXmTt4Z0c07uUZMKyDJ4aGecFeq663eyZn52FN5HX8dSUZrjaWWDG2A7q2cK63+zUqV/cCES8BualA4EAgbLRpRFWop4YoBj1wfA1wZjNg5QA8/m8gZEj93IuoapiwUqPEhJWoHBk5esz9zxn8ci4BgW42mNC9Odp4OEBppln9525n4JN9V5CUpcPrA1thSu8WdTaa2+gYjcD+j4E/3gec/IEuUwDP9jV/9V9dSReB398H0m8CbUYAQ/8NWNqa595ExfEPCWqUmLASleFmyj1MXHcMsSk5GN7BC0+0q99R1fLk5Bfg4z2XcPxmGh4NccMno9vDVmP+OGQtJxWIeBG4+qtpFn+XKYBtM/PHYdADx74xlQm4BAHj/2OqnSUyLyas1CgxYSUq4diNVLz43XEYjCJe7hOArv7OUCmlW2NDFEVsPhGH74/EwsvBCt9MeAhBzTh6B8DUS3XTeFObqvAJQJthprpVKd08COx7D7C0AcZuAbw6SBsPNTVMWKlR4kpXREXsOHUbY78+DK2FCm/8oxW6BbhImqwCgCAIGN3JFwuHtkFGrh4jVx7EyZtpksYkCzcOAOseM/VO7fcu0G609MkqADTvBgxdYfrntQOB8zuljYeIqBFgwkp032e/XcGMH08hqJkt5v+jFUI9HaCQUc1oOx9H/N+T7aBRKzF29REcikmROiTpXNkLbBhpavzffyHQ/GHTClVy4dISGPYlYO8LbH4WOLhC6oiIiBo0JqxEAP697wr+73+X0aOlC+YMCIavkwxG6srgYW+Fj0eFwUGrxoS1R7Hvwl2pQzK/8zuAH8YA9j6mHqiuwfJs3m/tDDzxmWlVrP/9C4hcLnVEREQNFhNWavJW/nEVy369jF6BLnilTwCcbCylDqlCzjaW+HhUOzSz0+Cl9Sew63S81CGZz6kfgC0TAJdAoN/bgHMLqSOqmMoSGPAe4NsN+PVt4ODnUkdERNQgMWGlJm3V/hh89Msl9GjpjJd6BTSYGfj2Vmp8ODIMfi7WmPHjSfz37B2pQ6p/Z7YA26cA7mHAI/8CHHykjqhqFCpgwELApyvwv/nAkVVSR0RE1OAwYaUm65sD17HkvxfxcAtnTOnVEnZWDSNZLWRjqcKSYW3h52yNmT+ewuGYZKlDqj9X991PVtsBj7wJ2HtKHVH1KFSmkVbvh4DdrwPH1kgdERFRg8KElZqkzcfj8N5P59HZ3wkv92khSY/VumBlocTCJ0LhZGOB5787jgvxmVKHVPdunwA2PQM4NAd6zQZsXKWOqGaUamDgYsArHPh5lml1LCIiqhImrNTkHLyajPkRZxHmbY+pfQLgqJV3zWpl7K3UWPxEKNRKBZ5ZcwS30nKkDqnuJF8FNj4JaOyB3vMaThlAeZQWwMD3AbcQYMcrwM1DUkdERNQgMGGlJuVqYhZe2nACng5WpglW1g07WS3kZqfBwqFtcC/PgLFfH0HavXypQ6q9rARgw3BANAJ95plaRTUGKktg0BLAytnU7SD1utQRERHJHhNWajKSsvIwYe0xqBQCpvVtCS8HrdQh1Sl/Fxu8NSQE8em5eG7NUeQXGKUOqebyc0wjq/eSgF5zAPe28mxdVVMaB+AfHwNGPbB+GJDLhSCIiCrChJWaBJ3egBe+O46krDy82rclghvp0qZtvRwws18gztzOwBsRZ6QOp2ZEEdg5DUg4C3SbAfh0blzJaiEHH2DAIiAjDvh+NFDQCEbFiYjqCRNWavREUcQ/t5zG6bh0TO7hj46+ThAaYwJ0X+9gNwzv4IX/RN3G2sgG+Lr54L+B6K1Au7FAy76AQil1RPXHswPQczYQdwTYMdWUrBMRUSlMWKnRW/3Xdfx85g6e7OSNfiHNoFQ03mS10HMP+6G9jwMW/XwBhxpSu6ure4G9C4DmPYD2Y02TlBq74MGm5PzsZuAoe7QSEZWFCSs1akeupeCDXy6is78TRnTwglrZNP6TVyoEzB3UCi42Fnh5QxTi0xtA54CUGGDrJMDRH3h4KmApz+Vx68VDzwOe4cCeN4G4o1JHQ0QkO03jb29qkhIzdZj6fRQ87DSY3N0f1pYNs9dqTdlYqvDOkDbQFRgwce1x6PQGqUMqX14W8OM4QATQYxZg5yF1ROalUAKPvgNY2QObxgP3UqSOiIhIVpiwUqOkNxgx9fsoZOcV4KXeLeDhYCV1SJLwcdLitUeDcOluFt7cdlbqcMomisCumUDyJaDbdKBZiNQRSUNjD/R/D8hJATY/Axhl/AsGEZGZMWGlRumD3Rdx7EYanu3aHG29HKQOR1LdAlzwRDtP/CfqNrafvC11OKWdXG+aZNX2KaBFL0Bown8suYUA3V4FbkYCe9+VOhoiItlown8zUGO1++wdfHPgOga0boYBbdybxCSrykzo5oeWbjaYv+0sbiTfkzqcvyVeAP47B/BoD7QbY1q+tKkLeQJo+ShwcDlwabfU0RARyQITVmpUbqfnYu5/ziDQzQbjOvvCUtWIWyJVg0qpwLxBrQAAL60/IY9FBfJzgC0TALXGNMnKykHqiORBEIBeswEHX2D7y0DWXakjIiKSHBNWajQMRhGv/XgKeoOIF3r6w8mmcSy7Wlea2WkwvW8gLt3Nwru7zkkdDrB7DpB0CegyFXBuJMuu1hWVBnh0AZCXDWydCBhl8AsGEZGEmLBSo/H571dx9EYqxnXxRVAzO6nDkaXuLV0wONQdG4/E4r9n7kgXyJnNptrV0BFAQK/GuZJVbTm1ALpOMdWzRi6XOhoiIkkxYaVG4cTNVCzfdwXdWzpjIOtWKzS5Rws0d9ZibsQZJKTnmj+AtBvAT68BzdoAHZ5pGosD1FSbEYBPF+D3RcDtk1JHQ0QkGSas1OBl6vSY/sMpuNpY4tmuzaFRs261IhYqBeYObAWd3oCZm09BNOdyoEYDsG2KqZVVV9atVkoQgEfeACxtgS3Pmep+iYiaICas1KCJooh/bYtGQoYOz/fwg4d90+y3Wl0+TlpM6OaPw9dS8c2B6+a78cF/A7GHgI4TANdg8923IdM4AI+8CaTHAj/NlDoaIiJJMGGlBm3XmTvYeToeQ9t7omNzJwishayyIWEeCPO2x0d7LuHK3az6v+GdM8Bvi4Hm3YFW/zCt7kRV490JaDsKOLMJOL9T6miIiMyOCSs1WHczdXhrezSCmtlgZAdvqJX8z7k6FIKAWY8GQa0Q8Or3J6E31ONMdL0OiHjRtJpT5xcBC+v6u1dj1fkFU6urXTO4dCsRNTn8G54aJFEUMfc/Z6DTGzCxmz/stWw4XxPONpaY+khLXLqbhf/bc6n+brRvIZB0AXjoRVPSRdWntAD6vgXkZQA7pkodDRGRWTFhpQbpx2Nx+ONSEkZ19EaIB1tY1UbPQFf0CnTB139dw7HrqXV/g2t/Aoc/B4L/AQT0Zgur2nAJBDo8C1zeDZz6XupoiIjMhgkrNThxqTlY9NN5hHra4bG2HmxhVQde6dMSDloLzNp8Cjq9oe4unJdlGg209wbCnwNUXMyh1jqMA1yCTQsvZMRLHQ0RkVkwYaUGxWgU8c8tpwEAE7r5wVbDUoC6YG2pwoy+LRGXlovFP5+vuwv/7y0g4xbQ+SXAtlndXbcpU6iAvm8CBXnAthdNLcKIiBo5JqzUoKyJvI6j11Mx+iEftHSzlTqcRiW8uRP6hzTDhiOxOHwtufYXjPkNOLEWaD0U8O1S++vR3xx8gYdeAG78BRxdJXU0RET1jgkrNRjXk+/h4z2X0NHXEQNbczWr+jC5pz+ctBb45+YztSsN0GUCO6YB9r5Au7Fczao+tB0JuLcF9i4w9WglImrEmLBSg2A0ipiz9TTUSgXGP+wLraVK6pAaJa2FCjP6BeJ2ei4W/VSL0oD//QvIige6vMhSgPoiKIA+b9xfPexllgYQUaPGhJUahG8P3cCxG2kY3ckbLVxspA6nUevg64gBrZth45FYHIypQWnA1b1A1LdAyFDAh6UA9crO09Sf9eYB4PgaqaMhIqo3TFhJ9m6m3MNHv1xCBx8H9G/dDAq2Rap3z/fwh7ONBV7fUs3SAF0msHO6qcay/VhAyUlx9S50BNCsDfDrW0D6LamjISKqF0xYSdZMpQBnoBCA8Q/7wtqSCZA5FC0NqFbXgF/fBrLumFazsnGrvwDpb4WlAYZ8YDtLA4iocWLCSrK28chNHLmeiic7eaOlK7sCmFN7H0c8er9rwPEbVVhQ4Pp+U1eAkMdZCmBu9t5Ap+eBG/uBqO+kjoaIqM4xYSXZikvNwfu7L6Kdtz0GtXFnKYAEJvfwh4OVGrO3nEZ+gbH8A/PvATunAXbe97sCcCTc7No+Cbi2Ava8CWTekToaIqI6xYSVZEkURczfdhaiCIzv2pylABKxtlTh1Uda4kZKDv5vz8XyD/xtMZB2A3hoMrsCSEWhBB55AyjQAbumSR0NEVGdYsJKsrT1xC38dSUZIzt6IZALBEiqs78zega64JsDN3D2dnrpA+KOAYdXAkGDgeYPmz0+KsKhORD+DHDlV+DMFqmjISKqM0xYSXYSs3R476fzaOVui3+EenCBABmY0isA1pZKzNp0GgWGIqUBBXnAjqmmCVYdxgMqS+mCJJP2YwFHf2D360BOmtTREBHVCSasJDtvbz8Hnd6IZx9uDlsNSwHkwM5KjZd7B+BKYjZW/Hb17x37PwaSL5km/Nh5Shcg/U2hAh6ZD+gygJ9ekzoaIqI6wYSVZGX32Tv45VwChrb3RGsPe6nDoSJ6BLqii78TVv5xFVfuZgEJZ4EDnwABfYEWvQFOipMPl0Ag7Cng/Dbg0m6poyEiqjUmrCQb6Tn5eGtHNAJcrfFEO0+WAsjQK31aQq1U4PXNURB3TAUs7YDw51gKIEcdJwJ2XsCumUBeltTREBHVChNWko33frqAtBw9nu3qBwethdThUBmcrC3wfA9/dEn4AcKd06akyMFX6rCoLCpLoM88IPsusHue1NEQEdUKE1aShT8uJeI/Ubfwj1B3hHmzFEDOHvO8h1nq/+B/xocQ16wvSwHkzL0t0HoocGojcO1PqaMhIqoxJqwkuey8AszfdhbejlYY3sEbKiX/s5Qt0QjvA/OgUFlgoeE5zP0jByKXApW3Li8B1i7AzlcBvU7qaIiIaoSZAUnu418u4k66Ds8+3ByutqyFlDOnCxthk3AE6a1G4/G27jh4W48t5+5JHRZVRK0Fer0OpMcCexdIHQ0RUY0wYSVJHbuRim8P3UT/1s3QqbmT1OFQBdTZt+F+7H3kuIQhw3cAngiyRCsnBRYdyERitkHq8KgiPp2BwP7A0a+AW8eljoaIqNqYsJJkdHoD5m49g2Z2lniykzfULAWQL1GE14E3IBgLkBLyLEQLaygVAmZ3sUJOATD/Nzaol71u001dHba/DBj0UkdDRFQtzBBIMv/edwXXku9hfFdfNLPVSB0OVcDh6n9ge+sPpAWOhM4p8MH25vZKjGttgb038rHrUo6EEVKlLG2BHq8ByZeBPz+UOhoiomphwkqSiL6dga/2X0PvIFc83MIFAmeay5YqJxEeh95FrmMrpLcYAgjKYvvHtLaEv72Ad/7MQHqusZyrkCy06A00725a8CHhnNTREBFVGRNWMrv8AiNmbzkNeys1xjzkA0uVsvKTSDKeB9+CoiAXKa0nwGhhW2q/SiFgTlct0vNEvPU7SwNkr+csQKUBtr0EGFl7TEQNAxNWMrsv/ojBxYQsjOviCy8HK6nDoQrYXf8Z9jd2Iz1gGHKdQ8o9rqWjEk8FW2DX1Tz8ei3XjBFStWmdgW7TgLv3l9YlImoAmLCSWV1MyMRnv19BtwBn9Ap0ZSmAjCl1afCMfAs6+wCkBQwFFBWPhD/T1hJeNgLm/5aOzDyWBsha4ABT54A/PwKSLksdDRFRpZiwktkUGIx4fcsZWFuqMLazDzRqlgLImcehBVDlpSGl9XMwahwqPd5CKWBOVysk54hY+Gd6vcdHtSAIpt6sCuX90gD+gkFE8saElcxm1V/XcPZ2Bp5+yAe+TtZSh0MVsL35PzjGbEN6i6HIcQ2r8nmtXVQYFqjG1os6/HGDpQGyZu0KPDwViI8CDn0mdTRERBViwkpmcTUxC5/uvYLO/o7o26oZSwFkTKlLh9eBN5Bn54/UwOGAQlWt8ye108DTWsDcvenIYmmAvAX/A/AMB35fDKRelzoaIqJyMWGlemcwinh96xloVAqM79KcpQAy53HoHah0qUhuPQFGTfVXH9OoBLze1QqJOSLeZWmAvAkC0Geu6Z8jXmRpABHJFhNWqner9l/Dydh0jOnsi+bOLAWQs2KlAG7tanydUFeWBjQYNs2ALlOAW0eBw59LHQ0RUZmYsFK9upSQhWW/XkJnf0c82soNCpYCyJZSl1arUoCSWBrQgLR+wlQa8Nt7QPJVqaMhIiqFCSvVG73BiH9uOQVrCxXGd20OK4vaJUBUvzwOLahVKUBJLA1oQAQB6DPPtIpZxGSWBhCR7DBhpXqz8vcYRN/OxLiuvvBjVwBZs7u+u05KAUoqWhqwjwsKyJuNm2lBgfiTwIFlUkdDRFQME1aqF9G3M7DiN9MCAX2C3NgVQMZUOYnwOjAPOvsApAaOqHUpQEnPt9PA20bA3H3pSM/lyJ2sBQ0CfLoAf34I3L0gdTRERA8wYaU6l1dgwOwtp2Fnpca4Lr7sCiBnogivA/Og0N9DcujzMGoc6/wWlioB87ppkaYT8ca+1Dq/PtUhQQB6zwGUFqbSAEOB1BEREQFgwkr1YNmvl3ExIQvju/rCx1ErdThUAcfLm2AXuxdpgaOQ69ym3u4T7KTE2NYW2H0tH9sv5tTbfagOaJ2BHjOBu9HAH+9LHQ0REQAmrFTHDl9Lwar91/BIK1f0CnRlKYCMqbNi4XHoXeQ4hyI94HHTMp31aGwbSwQ6KvDWH+lIyDLU672ollo+CrToAxz4BIg7KnU0RERMWKnuZOTqMWvzKXjYaTC2sy8sVSwFkC2jAT5/zgIAJLeZBKOFbb3fUqUQ8MbDVtAVAK/9LxWiKNb7PakWev4TsHIAtj4P5N+TOhoiauKYsFKdeWdHNO5m5GFiDz+421lJHQ5VwCX6a1gnHEVq8NPIc2xptvv62CnxQjsNDt3WY93pbLPdl2rA0hZ45E0gIw74ebbU0RBRE8eElerEztPx2H4qHkPbe6Cjb+17eFL90SSfQbPjHyPbvTMy/AYAgnn/GHgiSI32bgp8EJmFS8l6s96bqskrHGg7Ejj9PXBhl9TREFETxoSVai0+PRf/2nYWQc1sMLKDD9RK/mclVwr9Pfj+Ng0GC3skh0yAqDb/pDiFIGDew1pYKoGpu1ORV8DSAFnr/CLg6AfsnAZkJ0kdDRE1UcwsqFYMRhH/3HwaeoMRk7r7w16rljokqoDHoQWwyLyB5DbPQ2/nI1kczlYKzO5ihatpBizcny5ZHFQFSgug39tAXrap1RVrj4lIAkxYqVZW/n4Vh66lYExnX7Ryt5M6HKqA/bVdcLq8CektHsc9z66mnpsSethLjccD1NgYnYv/xXAVLFlzamEaab32BxC5XOpoiKgJYsJKNXb0eio+2XsZ3QKcMbiNO5QKtrCSK3XWLdNqVo5BSA16EqJSHiPhU8I18LUTMGdvOhLvsdWVrLUdBfh0BX57D7h1XOpoiKiJYcJKNZJ2Lx8zfjwJdzsNJnTzg5VF3S7nSXXIWACfP2YARgOS2rxQL6tZ1ZSFUsBb3bW4pxcxfXcqjHzdLF+CADzyBqCxB7Y8B+RlSR0RETUhTFip2kRRxOwtp5GUlYfJPVvAw54trOSs2YllsL57DCmtxkPnHCR1OKX42SvxUnsNDsfrseJoptThUEU09kDft4HM28C2KaxnJSKzYcJK1bY28gb2XUzEk5280cHHQepwqAK2cb/B7fRnyPTug0y/AYAgz8Uchgaq0cNLheVH7yEyTid1OFQRz3ZAh+eAiz8Bx9dIHQ0RNRFMWKlaztxKx/u7L6CjryOeaOcJFVtYyZY6+za8/5iJPDs/pLSeAFGlkTqkcgmCgNldrdBMK2D67jTWs8pd+DOAexiw5w0gIVrqaIioCWC2QVWWei8fUzacgKPWApO6+8HaUh4Td6g0wZAP332vQDDkI7HtyyjQukgdUqWs1QIW9NQiM1/EK/9NgcHI182ypVACj74DqKyAH58GdCzlIKL6xYSVqsRgFDHjx5NIysrDS71bwMfJ/A3nqercjy6GNukkkttMhM6lldThVFkLByWmddTg+J0CfHgwQ+pwqCJaZ1PSmnEL2DqJ9axEVK+YsFKVLPv1Ev66koxxXZqjo68TBIl7eFL57K7/DJdza5HRfCCyfPvKtm61PIMDLNDfT4VVUTn49Rr7s8qaZwdTf9arvwJ/fih1NETUiDFhpUr971wCPv89Bn2CXfFYWw/2W5Uxy9RL8P5zNnSOwUgOGQdRaSl1SDUyvZMV/OwEzPpfGq6lFkgdDlUkbDTg18uUsF7ZK3U0RNRIMWGlCl1PvodZm0+jpZsNnu3qB426YY3WNSVKXRr8fp0EUaVBYtjLMGqcpA6pxjQqAQt7WUMUgUm7kpGVZ5Q6JCqPIACPzAPsvID/TALSYqWOiIgaISasVK4snR4vrT8OhQC81KsFXG0b5mhdk2AsgO9vr0B1LwGJYS8jzyFA6ohqzcNGgbe6axGbacSrXFRA3tRaYOBioCAf+GE0oGcpBxHVLSasVCaDUcT0H04iJvEeXujVAkHNbKUOiSrgcWQRbOIjkRLyLO65P2Qa9WoEwt1VeLGdJf6MzceHkZyEJWsOvsAj84HEC6ZJWEaOihNR3WHCSmVa/PMF/H4pCeO6+KJHgAsUjSQBaowcL2+Gy7k1yGg+CBn+g00thxqREcGmSVhfReVg+8UcqcOhivj3BB56Hrj0X2Dfu1JHQ0SNCBNWKuX7I7FYE3kdA9s0w+NcHEDWtAnH4HlgPnJcwpAc8oysFweoKUEQMPMhKwQ7KTB3XzpOJ+RLHRJVpP04ILA/EPkpcHKj1NEQUSPBTISKOXg1GW/viEYHHwdOspI5i/QYNP/1eRRoXZHY9mUYNfZSh1RvLJQCFvbUwtZCwMSdKYjLYOcA2RIEoNccoFkb4KcZwM1DUkdERI0AE1Z6ICYpG1M2nICngxVe7h0AOyuuZCVXqpwk+O95FgBwt/1M6O19JI6o/jlZKfBBby3yDCKe2Z6MDB1rJGVLqQYGLgGsnE0rYaXdkDoiImrgmLASACAhQ4dnvzkKhSBgWt+W8HCwkjokKoegz0Hz/02EKicJd9tPg8654axkVVvNHZR4t4cWt7KMmLAjGXkF7BwgWxp74B8fmjoHfDsUuJcsdURE1IAxYSVk5Ojx3JqjSL2Xj+n9WiKYHQHk6377KquUaCSGvYQc906NpiNAVbVrpsLrna1w8m4Bpv/Cdley5tAcGLQEyIoHvhsK5GVJHRERNVBMWJu43HwDnv/2GK4lZ2PqIwEI57Kr8iWK8Ip8E3ZxvyElZDyyfPo0uGVX60pfPzUmtbXAnmt5ePfPdIhMWuXLox3Q921Tu6uNT5lGXImIqokJaxOmNxjx6vdROHEzDS/08EePlq5cdlWuRBHuRxbB6dIPSGsxFOkthgKKpl1jPKa1JYYEqPHtmVwsPZwpdThUEf+eQI9ZQOxBYOtE9mglomoTGsLIRKdOncTjx49LHUaZ3t11DgDwzuNtJI6keoxGEa9vPYP/RN3Cs12bY1gHL6graF/19V/XAAAv9GxRr3FVdJ/y9pkrtqL3i76dgVCvv2flv9CzRblxeBxaAABYWPBssf1Fr1Ne7B6HFiA6PgOLMwahR95+zGx5F6KgxL44YKX6OYS5qQAAr4SXbmflcmYVACA57EUAwMooXaljXc6swpkkA7Y7TsQr4Zpi55Q8Pn7fFwCA7Y4Ty72nz2/TAQBxff9d5s9TnrJi/e1GPn4U3oCthYDcxz7DyigdhqWtRZirEslhL8LlzCrY3voDOXoRW/Xd0dHDAu2aqYFu0wAA7+7PwMC7q9HVy6LYNgB4p1eJjgoHV5i+3z+u3G3/mWz6PnJ1sesVekf1XfFzih5f1vVqoq6uU96175w2jYrWx/Wj1gPHvwHCnwMeX179cpbd80zfB39Q97E1Hhx1oEZJJXUADd35+IY3smM0ipi/7Sz+E3ULwzt4Ymh7zwqTVQC4lpRtltgquk95+8wVW9H7JWTqoLVQltpeFk3KedP+/OL7y7tOyXNd0lORkDsIZzXtkNzGGx5HFsGnwICEfCOs0wzlnmuZca3Y55gyjrXMuAanHCNiYCh1TsnjnXKum7aj/HuqcxLK3VeRsmLN0gNeFklQ6IHc+9uccq7DMkPx4BxlfjZsATysiUXGXRF3dEp4dDNd43ySHqOyYoCUv0eizyfpyw4g5WrVtmXdKfax1PUsSpxT9PiyrlcTdXWd8q6ddQewsK6f63cYD+jSgahvAaUF8I+Pq5e0Jpytn7iISPZYEtDEFCarPx6Lw9B2Hnj6IV9YqppmHWRDoM5JgGDQQRRU0Nt4wmhhI3VIsuRrr4CtWsDNDAN+iDbvLzBUDYIAPPwqEDIUOPY18N/XgQbwlo+IpMcR1ibEaBTx5va/k9XxXZrDyoL/CciV0/n1sMy8AQgtICotIDbRCVZVIQDwc1DgWroB83/LhJITB+VLEIAer5n++djXps+DP2py3S6IqHqYrTQRhcnqD0eZrDYELqe/gMex92FQ28BosACM/Mu8MgIAf3sl2moVmLMvA83tmeDLVtGk9aipfplJKxFVhCUBTYDeYMTrW88wWW0IRBHNjn8Mj2PvI8ujG/JtvAGB/5tWlUIA3u9jjYfclbiZYUC+ga+bZaswaQ0Zakpad80AjOXXRxNR08a/CRu53HwDXlp/4sEEKyarMiYa4XF4AdxOrUCm9yNIbP8qRAX/XVWXhVLAwl5a2FsAeQbgZoaBfVrlqjBpbTPSNBFr8zNAQZ7UURGRDDFhbcTS7uVj7OrD+ONSIp7t2hzjmKzKl1EPr/2vw+XcWqT7DUZS2BQYLe2kjqrBUikE+NoroRKAO9kGvPlbOgxGJq2yJAhA92nAQ5OBiz8D3z0B6Bpe9xUiql/MXhqp+PRcPLvmKGJTcvBynwD0a9Ws0tZVJA1BNMD/l2dhEx+JtIDhSGn1NES1VuqwGjwBgKVKgIMK+P5cLuKzDTAaRSi4OIY8dRgPWDkBfy0FvhkAPLsDsG0mdVREJBPMYBqhs7cyMGLlQdxJz8Vr/QPRP8SdyapMCaIRmuRoWCccQVKb55HSejyT1TrmaaPAy+0t8OfNfEQn6VnXKmet/gEMeA9IvQZ8/Qhw97zUERGRTDCLaWR2nLqNUV8ehFEU8fqgVugW4MLlVmVKqc+CQn8PCqMeCR3+ifSAxyEqLaUOq1Ea0UqDhT2tkFcAnE3U4+xdrmcvW827AY8tA/LvAV/3Bc7vkDoiIpIBJqyNhMEoYsl/L2DGj6cQ4GqNtx4LQbiPAxRsEyNLDpe3QJNyAYCAPFtfZHt3BzjBql519VKjpaPpj7wntyZj1+VciSOicrm3AUasAuy9gM3PAr8tAoxGqaMiIgnxb8hGICNHj2k/nsT+y0l4NMQNzz7cHI5ajtTJkmiE1/7X4XR5E4yq92E0amBUq9i6ykw0KgFBjgoUGIFpv6Th6O08/KunPfh/iwzZuAFPfA78+QGw/2PgzhlALAAE/rVF1BTx//wG7viNVMz48RTuZuowoZsfHmvrAY2aDdPlSGHQwTLtMhzTNiMt4Ank320OZEgdVdOjUgj4v75afBmVi/VncxCVkI//aERoVHwbITsqS6Dv24BzIHBsNaBQA26tpI6KiCTAYZ0GqsBgxKd7L+Oprw7BYBQxd3AwnmjnyWRVplS6ZFglnYVgyEdi+GtIaT2BS61KSK0UMO0hLf71sAbX0wpwNlGP1Fy+cpYlQQDajzWthCUagDungb+WsUSAqInhCGsDdCstBzN/PIXjN9PQo6ULnu3aHO72GgisV5UdpS4NmvQrUOlSYFRZQm/tiUwfp/slAJz4I7XezS3Q0kkF9W/A5dQCbP01DW/3soedJX+Xlx3vToBrKyD9JrDvXSBmHzDyG8DWXerIiMgM+KdyA2I0ivjhaCwGL/8L5+9k4oWe/pjRLxAeDlZMVmXI9ub/EPifR6HUpUCvcUG+bXMYVVasV5UZL1sFAp2UcLESEHEhF/03JOLPmzqpw6KyKJSAUwug23Qg7ijweRfg1A8AVzIjavT4N2cDcTUxG2O+Pow3Is7Cx1GLBUNb47G2LAGQI6UuDd6/z4Dfr5NhVNsg384PeutmXGZVxgSY+rUu66eFShDx3I5UxKQVoIBvneUpdAQw7EvA2hX/396dx8lRlgkc/z19T8+dyUxmcpBkQgJJOBLCfSwQBbkWREXwWkXQRddjFXVBFllEV5HPqqyccgRYUQQXFBDC/ZEVAgRC7gTISe65z+7pq579o3qSyZDJQXqmuifP9/Op1NTbb73z1PSk5um33nqLP1/hPh2rZa3XURljBpElrHkukc5w8/PvcfbNL7NicweXnjSBH503lam15Ta/ah4KxrYx5ZHTqFjzOK2TLmTjideTCZbgpkQm302vDvDbs0u4cHKAxpjDwm1JHlkew7EevPxTVQ8X3gHHXeH2tt52vDu2NZPyOjJjzCCwLp88parMXbqVG+euZF1zjBMnVXHxMeOYUFVsc6vmoejW1ylqWoUv1U2qro6ts75LfOThNrdqAQr5ha/PilLZ7mNTl8P3n2/jgcVd3HBaBTO8Ds7szOeHIy+B+tPcZPWF62Hh7+FjP4XJZ7o3bBljhgXrYc1Db61v5VN3zONrDy4go8p3z5jMdz46hfqRJZas5hlfpodw2yomPXkR4iRJlYxh0/HXEa+ZaclqgYsGhcmVfr4zK8SG9jQff7iJeEpxrLM1/5TWwjm/cKfASnTC7z8N950Dmxd6HZkxJkfsL2oeeWdrJze/8C5PLdnKiGiQL504no8eWkt5NOh1aKYfUYdw+xoC8UZAaD30AhLb6lFfCA0Vex2eyaFzJkc4dXyYOYvjpDZAKqnc+FIb3zym1OvQTH8Hz4aJp8Dih2HxQ/DbU+GwT8JpV8PIyV5HZ4zZD5aw5oGFG9q49aVVPLd8G5Ggj4/PGM15R9RRU2pTVeWbYNcmwu1r8aUiBDKNZMKVpCNVNB02C220O8uHq+KQ8I2jo0S3QNKBh5bGeGRZjMoiH2NK/QT99v80b/iDMPNzMO0f4c37YMUTsPRRmHY+nPI9qDvC6wiNMR+CJawecRzlldVN3PG31byyqpnScIALZozmzGmjGFMRtRuq8okq0Ya3qFp6D+Xr5hJI/BD11ZMon4QTiADijqUzw54IhP1wz9nF/M/SHl54P8O2boeqIqGuxA8hryM024XL4KRvwZGfgUV/gHfnwvK/wKTZcPJ3YcLJNsbVmAJiCesQa+1O8qe3NvLg6+tZ1xyjMhrk4qPH8ZGpNYwqi9gY1Twi6R7K1z1N1bJ7iTYuIhMsoX38WSS2TUK7QjgBS1IPVGPK/Fx1YjEbu7pojDm0xpXmeJqlsRR1pX5KM0rIel3zQ0m1m7ge9UVY8kdY8STcf577uNdjLndv2iqq8DpKY8weWMI6BDKOMm91M48u2MiTS7aQTDscWlvKV0+p58RJVYwoDtml/3yhSrRhARXv/YmKNU/gT3aQLBlD09R/onPsqaSj1TgvxoGM15GaPBD2C2NL/dSVQHPcIZFR3mtJc8XdWzl/SpTrHPD5bFKzvFBUDsd+FWZ+Ht552l3m/hs8f507znXGZ+GgE903zBiTdyxhHSSqyoL323hi0WaeXLyZpq4k0ZCfUw4eyamHVDO1tswm/c8XqoRb36F83VwqVv+ZcPsaHH+Y7lGz6Ko7idioo3BCdoONGZhfoCbqY2rIT3tCmRry8Ydl3XwvqPgE7p/fyZdTSlHQUlfPBaNugnrYJ2HLEneYwJI/wcIHoaTWLT/8UzB6pg0ZMCaPWMKaQ/FkhldXN/HiygZeWtnA5vYegn5hxrgKLj5mLLMOcntTbXxqHnDSRBvepmz9s5Stm0u4cz0A8RFTaZx+GV11x5GOjrKxqWaflYeFH59STFvcIfQcpB24aV4ns0IpIn7hmVc6OHNShCNqgnYu8Frd4e6S6ITVL8Hav8Ebd8Jrt0L5QTDlYzDlLHe8azDidbTGHNAsYd0P6YxDdyJNR0+KL977BvPWNJNMO0SCPqaPLuecw+s4ZkIlNaURAn67zOQpVcJtqyjZ/ArFm/9OyZZ5+JOdqASIj5xO07jT6a6ZRap0DOq3O2fM/qso8hH0QdAHD5wTpW6e0JFQ7lzQxe1vdVEeFk4YG+bkg8KcNC7MhHK/DQ3ySrjUnUVg2vkQa4FVz8OG12HB/TD/LghEYOI/QP3pMP5EqD3cPswaM8QsYd0HLd1JlmxqZ/GGNuavb+WtdS10J92xjBkHTptSzZHjyjlsdDnlRdaT6iVJxYg2LaaoYQHRhreJNiwgGG8EIBUdRXfNLHpGTKO7ZgbpaDX4bK5bM3jqygJUR31UR+GhM0p4dWOSRY0Z3tycYO5qdzq0UcU+ZtaGOKouxMzaIIfXhIgE7Bwy5KIj4IhPu0sy5iauG96AzW/De8+6dUKlcNDxMP4EGH0UjJ4BRZWehm3McGcJ6y6kMg7rm2O8t62T9xq6WLGlgyWb2tnYGt9eZ1xlEcfVV7FySwfRkJ+ffeIIwgGf9ZAMNXWQTAJfOk71wluItKwk0rKCcPsaRN0PE8ni0fRUTqFt4nnEq6aTLB+P+iM2Ps14orLIx7mTI5w72Z3ebl1bmje3plnR7PD2lh0JbECgvjLA1JFBDh3prg+pClKL3cQ1ZEJRmHS6u6hC+wbYtAC2LYOG5bDquR11K8bDmKNg1HSoPhSqp0LlBPAP3p/ZmTNncsIJJ3DttddSV1c3aN/HmHxwwCasnT0ptrb3sLEtzvvNMdY3x3i/Jcb7Ld2sbeomlXGfvyjAqLII46uinDSpivFVxRxcU0JlNEQ44OOHjy0BsBuoBpGk44S6NuFPtCOZBKPm/5xw+1pC7WsJd6wlGvs+ALVv/oJUtIZkyVjaJp5HT0U9icoppIuqUX/YElSTd3w+oX5EkPoRbg+/KjR2p1namOadFod1HQ6vbOjhL+/ueB7sw+EkEb9w+1Mt1FcEGFfmZ0yZn9GlAcaU+q1XdrCIQMVB7jL94+6bFWuCbcuh6T1oWQXrX4Vlj+3Yxx+GERNhRL27VE5wt8sPgrI6dyjCfli4cCHLly9nzpw5XHrppZa4mmFt2CSsybRDezyVXZK0xVI0dydp7krS1JWguStBU1eSre09bOmI053YeVqicMBHTWmY6tJwdvL+IkZXFDGuMkpZUZBQwGdzpOaKk8Kf7MSfaMefaCXQ04q/p4VAoo1QZx3ipJjw9I0EYw0EYg0EEq0ARBL/DkD14jtJRWtJRWvoHDeb1MZa1B9m7Sl3kYlUumNQxcYMm8IjAjUlAWaXBJg90S1ThZZ4hlUtadZ1ZKhcJyQyysItCZ5Z3YOjO7cxIiLUFPupLvZTU+yjJuqnKuqjMuJjRJGPioj7dVlYKA357CldH5YIFFdD/anuAqAOxNugZQ20roW2jdCxCbYuccfFZpI7txEuhdLRUDYaSmrc9rYvI91hBpEKd11U4T7Fq59k0m3znnvuYc6cOfT09NwG3KCqWwbz8I0ZagWZsF758CLe3dZBZ0+aroS79KScAetHgj7KIkHKioJURIPUVxdTGQ0yojjEiOIQdeVFjCwJEQr4CfikcC7rqwK6fS3q7LSN9pY52a8VNOOWacbdx8kgmnYvrTsZ93UnjWgGcVLguGt/wm27bO1fkUwKn5NE0gnESeDLJJC0e1m+9/K8L9ODL9W90+JPdeFPduBLxwc8pGDyWlQCBLs2kwmXk6yZSToygnSkiuS741FfkDWnPoAGI6gvCOIj3dwNQLqkeNB/5MYMNRGoivqpivo5DhjT7H4Yu//kUnpS7pO2tnU7bIs5NMWUprjSmnDY0pFhRSO09igZHbj9aEAoCQslIeHmZIr6jLKxOc2tc1spCgiRgBAJChE/RAJCyC+Ee9d+CPqEoF8I+CCUXQd87tovQtDvTvvlE7fMJ265T9wpT32A3yeIuFe0fOIuAohIdp3dT3FfUM3P87T43DGw0REw9ugd5aqQSUBXA7RtcNexFog3Q6wZ2jdCwzKIt7v1BhKIuElu79JHb+IKfBm4VETmYImrGUYKMmH91tqvUdOzBnDPXhIAAjuu+Er2H3FfdcsViGWX5r6t7eZMvheKs5ejp8351E7lsjftat86utNatP/+un2fvWo7x3p7N8e/8JMB6zi+II4/jPau/WEcfxjHHyETriQRrcPxR3ACRWQCUXcdLCYTLCETKiMdKif+ZgUqwtrjf4qKb6fL+Om1aQASEnDnCsI9sau6H1ZSiYET4VxSdbKfB3Z8SEol4gPG0Vve//W+7QwUu1tHd6q3/fvq7o99l9+vX93e9vu3vcvjyf5O7s3Pe1/fi4F+Nn3bGyjWvvt/sA3d6ee7p/foAz+bAY6lf3t72md3vx/7qreddDJOABgTdReX0HeEq+MoyYxDLAVtCaW9R2lPQlcKYmnoTun2dSIDklQUiKcc5m+Kk8xAIgPJDCQH7hMYMg+F3KTskt/syMF6E9rtX/fdQXZa7fRa/3x3l+nv3tTZJwFgIsLEAWtE6aGKdiropFy6KKObcroop5sSJ05JMk5xZ5wSBvw9CmfXlriaYUX0A4lR/hGRRmB937KSolBlfn7E3iEHweX8+D7QoOzVeXqfqftPAfx2GWMGJgOcIfaR7morJ2eHA/oU0xlLDJz57uAAf1fVUwc7HmMGU0H0sKpqtdcxGGOMMflERHaXsCdwk9U5wA1DE5Exg6cgElZjjDHG7JWdElVV3epxPMbkhCWsxhhjTOGzRNUMa5awGmOMMYXLElVzQLCE1RhjjClMC4FXsUTVHAAKYpYAY4wxxhhz4DqgHgckIleKiIrISK9jyQURuUFEFovIQhF5VkRGex1TLojITSKyMntsj4lIhdcx7S8RuUhElomIIyJH73mP/CUiZ4nIOyKySkSu8jqeXBCRe0WkQUSWeh1LLonIOBF5SUSWZ3//vu11TPtLRCIi8oaILMoe0/Vex5RLIuIXkbdF5EmvYzEmnxwwCauIjAPOBN73OpYcuklVj1DVGcCTwI88jidXngMOU9UjgHeBqz2OJxeWAp8AXvY6kP0hIn7gVuBsYBrwGRGZ5m1UOXEfcJbXQQyCNHClqk4Djgf+ZRi8XwlgtqoeCcwAzhKR470NKae+DazwOghj8s0Bk7ACvwJ+wDCaaFpVO/psFjNMjk1Vn1XVdHbzNWCsl/HkgqquUNV3vI4jB44FVqnqGlVNAg8BF3gc035T1ZeBFq/jyDVV3aKqC7Jfd+ImQmO8jWr/qKsruxnMLsPi3CciY4Fzgbu9jsWYfHNAJKwicgGwSVUXeR1LronIT0VkA/A5hk8Pa19fBp72Ogiz3RhgQ5/tjRR4AnSgEJEJwEzgdY9D2W/Zy+YLgQbgOVUt+GPK+jVux0oePAjXmPwybGYJEJHngdpdvHQN8EPc4QAFZ3fHpap/UdVrgGtE5GrgG8B1Qxrgh7Sn48rWuQb3kuaDQxnbh7U3x2SMF0SkBPhf4F/7XZkpSKqaAWZkx7c/JiKHqWpBjz8WkfOABlV9S0RO8zgcY/LOsElYVfWjuyoXkcOBicAiEQH38vICETm2EKYBGei4duFB4CkKJGHd03GJyJeA84CPaIFMZbEP71Uh2wSM67M9Nltm8pSIBHGT1QdV9VGv48klVW0TkZdwxx8XdMIKnAScLyLnABGgTER+p6qf9zguY/LCsB8SoKpLVLVGVSeo6gTcS5hHFUKyuiciMrnP5gXASq9iySUROQv3stj5qhrzOh6zk/nAZBGZKCIh4BLgcY9jMgMQ91P6PcAKVf2l1/HkgohU984cIiJFwBkMg3Ofql6tqmOzf6cuAV60ZNWYHYZ9wjrM/VxElorIYtwhDwU/ZU3WLUAp8Fx2yq47vA5of4nIhSKyETgB+KuIPON1TB9G9ma4bwDP4N7A87CqLvM2qv0nIn8A5gGHiMhGEbnM65hy5CTgC8Ds7P+lhdkevEJWB7yUPe/Nxx3DalNAGTPM2YMDjDHGGGNMXrMeVmOMMcYYk9csYTXGGGOMMXnNElZjjDHGGJPXLGE1xhhjjDF5zRJWY4wxxhiT1yxhNcYYY4wxec0SVmMKkIj8UURaRKS2X7lfROaLyHvZSdWHIpb7RET7LI0i8nL2ARD9684SkZiIjMlu14lIa3Ye1F3VTYvIT7LbIiJvi8gPBv+ojDHG5BNLWI0pTN8EFLitX/n3gFnA5aoaH8J4GnEfinAC8BVAgKdE5CP96t0E3KuqmwBUdQvwfeASETm/t1L2caL3Au8CN2TrKvBj4GoRGTG4h2OMMSafWMJqTAFS1QbgO8CFInIRgIhMAf4DuFNV/zbEISVV9bXs8mfgXKCDPk9fE5FZwOnA7X13VNW7gReB20WkPFt8NXAYcJmqJvpUfxzoAS4frAMxxhiTfyxhNaZAqeoDwFzgFhEZifvM+EbA80vmqtqB2zt6cJ/iy4HFAzzK9StABfBfIjIduAb4jarO69duBngES1iNMeaAYgmrMYXtn4Eo8DpwMnCFqnZ6GxKISAAYB7T1KT4L+L9d1VfVNcC1wGXAE8Am3KR1V14GJotIfa7iNcYYk98sYTWmgKnq+8AtQD3wqKo+lcv2RSS8D3UD2WUscCtQCzycfW0UMAFYtJsmfo2bqE4Evq2q3QPUezu7Pn5vYzPGGFPYLGE1poCJSBnwBdwbsI4RkdJ92Pd8Efm7iCSzMw48JiKfE5FaEakUkc8Dz+5lc2OAVHbZAHwW+BHw39nXR2fXjbtp48JsOwr0v1mrr942Ru+mjjHGmGHEElZjCttNQCXuTU41wM/2ZicROQi4A3gUmA1cAXTh9oxuAVpw787vPwvBQBqAY4CjcXtIK1T1BlV1sq9HsuvErnYWkUrcnuKncWcC+KaIHDvA9+qd/WBIpu0yxhjjvYDXARhjPhwROQ33ZqUrVfXp7Hyl14vI71X11T3s3gxMV9XWPmUPZ8eeHgykgdXZqaT2RkpV39zD9wM3ud6VXwIluInzVuAi4G4RmaWqqX51e6e0atrL2IwxxhQ462E1pgBlHwpwFzAfuDlbfCOwDDfRC+1uf1Xt7pes9panVXWlqq7ah2R1b6zDnY7qAzdKicgZwJeAq1X1fVVN4s4CMB24ahdtTcyu38lhfMYYY/KYJazGFKYfA+NxHxDgAGR7Ii8HDmHgO+w9kU1CXwd2uswvIsXAb4FX6TP8IDud1W3ANSJyaL/mjsMdK/vaYMZsjDEmf1jCakyBEZGjcR8a8HNVXdL3NVV9A7fH9arsfKb55I/A7GyS2us/cW+eurzPeNdeVwPbgLtERPqUnwc8rqqxQY3WGGNM3pDcXvUzxphdy85osBH4uqr+7kO2MRp3FoIzVfWFXMZnjDEmf1nCaowZMiJyDXAxcOSHGSMrIr/K7js758EZY4zJWzZLgDFmKP0S8AN1wOZ92TE7LGAr7phXY4wxBxDrYTXGGGOMMXnNbroyxhhjjDF5zRJWY4wxxhiT1yxhNcYYY4wxec0SVmOMMcYYk9csYTXGGGOMMXnt/wH+IEjKjtFkSwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(Xs, Xt, ys, yt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eedea5b-9e4f-438d-903b-ceaaad52709e",
   "metadata": {},
   "source": [
    "### Question 1.1 :\n",
    "\n",
    "Use `scikit-klearn` to fit a linear model on the source labeled data `(Xs, ys)` and show the model predictions on the plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "371270ee-3077-4b27-a17c-d2dbdd9558f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "model = #...\n",
    "model.fit(#...); \n",
    "# Be carefull! You need to reshape Xs as a column vector : Xs.reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc36359d-48a5-45da-bff9-88f62bf1d5b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(Xs, Xt, ys, yt, model=model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5079fbba-1ad5-4f6f-8786-5142b54e63ef",
   "metadata": {},
   "source": [
    "### Question 1.2 :\n",
    "\n",
    "Compute the source and target empirical average squared errors : \n",
    "\n",
    "$$\\mathbb{E}_{(x, y) \\sim \\widehat{P_S}(X, Y)}[(h(x)-y)^2]$$\n",
    "\n",
    "$$\\mathbb{E}_{(x, y) \\sim \\widehat{P_T}(X, Y)}[(h(x)-y)^2]$$\n",
    "\n",
    "What do you observe ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6355162a-0131-4091-a5c5-5fda21ae80eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ! Do not forget to reshape Xs and Xt with Xs.reshape(-1, 1) to make the prediction\n",
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d3e12d-fb7d-4a87-be82-f73068faf14b",
   "metadata": {},
   "source": [
    "### Question 1.3 :\n",
    "\n",
    "Assuming that the model error $\\epsilon$ is a centered gaussian ($Y = X \\beta + \\epsilon$ with $\\epsilon \\sim \\mathcal{N}(0, \\sigma)$), compute the pvalue of the F-test that the source and target computed errors are drawn from the same law $\\epsilon$.\n",
    "\n",
    "For a risk $\\alpha = 0.01$, can we reject the hypothesis that $\\epsilon_s \\sim \\epsilon_t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a0903a1-2055-478a-abda-18e6ddf78e94",
   "metadata": {},
   "outputs": [],
   "source": [
    "F = err_src / err_tgt\n",
    "pvalue = scipy.stats.f.cdf(F, len(ys)-1, len(yt)-1)\n",
    "print(\"pvalue = %f\"%pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdce2dde-60fb-44ad-8417-ce299df0ade8",
   "metadata": {},
   "source": [
    "We now want to reweight our source risk to correct the shift between domains:\n",
    "\n",
    "### Question 1.4 :\n",
    "\n",
    "Show that we have:\n",
    "\n",
    "$$\\mathbb{E}_{(x, y) \\sim P_T(X, Y)}[(h(x)-y)^2] = \\mathbb{E}_{(x, y) \\sim P_S(X, Y)}[w(x) (h(x)-y)^2]$$\n",
    "\n",
    "With:\n",
    "\n",
    "$$w(x) = \\frac{p_s(t)}{p_s(x)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4909934-7899-4819-b368-c5a03fb2d292",
   "metadata": {},
   "source": [
    "### Question 1.5\n",
    "\n",
    "Express the density ratio $w(x)$ as a function of $x$ and write a python function for it.\n",
    "\n",
    "Then plot a comparison of the empirical source distribution with its reweighted counterpart."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f183000d-4003-4607-a1f9-fdc3fa08046a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def density_ratio(x):\n",
    "    return #..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ae2badf-a636-40da-bf71-166a8f9ad790",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(Xs, Xt)\n",
    "show(Xs, Xt, weights=density_ratio(Xs))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4104066-612a-4678-8061-98a66bb35c3b",
   "metadata": {},
   "source": [
    "### Question 1.6\n",
    "\n",
    "Now that we have computed $w(x)$, we can give a \"debiased\" estimation of the target risk by reweighting the source risk.\n",
    "\n",
    "Compute a reweighted average of the source errors and compare it to the target average error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98a57391-c94f-433b-bd02-8b83b1b751a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source reweighted average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb44ad08-f624-4262-ad7a-32af5859985d",
   "metadata": {},
   "source": [
    "### Question 1.7\n",
    "\n",
    "Compute the F statistic. Can we reject the hypothesis that $\\epsilon_s \\sim \\epsilon_t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f91e7fd9-6663-4a4e-8c74-3353c16c3f6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "F = err_src / err_tgt\n",
    "pvalue = scipy.stats.f.cdf(F, len(ys)-1, len(yt)-1)\n",
    "print(\"pvalue = %f\"%pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de9cf18c-30fb-417c-ad57-276c5edfb8de",
   "metadata": {},
   "source": [
    "### Question 1.8\n",
    "\n",
    "We know want to fit the linear regression model based on the reweighted source risk.\n",
    "\n",
    "Compute analyticaly the linear regression solution $\\beta^*$\n",
    "\n",
    "**Hint** : Use the `X, y` arrays defined below. The \"full ones\" vector added to the input matrix is usefull to fit the intecept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91ddb5a7-225e-4e34-892d-46cdfe05bd53",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.concatenate((Xs.reshape(-1, 1), np.ones((Xs.shape[0], 1))), axis=1)\n",
    "y = ys.reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "653a09a5-2336-4fc4-8690-6cd65d1c9544",
   "metadata": {},
   "outputs": [],
   "source": [
    "Delta = np.diag(density_ratio(Xs)) \n",
    "\n",
    "beta = np.linalg.inv(X.transpose() @ Delta @ X) @ X.transpose() @ Delta @ y\n",
    "\n",
    "print(\"Linear Regression coeficient\", beta.ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12d11d28-6c1b-485e-a241-373665167820",
   "metadata": {},
   "source": [
    "Define a debiased_model function which multiply the input matrix `x` with the $\\beta$ coeficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c63c276a-d036-498e-89b0-4830ec4022f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def debiased_model(x):\n",
    "    x = np.concatenate((x, np.ones((x.shape[0], 1))), axis=1)\n",
    "    return x.dot(beta).ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75cd7ade-0115-4648-88b7-220b338e056e",
   "metadata": {},
   "source": [
    "### Question 1.9\n",
    "\n",
    "Plot the model predictions and compute the new target risk. What do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8306edf2-81c0-4912-abf7-c39083678505",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(Xs, Xt, ys, yt, model=debiased_model, weights=density_ratio(Xs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4630d72b-1bf9-4091-a6e5-d0abe2ea459a",
   "metadata": {},
   "outputs": [],
   "source": [
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source reweighted average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48532549-f3c2-4600-83d0-94eb1ac39205",
   "metadata": {},
   "source": [
    "### Question 1.10\n",
    "\n",
    "Perform a weighted regression by with the `sample_weight` argument from the `fit` method of the `LinearRegression` object from `scikit-learn`. Verify that you get the same results than previously"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "328db9af-8370-455e-82d4-6fc145eebddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_2 = #...\n",
    "model_2.fit(#...);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7549264-a5b7-4e5c-8f37-6c29ace8be33",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(Xs, Xt, ys, yt, model_2, density_ratio(Xs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d90675c-8e00-4297-b0eb-de413cc8702f",
   "metadata": {},
   "outputs": [],
   "source": [
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source reweighted average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75921f99-442d-4d19-8586-10e9e5d8a1f6",
   "metadata": {},
   "source": [
    "## 2. Instance-based Approaches"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ae5496c-f0cc-471e-9bd2-5606712947e3",
   "metadata": {},
   "source": [
    "### Question 2.1 : PDE\n",
    "\n",
    "Implement the Paramatric Density-Estimation (PDE) method and compare the computed weights with the theoretical density ratio.\n",
    "\n",
    "**Hints** Use the `gaussian` function defined above with modified arguments `mu` and `s` to compute the estimated density functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab19c6e3-3ce6-45cb-ae2a-b4a3b76d6a92",
   "metadata": {},
   "outputs": [],
   "source": [
    "ps_Xs = #...\n",
    "pt_Xs = #...\n",
    "\n",
    "weights = pt_Xs / ps_Xs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7cf9cb5-4270-4d56-b5f9-d2b3da04b999",
   "metadata": {},
   "outputs": [],
   "source": [
    "args = np.argsort(Xs).ravel()\n",
    "plt.plot(Xs[args], density_ratio(Xs[args]), label=r\"$w(x)$\")\n",
    "plt.plot(Xs, weights, \".\", c=\"black\", label=r'$\\widehat{w}(x)$')\n",
    "plt.xlabel(\"X\"); plt.ylabel(\"Weight\"); plt.legend(fontsize=16); plt.show()\n",
    "\n",
    "show(Xs, Xt, weights=weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f2df112-50f4-4443-b14b-177b94353c69",
   "metadata": {},
   "source": [
    "### Question 2.2 : KDE (1)\n",
    "\n",
    "In many cases, the source and target input distributions are not gaussians, to estimate the densities, one can then use a kernel density approach.\n",
    "\n",
    "Compute the Kernel Density Estimation (KDE) with gaussians kernels for the source and target densities, plot the estimated densities and compare the obtained weights to the theoretical one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99da4ddd-a3b4-46fd-9916-3c8245729cb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics.pairwise import rbf_kernel\n",
    "\n",
    "def kde(centers, x, sigma=0.5):\n",
    "    # dist is a pairwise distance matrix of shape (len(centers), len(x))\n",
    "    dist = (1./np.sqrt( 2. * np.pi * sigma**2 )) * rbf_kernel(centers, x, gamma=1/(2 * sigma))\n",
    "    return dist.mean(0) # Average over the centers\n",
    "\n",
    "def log_likelihood(p):\n",
    "    return np.mean(np.log(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b14df28-6ece-4269-a135-63845482aff8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 0.5  # kernel bandwidth\n",
    "\n",
    "x = np.linspace(-4.2, 4.2, 100)\n",
    "ps_x = # ... ! You need to give Xs.reshape(-1, 1) and x.reshape(-1, 1) to kde\n",
    "pt_x = # ... ! You need to give Xt.reshape(-1, 1) and x.reshape(-1, 1) to kde\n",
    "\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.plot(x, ps_x); plt.plot(x, pt_x)\n",
    "plt.fill_between(x, ps_x, alpha=0.2)\n",
    "plt.fill_between(x, pt_x, alpha=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b33cc0a7-3438-4678-a020-f8571dadce56",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = #...\n",
    "\n",
    "args = np.argsort(Xs).ravel()\n",
    "plt.plot(Xs[args], density_ratio(Xs[args]), label=r\"$w(x)$\")\n",
    "plt.plot(Xs, weights, \".\", c=\"black\", label=r'$\\widehat{w}(x)$')\n",
    "plt.xlabel(\"X\"); plt.ylabel(\"Weight\"); plt.legend(fontsize=16); plt.show()\n",
    "show(Xs, Xt, weights=weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bcda40b-38eb-4d95-93c3-6a7476225468",
   "metadata": {},
   "outputs": [],
   "source": [
    "ll_src = log_likelihood(kde(Xs.reshape(-1, 1), Xs.reshape(-1, 1), sigma=sigma))\n",
    "\n",
    "print(\"Source KDE Log likelihood: %.3f\"%ll_src)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe42ff78-7ef8-40a0-9d6b-448d8b2dd7e7",
   "metadata": {},
   "source": [
    "### Question 2.3 : KDE (2)\n",
    "\n",
    "Try different value of kernel bandwidth `sigma`, what do you observe?\n",
    "\n",
    "**(BONUS)** Compute the log likelihood for multiple `sigma` values, in which case the log likelihood is maximized?\n",
    "Can we use the log likelihood to select a relevant value of `sigma`? What improvement can you propose to find the best `sigma` ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "493bcc31-9fac-481c-a2cc-dbdf1bef520d",
   "metadata": {},
   "source": [
    "### Question 2.4 NNW (1)\n",
    "\n",
    "We now consider the instance-based method \"Nearest-Neighbors-Weighting\" which derives the source weights based on the pairwise distances between source and target instances.\n",
    "\n",
    "Import the `NearestNeighborsWeighting` object from the [adapt.instance_based](https://adapt-python.github.io/adapt/contents.html#adapt-instance-based) module. \n",
    "You can access the documentation [here](https://adapt-python.github.io/adapt/generated/adapt.instance_based.NearestNeighborsWeighting.html). Then fit the weights with the `fit_weights` method and plot the results\n",
    "\n",
    "Try different value of `n_neighbors` in the range [1, 30], what do you observe? How do you explain the resulting weights for `n_neighbors=30` ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f466f83-8f1d-402f-8608-b08bd3dbd088",
   "metadata": {},
   "outputs": [],
   "source": [
    "from adapt.instance_based import NearestNeighborsWeighting\n",
    "\n",
    "nnw = NearestNeighborsWeighting(n_neighbors=5)\n",
    "weights = nnw.fit_weights(Xs.reshape(-1, 1), Xt.reshape(-1, 1))\n",
    "weights /= weights.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56113383-22d7-4ea2-8628-8cb77577125c",
   "metadata": {},
   "outputs": [],
   "source": [
    "args = np.argsort(Xs).ravel()\n",
    "plt.plot(Xs[args], density_ratio(Xs[args]), label=r\"$w(x)$\")\n",
    "plt.plot(Xs, weights, \".\", c=\"black\", label=r'$\\widehat{w}(x)$')\n",
    "plt.xlabel(\"X\"); plt.ylabel(\"Weight\"); plt.legend(fontsize=16); plt.show()\n",
    "\n",
    "show(Xs, Xt, weights=weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d33f691b-ce3b-4b3e-98d7-a988f378dff2",
   "metadata": {},
   "source": [
    "### Question 2.5 NNW (2)\n",
    "\n",
    "Every instance-based methods implemented in `adapt` can fit directly an estimator with the learned weights.\n",
    "\n",
    "Instantiate a `NearestNeighborsWeighting` object with a `LinearRegression` estimator and the target input data `Xt`, then use the `fit` method of the model on `Xs, ys`.\n",
    "\n",
    "Compute the target error, what do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdfb860c-9cdb-44be-8f99-3b4f56d86355",
   "metadata": {},
   "outputs": [],
   "source": [
    "nnw = NearestNeighborsWeighting(estimator=LinearRegression(),\n",
    "                                Xt=Xt.reshape(-1, 1),\n",
    "                                n_neighbors=1)\n",
    "\n",
    "nnw.fit(Xs.reshape(-1, 1), ys);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d7eb7d0-778c-4428-b37f-191d003735f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(Xs, Xt, ys, yt, nnw, nnw.weights_/nnw.weights_.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7761d881-bbe4-466a-b9ae-b99cb0dde133",
   "metadata": {},
   "outputs": [],
   "source": [
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source reweighted average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c431c4f7-8e33-4b3b-b2ab-f8f888fee118",
   "metadata": {},
   "source": [
    "### Question 2.6 KMM (1)\n",
    "\n",
    "We now consider the Kernel Mean matching method.\n",
    "\n",
    "Import `KMM` from `adapt.instance_based` and fit the weights. What do you observe? How do you explain that the estimated density ratio $\\widehat{x}$ is not continuous with respect to $x$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "972610de-e74b-4f5d-8445-fa5dd1d8e328",
   "metadata": {},
   "outputs": [],
   "source": [
    "from adapt.instance_based import KMM\n",
    "\n",
    "kmm = KMM()\n",
    "weights = #... Use the fit_weights method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c986be55-5465-4b8e-93d6-58d9494367a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "args = np.argsort(Xs).ravel()\n",
    "plt.plot(Xs[args], density_ratio(Xs[args]), label=r\"$w(x)$\")\n",
    "plt.plot(Xs, weights, \".\", c=\"black\", label=r'$\\widehat{w}(x)$')\n",
    "plt.xlabel(\"X\"); plt.ylabel(\"Weight\"); plt.legend(fontsize=16); plt.show()\n",
    "\n",
    "show(Xs, Xt, weights=weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "938adc03-40a5-485a-8e87-21a2f0783174",
   "metadata": {},
   "source": [
    "### Question 2.7\n",
    "\n",
    "Use the `fit` method to fit a linear model with the learned weights. Compute the target error, what do you observe?\n",
    "\n",
    "**(BONUS)** How do you explain that KMM and NNW provide lower target risk than the weighted regression using the optimal density ratio ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ffe73b4-dfba-4101-9f42-63ac83041fc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "kmm = KMM(estimator=LinearRegression(), Xt=Xt.reshape(-1, 1), verbose=0, random_state=0)\n",
    "kmm.fit(#...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c29488-1654-456c-be63-906812a524b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(Xs, Xt, ys, yt, kmm, kmm.weights_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c07d6048-80f7-48a6-8670-a6c717089a24",
   "metadata": {},
   "outputs": [],
   "source": [
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source reweighted average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b086cc8-59ac-438b-93ac-dfd2cd8c76cd",
   "metadata": {},
   "source": [
    "# Exercise 2 : CT-Scan Dataset\n",
    "\n",
    "We now consider a real-life transfer scenario using the [CT-Scan dataset](https://archive.ics.uci.edu/ml/datasets/Relative+location+of+CT+slices+on+axial+axis) from the UCI machine leanring repository.\n",
    "\n",
    "This dataset is composed of more than $50000$ Ct-Scan slices coming from around $100$ different patients. The goal is to predict the relative location of the slice in the patient based on features extracted from the slice (cf Figure below)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e34b98ee-253c-448c-8833-bd139fd30fba",
   "metadata": {},
   "source": [
    "![ctscan](ctscan.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c90f0d4b-6c4b-47f7-ae9b-b327d8887079",
   "metadata": {},
   "source": [
    "We consider the following transfer learning scenario, where the learner has access to labeled data for every patients except one:  discard the patient with the more data and consider it as a target patient. We are in the realistic situation where an expert (radiologist) has labeled the data for many patients, but the model is applied on new coming patient for which no labels is available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5eee38d7-cc7a-43b9-a9c9-e247f32b1463",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00858363-bc1b-451d-a825-62fa78f287eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"https://archive.ics.uci.edu/ml/machine-learning-databases/00206/slice_localization_data.zip\")\n",
    "\n",
    "print(df.shape)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b13ac7d1-3bf5-420c-80e3-335e2a5b2426",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df.drop([\"patientId\", \"reference\"], axis=1)\n",
    "y = df[\"reference\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8b1d02f-d613-412e-b7f6-2c298df3e68f",
   "metadata": {},
   "outputs": [],
   "source": [
    "source_ids = [i for i in range(97) if i != 46]\n",
    "target_ids = [46]\n",
    "\n",
    "X_src = X.loc[df.patientId.isin(source_ids)]\n",
    "X_tgt = X.loc[df.patientId.isin(target_ids)]\n",
    "\n",
    "y_src = y.loc[df.patientId.isin(source_ids)]\n",
    "y_tgt = y.loc[df.patientId.isin(target_ids)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81116b3b-cce5-4471-81fd-c4f8d5e25a05",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Shape of Xs\", X_src.shape)\n",
    "print(\"Shape of Xt\", X_tgt.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55f5846d-d48f-4033-9041-b5370c47618a",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X_src, y_src, train_size=0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fb64e66-4bcb-404b-a90d-831cbcbfa3cc",
   "metadata": {},
   "source": [
    "### Question 1\n",
    "\n",
    "The plots below show the two first components of the PCA of the 384 input features. We compare the difference between the train and test set which are uniform splits of the source dataset and the difference between the source and target dataset.\n",
    "\n",
    "What do we observe? Can we consider this problem as a transfer learning problem?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1711aba-653c-4873-b1f2-a202176c6a83",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "X_pca = PCA(2).fit_transform(X)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))\n",
    "ax1.plot(X_pca[X_train.index, 0], X_pca[X_train.index, 1], 'o', label=\"Train\", markeredgecolor=\"w\", ms=5)\n",
    "ax1.plot(X_pca[X_test.index, 0], X_pca[X_test.index, 1], 'o', label=\"Test\", markeredgecolor=\"w\", ms=5)\n",
    "\n",
    "ax2.plot(X_pca[X_src.index, 0], X_pca[X_src.index, 1], 'o', label=\"Source\", markeredgecolor=\"w\", ms=5)\n",
    "ax2.plot(X_pca[X_tgt.index, 0], X_pca[X_tgt.index, 1], 'o', label=\"Target\", markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "\n",
    "ax1.legend(); ax2.legend()\n",
    "ax1.set_title(\"PCA Train-Test\"); ax2.set_title(\"PCA Source-Target\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "568f725a-6368-49d7-8f5c-6a999b57d3db",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))\n",
    "ax1.hist(y_train, density=True, bins=30, label=\"Train\")\n",
    "ax1.hist(y_test, alpha=0.5, density=True, bins=30, label=\"Test\")\n",
    "\n",
    "ax2.hist(y_src, label=\"Source\", density=True, bins=30)\n",
    "ax2.hist(y_tgt, label=\"Target\", density=True, bins=30, alpha=0.5, color=\"red\")\n",
    "\n",
    "ax1.legend(); ax2.legend()\n",
    "ax1.set_title(\"Output Distribution Train-Test\"); ax2.set_title(\"Output Distribution Source-Target\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "906226c1-b64a-4f98-83e4-9db3f7777628",
   "metadata": {},
   "source": [
    "### Question 2\n",
    "\n",
    "Fir a `Ridge` regression model on the training data and compute the source risk (on the test set) and the target risk.\n",
    "\n",
    "What can be observed?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37719d19-d45d-4fdd-90b4-63c025d04e9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "\n",
    "model = #...\n",
    "model.fit(#...);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f2d267f-4253-4840-995f-54ce4f2f4449",
   "metadata": {},
   "outputs": [],
   "source": [
    "yps = # prediction test\n",
    "ypt = # prediction tgt\n",
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)\n",
    "\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(y_test, yps, 'o', label=\"Source - MSE : %.3f\"%err_src, markeredgecolor=\"w\", ms=5)\n",
    "plt.plot(y_tgt, ypt, 'o', label=\"Target - MSE : %.3f\"%err_tgt, markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "plt.plot([np.min(y_src), np.max(y_src)], [np.min(y_src), np.max(y_src)], c=\"k\", alpha=0.5)\n",
    "plt.ylabel(r\"$\\widehat{Y}$\", fontsize=14); plt.xlabel(r\"Y\", fontsize=14); plt.legend(fontsize=14); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a259582-0206-4db7-a2e3-f99b51f72271",
   "metadata": {},
   "source": [
    "### Question 3\n",
    "\n",
    "Can we consider that the source risk is a good estimation of the target risk? is the covariate-shift assumption realistic in this case?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48fc2573-f8f0-4ad2-8cc0-4aafabb5d556",
   "metadata": {},
   "source": [
    "### Question 4 : KDE\n",
    "\n",
    "We know try to apply KDE to the problem.\n",
    "\n",
    "Try to compute an estimation of $p_s(x)$ with different source sample sizes. Compute the computational time for each size. What is the complexity of the algorithm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee795a1e-6c35-4c26-8ad4-b68915110dd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "n = 100\n",
    "t0 = time.time()\n",
    "X_s = X_src.sample(n)\n",
    "ps_xs = kde(X_s, X_s)\n",
    "t1 = time.time()\n",
    "print(\"Time %.3f\"%(t1-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7bea247-3150-45ed-b53f-73da5188bfd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "times = []\n",
    "n_samples = [100, 500, 1000, 2000, 5000, 10000]\n",
    "\n",
    "t0 = time.time()\n",
    "\n",
    "for n in n_samples:\n",
    "    # ...compute time and append in times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f98a822c-66f2-4396-9c8d-5acf922f590b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(times, \".\")\n",
    "plt.plot(times[-1] * np.array(n_samples)**2 / n_samples[-1]**2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c39ce4e6-85c2-490b-896b-f6bc7474bc90",
   "metadata": {},
   "source": [
    "### Question 5 : NNW (1)\n",
    "\n",
    "What is the complexity of the NNW algorithm, is it more adapted than KDE to our problem?\n",
    "\n",
    "Fit a reweighted Ridge regression using the weights from NNW. Plot the weights on the PCA. What do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cc89074-2cb2-4fd9-9969-a7cb61c0e1e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from adapt.instance_based import NearestNeighborsWeighting\n",
    "\n",
    "nnw = NearestNeighborsWeighting(estimator=Ridge(), Xt=X_tgt, n_neighbors=5)\n",
    "\n",
    "nnw.fit(X_train, y_train);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a857dac-ae66-47db-8a62-102117b4a606",
   "metadata": {},
   "outputs": [],
   "source": [
    "yps = # prediction test\n",
    "ypt = # prediction tgt\n",
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)\n",
    "\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(y_test, yps, 'o', label=\"Source - MSE : %.3f\"%err_src, markeredgecolor=\"w\", ms=5)\n",
    "plt.plot(y_tgt, ypt, 'o', label=\"Target - MSE : %.3f\"%err_tgt, markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "plt.plot([np.min(y_src), np.max(y_src)], [np.min(y_src), np.max(y_src)], c=\"k\", alpha=0.5)\n",
    "plt.ylabel(r\"$\\widehat{Y}$\", fontsize=14); plt.xlabel(r\"Y\", fontsize=14); plt.legend(fontsize=14); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f06ea6c-9ee6-44f3-9ef5-23302b8d6fff",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = nnw.weights_\n",
    "weights /= weights.mean()\n",
    "\n",
    "X_pca = PCA(2).fit_transform(X)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(7, 6))\n",
    "\n",
    "ax.scatter(X_pca[X_train.index, 0], X_pca[X_train.index, 1], label=\"Source\", edgecolor=\"w\", s=5 * weights)\n",
    "ax.plot(X_pca[X_tgt.index, 0], X_pca[X_tgt.index, 1], 'o', label=\"Target\", markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "\n",
    "ax1.legend(); ax2.legend()\n",
    "ax1.set_title(\"PCA Train-Test\"); ax2.set_title(\"PCA Source-Target\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1232a67d-f911-4f73-8196-39f547031992",
   "metadata": {},
   "source": [
    "### Question 6 : NNW (2)\n",
    "\n",
    "Try different vlaues of `n_neighbors`. What is the best value?\n",
    "\n",
    "**(BONUS)** Propose some ideas to select the best value `n_neighbors`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33cc1094-642a-4550-a1fc-9f7b2375684e",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Question 7 : KLIEP (1)\n",
    "\n",
    "Try the same thing with KLIEP. Use different values of the `gamma` argument, what do you observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b21fde9e-4e5d-45b5-9934-aa1d3ee23a6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from adapt.instance_based import KLIEP\n",
    "\n",
    "kliep = KLIEP(estimator=Ridge(), Xt=X_tgt, gamma=1.)\n",
    "\n",
    "kliep.fit(X_train, y_train);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1428df53-b901-4b1f-8ae9-f7d0bed9e488",
   "metadata": {},
   "outputs": [],
   "source": [
    "yps = # prediction test\n",
    "ypt = # prediction tgt\n",
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)\n",
    "\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(y_test, yps, 'o', label=\"Source - MSE : %.3f\"%err_src, markeredgecolor=\"w\", ms=5)\n",
    "plt.plot(y_tgt, ypt, 'o', label=\"Target - MSE : %.3f\"%err_tgt, markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "plt.plot([np.min(y_src), np.max(y_src)], [np.min(y_src), np.max(y_src)], c=\"k\", alpha=0.5)\n",
    "plt.ylabel(r\"$\\widehat{Y}$\", fontsize=14); plt.xlabel(r\"Y\", fontsize=14); plt.legend(fontsize=14); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99fbf82b-d241-440a-bf43-bfbbc1c9ff10",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = kliep.weights_\n",
    "weights /= weights.mean()\n",
    "\n",
    "X_pca = PCA(2).fit_transform(X)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(7, 6))\n",
    "\n",
    "ax.scatter(X_pca[X_train.index, 0], X_pca[X_train.index, 1], label=\"Source\", edgecolor=\"w\", s=5 * weights)\n",
    "ax.plot(X_pca[X_tgt.index, 0], X_pca[X_tgt.index, 1], 'o', label=\"Target\", markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "\n",
    "ax1.legend(); ax2.legend()\n",
    "ax1.set_title(\"PCA Train-Test\"); ax2.set_title(\"PCA Source-Target\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5639527b-177c-40fb-86ff-34941735712c",
   "metadata": {},
   "source": [
    "### Question 8 : KLIEP (2)\n",
    "\n",
    "Put a list of possible values for gamma, explain what is happening?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "290b4915-1e3b-4dbb-9e22-f66886bfa8a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "kliep = KLIEP(estimator=Ridge(), Xt=X_tgt, gamma=#...)\n",
    "\n",
    "kliep.fit(X_train, y_train);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0242de52-13f9-4797-bc61-b2988b488bf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "yps = # prediction test\n",
    "ypt = # prediction tgt\n",
    "err_src = #...\n",
    "err_tgt = #...\n",
    "\n",
    "print(\"Source average squared error : %.4f\"%err_src)\n",
    "print(\"Target average squared error : %.4f\"%err_tgt)\n",
    "\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(y_test, yps, 'o', label=\"Source - MSE : %.3f\"%err_src, markeredgecolor=\"w\", ms=5)\n",
    "plt.plot(y_tgt, ypt, 'o', label=\"Target - MSE : %.3f\"%err_tgt, markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "plt.plot([np.min(y_src), np.max(y_src)], [np.min(y_src), np.max(y_src)], c=\"k\", alpha=0.5)\n",
    "plt.ylabel(r\"$\\widehat{Y}$\", fontsize=14); plt.xlabel(r\"Y\", fontsize=14); plt.legend(fontsize=14); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ce896a4-d029-4914-91d9-a54264c979e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = kliep.weights_\n",
    "weights /= weights.mean()\n",
    "\n",
    "X_pca = PCA(2).fit_transform(X)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(7, 6))\n",
    "\n",
    "ax.scatter(X_pca[X_train.index, 0], X_pca[X_train.index, 1], label=\"Source\", edgecolor=\"w\", s=5 * weights)\n",
    "ax.plot(X_pca[X_tgt.index, 0], X_pca[X_tgt.index, 1], 'o', label=\"Target\", markeredgecolor=\"w\", ms=5, c=\"red\")\n",
    "\n",
    "ax1.legend(); ax2.legend()\n",
    "ax1.set_title(\"PCA Train-Test\"); ax2.set_title(\"PCA Source-Target\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "statmath38",
   "language": "python",
   "name": "statmath38"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}