{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Before we begin, we will change a few settings to make the notebook look a bit prettier" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%html\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Understanding VERTIGO\n", "In this notebook we provide a step-by-step explanation of how VERTIGO works and how we implemented them. \n", "\n", "For ease of use, you can find an object version of VERTIGO in [`vertigo.py`](./vertigo.py). For a (simulated) distributed scenario where we actually use it, see the notebook [`demo_vertigo_local`](https://nbviewer.jupyter.org/github/IKNL/vertigo/blob/master/scripts/demo_vertigo_local.ipynb). For a real-life distributed implementation of VERTIGO using our [privacy preserving DL infrastructure](https://github.com/IKNL/ppDLI), see our repository, `d_vertigo`.\n", "\n", "All right. Let's get started." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Preliminaries\n", "Display plots." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import packages" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import sys\n", "import os\n", "import pathlib\n", "import pandas as pd\n", "import numpy as np\n", "import numpy.matlib as npm\n", "import seaborn as sns\n", "import warnings\n", "import matplotlib as mpl\n", "from matplotlib import pyplot as plt\n", "import statsmodels as sm\n", "from sklearn.model_selection import train_test_split\n", "\n", "import auxiliaries" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set (default) plotting parameters." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "mpl.rcParams['font.sans-serif'] = \"Calibri\"\n", "mpl.rcParams['font.family'] = \"sans-serif\"\n", "sns.set(font_scale=1.75)\n", "sns.set_style('ticks')\n", "plt.rc('axes.spines', top=False, right=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Setup paths." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Setup working path\n", "os.chdir(os.getcwd())\n", "\n", "PATH_DATA = pathlib.Path(r'../data/')\n", "PATH_RESULTS = pathlib.Path(r'../results/')\n", "\n", "if not PATH_DATA.exists():\n", " raise ValueError(\"Data directory not found.\")\n", " \n", "if not PATH_RESULTS.exists():\n", " path_results.mkdir()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Setup execution parameters." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Test/train partitioning.\n", "TEST_PROP = 0.2\n", "SEED = 31416 # For reproducibility.\n", "\n", "\n", "# VERTIGO execution parameters.\n", "LAMBDA_ = np.array([1000])\n", "TOLERANCE = 1e-8\n", "ITERATIONS = 500\n", "\n", "# We know alpha needs to be bounded between 0 and 1 (see Fig. 4, 3.ii).\n", "# Thus we will initialize it at 0.5.\n", "ALPHA_INIT_VAL = 0.5 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The data\n", "We will be using the [Breast Cancer Wisconsin Diagnostic Data Set](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+%28Diagnostic%29). The data consists of real-valued features computed from a digitalized image of breast biposy. The features describe the characteristics of the cell nuclei present in the image: \n", "* `radius` - mean of distances from center to points on the perimeter\n", "* `texture` - standard deviation of gray-scale values\n", "* `perimeter`\n", "* `area`\n", "* `smoothness` - local variation in radii lengths\n", "* `compactness` - perimeter^2 / area - 1.0\n", "* `concavity` - severity of concave portions of the contour\n", "* `concave points` - number of concave portions of the contour\n", "* `symmetry`\n", "* `fractal dimension` - \"coastline approximation\" - 1\n", "\n", "The mean, standard error, and \"worst\" or largest (mean of the three largest values) of these features were computed for each image, resulting in 30 features. For the sake of simplicity, we will only consider the mean values of the features, resulting in 10 features. The output variable is `diagnosis`, which states if result of the biopsy was a benign (`B`) or a malignant (`M`) tumor." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
diagnosisradius_meantexture_meanperimeter_meanarea_meansmoothness_meancompactness_meanconcavity_meanconcave points_meansymmetry_meanfractal_dimension_mean
id
842302M17.9910.38122.801001.00.118400.277600.30010.147100.24190.07871
842517M20.5717.77132.901326.00.084740.078640.08690.070170.18120.05667
84300903M19.6921.25130.001203.00.109600.159900.19740.127900.20690.05999
84348301M11.4220.3877.58386.10.142500.283900.24140.105200.25970.09744
84358402M20.2914.34135.101297.00.100300.132800.19800.104300.18090.05883
\n", "
" ], "text/plain": [ " diagnosis radius_mean ... symmetry_mean fractal_dimension_mean\n", "id ... \n", "842302 M 17.99 ... 0.2419 0.07871\n", "842517 M 20.57 ... 0.1812 0.05667\n", "84300903 M 19.69 ... 0.2069 0.05999\n", "84348301 M 11.42 ... 0.2597 0.09744\n", "84358402 M 20.29 ... 0.1809 0.05883\n", "\n", "[5 rows x 11 columns]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv(PATH_DATA/'data.csv', index_col=0, usecols=range(0, 12))\n", "df.head(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice how in this case, all variables are continuous. If the data had categorical values, they would need to be encoded, preferably using [one-hot encoding](https://en.wikipedia.org/wiki/One-hot)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data pre-processing\n", "Extract the features only" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "column_names = list(df.columns)[1:]\n", "X_centralized = df[column_names]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to avoid having features that are dominant just because of their range, we will normalize them." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "X_centralized = (X_centralized - np.min(X_centralized))/(np.max(X_centralized) - np.min(X_centralized))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Re-structuring of the data\n", "For the sake of this example, we will simulate two different parties. Each of these parties will have half of the features." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
radius_meantexture_meanperimeter_meanarea_meansmoothness_mean
id
8423020.5210370.0226580.5459890.3637330.593753
8425170.6431440.2725740.6157830.5015910.289880
843009030.6014960.3902600.5957430.4494170.514309
843483010.2100900.3608390.2335010.1029060.811321
843584020.6298930.1565780.6309860.4892900.430351
\n", "
" ], "text/plain": [ " radius_mean texture_mean perimeter_mean area_mean smoothness_mean\n", "id \n", "842302 0.521037 0.022658 0.545989 0.363733 0.593753\n", "842517 0.643144 0.272574 0.615783 0.501591 0.289880\n", "84300903 0.601496 0.390260 0.595743 0.449417 0.514309\n", "84348301 0.210090 0.360839 0.233501 0.102906 0.811321\n", "84358402 0.629893 0.156578 0.630986 0.489290 0.430351" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X1 = X_centralized[column_names[0:round(len(column_names)/2)]]\n", "X1.head(5)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
compactness_meanconcavity_meanconcave points_meansymmetry_meanfractal_dimension_mean
id
8423020.7920370.7031400.7311130.6863640.605518
8425170.1817680.2036080.3487570.3797980.141323
843009030.4310170.4625120.6356860.5095960.211247
843483010.8113610.5656040.5228630.7762631.000000
843584020.3478930.4639180.5183900.3782830.186816
\n", "
" ], "text/plain": [ " compactness_mean ... fractal_dimension_mean\n", "id ... \n", "842302 0.792037 ... 0.605518\n", "842517 0.181768 ... 0.141323\n", "84300903 0.431017 ... 0.211247\n", "84348301 0.811361 ... 1.000000\n", "84358402 0.347893 ... 0.186816\n", "\n", "[5 rows x 5 columns]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X2 = X_centralized[column_names[round(len(column_names)/2):]]\n", "X2.head(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Output\n", "VERTIGO uses an output of -1 or 1." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "y = df[['diagnosis']]\n", "y = y['diagnosis'].map(lambda x: 1 if x=='M' else -1)\n", "y = pd.DataFrame(y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Class balancing\n", "This is particularly important for any logistic regression, including VERTIGO. If a class is underrepresented, the model will apparently perform well, but it will likely just predict one class all the time. Thus we need to balance it. We will do so using downsampling." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "X1_balanced, y_balanced = auxiliaries.downsample(X1, y, random_state=0)\n", "X2_balanced, _ = auxiliaries.downsample(X2, y, random_state=0)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAFlCAYAAAB7rFmPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXyM5/7/8dfIHqRExJJobY1WI6LWOMjpsW8tdag0GEqLFj1FD7pr+22bOqqHKA2Nk6aEtF9Li6KUo5SgKoKS1lKNip1YsgyZ3x9+c38zMtEJk1ryfj4eeZjc1zXXfG6z5DPXdpusVqsVERERESk1ytzqAERERETkz6UEUERERKSUUQIoIiIiUsooARQREREpZZQAioiIiJQySgBFREREShklgHJb2blzJ6+99hqdOnWicePGNGjQgMjISEaOHMnXX3/Njexa9Le//Y169epRr149Tp8+fcOxZWRkGO1069bthtu5Ef379zceOy0t7U99bLnzTJs2zXi9XPsTGhpKREQE/fr148svv7zpxyqNr03b+TZq1Mjp+1y5coXHH3+cevXqkZKS8of1U1JSjMcZOnTozYQrLlbw/fXJJ58Yx8ePH28cX7FixZ8a04QJE6hXrx5z5sxx+j7uJRiPiNPy8/OZNGkS8fHxhcoyMzPJzMxk1apVhIeHM3XqVKpUqXILohS581ksFk6fPs3p06fZunUrBw8e5Pnnn7/VYd31/vOf/7B7925CQkJo3rz5rQ5H7jIDBgxg4cKFTJ06lU6dOlGtWrU/vI8SQLktvP/++8Y3F5PJRMeOHWndujU+Pj78/PPPJCcnc+rUKXbs2MGAAQP44osvKF++vFNtv/nmm+Tk5ABQrly5G46xUqVKTJ8+HYCyZcvecDsif6YuXbrQtWtXAKxWK3l5eaSkpLBgwQIAPv74Y3r37k316tVvZZh3tQsXLjBjxgzgao+p3J369+9Pu3btAGjQoMGf+tgPPvggTZo0Ydu2bXz44YfExMT84X2UAMott337diP5K1OmDLGxsbRt29auzoABAzCbzaSnp3Po0CGmTJnCa6+95lT7rVq1ckmcPj4+xptb5E5Ru3btQq/brl278vvvv/Pdd99x5coVUlNTlQCWoHnz5nH+/Hk8PDzo1KnTrQ5HSshDDz3EQw89dMsev2vXrmzbto1ly5YxduxYKleufN36mgMot9z8+fON23369CmU/AH4+/vzP//zP8bvX3zxBRcuXLCbJzNt2jRiY2Np1qwZ4eHhvPTSS8D15wBu2bKFAQMGEB4eTrNmzfjnP//JyZMnGThwoHEfm6LmAF4bw8GDBxk5cqQRx4ABA/jhhx8KndNvv/3GuHHjiIyMJDQ0lEaNGtGhQwfeeuutm5qrKOKMe+65x7hdoUIFuzJXvTYXLFhA7969efjhh3nooYdo0aIFZrO50PyoG30PAaxcuRKz2UzTpk1p2LAhHTt25M033+To0aOF6h48eJCxY8fSsmVLQkND+dvf/sbEiRM5duyYw7aXLl3K3//+dxo2bEhERARvvPEGWVlZTp+/je0zrnnz5vj5+dmV5eXlMX36dDp27EhoaCjt27fnk08+IT8//7pt/vjjj4wePdp4jlq1asWoUaP48ccf7erZPv9atWplN4d6zpw5xv/566+/bnef3r17G/NFs7Kybvj5Wb16NYMGDaJ58+bUr1+fRo0a8fjjjzNr1iwuX75cqH5KSgqDBw826jdt2pSePXsya9YsLBaLXT1bPHFxcezYsQOz2Ux4eDgtWrTgtdde48KFC1y6dImYmBhat25NWFgYjz/+OCtXrrR7zIULFxptLVy4kO+++46oqCjCw8Np2bIlr732mtOv+aLmABacJ5ufn8+nn35Kly5daNCgAY888giTJ08mNze3UHt79+5l2LBhNG7cmMaNGzN8+HB+/fVXXn75ZaO9jIwMo37btm0xmUxYLBaSk5P/MF71AMot99133xm3H3300SLrhYWFUbNmTQ4dOkRubi4//PAD3t7eRvnChQv5/fffjd9DQkKu+7hff/01Y8aM4cqVKwBkZ2ezZMkStm/ffsNDxXv27GHOnDlcvHjROJaSksLAgQNZvnw5NWrUAODYsWNER0fb/eGxWCz8+uuv/Prrr3z//fcsWbIET0/PG4pDxBGr1UpOTg7r1q3jm2++AaBmzZo0bdrUqOOq12ZsbCzTpk2zO3bmzBk2b97M5s2biYmJoUePHoXu5+x7CODtt98mMTHR7v6HDh3i0KFDrFq1irlz53LfffcBsG3bNp555hm7do8cOcK8efNYtWoViYmJ1K5d2yizfaG0ycnJISkpqchEtCi7du3iyJEjALRs2dKu7MqVKzzzzDNs2rTJOHb48GHef//9635+zZo1iw8++MAuSTxx4gQrV67km2++4cUXX+Spp54CriaAiYmJnDhxgr179/Lggw8CsHnzZuO+BRelnD59ml27dgGOE1Znn5+5c+fy5ptv2t330qVL7N69m927d7Nv3z7+9a9/GWXff/89zzzzjF2il5WVxZ49e9izZw8HDx7knXfeKfR/8d133zF16lTjftnZ2SxYsICMjAyysrLsFibt3r2bUaNGMWvWLNq0aVOoreXLl7NhwwYjUba19f3335OcnIy/v3+h+xTXuHHj7BZf/f7778TFxXH27Fneeust4/i2bdsYMmQI2dnZxrFvv/2W7du3F/naqFKlCnXq1OGXX35h/fr1PPfcc9eNRT2Ackvl5ubafbv6o6StTp06xu3ffvvNruz333+nUaNGTJo0CbPZTJcuXYps58KFC7z66qtG8hcZGcl7773H8OHDOX78OD/99NONnA7ffvstgYGBTJw4kddff92YK5iXl2fMuQL45JNPjD+w3bt3Z8qUKbz55pvGH6sDBw6wY8eOG4pBpKDY2Fijt+CBBx4gPDycf/zjH+Tm5hISEsKMGTNwd/+/vgBXvDZPnz7Nxx9/DEDlypV57bXXmDJlCtHR0Uad//3f/3V4X2ffQ6tXrzaSv7JlyzJq1Cj+9a9/GSMIJ06cMEYN8vLyGDNmDBcvXsTT05OhQ4fywQcfMGzYMDw9PTl58iQTJkww2v7555+N+b4A3bp1IyYmhqeeeooDBw5c99yvtW3bNuP2Aw88YFe2YMECI/nz8PDgmWeeISYmho4dO5Kenu6wve+++47JkyeTn5+Pu7s7/fv3JyYmhgEDBuDu7k5+fj4xMTGsX78ewG7433bs8uXLbN261Th+8OBBjh8/btSxJZaOprw4+/zY5jzec889/POf/+TDDz9kwoQJBAQEAPDVV1/ZvY4mTZqExWLB3d2d5557jilTpjB27Fh8fHwA+PLLL+2SQ5stW7YQEhLCe++9ZyS9ABs3bmTXrl0MGDCAmJgYQkNDjbIvvviiyP/b2rVrM3HiRF566SVjCPW3335j0qRJDu9TXF999RX9+vVj8uTJdqNdixYt4sKFC8DVRZETJkwwkr/w8HD+53/+hzFjxmC1WtmyZUuR7dteY2lpaZw/f/66sagHUG6ps2fP2v3u6+t73foFy699cbu5ufHvf/+bKlWqXLcnEa7+8bDdPzw8nI8//hiTyQRc7REZN26c0+dQkLe3N4mJicYHR25uLu+99x5w9Zu9zZNPPklYWBhHjhzhmWeeMR778uXLxrdm2weySEk5d+4cW7Zssev5csVr08fHh7i4OPbt20doaChNmjQBri5IWb58OWfOnCmyDWffQwWnjrz11lvGQpeuXbvSv39/ypUrxwMPPMDly5f59ttvyczMBODpp59m1KhRRl03NzemT5/Ojh072LNnD/Xr12fZsmVGL1C3bt2YPHmy8Vj+/v52PVd/5JdffjFu16xZ065s6dKlxu1XX32VJ554AoAePXrw9NNPGwlbQR9//LER2xtvvEHv3r2N+4SEhPDKK68A8NFHH9GmTRuaNm1KhQoVOHv2LOvXr2fo0KGkpaXZ9eDB1USqW7duxmOaTCaHCaCzz4+t/YCAANq1a2d8gWjatCm//vortWvXNl53VquViRMnsnfvXjw9Pe16htPS0li5ciUWi4UzZ84QGBhoF4+Pjw/x8fHGNIbvv/+evXv3AvD3v/+dl19+GYCgoCD69esHYDdSVJCfnx/z5s0z2mrRogU9evQgPz+fpUuXMnHixJselTGbzcaXjfbt2/PXv/6V06dPY7FYOHr0KPfffz/bt283/i+rV69OYmKi8bgPP/yw3Repa9leY1euXOHgwYOEhYUVWVcJoNxS1yZ8Fy5csJubdC3bNySg0CrgmjVrOr09zL59+4zbnTt3Nv7IwdUP/JdeesnoHSyOBx54wG7ibcEey7y8PLtYa9asyZkzZ1i3bh1paWn8+OOPdsNLjr7tihRXwVXA+fn5ZGdnc/DgQebNm8exY8d4/fXXsVgsxupUV7w2fXx8iIiIoHnz5uzbt4/58+eTlpZGSkoKZ86cuW4bzr6HCg7t/fWvfzVulylThrlz59q1uXPnTuP29OnT7Xr3Cvrxxx+pX7++Xe/btSMJXbp0KVYCeOLECeP2tcOpBR/H9hwVfJxrE8C8vDzjefD29qZXr1525b169eLtt98mJyeH1NRUcnJy8Pb25q9//SuLFy9mx44dnD9/3uh1DAgIoHr16uzcuZMtW7bQuXNnNm7cCFz9YuxoEYGzz8+jjz7K/Pnz2b9/Px06dODee+815rJFRkbaJXImk4mwsDDCwsLIyMjgq6++Ii0tjR9++IHdu3cb9Ry9ZurVq2c3h7Vy5cpGAlhwyL3g3wZH8+0AWrdubddWvXr1qFu3Lunp6eTl5XHgwIFCvbjFVXBRopeXFzVq1DBGwWxx2eKHq/P6CiadTZo0ISgoyJhWcK2CfxdPnTp13ViUAMotVb58eePbKVx94V9vj6yCH5gF5wLB1W1anFUwkbx2Xoe7uzsVK1bk5MmTTrdnc21S6uXlZdwuOAH75MmTvPnmm6xevdpINIOCgqhTpw579uwpVF/kRjlaBQxX53cNHDgQuJoU2RJAV702Fy9ezJQpU4yeN29vbx566CFOnTrFpUuXiryfs+8h23vY09PzD7dlcnbhhi3WgvOuKlasaFfHNoTprILJxrVfeG2P4+HhUWjesaPPs3PnzhnDs4GBgZQpYz+Lq0yZMlSuXJnffvuN/Px8srKy8Pb2pl27dixevJjLly/z/fffGwlgixYtCA4OZufOnaSkpLBjxw7js7hDhw4Oz8fZ5+eVV16hcuXKJCUlcfLkSQ4fPszhw4dZtGgRbm5udO3albffftu4/86dO5k4caIx/9DNzY06depQvXp1I9lx9Lq7Nh43NzfjdsFkruD/VVGv32sXQ8HV59v2d+faXtMbce2XAEf/fwX/Pjl6HQQEBBSZABZ8LxRsxxElgHLLtWjRwlgxtXDhwiITwC1bthgvek9PTxo3bmz37bDggpA/UvBD/NpvSbahhhtR8MPnekaPHk1KSgplypRhwoQJdOnShcDAQD7//HNjCEekJD388MPG7TNnznD69Gn8/f1d8trcsmUL48ePx2q10rBhQ8aNG0dYWBgeHh5ERkZeNwF09j1Urlw5zp49S15eHhcvXrT7w/fbb7/h4+NjvM8LlvXr14+IiAiHbdq+VF6vF6W4XwwLjmjk5ubaJYG2c7BYLJw7d86urqPem4oVK1KmTBny8/M5fvw4+fn5dolNfn6+0ePo5uZmJDStWrXC29ubnJwcVq1aZcy9syWAM2fO5NChQ3Zz49q3b+/wfJx9fjw8PBgxYgTPPvssu3btYsuWLaSmprJp0ybOnz/Pl19+SUBAAOPGjSMrK4unn36as2fP4u/vz5tvvklERATlypXjxRdfLDLZ+aN4nI3VxtHK8YLP97XJ2424Nml35Hp/n8C+V/laBZPU642mgRaByG1gwIABxu0lS5Y4vDzVsWPH7P749OrVq9A3ZmfeWDZ169Y1bl97ibnFixff0PCvs86dO2esuqtcuTIDBw40hkNsPSwiJa3ga83WA+Wq1+aqVauM91Tv3r1p3LgxHh4eHD9+/IZ61h0puNHut99+a9zOz8/n2Wef5S9/+QvNmzfnxIkTdnuznT17lnbt2hk/hw8fZuvWrUYCDFC/fn2j/vLly+0ed8mSJcWKs+Af82u/WBaMy5nHcXd3p2HDhsDVVckLFy60K//f//1fY9P7hx9+2Bg69PHxMYZDly9fbgzVtmjRgsaNGxtfnhcvXgxcHea9doSlOH777TemTJnCqFGjmDx5MmFhYQwZMoRp06bZXTrNNpydkpJi9Dz+5S9/oX379pQrV478/Hy74dCStnHjRqMXGK4ueLLN4fT19aVWrVp/ShwF/z6tXr3abmh906ZNRc5hBPu58dfOl7yWegDllmvcuDH9+/cnMTERq9XKiy++yIoVK3jkkUcoW7YsP/30EwsWLODcuXMA3Hvvvbzwwgs39Zjt27encuXKnDhxgtTUVIYNG0bnzp3Zv39/sa6leCN8fHzw8PDAYrFw7NgxXn31VZo3b87mzZvtvoEXfNOL3KgDBw6wevVq43eLxcKRI0f4z3/+Yxz729/+ZiQLrnhtFuxB++ijj4CrCcunn35q7P92s6/vvn37GltIvfLKKxw4cIDatWvzzTffGEN2ISEhVK5cmXbt2uHv78/p06dZunQpZcqUoXXr1hw4cIC4uDiuXLlC+fLl6dixI3B1/lpsbCwWi4Xly5fj5eVFq1at2LlzJ5999lmx4rRtuwJXF0kEBQUZv//973835ty99957nDx5knvvvZcVK1bYbQ1TkNlsNvb6e/3110lPT+ehhx5iz549Rmwmk6nQFiDt2rXj22+/NYaQg4ODjSSvcePGbNy40SgrqvfPWffccw+ffvqp0dN7/vx5mjdvTnZ2NosWLTLq2R6/4Jf5lStXGs/bokWL7Kb9lPRnYm5uLk8++SSDBw8GYPbs2cb/yWOPPWa3Wr4khYeH8+CDD/LTTz9x9OhRBgwYQO/evTl+/DizZ8++7n0PHToEXB0RKzg/0xElgHJbmDBhAh4eHsyZMwer1cqaNWtYs2ZNoXqhoaHExsb+Ydf2H/H09OSll15i9OjRWK1W1q1bx7p164Cr38p//fVXLly4YLc4xFVsq9w+//xzAJKTk41NO21/fMHxcIRIcS1fvrxQ71JBlSpV4sUXXwRc99p89NFH+c9//sOlS5f4/fff7Xrvbe2cOXOG3NxcuzlQxdGuXTvji2NOTo6RaNoU3Dze19eX999/n+HDh2OxWPjyyy/tRhrc3d156623jM+V6tWr8/LLL/PGG28AV7fosCUutWvXxmq1cvDgQafiLDilJTU11W74uUuXLqxevZply5aRk5Njt+9g27ZtHX4Gdu7cmfT0dGbMmMHly5dJSEiwKy9Tpgwvv/xyoWHuRx55BDc3N2N0o2B5RESEkYjCzSeAfn5+xMTEMHr0aCwWCwsWLLDbIgau9oz+4x//AK4moCEhIcZii4Krrgu+7jIzM+1WrLva/fffz6FDhwrtX1i7dm1Gjx5dYo/ryCuvvMKgQYPIy8vjxx9/NJL+oKAgKlSoYGwAfe3Il20LswYNGuDh4XHdx9AQsNwW3NzcGDduHEuWLCEqKoo6derg6+uLh4cHgYGB/O1vf2Py5Ml8/vnnTl3k2hldunRh1qxZNGrUCG9vb/z9/YmKiiIhIcH4pmfbg8rVXn31VUaMGMF9992Hl5cXwcHBdO3alcWLFxt/hAr22oi4SpkyZYzXXO/evfn888/thvtc8dqsVasWCQkJtGzZknvuuQc/Pz9CQ0N57bXXjN77y5cvO0xwiuOVV15hypQpNGvWjHLlyuHj40Pt2rXp378/ixYt4t577zXqtm7dmkWLFtG9e3cCAwONz5Z27doxd+5cOnfubNd2VFQUH3/8MY0aNcLLy4uAgACio6OZP39+seaC1ahRwxjS+/777wuVT5o0ifHjx1OzZk08PDyMbaiut9L4+eefZ968eXTr1o2qVavi4eFBQEAAnTp1YsGCBcZ2JwX5+/vTuHFj4/cWLVoYtwuulr333nvtroB0ozp06MAXX3zB448/Ts2aNfHx8cHLy4tatWoZz4/tdefp6cns2bPp0aMHgYGB+Pj4UKdOHQYNGmQ3ZLxq1aqbjut6mjVrxn/+8x/Cw8Px9PQkICCAqKgokpKSXDL/rziaNGnC3Llz+ctf/oKvry9+fn5069aNpKQkuw6Qgn+jsrKyjB7Aa1/PjpisWmoopdCxY8fYsGEDlStXpnLlyjzwwANGb9/Fixdp3rw5FouF++67r8Q/dETk7vbZZ5/x1ltv4enpyZYtW0rsi6UU38KFC419+aKjo52+xnxJunDhAl9//TWVK1cmICCAkJAQY4pGfn4+jzzyCJmZmXh5efHjjz8ai12WL1/OCy+8gKenJxs2bPjDkTINAUupdPHiReNawXD1jd+0aVMuXbrEokWLjCGHZs2a3aoQReQu0aNHDz788EPOnz/PunXrnOqdkdLLZDLxxhtvGPNlu3TpQvv27cnLy+Obb74xFqo0btzYbqWzbWPxbt26OTVNSgmglEq1a9emQ4cORu/e3LlzC20ee8899/DMM8/civBE5C5Srlw5RowYwbvvvsu8efOUAMp1lS1blv79+xsLEh3N4/X09DSuaANXr2u9bt06fH19jbmVf0RzAKXU+uCDD5gwYQKNGjWiQoUKuLm54eXlxX333Uffvn1ZuHCh3RwiEZEb1a9fP0JCQtiyZYvdVUxEHBk3bhzvvPMOLVq0oFKlSri7u+Ph4UFQUBCPPfYYX3zxBY0aNTLqf/LJJ1y5coWhQ4c6fUUszQEUERERKWXUA3gdly9fJiMjwxiHFxG5U+jzS0SuRwngdWRmZtK2bVu7ncFFRO4E+vwSketRAigiIiJSyigBFBERESlllACKiIiIlDJKAEVERERKGSWAIiIiIqWMEkARERGRUkYJoIiIiEgpowRQREREpJRRAigiIiJSyigBFBERESlllACKiIiIlDJKAF3Mqguv35b0vIg4R++V28+f9Zzk6bm/7ZTkc+JeYi2XUiZ3d46+GXurw5BrVHttxK0OQeSOoM+w28+f9fnl6e5Ohzkz/5THEuesGjSsxNpWD6CIiIhIKaMEUERERKSUUQIoIiIiUsooARQREREpZZQAioiIiJQySgBFREREShklgCIiIiKljBJAEZGbcOjQIcaMGUPLli0JDQ3lkUce4e233+b8+fN29fbu3cvQoUNp3rw5Dz/8ME8//TT79u1z2GZycjLdu3enYcOGREZGMmnSJLKzs/+M0xGRUkIbQYuI3KBjx47xxBNPYLFYePLJJwkODiY1NZW5c+eSkpJCcnIyPj4+pKenEx0dTfny5Rk8eDBubm4kJCTQt29fkpOTuf/++402Y2NjmTZtGm3atCEqKor09HTi4+PZvXs38fHxlCmj7+0icvOUAIqI3KD333+frKwsFixYQFhYGAB9+/blgQce4J133mH+/PkMGjSImJgY8vPzSUpKolq1agB07NiR7t27ExMTw+zZswE4evQoM2fOpF27dsTGxmIymQAIDg5m0qRJrFixgi5dutyakxWRu4q+SoqI3KDNmzfz4IMPGsmfTY8ePQDYunUrJ0+eZMOGDXTq1MlI/uBqUtepUyc2btzIqVOnAFi6dCkWiwWz2WwkfwD9+/fHy8uLJUuW/AlnJSKlgRJAEZEblJyczOTJkwsdP336NADu7u6kpqYCFEoSAUJDQ8nPz2fXrl0ARdb18vIiJCSEnTt3ujR+ESm9NAQsInKDgoKCHB6fM2cOAM2bNyczMxPArvfPpkqVKgBkZGQAkJmZScWKFfH29nZYNy0tjezsbHx8fOzKsrKyyMrKsjtme1wREUeUAIqIuNCyZctITk4mKCiIXr16kZCQAFAoaSt4zLbC9/z58/j6+jps15YUOkoAExISiI2Nddk5iMjdTwmgiIiLLF26lHHjxuHj48PUqVPx9vbGarUC2M3ps7Edc1RWVF1HzGYzPXv2tDuWmZlJdHR0ccIXkVJECaCIiAvExcXxwQcf4OvrS1xcHKGhoQBGj56jffxsx8qXL2/UPXbsmMP2c3Jy7OoW5Ofnh5+f382fhIiUGloEIiJyE/Lz83njjTeYPHky/v7+fPrppzRp0sQoDw4OBuD48eOF7mubp2ebCxgcHMzp06fJy8tzWLdSpUp4eHiUxGmISCmjBFBE5AZZrVZefvllkpKSqFWrFgsWLDB6/mxCQ0MxmUykpaUVun9aWhomk8lY9dugQQOsViu7d++2q5ebm0t6errDlcQiIjdCCaCIyA2Kj49n4cKF1KtXj6SkJGrUqFGoTmBgIM2aNWPZsmV2K3MzMjJYuXIlkZGRVKxYEYDOnTvj5uZGfHy8XRuJiYnk5uYa+wuKiNwszQEUEbkB586dM1betm/fnvXr1xeqExgYSEREBOPHjycqKoqoqCjMZjP5+fkkJCTg5ubGiy++aNSvUaMGgwcPJi4ujmHDhtG2bVv27NnD/Pnzad26NR06dPjTzk9E7m5KAEVEbsCOHTu4dOkSQJFbsERERBAREUH9+vVJTEzkww8/5N///jeenp6Eh4czevRo6tata3ef0aNHExAQQFJSEhMnTiQwMJAhQ4YwfPhwXQdYRFxGCaCIyA2IjIxk3759TtcPCwsrNLTriMlkwmw2YzabbyY8EZHr0tdJERERkVJGCaCIiIhIKaMEUERERKSUUQIoIiIiUsooARQREREpZZQAioiIiJQySgBFREREShklgCIiIiKljBJAERERkVJGCaCIiIhIKaMEUERERKSUuS0SwF9++YURI0bwl7/8hYcffpgBAwawadOmQvW2bt1K//79ady4Mc2aNeOFF17gyJEjhepdvnyZuLg4OnbsSFhYGO3btycuLo7Lly//GacjIiIiclu75Qng/v376dOnDz/++CP9+vXjhRdeICsri0GDBrFixQqj3qZNmxg0aBBnzpxh5MiRmM1mNmzYQN++fTlx4oRdm6+++iqTJ08mNDSUl156ibCwMCZPnsyrr776Z5+eiIiIyG3H/VYHMGnSJJgWCpoAACAASURBVPLy8vj888+pU6cOAI8//jidO3fm/fffp1OnTlitVt566y0qV65MUlIS5cuXB6BVq1Y88cQTxMbGMnHiRAB27NjBwoUL6d+/P6+88goAffv2xc/Pj3nz5vHEE08QHh5+a05WRERE5DZwS3sArVYrnp6edOrUyUj+AMqWLUt4eDhHjhzh7NmzpKWlsX//fnr16mUkfwANGzakadOmLF++3BjeXbx4MQADBw60e6whQ4YA8OWXX5bwWYmIiIjc3m5pD6DJZGLq1KmFjlssFvbt24efnx/ly5dnx44dAISFhRWq26BBA7Zs2cKhQ4eoW7cuqampBAQEEBwcbFcvKCgIf39/du7cWTInIyIiInKHuOVzAAs6e/Ys27ZtY/jw4Rw6dIjnn38eNzc3MjMzAahWrVqh+1SpUgXAWAySmZlJ1apVHbZfpUoVh4tGREREREqTWz4HsKBhw4bx448/AtCuXTt69uwJwIULFwDw8fEpdB/bsUuXLgFw/vx5fH19Hbbv7e1Ndna2w7KsrCyysrLsjtkSTxEREZG7yW2VAA4cOJAhQ4awZcsW5s6dS58+fZg3bx5WqxW4OmR8LdsxR2WO6hZVLyEhgdjY2JuIXkREROTOcFslgJ06dQKu9v7VrFmTiRMnkpiYaPTo5eTkFLqPrUfPtjjE19fXYT3b/cuVK+ewzGw2Gz2ONpmZmURHR9/YyYiIiIjcpm6rOYAFPfroowDs2rXLWNBx7NixQvVsw7S2eX/BwcEO69nqFjU/0M/Pj+DgYLufouqKiIiI3MluaQKYkZFBhw4deP311wuVXbx4EQAvLy8aNGgAQFpaWqF6aWlplC9fntq1awMQGhrKsWPHCiWBR44c4fTp0zRs2NDVpyEiIiJyR7mlCWC1atXIy8tj6dKlZGRk2JXNnDkTgPbt29OwYUNq1KjBggULjAUhADt37mTr1q08+uijxty+7t27AzBnzhy79mbPng3AY489VmLnIyIiInInuKVzAN3c3HjjjTd49tlniYqKIjo6mnLlyrF27Vo2bNhAp06d6Nq1KyaTiVdffZVhw4YRFRVF3759OXfuHHPmzKFatWoMHz7caLNp06Z07dqVOXPmcPr0aZo2bcrmzZtZunQpffr0MXoTRUREREqrW74I5K9//SuJiYl89NFHxMXFYbFYqFWrFq+88grR0dFGz15kZCRxcXFMnz6dmJgYypUrR+vWrRkzZgyVK1e2a/O9996jZs2aLF68mK+//prq1aszduxYnnrqqVtxiiIiIiK3lVueAAI0btyYTz755A/rtW7dmtatW/9hPU9PT0aNGsWoUaNcEZ6IiIjIXeW2XQUsIiIiIiVDCaCIiIhIKaMEUERERKSUKVYCmJ+fz/Lly40rbZw9e5Z//OMftGrViujoaHbt2lUiQYqIiIiI6zidAGZnZ2M2mxkzZgz79+8HYNy4caxcuZJTp07xww8/YDabOXz4cIkFKyIiIiI3z+kEcObMmWzduhWAAwcO8Ntvv/Hf//4XAKvVCsClS5eYNWtWCYQpIiIiIq7idAL47bffAnDvvfcSHh7O+vXrAfD392fjxo088sgjWK1WNm/eXDKRioiIiIhLOJ0AZmRkYDKZGDRoEDVq1OCHH34AICIigkqVKtGxY0eAQtfgFREREZHbi9MJoMViAaBChQoA/PDDD5hMJho1agRcvawbXN2EWURERERuX04ngLbE7+DBg+zdu9fo6WvatCkAa9asAaBq1aqujlFEREREXMjpBLBBgwZYrVZmzpzJ4MGDAQgKCiIkJIQJEyawYsUKTCYTERERJRasiIiIiNw8pxPAQYMGUaZMGfLy8jh16hQmk4mnnnoKgLy8PADKli2L2WwumUhFRERExCWcTgCbNWvGtGnTaNiwIffffz+jR4/mySefBKBWrVoEBQUxa9YsgoODSyxYEREREbl57sWp3LZtW9q2bVvoeHR0NM8++yxlyujKciIiIiK3O6cztgEDBmA2m9m+fXuhsooVK/LVV1/RsWNHRowY4dIARURERMS1nO4B3LJlCyaTiTNnzjgsz83N5ddff+Xs2bMuC05EREREXK/IBHDmzJl8/vnnhY6/9tprvPPOO3bHrFYrJ06cAODy5csuDlFEREREXKnIBLBPnz7Mnj2bixcvAmAymQA4ffp0kY2ZTCbq16/v4hBFRERExJWKnAPo7+/P0KFDsVqtWK1W47jtd0c/derU4Z///OefEriIiIiI3JjrzgEcPHgwf//737FarbRs2RKTycT7779P69at7eqVKVMGHx8fXQZORERE5A5w3QSwTJkyVKxYEYB3330XuLofoO2YiIj8n9TUVPr27UtiYiJNmjSxK2vRooXDRXRubm7s2bPH7lhycjKJiYkcPnyYChUq0K1bN0aMGIGPj0+Jxi8ipYfTq4B79uxZknGIiNzRDh8+zIgRI8jPzy9UduzYMc6cOUPXrl2JjIy0K7t2/9TY2FimTZtGmzZtiIqKIj09nfj4eHbv3k18fLz2WxURlyjWRtDbt29n7ty5HD58mNzcXLu5gTYmk4kvv/zSZQGKiNzu1q5dy/jx44vcBmvv3r0AdOrUiQ4dOhTZztGjR5k5cybt2rUjNjbWWHwXHBzMpEmTWLFiBV26dHH9CYhIqeN0Arh27Vqee+45h0mfjdVqNT6wRERKgzFjxrB06VLq1q1Lq1atWLp0aaE6+/btA+D++++/bltLly7FYrFgNpvtPkv79+/P1KlTWbJkiRJAEXEJp8cSZsyYQX5+/nVXAYuIlDb79+/nhRdeYNGiRdSsWdNhnX379uHp6cm9994LYGyvda3U1FQAwsLC7I57eXkREhLCzp07XRe4iJRqTvcA7tu3D5PJRJUqVXj55ZepU6cOXl5eJRmbiMhtLzk5+Q93QNi7dy/33HMPL730Et988w0XL16kUqVKREVFMXz4cNzdr34UZ2ZmUrFiRby9vQu1UaVKFdLS0sjOzi60GCQrK4usrCy7Y5mZmTd5ZiJyN3M6ASxbtix5eXmMGjWK9u3bl2RMIiJ3jD9K/vLy8jh06BCXL1/GYrHw3nvvcenSJb766itiY2PZv38/H374IQDnz5/H19fXYTu2pNBRApiQkEBsbKwLzkZESgunE8DIyEgWL17MpUuXSjIeEZG7Sl5eHv/4xz+oVKkSjz/+uHG8R48ePPfcc3z99dc88cQTREREXLed682vNpvNhXZqyMzMJDo6+uaCF5G7ltNzAEePHk1gYCCzZs0iPT29JGMSEblrlCtXjqefftou+bPp168fABs3bgTA19eXnJwch+3YjpcvX75QmZ+fH8HBwXY/VatWddUpiMhdyOkewA8++ICaNWuSkpLCY489RkBAAH5+foX2pNI2MCIizgkICAD+b1FIcHAwP/30E3l5eYWGljMzM6lUqRIeHh5/epwicvdxugdw0aJFbNmyBZPJhNVq5eTJkxw4cIBffvnF+Pn555/5+eefSzJeEZE7yurVq+nYsSNJSUmFyvbv3w/AfffdB0CDBg2wWq3s3r3brl5ubi7p6emFVgeLiNyoYm0pX3C7F20DIyLyx+rWrcvhw4eZO3cuubm5xvHs7GxmzJiBl5cXnTt3BqBz5864ubkRHx9v10ZiYiK5ubn06NHjT41dRO5eTg8Br1mzpiTjEBG5K9WsWZMhQ4YQFxdHnz59ePzxx8nNzWXRokUcPHiQt956iypVqgBQo0YNBg8eTFxcHMOGDaNt27bs2bOH+fPn07p16+teRUREpDicTgCDgoJKMg4RkbvWmDFjqFWrFp999hn/+te/cHd3p0GDBrz88su0atXKru7o0aMJCAggKSmJiRMnEhgYyJAhQxg+fLiuAywiLlOsawHb7N27l+3bt3Ps2DG6deuGv78/np6eDleniYiUFiNHjmTkyJEOyx5//HGHK4GvZTKZMJvNmM1mV4cnImIoVgJ46NAhxo8fb1yuCK5esmjHjh1MmjSJt99+W0MUIiIiIrc5p8cTjh8/Tr9+/UhNTS206GP//v1kZWUxevRoXatSRERE5DbndAI4c+ZMTp48iaenJ88995xdWZUqVXBzc+PKlSvMnj3b5UGKiIiIiOs4nQCuW7cOk8nEoEGDCs1xGTRoEE8//TRWq9VueFhEREREbj9OJ4AnTpwAICQkxGF5nTp1ADhz5owLwhIRERGRkuJ0Aujn5wdQ5HWAv/vuOwAqVqzogrBEREREpKQ4vQo4IiKCpUuXMmfOHNzd/+9umzZtYsWKFXz11VeYTCaaN29eIoGKiIiIiGs4nQCOGDGCtWvXcunSJaZPn47JZAJg7ty5Rh0PDw+GDh3q+ihFRERExGWcHgKuWbMmH3/8MYGBgQ6vA+zv789HH31kzAUUERERkdtTsTaCbtKkCd988w1r165l586dZGVlUbZsWR566CHatWuHj49PScUpIiIiIi5S7EvBeXp60rFjRzp27FgS8YiIiIhICSsyAVy1ahUAjRo1onLlysbvztDl4ERERERuX0UmgKNGjcJkMhEbG0vbtm2N3/+IyWRiz549Lg1SRERERFynWEPABa//KyIiIiJ3piITwB49emAymahevbrd7yIiIiJyZysyAXzvvfeu+7uIiIiI3JmcHgJevHgxAC1atKBq1aqFynfu3MmmTZsICAigV69erotQRERERFzK6QRw/PjxxqIQRwnggQMHmDJlCsHBwUoARURERG5jRSaAK1eu5Oeffy50fPny5fz000+Fjn/33XcAnDhxwoXhiYiIiIirFZkABgQE8PzzzxsLP2z/Ll++vMjGTCYT1apVc3GIIiIiIuJKRV4LuHHjxnTo0MG41q+No+sAF/zp3bv3nxK4iIiIiNyY684BHD9+POHh4QC8//77mEwmevbsyf33329Xr0yZMnh7e/PQQw/RoEGDkotW5DaWd/kynu7FvrqilDA9LyIihV33U7F69eo89dRTAKxduxaA3r178/DDD5d8ZCJ3GE93dzrMmXmrw5BrrBo07FaHICJy23H6a3FiYqJT9SwWCx4eHsUK4tChQ0ybNo1NmzaRlZVF5cqVadu2Lc8//zzly5c36u3du5cpU6awY8cOLBYLjRs3ZuzYsdSrV69Qm8nJySQmJnL48GEqVKhAt27dGDFiBD4+PsWKTURERORuU6xxkbNnz7Ju3TqOHz+OxWKxmxtosVg4evQo69evZ/PmzU63eezYMZ544gksFgtPPvkkwcHBpKamMnfuXFJSUkhOTsbHx4f09HSio6MpX748gwcPxs3NjYSEBPr27UtycrLdsHRsbCzTpk2jTZs2REVFkZ6eTnx8PLt37yY+Pp4yZYqc+igiIiJy13M6AczIyKBv376cOnWqyDpWq7XYl4t7//33ycrKYsGCBYSFhQHQt29fHnjgAd555x3mz5/PoEGDiImJIT8/n6SkJGOlcceOHenevTsxMTHMnj0bgKNHjzJz5kzatWtHbGysEU9wcDCTJk1ixYoVdOnSpVgxioiIiNxNnO4Kmz59OidPnrzuCmCg0AKRP7J582YefPBBI/mz6dGjBwBbt27l5MmTbNiwgU6dOtltMxMcHEynTp3YuHGjkZguXboUi8WC2Wy2S0b79++Pl5cXS5YsKVZ8IiIiIncbpxPAlJQUTCYToaGhfPzxx5QtW5bw8HDmzJnD2LFjcXd3x8PDg3fffbdYASQnJzN58uRCx0+fPg2Au7s7qampAIWSRIDQ0FDy8/PZtWsXQJF1vby8CAkJYefOncWKT0RERORu43QCaEvIoqOjiYyMpGnTpvz+++9EREQwZMgQ+vXrh8Vi4aOPPipWAEFBQdSqVavQ8Tlz5gDQvHlzMjMzARxuMl2lShXg6hA1QGZmJhUrVsTb29th3dOnT5OdnV2oLCsri4yMDLsf2+OKiIiI3E2cngN45coV4GpPGkDDhg3573//y9GjR6lWrZqx/5+jy8QV17Jly0hOTiYoKIhevXqRkJAA4HAFr+2YLak7f/48vr6+Dtu1JYXZ2dmF2kpISCA2NvamYxcRERG53TndA2jraVuyZAl5eXk0bNgQq9XK/PnzsVgsrFixAoAzZ87cVEBLly7ln//8Jz4+PkydOhVvb29jfqGjBSbXXqrueq5Xx2w2s2bNGrufuXPn3uBZiIiIiNy+nE4A27Rpg9VqZf369cyaNYuGDRvi6elJXFwc4eHhrF69GpPJRFBQ0A0HExcXx9ixY/Hy8mLWrFmEhoYCGD16joZubcds+wX6+vqSk5PjsH3b8YJ7C9r4+fkRHBxs91O1atUbPhcRERGR25XTCeCIESOMFb4hISGULVuW7t27Y7VauXLlitFL17dv32IHkZ+fzxtvvMHkyZPx9/fn008/pUmTJkZ5cHAwAMePHy90X9s8PVsPZXBwMKdPnyYvL89h3UqVKhV7o2oRERGRu4nTCaC/vz8LFixgwoQJxgrbV155hZ49e1KuXDmqVavG888/T79+/YoVgNVq5eWXXyYpKYlatWqxYMECo+fPJjQ0FJPJRFpaWqH7p6WlYTKZjJgaNGiA1Wpl9+7ddvVyc3NJT093uJJYREREpDQp1iUxfH19GTBggNHb5uPjw7vvvsu2bdtYu3Ytw4YV/5qb8fHxLFy4kHr16pGUlESNGjUK1QkMDKRZs2YsW7bMbmVuRkYGK1euJDIykooVKwLQuXNn3NzciI+Pt2sjMTGR3NxcY39BERERkdKqWJeCs9m2bRt79+7l4sWLlC9fntDQ0BvqWTt37pyx8rZ9+/asX7++UJ3AwEAiIiIYP348UVFRREVFYTabyc/PJyEhATc3N1588UWjfo0aNRg8eDBxcXEMGzaMtm3bsmfPHubPn0/r1q3p0KHDjZyyiIiIyF2jWAng9u3bmTBhAocPHy5UVrduXWJiYqhfv77T7e3YsYNLly4BFLkFS0REBBEREdSvX5/ExEQ+/PBD/v3vf+Pp6Ul4eDijR4+mbt26dvcZPXo0AQEBJCUlMXHiRAIDAxkyZAjDhw/XdYBFRESk1HM6Ady7dy+DBw8mJyfHWPBR0M8//0z//v1JTk6mTp06TrUZGRnJvn37nA42LCys0NCuIyaTCbPZjNlsdrptERERkdLC6QTwww8/NLZcqVq1Ku3atSMgIICTJ0+ydu1ajhw5wsWLF/nggw+YPn16iQUsIiIiIjfH6QRw69atmEwmGjZsyKeffoqnp6dRNn78eAYOHMi2bdtISUkpkUBFRERExDWcnhDn5uYGQJ8+feySPwB3d3d69erl2shEREREpEQ4nQA2b94cgFOnTjksP3LkCAAtW7Z0QVgiIiIiUlKcTgDHjh1LuXLlmDVrFlu3brUrW7VqFZ988gmBgYG89NJLLg9SRERERFzH6TmAM2fOpE6dOuzYsYMBAwYQEBBAhQoVOHbsGOfPnwegbNmyPP3003b3M5lMfPnll66NWkRERERumNMJ4KJFizCZTJhMJqxWKydPnuTkyZNYrVZMJhNwdXi44BBxwTIRERERuT0UayPogvv/FXVbRERERG5vTieAa9asKck4RETueKmpqfTt25fExESaNGliV7Z3716mTJnCjh07sFgsNG7cmLFjx1KvXr1C7SQnJ5OYmMjhw4epUKEC3bp1Y8SIEfj4+PxZpyIidzmnE8CgoKCSjENE5I52+PBhRowYQX5+fqGy9PR0oqOjKV++PIMHD8bNzY2EhAT69u1LcnIy999/v1E3NjaWadOm0aZNG6KiokhPTyc+Pp7du3cTHx+vy1mKiEsUawhYREQKW7t2LePHj+fs2bMOy2NiYsjPzycpKYlq1aoB0LFjR7p3705MTAyzZ88G4OjRo8ycOZN27doRGxtrzKEODg5m0qRJrFixgi5duvw5JyUidzV9lRQRuQljxoxh2LBhBAQE0K1bt0LlJ0+eZMOGDXTq1MlI/uBqUtepUyc2btxoLJ5bunQpFosFs9lst4Cuf//+eHl5sWTJkpI/IREpFZQAiojchP379/PCCy+waNEiatasWag8NTUVgLCwsEJloaGh5Ofns2vXruvW9fLyIiQkhJ07d7o4ehEprTQELCJyE5KTkwtdHrOgzMxMALveP5sqVaoAkJGRYdStWLEi3t7eDuumpaWRnZ1daDFIVlYWWVlZDh9XRMSRIhPAY8eOUbFixet+sImIlHZ/9Bl54cIFAIcreG3HsrOzATh//jy+vr4O27ElhY4SwISEBGJjY4sXuIiUakUmgCNHjmTPnj3ExMTQtWtX48Ola9eu1KpV608LUETkTmbbJ9XRpvi2Y85smH+9OmazmZ49e9ody8zMJDo6ujihikgpUmQCeODAAa5cuYKbmxuAsSLtwQcfVAIoIuIkW4+erZevINux8uXLG3WPHTvmsJ2cnBy7ugX5+fnh5+fnknhFpHQoMgG0WCwAJCUlcenSJeP45s2bjWv/FqVHjx4uCk9E5M4WHBwMwPHjxwuV2ebp2eYCBgcH89NPP5GXl1doaDkzM5NKlSrh4eFRwhGLSGlQZAJ477338ssvv7Blyxa2bNliDD989tln123QZDIpARQR+f9CQ0MxmUykpaXRu3dvu7K0tDRMJpOx6rdBgwasWrWK3bt306hRI6Nebm4u6enptGzZ8k+NXUTuXkVuA9O/f3+sVusN/YiIyFWBgYE0a9aMZcuW2a3MzcjIYOXKlURGRlKxYkUAOnfujJubG/Hx8XZtJCYmkpubqy/XIuIyRfYA9unTh+rVq7N161YuXrzIZ599hslkok2bNsaQhoiI/LHx48cTFRVFVFQUZrOZ/Px8EhIScHNz48UXXzTq1ahRg8GDBxMXF8ewYcNo27Yte/bsYf78+bRu3ZoOHTrcwrMQkbvJdfcBbNWqFa1atQKuDv1arVb69OlD27Zt/5TgRETuBvXr1ycxMZEPP/yQf//733h6ehIeHs7o0aOpW7euXd3Ro0cTEBBAUlISEydOJDAwkCFDhjB8+HBdB1hEXMbpjaDXrFkDQKVKlUosGBGRO9nIkSMZOXKkw7KwsLBCQ7uOmEwmzGYzZrPZ1eGJiBicTgCDgoIAOHXqFDNmzGDLli2cPXsWf39/IiIi6NevHxUqVCixQEVERETENYp1KbiUlBRGjRpld8mhQ4cOsX37dubNm8dHH31EeHi4y4MUEREREddxekLJ0aNHGTFiBFlZWQ5X/p4+fZphw4Y53OtKRERERG4fTieAs2fPNjaA7t69OwsWLGDdunXMnz+f7t27A3Du3Dlmz55dMpGKiIiIiEs4PQT83XffYTKZaNu2LZMmTTKOV61alfDwcLKzs1m9ejXr1q3jpZdeKpFgRUREROTmOd0DaNvAtF27dg7LbccLbnQqIiIiIrcfpxPAcuXKAXDixAmH5ba5f44uVC4iIiIitw+nE8B69ephtVqJj4/nwIEDdmUHDx5kzpw5mEwmHnzwQZcHKSIiIiKu4/QcwF69erFp0ybOnj3LY489RosWLahatSqZmZls3rwZi8WCyWSiZ8+eJRmviIiIiNwkpxPAbt268fXXX7NmzRouX77Mhg0bjDKr1QpcvXRc165dXR+liIiIiLhMsS4sOXXqVJ555hl8fHzs9gD08vLiySefZPr06SUVp4iIiIi4SLGuBOLm5sbo0aMZOnQoO3fu5MyZM1SqVImHHnrIWCQiIiIiIre3YiWANmXLliUiIsLVsYiIiIjIn6BYQ8AiIiIicudTAigiIiJSyigBFBERESlllACKiIiIlDJKAEVERERKGSWAIiIiIqXMDW0D48jhw4dZs2YNAIMGDXJVsyIiIiLiYi5LANPT04mJiaFMmTJKAEVERERuYy4fArZdF1hEREREbk8u6wGsXbs2I0aMcFVzIiIiIlJClACKiIiIlDI3nAD++uuvnD17Fn9/f2rUqOHKmERERESkBBU7AUxOTmb69OkcP37cOFa9enVGjhxJjx49XBqciIiIiLhesRLAd955h8TERMB+sceRI0eYMGEC+/fvZ8yYMa6NUERERERcyulVwOvXr+fTTz8FriZ/tWrVokWLFtSsWdM4Nnv2bDZt2lQigYqIiIiIazjdAzh37lwAvL29mT59Oi1btjTKNm3axLPPPktOTg4JCQlERES4PlIRERERcQmnewB37tyJyWRi0KBBdskfQEREBAMHDsRqtZKamuryIEVERETEdZxOAC9evAhgDPley3b80qVLNx2UiIiIiJQcpxPAypUrA7B582aH5SkpKQAEBga6ICwRERERKSlOJ4ARERFYrVYWLVrE1KlTOXXqFACnTp1i2rRpLFq0CJPJVGh4WERERERuL04vAhk4cCBffvklFouFGTNmMGPGDNzc3Lhy5QpwdRWwp6cn/fv3L7FgRUREROTmOd0DWLduXd5++208PDywWq1YrVYuX75s3HZzc+PVV1+lbt26NxxMamoqDz74INu2bStUtnfvXoYOHUrz5s15+OGHefrpp9m3b5/DdpKTk+nevTsNGzYkMjKSSZMmkZ2dfcNxiYiIiNxNirUR9KOPPsoDDzzAJ598wubNmzlz5gwBAQE0atSIAQMG0LBhwxsO5PDhw4wYMYL8/PxCZenp6URHR1O+fHkGDx6Mm5sbCQkJ9O3bl+TkZO6//36jbmxsLNOmTaNNmzZERUWRnp5OfHw8u3fvJj4+njJlnM55RURERO5Kxb4UXEhICDExMS4NYu3atYwfP56zZ886LI+JiSE/P5+kpCSqVasGQMeOHenevTsxMTHMnj0bgKNHjzJz5kzatWtHbGwsJpMJgODgYCZNmsSKFSvo0qWLS2MXERERudPc8u6wMWPGMGzYMAICAujWrVuh8pMnT7JhwwY6depkJH9wNanr1KkTGzduNBakLF26FIvFgtlsNpI/gP79++Pl5cWSJUtK/oREREREbnNF9gDGxsbecKMj2mR2WwAAFvZJREFURoxwuu7+/ft54YUXeOqpp/j4448Llds2lg4LCytUFhoaysKFC9m1axeRkZFF1vXy8iIkJISdO3cW5zRERERE7krXTQAL9qIVR3ESwOTkZDw9PYssz8zMBLDr/bOpUqUKABkZGUbdihUr4u3t7bBuWloa2dnZ+Pj4OB2fiIiIyN3munMArVZrsRssbtJ4veQP4MKFCwAOkzbbMdsK3/Pnz+Pr6+uwHVtSWFQCmJWVRVZWlt0xW/IpIiIicjcpMgF89913//DOa9asYfXq1UbSZ7VaqVGjhuui4/+SUEeJpe2YM0nnH9VJSEi4qWFvERERkTtFkQlgz549i7zT77//zptvvsl///tfTCYTVqsVDw8PnnrqKZ599lmXBmjr0XO0j5/tWPny5Y26x44dc9hOTk6OXd1rmc3mQuecmZlJdHT0jQUuIiIicpsq1jYwly9fJj4+nhkzZpCTk2P0zjVr1ozXX3+dOnXquDzA4OBgAI4fP16ozDZEa5sLGBwczE8//UReXl6hoeXMzEwqVaqEh4eHw8fx8/PDz8/PlaGLiIiI3Jac3gYmJSWFxx57jClTppCdnY3VasXf35/33nuPTz/9tESSP7i60tdkMpGWllaoLC0tDZPJZKz6bdCgAVarld27d9vVy83NJT093eFKYhEREZHS5g8TwFOnTjF27FgGDhzIgQMHsFqtmEwm+vTpw9dff02PHj1KNMDAwECaNWvGsmXL7BZlZGRksHLlSiIjI6lYsSIAnTt3xs3Njfj4eLs2EhMTyc3NLfFYRURERP5fe/cf1VV9x3H8hQgoQ0PmD0ppKvhFUOaaWpmgWcDWoTFRK38hbmVgQR21ztGalmuLlZumc00tww4apZI/5q9EBZ1CmZulQoIxwR8VeUR+69C4+8PDdyGg/JDv1773+Tinc+xz7/1+39f79X1e33s/33t/CBq9BGwYhtasWaPFixeroqLCerm3f//+mj9/fqse+9Zcs2fP1oQJEzRhwgTFxMSopqZG7777rpydnfX8889b1/Px8dHjjz+uFStWKC4uTg8++KBycnL0/vvvKyQkROHh4TarGQCu9cwzz+ijjz5qcNmaNWs0ZMgQSdKnn36qJUuWKCcnR87Ozho+fLiee+459ezZ05blAnBgjQbAsWPH6osvvqjzK9yAgABNnDhRJ0+e1MmTJxt90Zt9pi0wMFDJycl64403tHjxYrm6uupnP/uZZs6cKT8/vzrrzpw5U127dlVKSormz5+v7t2764knntD06dN5DjAAu8rNzVX//v3129/+tt6yvn37SpKysrI0bdo09e7dWwkJCaqsrNSqVat06NAhffjhh+rWrZutywbggBoNgDk5OXJycqpz+5QvvvhCc+fOve4LOjk5tTgAJiQkKCEhocFlP/3pT+td2m3s/WNiYhQTE9OiGgCgLVy8eFGnTp3SpEmT9Otf/7rBdQzD0CuvvKJu3bopJSXFeteC4OBgPfbYY1q6dKnmz59vy7IBOKjrnhIzDKNF/wEA6srLy1NNTY369evX6DpHjx5Vfn6+xo4dW+eWVYMGDdLQoUO1bds2XblyxRblAnBwjZ4BbM7j3AAA13f8+HFJsk5bqaqqUocOHepMTfnss88kNfzs86CgIB08eFAFBQX1pr4AQHMRAAHABmoD4LZt2/Tss8/q3Llz6tixo8LCwjRnzhx5eXk16dnnZ8+erRcAeZQlgOZq1o2gAQAtk5ubK+nqZd5Zs2bJ3d1dmZmZWrt2rY4cOaL169c36dnnVVVV9ZbxKEsAzUUABAAbiIqKUnBwsJ588km1b3+19f7iF79Qnz59lJiYqHfeeafFzz7nUZYAmosACAA28MgjjzQ4PnHiRC1YsEAHDhzQXXfdJen/zy7/vmufff59PMoSQHNxYzwAsCNXV1d16tRJlZWV1mefFxUV1Vuvdk6ft7e3TesD4JgIgADQxr7++mv96le/0uzZs+stKy4u1oULF3TnnXcqKChIkhp99nmnTp2sN4wGgNYgAAJAG+vRo4cqKiq0ffv2ek9RWrRokSRpzJgxGjRokHx8fPTBBx9YfxAiSUeOHNGnn36qyMjIBucAAkBzMQcQANpYu3bt9NJLL+mpp57SpEmTNHHiRHl6emrPnj06cOCARo8erbCwMEnS3LlzFRcXpwkTJmj8+PEqLS1VUlKSbr/9dk2fPt3OewLAURAAAcAG7r//fiUnJ+vNN99UUlKSqqur1adPH82bN08TJkywrjdy5EitWLFCf/vb3/Taa6/Jw8NDISEhmjVrFs8BBnDTEAABwEYGDx6slStX3nC9kJAQhYSE2KAiAGbFHEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEAAAACTIQACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMk4bAAsLi7WvHnzNHLkSA0aNEhjx47Vzp077V0WADQJPQxAW3LIAFhZWampU6dqw4YNioiI0OzZs+Xs7KyEhARt3LjR3uUBwHXRwwC0tfb2LqAtrFmzRrm5uVq6dKnCwsIkSWPHjtVjjz2mxMREhYeHy93d3c5VAkDD6GEA2ppDngHcuHGjevbsaW2ckuTq6qopU6aopKRE+/bts2N1AHB99DAAbc3hAmB5ebn+85//aNCgQfWWBQUFSZKOHDli67IAoEnoYQBsweEuARcVFckwDHl7e9db1qNHD0nS2bNn6y0rKytTWVlZnbHa9b755ptm1fBtRdmNV4JNfXfmjE3e50pJqU3eB013pgXH3tvbW+3b26c9tqSH3cz+JdHDbjW26l8SPexW05L+JTWthzlcACwvL5ekBufHdOjQQZJ08eLFesveffddLV26tMHXnDRp0k2sEHax6V17VwA7eXDZO83eZvfu3erVq1cbVHNjLelh9C8HR/8yrZb0L6lpPczhAqBhGJIkJyenestqxxpaFhMTo6ioqDpj1dXVOn36tHr37i1nZ+c2qPbW9c0332jSpElas2ZNg2ci4Lg49rLrfrekh9G/6uIzbF4c+6uasu8OFwBrvzU3dJbv0qVLkiQPD496yzp37qzOnTvXG+/bt+9NrvCHxdvb225nQmBfHHv7aEkPo381jM+weXHsb8zhfgRSe8CLiorqLaudC2PmbwUAbm30MAC24HAB0MPDQ3369NHRo0frLav95VxDv64DgFsBPQyALThcAJSkhx9+WAUFBdqzZ491rLq6WsnJyfLy8tKIESPsWB0AXB89DEBbc3755ZdftncRN9uAAQO0Y8cOpaam6tKlSyosLNSrr76qo0eP6ve//70GDBhg7xJ/ENzc3HTPPffIzc3N3qXAxjj29kUPaz0+w+bFsW8aJ6P2J2cO5ty5c/rLX/6i9PR0VVdXy9fXV3FxcQoNDbV3aQBwQ/QwAG3JYQMgAAAAGuaQcwABAADQOAIgmiQ2NlZTp061dxmwkc8//1wBAQE6dOiQvUsBWo3+ZS70r6YhAOKGFixYoIyMDHuXARs5deqU4uPjVVNTY+9SgFajf5kL/avpCIBoVFlZmWbMmKG3337b3qXARtLT0/XII4/o22+/tXcpQKvQv8yH/tU8BEA06PDhwwoLC9NHH32k+Ph4e5cDG5g1a5bi4uLUtWtXPfzww/YuB2gx+pf50L+ajwCIBhUUFMjf319r165VQkKCvcuBDeTn52vGjBnasGGDevfube9ygBajf5kP/av52tu7ANyaIiIiFBUVZe8yYENr166Vq6urvcsAWo3+ZT70r+bjDCAaxD8k8+GYw1HwWTYfjnnzEQABAABMhkvAJnbp0iWVl5fXGXN2dpaXl5edKgKApqF/Aa1DADSxbdu2ac6cOXXGevbsqT179tipIgBoGvoX0DoEQBMLDg5WUlJSnTE3Nzc7VQMATUf/AlqHAGhi3bt3V/fu3e1dBgA0G/0LaB1+BAIAAGAyBEAAAACTcTIMw7B3EQAAALAdzgACAACYDAEQAADAZAiAAAAAJkMABAAAMBkCIAAAgMkQAAEAAEyGAAgAAGAyBEDcMp577jn5+/vrk08+kSQ98MAD8vf315UrV+xcmW1FR0fL399fhYWF9i4FQBPRv66if/1w8Cxg3LKmTJmi8vJytWtnru8pUVFRuvvuu3XbbbfZuxQALUT/on/d6giAuGVNnTrV3iXYxZgxY+xdAoBWon/hVmeuryYAAADgDCBsr7i4WG+++aZ2796t8+fPy9fXV9OnT6+33gMPPKCzZ88qOztb7dtf/agahqENGzZo48aNys3NVUVFhTw8PBQYGKipU6dq5MiRdV6jurpaK1eu1KZNm/TVV1+pa9euioqK0pAhQzR16lTFx8crISFB0tW5K8eOHVNGRoaWLFmitLQ0FRcXq2fPnoqKitK0adPk7Oxc5/Xz8/O1fPlyZWZmqqSkRJ6enrrvvvsUFxenvn371ln3xIkT+utf/6pjx47p22+/laenp4YMGaLY2FgFBARY14uOjtbBgwe1c+dO/eQnP5EklZaWasmSJcrMzNTZs2fVoUMHBQQEaPLkyQoLC2v9QQHQJPQv+pejIADCpoqLizV+/HgVFhZq8ODBCg8P1/Hjx5WQkKBu3brdcPvf/e53Wr9+vfr166fIyEi5uLgoLy9P+/fvV1ZWllasWKERI0ZIkr777jvFxsYqMzNTfn5+Gj9+vC5cuKDly5frH//4R4Ovf+XKFUVHR6ukpESjRo2Ss7OztmzZokWLFunixYuaMWOGdd2srCzFxsaqurpaISEh8vX1VX5+vjZv3qy0tDT9/e9/17333itJOnXqlKKjo3X58mWFh4erW7duKiwsVFpamtLT07Vu3TpZLJYGa6qurta0adN07NgxjRo1SqGhoSopKdGOHTsUHx+vP/7xjxo3blxzDwWAZqJ/0b8cigHY0Lx58wyLxWIsWrSozviqVasMi8ViWCwW4+OPPzYMwzBGjRplWCwW4/Lly4ZhGEZ2drZhsViMiRMnWseu3X7mzJnWsdWrVxsWi8V4+umnjf/+97/W8UOHDhkBAQGGxWIxlixZYh2fPHmyYbFYjEcffdSoqKiwjp84ccIICAgw7r77buPKlSuGYRhGRUWFMWzYMGPAgAFGRkZGnVrS0tIMf39/Y9iwYUZlZaVhGIbxpz/9ybBYLMaBAwfqrPv+++8bFovFeOmll+rVUVBQYBiGYaSnpxsWi8VYuHBhnW2//PJLo3///sZDDz3U2F83gJuI/kX/ciTMAYTNXL58WVu3bpWXl5fi4+PrLIuJiZG/v/91t//xj3+s1157TfPmzbNeUql13333SZLOnz9vHUtNTVW7du30wgsvyNXV1To+ePBgRUZGNvo+v/nNb/SjH/3I+v9+fn7q06ePSkpKVFxcLElKT0/X+fPnFRUVVe+yTWhoqCIiInT+/HmlpaVJunrpR5I+++wz65+lq7+Y2717t+bOndtoPTU1NZKk3NxcXbx40Tru6+urnTt36sMPP2x0WwA3B/2L/uVoCICwmVOnTqm8vFwDBw6s1wAlaejQodfdvkePHho9erQsFovy8/OVlpamVatWae7cuXr22WclXb1sIl297JCTkyNvb2/dcccdzXqva+e+SFKnTp2srytJ2dnZkmS9RHKte+65R5KUk5MjSRo3bpzc3d21ePFiDR8+XLNmzVJqaqpKSkrUq1evenNzvi84OFi+vr5KT0/XsGHDNG3aNCUlJenEiRPy8fFRhw4dGt0WwM1B/6J/ORrmAMJmSktLJf2/GV3L09Pzhq+xefNmLV261HqTURcXF/n5+WngwIHKz8+3rnfhwgUZhtHovJwePXo0+h5ubm71xpycnCT9/5tweXm5JMnDw+O6r1/7jdfPz0/r16/X22+/rT179mjLli3asmWLnJycNHz4cL388svy8fFp8LVcXV2VkpKilStXavv27dq3b5/27dsnSerXr59efPFFDRs2rNH9AdB69C/6l6MhAMJmahtkbSO9VlVV1XW337Vrl55//nndcccdWrBggQYOHCgfHx/rROpNmzZZ1629BFLb6K5VWVnZkl2wqm2cRUVFDS6v3ccuXbpYx3x9fZWYmKiamhodP35cWVlZ2rp1q/bv36+nnnqq0YndknTbbbdp5syZmjlzpk6fPq1PPvlEu3btUnp6up588kmlpaXJ29u7VfsEoHH0L/qXo+ESMGzmzjvvlKenp44ePWq9FPF9hw8fvu72GzZskCT9+c9/VmRkpPr27SsXFxdJUl5enqT/f8P18PCQr6+vCgsLrfNemvNeNxIYGChJ1sc+Xevjjz+WJOu8oJSUFL3yyisyDEPt2rVTYGCgHn/8ca1du1a9evVSXl5eg3VK0j//+U/94Q9/sJ418PHx0bhx47Rs2TJFRkaqurq61fsD4ProX/QvR0MAhM20b99eUVFRKi0t1euvv26dHCxJmzZtumETqJ0r8vXXX9cZP336tBYuXChJdZ67+eijj+q7775TYmKiLl++bB3Pzs7WBx980Kp9CQ0NVZcuXbR161alp6fXWbZ3715t2rRJXl5euv/++yVJBw8e1OrVq+t8y5euftMuKytTp06d1Llz5wbf68yZM0pOTtby5cvrjNfU1Oirr76SJPXq1atV+wPg+uhf9C9HwyVg2NQzzzyjgwcPKjk5Wf/+9781dOhQnTx5Unv37lXv3r1VUFDQ6LZjxozR1q1b9cILL2jv3r3y9vbW6dOnlZGRITc3N7m4uOjChQvW9SdPnqxdu3Zp8+bNOn78uO69914VFxdr586d8vDwUFVVVYuf0+nu7q4FCxYoPj5e06dP14gRI6z30dq3b586duyohQsXyt3dXZIUHx+v/fv3a86cOdq+fbv69eunsrIy7dq1S2VlZZo/f36DE8slafTo0Vq3bp1SU1OVm5uroUOHqqamRllZWcrLy1NERISCgoJatB8Amo7+Rf9yJJwBhE25u7srOTlZsbGxKi0t1XvvvaczZ87o1Vdf1S9/+cvrbjt8+HAtW7ZMgYGBysjI0Hvvvacvv/xSY8aM0ebNm3XXXXepoKBAJ0+elHT1G/tbb72luLg4VVVVKSUlRYcPH9bTTz+t2NhYaz0tFRISovXr1ysiIkI5OTlKTk5WXl6exo0bp40bN9aZ2Ozr66vVq1crPDxc2dnZSkpK0o4dOxQQEKC33npL48ePb/R9OnbsqJUrVyo6OlplZWVavXq11q1bJzc3N82bN0+vv/56i/cBQNPRv+hfjsTJ+P5NfQAHcubMGXXp0qXOPbFqLVq0SMuWLdMbb7yhhx56yA7VAUDj6F9oa5wBhMNKTEzUz3/+c/3rX/+qM15UVKTU1FS5uLjc8N5dAGAP9C+0Nc4AwmFlZmbqiSeekJubm0JDQ3X77bfr3Llz2r17t0pLS/Xiiy9qypQp9i4TAOqhf6GtEQDh0I4cOaJ33nlHn3/+uc6dO6fOnTsrKChI0dHRCg4Otnd5ANAo+hfaEgEQAADAZJgDCAAAYDIEQAAAAJMhAAIAAJgMARAAAMBkCIAAAAAmQwAEAAAwmf8BKqkJCyHqUzkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5])\n", "sns.countplot(x='diagnosis', data=y, palette='husl', ax=ax1)\n", "sns.countplot(x='diagnosis', data=y_balanced, palette='husl', ax=ax2)\n", "ax1.set_title(\"Original\", weight='bold')\n", "ax2.set_title(\"Balanced (downsampling)\", weight='bold')\n", "ax1.set_ylabel(\"No. of patients\", weight='bold')\n", "ax2.set_ylabel(\"\")\n", "fig.savefig(PATH_RESULTS/('class_proportion.pdf'), dpi=1000, bbox_inches='tight')\n", "plt.show();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data partition\n", "Partition data into training and testing sets." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "X1_train, X1_test, y_train, y_test = train_test_split(X1_balanced, y_balanced, test_size=TEST_PROP, random_state=SEED)\n", "X2_train, X2_test, _, _ = train_test_split(X2_balanced, y_balanced, test_size=TEST_PROP, random_state=SEED)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## VERTIGO\n", "### 0. Preparations\n", "Before actually computing VERTIGO, we need to do a few things. These steps are shown here for didactic purposes. **Notice that the object `vertigo.vertigo` already does them internally**. \n", "\n", "First, we need to convert DataFrames to numpy arrays (for mathematical operations). Furthermore, we will rename the variables to match the paper's notation." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "X1 = X1_train.to_numpy()\n", "X2 = X2_train.to_numpy()\n", "\n", "Y = y_train.to_numpy()\n", "Y_diag = np.diag(Y.flatten())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add column of 1s to serve as constant term in the computation of the residual coefficient." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "X1 = sm.tools.add_constant(X1, prepend=False)\n", "X2 = sm.tools.add_constant(X2, prepend=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute some variables for ease of use." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "M = X1.shape[0] # Number of patients (same for both parties).\n", "N1 = X1.shape[1] # Number of party 1's features.\n", "N2 = X2.shape[1] # Number of party 2's features." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initialize alpha." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "ALPHA_INIT = ALPHA_INIT_VAL * np.ones((M, 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lastly, we will define a few useful functions. These are also part of `vertigo.vertigo` methods" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "# Check if matrix is diagonally dominant.\n", "def is_diagonally_dominant(X):\n", " D = np.diag(np.abs(X)) # Find diagonal coefficients\n", " S = np.sum(np.abs(X), axis=1) - D # Find row sum without diagonal\n", " if np.all(D > S):\n", " return True\n", " else:\n", " return False\n", "\n", "# Check if matrix is positive definite.\n", "def is_positive_definite(x):\n", " return np.all(np.linalg.eigvals(x) > 0)\n", "\n", "# Sigmoid function to bound a number between 0 and 1\n", "def sigmoid(x):\n", " return 1 / (1 + np.exp(-x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Local\n", "\n", "#### 1.1 Compute local gram matrix (i.e., linear kernel matrix)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "K1 = np.matmul(X1, X1.T)\n", "K2 = np.matmul(X2, X2.T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Server\n", "#### 2.1 Compute global gram matrix $K$" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "K = np.add(K1, K2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.2 Compute fixed Hessian matrix $\\tilde{H}$(Eq. 10)\n", "Originally, the paper states that $c$ must be set to a value that makes $\\tilde{H}$ diagonally dominant. We initialize it as \"the maximum of the elements in the original $H$\" (p. 573). Then, we increase $c$ by a factor in steps of 0.5 until the condition of diagonally dominance is satisfied.\n", "\n", "However, in our experiments we didn't find this condition to be necessary, getting the same results if we adjusted $c$ or not. We give the option to the user to do so (see line 14). \n", "\n", "Depending on the data, computing $\\tilde{H}$ can take quite long." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "not_diagonally_dominant = True\n", "c_factor = 1\n", "\n", "while not_diagonally_dominant:\n", "\n", " H = (1/LAMBDA_) * np.matmul(np.matmul(Y_diag, K), Y_diag) + np.diag(1/(ALPHA_INIT*(1-ALPHA_INIT)))\n", " c = c_factor * H.max()\n", " I = np.identity(len(Y_diag))\n", "\n", " H_wave = (1/LAMBDA_) * np.matmul(np.matmul(Y_diag, K), Y_diag) + c*I\n", " not_diagonally_dominant = not is_diagonally_dominant(H_wave)\n", " c_factor += 0.5\n", "\n", " # Comment break to make sure that H_wave is diagonally dominant.\n", " break" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2.3 Compute $\\tilde{H}^{-1}$" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "H_wave_1 = np.linalg.inv(H_wave)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. Local + server\n", "Compute $\\alpha*$" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "# We will store the evolution of the alpha coefficients here.\n", "alpha_all = np.zeros((Y.shape[0], ITERATIONS))\n", "\n", "for s in range(0, ITERATIONS):\n", "\n", " if s==0:\n", " alpha_old = ALPHA_INIT \n", " else:\n", " alpha_old = alpha\n", "\n", " alpha_all[:, s] = alpha_old.T\n", " \n", "\n", " # 3.1 LOCAL\n", " # Eq. 5\n", " beta1 = (1/LAMBDA_) * np.sum((npm.repmat(alpha_old, 1, N1) * \n", " npm.repmat(Y, 1, N1) *\n", " X1), 0)\n", " beta2 = (1/LAMBDA_) * np.sum((npm.repmat(alpha_old, 1, N2) * \n", " npm.repmat(Y, 1, N2) *\n", " X2), 0)\n", "\n", " # In order for the matrix multiplications to work, Y needs to be\n", " # transposed (which is not stated in the paper).\n", " E1 = Y.T * (np.matmul(beta1, X1.T))\n", " E2 = Y.T * (np.matmul(beta2, X2.T))\n", "\n", "\n", " # 3.2 SERVER\n", " E = E1 + E2\n", "\n", " # Catch incorrect cases (logarithm of a negative number)\n", " with warnings.catch_warnings():\n", " warnings.filterwarnings('error')\n", " try:\n", " J = E.T + np.log10(alpha_old/(1-alpha_old))\n", " except Warning as e:\n", " print('WARNING:', e)\n", "\n", " \n", " # 3.3 SERVER\n", " # Update alpha.\n", " alpha = alpha_old - np.matmul(H_wave_1, J)\n", "\n", " # Bound alpha (0 < alpha < 1)\n", " # This isn't defined in the original paper, but is done to avoid\n", " # problems when computing J (the logarithm of a negative number is\n", " # undefined).\n", " alpha = sigmoid(alpha)\n", "\n", " # Check for convergence.\n", " if max(abs(alpha - alpha_old)) < TOLERANCE:\n", " # Trim alpha_all.\n", " alpha_all = np.delete(alpha_all, np.s_[s:], axis=1)\n", " break\n", "\n", "# Final alpha value.\n", "alpha_star = alpha" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. Local\n", "Compute the $\\beta$ coefficients (and the intercept). This step isn't defined in Fig. 4, but it is the logical next step" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Party 1 coefficients are [0.02656502 0.01223053 0.02735187 0.02322579 0.01138798]\n", "Party 2 coefficients are [0.02350349 0.03031696 0.03433975 0.00996463 0.00073673]\n", "Intercept is -0.002826607481672041\n" ] } ], "source": [ "coeffs1 = (1/LAMBDA_) * np.sum((npm.repmat(alpha_star, 1, N1) * \n", " npm.repmat(Y, 1, N1) *\n", " X1[:, :]), 0)\n", "coeffs2 = (1/LAMBDA_) * np.sum((npm.repmat(alpha_star, 1, N2) * \n", " npm.repmat(Y, 1, N2) *\n", " X2[:, :]), 0)\n", "\n", "beta1 = coeffs1[0:-1]\n", "beta2 = coeffs2[0:-1]\n", "intercept = coeffs1[-1] # Notice how the intercept is the same for both parties.\n", "\n", "print(\"Party 1 coefficients are \", end='', flush=True)\n", "print(beta1)\n", "print(\"Party 2 coefficients are \", end='', flush=True)\n", "print(beta2)\n", "print(\"Intercept is \", end='', flush=True)\n", "print(intercept)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And that's it! These are the coefficients of our model. They have the same order as the input features. They can be found in a `vertigo.vertigo` object as `beta_coef_` and `intercept_`." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (Spyder)", "language": "python3", "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }