synthetic_experiment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fair Dummies: Synthetic Example\n",
    "\n",
    "This notebook implements a the fair dummies framework for learning predictive models that approximately satisfy the equalized odds notion of fairness.\n",
    "\n",
    "Paper: \"Achieving Equalized Odds by Resampling Sensitive Attributes,\" Y. Romano, S. Bates, and E. J. Candès, 2020\n",
    "\n",
    "### Proposed approach\n",
    "\n",
    "__Core idea__: fit a regression function, minimizing\n",
    "\n",
    "$$ \\text{loss = prediction error + distance to equalized odds}$$\n",
    "\n",
    "\n",
    "__Input__: $  \\{(X_i,A_i,Y_i)\\}_{i=1}^n \\sim P_{XAY}$ training data\n",
    "    \n",
    "\n",
    "__Step 1__: sample dummy protected attributes\n",
    "$$\n",
    "\\tilde{A}_i \\sim P_{A|Y}(A \\mid Y=y_{i}) \\quad \\forall \\ i=1,2,\\dots,n \\nonumber\n",
    "$$\n",
    "\n",
    "$A \\in \\{0,1\\}$? generate $\\tilde{A}$ using a biased coin-flip, with\n",
    "$$\n",
    "P\\{A=1|Y=y\\} = \\frac{P\\{y \\mid A=1\\}P\\{A=1\\}}{P\\{y \\mid A=1\\}P\\{A=1\\} + P\\{y \\mid A=0\\}P\\{A=0\\}}\n",
    "$$\n",
    "    \n",
    "__Step 2__: fit a regression function on $\\{(X_i, A_i, Y_i)\\}_{i=1}^n$\n",
    "\n",
    "$$\n",
    "        \\hat{f}(\\cdot) \\,= \\, \\underset{f \\in \\mathcal{F}}{\\mathrm{arg min}} \\, \\frac{1-\\lambda}{n} \\sum_{i=1}^n (Y_i - f(X_i))^2  + \\lambda \\mathcal{D}\\left( [\\hat{\\mathbf{Y}}, \\mathbf{A}, \\mathbf{Y}] , [\\hat{\\mathbf{Y}}, \\tilde{\\mathbf{A}}, \\mathbf{Y}] \\right) \\nonumber\n",
    "$$\n",
    "\n",
    "where\n",
    "    \n",
    "$$\n",
    "\\hat{\\mathbf{Y}} = \\left[f(X_{1}), f(X_{2}), \\dots, f(X_{n})\\right]^T \\ ; \\ \\mathbf{A} = \\left[A_{1}, A_{2}, \\dots, A_{n}\\right]^T \\ ; \\ \\tilde{\\mathbf{A}} = \\left[\\tilde{A}_{1}, \\tilde{A}_{2}, \\dots, \\tilde{A}_{n}\\right]^T \\ ; \\ \\mathbf{Y} = \\left[Y_{1}, Y_{2}, \\dots, Y_{n}\\right]^T\n",
    "$$\n",
    "\n",
    "and $\\mathcal{D}\\left( \\mathbf{Z}_1, \\mathbf{Z}_2 \\right)$ tests whether $P_{Z_1} = \\ P_{Z_2}$ given $\\mathbf{Z}_1, \\mathbf{Z}_2$,  here it implemented as a classifier two-sample test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "import matplotlib\n",
    "matplotlib.use('Agg')\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "np.warnings.filterwarnings('ignore')\n",
    "import random\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import pandas as pd\n",
    "\n",
    "import os\n",
    "import sys\n",
    "\n",
    "\n",
    "sys.path.append(os.path.abspath(os.path.join(os.getcwd() + '/others/third_party/cqr')))\n",
    "base_path = os.getcwd() + '/data/'\n",
    "    \n",
    "from fair_dummies import fair_dummies_learning\n",
    "from fair_dummies import utility_functions\n",
    "\n",
    "seed = 123\n",
    "\n",
    "random.seed(seed)\n",
    "np.random.seed(seed)\n",
    "torch.manual_seed(seed)\n",
    "\n",
    "# maximal number of testpoints to plot\n",
    "max_show = 100\n",
    "\n",
    "# save figures?\n",
    "save_figures = False\n",
    "\n",
    "# display results\n",
    "if save_figures:\n",
    "    SMALL_SIZE = 26\n",
    "else:\n",
    "    SMALL_SIZE = 18\n",
    "\n",
    "MEDIUM_SIZE = SMALL_SIZE\n",
    "BIGGER_SIZE = SMALL_SIZE\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "plt.rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_func(x,\n",
    "              y,\n",
    "              y_u=None,\n",
    "              y_l=None,\n",
    "              pred=None,\n",
    "              point_color=\"b.\",\n",
    "              shade_color=\"\",\n",
    "              method_name=\"\",\n",
    "              title=\"\",\n",
    "              filename=None,\n",
    "              save_figures=False,\n",
    "              show_legend=True):\n",
    "    \n",
    "    \"\"\" Scatter plot of (x,y) points along with the constructed prediction interval \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : numpy array, corresponding to the feature of each of the n samples\n",
    "    y : numpy array, target response variable (length n)\n",
    "    pred : numpy array, the estimated prediction. It may be the conditional mean,\n",
    "           or low and high conditional quantiles.\n",
    "    shade_color : string, desired color of the prediciton interval\n",
    "    method_name : string, name of the method\n",
    "    title : string, the title of the figure\n",
    "    filename : sting, name of the file to save the figure\n",
    "    save_figures : boolean, save the figure (True) or not (False)\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    inds = np.random.permutation(x.shape[0])[:max_show]\n",
    "    x_ = x[inds]\n",
    "    y_ = y[inds]\n",
    "    if y_u is not None:\n",
    "        y_u_ = y_u[inds]\n",
    "    if y_l is not None:\n",
    "        y_l_ = y_l[inds]\n",
    "    if pred is not None:\n",
    "        pred_ = pred[inds]\n",
    "\n",
    "    fig = plt.figure()\n",
    "    inds = np.argsort(np.squeeze(x_))\n",
    "    \n",
    "    if (y_u is not None) and (y_l is not None):\n",
    "        plt.fill_between(x_[inds].squeeze(),\n",
    "                         y_u_[inds],\n",
    "                         y_l_[inds],\n",
    "                         alpha=.2,\n",
    "                         color=shade_color,\n",
    "                         edgecolor=\"\",\n",
    "                         label = u'Prediction interval')\n",
    "    \n",
    "    if pred is not None:\n",
    "        if pred_.ndim == 2:\n",
    "            plt.plot(x_[inds,:], pred_[inds,0], 'k', lw=2, alpha=0.2,\n",
    "                     label=u'Predicted low and high quantiles')\n",
    "            plt.plot(x_[inds,:], pred_[inds,1], 'k', lw=2, alpha=0.2)\n",
    "        else:\n",
    "            plt.plot(x_[inds,:], pred_[inds], point_color, lw=2, alpha=0.8, markersize=15,\n",
    "                     fillstyle='none')\n",
    "    min_val = -35\n",
    "    max_val = 35 \n",
    "    \n",
    "    ident = [min_val, max_val]\n",
    "    plt.xlim(ident[0], ident[1])\n",
    "    plt.ylim(ident[0], ident[1])\n",
    "    \n",
    "    plt.plot(ident,ident, ls=\"--\", c=\".3\")\n",
    "    \n",
    "    plt.xlabel(r'$Y$')\n",
    "    plt.ylabel(r'$\\hat{Y}$')\n",
    "    \n",
    "    ax = plt.gca()\n",
    "    ax.set_aspect('equal', 'box')\n",
    "    \n",
    "    plt.title(title)\n",
    "\n",
    "\n",
    "    if save_figures and (filename is not None):\n",
    "        plt.savefig(filename, bbox_inches='tight', dpi=300)\n",
    "    \n",
    "    plt.show()\n",
    "    \n",
    "def plot_groups_pointwise(x_axis_0,x_axis_1,A,Y,Yhat,\n",
    "                          y_u=None,y_l=None,point_color=\"b.\",shade=\"\",\n",
    "                          filename_0=None,filename_1=None,save_figures=False):\n",
    "   \n",
    "    inds = A==0\n",
    "    \n",
    "    y_u0 = None\n",
    "    y_l0 = None\n",
    "    if (y_u is not None) and (y_l is not None):\n",
    "        y_u0=y_u[inds]\n",
    "        y_l0=y_l[inds]\n",
    "        \n",
    "    plot_func(x_axis_0.reshape(x_axis_0.shape[0],1),\n",
    "              Y[inds], y_u=y_u0, y_l=y_l0, pred=Yhat[inds],\n",
    "              point_color=point_color,shade_color=shade, method_name=\"\", title=r'$A=0$',\n",
    "              filename=filename_0, save_figures=save_figures,show_legend=True)\n",
    "\n",
    "    inds = A==1\n",
    "    \n",
    "    y_u1 = None\n",
    "    y_l1 = None\n",
    "    if (y_u is not None) and (y_l is not None):\n",
    "        y_u1=y_u[inds]\n",
    "        y_l1=y_l[inds]\n",
    "        \n",
    "    plot_func(x_axis_1.reshape(x_axis_1.shape[0],1),\n",
    "              Y[inds], y_u=y_u1, y_l=y_l1, pred=Yhat[inds],\n",
    "              point_color=point_color,shade_color=shade, method_name=\"\", title=r'$A=1$',\n",
    "              filename=filename_1, save_figures=save_figures,show_legend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6000, 2)\n",
      "(6000,)\n"
     ]
    }
   ],
   "source": [
    "def init_seed(seed = 0):\n",
    "    random.seed(seed)\n",
    "    np.random.seed(seed)\n",
    "    torch.manual_seed(seed)\n",
    "\n",
    "def synthetic_example(n = 6000):\n",
    "        \n",
    "    p0 = 0.9\n",
    "    p1 = 1 - p0\n",
    "\n",
    "    sigX_small = 1\n",
    "    sigX_large = 3\n",
    "\n",
    "    A = np.random.binomial(1,p1,n).T\n",
    "    X = np.random.randn(n,2)\n",
    "    X[A==0,0] = X[A==0,0] * sigX_small\n",
    "    X[A==0,1] = X[A==0,1] * sigX_large\n",
    "    X[A==1,0] = X[A==1,0] * sigX_large\n",
    "    X[A==1,1] = X[A==1,1] * sigX_small\n",
    "\n",
    "    beta0 = [0,sigX_large]\n",
    "    beta1 = [sigX_large,0]\n",
    "\n",
    "    Y = np.random.randn(n)\n",
    "    Y[A==0] = Y[A==0] + np.dot(X[A==0],beta0)\n",
    "    Y[A==1] = Y[A==1] + np.dot(X[A==1],beta1)\n",
    "\n",
    "    x_axis_0 = Y[A==0]\n",
    "    x_axis_1 = Y[A==1]\n",
    "    \n",
    "    return X, A, Y, x_axis_0, x_axis_1, beta0, beta1\n",
    "\n",
    "init_seed(seed)\n",
    "X, A, Y, x_axis_0, x_axis_1, beta0, beta1 = synthetic_example(n = 6000)\n",
    "X_cal, A_cal, Y_cal, x_axis_0, x_axis_1, beta0, beta1 = synthetic_example(n = 2000)\n",
    "X_test, A_test, Y_test, x_axis_0_test, x_axis_1_test, beta0, beta1 = synthetic_example(n = 2000)\n",
    "\n",
    "# compute input dimensions\n",
    "n = X.shape[0]\n",
    "in_shape = X.shape[1]\n",
    "    \n",
    "print(X.shape)\n",
    "print(Y.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate $P_{A|Y}$ using kernel density estimator, then compute $P\\{A=1|Y=y\\}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD+CAYAAADcWrmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAANEElEQVR4nO3dwXEbSZYG4Pc22gCI0bz2KtAeUPJg2B6IsxYM6cHIBsoD0oShPBiaQOm+h0bEnjXLwWHvbw8oSGg0iCoWIEJkfl+EIrqqkGAqW6o/KzMrlVUVALTpPw5dAQAORwgANEwIADRMCAA0TAgANOynQ1fgsX7++ed6/fr1oasB8Kx8+vTpX1V1vH7+2YXA69ev4+7u7tDVAHhWMvN/Np03HATQMCEA0DAhANAwIQDQMCEA0DAhANAwIQDQMCEA0LBBL4tl5nlE3HeH06r6sM8ymXlTVWdD6gLA/vSGwPJmXlUfu+NpZl5V1cU+ymTmSUS8G/07gDU3N4f5uWe6MTxDQ54ELqrqzfKgqmaZebrHMtMBdQDgO9g6J5CZk4g42XBp/tBN/TFlMvPd8mkBgKfXNzE8jYj5hvP38XAPflCZzJxGxGxAHQH4TvpC4Ci+Te6umkfEZMcyJ1X1ubeGsZhjyMy7zLz78uXLkCIADHCQJaLdsNDt0M9X1XVVva2qt8fHf9oOG4CRhoTA0YZzDz0F9Jbp5gyiqjYNGQHwhPpWB93F5hv+UUQ8NJTTV+Y84uvS0K8y8+8RMa+q6546AbAnW0OgquaZOcvMyVrPfVJVG4dzBpT5U7nMvBzyAhoA+zVkOOgyut57xNce/O3K8TQzb5bDPEPKAPBj6A2BbnhmnpmnmfkuIk7X3vydRsRprMwDDCgTEYsJ4sy86v77asBLaADs0aC9g7aN03dDPK8eU2at7G1EPLgFBQDfj11EARomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGjbon5cE+t3cHO5nn50d7mfzvHkSAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGjYoH9jODPPI+K+O5xW1YddymTmJCL+2h1OIuLXiLisqtnQigOwu94QWN7Mq+pjdzzNzKuqutihzGVEvK+qeXf9NCI+RcSr3X47ADzGkOGgi+XNPCKi662f7ljmbfdraRYRk+4JAYAnsjUEupvyyYZL8673PqpMVb2pqtuVa9OImC+fDAB4Gn1PAtOI2HRjvu+u7avM+4j4W09dANizvjmBo/g2ubtqHosJ3Z3KZOa7iPgtFpPCtxvKLD93HhHnERG//PJLT5UBGOqgS0S7eYP3EXHWBcJDn7uuqrdV9fb4+PjpKgjwwg1ZInq04VzfBO7gMt08wEVm/jszZ1X1eUCd+MHd3By6BsAQfU8Cd7H55n0UEQ/drLeWycxJN7yzbhYR/9VTHwD2aGsIdL302Yalm5OHxvAHlDmNxXsC6yYR8b/Dqg3APgyZE7iMblI2IiIzTyLiduV4mpk3azf9bWVuYzEPECvXp7F4Urh+7G8AgPF65wSq6jozz7s1/pNYbAGx+rbwNBa9+6PoloZuK1NV88y8zcy/d+XnEfEmIt54TwDgaQ3aO6iqHuyhd0M8f9ruoafMLCJ69x8C4PuyiyhAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMOEAEDDhABAw4QAQMN+GvKhzDyPiPvucFpVH3Yt012PiPg1IiYR8b6q5oNqDcBe9IbA8mZeVR+742lmXlXVxdgymXleVdcrn38XEZ9iEQgAPJEhw0EXy5t5RERVzSLidGyZzJzG2s2+++xRFwYAPJGtIZCZk4g42XBpnpkbg2BgmfMN1+8j4mhbfQDYr74ngWlEbBqnv++uPbpMVc2q6tUD5e566gPAHvWFwFF8m9xdNY/FZO5eynRzCLdV9fmh65l5l5l3X7586akyAEMdfIloN0dwUVW/PfSZqrquqrdV9fb4+PgJawfwsg0JgU3j9A89BYwpcxkRfxlQDwD2rC8E7mLzzfsoIjYO3TymTGZehvcDAA5mawh0N+dZt+Jn1aSqbncp080DXHXLR5fn+paeArBHQ4aDLmNlSWdmnkTE6s18mpk3azf9vjKnEXG3DIDMnAgAgKfX+8ZwVV13q3NOYzHMM117W3gaixfBjqJbGrqtTDcR/M/uv9d/3KalowB8J4P2Dlrd4mHDtdvYcPN+qEzX+//T3R+Ap3fwJaIAHI4QAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGiYEABo2KA3hoEf283NYX7u2dlhfi7740kAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYUIAoGFCAKBhQgCgYT8N+VBmnkfEfXc4raoPu5bJzElE/DUifquqs+FVBmBfekNgeTOvqo/d8TQzr6rqYmyZzDyJiGksQmK6h98HACMMeRK4qKo3y4OqmmXm6S5lqupzRHzuwgCAA9kaAt2QzaYb9TwzT6vqdh9l+H5ubg5dA+BH1jcxPI2I+Ybz24ZxxpQB4AD6QuAovk3urppHxGSPZbbKzPPMvMvMuy9fvoz5CgA2eBZLRKvquqreVtXb4+PjQ1cH4MUYEgJHG8719ejHlAHgifWFwF1svnkfRcTnPZYB4AC2hkBVzSNi1q34WTV5aJXPmDIAHMaQ4aDLiDhfHnRr+29XjqeZebN2099aZsWmYSMAnkjvy2JVdd2tzjmNxTDPdO1t4WlEnMbihj4fUiYzpxHxLiJ+i4iTzLyMiN+r6npfvzEA+g3aO2jbzbkb4nn1yDKziPjQ/QLgQJ7FElEAvg8hANAwIQDQMCEA0DAhANAwIQDQMCEA0DAhANAwIQDQMCEA0DAhANAwIQDQMCEA0DAhANAwIQDQMCEA0DAhANAwIQDQMCEA0DAhANAwIQDQMCEA0DAhANAwIQDQMCEA0LCfDl0B4Pm6uTnMzz07O8zPfYk8CQA0TAgANEwIADRMCAA0TAgANEwIADRMCAA0TAgANMzLYk/gUC/UAPTxJADQMCEA0DAhANCwQXMCmXkeEffd4bSqPuxaZsx3ArBfvSGwvFlX1cfueJqZV1V1MbbMmO8EYP+GPAlcVNWb5UFVzTLzdMcyY74TICIOu+LupW1jvXVOIDMnEXGy4dL8oZt2X5kx3wnA99H3JDCNiPmG8/fdtTFl7kd8515Yrw/wR30hcBTfJm9XzSNiMrLMo7+zm0M47w7/LzP/e0udifg5Iv516Eq8INpz/7Tpfg1pz//cdPJZvDFcVdcRcX3oejwXmXlXVW8PXY+XQnvunzbdr13ac8h7Akcbzj30FDC0zJjvBGDP+kLgLjbfnI8i4vPIMmO+E4DvYGsIVNU8Imbdip5Vk6q6HVNmzHfyaIbO9kt77p823a/R7TlkOOgyvk3KRmaeRMTtyvE0M2/Wbupbywy4zg66ORT2RHvunzbdr13aM6uq/0OL1TmzWAzj/GGLh25t/01EvKmq2ZAyQ64D8P0NCgEAXqZnsUQUDsEmh7vr2jAi4tdYPPW/7+YFV69r4xEy86aqztbOPbo9hQBsYJPD3WXm+epYdWa+i4hPsQgEbbyDbh713dq5Ue1pOOgF0vvaXWZ+Wt3ksDv3e1X9eqg6PSeZOY3FRpHv187/OyL+VlUftfF4XaDeVFWunBvVnv5RmRdm2fvqfr2PiH/Govf19Xp0vYWux/AxM68OVd8fkU0O9+Z8w7n7iDjSxuNl5rtlb3/l3Oj2FAIvSNf7+kPqd39YjrqeQ8Sid/Zx5fosIvyl+6MxGyeyoqpmVfVqw6VpLF4Y1cYjdH/HZxsujW5PIfDy6H3tbszGifTonkJvq+pzaOOxTrr2Wze6PU0MvyBdr17vix/OyhzBm94Ps1HXUdv7S7WeBF44va/RbHK4X5cR8Ze1c9p4oOWODKsLPDYY1Z6eBF4wva/RbHK4R5l5GWsr1EIbP9Z5xNeloV9l5t9j0Yn7R4xsTyHwg+oe/d73fnDhYnXLjhV6XyNU1TwzZ5k5Wbtx2eTwkbon0au1LWVOq+pWGw+3aRl3Zl6ubeEzqj2FwA+q+x83+i+D3tfOlpscfoiwyeEYXUfmbhkA3ZDG6j98oo33a1R7elnsBVqZB9jU+/o9Fpv9rb485gWdDWxyOF43FPn7A5dfLf/8aePH68L1LBY3/OtYvDR22117dHsKgRem+wNyv1xGttr76kLgPBaPiKu9hQuv6kObhMALovcFPJYQAGiY9wQAGiYEABomBAAaJgQAGiYEABomBAAaJgQAGvb/cVUMVpf2YsAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD+CAYAAADcWrmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYbUlEQVR4nO3d23IU57nG8ecD7MSbOLIK2bHjYjnjpFIp5yAFOknlKBW5KheAs64g4g7CNeA7QJew4DRVqVhxOXE2LkfgLQYDEmiHBAIxSEJog/Sug68HDcNI3TPT019v/r+qKWs2Pby0UT/d366dmQkAUE2HQhcAAAiHEACACiMEAKDCCAEAqDBCAAAq7EjoAjp19OhRe/vtt0OXAQCFcuHChbtmNtT6euFC4O2339bExEToMgCgUJxz0+1epzkIACqMEACACiMEAKDCCAEAqDBCAAAqjBAAgAojBACgwggBAKgwQgBoY3FR+u9/pfv3Q1cC9FfhZgwD/TY9LX32mf95bk76zW+k114LWxPQL1wJAE0uXdoLAEl6/Fj65BNpZSVcTUA/EQJA5M4d6dtvn319d1f66qvs6wGyQAgAkXYB0LCwIC0vZ1cLkBVCAJB07560tHTwZy5fzqYWIEuEACDpypX4z9y65cMCKBNCAJX38KE/wCfx5ZeSWX/rAbJECKDybt5M/tl796TZ2b6VAmSOEEClmfl5AZ344gtpfb0/9QBZIwRQaXfv+uagTmxuSn/5i/SPf/ifgSJLNGPYOTcqqTFArmZmH6S5jXPunJm9n6QWIE3dNu3s7Ei3b0vj435G8cDA0+8dOiQ5l06NQD/FhkDjYG5m56PnNefcWTM7lcY2zrnjkk52/TcAumSWvEN4P+vr0kcfST/7mfTCC37NocVF//PRo9K770ovv5xOvUA/JLkSOGVmJxpPzGzKOTeS4ja1BDUAqVtelh496v17dnaeHWK6vi7NzPg/4733pCOs0oWcOrBPwDk3IOl4m7fq+x3UO9nGOXeycbUAZO3Gjf7/GWtrB89EBkKL6xiuSaq3eX1Z+5/BJ9rGOVeTNJWgRiB1m5udjwrq1nffSfPz2fxZQKfiQmBQe527zeqSBtq83sk2x83sYmyF8n0MzrkJ59zEUtzcfiCBqSm/MFxWJiakra3s/jwgqSBDRKNmofGknzezMTMbNrPhoaGhPlaGKtjd9SGQpa0tViJFPiUJgcE2r+13FRC7TdRnIDNr12QE9N3UVJjJXjdu+HkJQJ7EjVmYUPsD/qCk/Zpy4rYZlZ4MDX3COfcnSXUzG4upCeja3bvS11+H+/O/+kr67W+ZQ4D8ODAEzKzunJtyzg20nLkPmFnb5pwE2zyznXPuTJIJaEC37t+Xrl71wzZDunfPBxGtmsiLJM1BZxSdvUtPzuDHm57XnHPnGs08SbYBsrKy4tf6+dvfwgdAw+Rk6AqAPbFTWMxsLBqdMyLfzFNrmflbkzQi39xTT7iNpCcdxO9HP5+VdG6/KwygU/Pz0qefZjsKKIlbt/y9i5lAhjxwVrDF0YeHh21iYiJ0Gcixx4/9DN483wns17+W3nordBWoEufcBTMbbn2dVURRKru70scf5zsApPw0TQGEAEplasp3AufdrVvSxkboKgBCACWytSVduhS6imTMslm7CIhDCKA0rl4t1tIMk5P567RG9RACKIXtbenatdBVdObRI/oGEB4hgFKYm/OjgormyhXfNASEQgigFObmQlfQndVVf+MZIBRCAIW3teXv91tUWd3XAGiHEEDhLSwUu0llbq7Y9aPYCAEUXtHv2rW5yRLTCIcQQKHt7kp37oSuone3boWuAFVFCKDQ6nU/PLToFhZCV4CqIgRQaGW4CpD8KKG1tdBVoIoIARRamc6gaRJCCIQACqtsHaqEAEIgBFBYi4uhK0jX0pIPNiBLhAAKq4xnzkWe9IZiIgRQSLu75TxglqWjG8VBCKCQlpfLMTS01e3bzB5GtggBFFKZRgU1W1/3w0WBrBACKByz4q4amkRZAw75RAigcOr1ck+sKtuoJ+QbIYDCKfqCcXGWlsrZ34F8IgRQOGUcGtrMzAcBkAVCAIWyuio9eBC6iv4r4/BX5BMhgEKpyo3ZuRJAVggBFIZZdW7F+OABS0ggG4QACuP2benhw9BVZKdMi+MhvwgBFMb166EryBYhgCwQAiiEpaXqTaKiXwBZIASQe2bSF1+EriJ79+/TL4D+IwSQe7OzfpZwFXE1gH47kuRDzrlRScvR05qZfdDLNs65AUl/iJ4OSHpH0hkzm0paOKrj2rXQFYRz+7b01luhq0CZxYZA42BuZuej5zXn3FkzO9XDNmcknTazevT+iKQLkl7t7a+Dslld9ctGVxVXAui3JM1BpxoHc0mKztZHetxmOHo0TEkaiK4QgCdmZ0NXENbqKv0C6K8DQyA6KB9v81Y9OnvvahszO2Fm403v1STVG1cGQEPVQ0CS7t0LXQHKLO5KoCap3YF5OXovrW1OS/pjTC2omHpdWlkJXUV4VW4OQ//F9QkMaq9zt1ldvkO3p22ccyclvSffKTzeZpvG50YljUrSsWPHYkpGWVRlnaA4TBpDPwUdIhr1G5yW9H4UCPt9bszMhs1seGhoKLsCEczubnXWCYqzvOz3B9APSYaIDrZ5La4DN/E2UT/AKefcfefclJldTFATSm5+XtrYCF1FPuzs+KaxwXa/VUCP4q4EJtT+4D0oab+D9YHbOOcGouadVlOS/jemHlRE1dYJikPnMPrlwBCIztKn2gzdHNivDT/BNiPy8wRaDUjinzq0skI7eCtCAP2SpE/gjKJOWUlyzh2XNN70vOacO9dy0D9om3H5fgA1vV+Tv1IY6/QvgPKhL+BZhAD6JbZPwMzGnHOj0Rj/AfklIJpnC9fkz+4HFQ0NPWgbM6s758adc3+Ktq9LOiHpBPMEsLsr3bwZuor8WV/3jxdfDF0JyibR2kFmtu8ZetTE88xyDzHbTEmKXX8I1bOwQIfwfu7elRghjbSxiihyZYolBPd1507oClBGhAByY3PTr5qJ9lhMDv1ACCA35uf9DWTQ3tqa7xcA0kQIIDdu3AhdQf5xpYS0EQLIhZUVFkpLomr3WUb/EQLIBRaLS2ZhQdreDl0FyoQQQC7Mz4euoBh2d6XFxdBVoEwIAQT38CH3DejErVuhK0CZEAIIjs7OziwssLQ00kMIIDjGv3dme5sF9pAeQgBBmTETths0CSEthACCeviQtYK6MTfHxDqkgxBAUDRrdOfRI5aXRjoIAQRFf0D35uZCV4AyIAQQFCHQPZqEkAZCAMGsr/s+AXTn0SOW2kDvCAEEw6ig3rHcBnpFCCAYJon1bmZG2tkJXQWKjBBAEKyBk46tLWl2NnQVKDJCAEEsLfkDGHp37RodxOgeIYAgaMtOT70u3b8fugoUFSGAzD1+zBj3tE1Nha4ARUUIIHPz8z4IkJ6ZGW42g+4QAsjczZuhKyifnR3u0YzuEALI1MOHzA/ol8lJOojROUIAmaJDuH/W1uggRucIAWSKMe39NT0dugIUDSGAzKysSA8ehK6i3KanmUGMzhACyAx3w+q/7W0/+gpIihBAZjg4ZYMmN3SCEEAmNjZY9jgri4vMGUByhAAywYqh2dndpekNyRECyAQrhmaLEEBSR5J8yDk3KqlxMV8zsw963SZ6X5LekTQg6bSZ1RNVjUIx40oga4uL/orgEKd5iBEbAo2DuZmdj57XnHNnzexUt9s450bNbKzp8yclXZAPBJTM3bvS5mboKqrl8WPfB3P0aOhKkHdJzhNONQ7mkmRmU5JGut3GOVdTy8E++uxgFAYoGUYFhcHVF5I4MASccwOSjrd5q+6caxsECbcZbfP+sqTBg+pB8ZjRPh0K+x1JxF0J1CS1a6dfjt7reBszmzKzV/fZbiKmHhTMgwd+0Thkr16XHj0KXQXyLi4EBrXXudusLt+Zm8o2UR/CuJld3O9959yEc25iaWkppmTkCWejYdEUhzjBxw5EfQSnzOy9/T5jZmNmNmxmw0NDQxlWh14RAmGxaiviJAmBdu30+10FdLPNGUm/S1AHCubRI5Y2Du3ePZrjcLC4EJhQ+4P3oKS2TTedbOOcOyPmB5QW9xHOB5aXxkEODIHo4DwVjfhpNmBm471sE/UDnI2GjzZeixt6igLh4JMPN29yxzHsL0lz0Bk1Del0zh2X1HwwrznnzrUc9OO2GZE00QgA59wAAVAu9TpNQXnx8KFvFgLaiZ0xbGZj0eicEflmnlrLbOGa/ESwQUVDQw/aJuoI/jD6ufWPazd0FAXEVUC+TE8zexjtJVo7qHmJhzbvjavNwXu/baKz/2eO/igPM0al5M3srPSrX0mHD4euBHkTfIgoyufePX//AOTH9jYd9WiPEEDqmKCUT9evh64AeUQIIHVMEMun5WU6iPEsQgCpWlvzD+TT1auhK0DeEAJI1cJC6ApwkPl5aX09dBXIE0IAqeI2kvlmJl27FroK5AkhgNTs7Eh37oSuAnFu3PCjhQCJEECK7t7197VFvm1vS1euhK4CeUEIIDXczrA4rl6lbwAeIYDU0B9QHLu70uefs7AcCAGkZHPT30oSxXHrFkEAQgApoUO4mCYnpU8+oaO4yggBpIKmoOK6fVu6cCF0FQiFEEDPzAiBopud9aO7UD2EAHpWr7NqaBl8+23oChACIYCesVREOdy+zd3gqogQQM9YOro8Ll8OXQGyRgigJw8f+uYglMP8vLS6GroKZIkQQE9mZ0NXgLRNTYWuAFkiBNA1M24oX0azs0wgqxJCAF1bWfEPlMujRwwXrRJCAF3jxuXlNTMTugJkhRBA1xgVVF5zcywLXhWEALry8CELxpXZ1hbzP6qCEEBXuAooP0Z+VQMhgK5wgCi/W7dYXbQKCAF0bG1NWl4OXQX6bWfHBwHKjRBAxxg5Uh1c8ZUfIYCOcWCojsVF30mM8iIE0JHVVSaIVYkZo4TKjhBAR5ggVj38Py83QgCJmdEfUEW3b/tOYpTTkSQfcs6NSmqMB6mZ2Qe9buOcG5D0B0nvmdn7yUtGKPfv0xRURTs70p070htvhK4E/RAbAo2DuZmdj57XnHNnzexUt9s4545LqsmHRC2FvwcycONG6AoQyuIiIVBWSa4ETpnZicYTM5tyzo30so2ZXZR0MQoDFMD2Nk1BVba4GLoC9MuBfQJRk027A3V9vyDoZhvk3/S09Phx6CoQytqaf6B84jqGa5La3TzwoGacbrZBju3uSlevhq4CoXE1UE5xITCovc7dZnVJAylucyDn3KhzbsI5N7G0tNTNV6AHMzN+1VBUGyFQToUYImpmY2Y2bGbDQ0NDocuplMePpW++CV0F8uDOHYaKllGSEBhs81rcGX032yCHrl3ztxsEdnake/dCV4G0xYXAhNofvAclXUxxG+TQ1pb03Xehq0Ce0CRUPgeGgJnVJU1FI36aDZjZeFrbIJ+uXmU9eTxtft7PHEd5JGkOOiNptPEkGts/3vS85pw713LQP3CbJu2ajZADGxuMCMKz1taYNV42sZPFzGwsGp0zIt/MU2uZLVyTNCJ/QK8n2cY5V5N0UtJ7ko47585ImjSzsbT+YujN5ct0AqK92Vnphz8MXQXS4qxg13bDw8M2MTERuoxS29iQ/vxnPz8AaPXyy9Lvfy85F7oSdMI5d8HMhltfL8QQUWTr8mUCAPtbW5Pq7aaDopAIATxlZUWanAxdBfJuejp0BUgLIYCnfP45oz8Qb3aWq8WyIATwxJ07/gHE2djwN5tB8RECkOTP/r/+OnQVKBKaDcuBEIAkfx/Z5XbL/gH7WFhgzkAZEALQ7i5XAegOy4oUHyEA3bzJUtHozswMCwwWHSFQcdvb0qVLoatAUe3u0jdQdIRAxV265Ed6AN26fp2FBouMEKiw1VX/Cwz0Ynubq4EiIwQq7KuvmBiGdFy96u9Ch+IhBCpqcVG6dSt0FSiLzU2uBoqKEKig7W3pwoXQVaBsLl/2YYBiIQQq6PPPpfX10FWgbLa3pS++CF0FOkUIVMzNm6wAif6ZmaGZsWgIgQq5c4dmIPTfhQvS1lboKpAUIVARq6vSv//N8r/ov40N6bPPGHlWFIRABWxsSP/8JxN6kJ2FBT9sFPlHCJTc2pr00Uf+v0CWvv5auns3dBWIQwiUlJl07Zr04YcsDocwzHwTJCcg+UYIlND6uvT3v/vhesziREibm9Inn7A+VZ4RAiVz964/+19aCl0J4K2tSf/6FyOG8ooQKJHFRX8FwC8b8mZ5mZOTvCIESmJuzo8AYggo8mp9Xfr4Y+nTT7kRTZ4QAiVw44b0n/8wLhvFMDsr/fWv/soV4R0JXQC6ZyZduSJ9803oSoDObG35DuMf/Uiq1aTXX5eOcDQKgt1eUGbSl1/6YaBAUS0u+sfhw9Jbb/nH66/758gGIVBAZtLEhF8MDiiDnR2/sOH0tPS970lvv+0fr7wSurLyIwQK6MsvCQCU1+am9N13/vHGG9KPfyy9+aYPB6SPECiYy5dpAkJ1LCz4x6FD0k9+Iv3859JLL4WuqlwIgYIw8wFw6VLoSoDs7e7621dOTvp+g1/8QhoYCF1VORACBbCxIV28KM3Ph64ECG9uzj+GhqR33/X/RfcShYBzblTScvS0ZmYf9LpNN99ZNZub0vXrfkle1gACnra05CefHTsm/fKXNBN1KzYEGgdrMzsfPa85586a2alut+nmO8vCzM+cXFnxB/ntbX+pe/iw9Pzz/udHj/w0+8VFZgADcWZm/AS0N9+UfvpTf2XgXOiqiiPJlcApMzvReGJmU865kR636eY7C2ttzV++3r4t3bvnh8MBSI+Zby6dn5e+//29UUWvvcacgzgHhoBzbkDS8TZv1Z1zI2Y23uk2kiY6/c6smfkz8EOHujuj2NqS7t/3K3ouLPifAWRjY8MvpXLjhv8dfv11/3jlFd9k9NJLXCk0i7sSqEmqt3l9OXqvm22Wu/jOnjUO7JubvrllY8M/trf9mXnzo7kJ5vBhP539uef2mmyee87/9/nn/T+ynZ29711dfXpxrB/8wD86qbOfslhfqAx/Rhn+Dln8GUX4Ozx+vHeVIPnf51de2Xu89JL/nd7c3DsGmPnf7cOH9/7b+Ln1ONH8vln740nzMcW5vc83f2/zdzdqaH2/H+EVFwKD2uu8bVaXtN8ArbhtOv7OqA9hVJKOHTt2cMX7aOz4F1/0DwBo9sILyT535Mj+6xw55w/Wzz3XXQ0HfXe/FGIVUTMbM7NhMxseYjwYAKQmSQgMtnktbppG3DbdfCcAIGVxITCh9gfnQUkXu9ymm+8EAPTBgSFgZnVJU9GIn2YD+43iidumm+8EAPRHkuagM4o6ZSXJOXdc0njT85pz7lzLQf3AbRK8DwDIQGw/tJmNOedGozH+A/JLPDTP7K1JGpFvzqkn2SbBdwIAMpBoMJKZjR3w3rikVzvZJsn7AID+K8QQUQBAfxACAFBhzrKY950i59ySpOnQdRTIUUl3QxdRAuzHdLAf09HNfvwfM3tmtm3hQgCdcc5NmNlw6DqKjv2YDvZjOtLcjzQHAUCFEQIAUGGEQPkxFDcd7Md0sB/Tkdp+pE8AACqMKwEAqDBCAAAqLON72ADFEt3VrnEnvJqZfRCynqKI9pskvSO/PtjpaAXh5vfZrx1wzp0zs/dbXut5PxICwD4av2Bmdj56XnPOnWWxw4M550ab1wZzzp2UdEE+ENivXYhWWj7Z8loq+5GO4RLjbKw3zrkLZnai5bVJM3snVE1555yrSTplZqdbXr8v6Y9mdp792rkoSM+ZmWt6LZX9SJ9ASTXOxqLHaUkfyp+NPXlf0VlEdCZx3jl3NlS9eRPdH+N4m7fq0RLo2N9om9eWJQ2yXzvnnDvZONtvei21/UgIlFB0NvbU2UD0j2gwOqOQ/Nna+ab3p+TvCwGvpuj+GC2Wo/fQhplNmdkzS8vL77MJsV87Ev0uT7V5K7X9SAiUF2djvRnUXlNZs7ra3yMb+4iuOsfN7KLYr506Hu23VqntRzqGSyg6q+dsDME19RGciP0wnhKdkPX9trtcCVQEZ2NdGWzzGvunM2ck/a7lNfZrjMY925sHcrSRyn7kSqACOBvryoTa/0INSmp3eY4WzrkzahmRJvZrUqPSk6GhTzjn/iR/svZ/Smk/EgI5F10Sno79oHcqagpqxdlYh8ys7pybcs4NtBzEBqL7auMA0ZXn2eZ/j865ETMbZ7/Gazdc2zl3pvn1tPYjIZBz0f/Qrn85OBvryRn5M7IPpCdnZRyoYkQnLhONAIiaNppvgMJ+TUcq+5HJYiXW1A/Q7mxsUtKJlsljTNhpEe3DKfnQZEJdjKjpcXKft19t/HtjvyYXher78gf8MflJY+PRez3vR0KgpKJ/OMuN4WXNZ2NRCIzKXzo2n0WcYuo+UC2EQAlxNgYgKUIAACqMeQIAUGGEAABUGCEAABVGCABAhRECAFBhhAAAVBghAAAV9v/fp5hvzHjfugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "def compute_density(Z,A,show_graphs=False,Z_test=[]):\n",
    "\n",
    "    bandwidth = np.sqrt(np.median(np.abs(Z)))\n",
    "\n",
    "    kde_0 = KernelDensity(kernel='linear', bandwidth=bandwidth).fit(Z[A==0][:, np.newaxis])\n",
    "    kde_1 = KernelDensity(kernel='linear', bandwidth=bandwidth).fit(Z[A==1][:, np.newaxis])\n",
    "    \n",
    "    if show_graphs:\n",
    "        plt.clf()\n",
    "        plt.hist(Z[A==0], fc='#AAAAFF', density=True)\n",
    "        plt.show()\n",
    "\n",
    "        X_plot = np.linspace(min(Z[A==0])-1, max(Z[A==0])+1, 1000)[:, np.newaxis]\n",
    "        log_dens = kde_0.score_samples(X_plot)\n",
    "\n",
    "        plt.clf()\n",
    "        plt.fill(X_plot[:, 0], np.exp(log_dens), fc='#AAAAFF')\n",
    "        plt.show()\n",
    "\n",
    "    log_dens_0 = np.exp(np.squeeze(kde_0.score_samples(Z[:, np.newaxis])))\n",
    "    log_dens_1 = np.exp(np.squeeze(kde_1.score_samples(Z[:, np.newaxis])))\n",
    "    p_0 = np.sum(A==0) / A.shape[0]\n",
    "    p_1 = 1 - p_0\n",
    "\n",
    "    # p(A=1|y) = p(y|A=1)p(A=1) / (p(y|A=1)p(A=1) + p(y|A=0)p(A=0))\n",
    "    p_success = (log_dens_1*p_1) / (log_dens_1*p_1 + log_dens_0*p_0)\n",
    "\n",
    "    p_success_test = []\n",
    "    if len(Z_test) > 0:\n",
    "        log_dens_0_test = np.exp(np.squeeze(kde_0.score_samples(Y_test[:, np.newaxis])))\n",
    "        log_dens_1_test = np.exp(np.squeeze(kde_1.score_samples(Y_test[:, np.newaxis])))\n",
    "        p_success_test = (log_dens_1_test*p_1) / (log_dens_1_test*p_1 + log_dens_0_test*p_0)\n",
    "    \n",
    "    return p_success, p_success_test\n",
    "\n",
    "p_success, p_success_test = compute_density(Y,A,True,Y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implementation of a fairness-unaware baseline algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Baseline (All): Test Error = 2.4019277050431986\n",
      "Baseline (A=0): Test Error = 1.7947208535163621\n",
      "Baseline (A=1): Test Error = 5.481289646272755\n"
     ]
    }
   ],
   "source": [
    "bhat, res, rank, s = np.linalg.lstsq(X, Y)\n",
    "\n",
    "Yhat_out_cal = np.dot(X_cal,bhat)\n",
    "Yhat_out_test = np.dot(X_test,bhat)\n",
    "\n",
    "print(\"Baseline (All): Test Error = \" + str(np.sqrt(np.mean((Yhat_out_test - Y_test)**2))))\n",
    "print(\"Baseline (A=0): Test Error = \" + str(np.sqrt(np.mean((Yhat_out_test[A_test==0] - Y_test[A_test==0])**2))))\n",
    "print(\"Baseline (A=1): Test Error = \" + str(np.sqrt(np.mean((Yhat_out_test[A_test==1] - Y_test[A_test==1])**2))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAElCAYAAABAu+8OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXjU1b3H8ffJwi6ZBJDFJELQuABVIaBWqbKoPBerAipF8SogQSrWDTAm1gVugii2al0KVqgPVlFBUFvRC+K1SgENYButVSugMxBQEgKVLSQ594/vb5jJMGFLMr/MzPf1PD5J5jeT/BKmn57zPZux1qKUUm5KcPsGlFJKg0gp5ToNIqWU6zSIlFKu0yBSSrlOg0gp5ToNIuU6Y8wyY4zH7ftQ7kly+wZUfDPGXA0MBtKAigj+3Fyg3Pkyy1r7SKR+tjqUtoiU20Y6H7Mi9QP9IWStXWitXQgsNMbMjtTPV4fSIFKuMcZMBe5xvoxk12yCE0AAWGs3IK0y5RINIuUKY0wWUOGEAESoReTUonqHuVRhjNEwcokGkXLLBGvtHOfzCqB7hH5uFuFrUeVEsHuoatNitYo4p+XxSdBDG5Bi9ZFed9R1HGvthDoupREoUgerILLdQxVEg0i54ZqQoDiq1shhwkVFOQ0iFVHOiFWaMWZm0MOR7hKFa31pa8hFGkQqYpxCcXdr7TUhj5cBM8O/qsEVEz500oB1EboHFUKDSEVSrrX2njCPH9VERmPMa0f5c8rr6sZZayuMMRuMMR5rbfDP9Vhrlx/l91cNTINIRYQzg3pDHZc3OM/JChrOP0RoS6oeZgK5wCPOz+0NaAi5SIfvVaMyxmQ5LZnXgENaKU7NyN9Kmh2JuTzOtIEKY8xg/xITLYS7y+ie1Uopt2mLSCnlOg0ipZTrNIiUUq7TIFJKuU6H7x3t27e3Xbt2dfs2lIo5lZWVbNu2je+//367tbZDuOdoEDm6du1KcXGx27ehVEzx+Xzk5+cDMH/+/G/rep52zZRSjSI4hAoLCw/7XA0ipVSDKysrqxVCGRkZh32+ds2UUg0uNTWVQYMGMXDgwCOGEGgQKaUakM/nIykpiU6dOnHjjTce9eu0a6aUahD+mtCsWbM41qVjGkRKqXoLLkzffvvtGGOO6fUaREqpegkdHTuamlAoDSKlVL388Y9/BI4/hECL1UqperrjjjvYuXMnJ5100nF/D20RKaWOmdfr5bHHHmP//v20adOmXiEE2iJSSh0jr9dLQUEBAOXl5XTu3Lne31NbREqpoxYcQoWFhQ0SQqBBpJQ6SqEhdLyF6XC0a6aUOio1NTWkpKQwderUBg0h0CBSSh1BRUUFKSkpnHzyyTzxxBMkJDR8R0q7ZkqpOnm9Xn71q1+xaNEigEYJIdAgUkrVIbgmdN555zXqz9IgUkodIjiEioqKSE9Pb9Sfp0GklKpl//793H///UBkQgi0WK2UCtG8eXNuvvlmTj755IiEEGgQKaUcPp+P0tJS+vbtywUXXBDRn61BpJQ6uJVHcnIyZ511Fs2aNYvoz9cakVJxLng/oYceeijiIQQaRErFteAQilRhOhwNIqXi2AcffAC4G0IA5lg3uY5VOTk5Vk96VfHCWosxBmstO3bsIC0trdF/pjFmrbU2J9w1bREpFWd8Ph933303W7ZswRgTkRA6Eh01UyqOBNeEqqurXb6bAG0RKRUnGuK0jcaiQaRUHNiyZUuTDSHQrplSccHj8XDmmWdy/fXXN7kQAg0ipWLali1bSEtLo1WrVuTl5bl9O3XSrplSMcrn85GXl8eTTz7p9q0ckQaRUjEouDA9atQol+/myDSIlIoxTXl0rC4aRErFEGstjz76KBA9IQRarFYqphhjmDJlCtbaqAkh0CBSKib4fD4++ugjRo4cSfqOHbB4MZSWQufOMGwY9Orl9i0elgaRUtGkpOSQkPGlph6sCV1WXk7qX/4Cv/gFDBkCX34JeXkwYgSMHevyzdct6oLIGJPrfNod8AD3WGsrQq6XO19mWWsfifAtKtU45s6FRYtqhYwvL4/8Vq2gZUuKbrqJ1Mceg/nzwb+QtV8/GDoURo+Gvn2bbMsoqoLIGJNrrZ0T9PXVwFoklA6GkLV2ofN1ljFmtrV2gis3rFRDKSmREAoKGV+XLuS/+y7s3EnRxImkf/SRhFToavq0NBg1SlpSTTSIombUzBiThRM4fk7gpDmBBDDBH0LO9Q3A4MjdpVKNZPHiQ0LG5/OR1KwZRUOGSAiVlsJpp4V/fXY2bN0aoZs9dlETRI7cMI+VI2HkAXqHuV5hjNEwUtEtKGT2798PyOmrzz77LOk5ORIynTtLTSicr76CTp0idbfHLGqCyFq7wVqbGuZSFlDsfKwIc73cuaZU9HJCxuv1MnHiRFavXg3IGWQHQ2bYMFiwAMrLa7+2vBxeflmuN1FRVSMK5dSElltr1zmtnvIwT6tAitpKRa9hw/Dm5VHw7ruQkBDYX9ofMjNnSv1nxAgpTI8aJd2xr76S61df3WTrQxDFQeTUjCZYa/vU43vk4nT3MjMzG+rWlGpwXo+HglatpDA9ZAjpmzfD++8fGjJjx8ro2OLFsGqVtJT8IdWERW0QATOBQSGPhdt8t87WkDMCNwdk8/yGuzWlGs6OHTsoKCiQIXr/6NgLL9QdMr16NfngCRWVQWSMmUnI/CGkThQudNKAdRG5MaUagcfjYejQoVxwwQXSJRs40O1banBRF0ROd2q2MzTvf2ywtXa5MWaDMcYTElAea+3yyN+pUvXj9XoByMjIYOTIkS7fTeOKmlEzkMABiv0hZIzxhAzNzyRoiN8Y0xvQEFJRx+v1UlBQwKxZs4iHswejpkXkFKeXOZ+HXk4FqfkYY3KdcPIgSzx0VrWKKv4QApg8eXK493vMiZogclpBR/wXCV4ColTUcBazer/7joLycmjRgsKZM6NqK4/6iKqumVIxado0mf/zwQcs2LULjKFw504yli1z+84iJmpaREo1eWG26DjiMPq0afDUU/DQQ9CnD7d9/jnlCxfS5bLLYOHCJr1iviFpi0iphjB3ruz7060bjBkjH/Py5PG6lJTAiy/inTKFoh072NOzJy3GjKHL/Pnwzjtw4YUSbHFAW0RK1VeYLTqOah+g3/8eb2IiBR9+CElJVKxZQ6sBAwLbdjiPxwNtESlVX2G26ABq7wMUau5cvMuWUdC9OyQlUThoEF1mzQq0oLKz4ZtvmvSK+YakQaRUfR3rPkAlJXgXLaKgZ09o3pzC/fvJGDVKWlQLF0oLa9068Hqb9Ir5hqRBpFR9Hes+QIsXkzh0KO1PPJHCffvIGDJEunB/+Yt06e69Fx54QB6Lg0I1aBApVX/HsA9QWVkZdssWuuTk8Nizz5IxYgQsXSqF6Q8/lDD66COYNAnuvz/Cv4h7NIiUqq/gfYDmz4c1a+Tj6NG1tujwer3ceeed/Mla+PJLmTE9dqysoK+qksJ0ZiZcf31chRDoqJlSDaNvX1i7Fp58EoyRr4O26Di4bKOqioushfvug88/l9Dxb9tRXi7hdcstLv8ykadBpFR9BR/zc+ONUi9asAA++QR69cL73nsUPP00VFZSuGkTGdddBzfcAG+9JfOFLr9cit1RsJNiY9EgUqo+jjCHqPK773hw7Vpo1ozCbdvI6NYNZs+Gc86R7ldxsbz28sujYifFxqI1IqXq43BziC68kGYvvcTEu+6isKqKjAMH4Npr4Q9/gI0b4Ykn4NRT4X/+Bzp2lNdNmwYTJ8rHkpLI/z4u0SBSqj7qmEPk9XpZ+a9/QUYGOe3bk/Hpp/Doo9Ilu+QSmDwZfvpTmTeUkCD7Tx/rEpEYol0zperDP4eoX7+DD/kL04nl5fTp1o0WixeD/+wxv+xs2dx+1CiYNw82b4bVq6PuqOiGoi0ipeojZA7RwdGxmhqmbdxIi6wsaTX17QvTp0trZ9o0eO89mejYsaPUiW666diWiMQYbREpVR9Bc4i8Q4dS8NFHUFVF4e7dZIwcCStXgrXS+vnZz+Cf/5Sv33hDNsFftAhSU2FwHYcR+1tOMU6DSKkjCbfPENR+bNw4Pl6+HPbupfCMM2SIvlcvaf08/jh06SKF6uJi6Yr16AHr14PHA//1X4d07w5q4kdFNxQNIqUOJ3iO0JAhEhjXXy/XpkyBIUOo+de/SJg7l+HDhzNoxAg8npBTrSZPho8/hvHjpVb04INQUwOzZslM6ltukcL00KG1u2fBp7jGOA0ipeoSbo5Qy5bSQqmpgbPPxuvx8Mj8+dx98810nToVz7vvQs+egd0ZS0ulLpSfH2hZrVwp32PqVPk8io+KbigaRErVJdwcocWL4eKLYe1avLm5FHTuDEDSH/4AgwZJQPmH3keMqD2qFnoC6/z5gW5XlB4V3VA0iJQKFlwPKi6GO++sff3996GyEu+FF1Jw4AAYQ9Enn5A+eTKce64cBX3DDYGh93HjpHt3NN2uKDwquqFoECnlF1oP+v57eOQR2LdPWiwlJbB5M6W5uRSsWQNJSRRlZ5M+eLCsG/vxx0ALxz/0/vnncd/tOhomHk6RPBo5OTm2uLjY7dtQbikpke5UcD2opATuvlu6W7/9rbSUkpLY/3//x9NJSVw7eTLpr70mNaA334Rnn5Xumb9GtGePtJCeeSbQ0tq6VcLqaE74iDHGmLXW2pxw17RFpBSErwf16iWPPfEEPPggm2tqSOnRgzalpdxljMyGrq6WovPf/y4BM3my1ITy8qQ+lJkZ+F5xFjzHQoNIKZCa0JAhhz4+diw0b47vqafIT0sju6SE+156Sa79/vewYgV8+60Ums85R4rS/frJOrLzz5dWkjoiDSKl4NA1Y0FFa9+2beR37AitWnHTzp1w0kmwZAls2gRnnSUr6FeuhLIy2fDMv7fQFVcEakTqsHStmVJQe81Y0GGJvlNPJX/PHti/n6Ldu0k/91y48krZifH226GiAnbulC7Z8OFSY1q9WkbDxo8/9AQPFZa2iJSCwKTCK6+E//wHHnkE+8orPL5hA6SkUDRqFOlFRbBjByQnQ/fu0iqyVjY1y8+X73PmmbLXUK9etecJqcPSIFLxKdz6sbFjZd/prVvh+ecx69cz+Te/oerss0lPT4fKSulqvfQSPPwwXHddYLStvFwK3f5FqnG0PKMhaBCp+DNtGrz4ImRkSMsmKSkwE7qmBt+4cby3ZAk33HcfnS6/PPA6f8hkZ8uuih9+KCF27rmBeULffAPbth1ygoc6PA0iFV+mTYOnnoKHHoI+fQIb3V92GSxciC89nfx58+DAAS7v0oV2wa/96ivYsgW2b5fJjsOHS+tpwQI47zzZ4uOll6SrdsstGkLHQCc0OnRCYxwoKZFWz513yr7Qfs4xPr4+fcgvKYHqagpPPZWMc86R5Rr+51x5pRwVNGoUvPuudNVGjZLNzaZMgbZtZXLj2LHu/H5NnE5oVAqkJpSRIS2hYGlp+IYOJX/ZMgAKBw8mY+FCmSPUsaN0tV5+Gdq1k5bTW28Faj/+RaqnnCKFaQ2h46JBpGJTuGJ0aanUhMJsQvZDWhrNqqp4oFcvMm6/XXZP/PWvA3sIXX+9DOvPnw833xzodvk/rlkjyznUcdF5RCr2BM0DqnUiRmkpdO1aa4/pvXv3AnBORQXPfv217KwIEjBLlsCf/yyTFv1bvl5+efhWT5zspNhYtEWkoltoy6dHj8BmZps3yzKMjz+WNWGlpeD1ykb1o0fjGzqU+1auZEz37lw0ezbJkyYdWmAOXiMWOlTvp0P19abFaocWq6NQ8LYdp50mXa7p02WT+p/+VPaKTkiAq64KPL+8HFq3xte9O/mpqQAUbdpE+sUXSz0ouCsXbtRr7lw5iyzclh5aHzqswxWrNYgcGkRRJty2HSBdsXXrZG5Qu3ZyrPOKFRIwf/87rF2Lr2tX8rOzoaqKojVrSLdWRrwmTJDTNPxD+iNGhA8X3dLjuOiomYo9dR313K0bfPYZ7N4tW3DMmCHPa9YMXnyRXQkJ5HftCpWVFG3ZQvpVV8Err8iyjQULAiNfhzvcULf0aHBarFbRqbQUEhMPPSt+2DCpA1VUwBdfSIvpwAH4zW+gqoq2ffsy/JtvKFq9mvTqaumq3XCDdK1OPFGO+ikpiavDDZsCDSIVnUpLZRvX0JGxTz6RbTq2b4fmzWHOHJgyBd++fXyTkwPNmnFVWRnpqamyBez//q/UeTZuhBtvlDDyh092tq6ejxDtmqnos2iRHE5YUyOLUM8+O7Bh/bXXyjB769bynOJifO3bk3/GGZzQrBlPLl1KYk2NtKauvlpmSG/aBD6f7CsEgfDRIfmI0RaRii7TpkFurgRNVpbUdS65RGZLjx8P330ntaFzzoF9+/CdcAL5p58OiYncu2oViS1ayFB+jx7SksrMlNbV44/LjOmdOyV8/EPy/lNdVaPSIFLRw79gtUcPWdt16qmwa5cMu2/aJC2g8nIZHTvxRHynnEL+2WdDTQ1FX31F+plnQosWMqT/r3/JNq5ffgmnny6b42/fLoXupCRdPR9hGkQqOpSUyNYdDz0kQ+ylpdJ1+u//luL0gQOyaZm1MpHx3XdZ1LWrnDu2eTPpKSlwzTWyaLVzZylmv/SSfH3eebLv9N690mWrqpLJiTovKGK0RqSiQ/CC1ZYtJYw6dJDicnW1HN2TkCDD8B06wLZt/PJvf6PsscfotGSJHIw4ZYq81hj5mJgoLaS1a6V7d9JJEkL33+/2bxt3jhhExpiu1tpNEbgXpcIrKZH6TVmZLM/IyJDtXDt0kPPE9uyREbLWrfG2b8/cHj2466OPOOH77+n0xhuyiX1ZmTzn3/+WFlSHDvDLXwZCx9kKRJdpuONoWkQzgZGNfSMqzoVbLd+rV2AZR02NtFZ275ZaUFKShEplpbw+MRHvjBkUvPMO7N/Pzl/8ghOef14WrW7cKEP9/u08BgyQUbJFi2Q1fk2NnrzqsiMu8TDG1ADnWGv/Hplbcocu8XBRuDVjCxbIFqxr1sC998p/LVtC+/bw2mvyurQ02SkR8LZrR8FPfwpJSRSuWkVG//6y2NXnkw3NOnaULhlI8CQkyM/ZtUtW1OsyjUZX3yUe7wG/MMYYa+2nDXtrDc8YkwuUO19mWWsfcfN+1BGUlARWy/uXa/TrJ3OCzj1XhthffllqQd99J/tEt2oVWKAK+Nq0oaBfP6ispPCDD8jweGQo/913pXh9xhky1B8adp9+KtMAtCbkurBBZIxpa63dBWCtvcR5bIQxJsta+3okb/BY+EPIWrvQ+TrLGDPbWjvB5VtTdalrzVhamszxef116YKdeqqESpcuMuK1ZYsUm888k+YbN9Jp925u++wzMpKTpQCdny/dtuxsePVVmRu0ZEntsPvxRxm2LynR1pDLDhm+N8YMAg5pPllrFwEbjTGTI3Fjx2mCP4QArLUbgMEu3o86ktJSaaGEs3Mn/OMfsjFZmzYy43nbNmktJSbyQ+vWVJ9xBh0SEpi5ciUZFRUyglZZKeHSrp2ctrFrl6yuD91D6K23pPit68lcVyuIjDE3A2XW2hXhnmytXQ88Z4x52BjTNhI3eLSMMR6gd5hLFcYYDaOmyn/Uc6iSEpmk2Ly5BMuWLVIfSkqCP/0JL3D3WWfxx82b4YwzMImJMnTv16aNTIAcMUKG5YuLpfu3Zo189E9YHDRI15M1AaFds41AFlBnLchauxPIM8ZMNsYsCy1iuzjcnwVUhHm83LmmmqJhw2Sx6tCh0mLxj54tXixdJ49HWkZ79sjM59278bZpQ8G550J1NZd6vfCTn8jzqqslWKqqam9k37UrpKbK6NmqVXJt5kw9jbUJqRVE1tr3jDHdjDHDD1cLMsYMBPoBM40xy4FlgD+UJgD3NuZN1yGNQJE6WAXgCfcCp6aUC5CZmdl4d6bq5j/qefRoqf988YU8VloqrRqQYLn0Uli6FG92NgXp6WCtjI6df75MVgRpMRUVwW23ybliIF2w77+XEbNJk3SL1ybqkBqRtXZjcAgZY2Y4HwcaY541xpQBy5Hay3PATuAWYL0xphqYGpE7bwDW2jnW2hxrbU6HDh3cvp34NXYsjBsHf/2r1Iu6dJFw8p/C6hSgq2bMYFpODiQmUrh+PRmVlfKarCyZaQ2ykVmXLtKC8nfBxow5uE912O6ZFqpddzTD9xOclkMq0rp4FXjNWvte8JOMMSnAJcDDDX6XRy8tzGNhW0Oqifn8czm+5+yzpVv2ww8yWtamjdSIli4lqU8fbisvJ3XdOpkn1L69dNtSUuBvf5ORsORkGW174YXaXTCQkPKfQxZ6TbnqaILIAywEZoeGTzCndrTQGBN2wlIEFBM+dNKAdRG+F3WsSkslcBYskOH8IUOkS/bmm/hOP52vqqsZWFjITzZvlqApLpaWT3W1fH3RRbIA9oUX6p4XpFu8NllHE0RzrLW3HMP3nH28N1Mf1toKY8wGY4zHWhtctPZYa5e7cU/qGCQkyHKMFSsCdZzXX8d3443kl5aSYAznrVxJq65dZZV9UpKcNX/VVbLtxw8/wHPPSSipqHM024AcUyXPWrvxOO+lIczEKT4DGGN6I/Us1ZSVlMgWrz6frAlbtAhuvRXfeeeRv307GMP01atp1aWLBM7OnTLT+sor4Z13pJ40b54sgO3Rw+3fRh2HmDtOyKlnbUC6aUe9xEPXmjWyIy1q3bcPTj4Zli2D77/Hl5JCfp8+MjpWXEzG7t1y7vy//y2jXdbK2WUnniijYnv3SlctM1OXbDRRcXWckLV2jtv3oEIEL2odMkQmMOblBRa1zp8vOy8mJcF770FiIp926ADJyRSedRYZH38sExvLyuDMM2HyZFkAu3ixdMu8Xnj+eRkp0/Pno1LMBZFqYoIXtW7eHGgVdeokR/ycdlpg+9e77qJ6714SH3uMy9u04aKvv+aEO+6AP/1JWjzNm8ukxi+/lM3y/RMSN27UyYlRTreKVY3Lv6h1yRJpBXXrJmvA1q6VUbHKSnns8cfx7t7NpJ49+WrhQtiwgRNWrpS5PjfdBCecENjQbMEC6Z4Fb3Cvm91HNQ0i1bj8ByH6W0Vnny0bm61YAcOHS7hs3Ii3VSsK+vRhd/PmtNyzR1bF/+MfgdnVJ50kz62qkhM6BgyQ/3r3lu08dHJiVNOumWpcnTsH6kNpadIN839++unw5pt4Z8+mICdH9hP6+GMyqqth/35ZP3bVVTJKVlwskxX9R/1ceKF8//Jy6Zrp5MSopkGkGod/lOyzz6QAfdll8nhpqRSsna7U1ubNKejZU0bHvviCjPJyWReWkiJ1oWnTpIhtrYyQ3XKLBk4M0iBSDccfPu+/L4Xpm26SEa7PPpOFqL/9reyIuH27bL1hDO1OO41+iYlcuWEDGZmZsin+nj1SxN6yRYbrO3SAqVOllrR4sQZRDNIakWoYc+dKMTopSQrQd94JK1fKIYbt2snJG7t2SZfr7bfxJSSwc+xYkktLmVRZScatt0o9aNAgCbBWrWRXxuxsqSXdf7+eRR/DtEWk6i94iP6pp6T75B9eHz5ctnv96it44gk45RS8F19MwVdfcfIbbzB9wAC5vnMnTJwos6traiSwtm6FggIJNNCz6GOYtohU/QXvOx289at/ZGzFCtnqY+pUvFu2UOD1QkICuVu3SguqQwdpLSUny9Yf8+bBBx/IKFhNjZ5FHwc0iFT9BYdP8NavpaWyNYfTnfJWVFDQuTO0bUvhFVeQsWeP7Cn96qtyKsfnn0the98+KVbv3y9dOz2LPuZpEKn6Cw6fYcMCEw4TEiRIPvkEm5fHM3/+M7RoQWFhIRlt2sh8oNGjZR7Qz34GS5fKRMfx4yV0fvMbmU2tZ9HHPK0RqeMTvIg1IQGWL5d9p/1bv150kQy/l5bCJZdg/vhHJhvD3n79SP/229pbtC5efOjcoGbNpHakLaC4oEGkjl24RawffggXXwxTpgRO3qipwXfiibz97beMu/BC2v3sZzBrlpxFf+WVgZDRsIl7GkTq2BzuZNYrr4TVq2VvoZ498Q0dSv7rr0NiIsO8Xjq8+KIMwZ98srSU9GBD5dAakTo2hzuZNTdXis59+uAbN478JUvAWoouvZQOv/udTGrcuxc2bIBRo/RgQ3WQtohU3cJtZuZfohFOdjasWoXP4yH/+efhwAGKbrqJ9JEj5XrwMc8JCTo5UR2kQaTCq2szs717pch84om1d1qEgxMOd55zDi3Xr+fXycmkV1YGvmfwMc+vvy5HSSuFBpGCQ1s+PXqErwNVVMADD8gExKeekhZNXp6Mkl11FbtffpnW48bRY/16nikpIbG6WsKpY0c5s/7ll2VYPj0dZs+GBx909ddWTYfWiOKdf41Yt25yEGG3bnDvvRJIwXWgkhLZqL6wEM44Q7pXALffDk88gW/YMH6ZnMz/zpsH3bqROHQo9O8v4TVqlGyMdv31Mnt6yhSZQ6SFauXQFlE8q2sEbMUKOaoneFTLX6TOzpbHJ048eFihr3t38quqoKaGMx9+GHr2lA3Q8vLk4MObbpJlGitXyozptm1h+nTXfm3V9GiLKJ7VNQLWrRv85Ce1R7X8yzj8C0979YL778eXn09+8+ZQXU3RkCGk9+wpz/dPbLztNunqHTgA558vm92PGaOtIVWLtojiWV0jYMOGyQr64EJz586y/OKtt+Sc+mnT+HHzZgq2bwegqKyM9L59a3+fsWPlmOff/U5aRunpupOiCktbRPEseI1YsF69pA700UfSbVuzRhaePvCAbN86dy5060abceMYlZ1N0SefkN6qVd3fq39/uPZa2VNIQ0iFoUEUz4IXqAYrL5fdEWfMkP2gX3hBFp4OGABvvIFv+3a+eOYZuO02hixcSPp110nX67nnwn8v3b5DHYF2zeKZv44zerSMbGVnSw3IP8yenS1bc1gL69bB+vX4evUiPzWVVtby9D//SWLfvlK8Pu002QK2ru+lLSF1GNoiindjx0rdxt/y8Z+IAYFh/f79YdMmfEB+Whq0a0fBc8+R+OmnMsHx3ntllK158/DfS7fvUEegLSIlrZXgFkvosP60afjGjyd/6VKorqZwzBgyMjLkuaNGyXD/T34iLdNxWZgAAAx0SURBVKDQ76XUUdAgUodavBguuEBmT5eWQnExf+7fHxISKNy0iYxly2SJx+LFsqnZpk1QXS3nlCl1HLRrpg71/vvwl78EZltnZjJ+xQpmGUNG9+5SlL7hBtlPGqSwXV4uW7sqdRw0iFRtJSVyJtno0fgGDKDgjTfYMWwYiT/+yImrVsmm9snJMv/ou+/k+Q8/LMXqsjL5WqljpF0zJfwLX996CzIz8b3yCvkrV8L+/ez+9FNSL75YWkrffSctn1WrZMj+ootk8/vrrpOA0gMQ1XHQFpGqvfA1MxNfTg75rVtDeTlFGzeSfvvtEjht2kDLlrIe7bPPZM1YZmZgZEwPQFTHSVtE8S5khMz3xRfkf/01pKZS9PHHpGdmysr5Tp1g5EjZJH/MGHntxo0yW9pPD0BUx0lbRPEuZOFr65//nJP37aMoL4/0jAzIyoJnnoFJk2TJx623ymzsjh1rt350BrWqB20RxTtn4eu2bdto164dqeefz/Rhw2DqVBkVW79eWkv+GdIjRsjx0FOmwCmnyDo0nUGt6klbRPGuc2e8n3zClClTmDNnjjzmn22dmAhffy3LPIJnSF91ldSHOnXSGdSqQWiLKM55L7iAgscfh5QUrrjiisCFk06SQw7vukuOCOrRA/bsCbR+xozR4FENRoMo3gTtT+31eCj49lto0YKi8nLS33//0MWqY8cGXrNqlbSCdE8h1cA0iOJJ0Mkc1ZdeStGTT8J//kPR4MGkDxxYd9jo+jHVyDSI4kXIMH0icLvHQ5uqKtKnToWBA2sPxSsVQVqsjhfOML13927efvttAE4//XTZY1pPXVUu0yCKF6WleNu1o6CggAULFvDjjz8GrumMaOUyDaI44fV4KJg3D4DCwkLatGkTuKgzopXLtEYUB7xer4yO7dtH4ejRZMybFzjVdeBAGSHz78qolAu0RRQHvvzyS0xyMoUeDxnjx8v+QX37ysfhw2W5ho6KKRdpEMWwqqoqAAYPHswzkyaRUV4Or78OXbpAcbF8fP11OZde9xFSLtKuWYzyer1Mnz6dX/3qV/Ts2ZPWS5fK4tYLL5T/gvlHzbRVpFyiQRQrws2YTkoiJSVFrtd1qivIqNmqVZG7V6VCaNcsFgRtbOb9+c8p+Pe/4T//obBPn8BpG3Wd6go6aqZcp0EU7YJmTP8wZAgFCxZAy5YUzphBxrvvBmo/hzvVVfcRUi6Lqq6ZMSbX+bQ74AHusdZWhFz3/y8ty1r7SIRvMfKCNjZLq66mf//+DBkyRFpCwbWfI53qqvUh5aKoCSJjTK61dk7Q11cDa5FQOhhC1tqFztdZxpjZ1toJrtxwpJSW4svJoWVZGe3atWP8+PGBa6G1n7FjZdheV9KrJiYqgsgYk4UTOH7W2oXGmOeMMVc74TPBWtsn6PoGY8zgSN9rpPk8HvLnzaPLX//KjBkzMMYELoar/ehKetUERVONKDfMY+VAmjHGA/QOc70ilsPI5/OR78yYvnX06NohpLUfFUWiokVkrd0ApIa5lAUUOx8rwlwvd66F5XTncgEyMzPrf6MR5PP5yM/Ph6QkCgcPJmPqVK39qKgVFUEUjhMiy62165xWT3mYp1UgRe2wnJrTHICcnBzbKDfaSJ577jlAFrBmZGTImjGt/agoFZVB5NSMatWE4s3dd9/Nrl27SE9Plwe09qOiWMSDyGm93HOUT5/gdMtCzQQGhTyWFuZ5dbaGopHX62XJkiXccssttG3blrZt27p9S0o1iIgHkbV2ObD8eF9vjJlJyPwhpE4ULnTSgHXH+7OaEq/XS0FBAQDXXHMNnXQmtIoh0TRq5q8LzQ5uJRljBjuhtMEZPQvmcYIvqgWHUGFhoYaQijlRE0ROl67YH0LGGE/I0PxMgob4jTG9qUfLq6kIDaGDa8eUiiFRUax2itPLnM9DL6eCjIAZY3KdcPIgSzyiflb1/v37adOmDffee6+GkIpZxtqoGrVuNDk5Oba4uNjt2zho165dB4vR1dXVJCYmunxHStWPMWattTYn3LWo6ZrFE6/Xy6RJk3jzzTcBNIRUzNMgamKCa0K9e4dbtaJU7ImKGlFMCdpJkc6dZS2YMxExOISKiooCkxWVinHaIoqkoJ0UGTNGPublwdy57Nmzh1//+teAhpCKP9oiipSQs+cB6NcPhg6F0aNp1bcvN954I6eeeqqGkIo72iKKlKCdFIP59uzh786C1QEDBmgIqbikLaJICXOKhn8rj+bAM61bk+zOnSnlOg2iSPGfotGvH5SU4HvpJfK/+AISEnigRw+SkzWGVPzSIIqUYcOkMF1Rge+dd8hv1QpatqSod2/SH30UJk1y+w6Vco3WiCKlVy8491yYPp1lp5wCSUkUXXgh6e+/D/fdB6tX67HPKm5piyiCrLWYO+7gxgMH+PnevbQvKwvspJiaqsc+q7ilQRQhPp+PJ0tKmJyby4mXXkr70Cfosc8qjmnXLAL8o2Nba2rY//XX4Z+kxz6rOKZB1MgOnrYBFN56Kxlvv63HPisVQrtmjWjz5s2BEPJvavbtt3rss1IhNIgaUUpKCqeccgpjxowJbGqmxz4rdQjdGM3RkBujbd26lbS0NJo1a9Yg30+pWKAbo0WQz+dj6tSpPP30027filJRQ4OoAQUXpq+++mqX70ap6KFB1EBqjY7paRtKHRMNogZQU1PDww8/DGgIKXU8dNSsASQkJHD33XeTlJSkIaTUcdAWUT34fD4WL14MQLdu3TSElDpO2iI6TsE1oYEDB5KSkuLyHSkVvbRFdByCQ6ioqEhDSKl60iA6RqEhpHtMK1V/GkTHaOPGjSQkJGgIKdWAtEZ0lA4cOEBycjL9+/enb9++tGjRwu1bUipmaIvoKHi9XiZOnMjatWsBNISUamDaIjqC4GOgO3bs6PLdKBWbtEV0GHoWvVKRoUFUh7KyMg0hpSJEu2Z1SE1N5ZJLLtFjoJWKAA2iEF6vl6SkJDp37swNN9zg9u0oFRe0axbEXxOaNWsWunOlUpGjQeSorKw8WBO64447MMa4fEdKxQ8NIsfmzZsB3U9IKTfo5vmOdu3a2U8//VRDSKlGcrjN8zWIHMaYH4Bv3b6P49Ae2O72TbhAf+/oc7K1tkO4CxpEUc4YU1zX/8vEMv29Y4vWiJRSrtMgUkq5ToMo+s1x+wZcor93DNEakVLKddoiUkq5ToNIKeU6XfSqoooxJhcod77MstY+4ub9NBbn9wToDniAe6y1FSHXY+bvoEEU5eLpDev/Xay1C52vs4wxs621E1y+tQZljMm11s4J+vpqYC3ybxyTfwctVkexOt6wM621db5hkaCKyjess0SgT8hj3/h/31jg/BtNsNbeE/L4DmC8tXZhLP4dtEYUpZw3bK03nhM4aU4ggbyhFwZd3wAMjtxdNhxjjAfoHeZShTEmKn+nw8gN81g58m8bk38HDaLoFk9v2CygIszj5c61mGCt3WCtTQ1zKQsoJkb/DhpEUSoO37BpBGpdwSqQ2ljMcrrYy62164jRv4MGUQyJhzdsvAmqGV3i9r00Jh01ixFBb9g+R3xy9EoL81ish+pMYFDIYzH3d9AgagKcms09R3yimOAUnUPF+hu2mPD3ngasi/C9RIQxZiYh0zGI0b+DBlETYK1dDiw/3tfHwxvWWlthjNlgjPGE/J4e5+8XU5xu9uzg/9Mxxgy21i6Pxb+D1oii3GHesBXABmf0LFg0v2FnEjRSaIzpTT0CvKlyWsjF/n9TY4wnZKQz5v4OOqExijlvznKnOO2fa5MD0spyQsrjn03tvGEnROuERjgYvBuQ1l5UzxQPx6n1fVPH5VR/KyjW/g4aRFEqXt+wKjZpECmlXKc1IqWU6zSIlFKu0yBSSrlOg0gp5ToNItXk+Tf+MsZYY8wOZwJnuOctc57zmjOqqKKEjpqpqGGMeQ0YXMeuA/4Z5q/451Wp6KEtIhVNZgOhs4yBg7tTaghFKQ0iFTWcpSkVwDXBjzszxtEQil4aRCravApc6//CWdYyOHhLXBV9NIhUtHmN2t2ze3XZSvTTYrWKOs6JFsuBT4A5IdthqCikLSIVjV4Frka2xdUQigEaRCoapRHYm1vFAA0iFY16A8vcvgnVcDSIVFRxRsmyiPIdCVVtGkQq2gwGnTMUazSIVFRw1pvNRPZrxhgzU9eTxQ4dvldKuU5bREop12kQKaVcp0GklHKdBpFSynUaREop12kQKaVcp0GklHKdBpFSynUaREop12kQKaVc9/9ATsK0YEFppQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAElCAYAAABAu+8OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXxU5bnHf282IIkkGbYEEiEBBUVAMVhU1LK4FWtFNqlwUW4bwA0oyJJYRL2JRvCi3oqCJdpyrVjBYNVu0HqtpQoERAO20ELQGQiLGQJSNGR57x/POczJMFnIMidn5vf9fPhM5pwzc96Jk5/P+3uf53mV1hqEEGInEXYPgBBCKESEENuhEBFCbIdCRAixHQoRIcR2KESEENuhEBHbUUptVEol2j0OYh8UImIrSqlxAEYBcAX5volKqSyl1JvBvC8JDIWI2M1E4zEjWDdUSg2GiJ83mPcldUMhIrahlJoPYIHxNGhTM631Dq31OgD7g3VPUj8UImILSqkMAOVaa1MMGJmEMRQiYhfTtdarjJ/LAfS2czDEXqLsHgAJP5RSowBssxzaj0aY1UqplY29h9Z6ehOGRmyCQkTsYLyfUDTKNKa4hC4UIhJUlFJZAFxKqXzLYfpDYQ6FiAQNI2mxt9Z6vN/xMgD5gV9FwgEKEQkmWVrrBQGOlzfmxeeRfOjlNM5ZUIhIUDAyqOvK29lvXJNhWc4/B/9IqgUIajY3qRsu35NWRSmVYUQybwI4J0oxPCMzSlpprKgFY0xmMuVgpVS+MQ5iE4o9qwkhdsOIiBBiOxQiQojtUIgIIbZDISKE2A6X7w06d+6se/XqZfcwCAk5zpw5gyNHjuDo0aNfaa27BLqGQmTQq1cvFBUV2T0MQkIKj8eD7OxsAMCaNWu+qOs6Ts0IIa2CVYRyc3PrvZZCRAhpccrKymqJUFpaWr3Xc2pGCGlxkpKSMHLkSIwYMaJBEQIoRISQFsTj8SAqKgrJycmYOnVqo1/HqRkhpEUwPaFly5bhfEvHKESEkGZjNaZnzZoFpdR5vZ5CRAhpFv6rY43xhPyhEBFCmsWrr74KoOkiBNCsJoQ0k9mzZ+PEiRPo0aNHk9+DEREh5Lxxu9145plnUFFRgfj4+GaJEMCIiBBynrjdbuTk5AAAvF4vUlJSmv2ejIgIIY3GKkK5ubktIkIAhYgQ0kj8RaipxnQgODUjhDSKmpoaJCQkYP78+S0qQgCFiBDSAOXl5UhISEDPnj3x3HPPISKi5SdSnJoRQurE7XbjoYcewvr16wGgVUQIoBARQurA6gkNHTq0Ve9FISKEnINVhPLy8pCamtqq96MQEUJqUVFRgcWLFwMIjggBNKsJIX60a9cOP/rRj9CzZ8+giBBAISKEGHg8HpSWlmLIkCG49tprg3pvChEh5Gwrj+joaAwaNAgxMTFBvT89IkLCHGs/occeeyzoIgRQiAgJa6wiFCxjOhAUIkLCmA8++ACAvSIEAOp8m1yHKpmZmZo7vZJwQWsNpRS01jh+/DhcLler31MptV1rnRnoHCMiQsIMj8eDuXPn4tChQ1BKBUWEGoKrZoSEEVZPqLq62ubR+GBEREiY0BK7bbQWFCJCwoBDhw61WRECODUjJCxITEzEpZdeirvvvrvNiRBAISIkpDl06BBcLhdiY2OxcOFCu4dTJ5yaERKieDweLFy4EM8//7zdQ2kQChEhIYjVmJ40aZLNo2kYChEhIUZbXh2rCwoRISGE1hpLly4F4BwRAmhWExJSKKXw8MMPQ2vtGBECGBEREhJ4PB6sXbsWWmukpqY6SoQARkSEOB6rJ3TzzTcjKSnJ5hGdP44TIqVUlvFjbwCJABZorcv9znuNpxla66eDPERCgoZ/PyEnihDgMCFSSmVprVdZno8DsB0iSmdFSGu9znieoZRaqbWebsuACWlF2kpTs5bAMR6RUioDhuCYGILjMgQJAKabImSc3w9gVPBGSUjw8Hg8iIqKcrwIAQ6LiABkAVjgd8wLEaNEAIMDvKZcKTVKa72p1UdHSBCoqKhAu3btMHToUFxxxRVo166d3UNqNo6JiLTW+7XWgSbAGQCKjMfyAOe9xjlCHI/b7cbMmTPx8ccfA0BIiBDgICEKhOEJbdJa7wDggs+ktlIOMbUJcTTmNtBVVVWOn4r541ghMjyj6VrrG5vxHllKqSKlVNGxY8dacHSEtCzB3os+2DhWiADkAxjpdyxQ8906oyGt9SqtdabWOrNLly4tOjhCWorjx4+HtAgBzjOrAQBKqXz45Q9BfKJAouMCsCMoAyOkFUhMTMTo0aNx7bXXhqQIAQ4UIsMXWmkszZvHRmmtNyml9iulEv0EKpErZsSJuN1uAEBaWhomTpxo82haF0dNzZRSowAUmSKklEo0jpnkQ5b4zesHA6AIEcdhekLLli1DOOw96JiIyDCnNxo/+59OAsTzMQzoUZBpWgazqonTsBrT8+bNC/R9DzkcI0RGFNTgfxFrCQghTsMqQk7qJ9RcHDU1IyTUWbt2LYDwEiHAQRERIeHAgw8+CK/Xi+7du9s9lKDCiIgQm3G73cjLy8Pp06fRvn37sBMhgBERIbZi9YTKy8sRGxtr84jsgRERITbhb0yHYyRkQiEixAbCdXWsLjg1I8QGIiMj0blzZ8yZMyfsRQigEBESVMrKyuByudC9e3c888wzYZGs2Bg4NSMkSLjdbsyZMwevvfYagIAVAmELhYiQIGD1hG644QabR9P2oBAR0srQmG4YChEhrciZM2ewZMkSABSh+qBZTUgrEhMTg5kzZ6Jbt24UoXpgRERIK+B2u7F582YAQGZmJkWoARgREdLCmJ5QZGQkrrzySrRv397uIbV5GBER0oJYjenHH3+cItRIKESEtBBcHWs6FCJCWoitW7cCoAg1BXpEhDSTmpoaRERE4M4778TIkSORmMiNhc8XRkSENAO3241Zs2bhwIEDUEpRhJoIIyJCmojVE4qKauBPqbgYKCwEdu0CvvoK6NwZuOwyYMwYYMCAIIy2bUMhIuR8MATF/eWXyPF6gfbtkff00/XvwFpQAKxfD6SkAAcPivAUFwMdOwILFwJjxwLTpgXvM7RBODUjpLEUFAALF6K0Sxfk/PvfQGQk8k6cQOof/1j3a4qLRYQWLQJKS4F33gFeekkeDx2S4+vWyXVhDCMiQhqDKShr1sAVF4fLv/kGEyZMQGpsLDB5MjBkSOApVmEhcNddwJ//LI8ulxx3uYBJk+T4pElyXRhP0RgREdIYCgtxcPRonIqJQbt27fCTn/xEpmOmoBQWBn5daSnQt6/v0crFFwOHD/sewxhGRIQEwjSXS0uBlBR4PvsM2dHRuNjrxSOPPFL7mp07AY8nsPGckgLs2eN7vOoq37m9e4HkZN9jGMOIiBB/DC8I6enAvffC06kTsk+fBr75Bvfcc8+51wwcCKSlyfOCgtrvNWYMsHYtMGKEPHq9ctzrBV5/XY6//rpcF8YwIiLEGv1ERMjzDRsAlwsejwfZW7YACQnI27wZqQcOAMePn/WLAIhHlJ8P9Ohxrl80YICsiuXlAd27A7fdJsL12WfApZfK8XHjwtofAhgRkXDHL/rB4cPAyZPAhg3QWuPZZ58FAOQtXYrUm26SZfYlS0Rs3ntPhMcUkrr8omnTRKguvFAip7175fHCC+V4mC/dA4yISDhjWQk7u5rVtSswYwawfDnUkCGYN28eqj7/HKkFBTKdiokRrweQ6Ck/v3Y0c/HFwEcfnXuvAQPCPuqpDwoRCV/MpXVThAAxpvfuxZ/698eUt95CclqaiNVddwE1NYBSImBJScDixfLz44+fNbURFRX2xnNToBCR8GXXLlnt+tvfRETGjIFn2DBkL18OxMXhtpoadNq69VwvaMcO4L77gO3bgX//G7jnHpnWbd8OPPoo8MADge/ntxLH8g4fSmtt9xjaBJmZmbqoqMjuYZDWxhSD998Xw3jECODhh4E9e+BZuxbZsbEAgNz9+5H2zTfADTcAR48C27YBgwfLvy1bgCNHgN27gYkTgX/+E7jkEsmUvuUW4Pe/P3fKZpZ53HWX5BPt2SOraGFU3qGU2q61zgx0jmY1CR9MYzoqCjhzBpg5E/jgA+Djj+EZPhzZSUnAt98id/JkpGkt0dIvfgH8/e8SAV11FfCznwE33yyCNHAgEBcHfPMN8OabUjs2fPi5hvX69cBTT4n/VFICdOgAXH65vH7JEuD++1niYfcACAkKxcXAq68C11wjInHhhRLNdOsGPPEEjnk8iFEKj3btirSsLCAhAYiNBaZOlSr5tWuB6GjgscekTqysDGjfHjhwAJgzB9i3T+6xcCFwxRW+fKGCAhGhESMk8tmzB7j7bjn38MNyzHxdGEVH/jAiIuHBT38qy/L9+4sIDRgALFyIb9q1A2bPxhVFRXgxPh5pJ04AkZEiCg89JD7S6NHiE+3YIQL1/e+LoBw+LMenTAFuvVVEbc0aESrT1F6/Hrj2WhGYq66SSCg5WaKjyy+v/bowLn5lREScQ0Nmb13ni4uBTz4Bfv5z4MYbZXqUng7PuHF4ZNEi3Hv11bihXz9Er1ghAnT11UB5ubw+OVlM6kmTgN69geXLJTKKi5MplkmgMg1zVa6kxFfeUVgowmWeT0+X11lzkMLQwGZERJyBf+Jhenrtkor6zhcWApmZvsLS/v3hyc9Hdk4OapRC7/ffFzHwesWUHjJEhGPRIvGJ/vM/RUwiI8WYHjHClyU9eTLw4ovyLypKnn//+4DWvkJXs8zD6/Udu/himdZZyzvCuPiVERFp+wRKPLzqKpkyTZ4s0ctLL8kUqKREpjxTpvjOx8dLpPPaa0B5OTy//z2y+/QBqqqQpzVS//IXMZonTwY6dRLRmDBBMqe7dBExGztW/KKYGGDVKjGao6KAYcOA//5vIDUVqKqS1bKdO2UcgAjalCny+smTJYr67W9FcP76VxFLMwIK4+JXChFpe/hPsY4cOTfxEJDn3bsDubnA+PE+M9hq/E6aBLzyClBdDXznOzj55JPIvuYaIDoaeX36IPXll+WP/7PPZBWtsFDEYsYMOb5unUQ4L74oHlPXrpIndOSIvO+kScBbb/nExOsFFiwQQQJkLKNHy1iGDBHBXL1apndLl8o4zde9/rrvdWEGhYi0Laz5NrfcIsKyerVkMvtTXCxL6126+MzgDh2Azz+XZfHt22WVrEsXiXIGDkTHBx7AnR4PMisrkVpVBVx0EbBihUQxr7wCZGUBlZU+X+juu0VoSkrkfVas8InOlVeKUPXvD5w+LRHN66/XLmI1I6FJk2TqNXQo8Je/iJl9+rREaW+95ctTClOY0GjAhMY2gLmMbZ2CAeLV/Pa3wP/+b20j9/HHJYlw717xhSorZYm+Sxfgyy8l67m6GrjkEniGD0fFihXoPWKEHNu3T66ZMkVKNbZsETF75RURNDMqO3xYIqPevYHNm0WI/MdsvS5QtnSgawBZyfvkE/Gvxo6VcYVwkmN9CY2MiEjbIVDtFwD88IfA734n0cNTT/mOm6tRc+cCzz4r05sLLpBz8fGSEX3iBDzR0ciOjcUFffrg+cJCRPbvL8IycqSIS0GBbyUsP1+mX/6rcmvWiIgEWplbvLj+zxWo4LW4WITzk09qf17T16qr9WyIwlUz0nYI1E4VkD/I226Tac+aNRK9rFkjZu+ll0oEER0tUc6pU3J9UhJQUQFPr17ITksDTp3ColOnEDlwoGRBr1oFZGfL+/zqVxJdmZHSnXfWXnUz/Rug/pW786Eu0W2o9WyIwoiI2I8ZZRQVSRSzZMm50UDfviJGJSXSZiM5WaIjUyiOHhXxSUqS6c2xY/DMno3sPXsArxd5O3ciNTdXXj90qPQD6t0b6NVLluTLymQF7PhxEaWsLGDWLMl+fuUVyT8yBTDQyt35RjClpeKBBaKuViIhDIWI2IvVnJ4zR0oo7rhDsp+HD5epT48evhUl/z/2EydECEpKgHbtxDj+7DNg2DCs/+tfgXbtkLd/P1IPHwY+/FD2FbvnHsnh2bcPePtt4Ntv5X4zZ8pUb9UqGVP37kCfPr4l9YYimPMRokA9rE3CcBmfUzNiH9b8oClTRBDatZOVL49HooYpU4DvfrfudqrTpknCYWWlTK0SE4GVK4FrrsF9WuPpL75A6vHjMnVbu1YKVrOzRWyefRYYNEgSFS+9VERhyhQRJ7Modv782smJgWhKIqI1ydGKOQ0Msx7WDQqRUqpX6w+DhCVWn8QsSv3e90QUoqMlgrn1VqnvGjKk7vfZvRu4/nqgY0e4d+zAY0uX4uuCAkSPGYPk6moRuIwMiXp+8xtfPVdhofzBx8RIUqSJNcoxoxMzgglEUyIYs5f15Mm1fS9r69kwojERUXhmWJHWxxplWItS582TzOavv5Y8n6ys+s3b0lJg2jS4u3ZFTnQ09n31FU5UVUm0dfQoUFEhUc011wDHjokHtWWLtABZswa4/XbxiazRiX8JRmtEMGYv65ISaTdSUhK2PawbzCNSStUAuEJr/WlwhmQPzCMKMsXFMqX66itJ7vvyS/FkunaVKVKPHiJMn30mSYX15fC88w7caWnIAYAzZ5C7bx/Svv1Wsp8rK6VdR5cuEmXdcIMkL2ZmAv/4BzBqlBjX69fLdMxMPHzlFdmFdeFCnzAUFEgCo3mNNYGxIfFgd8Z684gaI0QbARQBeENrvbMVxteiKKWyAJj/28rQWj/dmNdRiIJIQYH4MxUVMiU6dUqmT3FxYlCbyYbz5skS+YYN4uVY83UsJrfnn/9E9tatgFLILS9H2u7dEk1FRYn49O8vkVBZmdSU3XEH8MILIg5z50ov6uXL5X0LCyUS+vBDWZUzSzBMGpPAGOjzhnl3RqAJCY1KqY5a65MAoLW+0Tg2VimVobV+q/WG2jxMEdJarzOeZyilVmqtp9s8NGJSXCzRRnKyCMu0abLydfPNEoG8/754Qp07Axs3ijD96U9iKD/+uG8KZJrcGzag3fbtSI6Lw4O7dyOtpEQiG0DeIy5OROyLL0TsDh+W0o8tWySiKS2ViGznToly0tOBjz+WbG5/EQLOfzeOhgp2wyxxsS7OESKl1EgAGsCfrce11uuVUlcopeZprZcFa4DnyXSt9ZXmE631fqXUKDsHFPI0ZsphvWbPHlkZmzJFhGfAAIlctm6V6OX0aTmvlNSRbd0qHRLnzfMVtEZHA3fdhWO7d8O1fj26rFmDfI8H6oEHRGy+/FKmY4mJsvz+7rsSiXz9tTQhc7vFk0lOlmRGQMb37rsSOXXuLKt2xcXNF4nGJC5SiGqb1UqpHwEo01r/OdDFWutPALyslHpKKdUxGANsLEqpRACBqgbLKUatREM9goqLZRp0223Ap58C113niz5+/nPgjTekxMHrleP33SdtVi+4QOrEvF6JaF591be0bnRKdJ86hbnLl+PVfv0Alwtqwwbgxz+W9q4dOgD9+omoxcSIB2QKzOnT4j+tWCERmRnhpKaKUCYlSTuPTz+VXtJNyZq20tLL/iGKf0RUAiADQJ1ekNb6BICFSql5SqmN/ia2UqqX1vpAi4+0YTIAlAc47jXOkZakoSmHxyNTq6+/FtE5fFhqxbxe2ZiwpESu/8tf5JoxY2T3i5gYMZC3bRP/KCGhdsTgcsGdmYmc3/0OiIjATSNHynFzayCPR55//rn8oXftKj2GTpyQaOiFF2Qp3/+zPPusnL/jDp+P88tfyvHmTJ+YuNgoakVEWus/AfhEKXVnfS9SSo0AcBWAHUqpPxiiNMg4bZcf44LPpLZSDiAx0AuUUllKqSKlVNGxY8dadXAhR31TjmHDpHUHICLw0UfSrGzRIvFqUlLENK6uFrGJiJAmZAcOyDL6p5+KkERGyvL7hAniDxUXw+12IycqCqisRG5GBtLKyiRq2bJFxKuqSqZnFRXyft26STT0/vuSxNihg4zJyksvyRjeeEOiLjP6euMNidReeqnpvycmLjaKczwirXUJJDICACilntRaLzLEZzyACQCSIH/gL0MEYAaAp5VS5hLcotYeeEugtV4FYBUgq2Y2D8dZ1Fcr9Yc/yNQqJkaio6gon7czYYJER3/7m0QESsl2PCaJiWImR0TI9bGxct2+fahavRqPDxoEREcjNyYGadu2iYBVVvqKXkeOlJYhSokgvf02sH+/vF9lpURs/ivFW7dKJBRIVMeMkfdrKtbExUDL/vSHADQuoXG6UqoMwCYAEwG8CeBGrbVLaz1Daz1Ba90HIk4TYRExG3AFOBYwGiLNxD/T2Nx6ecIEiUyuuUZydqKiank76N5dzOikJCk6BUREevaU1h0nT4qHo7WsfqWlyfSqqgpROTl48OBB5H78MdKuvFJ8njNnZDp2ySUSoW3bJgKYkCDis2uXPO/YEfjBD+R+wZ4OMXGxQRpT9JoIYB2AlcbULSCGd7ROKRUwTyAIFCGw6LgA7AjyWEIP/9Wx/v1lSjR6tOT5mA3Jdu2SKddnn0nS4Nq1co3LJV7LO+/Iilb79jKFmzJFrqmululUbKxMrxITRZSKi+EZPx57+/bFiD/8AQOrqmQ5/uOPJVrp3Fkis5QUibQ6d5YoY/duMb1PnQIefFCuLS0VQ3vJktqf7aqr5DM89FDtqMjrlc88bFjzf3/nu+wfZjRGiFZprWecx3uubOpgmoPWulwptV8plai1tprWiVrrTXaMKWQI1L61oEAilR/8QLyYlBQRohMnJBoqLhaBuesu37TkootEiJKSRBDKysSYvuAC4JFHRIj69pWIITkZiI+HZ+xYZH/+OSJOnsTQY8cQO3++ZFlffbWv7KNbN1+y4+OPy+pdZSXw5JMSkd18s6zULVsmq3L+gjBjhiQwTpwI/Md/+KZPv/ylJDvOOJ+vP2kKjcmsTjd8ozaPkdCYaGZTK6UGQ3KLGjTQmVldB3W1b/V6RWBqauSP/Hvfk94+R45Iv5/yciAnR6KNe++VvswffCDTsJgYKWatqZHp05df+pINAfGMMjPhcbuRnZ4OVFUht7oaaZdfLsKyaZNETB6PrID94x+Sk2QWzy5cCDz/vNwjMVEis4QEaaRm7TltpaBAEi27dvUdO3pUxs4pVItQX2Z1gx6RU0QIOGs+lyulRimlxgEYxazqZtJQQl5xsUQMBQUSlcyYIZHQ3XfLv7//XSKYiy+W6vc775QcnzffBP71LxGQ7GyJhGbPltqzhx+GZ8sWZHfrJiI0dizStm8XgXr2WYmiBg4U/+jgQXnNxIkilqdPS3T2ne/Iz7NmSfT06aciKnVNj6ZNE5EaNEgirEGD5DlFKCiEXGM0Q4xIS9FQJ8HKytrHrKtE7dpJJHL11cDTT8tK1tChsj/YqlXSkqNfP+n5U1kpkU6/fsCSJdg5cKCI0N//LiIUFycpAY88Iu89ebIYvj16iCleViY+j8mFF8q/zZtlmheoqZo/9HFsI+SEiLQwDSXk9eolAjBqlEQ3ppk9caLs79W5sxjLSUly3cGDEmVFR8t+Yab4HD8OlJejOjkZkVdcgdsqK3FDVRUu0Fpe06OHCFxS0rk9e+bPl3t06+YrRg3Ubpa0WShEpH7GjPGVcbz+uuTcADI12rtXxGbWLMnDGTdOBGvrVvFo2reX1xQWSiS0YcO5pvf27ZJkOHEi3C+8gLz4eMy57DJcHBeHCzZskHM5OWJyHz3qW/q2iozZ47mh3TRIm4VCROpnwAAxcG+9VXyhG28Uo/rNN2XatXu3HI+NlQzm996T1/XpIx4NIFFSnz7Ayy/Litrf/iaCMmaMTNmmToW7vBw5l1wC1NSgQ1SUiM6tt4rBbGYfp6fLcr8/LJVwPBQiUj/FxWIqX3aZeDEnT0q5xk9+ItOf1auB6dNlmmXt07NkiaymFRbKVO2ZZyRXqH9/Xy3XwoVASop0VuzYEYiIQO6TTyItLc13/zVrfC1dze2b/Vfvwnir5lCBQkTqp7BQIqI77jg3Gpk6VerH3npLWmukpEjDeXPadPq0TJlGjpScntdek1UzQKZw11yDw9/9LnIGDgQA5Obm1hYhwDftYqlESEMhIvVTWiqPgVpZ7Nkj+UIdO8rSuBnlmJ0HzSnT7t3SEP+xx6QGzSIinW66CVf961/4QadO54oQUHvaNW2aZGcXFvr2NmvMahhp81CISP2kpIhfY66cmaUeu3ZJC49+/WS61qdP7TYgkZHAf/2XZDKXlADjx0vx6O7dwEcfwZOYiAsWL0bCyZN4YOpUybr2en1JiS+9JF7S/v3iFZlNyrjEHpI0mFkdLjCzug6Ki6VhWYcO4gP9/vfAtdcCv/61tOyoqpJC1upqiVbuvFNM6W3bZAqVlSVRy759IkpHjsA9ejRy/vpX9IyLwxPbtokJffvt0pi+e3fg//5P8opiYiSS+uILycKePZsJhg7mvHtWkzClrrav994rHk9OjojQq68Chw6JUCQmSjLhTTdJxfy8eVIR/53vSEKh2Z/60CHA7YZ70iTkvPceUFODLLMv0RNPyH0SEsQEr6nxZWYPGCCR0oQJ8l7s8RySUIjCDavYRBgVPjU18tzrlWJUM8fH6vds3y7/PvpI6rvS00VwPvpIGogtXy45Qr/6ley0um+fCJj5Xh9+CPexY8h5+20gLg65w4Yhbe1aaVK2YoUI35EjUhfmb4y7XGKMb9jAHs8hCrecDiesPaY7dZLyh82bZVp14ICUZFRW1u4PvW6diFdNjSQx/s//SF+h66+X6RoADB4sdV+FhVJbVlEhmdaWbof6/fexIiNDyjbi45G2caOUfIwd6+t3/e67Mo5AxvjFF8sjezyHJBSicMHaY/ryy6Vp/Z//7NvCZ8IEiWxM4QFq7zSRkiJRjikSNTXSFH/yZJkyXXCBVNd/+KFc065drdurTp0wr2dP5FVUIO2hh0TwNmyo3Zr1hz+Uco5NAbq27N0rj0xcDEkoROGCtYre+rPLJV7OgQO1hcfE3GlizBjJFdq+XY6npIjgZGfLvmOlpeIRZWZKM7KOssmLx+PBqlWrUH3sGDrt3o3Uzp3rruj/4Q/FJ1q9unaPZ69XOhsePcoezyEKPaJwwVpF719R37u3RDuAL4HQxMzjGTBAop9HH/W9Zve+kHcAAA0HSURBVOlSKWp98EFZTcvPlzYdH30E/OIX8Ozdi+yoKEBrjNmwAV0uvVREr66K/gEDZJn/xRdlY0VTdAoLReRmz6Y/FKJQiNoirbFPurWK3r+ivlcvWTL3emsnEPqXT5hFpcuXiyeUkCBR0urVsrT/3HMytVMKnn79kP3tt4DWyIuJQZfsbCAvTwpkCwvrrujv21eiJcDXtH7YMOlzRBEKWZhHZNBm8ohaa590a6fFgwd9PwMS6QwdKv7O0aPSVqOmxlc+4X9f6/7vSsmxr74SMcvIgKdvX2R/9RUQFYW8o0eR+uWXvk0Yp01ruOsjs6VDEuYROYXW3Cfdv1briitk+gNINJOeLhsipqc33EwsUHbz/fdLpPXeezhRWooOy5bhp8nJSHW5JDu6Rw+foLFujPjBiMigTUREZuP3QK0u1qyRUonm9twJFM1oLcLTnCngkCH49y23IO6JJwAA1dXViIyMlHNLlsg0y+xlFGgszb0/afMwInIKdZm4xcUybdq8WZ435w+2lWq1PJGRyNm9G3f/8Y+46aabfCJkw1iI8+DyfVvCf9NCwJeEGBEh5RVRUTKtGT787DbMduPxeJCdnAxdUYFLu3evfdLcG2zIEHsGRxwBhagt4b9PuukZPf+8bLmTmipR0Zw5kr1sbuVcUGDbkD0eD7Kzs4HYWOQdPozU2bNlGrllizxOmCBTP+4NRuqBU7O2hL+J++GHskz+4IOyqmX+cbtcsj1zSYk8b66R3UROnTqFnJwcAEDe008j9Y9/lCxr624aFRXM/yENQiFqa1ibf23eLNOx/Pxzs5HNxMODB6Xma9o0Wf1qLcM3QG5T/IABmDRpEi677DKkpqbWHjsNaHIeUIjaIlYTNz1dfl6xoraRvXevtNZYuNDnv6Sn166Ybyn8tpz2bNuGrxcuxCVjx+IW//vQgCZNgB5RW8bqGVmNbK9XNij0emVq1ru37EzqXzHfElhzm6ZMgad7d2Rv2YLnundHdUveh4Q1FKK2jNUzioqS9qkvvijPO3WS3kGAJAKadVmBClebg2VKeNaYBpCzeDEiW/I+JKyhELV1pk0Tj6iqSjoiLl8utVcREVKo6r/rKeCrmG8JSkuBvn1ridDZ3TZa8j4krKFH5ARM32XxYp9p7PFIZnSgMoyW3HDQmBK+u2MHAL8tf7ixIWkhGBE5DVOQVq8Gvv1WarismBXzLdW3x/Cpfjx2LJYtW+YToZa+DwlrGBE5lSAUjno8Hry4di3mjR6NpKlT0ZUFqqSVYNGrQZsoem0KrVQ4avWE8vLykHr8OPODSLNg0Wso0wp5O+eIUGqqlJdQeEgrQY+I1CKgCBHSyjAiIrWIi4tDz549MX36dIoQCRoUIgIAOHLkCDp16oSkpCQ8YTQ3IyRYcGpG4Ha78fDDD2PVqlV2D4WEKRSiMMftdp9t5XH77bfbPBoSrlCIwhirCNGYJnZCIQpTqqurkZeXB4AiROyHZnWYEhkZiVmzZiE+Pp4iRGyHEVGY4Xa78VtjB9V+/fpRhEibgBFRGGH1hK6//nrEx8fbPCJCBEZEYYJVhHJzcylCpE1BIQoD/EXobCsPQtoIFKIwYM+ePVBKUYRIm4UeUQhTVVWFqKgojBo1CldffTXi4uLsHhIhAWFEFKK43W7cd9992LVrFwBQhEibhhFRCGL1hBISEmweDSENw4goxKAxTZwIhSiEOHbsGEWIOBJHTc2UUlnGj70BJAJYoLUu9zvvNZ5maK2fDvIQbcXlcuG6667DLbfcQhEijsIxQqSUytJar7I8HwdgO0SUzoqQ1nqd8TxDKbVSaz3dlgEHEY/Hgw4dOqBTp074sbn7KyEOwhFTM6VUBgzBMTEEx2UIEgBMN0XIOL8fwKjgjdIezB7TS5cuBXdkIU7FEUJkkBXgmBciRokABgc4X66UClkxsja6v//++6GUsnlEhDQNRwiR1nq/1jopwKkMAEXGY3mA817jXECUUllKqSKlVNGxY8daZrBBIuBe9IQ4FEcIUSAMT2iT1noHABd8JrWVcoipHRCt9SqtdabWOrNLly6tNNLW4eWXXwZAESKhgWPMaiuGZzRda32l3WOxi7lz5+LkyZPsJ0RCgqALkeHZLGjk5dMN09mffAAj/Y65AlxXZzTkRNxuNzZs2IAZM2agY8eO6Nixo91DIqRFCLoQaa03AdjU1NcrpfLhlz8E8YkCiY4LwI6m3qstYc2YHj9+PJKTk20eESEth6M8IsMXWmmNkpRSowxR2m+snllJNITP0fiXbVCESKjhGCEypnRFpggppRL9lubzYVniV0oNRjMir7YCa8dIOOAIs9owpzcaP/ufTgJkBcxYjh8FmaZlhEJWdUVFBeLj47Fo0SKKEAlZFLNxhczMTF1UVGT3MM5y8uTJs2Z0dXU1IiMjbR4RIc1DKbVda50Z6JxjpmbhhNvtxgMPPIDf/OY3AEARIiEPhaiNYfWEBg8OVLVCSOhBIWpDcC96Eq5QiNoIp0+fxk9/+lMAFCESfjhi1SwciI2NxdSpU3HRRRdRhEjYQSGyGY/Hg7KyMgwaNAjDhw+3eziE2AKnZjZitvL42c9+hsrKSruHQ4htUIhswtpP6NFHH0V0dLTNIyLEPihENmAVIRrThFCIbGHjxo0AKEKEmLDEwyAYJR5aayilUFNTA6/Xi86dO7fq/QhpS7DEow3g8XiwYMECHD16FBERERQhQixw+T4IWD2hiooKm0dDSNuDEVErw902CGkYClErcvDgQYoQIY2AU7NWJCEhAX369MG9995LESKkHihErcDhw4fhcrkQHx+PxYsX2z0cQto8nJq1MB6PB/Pnz8cLL7xg91AIcQwUohbEakyPGzfO5tEQ4hwoRC0EV8cIaToUohagpqYGTz31FACKECFNgWZ1CxAREYG5c+ciKiqKIkRIE2BE1Aw8Hg8KCwsBAOnp6RQhQpoII6ImYvWERowYgYSEBJtHRIhzYUTUBPz7CVGECGkeFKLzhE3NCGl5KETnSUlJCSIiIihChLQg9IgaSWVlJaKjo3HddddhyJAhaN++vd1DIiRkYETUCNxuN2bOnInt27cDAEWIkBaGEVEDWLeB7tatm82jISQ0YURUD9yLnpDgQCGqg7KyMooQIUGCU7M6SEpKwo033ojhw4dThAhpZShEfrjdbkRFRSElJQVTpkyxeziEhAWcmlkwPaFly5aB+70REjwoRAZnzpw56wnNnj0bSimbR0RI+EAhMjh48CAA9hMixA645bRBp06d9M6dOylChLQS9W05TSEyUEodA/CF3eNoAp0BfGX3IGyAn9t59NRadwl0gkLkcJRSRXX9XyaU4ecOLegREUJsh0JECLEdCpHzWWX3AGyCnzuEoEdECLEdRkSEENuhEBFCbIdFr8RRKKWyAHiNpxla66ftHE9rYXxOAOgNIBHAAq11ud/5kPk9UIgcTjh9Yc3PorVeZzzPUEqt1FpPt3loLYpSKktrvcryfByA7ZD/xiH5e6BZ7WDq+MLma63r/MJChMqRX1ijROBKv2P7zM8bChj/jaZrrRf4HT8O4Mda63Wh+HugR+RQjC9srS+eITguQ5AA+UKvs5zfD2BU8EbZciilEgEMDnCqXCnlyM9UD1kBjnkh/21D8vdAIXI24fSFzQBQHuC41zgXEmit92utkwKcygBQhBD9PVCIHEoYfmFd8HldVsoh3ljIYkyxN2mtdyBEfw8UohAiHL6w4YbFM7rR7rG0Jlw1CxEsX9grG7zYubgCHAt1Uc0HMNLvWMj9HihEbQDDs1nQ4IXCdMN09ifUv7BFCDx2F4AdQR5LUFBK5cMvHQMh+nugELUBtNabAGxq6uvD4QurtS5XSu1XSiX6fc5E4/cXUhjT7JXW/+kopUZprTeF4u+BHpHDqecLWw5gv7F6ZsXJX9h8WFYKlVKD0QwBb6sYEXKR+d9UKZXot9IZcr8HJjQ6GOPL6TXMaTPXJhOQKMsQqUQzm9r4wk53akIjcFZ490OiPUdnigfC8Pr21XE6yYyCQu33QCFyKOH6hSWhCYWIEGI79IgIIbZDISKE2A6FiBBiOxQiQojtUIhIm8ds/KWU0kqp40YCZ6DrNhrXvGmsKhKHwFUz4hiUUm8CGFVH1wEzw/wNM6+KOAdGRMRJrATgn2UM4Gx3SoqQQ6EQEcdglKaUAxhvPW5kjIMi5FwoRMRp/BrABPOJUdYyytoSlzgPChFxGm+i9vRsEctWnA/NauI4jB0tNgHYBmCVXzsM4kAYEREn8msA4yBtcSlCIQCFiDgRF3y9uUkIQCEiTmQwgI12D4K0HBQi4iiMVbIMOLwjIakNhYg4jVEAc4ZCDQoRcQRGvVk+pF8zlFL5rCcLHbh8TwixHUZEhBDboRARQmyHQkQIsR0KESHEdihEhBDboRARQmyHQkQIsR0KESHEdihEhBDboRARQmzn/wGur6XcIbHB+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Init Loss = 76.43615\n",
      "Final Loss = 3.7835722\n",
      "Fair dummies test (regression score), p-value: 0.000999000999000999\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAElCAYAAAAY6L5SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYn0lEQVR4nO3dPXAj533H8d9fJxUejSc4KCluVFjCKT3Fo1tnYuHcO0PqyjQWqEpqEjKXGanMiRenjohLipQnUq4zJuRMnDIkfak0GYdw3Pt4cJwy8j/F8yy5xC2AZxdvC/L7mcHcYRe7+yxf9sfnZZ81dxcAACleWXYBAACrg9AAACQjNAAAyQgNAEAyQgMAkIzQAAAkIzSwUGbWMTMf8TpbYhlemNmRmbUWUYY6MrPN+HWY9DpZdlmxPK8uuwC4se5L6ktq5JYNUjaMF61zd78/ZRluu/sg7nNd0kNJJ2a25e69Kfc9FzM89yI9Se8NLTuR1JW0P4fjYQURGliW4+yCXcFMLmD547v7qaQtM9uXdCDp9iyOMQdzu3jHr8dpfpmZSdJZ/PpMxcwaktrufjjtvrA8NE9h5bh71927c9r3tiSZ2d489j+teZ77AmwoBDJWGKEBvOwLSZvLLgRQR4QGlqVpZo38K1thZm0zO8l1UF/5q9/MDszsIPd+38z24nZHZrYzZdlOJF10iA8fLy5zM9vMvd+Prx0zO4vl3onndpR19Me+k2m2KSrLXq6Tei+3fCcuy/bTmfLrckWZ45rZkaSj3NfOzaw9y/JgMQgNLMuZpBf5V27kUkvSrqS7krYk7Uy44DUVagYHChemaZtv+pJUciRVU1JWxvuxDHuSvor/3lPo6D+YcpsrYoBsKnRgb0najBfpdtzHewpfx21J5yXOZ6yyx40d91uS5O4WX7UcbIDx6AjHstwe1RE+1GbfN7NDhQvOuDBoSbo3iw7buC+5e7/kdgN3fyxJZvZI0o5Ch38vLtvXyx3ZVbZRXNdSuHBfnLeZ7SqMAnukEDj9+HUuey4jLeu4qAdCA7UUm34eKFzA1xWGg45zOqPAkMJfyFVGdh1n/3H3QRx5lL+noegv/SrbZLLmnZO4XV4vbvvCzHqSjrJwmoFlHRc1QPMUaifei/BQ0lOFGkbKEM1Z/kXbUegML6soaCY1CVXZ5sr2ueae/Gvg7lnz3kDS3nBfyJSWdVwsGTUN1Ers9F2XdDdrHir4a3aex99TuOFwd2hVvqO+oXroSWqY2fqoWla8J+Iw1txmdfFe1nFRA4QG6iarMezGC/i6Qvv5XDtNY1htK9Qy7hf0t7RiW/65pCfzLEsqd++bWVfSgZltK3zt1hWa9Z5K+m78d6DLO/CXedzsj4DNuE50hq8emqdQK/Fi/Vjh4n2icNE5VLU+hrFyQz9dYcRSU6GGM3whexTXnSn81fxI4QI48zKVFW9GzDrLz3TZrHeq0B/0VVy+oTh6aVnHjbWSU4WvYS1vnsRktuhnhMeq/fsKf80l/RDH4ZZZO2+LjjUAWI6FNk/FJoCsip80Bj4LjGy+GjNrmdl+Nt0DAGBxFl7TkC7C44m730v47Mnw58zsLI7QAAAsUK37NGJT1nrBqgFTEADA4tU6NBSasEaNY7+xD8sBgGWp+5DbpopvdBro6sN7LsQ+kI4kfetb37r31ltvVT/6734n/f73kz/3yivSt79d/TgAUCNff/31b9z9j4rW1T00SovzFnUlaWNjw4+PjydsMWZff/I9ffPhn0/83K3P/0n2rz+vfBwAqBMz+/WodXVvnpJCbWNYXe7IBYAbpe6hcazigGhq6LGUAID5q3VoZFMrF8z102D6AQBYvGWFRlGTU3bj3sFQSOzp8kE12T0eBAYALMGi7wjPHt5yX9J6nJDuLPfQnZbCXP1NXU5o1s09EayhMI0Id4MDwBIsNDTiVNeP46tofU/S7YLl0z6+EwAwA7Xu0wAA1AuhAQBIRmgAAJIRGgCAZIQGACAZoQEASEZoAACSERoAgGSEBgAgGaEBAEhGaAAAkhEaAIBkhAYAIBmhAQBIRmgAAJIRGgCAZIQGACAZoQEASEZoAACSERoAgGSEBgAgGaEBAEhGaAAAkhEaAIBkhAYAIBmhAQBIRmgAAJIRGgCAZIQGACAZoQEASPbqMg5qZh1J5/Fty90fJ26TaUjquvtgHuUDABRbeGhkgeHuh/F9y8z23X17zDY7GgoJM9uXNHIbAMDsLaN5ajsLDEly976k9oRtvltQqxiYWWPmpQMAjLTQ0IgX+fWCVQMzGxccrYL1DZqnAGCxFl3TaEkqutCfx3Wj7Eo6MrM96aKJa3/2xQMAjLPo0GjqsgM8b6DQuV3I3XuS7knqmJlL6rv76XyKCAAYZSmjp8oys5akDUlvS3qoUOvYdvduwWc7kjqSdOfOHT179qz6gdtt+WuvTy5fuy1NcxwAWBHLCI1mwbJJHdq7udFVu2b2VNJXZtaPtZALMUi6krSxseFra2uVC+off6Rv3nlz4udu9XqyTz6tfBwAWBWLbp46VnFANCUVNjfFDvCj/LLYNLUl6f6sCwgAGG2hoRFHO/ULhso2hmsMCY4lPZ9NyQAAKZZxn8aeYp+DJJnZuqRe7n3LzA6yYIlh8qBgPx3FZigAwGIsvE/D3btm1onNTg2FaUTyd3a3FG72a+pyeO4HcbjtWXzfkHTIfRoAsFhLGT1VNOopt64n6fbQsoHCvRoAgCVillsAQDJCAwCQjNAAACQjNAAAyQgNAEAyQgMAkIzQAAAkIzQAAMkIDQBAMkIDAJCM0AAAJCM0AADJCA0AQDJCAwCQjNAAACQjNAAAyQgNAEAyQgMAkIzQAAAkIzQAAMkIDQBAMkIDAJCM0AAAJCM0AADJSoWGmX1/XgUBANRf2ZpGz8yem9nfm9mfzqVEAIDaKhsatyU9lHRX0lcxQB6Z2drsiwYAqJtSoeHuv3X3rrv/QFJT0l8pBMiBmf0yBsh35lFQAMDyVe4Id/eBpOeSziVZfN2T9Csz+2cze2sWBQQA1Efp0DCztdin8VzSgUKT1ba7v5Orgfxa0k9nW1QAwLK9WubDZvZfkt6W9AuFpqkv3P23+c+4+8DM9iR9MLNSAgBqoVRoSOpJ+tzdn0343HNJ90etNLOOQrOWJLXc/XHKwc1sR9Ig29bdD1O2AwDMRtnQOJH0QNKV0DCzn0p66u7/KIUOc0lfFe0gC4zsgm9mLTPbd/ftcQc2syNJW7EvRWb2wsx62XsAwPyV7dPYkfTvBcu7kj5M3Md2vobg7n1J7XEbxBrGwVBA3CMwAGCxyobGXUmnBctPJK1P2tjMGiM+NzCzccHxUKFp7EIMGwDAApUNjb6kdwuWr0v6VcL2LYU+iWHncd1LYtA04v83zaxtZjtxOQBggcr2aXQl/YOZDdz9X6SL+ai6kh4lbN/UZQd43kAxGApsZOtz/SDHCsN9X+psj30mHUm6c+eOnj2b1Gc/Rrstf+31iR+zdlua5jgAsCJKhYa7PzazNxSmEPHcqifu/uPZFu2KhkItJyvHwMyaZrbu7leay9y9qxBi2tjY8LW16jOc+Mcf6Zt33pz4uVu9nuyTTysfBwBWRdmahtx918z+RuGvfJd06u4pTVOZZsGycU1N/Xjc4Watc4UO9KI+FgDAHJQODeliSG2VeySOVRwQTY24+Lt738xG7Y/RUwCwQKVDI84p1VLBxd/dfzJu29is1DezxlDNoeHuvZEbSqdm1hoaMdVSCCEAwIKUnUbkLyV9pjA54TCXdCthN3sKHdWP4z7XlRtOa2at+JkPcsGyG1/buW36w/0ZAID5KlvTeCjpbyU9Gp5zKpW7d82sE+/LaChMI5K/G7yl0FfRVGx+cveemTXiTX6S9Ia7j5ymBAAwH2VD41xh7qlKgZGJI5xGrespzJw7vJx5pgBgycre3LcvaXMeBQEA1F/ZmsZdSR0ze6jQCX1l9JK7P5hVwQAA9VNlyG2+mWjkWFgAwPVT9o7w1JlsAQDXUJXHvf7IzI7N7P9yy56a2Q9nWzQAQN2UCg0z+0zhfolHQ9seSPrrGZYLAFBDZfs0OpI23f1nQxMWJj1PAwCw2ko3Tyk8/1u62gneUtrzNAAAK6xsaDyRdGhmawrThsjM3pX0eXwBAK6xsqOnduPcUKeSZGbPFaYC6c75eRoAgBqo8jyNrRgc78VFvZLP0wAArKiqz9PoK/ckPQDAzVB2avQ/G7N64O4/m7I8AIAaK1vTGDXTrCvUPP54uuIAAOqs1Ogpd39l+KXw3IsvJW3NpYQAgNqocp/GFfHpeh8oPG0PAHCNTR0aUVPhaXsAgGusbEf404LFDYU7wgcF6wAA10jZmoYVvH4r6ReipgEA117ZO8Lfn1dBAAD1N6s+DQDADVC2T+P7ZT7PzX4AcL2Uvbmvp3Aj3/CzwUctu1WxXACAGirbPLUl6StJd4du8PtS4eFM+Rv/CAwAuGbKhsZnknYKZrXdlfRwNkUCANRV2dB4Q9K9guUt8bhXALj2yobGgaSumf0wWxA7x79Q6O8AAFxjZe/T2DazpqQvzczjYlN4kh8TFgLANVf1yX1v67I56pQn9wHAzVA6NMzsR5I+lLTm7q/GZU8lPXX3n8y4fACAGinVp2FmnymMlHqkq/dlHKjE6Ckz65jZZnztlClD3P6g7DYAgOmVrWl0FO7H+JnZlXv5TpQ4esrMOpLO3f0wvm+Z2b67byduvy5ps1yxAQCzULp5StLz+G8+NVqSUvs1tt39Ytiuu/fNrMwMua0SnwUAzFDZIbdPJB2a2ZrCNCEys3clfR5fY5lZQ8U1kkFKcJjZZlZDAQAsXtkht7tm1lIYYisze67wEKauu/84YRejHtZ0rgk1iHjcfpnyAgBmq+qQ25ak9+KiXokht02FgBg2UAifcdapZQDAcpWdGv0vFGoVfS3wr/7YdJV0x3nsaO9I0p07d/Ts2bPqB2635a+9PvmY7bY0zXEAYEWUrWn8QKEv4++mOGazYNnIWkbsB5G7Jz2D3N27krqStLGx4Wtra1XKGPb18Uf65p03J37uVq8n++TTyscBgFVRNjT2JH1hZj13/48KxztWcUA0FftJCnSki6G2F+L9HYMYEgCABSgbGvclvZB0amanGmqicvcH4zZ294GZ9c2sMVRzaLh7YfOTuz8eXmZme0XLAQDzVXbIbUOhRvClwn0ZNvRKsadYe5AuahC93PuWmR1kzVIAgPooO+T2w2kP6O7dOI1IWyGEWkN3g7cktRWarK70Y8RttuL/9yUdjKqhAABmb2RomNkvJd139/8eWv7W8LKyxvVDxBC4PWZdT1LSlCMAgNka1zx1V0Od1mb2B5LO5loiAEBtle3TkNL7LgAA10yV0PDJHwEAXEeTQiPphjoAwM0wafTU22aWnyuqKUlm9m0VNFO5+//MsGwAgJqZFBpFw1lNo2sgt6YrDgCgzsaFBsNaAQBXjAwNd3+yyIIAAOqvyugpAMANRWgAAJIRGgCAZIQGACAZoQEASEZoAACSERoAgGSEBgAgGaEBAEhGaAAAkhEaAIBkhAYAIBmhAQBIRmgAAJIRGgCAZIQGACAZoQEASEZoAACSERoAgGSEBgAgGaEBAEhGaAAAkhEaAIBkry7joGbWkXQe37bc/XHiNpJ0V1JD0q67D+ZURABAgYWHRhYY7n4Y37fMbN/dt8dt4+7d3PtNSScKAQIAWJBlNE9tZ4EhSe7el9Qe9WEza2koHOL2zRgeAIAFWWhomFlD0nrBqoGZjQwOSZ2CZeeSmjMpGAAgyaKbp1qSivohzuO6l8SayO0R+zqeXdEAAJMsOjSauuwAzxsodG4nif0iPXc/HbGuI0l37tzRs2fPKhZVUrstf+31yeVpt6VpjgMAK2Ipo6emEfs4tt39XtH62GHelaSNjQ1fW1urfCz/+CN9886bEz93q9eTffJp5eMAwKpYRkd4UT9Eci1D0p6k92ZUFgBACYsOjWMVB0RT0ktNTcPMbE/cnwEAS7PQ0IgX+34cRZXXcPfeuG1jX8V+7BjPlo0bcQUAmLFlNE/tKTeE1szWJfVy71tmdpAPlhgOx1lgmFmDwACAxVt4R7i7d82sEy/6DYVpRPJ3g7cUbvZrKty/0ZJ0JElmNry7oqG4AIA5WcroqfyUIAXresqFQaxdvJQWAIDFY5ZbAEAyQgMAkIzQAAAkIzQAAMkIDQBAMkIDAJCM0AAAJCM0AADJCA0AQDJCAwCQjNAAACQjNAAAyQgNAEAyQgMAkIzQAAAkIzQAAMkIDQBAMkIDAJCM0AAAJCM0AADJCA0AQDJCAwCQjNAAACQjNAAAyQgNAEAyQgMAkIzQAAAkIzQAAMkIDQBAMkIDAJDs1WUc1Mw6ks7j25a7P57HNgCA2Vp4aGQXf3c/jO9bZrbv7tuz3AYAMHvLqGlsu/u97I27982sPYdtAAAzttA+DTNrSFovWDUYFQJVtgEAzMeiO8JbkgYFy8/julltAwCYg0WHRlOXndl5A0mNGW4DAJiDpYyemqfYad6Jb//XzP5zqh3+/N+Gl/yhpN8UHHiqw9Rc8TlfbzftnG/a+Uqc8zjfGbViGaHRLFg2qcaQvI27dyV1yxYqlZkdu/vGvPZfR5zz9XfTzlfinKtadPPUsYov9k1JpzPcBgAwBwsNDXcfSOrHEVF5DXfvzWobAMB8LGMakT1d9jnIzNYl9XLvW2Z2MBQSY7dZsLk1fdUY53z93bTzlTjnSszdZ1GQcgcNndV9hWanK1OCxHsvDiTdc/d+yjYAgMVYSmgAAFYTs9wCAJJdu/s0pnETZ9+d4pwl6a5Cc+FuHLCwEqb9npnZgbtvzb5k81P1nM1sR+FG2nNJyiYNXQVT/mxL4We7uwo/27EP+H1J91N/Niv/Hrg7r9BE15G0mXvfkrQ/623q9Kp6zkPvNyWdLftc5nnOQ9uvh1+b5Z/LvM9Z0pHCKMXs/Yv8+zq/Kv5s7wyf3yr8Psefyc34Opnnz4S7Exq5L9pLX+xJF8Mq29TpVbb88Qdrr2D5i/wPYJ1f037P4i+mL/s85n3O8QI6/AdCa9nnMudzPihYtrdCQbleIjQq/x7Qp6GbOfvuFOXvFCw7V/Fd+7Uy7ffMzDZ9hZpnpKnO+aGGhrV7bjRjnU1xzq2C9Q1fgeapMqb9PSA0gps4+27p8rt7391vj9jX8QzLNi+Vv2dm1lIY8r1qSp9zvKg04v83zaxtZjsFN9jWVdXv866kIzPbky7a/PdnX7ylm+raRWgEN3H23ZmUP/5i9dx9FaZ0meac11fkHIdVOeeNbL27H3qYeaGrcP/UKqj0fY7neU9Sx8xcUn9Fv+eTTPW7T2igsvjX97a73192WeYpVtlv2pQ1DeVqVrGJphlnY7iW4s/zhqS3JT1WqHUUNcfeaITGpbnOvltT05Z/T9J7MyrLopQ656xJZsXbtct+n/tS4TmfS6p9f11U5Wd719277j5w912FWsfeKvRRVlD5d5/7NIKbOPvuVOWP7b4rdX+Gqp1zR7qY7+xCdv+Ch6n466z0Obt730Y/H2YVvt+lzzkGw1F+mbufmtmWpPu6XjXNqX73CQ2Fv6jMrG9mwyMlxs6+W3abOpmm/FkHoV+dG6xd9/Ou+H1+6YYnM9srWl5HU3yfT82sNTRiaiUGPMz4d/NYxSONVta0Xx+apy6t+uy7VZQ+5/gX2XF2MTGzxopV36t8n1ddlXPeja/8NqvUMVzqnOPF8kHBfjpandlwC4e9z/raxYSFOTdx9t0y5xw7Cs9G7Or2qjRVVfk+59Zt6fJCclD32lWm4s/2pi6HYL4R2/lXRtlzjhfVh7r8GW9IOqz7/Snx93JToRmtrdCJf5Y1nc762kVoAACS0TwFAEhGaAAAkhEaAIBkhAYAIBmhAQBIRmgAicxsP05kN+4zB2Z2kri/AzNblUkAAUmEBlDGgXRx/8Iom7qe02kDkggNIFm8kW+g4juH82HyxcIKBSwYoQGU01WoTRR5oPBskZW4Mx6ogtAAynkqjWyiomkK1x6hAZQQJ+zrS9rOL89CJP8M8fiY1BMzczN7kT1GdJSijvG47ebQsr24v4n7BGaNqdGB8vYVHs6Tn1p6W9Lh0OdaCjPF9uP/j8zsYiK5KmKorCs8/KopaX/afQJlUNMAysvC4f3csraGmqbiU+B67t6PneiHCrPkVpKbzXTL3U/jPnc1VOsB5omaBlBSnCb+VCEAunGK6UHRNOmxaemBQk1j2uetZM8tORnzZD1grggNoJp9XdYstlQwzDZ3k98jhcdo7mn6Z8gP3P32lPsAKqN5CqjmC+miJpE95OZCfBLaukJTUpkH+eSfkjgcMD1JjeHnlQOLRGgAFcQO8J6kJypumspCYjc+bnNTo+/vyGvFzzfivvPH7Cs+MTCOzGqZ2SZTkWCRCA2gugOFmsFLI5diqDxWeDTsicKjOA8V7igf5ZHCiKizuO9HCuFzsY27b+uyaexM4fGkT6c/FSANj3sFACSjpgEASEZoAACSERoAgGSEBgAgGaEBAEhGaAAAkhEaAIBkhAYAIBmhAQBIRmgAAJL9P1bRZ+u69scCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_groups_pointwise(x_axis_0_test,\n",
    "                      x_axis_1_test,\n",
    "                      A_test,\n",
    "                      Y_test,\n",
    "                      Yhat_out_test,\n",
    "                      y_u=None,\n",
    "                      y_l=None,\n",
    "                      point_color='r.',\n",
    "                      shade=None,\n",
    "                      filename_0=\"synth_baseline_0.png\",\n",
    "                      filename_1=\"synth_baseline_1.png\",\n",
    "                      save_figures=save_figures)\n",
    "\n",
    "p_val_base = utility_functions.fair_dummies_test_regression(Yhat_out_cal,\n",
    "                                            A_cal,\n",
    "                                            Y_cal,\n",
    "                                            Yhat_out_test,\n",
    "                                            A_test,\n",
    "                                            Y_test,\n",
    "                                            num_reps = 1,\n",
    "                                            num_p_val_rep=1000,\n",
    "                                            reg_func_name=\"Net\",\n",
    "                                            lr = 0.01,\n",
    "                                            return_vec=True)\n",
    "\n",
    "\n",
    "plt.clf()\n",
    "x, bins, p=plt.hist(x=p_val_base, bins=30, color='salmon', alpha=0.7, range=(0, 1), density=True, ec='red', align='left', label=u'Baseline')\n",
    "plt.grid(axis='y', alpha=0.75)\n",
    "plt.xlabel('Value')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Fair Dummies Test')\n",
    "for item in p:\n",
    "    item.set_height(item.get_height()/sum(x))\n",
    "plt.ylim(0, 1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# step size\n",
    "lr = 0.01\n",
    "\n",
    "# inner epochs to fit adversary\n",
    "dis_steps = 40\n",
    "\n",
    "# inner epochs to fit loss\n",
    "loss_steps = 40\n",
    "\n",
    "batch_size = 6000\n",
    "\n",
    "# utility loss\n",
    "cost_pred = torch.nn.MSELoss()\n",
    "in_shape = X.shape[1]\n",
    "out_shape = 1\n",
    "\n",
    "model_type = \"linear_model\"\n",
    "\n",
    "# equalized odds penalty\n",
    "mu = 0.99\n",
    "second_moment_scaling = 0\n",
    "\n",
    "# total number of epochs\n",
    "epochs = 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "fair_reg = fair_dummies_learning.EquiRegLearner(lr=lr,\n",
    "                                                pretrain_pred_epochs=0,\n",
    "                                                pretrain_dis_epochs=0,\n",
    "                                                epochs=epochs,\n",
    "                                                loss_steps=loss_steps,\n",
    "                                                dis_steps=dis_steps,\n",
    "                                                cost_pred=cost_pred,\n",
    "                                                in_shape=in_shape,\n",
    "                                                batch_size=batch_size,\n",
    "                                                model_type=model_type,\n",
    "                                                lambda_vec=mu,\n",
    "                                                second_moment_scaling=second_moment_scaling,\n",
    "                                                out_shape=out_shape)\n",
    "\n",
    "init_seed(seed)\n",
    "\n",
    "input_data_train = np.concatenate((A[:,np.newaxis],X),1)\n",
    "fair_reg.fit(input_data_train, Y)\n",
    "\n",
    "input_data_cal = np.concatenate((A_cal[:,np.newaxis],X_cal),1)\n",
    "Yhat_out_cal = fair_reg.predict(input_data_cal)\n",
    "\n",
    "input_data_test = np.concatenate((A_test[:,np.newaxis],X_test),1)\n",
    "Yhat_out_test = Yhat_test = fair_reg.predict(input_data_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fair Dummies (All): Test Error = 3.400726861346923\n",
      "Fair Dummies (A=0): Test Error = 3.4092725906435017\n",
      "Fair Dummies (A=1): Test Error = 3.318214246942538\n"
     ]
    }
   ],
   "source": [
    "print(\"Fair Dummies (All): Test Error = \" + str(np.sqrt(np.mean((Yhat_out_test - Y_test)**2))))\n",
    "print(\"Fair Dummies (A=0): Test Error = \" + str(np.sqrt(np.mean((Yhat_out_test[A_test==0] - Y_test[A_test==0])**2))))\n",
    "print(\"Fair Dummies (A=1): Test Error = \" + str(np.sqrt(np.mean((Yhat_out_test[A_test==1] - Y_test[A_test==1])**2))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAElCAYAAABAu+8OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXxU5b3GnzcbqyQZ1kDCkiguJSIIiAIiGCpqq4JIpcJFaWUrVhAUCJUqNsHgUqFaJbap3lyUewtC1dZWopZbtKBhu6EiVIIyA2GRkCCLgSTv/eN3DnMymSwkkzk5M8/38+GTzDln5rwzDA+//VVaaxBCiJ1E2L0AQgihEBFCbIdCRAixHQoRIcR2KESEENuhEBFCbIdCRGxHKbVBKRVn9zqIfUTZvQAS3iilxgFIA+ACUBLE+04FUGw8TNZaLwvWvUl1aBERu/mR8TM5WDc0RUhrvUZrvQbAGqXUymDdn1SHQkRsQyn1GID5xsNgumbTDAECAGitCyFWGbEJChGxBaVUMoASQwSAIFlERiyqv59TJUopipFNUIiIXUzTWmcbv5cASAnSfZPhPxZVjCC6h6QqDFaToGNYHp9ZDhVCgtV1Pa/ecRyt9bQaTrngDVJbKUFw3UNigUJE7OAeH6GolzVSi7gQh0MhIkHFyFi5lFJZlsPBdon8WV+0hmyEQkSChhEoTtFa3+Nz/DiALP/PCjj58C86LgDbgrQG4gOFiASTqVrr+X6O16uQUSn1x3rep7gmN05rXaKUKlRKxWmtrfeN01rn1fP1SYChEJGgYFRQF9ZwutC4JtmSzq+GryXVCLIATAWwzLhvfwAUIRth+p40KUqpZMOS+SOAalaKETMyraSVwajlMcoGSpRSaWaLCQPh9qI4s5oQYje0iAghtkMhIoTYDoWIEGI7FCJCiO0wfW/QoUMH3bNnT7uXQUjIce7cORw5cgRHjx79Rmvd0d81FCKDnj17Ij8/3+5lEBJSeDwepKenAwByc3O/ruk6umaEkCbBKkIZGRm1XkshIoQEnOPHj1cRoaSkpFqvp2tGCAk48fHxuPnmmzFy5Mg6RQigEBFCAojH40FUVBS6dOmCyZMn1/t5dM0IIQHBjAk9++yzuNjWMQoRIaTRWAPTDz/8MJRSF/V8ChEhpFH4ZsfqExPyhUJECGkUr732GoCGixDAYDUhpJHMnj0bpaWl6NatW4NfgxYRIeSicbvdeO6551BWVoa2bds2SoQAWkSEkIvE7XZj0aJFAIDi4mIkJCQ0+jVpERFC6o1VhDIyMgIiQgCFiBBST3xFqKGBaX/QNSOE1IvKykrExsbiscceC6gIARQiQkgdlJSUIDY2Fj169MDy5csRERF4R4quGSGkRtxuN37+859j7dq1ANAkIgRQiAghNWCNCQ0ePLhJ70UhIoRUwypCmZmZSExMbNL7UYgIIVUoKyvD4sWLAQRHhAAGqwkhPrRo0QI//elP0aNHj6CIEEAhIoQYeDweFBUVYeDAgRgyZEhQ700hIoRcGOURHR2Nvn37IiYmJqj3Z4yIkDDHOk/oySefDLoIARQiQsIaqwgFKzDtDwoRIWHMxo0bAdgrQgCgLnbIdagyYMAAzZ1eSbigtYZSClprnDhxAi6Xq8nvqZTaqrUe4O8cLSJCwgyPx4O5c+fi0KFDUEoFRYTqglkzQsIIa0yooqLC5tV4oUVESJgQiN02mgoKESFhwKFDh5qtCAF0zQgJC+Li4nDVVVfhvvvua3YiBFCICAlpDh06BJfLhdatW2PBggV2L6dG6JoREqJ4PB4sWLAAK1assHspdUIhIiQEsQamJ0yYYPNq6oZCREiI0ZyzYzVBISIkhNBa45lnngHgHBECGKwmJKRQSuHRRx+F1toxIgTQIiIkJPB4PFi9ejW01khMTHSUCAG0iAhxPNaY0C233IL4+HibV3TxOE6IlFJTjV9TAMQBmK+1LvE5X2w8TNZaLwvyEgkJGr7zhJwoQoDDhEgpNVVrnW15PA7AVogoXRAhrfUa43GyUmql1nqaLQsmpAlpLkPNAoFjYkRKqWQYgmNiCI7LECQAmGaKkHG+EEBa8FZJSPDweDyIiopyvAgBDrOIAEwFMN/nWDFEjOIA9PfznBKlVJrWOq/JV0dIECgrK0OLFi0wePBg9OvXDy1atLB7SY3GMRaR1rpQa+3PAU4GkG/8LPFzvtg4R4jjcbvdmDFjBjZv3gwAISFCgIOEyB9GTChPa70NgAveILWVEkhQmxBHY24DXV5e7nhXzBfHCpERM5qmtR7ViNeYqpTKV0rlHzt2LICrIySwBHsv+mDjWCECkAXgZp9j/obv1mgNaa2ztdYDtNYDOnbsGNDFERIoTpw4EdIiBDgvWA0AUEplwad+CBIn8ic6LgDbgrIwQpqAuLg43H777RgyZEhIihDgQCEy4kIrjdS8eSxNa52nlCpUSsX5CFQcM2bEibjdbgBAUlISfvSjH9m8mqbFUa6ZUioNQL4pQkqpOOOYSRYkxW9e3x8ARYg4DjMm9OyzzyIc9h50jEVkBKc3GL/7no4HJOZjBKDTIG5aMquqidOwBqbnzZvn7/secjhGiAwrqM6/EWsLCCFOwypCTpon1Fgc5ZoREuqsXr0aQHiJEOAgi4iQcOChhx5CcXExunbtavdSggotIkJsxu12IzMzE2fOnEHLli3DToQAWkSE2Io1JlRSUoLWrVvbvCJ7oEVEiE34BqbD0RIyoRARYgPhmh2rCbpmhNhAZGQkOnTogDlz5oS9CAEUIkKCyvHjx+FyudC1a1c899xzUEqhoABYtw4oKgISEoAxY4DUVLtXGlzomhESAAoKgCVLgBkz5GdBQfVr3G435syZg1WrVgGQDoGcHGDBAqBXL+CBB+TnggVATk6Q34DN0CIipJHk5ABr1wL33guMHg3s2SNicvfdwJQpco01JjR8+HAAIlZr1wK5uYDLGGAzaBBw++3AxInAwIHhYxlRiAhpBPURk7g4/4HpdetEvFw+U7RcLmDCBDlPISKE1EldYrJ27Tm43U8AqJ4dKyoSC8ofvXsD774rbl44xI4oRIQ0gl27AI8H+OST6mLRuzfwz3/GYMaMGejcuXO17FhEBJCVBXTqVP25r74KbN8O/OAHNbt7oQSFiJAGkpMD7NgBjBwpgWarWIwa5caHHx6AUkPwl78MqGbV5OQA//d/wLffAtOnA4cPe5/buzfw9tvAW28BQ4fKvUI9dkQhIqQBmLGhnBxg6VLg0kutYuHGhg2LUFoaiVOnrsXgwS2rWDXXXQds2QL86U/A+vXAr38tbtzDDwOPPgqUlAB33OEVIZNQjh1RiAhpAGZsaOhQsWImThSRcLnciIlZhJIS4OjRJdiwoWW1IPZ114kF5XKJmzVwoLzeP/8pgrZnDzB1qv/7irsXvPcZLChEhDQAa6DZFJM33nBj927Jjh0+nIFHHknyG8Tu3h346ivvsdRUr4WzZQvw0EMiRoMGVb/v3r1Aly6Bfz92w4JGQnyoT3FiQoKIhUlqKjBgwKdo1w4YNSoDUVFJuPxy/6+fkgLs2+f/3N69ImqrVwPFPtuFFhcDb74pcaZQg0JEiIX6VjqPGeMVi8rKSgDA2LFjsWTJCvz5z0kYNKiqUFnp2RNwu2sWmunTve5ebq5YSbm58njcuNCLDwF0zQi5wMVUOqemmmLhRvv2yzB+/FyUlPTEm2/GYdw4uXbBAnmu1T0rLgY2bZLXM+NKvXuLJfTmm16hSU2tGjvq0kVS/aEoQgCFiJALXGyl86hRbuTlLcLp08Dbb0chKamqWFiD2L5iM2UKLjS71iQ01thRqEMhIsSgrkpna7bK7B2LigJWrPC/DfTAgcDWrcCKFYBS8tgqNuEkNHXBGBEhBr4BaCvWbFVRUVGde9GbsabBg4GXXpJM2FdfAZ991kSLdzi0iAgxGDOm5rjOm2+KNQMALpcL11xzDcaPH+9XhNhVf/FQiAgx8Aag/cd1XK6D+PTTWPz1r21RVPQIcnL8N6Kyq/7ioRARYsG30tkMIsfHezB7djq++643xo//Ra2NqBcTayIChYgQH3yDyB6PB489lo7vvgOWLr0fffrI8ZrcLTPWFE6V0Y2FwWpCasHj8SA9XURo9OhM9OlTNSZkdbdMrMWOVkK5Mrqx0CIiIUMghtAXFACvvAJ8+ikAaPTo8QJiYoCIiEzs3JmIGTNkjhAAVFbKfVJSZIyHSV2xJsaHqkMhIiFBfeZG1+c1XnhBhOauuwBA4Z135qGkpBynTiVizBigfXuZnAjI0LJevYBnnpG2DSs1xZooQv6hEBHHE4h0eUEB8Ic/iGA8/7wHW7d+gEmTJiEtrQvuugvo2lVqgNq3Bz78UJ4zcaKIXocOwPHj8hrhWhndWBgjIo6nPuny+rxGp07AnXd6sGxZOj744AN88skJLF4sx8vLJcZz5Ajw5z8D//430KYNcOedQEWFXPPKK03z/sIBChFxPEVFqHHkRu/eVeM3tb2G1h5s2pQOALj22gwsXdoeSgFjx4rofPONDDXbv1/cvi1bZJCZ6XK9+2747UcWKChExPHUtzWjNuLiPCgvT0d5OTB5cgbefz8JubnA8OEiQN27A9dcI+Lzve8BrVsDCxeK+zdokASsJ0wA1qzxP7+I1A6FiDie2tLl2dniTtU25AwAEhKO4dy5GGzfnoEnnkjCkCHi2o0cKeJy8CAwezbw3XfAI4/I0PvnnxdRMtPy991Xf1eQVIXBauJ4akqXL1smXe+DB4vr5i+TdvbsWbz5Ziv87W/90K3by/jyy2hs3y4Nqps3A19/DbRsCZw5A7z3HnD6tAhbZSXQrx/w1FPAk0+KSKWmynWsnL54KETEkfirGbKmy5UC4uNll4yaMmnx8R7Mn/8LlJY+gNzc4XC5olFQAPzkJ8DRo8DGjcCttwKLFomwLVwIJCcDWouwJSaKRbRpk7wewMrphkLXjDiOmsa5fvYZsHgx8NvfAp07Aw8+WHMm7Y03pGL67NlK3HpryoXrUlOB3/9eBObpp4GrrpJj//oXMGcO0LatzBdq0QL41a9kK6GpU0UAWTndcGgREUdR35qhmhpPCwqAjz7y4PDhdLRsKRXTAwdWbdswXb3f/Q5o1UoC1Bs3AmfPAj/9qZwrLfW6gufOyfnNm1k53VAoRMRR+NYMWV20Vq2kluell/w3nkr19UlER6cjJga4+eZMLF+eiOzsqtcVFMg20hUVQGGhBKW1lkpqM7ZkrZzeuFHOs3K64dA1I47CWjNkddF69AC2bZMs2ZVXigtlzaSZltSKFe1w6tRYdOqUiV27EjFsmGztvGlT1dfs2BGIjBT37tQp4KabgI8/rpqZS00FZs0SN235copQY6BFRByFaem0auV10ebOlWrnPn1EhHbtkvjNVVd53aePPvKgS5cy3HVXCpS6C5MnezNp778v7tbDD8vrTJwIvPOOuGFTpoj4TJwo2Tc2sjYNSmtt9xqaBQMGDND5+fl2L4PUQUGBWCxXXy1C07q11Ai99x7w+OPiHu3YAeTni0gsWAAcO+bBzp3piI6+BCdOrMD69ZFVgtjvvw9MngwkJQGXXCJFjL6d+7m5UlE9Zoy4Y4cPS3asIR3+4YpSaqvWeoC/c7SIiGMw40FnzgAvvyyWyWefiSX0+ONey8Ss5Rk+HMjO9uCKK9IREQGcObMQo0dH4sUXq6b9jxwBbrgBOHBAsmD+BpqZkxXZyNo0MEZEHMGSJeI+bdwIXHYZcMUVIkq7dwPt2oklZAaS9+4FDh0CCgs96NUrHa1aAf36ZWLjxkS8+mrVtP8jj8gYj7Fjpfaosa0ipGFQiEizZ8kS4MUXpY5n6VJg2DAgOlpGckRGSif8b38r123aJAHr4mKgV6+1iIwEnnkmE+PHJyIuTvrG3G7JcploLRXU4bjnfHOBMSIDxoiaJwUFYgnNmSOxIJNNm4DbbpOWC0CyV2VlkgHr1k36vp555jyyso7j/vu7YMkSICoKeP11GdmRmioWkHmPf/1LBKy0VHrLrAHp7GwRvYSEhk9+JLXHiOq0iJRSPQO9IELqy7p1EkS+9lrvsZwcabeIj5fYTuvWYjG9/TZw5ZVutGjxJJ544ltcf3007r9f/KmiIiAtDbjnHokFnTsnQe2KCuDECREnc4RHVpYEpl9/XTJzSokYWqu4Oe4jsNQnWJ0F4EdNvRBC/LFrl2So7r1XrJKUFIn/DB0KlJQAf/ubxIuk+NCNVq0WoVUrICmpFOfPX4KcHIkdmWn/yy8Hrr8e+N//ldfo2lWsn/JysarM6uzFi70Zutr61WgZBYb6xIjuUUr1bfKVEOJDTo4Ep48fF7dq5Ejg738Xa2XdOhGhM2fEWomIcOOGGxYhJga4884MnDmTiBUrvPOBxoyRFPzrr4ugPf64vP6sWeLmjRlTfaJjICY/kvpRHyH6AMC9SqlrmnoxgUApNVUpNc7485jd6yENw5whnZwsNUMPPSTHevcWUdq3TyylmTOB0lIPevVahHbtgC++yMCJE0no10+e07078JvfiBtWVCRd9V98AZw/L8I0cWLVgkTrRMdATH4k9cOvECml2pm/a61Haa0XAkhRSo0N2soagFJqKoBirfUarfUaAGuUUivtXhe5eMwZ0jNnilXy3nviSp05IxYSINXVs2YB33zTAh07dsGwYRlISEjCc88BMTESz/nmG+CTT8SKeuMNYNUqyZKtWiXHrGl/oGqaPhCTH0n9qCZESqmbAVSLbGut1wLYr5SaF4yFNZBphgABALTWhQDSbFwPuQgKCiQFP2OGtFiUlopFMmWKCEZ5ucSIKiuliLFv32OYOLECxcUd0aZNFp5+Ognbt4uAzZsnLRlffilzpRcv9hYjLl8uw85mzaoa4/FN03OjxOBRJVitlPopgHyt9Q5/F2uttyulCpVSTwPI1FqfDMYi64NSKg5Afz+nSpRSaVrrvGCviVTHrI7etUuslQ4dcGEL5y1bvPuS/fvfMlZj8mSJx4wZI2ICSAbN7XYjIWERtB6O7dt/gnfeUYiLA37+c7GGBg0SAXr99epb/dR3A0RulBg8fLNm+wEkA/ArRACgtS4FsEApNU8ptUFrvdN6XinVU2v9VcBXWjfJAEr8HC82zhGbMTdBTEiQGdCpqSIQ5eUSMP7FL4BJk+S6sjJxfYqKJNVujni96y7gzBk3+vVbhIgIICbm++jXT6yn6GixpB57TOI/b74J/PjHcnzduqrCUd8NELlRYnCoIkRa6w+UUr2UUmO11m/V9CSl1EgAgwBkKaXyAGwAYIrSNAALm3LRNeCCiI4vJQDi/D3BiClNBYDu3bs33crIhTEcCxdKdfQ770j2qbhYtuiZOhX4619FlF55BRgyRLJbmzfLgLL/+A95XlaWG716LUKLFsD58xkYMyYJvXsDx45JQPryy2Vch1UwtmzxP0e6vn1j7C9reqrVEWmt90MsIwCAUmqp1nqhIT73ABgPIB7yD/xViABMB7BMKWWWadshRBeN1jobQDYgldU2LyekMVPhH35YNSXucklm6+hRqenJyJCiwylTJFD81VdSN7RqFXDqVDkGDVqCli2B3/wmAyUlSVi3TvYTc7tlhhAgMSarcDCw3PypT/p+mlLqOIA8SGHjHwGM0lq7tNbTtdbjtdaXQsTpR7CImA24/Bzzaw2R4GKmwv2lxFNSgJ07pYG1Y0cRoUGDZB8xraXy+corgblzo+B2P4Ti4gxs2JCE1FQZYH/qlAwwa9NGLCpr5TMDy86gPpXVcQDWAFiptf6gpouM2NEapZTfXpIgkA//ouMCsC3Iawk7/O2qYVolBQVi3SxcKMKSl1d11EbPnuK2/fCHMu1wzx4JZr/6qqTYJ0zw4Ntv9+Ktt0ZCqasxYYLUGL3/vrhcEybIcwsLpWJ6wgRx5w4fltgTA8vNnzqbXpVSr2itp9f7BZXqZbh3QUcptQ/AtVrrEusxrXVKXc9l02vDMYPQ997rnXq4erUElwE5N2SITD8cMUI65QcPlvGucXEiFqdPi+hMmiRiUl4uNUM33uhBWVk6zp6NwKlTv0VxcWucOiWZscpKYPx4sahyc0UEn3hCGljNOdIc4dp8aOxgtKyLuZldImSQBQk+LwMApVR/iEtJmojadtW4805pGDV7tbp0AX79aylE3LhRxOjAAdk99dw5aeHo0EEyYGfPAsnJHpw+nY7z54Evv3wKjzzSGocPy3Y+SUkiPGaOobJS7tW7twjgli2SuqcIOYM6Y0Q2C8tFYQSfS5RSaUqpcQDStNbT7F5XKFNbP1anThLzMc/FxgInT8rx+HhptejfX+p0unWTzNesWVJH9MQTHrhc6Th7Fvj3vzPw5JNJmDJFgtff+564ftHR3grpX/8aePRRb/8XA9TOIuQGo2mts7XWeUabxzK71xPq1NaPZfX6zdEdt94qgnHZZSJKZWVS6zNnjgjWxImSQfv44x1QCti8OQM9ekiKHpAes27dxJ179tmqFdJm/xcD1M6Dg9EMGCNqGEuWSE/XpEnVz6WlSeo9KUmG2bdsKVZMWZlkuxISZDJibCxw440iINOmVWDOnEj8859AWtq3SEi4BFu2iNs1cqT0i8XGAvPnS3YtJ8c7yGzfPnEVz56VALW1h4zYT20xIgqRAYWoYZgze6wxIkA63n/5Swkqay0xnOJiGe2akiK/t2wpQ806dJAB9tdc40Z5eSY+/3wOkpN749tvpdAxIkKC33l5Yvm0bFl1RlBBgQjUG2/IJojTpzM21BzhLh6kyfDXj/XBB8ALL4iIrFolrlZ8vLRKbNsmcR6tpfbnT3+S7FhxsRsnT5rzhFqhTx8Z0bp2rRQ6DhwoA83Ky8Wa8u3/2rlThI9WkDOhEJFG49uP9cUXwOzZ0k/Wvr1319Rvv5X4zuHDcvyLL+Tc7t1uDBy4CAAQE5OBu+9OwqBB3t02YmIk8Nyjh7RvLF7M/q9Qg0JEAoK1H2vGDIkP/exn0qT6+eciSmlp0jvWqZO4au+9B7RqdRgDBixCdDQQFZWB4uIk7NkjfWf/9V8SUyoslCrrZ56R4kff+xHnQyEiAcd3oFifPrL32KZN0pjas6ek6rUGbrihPY4fH4T9++/ED36QhOho4KmnpK1jyRKptF61SkbElpVVH+lBQgMKURhSWztGQ58PyLGPP5b2jBMnJLv16qtSF5SXJ/uRffaZ9JTFx3uQkHAJPv44FhUVszB0qBQhXnmlFDR26gS0bSuZsLg475iPW2+tPtKDOB8KUZhhbccwBeLuuyV2M2JE3aLk+/w9e2QPMUAEo6AA6NdPAsjHjknF9H/+p2TG1q4VK6hNG5kndPZsDyQmPoW+fSXj1aKFuGF9+gA/+Yk3G/aDH0jrxo4dwD/+ITOrSWgRcgWNpGas7RiTJonl8vHHUkxYVib/wGvbs8v3+a1aScznxAlJv+/cKfVCH34o5ysrxQ0rL5fXb9ECaNvWjeuuk6FmN900FcuWSeC5UyfZHujrr4FLLpF+MzMT9tJL3oLFfftYMR2K8P+WMMLajuHbI9a2rbRKmDtbxMZK86jV/bI+37SMyspkJOunn4rgrFol165dK5ZScrLcKz8fGDbMjSNHFiEyEkhMzEB0dBJef12E5b//W2JCeXlyz88/l3sMHepd/7ZtMneIFdOhB4UojCgqEncKqN4j1ru3pMJdLqnbeeQRYIBRenb0KLBhg1hMo0d7RWz0aCAzUzY4BETMXnhBAstJSVItnZ8v1+3erdGt229x/Djwr39loEuXpAszqAEpcDx4UBpmo6JEDDMzvbVC27aJdeQ78J6EBhSiMMLMZg0aVFWUAG+TaEGBCEllpaTezbEeubniKuXlieUzZIgEolu3lmu++kraMFq3FnftttvESuraVQLUnTopzJs3D48+ehbbtydi0yZ5TetA+quukm76xERp2xg6VGJCr70mXfqzZqGKeJHQgTGiMMK6PY41xW5tEn3lFdl88L77JM4zaJD8XL1a3LWVK8WF27ZN3LKbbpI5Q+Z+Y3Pnymt+9RXw1ltAly4ebN2ajRkzKtC+fXt07ZqI0lJg0SLv/vL79wPp6dKXNmZM1e2DoqKA4cPltShCoQstojDC2o4xdKiIzqlTMsjenGL4ySdSyWxmwqyp+m7dxHrZtEmqo+PiJKXeooVce/689HnFxADbtwOXX+7B2rXp6NwZiI0dg5df7ojsbGDUKIkJTZgA3HKLvGZmZvWtfOiChQ+0iMIMq7UREyNzfIYOlRk/ubmSPk9MFPEZMUKEKypKWi3i4mSSYkqK/Dx9WkQkKwv4y19kp43KSvlzySUedO+ejpYtgU6dMpGV1RHPPy/u1fr18hyrReS74yoJL9h9bxCO3fcFBWIVffqpxHcGDpR4zhdfyEaFH30k1tO774rl8vvfS4V0fr50xJ87JzU/Wktg+9Ah6SH7xz88uOaadLRuDfTokYnS0kR06SJi55uJo9UTPrD7nlTDWpg4ebLEi7KzpZYnKUkKCb//fZmgCACPPy67ZJw9K/VDKSlSO3TwoFhXBw6IaG3eDChVivj4VnjqqceRmJh44X45OVULIc1NE2kJEVpEBuFkEdU0Q2jhQnGx7r4bePllScub53ftEgsoNlYyW7t2STd8YaEULe7YAXTseBoREW2QmAgMH16Bw4cjkZAgllBOTvX7FReLeLFzPjyozSJijCgM8S1sXLJEOubz8oAbbpBrxo6V8507A337Ssbq9GkRnz17JEuWng5Mm2a6Wh5ceulMxMa+j/JyIDk5Eg88IKM8Fi4UV8zfXOsJE7xzpkn4QiEKQ8w50zk5YhmZc3+6dxcx+ugjyX599x3wq1+JCFVWynN37RLXbeZMCXQnJQFPPeVBz57p0FqjsvIqrF9fNfU/ZIhUShcUVF+LOWeahDeMEYUhCQkiOB9/XNVd+vGPJTW/dasEsa+7TlynH/5QrouJkYD0ww9LMFuaUj34/PN0AMDu3ZnIyEisZvn06iUpfn9d89xtgwAUorBkzBiJA5k7ZwBiHa1eLWLTooVk0mJixDV7/nlJ6cfFSbbs7bdl9Gv37qewb98iVFYCBw5kIiYm0e+OHmPGiIV17lzV41+pvvwAAA12SURBVGYhZdZF7ZxHQhEKURiSmirFibm50h8WESGbFnboADz5pAjR8uXivv3jH8DTT4twLVki2bE9e6QWCGiLiooJOHmyD264IRGFhd4WEt/7XXmldNT7tnVwO2gCMEYUtowYIXN+9u8X0UlJkZjPlCkyigMQ16uiQiyWggKxbIqKgHnzPEhJ2Y3OnYH+/UcjIyMRhw7JaFizhcRKcbG4dEuXsoiR+IcWUQhS2wRG89yuXZJyz8mR6x54QK75zW9EMLp2lYzWwYNy/X33yUD873/fg9/9Lh2XXNIaDz30Er78MvJCe8bdd8skRd8dNkzL5+675Q8hvrCOyCBU6oishYpm5/zq1V4BsJ7LzpZ4T1KSBKRjYmT31IQEcclcLnGl/v53yaRp7UHfvumIiQFSUjJQWpqELl2qV0ibYnf4MPyeJ+EJN1isB04SoposHrNQceFCmZJonu/RQyqjz5+XKur77vMKw6ZNki0rLxf3rEULuWbSJHGphg+XwPRll3lw5Eg6IiKAigrZ8oduFbkYWNAYQvjW/vTq5R3vum6dCM/Spd7zBw7IXvMREVIp/b3vVR0HO3QokJEhAeavv5YWj/PnxRK6807pQVu/Huje/V1ERACdO2cgNzcJa9b4rwsipCEwRuQgfMe7AiIgt98ucZmKCu9YD7NquqhIJi8OGSKbHE6a5L1+4EDvLOiuXSWT9tFHMu51+HDg6quBwYPltXr2fBA7d45FUlKnKhXRdLlIIKBF5CB8x7uamMKwd68Ig3nevD4lReJAhw5Vvd5srTCLCqdPl+75yEgZ11FW5sGnny5CYeEJ5OZG4tixThfmRbMimgQSCpGDMFsz/NG7twhIQYE3fW5eX1wsxYQlJd5zppBYpzOmpoo7V1QEjBrlweHD6fj6azfuuec0jh71ZtYAVkSTwELXzEFYZ077snev9IolJ3vT5xUVEgs6cEDcsG++8Z7bt0/GeEycWLWocMoUoGdPD1asSMe5c8DWrZm47bZEzJ7tvYYV0STQ0CJyENaZ01ZMYfjZz8SaSU+XgsHSUsmezZolbtlTT4l4fP458Ic/iEXjW1To8Xjw2mvpaNcOyM6W3rGiIqk52rLFu90QK6JJIKFF5CCsM6drKhgsLfVuw3PLLbLTxv33A3fcIaM79u71blzoL/3epk0b9OjRA9OmTUNiYiKmTBFrat06CXqb4kURIoGEdUQGTqwjqqlg0Pe8OaK1tgLDI0eOoH379ojifs6kiWBBYz1wkhAFGrfbjUWLFmHw4MGYOXOm3cshIQoLGkmNmCIEAHfccYfNqyHhCoUojLGKUGZm5oVB94QEGwYEHEJtHfUNoaKiApmZmQAoQsR+KEQOwNpRH6iteCIjI/Hwww+jbdu2FCFiOxSiZoQ/qweovb/M7BerL263GwUFBbjttttwxRVXBP5NENIAKETNhJqsnujo2vvLLqbx1BoTuvHGG9G2bdsAvwtCGgaFqBlQW1f9NdcA48f7f17v3lJkWB+sIpSRkUERIs0KZs2aAbV11Q8cKCLlj/o2nvqKUFJSUiNXTEhgoRA1A2rrqh87FsjPr7m/zIwj1caePXuglKIIkWYLXbNmQG1d9ZWVQL9+NfeX1RYfKi8vR1RUFNLS0nD99dejTZs2TfcmCGkEtIiaAXV11Ztd8+ZWPJs3yyjYzz6Tvcb8jWx1u92YOXMmdu3aBQAUIdKsoUXUDKirq960elJT/e/S4VtTZI0JxcbG2vSuCKk/FKJmQn3GbdQ1s3rgQCAujoFp4jwoRM2I1NTaYz51zaxevfoY3G6KEHEejhIipdRU49cUAHEA5mutS3zOm5GWZK31siAvsUkpKpJiR3/07g188okLw4YNw+jRoylCxFE4RoiUUlO11tmWx+MAbIWI0gUR0lqvMR4nK6VWaq2n2bLgJqCm7JrH48HOna2QkNAeDz74oD2LI6QROCJrppRKhiE4JobguAxBAoBppggZ5wsBpAVvlU2Pv+yax+PBggXpeOedZ3DXXRxyR5yJI4TIYKqfY8UQMYoD0N/P+RKlVMiIkTW7lpsL/PnPHsyZk47SUiAt7We4+mpl9xIJaRCOcM0M6ybez6lkAPnGzxI/54uNc34x3LmpANC9e/fGLzQImNm1N97wYPdu2Yt+9uwM3HwzY0LEuTjJIqqCISJ5WuttAFzwBqmtlECC2n7RWmdrrQdorQd07NixiVYaeFJTgZiYV9GuHbB8OUWIOB9HWES+GDGjaVrra+1ei13MnTsXJ0+e5FAzEhIEXYiMmM38el4+zXDLfMkCcLPPMZef62q0hpyI2+3G+vXrMX36dLRr1w7t2rWze0mEBISgC5HWOg9AXkOfr5TKgk/9ECRO5E90XAC2NfRezQlr28Y999yDLtx4noQQjooRGXGhlVYrSSmVZohSoZE9sxJnCJ+j8Z0nRBEioYZjhMhw6fJNEVJKxfmk5rNgSfErpfqjEZZXc4FDzUg44IhgtRGc3mD87ns6HpAMmFJqqiFOcZAWD8dXVZeVlaFt27ZYuHAhRYiELNxy2qC5bTl98uTJC8HoiooKREZG2rwiQhoHt5x2GG63G7NmzcLbb78NABQhEvJQiJoZ1phQ//7+ulYICT0oRM0I7kVPwhUKUTPhzJkzePzxxwFQhEj44YisWTjQunVrTJ48GZdddhlFiIQdFCKb8Xg8OH78OPr27YsRI0bYvRxCbIGumY14PB6kp6fjxRdfxPnz5+1eDiG2QSGyCVOEAOCXv/wloqOjbV4RIfZBIbIBqwgxME0IhcgWNmzYAIAiRIgJWzwMgtHiobWGUgqVlZUoLi5Ghw4dmvR+hDQn2OLRDPB4PJg/fz6OHj2KiIgIihAhFpi+DwLWmFBZWZnNqyGk+UGLqImxihDnCRHiHwpRE3Lw4EGKECH1gK5ZExIbG4tLL70UDzzwAEWIkFqgEDUBhw8fhsvlQtu2bbF48WK7l0NIs4euWYDxeDx47LHH8NJLL9m9FEIcA4UogFgD0+PGjbN5NYQ4BwpRgGB2jJCGQyEKAJWVlXj66acBUIQIaQgMVgeAiIgIzJ07F1FRURQhQhoALaJG4PF4sG7dOgBAr169KEKENBBaRA3EGhMaOXIkYmNjbV4RIc6FFlED8J0nRBEipHFQiC4SDjUjJPBQiC6S/fv3IyIigiJESABhjKienD9/HtHR0Rg2bBgGDhyIli1b2r0kQkIGWkT1wO12Y8aMGdi6dSsAUIQICTC0iOrAug10586dbV4NIaEJLaJa4F70hAQHClENHD9+nCJESJCga1YD8fHxGDVqFEaMGEERIqSJoRD54Ha7ERUVhYSEBEyaNMnu5RASFtA1s2DGhJ599llwvzdCggeFyODcuXMXYkKzZ8+GUsrmFRESPlCIDA4ePAiA84QIsQNuOW3Qvn17vWPHDooQIU1EbVtOU4gMlFLHAHxt9zoaQAcA39i9CBvg+3YePbTWHf2doBA5HKVUfk3/y4QyfN+hBWNEhBDboRARQmyHQuR8su1egE3wfYcQjBERQmyHFhEhxHYoRIQQ22HTK3EUSqmpAIqNh8la62V2rqepMN4nAKQAiAMwX2td4nM+ZD4HCpHDCacvrPletNZrjMfJSqmVWutpNi8toCilpmqtsy2PxwHYCvk7DsnPgcFqB1PDFzZLa13jFxYiVI78whotAtf6HNtnvt9QwPg7mqa1nu9z/ASAB7XWa0Lxc2CMyKEYX9gqXzxDcFyGIAHyhV5jOV8IIC14qwwcSqk4AP39nCpRSjnyPdXCVD/HiiF/tyH5OVCInE04fWGTAZT4OV5snAsJtNaFWut4P6eSAeQjRD8HCpFDCcMvrAveWJeVEkhsLGQxXOw8rfU2hOjnQCEKIcLhCxtuWGJGo+xeS1PCrFmIYPnCXlvnxc7F5edYqItqFoCbfY6F3OdAIWoGGDGb+XVeKEwzgs6+hPoXNh/+1+4CsC3IawkKSqks+JRjIEQ/BwpRM0BrnQcgr6HPD4cvrNa6RClVqJSK83mfccbnF1IYbvZK6386Sqk0rXVeKH4OjBE5nFq+sCUACo3smRUnf2GzYMkUKqX6oxEC3lwxLOR88+9UKRXnk+kMuc+BBY0OxvhyFhvBabPWZgAgVpYhUnFmNbXxhZ3m1IJG4ILwFkKsPUdXivvDiPXtq+F0vGkFhdrnQCFyKOH6hSWhCYWIEGI7jBERQmyHQkQIsR0KESHEdihEhBDboRCRZo85+EsppZVSJ4wCTn/XbTCu+aORVSQOgVkz4hiUUn8EkFbD1AGzwvy/zboq4hxoEREnsRKAb5UxgAvTKSlCDoVCRByD0ZpSAuAe63GjYhwUIedCISJO438AjDcfGG0tadaRuMR5UIiI0/gjqrpnC9m24nwYrCaOw9jRIg/AZwCyfcZhEAdCi4g4kf8BMA4yFpciFAJQiIgTccE7m5uEABQi4kT6A9hg9yJI4KAQEUdhZMmS4fCJhKQqFCLiNNIA1gyFGhQi4giMfrMsyLxmKKWy2E8WOjB9TwixHVpEhBDboRARQmyHQkQIsR0KESHEdihEhBDboRARQmyHQkQIsR0KESHEdihEhBDboRARQmzn/wHaNA8jz0F5EAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAElCAYAAABAu+8OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXxU5b3GnzcbECLZoBBIWIIiIsHIVhSQG5ZbKloBwaqACFdZFFBEtqRFhSYaRKsoKtiiFrW0olAuoq1otYKghkUjKHoh6Ay7JCFAIBDy3j9+55CTYbKSzMmZeb6fj59kzjkz553J+PDbX6W1BiGE2EmQ3QsghBAKESHEdihEhBDboRARQmyHQkQIsR0KESHEdihExHaUUh8opaLsXgexDwoRsRWl1AgAAwHE+Pi+UUqpCUqpt3x5X+IdChGxm98aPxN9dUOlVFeI+OX68r6kfChExDaUUrMAzDYe+sw101pv01qvArDXV/ckFUMhIraglEoEkK+1NsWAlkkAQyEidjFRa73M+D0fQHs7F0PsJcTuBZDAQyk1EMCXlkN7UYVgtVJqaVXvobWeWIOlEZugEBE7GOkhFFUKGlNc/BcKEfEpSqkJAGKUUpmWw4wPBTgUIuIzjKLF9lrrkR7HjwHI9P4sEghQiIgvmaC1nu3leH5VnlyN4sNcunHOgkJEfIJRQV1e3c5e45pESzr/IjwtqVrAp9XcpHyYvid1ilIq0bBk3gJwkZVixIxMK2mpkVHzxZrMYsquSqlMYx3EJhRnVhNC7IYWESHEdihEhBDboRARQmyHQkQIsR2m7w2aNm2q27Zta/cyCPE7zp49i8OHD+PIkSM/a62bebuGQmTQtm1bZGVl2b0MQvwKt9uN1NRUAMCKFSt+LO86umaEkDrBKkLp6ekVXkshIoTUOseOHSsjQgkJCRVeT9eMEFLrREdHY8CAAejfv3+lIgRQiAghtYjb7UZISAhatGiBsWPHVvl5dM0IIbWCGRNatGgRqts6RiEihFwy1sD0Aw88AKVUtZ5PISKEXBKe2bGqxIQ8oRARQi6JV199FUDNRQhgsJoQcok8+OCDOH78OFq1alXj16BFRAipNi6XC0899RSKiooQERFxSSIE0CIihFQTl8uFtLQ0AEBubi7i4uIu+TVpERFCqoxVhNLT02tFhAAKESGkiniKUE0D096ga0YIqRIlJSWIjIzErFmzalWEAAoRIaQS8vPzERkZiTZt2uDZZ59FUFDtO1J0zQgh5eJyuTBt2jS8/fbbAFAnIgRQiAgh5WCNCfXq1atO70UhIoRchFWEMjIyEB8fX6f3oxARQspQVFSEefPmAfCNCAEMVhNCPGjQoAHuuecetGnTxiciBFCICCEGbrcbBw8eRI8ePdC7d2+f3ptCRAi5MMojNDQU11xzDcLCwnx6f8aICAlwrPOEHnvsMZ+LEEAhIiSgsYqQrwLT3qAQERLAfPLJJwDsFSEAUNUdcu2vdO/eXXOnVxIoaK2hlILWGnl5eYiJianzeyqltmqtu3s7R4uIkADD7XZjxowZOHDgAJRSPhGhymDWjJAAwhoTOn/+vM2rKYUWESEBQm3stlFXUIgICQAOHDhQb0UIoGtGSEAQFRWFTp06YdSoUfVOhAAKESF+zYEDBxATE4Pw8HDMmTPH7uWUC10zQvwUt9uNOXPmYPHixXYvpVIoRIT4IdbA9B133GHzaiqHQkSIn1Gfs2PlQSEixI/QWuPJJ58E4BwRAhisJsSvUEph5syZ0Fo7RoQAWkSE+AVutxsrV66E1hrx8fGOEiGAFhEhjscaE/rVr36F6Ohom1dUfRwnREqpCcav7QFEAZittc73OJ9rPEzUWi/08RIJ8Rme84ScKEKAw4RIKTVBa73M8ngEgK0QUbogQlrrVcbjRKXUUq31RFsWTEgdUl+GmtUGjokRKaUSYQiOiSE4MYYgAcBEU4SM83sBDPTdKgnxHW63GyEhIY4XIcBhFhGACQBmexzLhYhRFICuXp6Tr5QaqLXeUOerI8QHFBUVoUGDBujVqxeuvfZaNGjQwO4lXTKOsYi01nu11t4c4EQAWcbPfC/nc41zhDgel8uFyZMnY8uWLQDgFyIEOEiIvGHEhDZorbcBiEFpkNpKPiSoTYijMbeBLi4udrwr5oljhciIGU3UWg+6hNeYoJTKUkplHT16tBZXR0jt4uu96H2NY4UIQCaAAR7HvA3fLdca0lov01p311p3b9asWa0ujpDaIi8vz69FCHBesBoAoJTKhEf9ECRO5E10YgBs88nCCKkDoqKiMGTIEPTu3dsvRQhwoBAZcaGlRmrePDZQa71BKbVXKRXlIVBRzJgRJ+JyuQAACQkJ+O1vf2vzauoWR7lmSqmBALJMEVJKRRnHTDIhKX7z+q4AKELEcZgxoUWLFiEQ9h50jEVkBKc/MH73PB0NSMzHCEAPhLhpiayqJk7DGph++OGHvX3f/Q7HCJFhBVX6F7G2gBDiNKwi5KR5QpeKo1wzQvydlStXAggsEQIcZBEREghMnToVubm5aNmypd1L8Sm0iAixGZfLhYyMDBQWFqJhw4YBJ0IALSJCbMUaE8rPz0d4eLjNK7IHWkSE2IRnYDoQLSETChEhNhCo2bHyoGtGiA0EBwejadOmmD59esCLEEAhIsSnHDt2DDExMWjZsiWeeuqpgChWrAp0zQjxES6XC9OnT8cbb7wBwGuHQMBCISLEB1hjQv369bN5NfUPChEhdQwD05VDISKkDjl79iweffRRABShimCwmpA6JCwsDJMnT0bz5s0pQhVAi4iQOsDlcmHTpk0AgO7du1OEKoEWESG1jMvlwuzZaSgqCsZf/tINrVo1xLBhQFKS3Surv1CICKlFXC4XHnooDWfOAIMHz0ePHg2xezcwZw5w663A+PHen5edDaxeDRw8CMTFIeCEi0JESC1hWkJnzgCPP56Ozp3FHevZExgyBBg9GujR42KBWb4cePtt4PbbgcGDUSXh8jcoRIRUA6vlEmREWEtKxIqJjf3CsIRKRcgkJga44w55rlWIsrNFhFaskGuAyoXLH2GwmpAqsny5WCrt2gGxscCmTfJfTEwJ2rUD1q8fjvz8xejRw3tgukMH4NChssdWrxZLKMZjRz6rcAUCFCJCqoDVcklOBrZvBz76CHj9dRe+/voBtGixDytWKLhcUdhQzr4xH34IfPcdMHkyMH++vObBg8CVV3q/3ptw+SsUIkKqgNVyMX8/dcqFhQvT0LDhcXz0UQhiYoBx44BXXwVyc8s+/7nngGeeAQYOlGtCQiQG9P77QGamiJIn338PtGjhk7dnOxQiQryQnS1Wi2m9fPNNqeVy8CAQG1vatnH33Rk4flx2YB0wAGjVSuI7K1YAn38OTJsmLl2HDkBxsYjPpk3A9OlAkyby2g89JK6fSW4u8Ne/SvYsEKAQEWIhOxsYOhS46Sbgq6+Avn0lJrRjB7DM2KgqOvogXnmldC/63Nz4C5bL998DKSli5eTkiAC98QZw443A00+LJfT888CvfiUi9/DDwFVXyXOffVauXbFChGzEiMAIVAMUIkIusHw5cN99wL59wJ/+JIL0xhvAuXNybu1aYONGYMSIGJw+nYxZszKwb1/8BYtp7lwRK7MGaNgwIDwc+M1vgJkzJRtWXAw89phYRdnZYiW1bCkidfnlIkY5OSJkgZK6B5i+JwRAaTC6Tx+gUydg0CA5bqbRMzOBG2/cj3vuiURaWgSSkx/Cr38NHD8uFlCnThI70hr48ksRIjOWlJMjtUE9e4pbN3gwEBEh59u1kzhQUhIwaxbw2mvAvHn2fhZ2QIuIEJSKRn5+2SyWmUZ/8003iormomPHZ5CTA/zwg9QRjRwpMaGgIOD114GPPwZWrSqbERs2DFi5UuI+cXEiSh06iOW1bBlw+LC4aQsXAoE6K40WESEotVSs1otJTIwb336bitBQ4Kqr7sa8eRLAvukmYMyYi1/rjjuAl16S15k7F+jXD/jlL8Wy6tNHzl1xBfDee0Dz5kCvXmIVPfwwcPSouIGB5JYBtIgIAVBqqVx9NbBggaTY588HPvrIjVdfTYVSwKlTGbjzTsmOVVT/s3s3sG6dpOrPnAGaNZPsWa9eEiNSCvjb34BGjYBFi+Q5f/wj8MADwJo1pRZVIEGLiBCI+zRqlFgmN9wA7NoFhIVpLFz4DJQCCgszMHZs/IUslilcVssJEAFZt06sotRUeb1Vq4Cbb5Zs2E03Scr+uuuA6GgRnhYtJAZlvra3VhB/h0JEiIWSEnGlUlKAd95R+P77h3H+fDGefjoet95aet2wYZKaHzKkbHvGm29KwHrUKHk8frz0i61eLRbQmjXy+OOPJUs2efLFgtOhA7B5c52/1XoFhYgEFOWN21i9WlLsyckSmP7ppw/RpcsYPPpoC+zYAezciTJClJQkj0ePFgumQwepIXrjDeCuu8qKS1KS/Pe734lV1KuXlAQEBXnvsg+kimoTxohIwGBtWh03Tn7OmSPHzZhPdLQb+/enIjj4Q0yZkoekpPJ7vsaPLy1cfO01+Xnzzd5jR1aXbcwYYOpUwOUCFi8uGxMKtIpqE1pExG/xHNmRnS2ukbdxG+3aAV9+6cbnn6cCkEH3sbGxACq2UExrx3rPqrhspkU1dSrQurX0ovXtKyIUSBXVJhQi4pd4DhvLzAQKCkSIrG6QWSf06aduvP9+KiIjgZSUdIwdm4Aff5QsV2EhcNttpVZLRZMUq+OymfGj554DPvsMiI8vG7QOJChExO/wNmzsF78AJk2SNHlkpMR8TDFp3x4oLj6K8PAwbNnyCNauTcDll0sc59QpScFnZ5daMzNnVjxJ0Rqg3rxZrKnyXLakJLGE4uMDs6LahDEi4nd4GzYWFydxnpYtpcjQjBMlJJzGk08CubnXIj7+Rfz4YwKUkvR9t26S3fr6a+klU0oELTlZ3LoxY0TsvNX9JCWJsLzwgvycNKm0utpKoMaEPKEQEb/DW7HhsGHSTpGdLdXNY8YALVu68e9/T0bz5p9g82Zg8eJQXHEF0L8/cP31Mqrj6afFYmrfXiyksWPLTk2s6iRFq8tmjgcJxC778qBrRvwOb8WGSUky3nXbNnGD3n1XKqbPnAEiI9sjPx8IDpaYTmgo8MQTcv1114kY/fwzcOwY8Omn0uhqpap1P95ctkCNCXlCISJ+R3nFhtHRUlR44IAbS5akIigIaNQoA+vWxaNJEwlKHzggja/LlskokE6dgPR0qYaOjZXs2/btZfvBqlP345llIwKFiDiGqu79VV7mav16oFGjAsTHpyIkBDhzJgNr18YjLKxUZLZuBYYPl9lDr70mMZzISBGnhATg//5PRCgjQ6ybVq0kxpOZ6fvPw59QWmu711Av6N69u87KyrJ7GaQcrOn4K68U12vlyqptWnjokASat2wRwbn77jVo0aI7Jk2KR3KyXNO4sQSWn3tO4jexscDp0+KmnTghMaO8PHmdmTOBPXvk9U+flhhPoHXL1wSl1FatdXev5yhEAoWo/mIWCVrT8YBYK6NHA//zP5KO/+YbyXYdOSLn27YFZs8WsZo71424uCJERLTHqlXion33nQSgT5+W5+zaJUHpCROAv/xFaohKSoAuXcRKMjNbq1fLGFmXC/jzn+lqVZWKhIiuGan3VLT3l5mOv+EGSbWfOCHxnCuvlMFj48YB//iHGwUFqTh8+DK8/PJi9OgRjPHjRXzatgUeeURGf1x3nQhOs2YiVAUFkmH7z3/K3jcpSUQxJ4ciVFswfU/qPeXN/snOBr79VgLK27ZJC8U338ixiAhpq+jWzY0jR1JRUgJ06zYXwcHBSEqSAsPERCAsTATnzjuBF1+UWNIrr4iV1LkzcNllrP3xBRQiUu/w3MonKEhiQp6sXi0W0VdfiXXStq1kvszannffdaNp01QEBwNHj2Zg/fr4C6IybJjEfFq1EotIgtcSK8rLE5Hq04e1P76CQkTqFd465L/+Gkb1c9lr//UvYMMGqXYOD5eKZ7ObvkMHICfn7QvZsRMn4suISmGhWDz/+IfUG/3znzIh8YsvZA51aGhptbRnh32g7bDhCxgjIvUGbz1iZod8v37ALbdIILlDB9m+eft2ielccYVk0IKDSy2WQ4cAre/Drl3H8MMPLdCwoaThMzNLCwqTk6X+5913pafsxhuBjh1ly6AHHpDthEaPljR9IPeB+YJKhUgp1VZrvc8HayEBTkVB6VmzJP2ekyMi8t13ki3bs0daMtasERG76y4X4uKWY9Gih6D1ZWjYsAWCgkRMJk8WUXvllbKvf//9IlwREZIpC/SxrXZQFdeMpVrEJ1Q0kL5DBwlGz5sngnLihKTb9+4VV0opICzMhcceS4PbvQfFxccRHCzuXEKCxHtefVWsn7ffLvvaJSWS5jcbVK2iU95QNFK7VEWIRiqlrqnzlZCAJy5OYj7WQLUZpzHbKMwY0qlT4popJQJVXOxCly5pCAoCtmxJx8mT8SguBiZOlJqgdu2AJUskw7ZkycX39RYMt96X1C1VEaIPAdyulEqu68XUBkqpCUqpEcZ/s+xeD6kezzwjdTzWUa6LF0u6/OqrxZrp0kWslLg42aJ5xgw34uLS0LAhUFycjqKiBCQkyOiNc+dke+fkZIkfHTsG/PRT2XtaN0C0wjS97/AaI1JKNdFaFwCA1nqQcexWpVSi1vodXy6wOiilJgDI1VqvMh4nKqWWaq0n2rw0UgnZ2WINde8OPP64xHSGD5ftdx55BJgyRaqne/cWFys9XTYjvP9+oHXrBggNbYGsrKkID0+A1sDJkyJcZjuIOcCsWTOJM1kprzctUMe22sFFLR5KqQEAtNb6o4suVupaAAO01ot8tL5qYZSQd/M4tkdr3b6y57LFw16GDpVK6BkzJPu1fDmQlSV1PQ0aSHq+USMgKkoC1FFRwJkzRzFgQAyKioKxfbvGzp0K7drJ+YYNpYG1vfGXz80Vy2b7dukbe//9i5tor75axO7QIXHHymuqJTWjyi0eSql7AGRprXd4u1hrvV0ptVcp9QSADNNqqg8opaIAdPVyKl8pNVBrvcHXayJVIztbBOJPfwIGDRIRatBAqp/XrZOU+s6dcuyrr6QI8exZFzp1SkNOTj8MH/4/aNtW4dw5qay+7DIpSOzRQ4SocWOJKe3aBTRtKkLkOdN69245duutTNXbgadrlgMgEYBXIQIArfVxAHOUUg8rpT7QWn9lPW9juj8RQL6X47nGOWIzFe0pFhkpbRZnz0pcZ/x4caEyMsRlGzJErmvWTEToqqvSEBoKtG7935gxQ6qie/WSCumWLYEffpCes1OnJC7UsKHcMy9PLJ/ly73XK5l1Q7SEfEuZYLXW+kMA25VSwyt6klKqP4CeALYppf5piJKZWbMrHhMDER1P8gFEeXuCEdjOUkplHT16tE4XF+h4VkyHhIj1kZICLFwoFk9JiVQ7t20LvPSSiEn79uKiNWwoO2l07OhC585pKCkBPvkkHRs3JuDMGRGUzz4TlyoiQlL2Z86IsPz1r9IY27Ch3Pejj8qvV6rK2FdS+1yUNdNa51gD0kqpx42f/ZVSLyqljgHYAGAggJcBHAcwCSJg5wE4JlOltV6mte6ute7erFkzu5fjt1grpseMEfdp0yZg+nSpAyosFFfqN7+R0Rq33Sbd9Hv2AHffLWL0xRfAd98Vo0GD+QgKAjZvTkdcXALcbhGQU6fk9cLCxM2aMgV46CGxpl57Tdyy6OjSNo6K6pVYN+R7qpK+n2gRn98CeAvAIK11jNZ6ktb6Nq315QCijfM5FbxWXRPj5ZhXa4j4DmvFtFWUJk8WEWrZUiwll0viNxs3ipiEhoooHD8umxBGR4egQYOp+OabdJw/n4BBg0RcOnSQrX5mzZKAduPG8vpnz8r5666T2NK4ccCAAVJ7xLqh+kVVhCgKUktkFZ8PPS/SWh830uaranuRVSQL3kUnBsA2H6+FWPjmGxk6P3myWCS9e5e6RUVFElwODwd27JDm040bRYTOnZNrYmPdaNbsI2zaBAwc2AXHjycAEJeraVPA7ZYA9ujREuBu1UrS/tnZIm7WRtXvv5cYEOuG6hdVaXpdprWeVI3XXFrTxVwKWut8I6MXpbW2Bq2jmDGzj+XLRWD69xeLZO5cabNo0UKEQWtJ26emylygTz6RWE5BgcRzwsLcCAtLRU5OEBo37oWHHgpHXBzQpo3srPHFF5LuLygQsdmxQwLSH38sMSfrBEVTaDIzgS+/ZN1QfaIqQlStXjOttZ2uWSaACQAWAoBSqivEpSQ2YLphaWlSgFhUJMKTkiKbEgYHm+0ZEtvZtQv49a+lJSMlBdi7140+fVJRXAzs378ABQXhWLZMxOftt8Wl++ADsbQ2bRJxWrFCrKH//EdcsB07xP3zFJqkJG7tU5/wu5nVRnX1Xoiblqi1XliV57GgsWZUtLPG/PnSTmGe27VL3KYNG4CrrpI4zbBh0km/ebMMQLvxRrn+22/duP76VDRoAFxzTTreey8BkZEy/gMA7rlHxOb3v5eg9saNwI8/yp5lKSll50uzQLF+EFAzq7XWy+xeQ6DgrSjQbKXo0QP429+khufuu0v3jTdFa+tWiQtdf73MGtq5U8azvv++WEr9+u1ASAhw4EA6Dh1KwJ490mUPiIt14IC0guzcKSL07bcSe+rcuazgUHicgd9ZRDWFFlH1qGhnjX79JFv188+SJh8+/OKtf0aOFMuoa1e5/vBhEaCVK89j8+ZgfPYZ8MMPJ9C27WV48kmpKRo+XILRpos1fryI4eLFkuKfPbtq2wwRewgoi4j4hvKGmO3fL4Jy/fVSI7R/v1QsW6uWW7WS4z17SkB6yhR5fMcdLlx/fQYmT56Ow4c7YNasy1BQICLVurWMhk1MLI3lZGeL6DRtCjz6qBxjhbQzoRCRGnHwoLhjnqxeLQ2sR46Im9SkSWl2qmdPEYzTp2UuUOvWEtO55RagWTMX4uPTcOYM8NVXjVBYKEHooCBx6+66Syywr7+WFo3CQqnIPnxYSgKsgmOtkKYQOQMOzyc1orxhYgcPyk8zOHzwoKTmc3IklvPpp/L7W29J8NrtBoKCXCguToNSwBdfpCMvLwGxsUC3bkDz5vKcl1+WILR1kP1PP4k75s0FY4W0s6AQkRpR3jCxqCixRMyA8a23SuPqTz+JK3X2rKTr/+u/pOP++ecPoVmzNEREAMeOpSM4OAGbN4vFtHixZNOio4F33hFLKClJuuNfeEHiRefPe18fK6SdBV0zUiNMkbnlFtnOR2up28nJkd9btZLrxo+XzvoHH5RWjVGj5PekJCluXLs2FidO9ESbNrcgISEBXbpISr51a9mHvm9fID9fWkB27pR7mgwbJu7akCEXB8zNwkXiDChEpMp41gxZCQoSAQIkWJ2cXDpl8Z13pF1j9GiZOQQAbrcbx45dhpEjI5GdPQXbt4v4TJgg/WLPPSfd9PHx0kP2zjsXu1qcrOg/UIhIlZg/H3j9dRmv0b69pOZXrgR+9zuJ3aSnS5/YqVOSCRs+XILSqakiUp07S5uHFDm6kJubhuDgNmjZcgEAsaYAiTuNGSOWUHy8uGErVsg5b67W+PGskPYHWEdkwDqi8pk/H3j+eeCxxySAvHs3sGAB8MtfSgvFkSMyI/qyy2TC4vr18vjOO4GnnpLrQkOljePee13YuDENxcXAkSPp2Ls3AcnJYg19/bVYQy++CEydKoLSqpWMBSkqkrgQBca5VFRHRCEyoBB5Jztb3J/p02VvMNM1+/xzsX5OnJAYzoQJIlJjxki9z6hRIkY33CCxoe3bgT59XIiLky1/0tPTkZ2dgKFD5T433CA1Qu+9J9enpADXXiv301pcLcB7KwlxBhUJEbNmpEJWrxZ3bP/+i/ekP3RIRKhtW7FYzGFjf/+7xIni4qTFo1EjIDRUo6TkBeTlAXFx6ZgyJQG33y7TFLt0kYD2P/8p20eb/Wbr14v4jRghwme9t7nHPfEPGCMiFXLwoKTP162TEatmdurRR6Uu6ORJicuYdUVnzwJr14rrtnKlbPkTGQkMH66wb9/DaNz4NN59Nx5ut7R59O8vnfMvvCDZrtGjpVvftHa8tZKwetr/oEVEKiQuTqweT5KSRIBKSkSA2reXeNBdd8nI1/BwiR2dPetGXt4yHDhwHoMHx6Jz53jcdptkxf78Z3m+GYT2NjO6vFYSzpf2LyhEpEKGDZOG0wEDxAJZsULcpBdflOmHERHioj36qLhnJ04AsbFi1bhcbrRvn4rWrTciLy8Xa9YA//536R733iYielZEm9d6g9XT/gOFiFRIUpIEjTdskHjNp58C06aJu6S1iI9S0qoRHi4p+6efBvLyZKhZdDTw/PMZuP/+ZigslHEdR49KfGf06IvrfTwrorkvfWBAISKVsmCBxHkKCsQ62rdP4kYxMVLrEx8vqfm9eyVN36iRGykpqWjaFFi4MAPx8fEYOlSe07mzpPzfekue43ZLHAjwbiFxX/rAgMFqUilJSTLc7JVXJO7Ttq20dfz0k3S+Dx1aWuvz4IPAzJnHUVjYCCEhv8cnn8SjpEREY9w4EasffhAB2bVLBpxNmiQTGw8cuNhCYvV0YMA6IgPWEVXO/fdLjCcuTqyj6dOlaBGQ2NGzz55CQUFjtGoFHDlyHh07BiMrS1y7BVJAjTFjZHzIE0+Utozs2ycu3xNPlO0ls2Jey7GvzoWD0YhXKpo37Y2SErFcZs+Wka7WzveTJ91o2jQNHTqMQufO/41XXw3Gyy/LudGj5eebb0pcyRwbaw6xB0TIPJtarVivJf4HY0QBiucW0FUpEjQbXc2h92bsxu124733UhEcrNG8eSds3CjiM3q0bB1kDkR74w0Z3eFNUJgBC2xoEQUg1t1Wq1MkOGyYbN+zYoVcK7EbN8LDU3H2LLBnTwYKCuIxbpw0o5oWl7nR4c03l5+KZwYssKFFFIDUtEgwKUmsp4MHzX3HTiIsLA0nTwJZWRlo0iQeL7xQOjHRHGI2dKjsNTZpEjNgxDu0iAKM7Gzgf/9X+sdyci6OC3XoIOM0ysMcu/HSS8C//hWBRo3uwJVXdsa99zNzjRgAAAvFSURBVMZj+fLSgWgm1iFlzICR8qAQBRDmPmStW8v/9GZcyLr1TlVcpOhoN6ZMOYGrrroKQOkE/ePHKxcZzg8i3qAQBQjWuJDZST9t2sXb/FQ2YtXtdiM1NRXh4eFYsmQJgoODL5yrqsgwA0Y8oRAFCNa4UEyMDCvr1UtctJISqQkKC6vYRTJFCADS0tLKiJAJRYbUBAarAwRr8+jy5dK4evfd0jUvwWYRpvJ2R7WKUHp6OhISEnyzcBIQ0CIKEMzmUZdLKph795ZtfaZOld6vnTuBLVvEhfNm0axbtw4ARYjUDRSiAGHYMKloPnVK5kqPHy/CNGOGWEtvvin7hpW3O2qvXvfC5RqOjIxfcFQrqXXomgUYRUUSRDZbDEtKSnfQ8KxudrvdSEtLwwsv5CEtLRgdO/6Co1pJnUCLKEBYvRqYObPU/TpyRATpj38U12z1ahEYM3VvxoSKi4EdO05hxYpojmoldQaFKEA4eFC63pOTxZr5wx9KK6sLC6WrfssWSbdbA9Nt2mSgU6f4CquwKUTkUqFrFiCYwWprdbM59vWVV2QQ/ogRUqxoilBGRgby8+M5qpXUORSiAMHaLT9+vFg+OTnA0qXAhx8Cjz8uxxs3bow2bdogI0MmK3JUK/EFFKIAwdMSKiws3Zts7lygT5/DKC4uRnR0NBYsWID4+HgAHNVKfANjRA6nOsPNymvBiIpyYebMNPTq1Qv33XdfmeewUZX4Ao6KNXDiqFizifX226VqevdusV6sTayV4XK5kJaWBgAX3DFvcFQruVQ4KtbheLN6gJoNN7NSVREC2ENG6hYKUT3HavUMHixWz5w5sm1PZcPNKhKO8+fPIyMjA0DlIkRIXUMhqsdUNNI1OVm28PFGZcPNACA4OBgPPPAAIiIiKELEdpg1q8dUNNK1Rw8RKW9UlFZ3uVxYv349AKBjx44UIVIvoBDVYyra9334cBndUZ20uhkTWrlyJU6ePFn7CyakhtA1q8eYxYQ9e158rqRENi6salrdGphOT09HRESED94BIVWDQlSPGTZMAtNDhpR1z6wD6YHKR7N6ihDnCZH6BoWoHlPVYsLK0uq7d++GUgp/+MMfKEKkXsKCRoP6XNBY02LC4uJihITIvzWnTp1C48aN63ilhJQPCxodTk2KCV0uFxYsWIBp06ahc+fOFCFSr6EQ+SHWmFBkZKTNqyGkcpi+9zMYmCZOhELkRxw9epQiRByJo1wzpdQE49f2AKIAzNZa53ucN0v8ErXWC328RFuJiYlB3759MXjwYIoQcRSOESKl1ASt9TLL4xEAtkJE6YIIaa1XGY8TlVJLtdYTbVmwD3G73WjUqBFiY2Nx77332r0cQqqNI1wzpVQiDMExMQQnxhAkAJhoipBxfi+Agb5bpT2Yg+6ffPJJsBSDOBXHWEQAJgCY7XEsFyJGUQC6enlOvlJqoNZ6Q10vrjqTEmsL624b999/P5S5QRkhDsMRFpHWeq/WOtrLqUQAWcbPfC/nc41zXlFKTVBKZSmlso4ePVrj9S1fLq0Y7drBZxsQci964k84ySIqgxET2qC13qaUGojSILWVfEhQ2ytGzGkZIJXVNVlHRTOD6nIDwpdffhkARYj4B44UIiNmNFFr3c3utVQ0M6guNyCcMWMGCgoKOE+I+AU+FyLDevGM9ZTHRCPo7EkmgAEex2K8XFeuNVRbmDuoeqMqkxKrg8vlwpo1azBp0iQ0adIETZo0qb0XJ8RGfC5ERuC4xsFjpVQmPOqHIHEib6ITA2BbTe9VFSqaGVSbGxBaK6ZHjhyJFtzZkPgRjghWmxhxoaVWK8nIiuUD2Gtkz6xE1XXGzBcbEHq2bVCEiL/hmBiR4dJlmSJkiI51pEAmJMW/0DjfFZdgeVWVut6AkL1jJBBwhBAZwekPjN89T0cDkgEz0vEDIW5aoq+qqsvbQbU2gtRFRUWIiIjA3LlzKULEb+FgNIP6NhitoKDgQjD6/PnzCA4OtnlFhFwaFQ1Gc1SMKFBwuVyYMmUK1q5dCwAUIeL3UIjqGdaYUNeu3rpWCPE/KET1iOrsRU+IP0EhqicUFhbi97//PQCKEAk8HJE1CwTCw8MxduxYXHHFFRQhEnBQiGzG7Xbj2LFjuOaaa5CSkmL3cgixBbpmNmKO8nj++edx7tw5u5dDiG1QiGzCOk/okUceQWhoqM0rIsQ+KEQ2YBUhBqYJoRDZwgcffACAIkSICVs8DHzR4qG1hlIKJSUlyM3NRdOmTev0foTUJ9jiUQ9wu92YPXs2jhw5gqCgIIoQIRaYvvcB1phQUVGRzashpP5Bi6iO4W4bhFQOhagO2b9/P0WIkCpA16wOiYyMxOWXX45x48ZRhAipAApRHXDo0CHExMQgIiIC8+bNs3s5hNR76JrVMm63G7NmzcKSJUvsXgohjoFCVItYA9MjRoyweTWEOAcKUS3B7BghNYdCVAuUlJTgiSeeAEARIqQmMFhdCwQFBWHGjBkICQmhCBFSA2gRXQJutxurV68GALRr144iREgNoUVUQ6wxof79+yMyMtLmFRHiXGgR1QDPeUIUIUIuDQpRNeFQM0JqHwpRNcnJyUFQUBBFiJBahDGiKnLu3DmEhoaib9++6NGjBxo2bGj3kgjxG2gRVQGXy4XJkydj69atAEARIqSWoUVUCdZtoJs3b27zagjxT2gRVQD3oifEN1CIyuHYsWMUIUJ8BF2zcoiOjsagQYOQkpJCESKkjqEQeeByuRASEoK4uDiMGTPG7uUQEhDQNbNgxoQWLVoE7vdGiO+gEBmcPXv2QkzowQcfhFLK5hUREjhQiAz2798PgPOECLEDbjltEBsbq3fs2EERIqSOqGjLaQqRgVLqKIAf7V5HDWgK4Ge7F2EDfN/Oo43Wupm3ExQih6OUyirvXxl/hu/bv2CMiBBiOxQiQojtUIiczzK7F2ATfN9+BGNEhBDboUVECLEdChEhxHbY9EochVJqAoBc42Gi1nqhneupK4z3CQDtAUQBmK21zvc47zefA4XI4QTSF9Z8L1rrVcbjRKXUUq31RJuXVqsopSZorZdZHo8AsBXyN/bLz4HBagdTzhc2U2td7hcWIlSO/MIaLQLdPI7tMd+vP2D8jSZqrWd7HM8DcK/WepU/fg6METkU4wtb5otnCE6MIUiAfKFXWc7vBTDQd6usPZRSUQC6ejmVr5Ry5HuqgAlejuVC/rZ++TlQiJxNIH1hEwHkezmea5zzC7TWe7XW0V5OJQLIgp9+DhQihxKAX9gYlMa6rORDYmN+i+Fib9Bab4Offg4UIj8iEL6wgYYlZjTI7rXUJcya+QmWL2y3Si92LjFejvm7qGYCGOBxzO8+BwpRPcCI2cyu9EJhohF09sTfv7BZ8L72GADbfLwWn6CUyoRHOQb89HOgENUDtNYbAGyo6fMD4Qurtc5XSu1VSkV5vM8o4/PzKww3e6n1Hx2l1ECt9QZ//BwYI3I4FXxh8wHsNbJnVpz8hc2EJVOolOqKSxDw+ophIWeZf1OlVJRHptPvPgcWNDoY48uZawSnzVqb7oBYWYZIRZnV1MYXdqJTCxqBC8K7F2LtObpS3BtGrG9POaejTSvI3z4HCpFDCdQvLPFPKESEENthjIgQYjsUIkKI7VCICCG2QyEihNgOhYjUe8zBX0oprZTKMwo4vV33gXHNW0ZWkTgEZs2IY1BKvQVgYDlTB8wK87+ZdVXEOdAiIk5iKQDPKmMAF6ZTUoQcCoWIOAajNSUfwEjrcaNiHBQh50IhIk7j7wBuMx8YbS0DrSNxifOgEBGn8RbKumdz2bbifBisJo7D2NFiA4AvASzzGIdBHAgtIuJE/g5gBGQsLkXID6AQEScSg9LZ3MQPoBARJ9IVwAd2L4LUHhQi4iiMLFkiHD6RkJSFQkScxkCANUP+BoWIOAKj3ywTMq8ZSqlM9pP5D0zfE0JshxYRIcR2KESEENuhEBFCbIdCRAixHQoRIcR2KESEENuhEBFCbIdCRAixHQoRIcR2KESEENv5f6C4VEKW6Yu9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Init Loss = 50.116447\n",
      "Final Loss = 5.8048\n",
      "Fair dummies test (regression score), p-value: 0.47652347652347654\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAElCAYAAADX8PdgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3db5BV9Z3n8c/XRtjB6uZ22w4SsTCXhNJuLEg3UxYOaJk0E5pKWE1BspAaZ7HWxp4qyIMQWFdipSawBsI+GH0w0GZhZWqhpIm6CQU4tGNFULGG7sFKaGuN3EjUiY4t3cBCRaTz3QfnXLhcb/+5f87907xfVbfoe+758zu3m9uf/v2+53fM3QUAAIDCuq7UDQAAABiLCFkAAAARIGQBAABEgJAFAAAQAUIWAABABAhZAAAAESBkoajMrM3MfIjHyRK2od/MDplZvBhtKEdmtiR8H0Z6dJe6rQBQCcaVugG4Zi2QlJAUS1k2MJoNw1/yp919QZ5tqHX3gXCfTZIeldRtZkvdvSvPfUeigOeeSZekr6Ut65bUIWlbBMcDgDGNkIVSOZYMODkoyC/81OO7e4+kpWa2TVKnpNpCHCMCkYWd8P3oSV1mZpJ0Mnx/8mJmMUkt7r43330BQCVguBAVx9073L0jon2vlCQz2xTF/vMV5bkXwRwFARYArgmELODz9khaUupGAAAqGyELpVJnZrHUR/IFM2sxs+6UgvSrepXMrNPMOlOebzOzTeF2h8xsbZ5t65Z0uQA+/XjhMjezJSnPt4WPtWZ2Mmz32vDcDiUL+8Par3y2ydSWTSlF6ZtSlq8NlyX305bn+3KVbI5rZockHUp579zMWgrZHgAoN4QslMpJSf2pj5Qr++KS1kmaLmmppLUjBIQ6BT1PnQp+kec7nJaQpCyvNKyTlGzjgrANmyS9FP7brKCwvzPPba4SBq4lCgrWl0paEoaalnAfX1PwPq6UdDqL8xlWtscNC/WXSpK7W/goy4sLAKBQKHxHqdQOVfieVnOUMLO9Cn5BDxee4pKaC1GgHe5L7p7IcrsBd98sSWb2hKS1Cgr8u8Jl2/T5wvVctlH4WlxB0Ll83ma2TsFVkk8oCGiJ8H3O9lyGVKrjAkClIWShLIVDcd9REHiaFEwvMJyeAgUsKeiByeXKx2PJL9x9ILwyL3VOqUw9Sblsk5QcbusOt0vVFW7bb2Zdkg4lw1wBlOq4AFBRGC5E2QnngnpU0rMKerBGc8l/IXtM2hQUv2crUzAbaYgul22u2j5l+C31MeDuyeHWAUmb0mu58lSq4wJAxaAnC2UlLPJukjQ9OVyXobckyuNvUjBB6rq0l1IL82MqD12SYmbWNFQvXjgn1d6wZ7BQYadUxwWAikLIQrlJ9kitCwNPk4L6n0iLpMNwt1JBL9aCDPVi8bAW6bSkp6Nsy2i5e8LMOiR1mtlKBe9dk4Jh1mcl/UX474CuzLBfyuMmQ/OS8DVR/A5gLGO4EGUlDDebFYSdbgW/pPcqtxqpYaVMJeAKruirU9CDlv6L/4nwtZMKemWeUBAYCt6mbIWTpyaL40/qyjBrj4J6tpfC5XMUXt1XquOGvV49Ct7DspzsFQAKydy9uAcMhlq+raC3YFQf+uHl+8k6lTiFtAAAoNwVdbgwHJJJDrmMag6iZMBK3u/MzOJmti15+xMAAIByVPSeLOly2Hra3ZtHsW53+npmdjK8ggkAAKAslXVNVji02JThpQFuyQEAAMpZWYcsBUOKQ80jlM0tTwAAAIqq3KdwqFPmiRkHlDJvUaqwhqtNkv7sz/6s+bbbbsv96OfOSX/608jrXXedVF2d+3Ewppw7Jw2O4semih8blKm33nqrz91vKnU7gEpX7iEra+F97zokac6cOX7s2LERthhmX/feo8FH/mbE9aq2PiP71Ss5Hwdjy9z5l3TfspEnbX95d51ePzzm/gtiDDCzU6VuAzAWlPtwoRT0ZqUrlxm3AQAAMir3kHVMmQNVnYJJDQEAAMpSWYescPbvRIZ7xcW4HQcAAChnpQpZmYYAkxONdqaFqk0KC9nDdZoU8X3sAAAA8lXsGd/jCm72u0BSU3gD4JNhsboUTMvQoiCEJW8g22FmbeG8WDEFt9VhtncAAFDWihqy3D2h4Oa/Ge89GA4B1mZY3pFhdQAAgLJV1jVZAAAAlYpJegAAJdPd3X1bVVVV23XXXdfq7p8byQDKlZn1/+lPfzowODjY0dzc/G6mdQhZAICS6O7uvu36669/bvLkybFYLHZu/PjxfWZW6mYBI3J3Xbx48fqBgYH/9NFHHy3s7u7+VqagxXAhAKAkqqqq2iZPnhybPHny6QkTJnxGwEKlMDNNmDDhs8mTJ5+ePHlyrKqqqi3TeoQsAEBJXHfdda2xWOxcqdsB5CMWi5277rrrWjO9RsgCAJSEu9eOHz/+s1K3A8jH+PHjPxuqnpCQBQAoGYYIUemG+xkmZAEAAESAkAUAABABQhYAABHasmVLvZk1D/Xo6+urKnUbM2ltbY23trbG812nnJhZ8wsvvFA9mnULcW7MkwUAKDuLvjk4/Q8f+vhStyNpys12cf8vq07ms4/nn3/+7RkzZnx6+vTpq3731tfXD+bXuqs1NjbeMWnSpEuvvfbabwu5X2SPkAUAKDt/+NDH3/fdT8rmysOX//eNeQe+efPmXQgD1cUCNGlIK1as+Dj1eV9fX9Uvf/nLmhUrVvRHeVx8HiELAIAxZM2aNX2pz48cOTLxoYceiq9YsaK7VG26VlGTBQAAEAFCFgAARfDv//7vVX19fVc90tdZvnz5tJqamtm33nrrzB07dtSmF2pnKsY2s+YdO3bUZlrn7rvv/vIDDzwwI7le6v5eeOGF6sbGxjvMrLmmpmZ2e3v7LentOXPmTFV7e/stt95668yh1knX3t5+S01NzezRrL98+fJpy5cvn7Z+/frJyWOsX79+cl9fX9Xdd9/9ZTNrvvXWW2ceOXJkYup2qeu3trbGM72XfX19Va2trXEza25sbLxjy5Yt9fm2N1uELAAAiqCxsfHOm266aXbqIzUwtba2xnfv3l3/+OOPv7979+7E008/nTEUZOO111777fbt2xOS5O7d7t59//33n5Okd955Z8LGjRvfP3HixK937tx5cuvWrTenB5HXX3+9Zvr06Z+++OKLb69evfoPmdZJ1draGt+3b1/t/v373965c+fJffv21Q63fn9/f9Xu3bvrJenFF198+7vf/e7HGzdunHrvvffOWLt27YeHDx9+q6amZnDZsmWX36fly5dPe+aZZ2566qmnTh09erS3trZ2MB6P35ketBYvXhx/9dVXa7Zv357o7Ow82dPTc0O+7c0WNVkAABTBxx9/fHyoKwl7e3vHHzx4sHb79u2JZIH6z372s1ONjY13RtWe1NqthoaGiwsXLux/7rnnalOXz50792zy+YYNGz5KJBL/4e///u9vTq/7Sj2Hw4cPvzVv3rwLktTf3//Bli1bMq6fVF1dPbhhw4aPJOnHP/7xh1u3br151qxZF5JhcMWKFR//4Ac/mCYF9WW7d++uP3HixK8bGhouStKuXbtONTY2Tly9evXUXbt2nUq25fXXX695/vnn307uZ9euXaeSgS6f9maDkAUAQIm98cYbN0hS6hWAyRARpR07dtTu2bOn9ve///2E3t7eiXPnzj073PoLFiw4mxpUUu3fv79GkubPn39HNm2YOXPm+eTXyRDa1NR0edmNN954OZgePXp0YnV19WD6e3PPPfecfeWVV2qSz5PvZzJgFbK92SBkAQBwDWpsbLxDktasWfPhXXfddf773//+1DNnzgw7Meonn3wy7OvV1dWDZ8+ePZ5NOyZNmvS53r3UYJVqYGBgyOOfPXs260ldc2lvNqjJAgCgxO66667zUlCMPtK6qUEo19nijxw5MrG3t3diZ2fnyRUrVvSPttfspZdeqmloaLiQ6bVFixadPXfuXFV6kXohLVy48FymY+zbt692/vz5l3utMr2f6e9VMdpLTxYAAEVw5MiRiZlqsm6//fZPGxoaLjY0NFx48MEHp+/cufPkjBkzPv3+978/NdN+3nvvvQm9vb3j//zP/3zwr//6r6eNdNwvf/nLn0rB0GBtbe0lKZgYVZI2bNgwZf369X944403bjh48GBt+nDh66+/XrNly5b6RYsWnd21a1dtsm4s03EaGhouLlu2rG/ZsmXxp5566tSMGTM+feONN27Ys2dP7YEDBzJuk6158+ZdWLhwYX/yGPX19YMbN268+cyZM+OefPLJ91PbMnfu3LOrVq2aJulUfX394Nq1a6+6crAY7SVkAQBQBMmpFNIli91/9atfvb148eL4Aw88MKOhoeHCmjVrPjx48GBt6rqPPfbYh4sWLZrR2Nh459y5c89u3rz5g9/85jcTk+Epk3nz5l1oaGi48NBDD8UbGhoubNy48f36+vrBRx555MOtW7fevG/fvtpvfOMb/QsXLrxqRvja2trBhQsX9vf09Nzwd3/3d1MnTZp06ac//emp4WaO37Vr16n169dPXrVq1bT3339/QvI8sn2vhnPgwIFEe3v7LatWrZp25syZcTNnzjx/9OjR3vQA+4tf/CKxePHi+IMPPjh95syZ5x9++OG+9957b0Ix22vuXqh9lZ05c+b4sWPHct7e771Hg4/8zYjrVW19RvarV3I+DsaWufMv6b5lp0dc7+XddXr9MH/noPyYWbe7z4n6OG+++ea7s2bNyngV11f+4tIdZXZbnev/9V/GvVXs45pZc+oVcihPb775Zv2sWbNuS1/OJzwAoOxMudkuFuJ+gYUy5WaL/Eo/jD2ELABA2dn/y6qTpW4DkC+uLgQAAIgAPVkAAJQpd+8udRuQO3qyAAAAIkDIAgAAiAAhCwAAIAIlqckyszZJyYmE4u6+eZTbJMUkdbj7QBTtAwAAyFfRQ1YyYLn73vB53My2ufvKYbZZq7RQZWbbJA25DQAAQCmVYrhwZTJgSZK7JyS1jLDNX2TotRows1jBWwcAAFAARQ1ZYShqyvDSgJkNF7TiGV6PMVwIAADKVbF7suKSMgWj0+FrQ1kn6ZCZbZIuDzluK3zzAAAACqPYIatOVwreUw0oKGbPyN27JDVLajMzl5Rw955omggAAJC/ipjx3czikuZI+qKkRxX0aq10944M67ZJapOkKVOm6Pjx47kfuKVFfv0NI7evpUXK5zgYU76+wPWF2B9Hsd4HOn7citAiAEAplCJk1WVYNlIB+7qUqw/Xmdmzkl4ys0TYy3VZGLw6JGnOnDk+e/bsnBvq31utwS/dMuJ6VV1dsh8+nvNxMLa0r7qk++ozddhe7eVDdfrR4xXxdw5QdP4fF0/XRx+OL3U7Lpt880X7P7/I6abVW7Zsqf/BD34wLXXZ1KlTP50/f/65J5988v36+vrBwjQys9bW1rgkHThwIBHlcUrNzJqff/75t++///5zI61brPek2J/wx5Q5UNVJyjj8Fxa8H0pd5u49ZrZU0gJJXZm2AwBUsI8+HD/Y/p8/K3Uzkqr+4X/lHfiS9yHs6+urOnLkyMTHHntsajwevzORSPw66qCVqrGx8Y5JkyZdeu21135brGNeq4pakxVeDZjIMPVCLL1HahSOSfqkMC0DAKA46uvrB++///5zJ06ceGvSpEmXVq9ePbWYx1+xYsXH3/rWt/qLecxrVSnGKjYpqJnaLElm1qSU3qiw/mqTpIfdfcDdu8ysU9LetP20KRwWBACgEs2cOfPCu+++W9Rh0TVr1vQV83jXsqKHLHfvMLO2cBgwpuC2Oqkzt8cVTE5apyvTPTwcTt+QHA+PSdrLPFkAgErV19dX9eqrr9asXr36D6VuC6JRkhtEu3uHu3e5+970+xaGy2vDmeCTywbcfV24XYe7b059HQCActfX11fV19dX1dvbO37Hjh21X/nKV+74y7/8y7MbNmz4KLnOCy+8UN3Y2HiHmTXX1NTMbm9vv+rqq/Xr10+uqamZbWbNt95668wtW7bUJ19rb2+/paamZnam7VK1trbGk4XfkrR8+fJpy5cvn7Z+/frJt95668yamprZ69evn5y+3Wj3n2qofff19VXdfffdX06ex5EjRyamn2dy/dbW1nhfX19V+nvZ2toaN7PmxsbGO1Lfh3zbXEglCVkAAFxrbrrpptk33XTT7MbGxjsfeuih+Pz588/94z/+46nUdd55550JGzdufP/EiRO/3rlz58mtW7fenAwQL7zwQvXGjRun7t+//+0TJ078+qmnnjp14403DkpBcNq3b1/t/v373965c+fJffv21Q4VPNL19/dX7d69u16SXnzxxbe/+93vfrxx48apvb29l4cxc93/UPu+9957Z6xdu/bDw4cPv1VTUzO4bNmyq0LfM888c9NTTz116ujRo721tbWD8Xj8ztSgtXjx4virr75as3379kRnZ+fJnp6ez823lM97UihcPw4AQBEkry6UpN7e3vEbNmyYEo/H79y5c+fJ5LQDqfVSDQ0NFxcuXNj/3HPP1a5Zs6avv79/XHV19eDtt9/+aX19/WBDQ8PF5L4OHjxYe/jw4bfmzZt3QZL6+/s/2LJly82jrb+qrq4eTPao/fjHP/5w69atN7/99tsTGhoaLua7/0z7njVr1oXkOa9YseLj5BQXR44cmbh79+76EydO/Dp5frt27TrV2Ng4cfXq1VN37dp1qre3d/zrr79ekzpdw65du04lw1yh3pNCIGQBAFBkDQ0NF3ft2nVq+fLlevDBB6efPXv28ozWO3bsqN2zZ0/t73//+wm9vb0T586de1aSvvnNb559/PHHL910002z586de/arX/3q2Q0bNny0f//+GkmaP3/+Hbm2Z+bMmeeTXyenk+jv7x8nSfnuP9O+m5qaLi9L9sZJ0tGjRydWV1dfDpBJ99xzz9lXXnmlRpLeeOONGyRpuPmwCvGeFALDhQAAlMjf/u3ffnzu3LmqZE1SWF9087e//e3+zs7OkwsXLrw81UJ9ff3ge++995vt27cnJk2aNLhx48apydqq6urqQXfvTn+Mth2TJk0adp6ufPafad+pwSrVwMBAVablknT27NkhX8sk3/ekEAhZAACUyMGDB6slad68eReOHDkysbe3d2JnZ+fJFStW9Kf35iStWLGi/8CBA4nt27cnDh48WLto0aKzqUGt0KLef6qFCxeey3Ssffv21c6fP/+cJN11113npaBGLfl6emF8Mds8HEIWAABF1tvbO769vf2WjRs3Tn3sscfel6Tbb7/9U0nasGHDlOQViAcPHqxNbrNjx47a9vb2W8IwNv7QoUM1U6dO/bShoeHismXL+pYtWxZ/4YUXqpPbpl5BmI+o959q3rx5FxYuXNifPNaRI0cmtra2xs+cOTPuySeffD/Znrlz555dtWrVtOQ6ixcvvqotxWzzcKjJAgCgCMysOfl1dXX14MyZM8+nFm/X19cPPvLIIx9u3br15n379tV+4xvf6E8dLrzrrrvO79mzp3bRokUzzp07V9XQ0HBh9+7dCSko/F6/fv3kVatWTXv//fcnNDQ0XFizZs2HhWp71PtPdeDAgUR7e/stq1atmnbmzJlxM2fOPH/06NHe1FsP/eIXv0gsXrw4/uCDD06fOXPm+Ycffrjvvffem1CqNg/F3L2YxyuqOXPm+LFjx3Le3u+9R4OP/M2I61VtfUb2q1dyPg7GlrnzL+m+ZaO4QfTuOr1+mL9zUH7MrNvd50R9nDfffPPdWbNmZbzSayzdIBpj35tvvlk/a9as29KX8wkPACg7BBqMBdRkAQAARICQBQAAEAFCFgAAQAQIWQAAABEgZAEAAESAkAUAKBUfy9MI4doQ/gxn/EEmZAEASsLM3j1//nxJb3sC5Ov8+fMTzex3mV4jZAEASuKzzz7b/O67747r6+uLXbx4cRy9WqgU7q6LFy+O6+vri7377rvjPvvss82Z1mMyUgBASTQ3Nx/o7u5+54MPPnj03/7t32a5e5344x+V4U9m1u/ur1y6dOmJ5ubm32ZaiZAFACiZ8JfTQ6VuBxAF/mIAAACIACELAAAgAoQsAACACBCyAAAAIkDIAgAAiEBWIcvMvhpVQwAAAMaSbHuyuszsEzP7BzO7L5IWAQAAjAHZhqxaSY9Kmi7ppTBwPWFmswvfNAAAgMqVVchy9zPu3uHufyWpTtJ/VRC4Os3st2HgmhZFQwEAACpJzoXv7j4g6RNJpyVZ+GiW9DszO2hmtxWigQAAAJUo65BlZrPDmqxPJHUqGEJc6e5fSunhOiXpnwrbVAAAgMqR1b0LzewdSV+U9K8Khgr3uPuZ1HXcfcDMNkl6uGCtBAAAqDDZ3iC6S9JWdz8+wnqfSFow1Itm1qZgmFGS4u6+eTQHN7O1kgaS27r73tFsBwAAUGzZhqxuSd+RdFXIMrN/kvSsu/9PKSiQl/RSph0kA1YyIJlZ3My2ufvK4Q5sZockLQ1rwWRm/WbWlXwOAABQTrKtyVor6V8yLO+Q9Mgo97EytQfK3ROSWobbIOzB6kwLVM0ELAAAUK6yDVnTJfVkWN4tqWmkjc0sNsR6A2Y2XNB6VMFQ5WVhOAMAAChL2YashKSvZFjeJOl3o9g+rqCmKt3p8LXPCYNZLPx6iZm1mNnacDkAAEBZyrYmq0PSz8xswN1fli7fz7BD0hOj2L5OVwreUw0oDFIZzEm+nlLHdUzB9BGfK64Pa77aJGnKlCk6fnykGv1htLTIr79hxNWspUXK5zgYU76+wPWF2B9Hsd4HOn7citAiAEApZBWy3H2zmd2o4JY6nvLS0+6+pbBNu0pMQS9ash0DZlZnZk3uftXwpbt3KAh9mjNnjs+enfsdf/x7qzX4pVtGXK+qq0v2w8dzPg7GlvZVl3Rffaa/Ja728qE6/ejxbP/OAQBUiqw/4d19nZn9dwW9SC6px91HM1SYVJdh2XBDf4nwuOnDjKcVFMxnqhEDAAAoqZz+jA6naMhljqpjyhyo6jREWHL3hNmQQypcXQgAAMpS1iErvCdhXBnCkrs/N9y24TBfwsxiaT1TMXfvGnJDqcfM4mlXFMYVhDYAAICyk+1tdX4g6ScKbgadziVVjWI3mxQUpm8O99mklOkZzCwervNwShBbFz5WpmyTSK/HAgAAKBfZ9mQ9Kumnkp5Iv2fhaLl7h5m1hfNixRTcVid1tve4glqrOoXDge7eZWaxcFJSSbrR3Ye8bQ8AAECpZRuyTiu4d2FOASspvAJwqNe6JNVmWM59CgEAQMXIdjLSbZKWRNEQAACAsSTbnqzpktrM7FEFRedXXd3n7t8pVMMAAAAqWS5TOKQO2zFdNQAAQAbZzvj+SFQNAQAAGEuyrcmSmf0XMztmZpdSlj1rZg8UtmkAAACVK6uQZWY/UTBf1RNp23ZK+m8FbBcAAEBFy7Ymq03SEnf/57QbRHdLaipcswAAACpb1sOFkj4J/00teo9LyuYm0QAAAGNatiHraUl7zWy2gtvoyMy+Imlr+AAAAICyv7pwXXhvwR5JMrNPFNwap8Pdt0TQPgAAgIqU9TxZ7r40DFpfCxd1uTtDhQAAAClymYxU7p6QlChwWwAAAMaMrEKWmX1rmJcH3P2f82wPAADAmJBtT9beIZa7gp6tL+fXHAAAgLEhq6sL3f269IekOkk/l7Q0khYCAABUoFzmybqKuw9IeljSpvybAwAAMDbkHbJCdZJaCrQvAACAipdt4fuzGRbHFMz4PlCQFgEAAIwB2fZkWYbHGUn/KnqyAAAALst2xvdvR9UQAACAsaRQNVkAAABIkW1N1lezWZ/JSQEAwLUq28lIuxRMPGppy4daVpVjuwAAACpatsOFSyW9JGl62oSkP5e0JG2iUgIWAAC4ZmUbsn4iaa27/y5t+TpJjxamSQAAAJUv25B1o6TmDMvjkprybw4AAMDYkG3I6pTUYWYPJBeExfB7FNRrAQAAQNnPk7XSzOok/dzMPFxsknrEDaIBAAAuy/bqQrn7UjP7oq4MD/ZkqNECAAC4pmUdsszsv0h6RNJsdx8XLntW0rPu/lyB2wcAAFCRsqrJMrOfKLiS8AldPS9Wp7K4utDM2sxsSfhYm00bwu07s90GAACgmLLtyWpTMB/WP5tdNfdot0Z5daGZtUk67e57w+dxM9vm7itHuX2TpCXZNRsAAKC4sh4ulPRJ+G9qyopLGm1d1kp3vzwNhLsnzKwli+PHs1gXAACgJLKdwuFpSXvNbLaC2+bIzL4iaWv4GJaZxZS5x2tgNEHLzJYke8AAAADKWbZTOKwzs7iCKRtkZp9IiknqcPcto9hFXNJAhuWnNUIPVXjcRDbtBQAAKJVcp3CIS/pauKgriykc6hQEqnQDCsLacJroxQIAAJUiq5BlZmsU9FolVMRepXAocVQzyoeF9W2SNGXKFB0/fjz3A7e0yK+/YeRjtrRI+RwHY8rXF7i+EPvjKNb7QMeP24jrAQAqU7Y9WX+loBbrf+RxzLoMy4bsxQrruOTumYYZP8fdOyR1SNKcOXN89uzZubQx2Nf3VmvwS7eMuF5VV5fsh4/nfByMLe2rLum++kwdtld7+VCdfvR4LteeAAAqQbaf8Jsk7TGzLnd/M4fjHVPmQFWnsM4rgzbp8tQNl4Xzaw2EoQoAAKCsZBuyFkjql9RjZj1KGzJ09+8Mt7G7D5hZwsxiaT1TMXfPOBzo7pvTl5nZpkzLAQAAykW2UzjEFPQ4/VzBvFiW9hiNTQp7p6TLPVRdKc/jZtaZHCYEAACoRNlO4fBIvgd0947wtjotCkJbPG2297ikFgVDiFfVYYXbLA2/3iapc6geMAAAgFIaMmSZ2W8lLXD3d9OW35a+LFvD1VGFoal2mNe6JI3qFjwAAAClMtxw4XSlFamb2SRJJyNtEQAAwBiQbU2WNPraKwAAgGtWLiHLC94KAACAMWakkDWqCUABAABwtZGuLvyimaVOXV0nSWZWrQzDhu5+toBtAwAAqFgjhaxM0yOYhu7hqsqvOQAAAGPDcCGLaRIAAAByNGTIcveni9kQAACAsSSXqwsBAAAwAkIWAABABAhZAAAAESBkAQAARICQBQAAEAFCFmJU79gAAArFSURBVAAAQAQIWQAAABEgZAEAAESAkAUAABABQhYAAEAECFkAAAARIGQBAABEgJAFAAAQAUIWAABABAhZAAAAESBkAQAARICQBQAAEAFCFgAAQAQIWQAAABEgZAEAAESAkAUAABABQhYAAEAECFkAAAARGFeKg5pZm6TT4dO4u28e5TaSNF1STNI6dx+IqIkAAAB5KXrISgYsd98bPo+b2TZ3XzncNu7ekfJ8iaRuBYELAACg7JRiuHBlMmBJkrsnJLUMtbKZxZUWpsLt68KwBQAAUHaKGrLMLCapKcNLA2Y2ZNCS1JZh2WlJdQVpGAAAQIEVe7gwLilTHdXp8LXPCXu6aofY17HCNQ0AAKBwih2y6nSl4D3VgIJi9lEJ67q63L1niNfaJGnKlCk6fvx4jk2V1NIiv/6GkdvT0iLlcxyMKV9f4PpC7I+jWO8DHT9uRWgRAKAUSnJ1YT7CGq2V7t6c6fWwQL5DkubMmeOzZ8/O+Vj+vdUa/NItI65X1dUl++HjOR8HY0v7qku6rz7T3xJXe/lQnX70eMX9FwQAjFIpCt8z1VGNuhdL0iZJXytQWwAAACJR7JB1TJkDVZ2kzw39pTOzTWJ+LAAAUAGKGrLCcJQIrzJMFXP3ruG2DWuttoWF8Mllw12RCAAAUDKlGC7cpJQpGcysSVJXyvO4mXWmBrEwTB1LBiwzixGwAABAOSt61a27d5hZWxiSYgpuq5M623tcweSkdQrmz4pLOiRJZp+7EivT1A4AAAAlV5JLm1JvkZPhtS6lhKew94rr3AEAQEUpxXAhAADAmEfIAgAAiAAhCwAAIAKELAAAgAgQsgAAACJAyAIAAIgAIQsAACAChCwAAIAIELIAAAAiQMgCAACIACELAAAgAoQsAACACBCyAAAAIkDIAgAAiAAhCwAAIAKELAAAgAgQsgAAACJAyAIAAIgAIQsAACAChCwAAIAIELIAAAAiQMgCAACIACELAAAgAoQsAACACBCyAAAAIkDIAgAAiAAhCwAAIAKELAAAgAgQsgAAACJAyAIAAIjAuFIc1MzaJJ0On8bdfXMU2wAAAJRK0UNWMiy5+97wedzMtrn7ykJuAwAAUEql6Mla6e7NySfunjCzlgi2AQAAKJmi1mSZWUxSU4aXBoYKTblsAwAAUGrFLnyPSxrIsPx0+FqhtgEAACipYoesOl0pXk81IClWwG0AAABKqiRXF0YpLJJvC5/+PzP7v3nt8JXD6UvqJfVlOHBehylzmc95bMvrnI8eGd16ZfRjw/f42jDac54WdUOAa0EpQlZdhmUj9UiNeht375DUkW2jRsvMjrn7nKj2X44457HvWjtfiXMGEL1iDxceU+ZwVCepp4DbAAAAlFRRQ5a7D0hKhFcMpoq5e1ehtgEAACi1UtxWZ5Ou1EzJzJokdaU8j5tZZ1qoGnabIotsKLKMcc5j37V2vhLnDCBi5u7FP2hQnJ5QMAx41S1ywrmvOiU1u3tiNNsAAACUm5KELAAAgLGuFMOFAAAAY96YmycrH8kbUYdPRzUkmcs25SSPc5ak6QqGb9eFFyhUhHy/Z2bW6e5LC9+y6OR6zma2VsHEv6clKXmT9kqQ58+2FPxsd1TCz3ZYw/ptSQtG+7NZ6Z9dQCUgZIWSHzjJXyJhAf42d19ZyG3KSa7nHM5Flny+RFK3gsBV9vL9noUXXSyJso2Flus5m9khSUuTIcPM+s2sq0JCRy4/22uVFqrMbJuksv7/HP5MxpXFrcYq/bMLqBTUZIXMrNvdm9OWnXT3IcNDLtuUk2zbb2ZxSSvdfV3a8n5JD1dCL0e+37MwVHa6e/nM1T6CHH+210oaSAvU8dSLUcpZjuf8uR5KM9sk6YkKCZZNkp5OP+8h1q3ozy6gUlCTpctd7U0ZXhoIr3YsyDblJI/2t2VYdlqZZ+UvK/l+z8xsSSUEyVR5nPOjSpsmpYICVq7nHM/weqwSAlY2Kv2zC6gkhKxAXEHdSbrhut9z2aacZN1+d0+4e+0Q+zpWwLZFJefvWdiLVxEhI03W5xz+Eo6FXy8xsxYzW5thQuBylev3eZ2kQ2HvVXJIbVvhm1dylf7ZBVQMQlagTlcKQFMNaOj7KuayTTkpSPvDX0Rd7l4JtzjK55ybKuQc0+VyznOSr7v73vDOCh0K5q+rBDl9n8PzbJbUZmYuKVGh3/ORVPpnF1AxCFnIWUqN1oJStyVK4RDKtXYLp5hSeu7CIbO6sO5nTAp/nudI+qKkzQp6tTINjwPAqBCyrshUUzTSX3W5bFNO8m3/JklfK1BbiiWrc04OkVV4XU623+eElPGcT0uqlJqdXH6217l7h7sPhBd3NEvaNEbrlCr9swuoCEzhEDimzB8wdZKGGi7IZZtyklf7w7qVipofS7mdc5t0+cqtyzJdfVemsj5nd0+YDXnxZCV8v7M+5zBIHUpd5u49ZrZU0gKNrZ7MSv/sAioGIUvBX+xmljCz9CuJYmGdRkG2KSf5tD9ZEJx2b8mWcj/vHL/Pn5ug0cw2VcrEjXl8n3syTNlQERc4FPj/5jFlvhKvYlX6ZxdQSRguvGKTUqYnCHsuulKex82sM+0Kq2G3qQBZn3P4F/+x5C9fM4tV2HBKLt/nSpfLOa8LH6nbVFIheFbnHIaL72TYT5uCov9KkHEalTH62QVUBCYjTRH20CQUdKVfdZuJMEh0SmpO68EZcptKkM05h4XBJ4fYVW2lDB3m8n1OeW2prvzi7ayUv/xz/NleoiuX9N+YPgltucv2nMMQ8qiu/IzHJO0t9/nBwv+XSxQMa7YoKNo/mRzKHqufXUAlIGQBAABEgOFCAACACBCyAAAAIkDIAgAAiAAhCwAAIAKELAAAgAgQsoAyE85pVCk3YwYADIGQBQAAEAFCFgAAQAQIWQAAABEgZAEAAESAkAVkYGabzGybma01s5Nm1m9m20ax3TYzO5m2LGZmnryRtpm1mFl3uKzfzDaNsM/PFcKH2y7J0Ob+0ewTABA9QhaQWVzBjaCl4KbQHZLazGztCNt1SoqHN+1N+rakgZSbScclrZM0Pdz32vBmvTkLQ9gSSV8L97kk330CAPJDyAKG1uXum929x93XSdoraeVwG4RBaiBtvZWS9qSs0+HuXe6eCNffqyAY5SQMdEskLQ3b2qUgxA3bVgBAtAhZwOgdUtALNZIOBaFHZhaT1CTpqqFGM1sSDgN2J9fNQ0v4b3II0hX0qDXluV8AQB4IWUAOwjotT3mk1kxtUzBk2KQrQ4U9Kdt2S3pU0rMKerD2FqBJA+5u6Y8C7BcAkKNxpW4AUEGWSuqRJHefPtRK7p4ws4SCHqYFCnq2JElh8GqSNN3dE+Gy0Rw7lrKPWNprXZJiZtaUGuYAAKVFTxYwtBYzazOzeHi1XoukJ0a57TZJ3wm3SR0qTIT/rgv3u0SjGy6Mh+vHJD2d+kIY1jokdYZXLsaTw5GjbCsAIAKELGBoPZKaJSXrpla6+2iH9joU9Fglkj1WkuTuA5I2K7hysVtBT9deBcXyQ3lCUp2kkwpqrZ5QENYub+PuKxWEuW3hesnhSABAiZi7l7oNQNkJe4Fi7r6g1G0BAFQmerIAAAAiQMgCAACIACELAAAgAtRkAQAARICeLAAAgAgQsgAAACJAyAIAAIgAIQsAACAChCwAAIAIELIAAAAi8P8Bh5wGrIdy36cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_groups_pointwise(x_axis_0_test,\n",
    "                      x_axis_1_test,\n",
    "                      A_test,\n",
    "                      Y_test,\n",
    "                      Yhat_out_test,\n",
    "                      y_u=None,\n",
    "                      y_l=None,\n",
    "                      point_color='b.',\n",
    "                      shade=None,\n",
    "                      filename_0=\"synth_equi_0.png\",\n",
    "                      filename_1=\"synth_equi_1.png\",\n",
    "                      save_figures=save_figures)\n",
    "\n",
    "\n",
    "p_val_equi = utility_functions.fair_dummies_test_regression(Yhat_out_cal,\n",
    "                                                            A_cal,\n",
    "                                                            Y_cal,\n",
    "                                                            Yhat_out_test,\n",
    "                                                            A_test,\n",
    "                                                            Y_test,\n",
    "                                                            num_reps = 1,\n",
    "                                                            num_p_val_rep=1000,\n",
    "                                                            reg_func_name=\"Net\",\n",
    "                                                            lr = 0.01,\n",
    "                                                            return_vec=True)\n",
    "\n",
    "\n",
    "plt.clf()\n",
    "x_equi, bins_equi, p_equi=plt.hist(x=p_val_equi, bins=30, color='royalblue', alpha=0.7, range=(0, 1), density=True, ec='blue', align='left', label=u'Equitable model')\n",
    "x_base, bins_base, p_base=plt.hist(x=p_val_base, bins=30, color='salmon', alpha=0.7, range=(0, 1), density=True, ec='red', align='left', label=u'Baseline model')\n",
    "\n",
    "plt.grid(axis='y', alpha=0.75)\n",
    "plt.xlabel('p-value')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('Fair Dummies Test')\n",
    "for item in p_equi:\n",
    "    item.set_height(item.get_height()/sum(x_equi))\n",
    "for item in p_base:\n",
    "    item.set_height(item.get_height()/sum(x_base))\n",
    "\n",
    "plt.ylim(0, 1)\n",
    "plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left') \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}