{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Secure Kaplan-Meier Survival Analysis Explained\n", "\n", "The MPyC demo [kmsurival.py](kmsurvival.py) implements privacy-preserving Kaplan-Meier [survival analysis](https://en.wikipedia.org/wiki/Survival_analysis), based on earlier work by Meilof Veeningen. The demo is built using the Python package [lifelines](https://pypi.org/project/lifelines/), which provides extensive support for survival\n", "analysis and includes several datasets. We use lifelines for plotting Kaplan-Meier survival curves, performing logrank tests to compare survival curves, and printing survival tables." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# use pip (or, conda) to make sure lifelines package is installed:\n", "!pip -q install lifelines " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import os, functools\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import lifelines.statistics\n", "from mpyc.runtime import mpc\n", "mpc.logging(False)\n", "from lifelines import KaplanMeierFitter\n", "from kmsurvival import plot_fits, events_to_table, events_from_table, logrank_test, aggregate, agg_logrank_test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Actual use of the lifelines package is hidden mostly inside the [kmsurival.py](kmsurvival.py) demo, except for the function `lifelines.statistics.logrank_test()` used below to validate the results of the secure computations.\n", "\n", "## Kaplan-Meier Analysis\n", "\n", "We analyze the aml (\"acute myelogenous leukemia\") dataset, which is also used as a [small example on Wikipedia](https://en.wikipedia.org/wiki/Survival_analysis#Example:_Acute_myelogenous_leukemia_survival_data). The file `aml.csv` (copied from the R package `KMsurv`) contains the raw data for 23 patients. Status 1 stands for the event \"recurrence of aml cancer\" and status 0 means no event (\"censored\")." ] }, { "cell_type": "code", "execution_count": 3, "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", " \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", "
observationtimestatusgroup
12512
13512
14812
15812
1911
161212
21311
31301
171602
41811
52311
182312
192712
62801
203012
73111
213312
83411
224312
94501
234512
104811
1116101
\n" ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv(os.path.join('data', 'surv', 'aml.csv')).rename(columns={'Unnamed: 0': 'observation', 'cens': 'status'})\n", "df.sort_values(['time', 'observation']).style.hide(axis='index')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Time is in weeks. The study compares the time until recurrence among two groups of patients. Patients in group 1 received maintenance chemotherapy, while patients in group 2 did not get any maintenance treatment. To plot the [Kaplan–Meier survival curve](https://en.wikipedia.org/wiki/Kaplan–Meier_estimator) for groups 1 and 2, we use the function `plot_fits()` as follows." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD+CAYAAABvEpGeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9eUlEQVR4nO2deZhUxdX/P18WQQFBBRQFBBdQEFlEwY2AKO6SV3FB8wqJSzTuedWfu2g0ajBGjUbjSlQ07sa4L4GouIIiiLigooKIgLIpIMv5/VHVcOnpmekeZqbvzJzP8/Qzt+vWrfu9Nd19btU9dY7MDMdxHMdJG/WKLcBxHMdxcuEGynEcx0klbqAcx3GcVOIGynEcx0klbqAcx3GcVOIGynEcx0klbqCcOoOk/pJmFFtHWUh6VtKwYusoFEmLJW1VCe2YpG0qQ5NT83ED5aQOSdMl7Z14f5SkHyT9opi6ykLS8Pjj+pes8sGxfFQ+7ZjZ/mb2jyoRWYWYWVMz+7zYOpzahRsoJ9XE0cTNwIFm9t9i6ymHz4AjJDVIlA0DPqnqE2eds8a0XQwk1S+2Bic/3EA5qUXSb4E/A/ua2euxbGtJ/5E0T9JcSaMltUgcM13S+ZI+jKOuuyU1LqX98yR9JmlRrP8/iX3DJb0m6drYzheS9i9H8rfAZGDf2MbGwG7Ak1nn7SvpdUnzJb0vqX9i31hJxyfe/0bS1KjheUlbJvaZpFMkfQp8muP6Gku6L/bVfEnvSNo00U/JUeoISffF7Q6x7eMkfQX8J049nprV/vuSDk1o2UZSH0nfJo2ApP+RNClu7yLpjahnlqSbJK1XTr9m2tk4/j+/if3xRCwfLum1rLqrpwoljZJ0i6RnJP0InF2OxnqJz8Y8SQ/F/2WZfepUPm6gnLRyMnA5MNDMxifKBVwFbA5sD7QDRmQdewzBSGwNdAIuKuUcnwF7As2By4D7JLVJ7O8DfAy0BP4E3ClJ5ei+Bzg2bh8F/AtYtlq8tAXwNHAFsDFwNvCopFbZDUkaDFwAHAq0Al4FHsiq9suos0sOLcPitbUDNgFOApaUoz/JLwh9vG8879CEti7AlvFaVmNmbwE/Anslio8G7o/bK4GzCH26KzAQ+F2eeu4FNgC6Aq2Bv5RdfS2OBq4EmgE3lKPxNEK//oLwOfuBMIqHde9TpwDcQDlpZR/gTcKIZDVmNs3MXjSzZWY2B7iO8EOS5CYz+9rMvif8KA0lB2b2sJl9Y2arzOxBwihkl0SVL83sdjNbCfwDaAOUd7f8ONBfUnOCobona/+vgGfM7Jl43heB8cABOdo6CbjKzKaa2Qrgj0CP5Cgq7v/ezHL9SC4n/IhuY2YrzWyCmS0sR3+SEWb2Y2z78axzHwM8ZmbLchy32phJahav7QGAqOFNM1thZtOBv1Py/1eCeOOwP3CSmf1gZssLnPL9l5mNi32+tCyNhH6/0MxmxOsbAQyJU53r2qdOAbiBctLKyYTRzx3JUYukTSX9U9JMSQuB+wh340m+Tmx/SbgLLoGkYyVNjFM184Edstr6NrNhZj/FzaaS9lTwWlssaUqyzfhj/jRh1LaJmY3LOu2WwOGZc8bz7kEwftlsCdyQqPc9YQS5RSnXms29wPPAP+O02J8kNSyjfjar2zazRfG6jopFQ4HRpRx3P3CopEaE0d+7ZvYlgKROkp6KU2wLCUY3+/+Xi3bA92b2QwH6k2T3U6kaCf3+eKLfpxJGfpuy7n3qFIAbKCetzCZM/+wJ/C1R/kfAgG5mtiFhRJI97dYusd0e+Ca78TgSuB04lWBIWgAf5GirBGb2avRaa2pmXXNUuQf4P4LxzOZr4F4za5F4NTGzq0up+9usuutnnsdl5JShc7mZXWZmXQjPwg5izfTjj4Tpsgyb5Woi6/0DwFBJuwKNgTGlnPdDwo3B/qw9dQZwC/ARsG38/11AHn1O6IuNlXjemGCta5FU7rWUo/FrYP+sfm9sZjPL6VOnknED5aQWM/uGYKT20xr37WbAYmBBfJ5zTo5DT5HUNj7YvhB4MEedJoQfrTkAkn5NGEFVBv8lTFH+Nce++4CDJe0rqX586N5fUtscdW8FzpfUNWpsLunwfEVIGiCpW3QGWEiYnloVd08EjpLUUFJvYEgeTT5DGF1cDjxoZqvKqHs/cAbQD3g4Ud4salksaTvCSLlczGwW8CzwN0kbRd394u73ga6Seig4xIzIp80yNN4KXJmZzpTUKj4PLK9PnUrGDZSTaszsK8LD7CGSriI4M/QCFhCmnB7Lcdj9wAvA5wRHiCtytPshwUPwDcJorRuQPR1XUc1mZi/HZ2DZ+74GMs4Pcwh36+eQ47toZo8D1xCmkxYSRnjleRIm2Qx4hPBDOpVgOO+N+y4mOJH8QOjT+3M1kKVnGaG/986j/gOEZ0v/MbO5ifKzCSOWRYQRbK6bh9L4X4JB+Aj4Djgz6vqEYDRfIjxHfK2U4/PVeAPB8/IFSYsIz0L7xH1l9alTycgTFjq1CUnTgePN7KVia3EcZ93wEZTjOI6TStxAOY7jOKnEp/gcx3GcVOIjKMdxHCeVuIFyHMdxUkmtilJcCJLuIiyy+87MSqx/idELbiCEQPkJGG5m75bXbsuWLa1Dhw6VrNZxHKf2MmHChLlmViIeZZ01UMAo4CZKxkrLsD+wbXz1IayA71NK3dV06NCB8ePHl1fNcRzHiUj6Mld5nTVQZvaKpA5lVBkM3GPBi+RNSS0ktYkr2iudN/92As3mT82r7rj1B/DyBrlii8LgHltwdJ/2lSnNcRynKNRZA5UHW7B2gMkZsayEgZJ0InAiQPv2FTcOq1asKLdOx1VfsmrFCh6vt0+JfV/O+4mFS5e7gXIcp1bgBqoSMLPbgNsAevfuXSG//b6/uz2/incfSLel83n25H4ldh359zf44aefeXLizFIPb7Z+QwZ0bl0RiY7jONWKG6jSmcnaUbHbxrJUs2qV0apZzgSyAMxZtLQa1ThOSZYvX86MGTNYutQ/i3WNxo0b07ZtWxo2zC9DiRuo0nkSOFXSPwnOEQuq6vmT49QlZsyYQbNmzejQoQPlJyh2agtmxrx585gxYwYdO3bM65g6a6AkPQD0B1pKmgFcCjQEMLNbCakFDgCmEdzMf10cpTmwlTD5kbXLGjcHmpZ76NIVK8ucAgSfBnSqlqVLl7pxqoNIYpNNNmHOnDl5H1NnDZSZ5UwDnthvwCnVJKcwVq2CplmZxxfPJh8D1W6jJuXW8WlAp6px41Q3KfT/7pEkHMepc/zmN7+hdevW7LBD/jkqx44diyTuuOOO1WUTJ05EEtdee22Zx956663cc09pSy7XtPXMM8+Uq2P8+PGcfvrp+Ykuh1GjRnHqqadWSltVQZ0dQdVWZs5fwuVPTSl1/+5bt2Tg9puWut9x6gLDhw/n1FNP5dhjC8vWvsMOO/DQQw9x/PHHA/DAAw/QvXv3co876aSTyq0zceJExo8fzwEH5F7jmKF379707t07P8E1HDdQtYjBPbZg3uJlpe7/ct5PwNxyDVQ+z6nAn1U5NZd+/foxffr0go/bcsstWbhwIbNnz6Z169Y899xzaxmU22+/ndtuu42ff/6ZbbbZhnvvvZcNNtiAESNG0LRpU84++2z69+9Pnz59GDNmDPPnz+fOO++kT58+XHLJJSxZsoTXXnuN888/n44dO3LGGWewdOlS1l9/fe6++246d+7M2LFjufbaa3nqqacYMWIEX331FZ9//jlfffUVZ5555urR1X333ceNN97Izz//TJ8+ffjb3/5G/fr1ufvuu7nqqqto0aIF3bt3p1GjRpXVrZWOG6haxNF92tO0Uf1S3cwvf2oKcxaVbsAy5POcCvxZlbPuXPbvKXz4zcJKbbPL5hty6cFdCz5u5MiRjB49ukR5v379uPHGG1e/HzJkCA8//DA9e/akV69ea/3AH3rooZxwwgkAXHTRRdx5552cdtppJdpcsWIFb7/9Ns888wyXXXYZL730Epdffjnjx4/npptuAmDhwoW8+uqrNGjQgJdeeokLLriARx99tERbH330EWPGjGHRokV07tyZk08+mWnTpvHggw8ybtw4GjZsyO9+9ztGjx7NPvvsw6WXXsqECRNo3rw5AwYMoGfPngX3VXXhBqqOMXfxz8WW4Dip5JxzzuGcc84pt94RRxzBkUceyUcffcTQoUN5/fXXV+/74IMPuOiii5g/fz6LFy9m3333zdnGoYceCsBOO+1U6khuwYIFDBs2jE8//RRJLF++PGe9Aw88kEaNGtGoUSNat27N7Nmzefnll5kwYQI777wzAEuWLKF169a89dZb9O/fn1atQlzWI488kk8++aTcay4WbqBqC8uXwuRH2PzrH9hw/ZKL4FastyHQsvp1OU4ZVGSkU1XkO4LabLPNaNiwIS+++CI33HDDWgZq+PDhPPHEE3Tv3p1Ro0YxduzYnOfKjLrq16/PilJCnF188cUMGDCAxx9/nOnTp9O/f/8y20q2Z2YMGzaMq666aq26TzzxRM420oobqJrIgq/gufPWLtuqP3Taj2WN67N8/fVKHNJwyRzcQDlO6eQ7ggK4/PLL+e6776hfv/5a5YsWLaJNmzYsX76c0aNHs8UWW+R9/mbNmrFo0aLV7xcsWLD6+FGjRuXdDsDAgQMZPHgwZ511Fq1bt+b7779n0aJF9OnThzPOOIN58+ax4YYb8vDDD+fl5FEs3EDVNLoNgR+zFrp9/0X422m/apWSrzNFLtzBwikmQ4cOZezYscydO5e2bdty2WWXcdxxx+V9/G677Zaz/A9/+AN9+vShVatW9OnTZy2DUx4DBgzg6quvpkePHpx//vmce+65DBs2jCuuuIIDDzww73YAunTpwhVXXMGgQYNYtWoVDRs25Oabb6Zv376MGDGCXXfdlRYtWtCjR4+C2q1uFNajOpVF7969rcrzQU1+ZO2FupnR1H5XM27aXFpskHsEdcaUrZg6axHbt2lWZvPV4Yo+Z9FSDumR/92lU3uYOnUq22+/fbFlOEUi1/9f0gQzK+E77yOoWsroj+H+tZ59tgLC3dzUWSXv6lo2XY9WzRrl7YruOI5T1biBqqUc0zm8MjRa/DWr6jdm0EuteGHv0mNhnT2+OSxZwiZf/Lvcc6xYb0MWbPGLypDrOI5TAjdQdYRlTddkDlm+fqtS61n98utkCI4XjuM4VYPH4nMcx3FSiRsox3EcJ5X4FJ9TFHK5qLvrueM4SXwEVcc4ulOxFQTabdSEVs0ar/VatCR3KBfHqUy+/vprBgwYQJcuXejatSs33HBD3sdmUm78+99rnIgOOuigUiNGZLj++uv56aefKiq5TEpbk1WR8x9//PF8+OGHlSGLDh06MHfu3HVqw0dQdYykZ19pfL4Aznu9/HoDW9bnAHJ7+7mHn5NWGjRowJ///Gd69erFokWL2Gmnndhnn33o0qVLXse3bduWK6+8koMPPjjvc15//fX86le/YoMNNqio7FJJhlpa1/Mnc12lAR9BOWvxiy1gq+bl1/t8Abw8d2OWr98q56vBz5UbodpxKos2bdrQq1cvIIQX2n777Zk5M/+IKN27d6d58+a8+OKLJfa9/PLL9OzZk27duvGb3/yGZcuWceONN/LNN98wYMAABgwYUOKYDh06cP7559OjRw969+7Nu+++y7777svWW2/NrbfeCsDixYsZOHAgvXr1olu3bvzrX/9afXzTpiGT9tixY+nfvz9Dhgxhu+2245hjjsHMcp7/5JNPpnfv3nTt2pVLL710dVv9+/cnE2igadOmXHjhhXTv3p2+ffsye/ZsAObMmcNhhx3GzjvvzM4778y4ceMAmDdvHoMGDaJr164cf/zxVEYQCB9BOWux/5bhVR75jLAcp1yePQ++nVy5bW7WDfa/Oq+q06dP57333qNPnz5A/gFjL7zwQi6++GL22Wef1WVLly5l+PDhvPzyy3Tq1Iljjz2WW265hTPPPJPrrruOMWPG0LJl7niY7du3Z+LEiZx11lkMHz6ccePGsXTpUnbYYQdOOukkGjduzOOPP86GG27I3Llz6du3L4ccckiJFOrvvfceU6ZMYfPNN2f33Xdn3LhxnH766SXOf+WVV7LxxhuzcuVKBg4cyKRJk9hxxx3XauvHH3+kb9++XHnllZx77rncfvvtXHTRRZxxxhmcddZZ7LHHHnz11Vfsu+++TJ06lcsuu4w99tiDSy65hKeffpo777wzr/9BWbiBqi18/wU8dx47LFlOg/olB8YLNtuNH9oOrHZZj0z4miE7tSu/ouNUM4sXL+awww7j+uuvZ8MNNwTyDxjbr18/AF577bXVZR9//DEdO3akU6fwoHfYsGHcfPPNnHnmmeW2d8ghhwDQrVs3Fi9eTLNmzWjWrBmNGjVi/vz5NGnShAsuuIBXXnmFevXqMXPmTGbPns1mm222Vju77LILbdu2BaBHjx5Mnz6dPfbYo8T5HnroIW677TZWrFjBrFmz+PDDD0sYqPXWW4+DDjoICGlBMiPGl156aa3nVAsXLmTx4sW88sorPPbYY0BIAbLRRhuVe93l4QaqNrBV/zJ3N170JUBRDNSj7850A+WUTp4jncpm+fLlHHbYYRxzzDGrczNB/iMoCKOoK664ggYN1v1nNJMyo169emulz6hXrx4rVqxg9OjRzJkzhwkTJtCwYUM6dOjA0qUlE4bmSr2RzRdffMG1117LO++8w0YbbcTw4cNzttWwYcPVI7RkW6tWreLNN9+kcePciVErEzdQtYFO+62OZP5BjmCxHcb/oczDW332CHO2HlJl8vJlXaKjF4q7tNddzIzjjjuO7bffnt///vdr7Ssk5cagQYO4+OKLmTVrFgCdO3dm+vTpTJs2bXW691/8IjgKZVJplDbFVx4LFiygdevWNGzYkDFjxvDll18WdHzy/AsXLqRJkyY0b96c2bNn8+yzz5aaayoXgwYN4q9//evqfpo4cSI9evSgX79+3H///Vx00UU8++yz/PDDDwVpzIUbKIfWnz9W6Qaq3sqlMZ5fq7Xi+pXl3ZdvqvnKwNPV113GjRvHvffeS7du3Vanm/jjH//IAQccUHBbF154IYMHDwagcePG3H333Rx++OGsWLGCnXfemZNOOgmAE088kf3224/NN9+cMWPGFHyeY445hoMPPphu3brRu3dvtttuu4KOzz5/z5492W677WjXrh277757QW3deOONnHLKKey4446sWLGCfv36ceutt3LppZcydOhQunbtym677Ub79u0LajcXnm6jkilKuo0EudJtZEZQ03tfnPOYri8ezZR97i9IQsZJ4upylmAc+G94OuGN23DJHOZ1zN89t6rwdB/Fw9Nt1G083YZTLeS7XipZRyubs3zKlGrJOeU4Ts3GDZRTIX6RNfiY/RN8tyR33cnzku/Wg/mLmLNomRsox3HKxA1UHaHxoi/LdJYoz5Eim5OBo9uX77qea4rvjClbFXQux3HqJm6g6gALNgsPihoumcN6S3PHxmryw9QSZT83bllqXqhiuq47NR8zK7HI1Kn9FOrz4AaqDvBD24FlGpKKOEkUOuJynAyNGzdm3rx5bLLJJm6k6hBmxrx58wpaP+UGyqlW6q1cSsMl4aFUxv3cA8vWLdq2bcuMGTOYM8czMtc1GjduvDrSRT7UaQMlaT/gBqA+cIeZXZ21fzgwEsisHr3JzNIV7reGsaxpuxJp5T11fN2iYcOGdOzYsdgynBpAnTVQkuoDNwP7ADOAdyQ9aWbZyVAeNLNTq11gLSEt+accx6l51FkDBewCTDOzzwEk/RMYDFROti4HyC//VDGozrBKjlPbqarQYXXZQG0BfJ14PwPok6PeYZL6AZ8AZ5nZ19kVJJ0InAhUSniPdWGDRg2Y/9PPpe7/eeUqWjdb+yHld1sdWkrtqiO5yLcYi3erM6yS49R2qip0WF02UPnwb+ABM1sm6bfAP4C9siuZ2W3AbRBCHVWvxLXp2a5FmfvHTSvpZl7ROHxlra0qK71H9iJfgC/n/QTM9cW7juOspi4bqJlAMg9EW9Y4QwBgZskYCHcAf6oGXeXTuDksnr122fKlsFEemQYriczaqlyUt0YqOyliwyULOGPKJpWqz3Gcmk9dNlDvANtK6kgwTEcBRycrSGpjZrPi20OAkqtZi8G2+5Qsm/xItUooa21VoWukcrmeVxR3WXec2kOdNVBmtkLSqcDzBDfzu8xsiqTLgfFm9iRwuqRDgBXA98DwogmuxeRyPa8o7rLuOLWHOmugAMzsGeCZrLJLEtvnA+dXt66qpDwnCsjtSFFT8ZTzjlNzqdMGqi5SnhMF5HakqKl4ynnHqbnUK7YAx3Ecx8mFj6CcKqE0F/Sy3M8dx3GSuIFyKp3SXNA9RYfjOIXgBsqpdEpzQS/P/TzfFPJlkYlKkeTyp6aUqOcp5x0n/biBclJBrugShbAm5XxIKZ9k6qxFJep7ynnHST9uoJwS5OOKDpXrjp4dXaKiNFr8Navqr9E06KVWPHBC37Xq5BpROY6TPtxA1RZyhT9KUkAopHxc0SGd7ujLmrpLuePUFtxA1RZyhT9KUs2hkBzHcdYVXwflOI7jpBIfQTnVSq71Ub42ynGcXLiBcqqNXOujqnpt1K86/lgl7TqOU/W4gaorlOdEUQE2+HEWbLBV3vVzrY8qNDVHoRy79U/MK7+a4zgpxA1UXaE8J4oKsN43d+Xljp6L2hQx3XHqOs9MnsUhPdZxMWMO3EA5Fabr5htC05YVOjaNLuqO41SM56ZU7uxMBvficxzHcVKJGyjHcRwnlfgUn1N0SkvNURHcZd1xag9uoJyKU5ZnYJ6hlUpLzVEhOQW4rH857yePyec4KccNlFNxyvIMzDO0UmmpOSpCvqOw3bduCbiThuMUypxFy5i7OLfnbofzni5RdsbAbTlrn04VPp8bKKco5IqYXl2u5wO339RTbThOJTL09jeZfvWBld6uGyinKOSKmO6u547jJHEvPsdxHCeVuIFyHMdxUolP8TlVQ2kefgUkTnQcp27jBsqpGkrz8KvixInZa6oWb7Ij8zoeXKXndJy6zn5dq8bpyA2UkxpyefYVwqyN+9Bq5SpYuQqAJj9+xarly5izaGmJuktXrKTdRk0qfC7HcdZwQLc2VdKuGygnNeTy7CuIbYYAQ9a8f+48mq1cnjPK8pMTZ67buRzHqXLcQDnVS65nU/5cynGcHLiBcqqXXM+mqvi5lOM4NZM67WYuaT9JH0uaJum8HPsbSXow7n9LUociyHQcx6mT1NkRlKT6wM3APsAM4B1JT5rZh4lqxwE/mNk2ko4CrgGOrH61tZzKTEfv04WOU2uoswYK2AWYZmafA0j6JzAYSBqowcCIuP0IcJMkmZlVp9BaT2Wmo/fpQsepNdRlA7UF8HXi/QygT2l1zGyFpAXAJmSFwpZ0InAiQPv27atKr5MPydFYs82gfqOc1Zqt3zCn+7njOIXTbP2GVdJuXTZQlYaZ3QbcBtC7d28fXRWT5Gis25BSqw3o3LoaxDiOsy7UZSeJmUC7xPu2sSxnHUkNgObAvGpR5ziOU8epywbqHWBbSR0lrQccBTyZVedJYFjcHgL8x58/OY7jVA+qy7+3kg4ArgfqA3eZ2ZWSLgfGm9mTkhoD9wI9ge+BozJOFWW0OQf4soKS0prq1XXlTxo1gesqhDRqgtqta0sza5VdWKcNVNqQNN7MehdbRzauK3/SqAlcVyGkURPUTV11eYrPcRzHSTFuoBzHcZxU4gYqXdxWbAGl4LryJ42awHUVQho1QR3U5c+gHMdxnFTiIyjHcRwnlbiBchzHcVKJG6iUUF7qj2rS0E7SGEkfSpoi6YxYvrGkFyV9Gv9uVCR99SW9J+mp+L5jTIMyLaZFWa8ImlpIekTSR5KmStq12P0l6az4//tA0gOSGhejryTdJek7SR8kynL2jQI3Rn2TJPWqZl0j4/9wkqTHJbVI7Ds/6vpY0r7VqSux7/8kmaSW8X219FdpmiSdFvtriqQ/Jcort6/MzF9FfhEWCn8GbAWsB7wPdCmCjjZAr7jdDPgE6AL8CTgvlp8HXFOkfvo9cD/wVHz/EGHxNMCtwMlF0PQP4Pi4vR7Qopj9RQhw/AWwfqKPhhejr4B+QC/gg0RZzr4BDgCeBQT0Bd6qZl2DgAZx+5qEri7x+9gI6Bi/p/WrS1csbwc8TwgA0LI6+6uUvhoAvAQ0iu9bV1Vf+QgqHaxO/WFmPwOZ1B/VipnNMrN34/YiYCrhB28w4YeY+PeX1a1NUlvgQOCO+F7AXoQ0KEXRJak54Qt8J4CZ/Wxm8yl+fzUA1o/xIzcAZlGEvjKzVwgRWJKU1jeDgXss8CbQQlKb6tJlZi+Y2Yr49k1CbM6Mrn+a2TIz+wKYRvi+VouuyF+Ac4GkR1u19Fcpmk4GrjazZbHOdwlNldpXbqDSQa7UH1sUSQsACtmDewJvAZua2ay461tg0yJIup7wJV0V328CzE/8qBSjzzoCc4C749TjHZKaUMT+MrOZwLXAVwTDtACYQPH7KkNpfZOm78BvCKMTKLIuSYOBmWb2ftauYurqBOwZp4z/K2nnqtLkBsopgaSmwKPAmWa2MLnPwli+WtcmSDoI+M7MJlTnefOgAWH64xYz6wn8SJi2Wk1191d8pjOYYDw3B5oA+1XX+QuhGJ+l8pB0IbACGJ0CLRsAFwCXFFtLFg2AjQlTi+cAD8UZjUrHDVQ6yCf1R7UgqSHBOI02s8di8ezM9EH8+11px1cRuwOHSJpOmP7cC7iBMK2RyWlWjD6bAcwws7fi+0cIBquY/bU38IWZzTGz5cBjhP4rdl9lKK1viv4dkDQcOAg4JhrPYuvamnCj8X787LcF3pW0WZF1zQAei9OLbxNmNVpWhSY3UOkgn9QfVU68C7oTmGpm1yV2JdOODAP+VZ26zOx8M2trZh0IffMfMzsGGENIg1IsXd8CX0vqHIsGAh9S3P76CugraYP4/8xoKmpfJSitb54Ejo3eaX2BBYmpwCpH0n6EKeRDzOynLL1HSWokqSOwLfB2dWgys8lm1trMOsTP/gyCE9O3FLe/niA4SiCpE8E5aC5V0VdV4fnhrwp5yxxA8Jr7DLiwSBr2IEy5TAImxtcBhOc9LwOfErx3Ni5iP/VnjRffVvELMA14mOhVVM16egDjY589AWxU7P4CLgM+Aj4gpItpVIy+Ah4gPAdbTvhxPa60viF4o90cP/+Tgd7VrGsa4flJ5nN/a6L+hVHXx8D+1akra/901njxVUt/ldJX6wH3xc/Xu8BeVdVXHurIcRzHSSU+xec4juOkEjdQjuM4TipxA+U4juOkEjdQjuM4TipxA+U4juOkEjdQjuM4TipxA+U464hCyo3fJd5vLumRso4poO0Rks6O25dL2ruS2m2jmLakKpC0uIC6L6lIKVycdOMGynHWnRbAagNlZt+Y2ZDSq1cMM7vEzF6qpOZ+D9xeSW2tK/eS6D/HyVA0A1VWcq4yjukfk3YdnyjrEcvOLufYkyQdW06dHpIOyENHb0k35qu7nLaGS7qpMtpyisbVwNaSJsbEdx0yn+v4/31CITnfdEmnSvp9jH7+pqSNY72tJT0naYKkVyVtl30SSaMkDYnb0yVdJuldSZMz9SU1id+tt+M5SkvbchjwXDzmaUk7xu33JF0Sty+XdELcPkfSOwrJ8S5LaPpVPNdESX+XVD9Lc0tJb0g6MI7aXol1P5C0Z6z2JDC0gn3v1GKKOYIaRcWiLH8AHJF4P5SQJKtMzOxWM7unnGo9CKF9ymtrvJmdXl49p85wHvCZmfUws3Ny7N8BOBTYGbgS+MlC9PM3gMxN023AaWa2E3A28Lc8zjvXzHoBt8RjIISa+Y+Z7UKIlzZSIQXIamKctB8s5vMBXiWkT2hOiOS9eyzfE3hF0iBCXLVdCN+RnST1k7Q9cCSwu5n1AFYCxyTOsynwNHCJmT0NHA08H+t2J4QUwsx+ABpJ2iSPa3bqEEUzUFZ6cq7y+BJoLGnTGAxzP9bkbkHSCfFO731JjyqErM+eyx8r6Zp45/eJpD0VgrReDhwZ7/COlLRLvPt7T9LrikFB40juqUS7d8U2P5d0ekJLzrtLSb+O532bNT8GTu1ljJktMrM5hPxM/47lk4EOCulNdgMeljQR+Dshu3F5ZKLNTwA6xO1BwHmxnbFAY6B91nFtCHmsMrxKSLy4O8GgNI3fm45m9nFscxDwHiH22nYEgzUQ2Al4J55vICHmH0BDQsy9c83sxVj2DvBrSSOAbhaSYmb4jpAexHFW06D8KtWHpHNI3IEleCVrxPIIcDhrvjDLEvseM7PbY3tXEIIb/jVHmw3MbJc4pXepme0dpzZ6m9mp8fgNgT3NbIXCw+k/EqZGstmOcLfaDPhY0i3ANqy5u1wu6W/AMZJeJATz3InwYzUmXodTe0l+Plcl3q8ifAfrERIK9qhguytZ810WcFg0LKWxhGC4MrwD9AY+B14kpE44gWD4Mm1eZWZ/TzYi6TTgH2Z2fo5zrIjH7wv8F8JNqaR+hMzIoyRdl5jVaBx1Oc5qUuUkYWYj4zRJ9it7Ou0hgoEaSoi2m2SHOIc/mWDsupZyulx3n9k0J9zVfkBIu1xaW09bSHM8l3AnuCml3132AcZayNfzM/BgKW06NYdFhJuTCmEhKeQXkg6HkPZEUvcKNvc8cFqcXUBSzxx1PiHxmY+fw68J36k3CCOqs4FXEm3+Jo70kLSFpNaEEdKQuI2kjSVtmWmWkJl2O0n/L+7fEpgdbyDvIOTOyqR52YwQrdtxVpMqAxUfxE7M8VrLIcFCPpTlwD6EL0mSUcCpZtaNMFJpTG5y3X1m8wfC9MwOwMF5tJVsT4S7y4yR7WxmI0o53qnBmNk8YFx88D+ygs0cAxwn6X1gCiErbkX4A2F6bZKkKfF9tt4fgc8kbZMofpWQtXhJ3G4b/2JmLwD3A2/EG79HgGZm9iFwEfCCpEmE0VebxHlWEm4i91Jww+9PSL73HmF24YZYdSfgTVuTkt5xgJRN8ZnZSCDfL/glQGszW6m1sw03A2YpZIY9hsIyOmbfCTdPHD+8gHYgGM5/SfqLmX2n4K3VDHgLuCE+EF5IuGst18nDSTdmdnRW0Q6xfBThpilTr0Nie/U+M/uCHE5DyZsaMxteSjvjCT/+RAPz2zwk30T4TF8Uj7sYuDhuf0O4wUrquIE1BiVZ/iA5ZgHMrGn8u4wwzZfhHzm0/C/5OYU4dYxiupk/QJhO6CxphqTjCjnezF43sydy7LqYYATGEZK2FcIYoEvGSQL4E3BVvOMryJiXdndpIevlCMK1jwOmFqjRcdYZM3uc9EypfWBm2TMhjuMJCx3HcZx0kqpnUI7jOI6TwQ2U4ziOk0rcQDmO4zippEYbKEkNJM2RdHVW+QUFtHGHpC5l7B8rqXcF9fVTiJW2QjGGWmLfc5LmqwojSheKpM5Z7v0LJZ2ZAl2NY0SO9yVNUSIWXDFRBeJJVgdp1JVGTeC6CqEYmmq0gSKsg/oEODyzMDGSl4GSVN/Mjo8ed1XBVwRX3vtz7BtJcK9NDWb2cWbdFmFtyk/A48VVBYR1ZnuZWXdCLLj9JPUtriSg4vEkq5pRpE/XKNKnCVxXIYyimjXVdAM1lLA24ytgV4A4mlo/jgBGZx8gabGkP8cFkbtmRkiS6itEi/5AITr0WVnH1Yv7r8hXnJlNN7NJhJA22fteJqy7SisDCQFQvyy2EAtk8gs1jK+iu5+uQzzJKiWNutKoCVxXIRRDU6oW6haCpMbA3oRFiS0Ixup1MztP0qllxDVrArxlZv8X28mU9wC2iFEjkNQicUwDYDRhvcaVlXoh6eUoSoaRKhoKgXYnEGIc3mxmbxVZkuM4VUxNHkEdRAhDtAR4FPilsnLRlMLKWD+bz4GtJP1V0n6EKA8Z/k4dMk4Kkd0PAR4utpYMZrYy3nS0BXaRtEORJTmOU8XUZAM1FNhb0nTCnfUmwF55HLc0xghbi5iTpjshRcFJhGCWGV4HBsRRW11gf+BdM5tdbCHZmNl8QsSPtM3PO45TydRIA6WYBgNob2YdYlyyU1iTlXN5jMVXSJstgXpm9ighRFGvxO47gWeAhyTV2GnRAsgVJb5oSGqVmXKVtD7BOabQMFaO49QwaqSBAv6HkDU0GUX8X8DBkhoRspNOyuUkUQZbAGMVUmPcB6yV48bMriPkbbpXUl79JmlnSTMIAWH/HqNLZ/a9SphCGxhjEe5bWjvViUL21X1Yk44kDbQBxsSYhu8AL5pZ0d3z1zWeZFWRRl1p1ASuK+2aPBaf4ziOk0pq6gjKcRzHqeW4gXIcx3FSSVENlKSGkq6W9GkMCfSGpP2Lqaks4kLdIeXXzKutYfG6P5U0rDLarAxSGmKlnaQxkj6MoY7OKLYmSGcIpjRqAtdV0zVBkXSZWdFewNWEDJuN4vtNgSOqWUODAuqOAoZUwjk3Jqy72hjYKG5vVMz/RUJbP4IH4wfF1pLQ1AboFbebEcJbdUmBLgFN43ZDQqLMvq7JddU2TcXSVcyMuhsAJwCnWfTGM7PZZvZQ3D8ojqjelfSwpKaxfLqky2L5ZEnbxfJfaE2Q0/ckNVNgpNaELzoy1u0v6VVJTwIfKoQ5GinpHUmTJP021pOkmyR9LOkloHUlXf6+BE+07y2sv3qRlKzrsXSGWJllZu/G7UWELMRbFFdVOkMwpVETuK5CSKMmKI6uYk7xbQN8ZWYLs3corEm6CNjbzHoB44HfJ6rMjeW3AGfHsrOBUyxEG9gTWAIcSghh1J0QFmmkpDaxfi/gDDPrBBwHLDCznYGdgRMkdSS4s3cGugDHArtVzqWzBfB14v0MUvCDWxOQ1AHoSbh7Kzrx5mYi8B3hpqPoutKoCVxXIaRRE1S/rrQ6SfQlGIVxsTOGAVsm9mfW6EwAOsTtccB1kk4HWpjZCmAP4AELYXJmA/8lGCCAt83si7g9CDg2nustQlSKbQnTXZnjvwH+U9kX6uRPHEU/CpyZ68amGFgKQzClURO4rkJIoyaofl3FNFDTgPYKUSGyEcE694ivLmaWXBSWWaC7khjw1syuBo4H1icYtu3KOf+PWec7LXG+jmb2QkUuKk9mAu0S79vGMqcUFCKDPAqMNrM0LSIG0hmCKY2awHUVQho1QfXpKpqBMrOfCCGEblAITpoJaXM48Cawu6RtYnkTSZ3Kak/S1mY22cyuIUQb2A54FTgyDktbEUZEb+c4/Hng5PgjiKROChEVXkkc3wYYUAmXnjnfIEkbSdqIMIJ7vpLarnVIEuGzMtVCRI9UoBSGYEqjpqjFddVgTVFLtesqdly5i4ArCI4KSwmjmkvMbI6k4cADCqGLMnU/KaOtMyUNIORemgI8C/xMyBP1PuFh3rlm9m2O0dUdhKnCd+OP4Rzgl4RkfXsBHxJyTr2xTlcbMbPvJf2BYEgBLjezVDgmKIQz6Q+0VAjTdKmZ3VlcVexOSO44OU7DAlxgZs8UTxIQvAv/oRBFvx7wkBU/BFMaNYHrqumaoAi6PNSR4ziOk0rS6iThOI7j1HHcQDmO4zipxA2U4ziOk0pqnIGSNELSzETUiIkZz5JKPMcFldleKedIXSw+pTQGGICkFpIekfSRpKmSdk2Bps5Zn8OFks50Ta6rtmkqlq4a5yQhaQSw2MyurcJzLDazplXY/saE6Bi9Cd6FE4CdYtijohE9GJuY2eLocv8aIdrGm8XUBSDpH8CrZnZHXJawQVyLkQqiZ9NMoI+ZfVlsPZBOTeC6CiGNmqD6dNW4EVRpSHpTUtfE+7GSeiusoborjgzekzQ47h8u6TFJz8VRzJ9i+dXA+vEOYXQ8/uk4qvhAMZ7fOpLKWHxpjQEmqTlhDdudAGb2c5qMU2Qg8FmafkRIpyZwXYWQRk1QTbpqqoE6KzHMHBPLHgSOAFBYVNvGzMYDFxLSw+9CWGg7UmERLoQ4fUcC3QgLctuZ2XnAkhhR4hiC4fjGzLqb2Q7Ac5WgP7Wx+JTOGGAdCWvT7o43GXck/odp4SjggWKLyCKNmsB1FUIaNUE16aqpBuovibBEmegODwGZXE1HAI/E7UHAefFHdyzQGGgf971sZgvMbClhMW4y3l+GycA+kq6RtKeZLaj8y0kPKY0B1oAQ3PcWM+tJWNB9XnElrSFOOR4CPFxsLRnSqAlcVyGkURNUr66aaqBKYGYzgXmSdiSMih6MuwQcljBo7c1saty3LNHE6rh+We1+QvhxnAxcIemSSpCb+lh8KYsBNgOYkRjNPUL4n6SF/YF3Y0DitJBGTeC6CiGNmqAaddUaAxV5EDgXaG5mk2LZ88Bp0QEAST3zaGe51sTl2xz4yczuA0ZSOT+MqYzFp5TGADOzb4GvJXWORQMJI960MJT0TcOkURO4rkJIoyaoRl3FjsVXUc6S9KvE+1+a2XTCnfUNwB8S+/4AXA9MklQP+AI4qJz2b4v13wXuITy3WgUsB05eV/EpjsWX1hhgAKcBo+P0wufAr4usBwiBjAmG/LfF1pIhjZrAdRVCGjVB9euqcW7mjuM4Tt2gtk3xOY7jOLUEN1CO4zhOKnED5TiO46QSN1CO4zhOKnED5TiO46QSN1CO4zhOKnED5TiO46QSN1CO4zhOKnED5TiO46QSN1CO4zhOKnED5TiO46QSN1CO4zhOKnED5TiO46QSN1CO4zhOKinTQElqJ2mMpA8lTZF0Rr4NS+ovySQdnCh7SlL/co47U9IG+Z6nECS9nkedvM4v6Q5JXSpJ13RJLSujLcdxnNpCeSOoFcD/mVkXoC9wSoE/yjOACwvUdCZQJQbKzHarrPOb2fFmlqasro7jOLWKMg2Umc0ys3fj9iJgKrBFAe2/DyyQtE/2DkkDJb0nabKkuyQ1knQ6sDkwRtKYHMdMl3SVpImSxkvqJel5SZ9JOinWaSrpZUnvxrYHJ45fHP/2lzRW0iOSPpI0WoES55d0SzzXFEmXJdoaK6l3pl1JV0p6X9KbkjaN5a0kPSrpnfjaPZZvIumF2OYdgAroU8dxnLqBmeX1AjoAXwEbxvfnABNzvG6M+/sDTwH9gP/GsqdieWPga6BTLL8HODNuTwdalqJhOnBy3P4LMAloBrQCZsfyBgmNLYFprMkcvDihbQHQlmCk3wD2yHV+YOP4tz4wFtgxvh8L9I7bBhwct/8EXBS370+02x6YGrdvBC6J2wfG43Nes7/85S9/1dVXg3LsFxBGJcCj0YgsBDCzkcDI8o41s1ckIWmPRHFn4Asz+yS+/wdwCnB9HnKejH8nA00tjOwWSVomqQXwI/BHSf2AVYQR36bAt1ntvG1mM+L1TSQY4NdynO8ISScSDF8boAvBMCb5mWB8ASYAmRHj3kAXafUAacPYl/2AQwHM7GlJP+Rx3Y7jOHWKcg2UpIYE4zTazB5LlJ8DHJPjkFfM7PSssiuBiwjPtNaVZfHvqsR25n2DqKkVsJOZLZc0nTBiK60dgJXk6AtJHYGzgZ3N7AdJo0ppa7mZWY626gF9zWxpVrulXpzjOI4TKM+LT8CdhKmp65L7zGykmfXI8co2TpjZC8BGwI6x6GOgg6Rt4vv/Bf4btxcRpu0qSnPgu2icBgBbFnh88vwbEkZkC+Jzpf0LbOsF4LTMG0k94uYrwNGxbH9C3ziO4zgJyvPi251gPPaKjgkTJR1QwXNdCbQDiCOKXwMPS5pMGP3cGuvdBjyXy0kiT0YDvWO7xwIfFXj86vOb2fvAe7GN+4FxBbZ1etQySdKHwEmx/DKgn6QphKm+rwps13Ecp9ajNTNTjuM4jpMePJKE4ziOk0rcQDmO4zipxA2U4ziOk0pqtIGS1EDSHElXZ5VfUEAbZcbUS0aMqIC+fjGixQpJQxLlPSS9ESNJTJJ0ZEXarwoktUhE2JgqadcUaKpwTMiqRtJ+kj6WNE3SecXWA+nUBK6rENKoCYqgq9grhdflRXD7Hgd8RnT4iOWL8zy+fh51xhIjRlRAXweCa/09wJBEeSdg27i9OTALaFHs/ox6/gEcH7fXS4MuwgLpXnG7GfAJ0CUFuurHz95Wsa/eL7auNGpyXTVfU7F01egRFDAUuIHgpr0rQBxNrR9d4kdnHxDj5v1Z0vvArpkRkqT6kkZJ+kAhht9ZWcfVi/uvyFecmU03s0kEN/pk+Sdm9mnc/gb4jrC4uKhIak6IcnEngJn9bGbziyqKSokJWVXsAkwzs8/N7Gfgn8Dgco6pi5rAddV0TVAEXTXWQElqTAgl9G/gAYKxwszOA5ZYWDScK9JFE+AtM+tuZsnQRj2ALcxsBzPrBtyd2NeAsL7qUzO7qJKvYxfC3chnldluBekIzAHuVgjke4ekJsUWlURSB6An8FaRpUAwkl8n3s+g+IYzjZrAdRVCGjVBEXTVWAMFHASMMbMlhFBMv5RUP4/jVsb62XwObCXpr5L2AxYm9v0d+MDMrlxX0UkktQHuBX5tZqvKq18NNAB6AbeYWU9CFI00zX+XiAnpOE7tpSYbqKHA3jHW3gRgE2CvPI5bamYrswvN7AegO+GZ00nAHYndrwMD4qitUpC0IfA0cKGZvVlZ7a4jM4AZZpYZnTxCMFhFp7SYkEVmJjE6SqRtLCsmadQErqsQ0qgJiqCrRhqo+OO+J9DezDqYWQdCNPShscry+INWSJstgXpm9ighsG3yh/lO4BngIUl5RYAv51zrAY8D95jZI+vaXmVhZt8CX0vqHIsGAkVPylhWTMgi8w6wraSO8X96FGui7bumtXFdNVsTFEHXOv/YFon/Af5jZsmI5P8C/iSpESGe3iRJ75byHCoXWxCevWSM9vnJnWZ2XXQiuFfSMflMyUnamWCINgIOlnSZmXUFjiA4I2wiaXisPtzMJuaptSo5DRgdP4CfE2ImFptMTMjJCqlRAC4ws2eKJwnMbIWkU4HnCR5Od5nZFNdUEtdVszVBcXR5LD7HcRwnldTIKT7HcRyn9uMGynEcx0klRTVQkhpKulrSpzEk0BsKCfxSSVyoO6T8mnm1NSxe96eShlVGm5WBpLskfSfpg2JryZBGTeC6CiGNmsB1FUIxNBV7BPUHQhibHcysF/BL1i2bbsFUhldeBc65MXAp0IewOvtSSWnJqjsK2K/YIrIYRfo0gesqhFGkTxO4rkIYRTVrKpqBkrQBcAJwWsYbz8xmm9lDcf+gOKJ6V9LDcZEmkqZLuiyWT5a0XSz/hdZk/X1PUjMFRibCFx0Z6/aX9KqkJ4EPFcIcjZT0jkLw1t/GepJ0k0JwxJeA1pV0+fsCL5rZ93H91Yuk5MNoZq8A3xdbR5I0agLXVQhp1ASuqxCKoamYI6htgK9yRQSIa5IuAvaOI6vxwO8TVebG8luAs2PZ2cApZtaDsEZqCSGdeg/CAty9gZEK0RsgrHM6w8w6AccBC8xsZ2Bn4ARJHQnu7J2BLoT08btVzqWnNpSJ4zhOakjrOqi+BKMwLqzRZD3gjcT+TCSBCQQjBCGq+XUKAWIfM7MZkvYAHoiRI2ZL+i/BAC0E3jazL+Kxg4AdE8+XmgPbEtYqZY7/RtJ/quBaHcdxnBwU00BNA9pL2jDHKEqEKbChOY4DyCzQXUm8BjO7WtLTwAEEw7ZvOef/Met8p5nZ82uJkA7I4zoqwkygf+J9W0KIJcdxHCdStCk+M/uJEL7mhhi1AEmtJB0OvAnsLmmbWN5EUqey2pO0tZlNNrNrCCE5tgNeBY6Mz5haEUZEb+c4/HngZMXwSJI6KUTxfiVxfBtgQCVceuZ8gyRtFJ0jBsUyx3EcJ1JsL76LCOkdPoyui08BC81sDjAceEDSJML03nbltHVmdIaYBCwHniWEGZpESKz1H+DcGG8umzsIMefejTr+ThiZPQ58Gvfdw9rTjBXGzL4neDC+E1+Xx7KiI+kBwnV2ljRD0nGuKTeuq2ZrAteVdk0e6shxHMdJJcUeQTmO4zhOTtxAOY7jOKnEDZTjOI6TSmqcgZI0QtLMRNSIiZJaVPI5LqjM9ko5R1pj8U2PUTcmShpfbD0ZJJ0RnWCmSDqzms9dIgaZpMOjllWSelennrRqcl01X1PadNU4AxX5i5n1SLzmV3L7VWqglO5YfAADYr8W5QuSjaQdCGGxdiFEBTkoswShmhhFyVBUHxAWib9SjTqSjCJ9msB1FcIo0qcJUqSrphqoEkh6U1LXxPuxknrHNVR3SXpbIUbf4Lh/uKTHJD0XRzF/iuVXA+vHEcToePzTkt6Pd/BHVoLc1MbiSynbA2+Z2U9mtgL4L2siiFQ5uWKQmdlUM/u4ujRkk0ZNUYPrypM0aooaUqOrphqosxLTe2Ni2YOEVOrERbVtzGw8cCEhPfwuhIW2I+MiXAhx+o4EuhEW5LYzs/OAJXEEcQzBcHxjZt3NbAfguUrQn+ZYfAa8IGmCpBOLLSbyAbCnpE0UggwfALQrsibHcaqYmmqgklN8megODwGZWHpHAI/E7UHAeZImEsIJNQbax30vm9kCM1tKWIy7ZY5zTQb2kXSNpD3NbEHlX06q2CMG4t0fOEVSv2ILMrOpwDXAC4QbhImEMFeO49RiaqqBKoGZzQTmSdqRMCp6MO4ScFjCoLWPP3iwJqYfJOL6ZbX7CSHy+WTgCkmXVILcmaw9Amgby4pO7EfM7DtCJI1diqsoYGZ3mtlOZtYP+AH4pNiaHMepWmqNgYo8CJwLNDezSbHseeA0KYRFl9Qzj3aWa01cvs2Bn8zsPmAkwVitK6mMxReftzXLbBN0pSKjp6TW8W97wvOn+4uryHGcKsfMatQLGEEYbUxMvDrEfZsCK4BLE/XXJ8TWmwxMAZ6K5cOBmxL1ngL6x+1rgKnAaIJDw6R4nneA3pV0Hb8hRHSfBvy62P0aNW1FiFv4fuyrC4utKaHtVcI07PvAwGo+9wPALEKMxxmE/GH/E7eXAbOB5+u6JtdV8zWlTZfH4nMcx3FSSW2b4nMcx3FqCW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJW6gHMdxnFTiBspxHMdJJf8frSGk2oB96lEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "T1, T2 = df['time'] [df['group'] == 1], df['time'] [df['group'] == 2]\n", "E1, E2 = df['status'][df['group'] == 1], df['status'][df['group'] == 2]\n", "kmf1 = KaplanMeierFitter(label='1=Maintained').fit(T1, E1)\n", "kmf2 = KaplanMeierFitter(label='2=Not maintained').fit(T2, E2)\n", "plot_fits(kmf1, kmf2, 'Kaplan-Meier survival curves', 'weeks')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The vertical ticks on the graphs indicate censored events. At the bottom, the number of patients \"at risk\" is shown. The study concerned 11 patients in group 1 with one patient censored after 161 weeks and 12 patients in group 2. A Kaplan-Meier curve estimates the survival probability as a function of time (duration). \n", "\n", "In this example, the two curves do not appear to differ very much. To analyze this more precisely one performs a [logrank test](https://en.wikipedia.org/wiki/Logrank_test):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0653393220405049" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lifelines.statistics.logrank_test(T1, T2, E1, E2).p_value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The null hypothesis is that survival is the same for both groups. Since $p=0.065$ is not particularly small (e.g., not below $\\alpha=0.05$), the null hypothesis is not strongly rejected. In other words, the logrank test also says that the curves do not differ significantly---with the usual caveat of small sample sizes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Privacy-Preserving Survival Analysis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In a multiparty setting, each party holds a private dataset and the goal is to perform a survival analysis on the *union* of all the datasets. To obtain the secure union of the datasets in an efficient way, each private dataset will be represented using two survival tables, one survival table per group of patients. \n", "\n", "The survival tables for the aml dataset are available from the Kaplan-Meier fitters `kmf1` and `kmf2`:" ] }, { "cell_type": "code", "execution_count": 6, "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", " \n", " \n", " \n", " \n", " \n", " \n", "
removedobservedcensoredentranceat_risk
event_at
0.00001111
9.0110011
13.0211010
18.011008
23.011007
28.010106
31.011005
34.011004
45.010103
48.011002
161.010101
\n", "
" ], "text/plain": [ " removed observed censored entrance at_risk\n", "event_at \n", "0.0 0 0 0 11 11\n", "9.0 1 1 0 0 11\n", "13.0 2 1 1 0 10\n", "18.0 1 1 0 0 8\n", "23.0 1 1 0 0 7\n", "28.0 1 0 1 0 6\n", "31.0 1 1 0 0 5\n", "34.0 1 1 0 0 4\n", "45.0 1 0 1 0 3\n", "48.0 1 1 0 0 2\n", "161.0 1 0 1 0 1" ] }, "metadata": {}, "output_type": "display_data" }, { "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", " \n", " \n", " \n", " \n", " \n", " \n", "
removedobservedcensoredentranceat_risk
event_at
0.00001212
5.0220012
8.0220010
12.011008
16.010107
23.011006
27.011005
30.011004
33.011003
43.011002
45.011001
\n", "
" ], "text/plain": [ " removed observed censored entrance at_risk\n", "event_at \n", "0.0 0 0 0 12 12\n", "5.0 2 2 0 0 12\n", "8.0 2 2 0 0 10\n", "12.0 1 1 0 0 8\n", "16.0 1 0 1 0 7\n", "23.0 1 1 0 0 6\n", "27.0 1 1 0 0 5\n", "30.0 1 1 0 0 4\n", "33.0 1 1 0 0 3\n", "43.0 1 1 0 0 2\n", "45.0 1 1 0 0 1" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(kmf1.event_table, kmf2.event_table)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using `mpc.input()`, each party will secret-share its pair of survival tables with all the other parties. To allow for a simple union of the survival tables, however, we modify the representation of the survival tables as follows. \n", "\n", "First, we determine when the last event happens, which is at $t=161$ in the example. We compute $maxT$ securely as the maximum over all parties, using the following MPyC one-liner:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "161" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "secfxp = mpc.SecFxp(64) # logrank tests will require fixed-point arithmetic\n", "maxT = int(await mpc.output(mpc.max(mpc.input(secfxp(int(df['time'].max()))))))\n", "timeline = range(1, maxT+1)\n", "max(timeline)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All parties use $1..maxT$ as the global timeline. Subsequently, each party pads its own survival tables to cover the entire timeline. This is accomplished by two calls to the function `events_to_table()`, which only keeps the essential information:" ] }, { "cell_type": "code", "execution_count": 8, "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", "
d1n1d2n2
1011012
2011012
3011012
4011012
5011212
6011010
7011010
8011210
911108
1001008
\n", "
" ], "text/plain": [ " d1 n1 d2 n2\n", "1 0 11 0 12\n", "2 0 11 0 12\n", "3 0 11 0 12\n", "4 0 11 0 12\n", "5 0 11 2 12\n", "6 0 11 0 10\n", "7 0 11 0 10\n", "8 0 11 2 10\n", "9 1 11 0 8\n", "10 0 10 0 8" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d1, n1 = events_to_table(maxT, T1, E1)\n", "d2, n2 = events_to_table(maxT, T2, E2)\n", "pd.DataFrame({'d1': d1, 'n1': n1, 'd2': d2, 'n2': n2}, index=timeline).head(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Column `d1` records the number of observed events (\"deaths\") for group 1 on the entire timeline, and column `n1` records the number of patients \"at risk\". Similarly, for group 2. \n", "\n", "To obtain the secure union of the privately held datasets, we let all parties secret-share their survival tables and simply add all of these elementwise:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "d1, n1, d2, n2 = (functools.reduce(mpc.vector_add, mpc.input(list(map(secfxp, _)))) for _ in (d1, n1, d2, n2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The joint dataset (union) is now represented by `d1, n1, d2, n2`. Note that these values are all secret-shared, for example:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d1[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Aggregate Kaplan-Meier Curves\n", "\n", "We now proceed to analyze the joint dataset in a privacy-preserving way by plotting *aggregated* versions of the Kaplan-Meier curves. The exact Kaplan-Meier curves would reveal too much information about the patients in the study. By aggregating the events over longer time intervals, the amount of information revealed by the curves is reduced. At the same time, however, the aggregated curves may still be helpful to see the overall results for the study---and in any case to serve as a sanity check.\n", "\n", "The function `aggregate()` securely adds all the observed (\"death\") events over intervals of a given length (stride). The aggregated values are all output publicly, and used to plot the curves via the function `events_from_table()`:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD+CAYAAABvEpGeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7eElEQVR4nO3dd7wU5dn/8c+XIigCFkBRRLAAUgQUEBUJRMESW+yG51FiiyY28lNji2JLNCRGjXk0VqKiYk2wVwwGRQVFELCgooCIiEpRUcr1++O+DwzLnnP2wJZZzvV+vfZ1Zqfcc+3snr1mZmeuW2aGc845lzZ1Sh2Ac845l40nKOecc6nkCco551wqeYJyzjmXSp6gnHPOpZInKOecc6nkCcq5SNJLkk4qwXpnSNqn2OvNlaRBkp4tdRw1JelmSb/PQzvDJV2Zj5hczXiCcsDKL+evJTUodSxrQ1IbSSapXoHaHyrpnsTzrSW9K+kGSSrEOvMhbpMvkttFUv04LqebIM1shJkNLFyUhWFmp5rZFaWOw609T1AOSW2AvQADDi7gegqSPIpN0rbAGGCUmZ1p6b/b/Wtg/8Tz/eO4gpNUt0DtStJ68/21vvxv5Nt68wa7dXIcMA4YDhyfnCBpc0mPSVoo6Q1JV0r6b2L6QEnvSVog6f8k/afiNJmkwZLGSvqrpPnAUEkNJP1Z0qeS5sbTMBsm2jtP0hxJn0k6KR4B7BCn/UzSWzGWmZKGJkIdE/9+I2mxpN3jMidImhaPDp+JyaViXQPiUdACSTcC1R4JSdo+rmuEmZ2XGH99jGmhpAmS9kpMGyrpIUkjJS2S9KakrpW030vSq5K+idvhRkkbJKabpFMlfRDn+XsOR3B3E97jCscBd2Wst6mk2+M6Z8f3uW6cNjjjPe8g6TlJX8X3/qjEtOGSbpL0pKRvgf5ZXuNgSR/FbfGxpEGJ7ZQ8Sl3tqFjhKP8qSWOB74BzJY3PaHuIpFGJWK6Mw9MkHZiYr56keZJ2ic8flPR5/CyMkdSpmm2aXOfJsf1FkqYm2lz52c0STz9JsyT9TtLnwJ05xNhb0ivxfX9bUr/qtmnZMzN/1PIHMB34NbArsBTYIjHt/vjYCOgIzAT+G6c1AxYChwH1gLPi8ifF6YOBZcAZcfqGwF+BUcBmQGPgMeCPcf79gM+BTnF99xCO6naI0/sBXQg7VjsDc4FD47Q2cd56idgPia9tp7j+i4FXErEvAo4A6gNDYqwnVbKNhgKvALOBC7NM/x9g87ie/xdfR8PEsksT6zoH+BioH6fPAPaJw7sCvWM7bYBpwNmJ9RjwOLAJ0BqYB+xXxXtrQOe4rTYBNo3DncO//8r5HgX+ATQCWgCvA79KvI8V73mj+Bn4ZYyxO/Al0DFOHw4sAPaM71PDjHgaET4z7ePzlkCnxHa6JzHvau8p8BLwKeHzUQ9oGt/DHRPLvAEck4jlyjh8CWGnomK+nwHTEs9PIHweGwDXARMT01a2k2X7Hhk/Ez0JOzg7ANsmtv0O2dohfJaXAdfEdW5YVYzA1sB84IC4XQfE582r2qbl/ih5AP4o8QcA+hC+PJvF5+8CQ+Jw3TitfWL+KxNfVscBryamKX55JRPUpxnTvwW2T4zbHfg4Dt9BTFbx+Q6Z/+QZsV8H/DUOr/ZlFsc9BZyYeF6HsOe9bYx9XEZss6g6QS0EvknGX8V2/Rromlh2XEYcc4C94vMZxASVpZ2zgUcTzw3ok3j+AHB+FXFY3I63Ab8CTgVurdi2cZ4tgB+ADRPLHQuMTryPFe/50cDLGev4B3BpHB4O3FVFPI3iNjw8ub7EdqouQV2escw9wCVxeEdCwtooEUtFQtghY9qIiuWyxLhJXG/TzHayzPsMcFZV2z7xPBlPP+BHEgm8qhiB3wF3Z1n38VVt03J/+Ck+dzzwrJl9GZ/fy6rTfM0Je6ozE/Mnh7dKPo/fdrMy2k/O35xwZDQhnqb4Bng6jl+jvYxhJO0maXQ87bGA8GXbrIrXti1wfWJdXxES0daVxD4zWyMJowhJ9MXkqcIY2znxFM2CuK6mGbEl17WCsJ22ylyBpHaSHo+nmxYCf8jyGj9PDH8HbByXnaJwenNx8hRjdBchKa9xeo+wneoDcxLb6h+EI6lM2wK7VcwX5x0EbJnttWYys28JSe7UuL4nJHWobP4sMtu+l5BMAX4B/MvMvsuy3umEo9GDJG1E+K31Xgi/k0m6WtKHcZvPiItV9dmqsA3wYQ3iT5pnZktyiZGw3Y/M2O59gJZ52Kap5T/M1WIKv/0cBdSN58EhnG7YROE3kncIpyFaAe/H6dskmpgTp1W0p+TzKHkBwZfA94TTD7OzhLRaexnrgvDPeiOwv5ktkXQdq75EjDXNBK4ysxGZEyTtmGw/xp65vjWY2W8VrnR8UVJfM5sdk8F5wN7AFDNbIelrVv9NK7muOvF1fpZlFTcBbwHHmtkiSWcTTg1Wy8yq+t3kZcKpHwP+C2yfmDaTcATVzMyWVbOamcB/zGxAVaFUE+czwDPx83cl4YhuL8LR9UaJWbfMtnjG8+eA5pK6ERLVkCpWfV+cpw4wNSYECIntEGAfQnJqSjgCzuXqzJmsvi2TvmPN15Pcgcu2nSqLcSbhCOrkbCuqYpuWNT+Cqt0OBZYTflvqFh87Eb7MjjOz5cAjhIsbNop7Zckf258Aukg6NP6Q/Ruyf6kAK48cbgX+KqkFrLxce984ywPALyXtFPcgM+9haQx8FZNTL8IXS4V5wApgu8S4m4ELKn7wVrgQ4MhE7J0kHRZjP7Oq2DOcDowGXpC0RYxrWYyhnqRLgCYZy+yaWNfZhIQwLkvbjQmnEhfH7X1ajjFVKR4hHgQcHIeT0+YAzwJ/kdREUh1J20v6SZamHgfaSfpfhcvV60vqKWmnXOKQtIWkQyQ1ImyDxYT3DWAi0FdSa0lNgQtyeF1LgQeBYYTfNZ+rYvb7gYGEbXpvYnzjGMt8QkL5Qy6vJboNOEfSrgp2SBxdTwR+EY/Q9gOybc9cY7yHcGS1b2yvYbzQolU127SseYKq3Y4H7jSzT83s84oH4ShlUPwyPZ2wR/k54Wqw+wj/BMTTgkcCfyL8c3cExldMr8TvCBcujIunU54H2sf2ngJuIHz5T2fVF3hFe78GLpe0iPCD8gMVjcbTOlcBY+MpkN5m9ijhR+j747reIV5unYj96hj7jsDYXDZa/II/hXAhwfPABMKpyveBT4AlrHkq6t+E0zBfA/8LHBa/XDOdQ0i8iwjJfGQuMeUY9xQzm1LJ5OOADYCpMcaHCEdcmW0sInyBHkM4AvycVT/056IO8Nu47FeEL+3TYtvPEV7vJMI2fTzHNu8lHP08WNURYEzErwJ7sPp2vYvwvs0mvP5sOw6Vtfkg4XN3L+E9+xchUUK4aOggwu9Dg+K06trLGqOZzSQc5V1I2BGaCZxL2J6VbtNyp4ydKeeqJOkaYEszOz7LtDqEUxiDzGx0Hta1EyGpNMjh1FNqKVwOv4OZ/U+pY3GunPgRlKuSwj0vO8fTF72AEwmXJFdM31fSJvF3mQsJ5+1z3gPNsr6fK9wrtSlhz/yxck5Ozrm15wnKVacx4XeobwmnHP5COF1VYXfCVUxfEk5nHGpm36/D+n4FfBHbXM56cqrCOVdzforPOedcKvkRlHPOuVTyBOWccy6Vau2NupLuAA4EvjCzzlmmC7ieUPvqO2Cwmb1ZXbvNmjWzNm3a5Dla55xbf02YMOFLM2ueOb7WJihCXawbWbPsS4X9CffG7AjsRrjDf7fqGm3Tpg3jx4+vbjbnnHORpE+yja+1CcrMxij0g1SZQwhFL41wU+kmklrGG+nybtz/nUzjb6YVoulKLd7x5+x25P8r6jqdcy5XtTZB5WBrVq8GMCuOWyNBSTqFUFmA1q1br/UKVywr3u0+2y7/hE+mPUjoGcI559LHE1QemNktwC0APXr0WKvr9nv/+ta8xlSdKX/oA0VMiM45V1OeoCo3m9WrW7eK45xz62Dp0qXMmjWLJUuWVD+zW680bNiQVq1aUb9+/Zzm9wRVuVHA6ZLuJ1wcsaBQvz85V5vMmjWLxo0b06ZNG1Rtb/VufWFmzJ8/n1mzZtG2bduclqm1CUrSfYReLZtJmgVcSui0DTO7GXiScIn5dMJl5r8sTaSFY8CoicU5KGy8YX36t8/W/52rbZYsWeLJqRaSxOabb868efNyXqbWJigzO7aa6Ubo32j9ZUbzxg2Lsqp5i/x0jlvFk1PtVNP33StJOOdqnRNOOIEWLVrQufMa9+hX6qWXXkISt91228pxEydORBJ//vOfq1z25ptv5q67KrvlclVbTz75ZLVxjB8/njPPPDO3oKsxfPhwTj/99Ly0VQieoJxztc7gwYN5+umna7xc586deeCBlf1kct9999G1a9dqlzv11FM57rjjqpwn1wTVo0cPbrjhhuqDXQ94gnLO1Tp9+/Zls802q37GDNtuuy1Llixh7ty5mBlPP/00+++//8rpt956Kz179qRr164cfvjhfPfddwAMHTp05VFWv379+N3vfkevXr1o164dL7/8Mj/++COXXHIJI0eOpFu3bowcOZLXX3+d3Xffne7du7PHHnvw3nvvAeFI7sADD1zZ7gknnEC/fv3YbrvtVktc99xzD7169aJbt2786le/Yvny5QDceeedtGvXjl69ejF2bE6dSJdMrf0NysFS898BXGld9tgUpn62MK9tdtyqCZce1KnGyw0bNowRI0asMb5v376rffEfccQRPPjgg3Tv3p1ddtmFBg1W9XZ/2GGHcfLJJwNw8cUXc/vtt3PGGWes0eayZct4/fXXefLJJ7nssst4/vnnufzyyxk/fjw33ngjAAsXLuTll1+mXr16PP/881x44YU8/PDDa7T17rvvMnr0aBYtWkT79u057bTTmD59OiNHjmTs2LHUr1+fX//614wYMYIBAwZw6aWXMmHCBJo2bUr//v3p3r17jbdVsXiCqsWWrih1BM6lx7nnnsu5555b7XxHHXUURx99NO+++y7HHnssr7zyyspp77zzDhdffDHffPMNixcvZt99983axmGHHQbArrvuyowZM7LOs2DBAo4//ng++OADJLF06dKs8/3sZz+jQYMGNGjQgBYtWjB37lxeeOEFJkyYQM+ePQH4/vvvadGiBa+99hr9+vWjefNQl/Xoo4/m/fffr/Y1l4onKOdcyazNkU6h5HoEteWWW1K/fn2ee+45rr/++tUS1ODBg/nXv/5F165dGT58OC+99FLWdVUcddWtW5dllVR0+f3vf0///v159NFHmTFjBv369auyrWR7Zsbxxx/PH//4x9Xm/de//pW1jbTyBOWcc+R+BAVw+eWX88UXX1C3bt3Vxi9atIiWLVuydOlSRowYwdZbb53z+hs3bsyiRYtWPl+wYMHK5YcPH55zOwB77703hxxyCEOGDKFFixZ89dVXLFq0iN12242zzjqL+fPn06RJEx588MGcLvIoFb9IwjlX6xx77LHsvvvuvPfee7Rq1Yrbb7+9RsvvscceHHrooWuMv+KKK9htt93Yc8896dChQ43a7N+/P1OnTl15kcR5553HBRdcQPfu3Ss9yqpMx44dufLKKxk4cCA777wzAwYMYM6cObRs2ZKhQ4ey++67s+eee7LTTjvVqN1iU7gf1eVLjx49rBz6g5ryhz4sWrKMyzb/U1HWt/PWTbnmiPTuqbnimTZtWuq/GF3hZHv/JU0wsx6Z8/opvlrgr8+9z/UvfLDauPs3CHtk0+YsWmP+ZhtvQPPGDdYYv7Y+mf8dS5f5FRnOuZrxBFULDBnQjiED2q0+8s7/Y9zH83l2n9zrYq2tc8Y3ZemPPxR8Pc659YsnqFpu6YbNC74Oqwta6gnKOVczfpGEc865VPIjKFcUhhWtaw/w7j2cWx94gnJFUb9OnaJ17QHevYdz6wM/xVeLNW+wvNQhOFd0M2fOpH///nTs2JFOnTpx/fXX57xsRZcbjz322MpxBx54YKUVIypcd911KwvH5tsee+xR7Ty5rv+kk05i6tSp+QiLNm3a8OWXX65TG56garEWDWp2859z64N69erxl7/8halTpzJu3Dj+/ve/1+hLuVWrVlx11VU1WmchE1Sy1NK6rv+2226jY8eO+QgrLzxBOedqlZYtW7LLLrsAobzQTjvtxOzZuf8+2rVrV5o2bcpzzz23xrQXXniB7t2706VLF0444QR++OEHbrjhBj777DP69+9P//7911imTZs2XHDBBXTr1o0ePXrw5ptvsu+++7L99ttz8803A7B48WL23ntvdtllF7p06cK///3vlctvvPHGQDi669evH0cccQQdOnRg0KBBmFnW9Z922mn06NGDTp06cemll65sq1+/flQUGth444256KKL6Nq1K71792bu3LkAzJs3j8MPP5yePXvSs2fPlV12zJ8/n4EDB9KpUydOOukk8lEEwn+Dcs6VzlPnw+eT89vmll1g/6tzmnXGjBm89dZb7LbbbkDuBWMvuugifv/73zNgwICV45YsWcLgwYN54YUXaNeuHccddxw33XQTZ599Ntdeey2jR4+mWbNmWeNo3bo1EydOZMiQIQwePJixY8eyZMkSOnfuzKmnnkrDhg159NFHadKkCV9++SW9e/fm4IMPXqML9bfeeospU6aw1VZbseeeezJ27FjOPPPMNdZ/1VVXsdlmm7F8+XL23ntvJk2axM4777xaW99++y29e/fmqquu4rzzzuPWW2/l4osv5qyzzmLIkCH06dOHTz/9lH333Zdp06Zx2WWX0adPHy655BKeeOKJGpePysYTlHOuVlq8eDGHH3441113HU2aNAFyLxjbt29fAP773/+uHPfee+/Rtm1b2rULN8Uff/zx/P3vf+fss8+utr2DDz4YgC5durB48WIaN25M48aNadCgAd988w2NGjXiwgsvZMyYMdSpU4fZs2czd+5cttxyy9Xa6dWrF61atQKgW7duzJgxgz59+qyxvgceeIBbbrmFZcuWMWfOHKZOnbpGgtpggw1Wdoy46667rjxifP7551c7Jbpw4UIWL17MmDFjeOSRR4DQBcimm25a7euujico51zp5Hikk29Lly7l8MMPZ9CgQSv7ZoLcj6AgHEVdeeWV1Ku37l+jFV1m1KlTZ7XuM+rUqcOyZcsYMWIE8+bNY8KECdSvX582bdqwZMmaV6pm63oj08cff8yf//xn3njjDTbddFMGDx6cta369euvPEJLtrVixQrGjRtHw4aFvyrXf4OqxTb48etSh+Bc0ZkZJ554IjvttBO//e1vV5t27rnnMnHixDUemckJYODAgXz99ddMmjQJgPbt2zNjxgymT58OwN13381PfvITYM2uNGpqwYIFtGjRgvr16zN69Gg++eSTGi2fXP/ChQtp1KgRTZs2Ze7cuTz11FM1amvgwIH87W9/W/l84sSJQEji9957LwBPPfUUX3+97t8vnqBqsQaeoFwtNHbsWO6++25efPFFunXrRrdu3XjyySfXqq2LLrqImTNnAtCwYUPuvPNOjjzySLp06UKdOnU49dRTATjllFPYb7/9sl4kkYtBgwYxfvx4unTpwl133VXjrjyS6+/atSvdu3enQ4cO/OIXv2DPPfesUVs33HAD48ePZ+edd6Zjx44rL+S49NJLGTNmDJ06deKRRx6hdevWNWo3G+9uI8/KpbsN7vwZfPJfpgy4t+CrOv8V0PIfueDIvgVfV4V5i5ZwcLfcO4tzxePdbdRuNeluw4+gnHPOpZInKOecc6nkV/HVcm3GX1HwdfzhexhTZzegeKf4nHPlzxNUbTD6j/Cf7JfzNvp62hrjfmzYLK/9RLVd8Qk/rjB+zFuLrtyZ2Ro3mbr1X02vefAEVRv0vyA8Mg1tWpSLJFa8cAVLveyfixo2bMj8+fPZfPPNPUnVImbG/Pnza3T/lCcoVzSbf/xY9TPlydIf6zOKfkVbn8tdHavL5gvn8cnszwG/inh9UEeiYf261c7XsGHDlZUuclGrE5Sk/YDrgbrAbWZ2dcb0wcAwoKKS5I1mdltRg1yPFKN7+QpbMo/6Rex/ytVUIz/lux6Zt2gJB++U/9s6am2CklQX+DswAJgFvCFplJll1t0faWanFz1A55yr5WptggJ6AdPN7CMASfcDhwD56a3LreH86rutyZufNm9I77bFW59zLv9qc4LaGpiZeD4L2C3LfIdL6gu8Dwwxs5mZM0g6BTgFyEt5j2L5rOWA6meqoRHvwb3vrz7u/g3C38nz15y/xYawxUb5jeGjBaDlDeid32adc0VWmxNULh4D7jOzHyT9Cvgn8NPMmczsFuAWCKWOihvi2vtsq33ZMM9tDmofHkltxofk9MRBeV5ZJc5/BfDe7J0re7W5ksRsYJvE81asuhgCADObb2Y/xKe3AbsWKTbnnKv1anOCegPYUVJbSRsAxwCjkjNIapl4ejCw5l2tzjnnCqLWnuIzs2WSTgeeIVxmfoeZTZF0OTDezEYBZ0o6GFgGfAUMLlnAzjlXy9TaBAVgZk8CT2aMuyQxfAGQpQTD+mHDDeryzXeFvxtl2fIV1O6Ddefc2qjVCaq267RVE9i4WeFXNL0+zZYs44fq53TOuZV8t9YVRfOGZXNxo3MuJTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RbL32+pPrePZ1z+fHk5DkFadcTlFsvfeEJyrmieXrK3IK065UkXFE0+vZT2oy/oijr+sP3cG/dPcB7hHKurHmCcoW3XT++/X5p0T5sbVd8wiF1wWtXOFfePEHVZg2bwuLCHJqvZqvufPT9lmzQfLvCrwtY8cIVsKIoq3LOFZAnqNpsx/x3+V6Z5TNuLdq6Klz++JSir9M5lz+eoFxZG/Ee3Pv+6uPu3yD8nTZn0RrzN9t4A5o3blCEyJxb/8xb9ANfLs7eRU+b859YY9xZe+/IkAHt1np9nqBcWRvUPjyS2oyHyfPhvpP9IgnniuHYW8cx4+qf5b1dv8zcOedcKnmCcs45l0qeoJxzzqWS/wblimLDDeryzXfZf1zNt2XLV+D7Xs6VP09Qrig6bdUENm5WnJVNr0/z73/ki+Kszblab79OWxSkXd/NdOulFg2WlToE52qNA7q0LEi7nqCcc86lkico55xzqeQJyjnnXCp5gnLOOZdKnqCcc86lkico55xzqeT3QbniKFbfUwDLl2LyLt+dK3eeoFxxFLHvKcbfif04r3jrc84VRK1OUJL2A64H6gK3mdnVGdMbAHcBuwLzgaPNbEax43Q1V6eOmLdoSanDcK5WaLxh/YK0W2sTlKS6wN+BAcAs4A1Jo8xsamK2E4GvzWwHSccA1wBHFz9aV1ONG9Tj4G5blzoM59w6qM0XSfQCppvZR2b2I3A/cEjGPIcA/4zDDwF7S1IRY3TOuVqrNieorYGZieez4ris85jZMmABsHlmQ5JOkTRe0vh58/y3j5Lbsgtstl2po3DOraNae4ovn8zsFuAWgB49eliJw3H7X139PM651KvNR1CzgW0Sz1vFcVnnkVQPaEq4WMI551yB1eYE9Qawo6S2kjYAjgFGZcwzCjg+Dh8BvGhmfoTknHNFoNr8fSvpAOA6wmXmd5jZVZIuB8ab2ShJDYG7ge7AV8AxZvZRNW3OAz5Zy5CaAV+u5bLF5rEWRjnFCuUVr8daGPmIdVsza545slYnqLSRNN7MepQ6jlx4rIVRTrFCecXrsRZGIWOtzaf4nHPOpZgnKOecc6nkCSpdbil1ADXgsRZGOcUK5RWvx1oYBYvVf4NyzjmXSn4E5ZxzLpU8QTnnnEslT1ApIWk/Se9Jmi7p/FLHkyRpG0mjJU2VNEXSWXH8UEmzJU2MjwNKHSuApBmSJseYxsdxm0l6TtIH8e+mKYizfWLbTZS0UNLZadmuku6Q9IWkdxLjsm5HBTfEz+8kSbukINZhkt6N8TwqaZM4vo2k7xPb9+YUxFrpey7pgrhd35O0bwpiHZmIc4akiXF8/rermfmjxA/CjcIfAtsBGwBvAx1LHVcivpbALnG4MfA+0BEYCpxT6viyxDsDaJYx7k/A+XH4fOCaUseZ5TPwObBtWrYr0BfYBXinuu0IHAA8BQjoDbyWglgHAvXi8DWJWNsk50vJds36nsf/s7eBBkDb+D1Rt5SxZkz/C3BJobarH0GlQy5df5SMmc0xszfj8CJgGmtWfk+7ZNcp/wQOLV0oWe0NfGhma1uFJO/MbAyhgkpSZdvxEOAuC8YBm0hqWZRAyR6rmT1roRcCgHGEepslV8l2rcwhwP1m9oOZfQxMJ3xfFEVVscauh44C7ivU+j1BpUMuXX+kgqQ2hNJPr8VRp8dTKHek4bRZZMCzkiZIOiWO28LM5sThz4EtShNapY5h9X/0NG5XqHw7pv0zfALhCK9CW0lvSfqPpL1KFVSGbO95mrfrXsBcM/sgMS6v29UTlMuZpI2Bh4GzzWwhcBOwPdANmEM43E+DPma2C7A/8BtJfZMTLZyPSM39FbFY8cHAg3FUWrfratK2HSsj6SJgGTAijpoDtDaz7sBvgXslNSlVfFFZvOcZjmX1naq8b1dPUOmQS9cfJSWpPiE5jTCzRwDMbK6ZLTezFcCtFPHUQ1XMbHb8+wXwKCGuuRWnnOLfL0oX4Rr2B940s7mQ3u0aVbYdU/kZljQYOBAYFBMq8XTZ/Dg8gfC7TruSBUmV73lat2s94DBgZMW4QmxXT1DpkEvXHyUTzzXfDkwzs2sT45O/MfwceCdz2WKT1EhS44phwg/l77B61ynHA/8uTYRZrbYnmsbtmlDZdhwFHBev5usNLEicCiwJSfsB5wEHm9l3ifHNJdWNw9sBOwJV9lJQaFW856OAYyQ1kNSWEOvrxY4vi32Ad81sVsWIgmzXYl0N4o9qr5Y5gHB13IfARaWOJyO2PoRTOZOAifFxAKErkslx/CigZQpi3Y5w1dPbwJSKbQlsDrwAfAA8D2xW6lhjXI0InWA2TYxLxXYlJM05wFLCbx8nVrYdCVfv/T1+ficDPVIQ63TC7zcVn9mb47yHx8/GROBN4KAUxFrpew5cFLfre8D+pY41jh8OnJoxb963q5c6cs45l0p+is8551wqeYJyzjmXSp6gnHPOpZInKOecc6nkCco551wqeYJyzjmXSp6gnFtHkjaR9OvE860kPZSntodKOicOXy5pnzy121LS4/loq5L2F9dg3udTVm/QpYQnKOfW3SbAygRlZp+Z2RH5XomZXWJmz+epud8SSuqkwd0ktp9zFUqWoLJ1hJXDMv0kmaSTEuO6xXHnVLPsqZKOq2aebsqhczhJPSTdkGvc1bQ1WNKN+WjLlczVwPaxk7ZhseO2d2Dl+/svhc79Zkg6XdJvY8XncZI2i/NtL+npWIH9ZUkdMlciabikI+LwDEmXSXpToXPGDnF8o/i/9XpcR2XdthwOPB2XeULSznH4LUmXxOHLJZ0ch8+V9Eastn1ZIqb/ieuaKOkfFaVuEtObSXpV0s/iUduYOO87iWrXowjlnpxbTSmPoIYD+63Fcu8Q+iCpcCyhrE2VzOxmM7urmtm6EUr4VNfWeDM7s7r5XK1xPqEvp25mdm6W6Z0JhTV7AlcB31mo+PwqULHTdAtwhpntCpwD/F8O6/3SQtX2m+IyEMrivGhmvYD+wLBYk3ClWNPtazP7IY56GdhLUlNC1e894/i9gDGSBhLqqvUi/I/sKqmvpJ2Ao4E9zawbsBwYlFjPFsAThA7tngB+ATwT5+1KKImDmX0NNJC0eQ6v2dUiJUtQVrNOu5I+ARpK2iIWMd2PRD8vkk6Oe3pvS3pY0kZxfPJc/kuSrol7fu9L2kuhSOvlwNFxD+9oSb3i3t9bkl6R1D4u36/i/H1s947Y5keSzkzEknXvUtIv43pfZ9WXgVt/jTazRWY2D1gAPBbHTwbaKHRjsgfwoEL32f8g9GJcnUfi3wmE3kwhFMc9P7bzEtAQaJ2xXEtgXuL5y4SeU/ckJJSN4/9NWzN7L7Y5EHiLUGOtAyFh7Q3sCrwR17c3oRYiQH1Czb7zzOy5OO4N4JeShgJdLHR+WeELYKscXrOrReqVOoAkSeeS2ANLGJNxxPIQcCSr/mF+SEx7xMxuje1dSSjE+LcsbdYzs17xlN6lZrZPPLXRw8xOj8s3AfYys2UKP07/gXBqJFMHwt5qY+A9STcBO7Bq73KppP8DBkl6DriM8I+9ABgdX4dbfyU/nysSz1cQ/gfrAN/EI4u1aXc5q/6XBRweE0tlvickrgpvAD0IlaefA5oBJxMSX0WbfzSzfyQbkXQG8E8zuyDLOpbF5fcF/gNhp1Shb66fAcMlXZs4q9EwxuXcSqm6SMLMhsXTJJmPzNNpDxASVGaHWQCd4zn8yYRk16mS1WXb+8zUlLBX+w7w1yraesJCXyhfEvYEt6DyvcvdgJfMbJ6F7t1HVtKmKx+LCDsna8VC548fSzoSQvcmkrquZXPPAGfEswtI6p5lnvdJfObj53Am4X/qVcIR1TnAmESbJ8QjPSRtLakF4QjpiDiMpM0kbVvRLKEX2w6Sfhenb0vogfVW4DZgl4rXC2wJzFjL1+zWU6lKUPGH2IlZHqtdkGBmnxPKvw8g/JMkDQdON7MuhCOVhmSXbe8z0xWE0zOdgYNyaCvZngh7lxVJtr2ZDa1keVfGLHTSNjb+8D9sLZsZBJwoqaKbkMoubqjOFYTTa5MkTYnPM+P9FvhQ0g6J0S8DX5jZ93G4VfyLmT0L3Au8Gnf8HgIam9lU4GLgWUmTCEdfLRPrWU7YifypwmX4/YC3Jb1FOLtwfZx1V2CcmS1by9fs1lOpOsVnZsOAXP/BLwFamNnyuLNYoTEwR6EH2EHUrPfJzD3hponlB9egHQiJ89+S/mpmXyhcrdUYeA24Pv4gvJCw11rtRR4u3czsFxmjOsfxwwk7TRXztUkMr5xmZh+T5aKh5E6NmQ2upJ3xhC9/YoL5VQ4h30j4TF8cl/s98Ps4/BlhBysZx/WsSijJ8SPJchbAzDaOf38gnOar8M8ssfwvuV0U4mqZUl5mfh/hdEJ7SbMknViT5c3sFTP7V5ZJvyckgbHAuzUMazTQseIiCeBPwB/jHl+Nknlle5cWehkdSnjtY4FpNYzRuXVmZo+SnlNq75hZ5pkQ57zDQuecc+mUqt+gnHPOuQqeoJxzzqWSJyjnnHOpVNYJSlI9SfMkXZ0x/sIatHGbpI5VTH9JUo+1jK+vQq20ZYo11BLTnpb0jQpYUTqfJA2RNCVeSn2fpMouuS8prUWNx1LxWAvDYy2cYsdb1gmKcB/U+8CRFTcmRjklKEl1zeykeMVdIXxKuJT33izThhEur009SVsDZxKqbHQG6gLHlDaqSg1n7Wo8lsJwPNZCGI7HWijDKWK85Z6gjiXcm/EpsDtAPJraMF4qPiJzAUmLJf0l3hC5e8URkqS6CtWi31GoDj0kY7k6cfqVuQZnZjPMbBKhpE3mtBcI912Vi3qE7VoP2Aj4rMTxZLUONR6LzmMtDI+1cIodb6pu1K2JeIppH8JNiZsQktUrZna+pNOrqGvWCHjNzP5fbKdifDdg63iEgKRNEsvUA0YQ7te4Kq8vpAyY2WxJfybsCHwPPBurCzjnXMGU8xHUgYQyRN8DDwOHKqMvmkosj/Nn+gjYTtLfJO1HqPJQ4R/U0uQEoNDb6SFAW0LF6UaS/qe0UTnn1nflnKCOBfaRNINQ8HVz4Kc5LLck1ghbTeyTpiuhi4JTCcUsK7wC9E/rhQFFsA/wcSxwu5RQaHePEsfknFvPlWWCUuwGA2htZm1iXbLfsKpXzqWxFl9N2mwG1DGzhwklinZJTL4deBJ4IP4GU9t8CvSWtFG8GGVvvESTc67AyjJBAT8n9BqarCL+b+AgSQ0IvZNOynaRRBW2Bl5S6BrjHmC1Pm7M7FpCv013S8ppu0nqKWkWoSDsP2J16YppLwMPAnvHWoT7VtZOqZnZa4QK1m8SOtmrQ9jGqbOuNR6LyWMtDI+1cIodr9fic845l0rlegTlnHNuPecJyjnnXCqVNEFJqi/pakkfxJJAr0rav5QxVSXeqHtE9XPm1Nbx8XV/IOn4fLRZKOVUjkVSQ0mvS3o7lma6rNQxVcZjLQyPtTBKEWupr0i7gtBFdGcz+0HSFsBPihmApHrF7mpaoXfdS4EegAETJI2Kl7qn0XBCD6x3lTiOXPwA/NTMFscrOf8r6SkzG1fqwLLwWAvDYy2Mosdayh51NwJOBs6ouBrPzOaa2QNx+sB4RPWmpAclbRzHz5B0WRw/WVKHOP4nCuWNJkp6S1JjBcO0qnzR0XHefpJeljQKmKpQ5miYpDckTZL0qzifJN0o6T1JzwMt8vTy9wWeM7OvYlJ6jhTX4yqnciwWLI5P68dHKq8E8lgLw2MtjFLEWspTfDsAn5rZwswJCvckXQzsY2a7AOOB3yZm+TKOvwk4J447B/hNLHG0F6Ekz2GEEkZdCTebDpPUMs6/C3CWmbUDTgQWmFlPoCdwsqS2hMvZ2wMdgePI382pWwMzE89nxXEuD+IOx0TgC8KOwGslDqlSHmtheKyFUexY03qRRG9CUhgbN8bxwLaJ6Y/EvxOANnF4LHCtpDOBTeJpuz7AfWa23MzmAv8hJCCA183s4zg8EDgurus1QlWKHYG+ieU/A17M9wt1+Rffr25AK6CXpM4lDqlSHmtheKyFUexYS5mgpgOtFapCZBIhO3eLj45mlrwhrOIG3eXE39HM7GrgJGBDQmLrUM36v81Y3xmJ9bUtcDHU2cA2ieet4jiXR2b2DTCaFJ8+reCxFobHWhjFirVkCcrMviOUELpe0gYAkppLOhIYB+wpaYc4vpGkdlW1J2l7M5tsZtcAbwAdgJeBo+NhaXPCEdHrWRZ/Bjgt/vCHpHaSGgFjEsu3BPrn4aVXrG+gpE0VCrEOjOPcOoqfoU3i8IaEPsPeLWlQlfBYC8NjLYxSxFrqq/guBq4kXKiwhHBUc4mZzZM0GLhPoXRRxbzvV9HW2ZL6E/pemgI8BfxI6CfqbcKPeeeZ2edZjq5uI5wqfFOSgHnAocCjhAK0Uwn16F5dp1cbmdlXkq4gJFKAy80stRchKJQ36Qc0UyjddKmZ3V7aqCrVEvinQmX7OsADZpbWXos91sLwWAuj6LF6qSPnnHOplNaLJJxzztVynqCcc86lkico55xzqVR2CUrSUEmzE1UjJlZcWZLHdVyYz/YqWUc51eLbRtJoSVMVanCdVeqYKiOpfcZnY6Gks0sdVzYea2F4rIVRiljL7iIJSUOBxWb25wKuY7GZbVzA9jcjVMdYWYsP2DWttfjiJfYtzexNSY0J8R5qZlNLHFqV4tVGs4HdzOyTUsdTFY+1MDzWwihWrGV3BFUZSeMkdUo8f0lSj3gP1R0KVXjfknRInD5Y0iOSno5HMX+K468GNox7CCPi8k8oVPB9R7Ge3zoqt1p8c8zszTi8iNDdezmUZtob+DDt/+yRx1oYHmthFCXWck1QQxKHmaPjuJHAUbDaHv944CJC9/C9CDfaDlO4CRdCnb6jgS6EG3K3MbPzge9jRYlBhMTxmZl1NbPOwNN5iL9sa/FJagN0J5SESrtjgPtKHUSOPNbC8FgLoyixlmuC+muiLFFFdYcHgIq+mo4CHorDA4HzFersvQQ0BFrHaS+Y2QIzW0K4GTdZ76/CZGCApGsk7WVmC/L/csqDQkX5h4GzsxX5TROF6iQHAw+WOpbqeKyF4bEWRjFjLdcEtQYzmw3Ml7Qz4ahoZJwk4PBEQmttZtPitB8STays65fR7vuEyueTgSslXZKHcMuuFp9CGaiHgRFm9kh186fA/sCbsUhw2nmsheGxFkbRYl1vElQ0EjgPaGpmk+K4Z4AzYgkjJHXPoZ2lWlWXbyvgOzO7BxhGSFbrqqxq8cVtdzswzcyuLXU8OTqW8jld4rEWhsdaGEWLtVwTVPI3qInxdxEIp/WOIZzuq3AFoWOtSZKmxOfVuSXOP4Lw+9Tr8RThpYTagesk1t2rqMX3BimvxQfsCfwv8NPENj+g1EFVJv7GOIBV3bKklsdaGB5rYRQ71rK7zNw551ztUK5HUM4559ZznqCcc86lkico55xzqeQJyjnnXCp5gnLOOZdKnqCcc86lkico55xzqeQJyjnnXCp5gnLOOZdKnqCcc86lkico55xzqeQJyjnnXCp5gnLOOZdKnqCcc86lUpUJStI2kkZLmippiqSzcm1YUj9JJumgxLjHJfWrZrmzJW2U63pqQtIrOcyT0/ol3SapY57imiGpWT7acs659UV1R1DLgP9nZh2B3sBvavilPAu4qIYxnQ0UJEGZ2R75Wr+ZnWRmU9c5KOecc1lVmaDMbI6ZvRmHFwHTgK1r0P7bwAJJAzInSNpb0luSJku6Q1IDSWcCWwGjJY3OsswMSX+MPbqOl7SLpGckfSjp1DjPxpJekPRmbPuQxPKL499+kl6S9JCkdyWNULDG+iXdFNc1RdJlibZektSjol1JV0l6W9I4SVvE8c0lPSzpjfjYM47fXNKzsc3bANVgmzrnXO1gZjk9gDbAp0CT+PxcYGKWxw1xej/gcaAv8J847vE4viEwE2gXx98FnB2HZwDNKolhBnBaHP4rMAloDDQH5sbx9RIxNgOms6rn4MWJ2BYArQhJ+lWgT7b1A5vFv3WBl4Cd4/OXgB5x2ICD4vCfgIvj8L2JdlsD0+LwDcAlcfhncfmsr9kf/vCHP2rro141+QsIRyXAwzGJLAQws2HAsOqWNbMxkpDUJzG6PfCxmb0fn/8T+A1wXQ7hjIp/JwMbWziyWyTpB0mbAN8Cf5DUF1hBOOLbAvg8o53XzWxWfH0TCQn4v1nWd5SkUwiJryXQkZAYk34kJF+ACUDFEeM+QEdp5QFSk7gt+wKHAZjZE5K+zuF1O+dcrVJtgpJUn5CcRpjZI4nx5wKDsiwyxszOzBh3FXAx4TetdfVD/LsiMVzxvF6MqTmwq5ktlTSDcMRWWTsAy8myLSS1Bc4BeprZ15KGV9LWUjOzLG3VAXqb2ZKMdit9cc4554LqruITcDvh1NS1yWlmNszMumV5ZCYnzOxZYFNg5zjqPaCNpB3i8/8F/hOHFxFO262tpsAXMTn1B7at4fLJ9TchHJEtiL8r7V/Dtp4Fzqh4IqlbHBwD/CKO25+wbZxzziVUdxXfnoTk8dN4YcJESQes5bquArYBiEcUvwQelDSZcPRzc5zvFuDpbBdJ5GgE0CO2exzwbg2XX7l+M3sbeCu2cS8wtoZtnRljmSRpKnBqHH8Z0FfSFMKpvk9r2K5zzq33tOrMlHPOOZceXknCOedcKnmCcs45l0qeoJxzzqVSWScoSfUkzZN0dcb4C2vQRpU19ZIVI9Yivr6xosUySUckxneT9GqsJDFJ0tFr034xSdokUXljmqTdSx1TZSTtJ+k9SdMlnV/qeKrisRaGx1oYRY+11HcKr8uDcNn3WOBD4gUfcfziHJevm8M8LxErRqxFfG0Il9bfBRyRGN8O2DEObwXMATYp9fas5rX8EzgpDm+Q1ngJFT8+BLaLcb4NdCx1XB6rx+qx1vxR1kdQwLHA9YTLtHcHiEdTG8ZL4kdkLhDr5v1F0tvA7hVHSJLqShou6Z1Yw29IxnJ14vQrcw3OzGaY2STCZfTJ8e+b2Qdx+DPgC8LNxakkqSmh+sXtAGb2o5l9U9KgKtcLmG5mH5nZj8D9wCHVLFMqHmtheKyFUfRYyzZBSWpIKCX0GHAfIVlhZucD31u4aThbpYtGwGtm1tXMkqWNugFbm1lnM+sC3JmYVo9wf9UHZnZxnl9HL8LeyIf5bDfP2gLzgDsVCvzeJqlRqYOqxNaEOo8VZlGzAsfF5LEWhsdaGEWPtWwTFHAgMNrMvieUYjpUUt0cllse58/0EbCdpL9J2g9YmJj2D+AdM7tqXYNOktQSuBv4pZmtqG7+EqoH7ALcZGbdCdU1Un2u3DlX/so5QR0L7BNr7U0ANgd+msNyS8xseeZIM/sa6Er4zelU4LbE5FeA/vGoLS8kNQGeAC4ys3H5ardAZgGzzOy1+PwhQsJKo9nEiiVRqzgujTzWwvBYC6PosZZlgopf7nsBrc2sjZm1IVRDPzbOsjQWua1Jm82AOmb2MKGwbfIL+HbgSeABSTlVgK9mXRsAjwJ3mdlD69peoZnZ58BMSe3jqL2BtHbW+Aawo6S2cTsfw6oK+GnjsRaGx1oYRY91nb9sS+TnwItmlqxI/m/gT5IaEOrpTZL0ZiW/Q2WzNeE3loqkfUFyopldGy8WuFvSoFxOyUnqSUhEmwIHSbrMzDoBRxEuOthc0uA4+2Azm5hjrKVwBjAifjA/ItRSTB0zWybpdOAZwlVHd5jZlBKHlZXHWhgea2GUIlavxeeccy6VyvIUn3POufWfJyjnnHOpVNIEJam+pKslfRBLAr2q0IFfKsUbdY+ofs6c2jo+vu4PJB2fjzYLRdIdkr6Q9E6pY6mOx1oYHmthlFOsUPx4S30EdQXQEuhsZrsAh7JuvenWWD6uyluLdW4GXArsRrg7+1JJae5VdziwX6mDyNFwPNZCGI7HWgjDKZ9YocjxlixBSdoIOBk4o+JqPDOba2YPxOkD4xHVm5IelLRxHD9D0mVx/GRJHeL4n2hVr79vSWqsYFiifNHRcd5+kl6WNAqYGsscDZP0hkLx1l/F+STpRoXiiM8DLfL08vcFnjOzr+L9V8+R4g+pmY0Bvip1HLnwWAvDYy2McooVih9vKY+gdgA+NbOFmRPiPUkXA/vEI6vxwG8Ts3wZx98EnBPHnQP8xsy6Ee6R+p7QnXo3wg24+wDDYvUGCPc5nWVm7YATgQVm1hPoCZwsqS3hcvb2QEdC9/F75Oell1V5E+ecK4m03gfVm5AUxkqCUKvu1cT0R+LfCYQkBKGq+bUKBWIfMbNZkvoA98XKEXMl/YeQgBYCr5vZx3HZgcDOid+XmgI7Eu5Vqlj+M0kvFuC1Ouecy6KUCWo60FpSkyxHUSKcAjs2y3IAFTfoLie+BjO7WtITwAGExLZvNev/NmN9Z5jZM6sFIR2Qw+tYG7OBfonnrQgllpxzzkUlO8VnZt8RSghdH6sTIKm5pCOBccCeknaI4xtJaldVe5K2N7PJZnYNoSRHB+Bl4Oj4G1NzwhHR61kWfwY4TbE8kqR2CtW6xySWbwn0z8NLr1jfQEmbxosjBsZxzjnnolJfxXcxoRuHqfGyxceBhWY2DxgM3CdpEuH0Xodq2jo7XgwxCVgKPEUoMzSJ0LHWi8B5sa5cptsIteXejHH8g3Bk9ijwQZx2F6ufZlxrZvYV4QrGN+Lj8jgulSTdR3jt7SXNknRiqWOqjMdaGB5rYZRTrFD8eL3UkXPOuVQq9RGUc845l5UnKOecc6nkCco551wqlV2CkjRU0uxE1YiJkjbJ8zouzGd7layjbGrxwcoKHpPj9h6fgnjWqAkm6UhJUyStkNSjlPEleayF4bEWRppiLbsEFf3VzLolHt/kuf2CJiiVXy2+Cv3j9k7DP9Nw1iwP9Q7hxu0xRY+masPxWAthOB5rIQwnJbGWa4Jag6Rxkjolnr8kqUe8h+oOSa8r1Og7JE4fLOkRSU/Ho5g/xfFXAxvGI4URcfknJL0dL2M/Og/hllUtvjTKVhPMzKaZ2XslCqlSHmtheKyFkaZYyzVBDUmc3hsdx40kdKVOvKm2pZmNBy4idA/fi3Cj7bB4Ey6EOn1HA10IN+RuY2bnA9/HI4VBhMTxmZl1NbPOwNN5iL8ca/EZ8KykCZJOKXUwzrn1X7kmqOQpvorqDg8AFbX0jgIeisMDgfMlTSSUE2oItI7TXjCzBWa2hHAz7rZZ1jUZGCDpGkl7mdmC/L+cstAnFujdH/iNpL6lDsg5t34r1wS1BjObDcyXtDPhqGhknCTg8ERCa21m0+K0HxJNrKzrl9Hu+4TK55OBKyVdkodwZwPbJJ63iuNSK25fzOwLQoWNXqWNyDm3vltvElQ0EjgPaGpmk+K4Z4AzpFAWXVL3HNpZqlV1+bYCvjOze4BhhGS1rsqqFl/8Ha5xxTAh3rLoAdQ5V8bMrKwewFDC0cbExKNNnLYFsAy4NDH/hoTaepOBKcDjcfxg4MbEfI8D/eLwNcA0YAThgoZJcT1vAD3y9DpOIFR0nw78stTbtZpYtyPUM3w7bsOLUhDTfcAcQt3FWYQ+vX4eh38A5gLPlDpOj9Vj9VjX/uG1+JxzzqXS+naKzznn3HrCE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulf4//r4tcV5c554AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "stride = 16\n", "agg_d1, agg_n1 = aggregate(d1, n1, stride)\n", "agg_d2, agg_n2 = aggregate(d2, n2, stride)\n", "agg_d1, agg_n1, agg_d2, agg_n2 = [list(map(int, await mpc.output(_))) for _ in (agg_d1, agg_n1, agg_d2, agg_n2)]\n", "T1_, E1_ = events_from_table(agg_d1, agg_n1)\n", "T2_, E2_ = events_from_table(agg_d2, agg_n2)\n", "T1_, T2_ = [t * stride for t in T1_], [t * stride for t in T2_]\n", "kmf1_ = KaplanMeierFitter(label='1=Maintained').fit(T1_, E1_)\n", "kmf2_ = KaplanMeierFitter(label='2=Not maintained').fit(T2_, E2_)\n", "plot_fits(kmf1_, kmf2_, 'Aggregated Kaplan-Meier survival curves', 'weeks')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Picking `stride = 16` achieves a reasonable balance between privacy and utility. To enhance both privacy and utility at the same time, one may look for differentially private randomization techniques, adding a suitable type of noise to the Kaplan-Meier curves. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Secure Logrank Tests\n", "\n", "The function `logrank_test()` performs a secure logrank test on a secret-shared dataset, similar to function `lifelines.statistics.logrank_test()` used above for a dataset in the clear. The input parameter `secfxp` specifies the secure type to be used for fixed-point arithmetic, and the output is an instance of `lifelines.statistics.StatisticalResult`:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.06533928754794803\n" ] } ], "source": [ "print((await logrank_test(secfxp, d1, d2, n1, n2)).p_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Relying solely on p-values is in general not a good idea, and this is especially true when handling otherwise (mostly) hidden data. Together with the aggregated curves, however, the p-value may lead to a useful conclusion for a study. \n", "\n", "The function `logrank_test()` uses one secure fixed-point division per time moment in $1..maxT$. Even though these divisions can all be done in parallel, the total effort is significant when $maxT$ is large. However, \"most of the time\" there is actually no event happening and no divisions need to be performed for these time moments. E.g., in the survival tables for the aml dataset above, there are only 7 time moments with nonzero `d1` entries on the entire timeline $1..161$, and only 9 time moments with nonzero `d2` entries.\n", "\n", "Therefore, it may be advantageous to first extract the nonzero rows of the survival tables, and then limit the computation of the logrank test to those rows. The extraction needs to be done obliviously, not leaking any information about (the location of) the nonzero entries of the survival tables. To prevent this oblivious extraction step from becoming a bottleneck, however, we will actually exploit the fact that the aggregate curves are revealed anyway. We may simply use `agg_d1` and `agg_d2` to bound the number of events per stride, and extract the nonzero rows obliviously and efficiently for each stride. \n", "This is basically what the function `agg_logrank_test()` does:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0653392245841535\n" ] } ], "source": [ "print((await agg_logrank_test(secfxp, d1, d2, n1, n2, agg_d1, agg_d2, stride)).p_value)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even for a small dataset like aml, the speedup is already noticeable. For larger datasets, the speedup gets really substantial, as can be noticed for some of the other datasets included with the [kmsurival.py](kmsurvival.py) demo. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "We end with two complete runs of the demo on the aml dataset, showing the Chi2 test statistic and p-value for each logrank test. \n", "\n", "The help message included with the demo shows the command line options:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Showing help message for kmsurvival.py, if available:\n", "\n", "usage: kmsurvival.py [-h] [-i I] [-s S] [-a A] [--collapse] [--print-tables]\n", " [--plot-curves]\n", "\n", "options:\n", " -h, --help show this help message and exit\n", " -i I, --dataset I dataset 0=btrial(default) 1=waltons 2=aml 3=lung 4=dd\n", " 5=stanford_heart_transplants 6=kidney_transplant\n", " -s S, --stride S interval length for aggregated events\n", " -a A, --accuracy A number of fractional bits\n", " --collapse days->weeks->month->years\n", " --print-tables print survival tables\n", " --plot-curves plot survival curves\n" ] } ], "source": [ "!python kmsurvival.py -h" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Complete Run: 5 logrank tests + survival curves\n", "To show the plots of the survival curves the `main()` function of the demo is called directly from a notebook cell:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using secure fixed-point numbers: SecFxp64:32\n", "Dataset: aml, with 1-party split, time 1 to 161 (stride 16) weeks\n", "Chi2=3.396389, p=0.065339 for all events in the clear\n", "Chi2=3.396389, p=0.065339 for own events in the clear\n", "Chi2=2.685357, p=0.101275 for aggregated events in the clear\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD+CAYAAABvEpGeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9yElEQVR4nO2deZgU1dWH3x87CoIKGmQRNAqyi4AoEVEUdzSGuCaCRJMYRU2ihqiJW0w0ms3oZ9zRBFGDS4i7RohEURFE1rijggQRAUEWWc73x70NzdDN9PT0TNXAeZ9nnqm6devWr6q769S9de45MjMcx3EcJ23USlqA4ziO4+TCDZTjOI6TStxAOY7jOKnEDZTjOI6TStxAOY7jOKnEDZTjOI6TStxAOds0koZK+k/SOgAkzZTUvwTtzJF0WJH7ni7p2SL3bSvJJNWJ609JGlLAfgdJemsL20dK+lUxmrakz0k/bqC2IeKNa6Wk5ZIWxB9+oyLbGi/prEpoaStpnKQVkv5b7A21qpH0vahvWbxmT0pqXBXHMrNOZja+KtqugIZRZjawRG0dZWb3FlBvgpm1L8UxaxppekBKI26gtj2OM7NGQA+gJ3B5RXZWoBTfm9HAG8DOwGXAGEnNS9BuyZB0MPBr4FQzawzsAzxYZFv+1O44FcQN1DaKmc0DngI6S9pR0uOSFkpaHJdbZerG3tK1kl4CVgB/BQ4Cbo69sZsl3SLpd9nHkDRW0o/LHlvS3gQDeYWZrTSzh4HpwLcK0S7pGElvSPpC0seSrszalhnGOTNuWyzph5J6SZomaYmkmwu8TL2AiWb2Rrxmn5vZvWa2LOu6bOhFln0ajjrOlfQO8I6kWyXdWOZc/iHpJ3F5jqTDJO0We7o7ZdXbV9JnkupK2lPSC5IWxbJRkpoWeE5bJM85/FDSO/Ha3SJJcVttSTdGDe8Dx5Rpa7yksyTVj/t2ztrWPJ7jLpL6S5pb5lynxF7rg0CDfPqyNH49Luf9bhRw7rtJejj+Dj6QdH5Wed7PI64PkzQ7ft+ekbR7eddQ0j7AX4AD4u9oSax/tKRZ8fznSbqo0HPY2nADtY0iqTVwNKEXUwu4B9gdaAOsBMrexL8LfB9oDAwFJgDnmVkjMzsPuBc4VbF3JakZcBhwf47DdwLez9zoI2/GciR9I/NjzcOXwBlAU8JN8RxJJ5Spsz+wF3Ay8EdCL+2weIyTFHpH5fEqcISkqyT1lVS/gH3KckLU0pHQazw56wa/IzAQeCB7BzP7BJjIpgb7NGCMma0BBPwG2I3Qq2sNXFmEtkI5lmCsuwInAUfE8rPjtn0JvfHBuXY2s9XAI8CpWcUnAf82s0+z60qqBzxGeAjaCfg7BT64RAr5bmxG/N7+k/A9bAkMAC6UdER5n4ek44FLgROB5oTfxugyh9jsGprZbOCHhIegRmbWNNa9C/hB7LV3Bl6owPlvVbiB2vZ4LN78/wP8G/i1mS0ys4fNbEU0GtcCZW/gI81sppmtjTfJTTCz14ClhB82wCnAeDNbkENDo1g3m6UE44eZ/Sfrx7oZZjbezKab2Xozm0a4GZTVe42ZrTKzZwk3rdFm9mnsOU4g3FS3iJlNINx0egBPAIsk/V5S7fL2zeI3see1Mh7XCL1PCDf0ifEGWJb7iTf0aNBOiWWY2btm9pyZrTazhcDvc5x/KbnOzJaY2UfAOKB7LD8J+KOZfWxmnxOMZj7uJ5xDhtPI/fDSB6gb211jZmOASYUKLfC7kYteQHMzu9rMvjKz94E7sjTn/TwIRuY3ZjbbzNYShoW7Z/eiyH8Nc7EG6ChpBzNbbGZTCjn3rRE3UNseJ5hZUzPb3cx+ZGYrJW0n6TZJH0r6AngRaFrmRvxxAW3fC3wnLn+H8BSci+XADmXKdgCW5ai7GZL2V3CwWChpKeEG0axMtWzDuDLHekHOIWb2lJkdR3iaP57Qe6yIc8iG62YhMvMDbOxJnAaMyrPfw4ShnxZAP2A9wcAhaVdJD8Thny+Av7H5+W+GpDZxKGm5pOUVOIf/ZS2vYOO1241NvxcfbqGNccB28bNrS7hBP5qj3m7APNs0ivWW2t2EAr8budgd2C0OwS2JD3GXArvG7Xk/j7jvn7L2+5zQy22Z1X6+a5iLbxFGNz6U9G9JBxSgf6vEDZQD8FOgPbC/me1A+AFC+JFlKBv2PlcY/L8Bx0vqRhh6eizP8WYCe2hTb7husbwQ7gfGAq3NrAlhHF9b3qVyxCfyfxGGWzLvUr4Etsuq9rVcu5ZZHw0Mjk/X+xNufLmOtxh4ljBEeRrwQNZN+9ex3S7x8/oOBZy/mX0Uh5IaRUeZyjKfMLyYoc0Wjr0OeIhgnE8FHi8zxJvdZsvMMGiOdje55pLKXvNivxsfAx/Eh7fMX2MzOzrq39Ln8TFhSC5734Zm9nIBx93sd2Rmk8zseGAXwm/ooQLa2SpxA+VAGFpbCSyJL4KvKGCfBcAe2QVmNpcwHPNX4OE4rLUZZvY2MBW4QlIDSd8kjM3nvFnn0fu5ma2S1Jtwwyg5ko6XdIqCE4nisQ4GXolVpgInxh7o14HvlddmdLj4DLgTeMbMlmyh+v2E9ymD2XQ4rDGhF7pUUkvg4oqdWcl4CDhfUqv4Pm1EOfXvJ9zgTyf38B6Edz1rY7t1JZ0I9M7a/ibQSVJ3SQ3Y/N1bsd+N14Blkn4mqaGCA0hnSb3K6M/1efwF+LmkzDvUJpK+XeBxFwCt4rs3JNVTmIvWJA6lf0HorW2TuIFyIDgRNCTcOF8Bni5gnz8RegKLJd2UVX4v0IX8w3sZTiG8WF8MXAcMju9TMhM3tzQE9SPgaknLgF9SdU+YiwmOAO8QbhR/A24ws8yw3B+Arwg3mXvJP1xXlvvJ70CSzViCo8f/zOzNrPKrCO/FlhLejT1S4HFLzR3AMwSjMaU8HWb2KqEHtBvBgzRXna8I7/2GEobKTs5uNz7cXA08T/hcys4hKuq7EXt4xxKGHj9g40NEk6xqOT8PM3sUuB54IA65zgCOKuS4hB75TOB/kj6LZd8F5sS2fkgw6NskMk9Y6JQQSf0IN/Ldzb9cjuNUAu9BOSVDYU7IBcCdbpwcx6ksbqCckhAnHS4BWhCGDB3HcSqFD/E5juM4qcR7UI7jOE4qcQPlOI7jpJIaHWFZ0t0E19BPzaxzju0iuEMfTZi9PbSQsCHNmjWztm3bllit4ziOU5bJkyd/ZmY5MxnUaAMFjCQENb0vz/ajCPMW9iLM2r81/t8ibdu25fXXXy+RRMdxHCcfkvKGsqrRBsrMXoxxvfJxPHBfdHl+RVJTSS3MbH5V6Hnl/86m8ZLZBdV9qeEh/Gu7o/NuP757S07bP2/kGMdxnK2eGm2gCqAlmwaznBvLNjNQkr5PSCdBmzbFG4b1a9eWW6fd+g9Zv3Ytj9Y6POf2Dxet4ItVa9xAOY6zTbO1G6iCMbPbgdsBevbsWZTvfZ8f3VFYxXuOocuqJTx1Tr+cm0++bSKLV3zF2Knz8jbRuGFdDmm/SzEyHcdxagRbu4Gax6bRllvFstSzfr3RvHGDvNsXLltVjWocJ3nWrFnD3LlzWbXKv/s1kQYNGtCqVSvq1q1b8D5bu4EaC5wn6QGCc8TSqnr/5DhO1TJ37lwaN25M27Zt2TQbh5N2zIxFixYxd+5c2rVrV/B+NdpASRoN9AeaSZpLSBNRF8DM/gI8SXAxf5fgZn5mMkpzYOtg+pjNyxs0oZBceqvWrtviEGAGHwp0thZWrVrlxqmGIomdd96ZhQsXVmi/Gm2gzOzUcrYbcG41yakY69dDo103L1++gEIMVOsdty/oMD4U6GxNuHGquRTz2XkkCcdxnAKRxHe+850N62vXrqV58+Yce+yxW9zv9ddf5/zzz99inSVLlvB///d/Bek48MADC6pXHnPmzKFz581iHKSGGt2D2pqZt2QlVz+ePwN63z2bMWCfHD0wx3GqjO23354ZM2awcuVKGjZsyHPPPUfLli3L3a9nz5707Nlzi3UyBupHP/pRue29/HIh2eRrPt6DSiHHd29Jy6YN827/cNEKXnrvs7zbs8m8qyrvb9xbn5ZKvuNs1Rx99NE88cQTAIwePZpTT934puG1117jgAMOYN999+XAAw/krbfeAmD8+PEbellXXnklw4YNo3///uyxxx7cdFNISD1ixAjee+89unfvzsUXX8zy5csZMGAAPXr0oEuXLvzjH//YcJxGjRptaLd///4MHjyYDh06cPrpp5PJUDF58mQOPvhg9ttvP4444gjmz5+/obxbt25069aNW265pYqvVuXwHlQKOW3/NjSqXzuvm/nVj89k4bLVBbXl76qcrZGr/jmTWZ98UdI2O+62A1cc16nceqeccgpXX301xx57LNOmTWPYsGFMmDABgA4dOjBhwgTq1KnD888/z6WXXsrDDz+8WRv//e9/GTduHMuWLaN9+/acc845XHfddcyYMYOpU6cCYfjw0UcfZYcdduCzzz6jT58+DBo0aLN3OW+88QYzZ85kt912o2/fvrz00kvsv//+DB8+nH/84x80b96cBx98kMsuu4y7776bM888k5tvvpl+/fpx8cUXV/7CVSFuoGoony3/KmkJjrNN0rVrV+bMmcPo0aM5+uhNw5UtXbqUIUOG8M477yCJNWvW5GzjmGOOoX79+tSvX59ddtmFBQsWbFbHzLj00kt58cUXqVWrFvPmzWPBggV87Wtf26Re7969adWqFQDdu3dnzpw5NG3alBkzZnD44SFazbp162jRogVLlixhyZIl9OsXggR897vf5amnnqr0Nakq3ECljTWrYPoYdvt4MTs03HxC29p6OwDNql+X46SIQno6VcmgQYO46KKLGD9+PIsWLdpQ/otf/IJDDjmERx99lDlz5tC/f/+c+9evX3/Dcu3atVmbI0TaqFGjWLhwIZMnT6Zu3bq0bds25yTlXG2ZGZ06dWLixImb1F2yZEkFzzRZ3EAlxdKP4OkRm5fv0R/2PpLVDWqzpmG9zTbXXbkQN1COkyzDhg2jadOmdOnShfHjx28oX7p06QaniZEjR1aozcaNG7Ns2bJN2tpll12oW7cu48aN48MP8wb93oz27duzcOFCJk6cyAEHHMCaNWt4++236dSpE02bNuU///kP3/jGNxg1alSFNFY3bqCSoMtg+DLHhLXPPwj/9z6yevVQ+MTfXPhkYGdbo1WrVjndxi+55BKGDBnCr371K4455pgKtbnzzjvTt29fOnfuzFFHHcXPfvYzjjvuOLp06ULPnj3p0KFDwW3Vq1ePMWPGcP7557N06VLWrl3LhRdeSKdOnbjnnnsYNmwYkhg4cGCFNFY3ynh8OBvp2bOnVXk+qOljNp+om+lRHXkdL737GU23y92DumDmHsyev4x9WjQu9zDV4Y6+cNkqBnUv39XWcSrD7Nmz2WeffZKW4VSCXJ+hpMlmltMH33tQKWfUW3D/29klzYEwDDB7/rLN6jdrVI/mjcOY9IeLVgCf+Xwpx3FqJG6gUs7p7cNfhvrLP2Z97QYMfL45zx625bhWF73eBFauZOcP/lnucdbW24GlLQ+urFzHcZyS4QaqhrG60cbsIWsaNt9iXatdWD3IOF84juOkB48k4TiO46QSN1CO4zhOKvEhPqfS5HNRd/dzx3Eqg/egaiin7Z20go203nF7mjdusNnfspW5w7w4Tk1FEj/96U83rN94441ceeWVW9znscceY9asWVWi56yzziq37UKP/5e//IX77ruvJLqGDh3KmDE5ErJWEO9B1VCyPfu2xPtLYUQBkfkHNKvN0eT39nMvP8cJYYUeeeQRfv7zn9OsWWERXR577DGOPfZYOnbsWHI9d955Z8mO/8Mf/rBUskqG96C2Yg5uCXs0Kb/e+0vhX5/txJqGzfP+1fmqtJGjHacmUqdOHb7//e/zhz/8YbNtc+bM4dBDD6Vr164MGDCAjz76iJdffpmxY8dy8cUX0717d957771N9hk6dCjnnHMOffr0YY899mD8+PEMGzaMffbZh6FDh26od84559CzZ086derEFVdcsaG8f//+ZIIKNGrUiMsuu4xu3brRp08fFixYkPP4d9xxB7169aJbt25861vfYsWKFUBIA3LjjTduaPdnP/sZvXv3Zu+9994QrX3dunVcfPHF9OrVi65du3LbbbcBIbDteeedR/v27TnssMP49NPSpO/xHtRWzFG7h7/yKKSH5Tip4qkR8L/ppW3za13gqOvKrXbuuefStWtXLrnkkk3Khw8fzpAhQxgyZAh33303559/Po899hiDBg3i2GOPZfDgwTnbW7x4MRMnTmTs2LEMGjSIl156iTvvvJNevXoxdepUunfvzrXXXstOO+3EunXrGDBgANOmTaNr166btPPll1/Sp08frr32Wi655BLuuOMOLr/88s2O37RpU84++2wALr/8cu666y6GDx++ma61a9fy2muv8eSTT3LVVVfx/PPPc9ddd9GkSRMmTZrE6tWr6du3LwMHDuSNN97grbfeYtasWSxYsICOHTsybNiwgi77lnADlTY+/wCeHkHnlWuoU3vzDu7Srx3I4lYDEhAWGDP5Ywbv17r8io6zlbLDDjtwxhlncNNNN9Gw4cbEohMnTuSRRx4BQhqLsgYsH8cddxyS6NKlC7vuuitdunQBoFOnTsyZM4fu3bvz0EMPcfvtt7N27Vrmz5/PrFmzNjNQ9erV25AUcb/99uO5557LebwZM2Zw+eWXs2TJEpYvX84RRxyRs96JJ564oa05c+YA8OyzzzJt2rQN75eWLl3KO++8w4svvsipp55K7dq12W233Tj00EMLOvfycAOVJvbov8XNDZaFaMZJGqiHp8xzA+UkTwE9narkwgsvpEePHpx55pmVbiuTLqNWrVqbpM6oVasWa9eu5YMPPuDGG29k0qRJ7LjjjgwdOjRn2o26detuSGaYL4UHhGHFxx57jG7dujFy5MhNorHn0pXdlpnx5z//eTOj9uSTT1bspAvEDVSa2PvIDZHMZ+QIFtv29WvKbaL5e2NYuGfuoYTqpjIR0ovB3dqd6mKnnXbipJNO4q677towlHXggQfywAMP8N3vfpdRo0Zx0EEHAZun0agoX3zxBdtvvz1NmjRhwYIFPPXUU3nzTOWi7PGXLVtGixYtWLNmDaNGjdqQHqQQjjjiCG699VYOPfRQ6taty9tvv03Lli3p168ft912G0OGDOHTTz9l3LhxnHbaaRU5zZy4gdrK2OX9R6rEQNVatyrG9Gu+IbZfeZ59haabLxWett6pTn76059y8803b1j/85//zJlnnskNN9xA8+bNueeee4CQIv7ss8/mpptuYsyYMey5554VOk63bt3Yd9996dChA61bt6Zv374V2r/s8a+55hr2339/mjdvzv77718h43nWWWcxZ84cevTogZnRvHlzHnvsMb75zW/ywgsv0LFjR9q0acMBBxxQIY358HQbOUgs3UYWudJtZHpQc3r+Iu9+nZ47jZmH318hKRkniesOLL/uMf+EJ44Ly3VXLmRRu+MqdKyqxNN+bN14uo2aj6fbcIqi0PlSsLGe1jVhzcyZ1ZJzynGcbQ83UA4H5+h0LFgBn67MXX/6osxSPViyjNnzl7F4xVfuPOE4TklxA1XDaLDsw3KdJQpxpsjmHOC0NoW5r5cd4rtg5h4Abpwcxyk5bqBqEEu/tvElUd2VC6m36rOc9bZfPHuzsq8aNMubFyoN7uuOUwhmtsGV2qlZFOPv4AaqBrG41YByjUgxThIV7XE5ThI0aNCARYsWsfPOO7uRqmGYGYsWLaJBgwYV2s8NlFM0tdatou7K8EIqO628B5Z1qoJWrVoxd+5cFi707M81kQYNGtCqVasK7VOjDZSkI4E/AbWBO83sujLbhwI3AJnZojebWfnhf52CWN2odc608p4+3qkK6tatS7t27ZKW4VQjNdZASaoN3AIcDswFJkkaa2ZlE588aGbnVbvArZQ05aFyHGfrpsYaKKA38K6ZvQ8g6QHgeKBqMoM5QOF5qJKgukMrOY5TtSHGarKBagl8nLU+F9g/R71vSeoHvA382Mw+zlEHSd8Hvg/Qpk2bEkutONvVr8OSFV/l3f7VuvXs0njzF46f7nFiVcrKSdlJvklN4K3u0EqO41RtiLGabKAK4Z/AaDNbLekHwL1AzjjwZnY7cDuEUEfVJzE3+7ZuusXtL72b28W82Dh85c2vypfmI9ckX4APF60APvMIE47jFE1NNlDzgOzZoa3Y6AwBgJktylq9E/htNegqjAZNYPmCzcvXrIIdC8gyWEKy51flYkvzpHIlRay7cikXzNy5ZPocx9k2qckGahKwl6R2BMN0CrBJfHdJLcxsflwdBGw+gzUp9jo8d/n0MdWrg/LnV1V0nlQ+9/Nicbd1x9k2qbEGyszWSjoPeIbgZn63mc2UdDXwupmNBc6XNAhYC3wODE1M8DZEPvfzYnG3dcfZNqmxBgrAzJ4EnixT9sus5Z8DP69uXdVBeU4UGfI5U9RkPO2842wb1GgDtS1TnhNFhnzOFDUZTzvvONsGtZIW4DiO4zi58B6UUxC53NDzuZ47juOUAjdQTrnkckP3FB2O41Q1bqCccsnlhl6I63lF0shviUxkimyufnzmZvU89bzjbF24gXKqhHwRJirCxrTzIbV8NrPnL9us/sJlq91AOc5WhBuorZyk3NFzRZgolvrLP2Z97Y3aBj7fnNFn99mkTq4eleM4NRs3UGkjXwikDBUMhbQ1uKOvbuQu5Y6zLeIGKm3kC4GUIYFQSI7jOEng86Acx3GcVOI9KKdo8qXo8PlRjuOUAjdQTlHkS9FRHfOjvtPuyypr23Gc9OAGqqZRnhNFkWz35XzYbo+C6+dL0VHR1BzFcMaeK1hUfjXHcWo4bqBqGuU5URRJvU/uLsgdPRdbY8R0x3EK48np8xnUvQQTH3PgBsoBoNNuO0CjZkXtm2YXdcdxqpanZ5Z+RCeDe/E5juM4qcQNlOM4jpNKfIjPKTn53M+Lxd3WHWfbxA2UEyhRiKV87ufFUhG39Q8XrfCYfI6zFeEGygmUKMRSPvfzYim0J9Z3z2aAO2s4TlWxcNlqPlue29O37YgnNiu7YMBe/PjwvSt1TDdQTqXJFzG9Ot3PB+yzq6facJwEOPWOV5hz3TFV0rYbKKfS5IuY7u7njuNUBvficxzHcVKJGyjHcRwnlfgQn1MYubz8Kpg80XEcpyK4gXIKI5eXXzUlTyw7r2r5zl1Z1O64ajm24zhb5shOVeec5AbKqTLyefdVhPk77U/zdeth3XoAtv/yI9avWc3CZas2q7tq7Tpa77h9pY7nOE7FOLpLiypr2w2UU2Xk8+6rEF8fDAzeuP70CBqvW5MzevLYqfMqfzzHcVKDGyinePJFn/B3U47jlAA3UE7x5Is+UU3vphzH2bqp0W7mko6U9JakdyWNyLG9vqQH4/ZXJbVNQKbjOI5TBDW2ByWpNnALcDgwF5gkaayZzcqq9j1gsZl9XdIpwPXAydWvdhuj1GnpfcjQcbZJaqyBAnoD75rZ+wCSHgCOB7IN1PHAlXF5DHCzJJmZVafQbY5Sp6X3IUPH2SapyQaqJfBx1vpcYP98dcxsraSlwM7kCHst6fvA9wHatGlTFXqdYsnukTX+GtSun7Na44Z1c7qfO45TdTRuWLfK2q7JBqqkmNntwO0APXv29B5WmsjukXUZnLfaIe13qQYxjuNUFzXZSWIe0DprvVUsy1lHUh2gCbCoWtQ5juM4laImG6hJwF6S2kmqB5wCjC1TZywwJC4PBl7w90+O4zg1A9Xk+7Wko4E/ArWBu83sWklXA6+b2VhJDYC/AvsCnwOnZJwqyml3IfBhkbLSnto1zfpcW/GkWZ9rK5406yuVtt3NrHmuDTXaQKURSa+bWc+kdeQjzfpcW/GkWZ9rK54066sObTV5iM9xHMfZinED5TiO46QSN1Cl5/akBZRDmvW5tuJJsz7XVjxp1lfl2vwdlOM4jpNKvAflOI7jpBI3UI7jOE4qcQNVQspL/1HNWlpLGidplqSZki6I5TtJek7SO/H/jglqrC3pDUmPx/V2MS3KuzFNSr0EtTWVNEbSfyXNlnRAWq6dpB/Hz3SGpNGSGiR57STdLelTSTOyynJeKwVuijqnSeqRgLYb4uc6TdKjkppmbft51PaWpCOqUls+fVnbfirJJDWL64lfu1g+PF6/mZJ+m1Ve+mtnZv5Xgj/CZOH3gD2AesCbQMcE9bQAesTlxsDbQEfgt8CIWD4CuD5BjT8B7gcej+sPESZTA/wFOCdBbfcCZ8XlekDTNFw7QgDkD4CGWddsaJLXDugH9ABmZJXlvFbA0cBTgIA+wKsJaBsI1InL12dp6xh/t/WBdvH3XLu69cXy1sAzhIABzVJ07Q4Bngfqx/VdqvLaeQ+qdGxI/2FmXwGZ9B+JYGbzzWxKXF4GzCbc3I4n3HyJ/09IQp+kVsAxwJ1xXcChhLQoSWtrQvhx3gVgZl+Z2RJScu0IQZ4bxviS2wHzSfDamdmLhEgt2eS7VscD91ngFaCppBbVqc3MnjWztXH1FUIcz4y2B8xstZl9ALxL+F1XGXmuHcAfgEuAbC+2xK8dcA5wnZmtjnU+zdJW8mvnBqp05Er/0TIhLZugkEl4X+BVYFczmx83/Q/YNSFZfyT8ANfH9Z2BJVk3jiSvXztgIXBPHIK8U9L2pODamdk84EbgI4JhWgpMJj3XLkO+a5W238kwQq8EUqJN0vHAPDN7s8ymNOjbGzgoDif/W1KvqtTmBmorR1Ij4GHgQjP7Inubhb55tc8zkHQs8KmZTa7uYxdIHcLQxq1mti/wJWGYagMJXrsdCU+r7YDdgO2BI6tbR0VI6lqVh6TLgLXAqKS1ZJC0HXAp8MukteShDrATYYjxYuChOPpRJbiBKh2FpP+oViTVJRinUWb2SCxekBkWiP8/zbd/FdIXGCRpDmEo9FDgT4Qhi0yOsiSv31xgrpm9GtfHEAxWGq7dYcAHZrbQzNYAjxCuZ1quXYZ81yoVvxNJQ4FjgdOjAYV0aNuT8PDxZvx9tAKmSPpaSvTNBR6Jw4yvEUZAmlWVNjdQpaOQ9B/VRnyquQuYbWa/z9qUnYJkCPCP6tZmZj83s1Zm1pZwnV4ws9OBcYS0KIlpi/r+B3wsqX0sGgDMIgXXjjC010fSdvEzzmhLxbXLIt+1GgucET3S+gBLs4YCqwVJRxKGlweZ2YqsTWOBUyTVl9QO2At4rTq1mdl0M9vFzNrG38dcgrPT/0jBtQMeIzhKIGlvggPRZ1TVtatKL5Bt7Y/gZfM2wYPlsoS1fIMwrDINmBr/jia86/kX8A7BG2enhHX2Z6MX3x7xS/0u8Heip1BCuroDr8fr9xiwY1quHXAV8F9gBiGdTP0krx0wmvA+bA3hhvq9fNeK4IF2S/yNTAd6JqDtXcL7kszv4i9Z9S+L2t4Cjkri2pXZPoeNXnxpuHb1gL/F794U4NCqvHYe6shxHMdJJT7E5ziO46QSN1CO4zhOKnED5TiO46QSN1CO4zhOKnED5TiO46QSN1CO4zhOKnED5TjloJB640dZ67tJGrOlfSrQ9pWSLorLV0s6rETttlBMY1IVSFpegbrPK8G0Lk7NxQ2U45RPU2CDgTKzT8xscP7qxWFmvzSz50vU3E+AO0rUVmX5K1nXz3EKJTEDFRNx/S1rvY6kheU99UnqKemmcups8sRbTt2XC1NcbjttcyUdc7YKrgP2lDQ1Jrvb8FlLGirpMYWkfHMknSfpJzEK+iuSdor19pT0tKTJkiZI6lD2IJJGShocl+dIukrSFEnTM/UlbR8Tyb0Wj5Evpcu3gKfjPk9I6hqX35D0y7h8taSz4/LFkiYpJMK7KkvTd+Kxpkq6TVLtMpqbSZoo6ZjYa3sx1p0h6aBYbSxwapHX3tmGSbIH9SXQWVLDuH44BQQXNLPXzez8cqo1pcAnNjM7sJB6zjbNCOA9M+tuZhfn2N4ZOBHoBVwLrLAQBX0icEasczsw3Mz2Ay4C/q+A435mZj2AW+M+EMLJvGBmvQkx0W5QSAWygRgLbbHFnD3ABEKKhCaE6N19Y/lBwIuSBhJip/UmhHjaT1I/SfsAJwN9zaw7sA44Pes4uwJPAL80syeA04BnYt1uhDBCmNlioL6knQs4Z8fZQNJDfE8SktZBeMIandkgqXd8MntD0suZwJ2S+mtjivAr49PkeEnvS8oYrrJPvI0k/SvrafT4rOMsz2p3vDam+R4lhTDykvZTyH0yWdIz2hileT9Jb0p6Ezi3ai+Vk2LGmdkyM1tIyM/0z1g+HWirkPLkQODvkqYCtxEyHpdHJgL9ZKBtXB4IjIjtjAcaAG3K7NeCkM8qwwRCAsa+BIPSSCGtQzszeyu2ORB4gxBfrQPBYA0A9gMmxeMNIMT8A6hLiLV3iZk9F8smAWdKuhLoYiFRZoZPCelBHKdg6pRfpUp5APhlNDhdgbsJT3UQgmEeZGZr44vjXxOGLcrSgfAk2Rh4S9KthCfezvFJDoU0BN80sy8kNQNekTTWNg9EuC/QCfgEeAnoK+lV4M/A8Wa2UNLJhKfkYcA9wHlm9qKkG0pxQZwayeqs5fVZ6+sJv7FahISC3Ytsdx0bf6sCvhUNSz5WEgxXhklAT+B94DlCeoSzCYYv0+ZvzOy27EYkDQfuNbOf5zjG2rj/EcC/IWRgldSP8NA5UtLvzey+WL9B1OU4BZNoD8rMphGeDE8l9KayaUJ44pxBSH/cKU8zT1hIM/wZ4SktV5ZTAb+WNI0QWbllnnqvmdlcM1tPGJ5oC7QnDOE8F58iLwdaSWoKNLWQFhnCi2Bn62QZ4QGoKCwkivxA0rchpEKR1K3I5p4Bhmf17vfNUedtNva4MLOvCNG7v00YdpxAGDLMfHefAYbFnh6SWkrahdBDGhyXkbSTpN0zzRIe0jpI+lncvjuwwMzuAO4k5NDKpH75GiEyt+MUTNI9KAgvUG8kpF3IHqO+hjB08k2FlOXj8+yf/fSa/aSZzelAc2A/M1ujkAisQY56udoSMNPMDsiuGA2Usw1gZoskvRQflp4ipDyoKKcDt0q6nDA89gBQNqV3IVwD/BGYJqkW8AEh8V623i8lvSfp62b2biyeAAwws5WSJhASyk2I9Z+N75smRru3HPiOmc2Kep+Nx1pDGMr+MO63TtKpwFhJywjvlS+WtCa2kXn/th/wim1MSe84BZEGA3U3YfhjuqT+WeVN2Og0MbSCbZZ94m1CSDG+RtIhwO65d8vJW0BzSQeY2USFLLV7m9lMSUskfcPM/kPWy2Nn68PMTitT1DmWjwRGZtVrm7W8YZuZfUCO1OxmdmXW8tA87bxOeIDDzFYCPyhA8s2E383lcb9fAL+Iy58QHryydfyJkNW4rL4HgQdzlDeK/1cThvky3JtDy3cpzCnEcTYhaScJ4pBaLrfx3wK/kfQGFTSkZrYIeCm6ut4AjAJ6SppOeKr7bwXa+oqQqfT66AwxlfDCG+BM4JY49KecDThOApjZo6RnSG2Gmf0raRFOzcMTFjqO4zipJPEelOM4juPkwg2U4ziOk0rcQDmO4zipJHEDpY0x+K4rU35pBdq4U1LHLWwfL6lnkfr6xQgUaxXjpGVtezp68lVZ1OhikdQ+RtLI/H0h6cKkdWWQ1EAhxtubkmYqK/5b0sToJJ8qpbEV06zPtRVPmvUlpS1xA0WIwfc28O3M5MNIQQZKUm0zO8vMZlWJOviI4K57f45tNxBcaFOHmb0VY8d1J8xDWQE8mqyqTVgNHGpm3Qjx346U1CdZSRsYSQ6X8BQxkvTqG4lrK5aRpFffSBLQlgYDdSph/sVHwAEAsTfVMD75jyq7g6Tlkn4X3b4PyPSQJNVWiAg9QyHm3o/L7Fcrbv9VoeLMbE6MeLE+x7Z/EeZcpZ0BhGCnHyYtJIMFMjmF6sa/VLiUxuggnyetIx9p1ufaiifN+pLSluhEXUkNgMMIEw+bEozVy2Y2QtJ5W4hdtj3wqpn9NLaTKe8OtDSzzrG8adY+dQjzoWaY2bUlPZH0cwpZgXjTgkLqhsnA14FbzOzVhCU5jpMiku5BHUsIZ7QSeBg4QWXyzeRhXaxflveBPST9WdKRwBdZ225jGzROkuoBg4C/J62lLGa2Lj6EtAJ6S+qcsCTHcVJE0gbqVOCwGBtvMiEW36EF7LfKzNaVLYx5Z7oR4vb9kBCwMsPLwCGx17YtcRQwxcwWJC0kH2a2BBhHesffHcdJgCQz6u5ASK3Rxszaxthj57Ix8+aaGPeuIm02A2qZ2cOEGGQ9sjbfRYiY/pBC+o1thU3ybKUFSc0zQ7AKSSsPpwIhqBzH2fpJsgf1TUJm0OwI4v8AjpNUn5CBdFouJ4kt0BIYH2Pj/Q3YJI+Nmf2ekJTtrzE6c7lI6iVpLiFVwW2SZmZtm0AYOhsgaa6kI/K1kwQKmVYPZ2PiuzTRAhinkAJlEvCcmaXCXV/SaEJaivbxc/1e0pqySbM+11Y8adaXlDaPxec4juOkkqTfQTmO4zhOTtxAOY7jOKmk0gZKUl1J10l6J4YEmijpqFKIqwriRN3B5dcsqK0h8bzfkTSkFG2WkpSHTmktaZykWTHU0QVJa8qQ8jBMqdUG6dbn2oonMX1mVqk/4DpCFs36cX1X4KTKtltBDXUqUHckMLgEx9yJMO9qJ2DHuLxjdZ53ARr7ETwZZyStJYe2FkCPuNyYEO6qY9K6oh4BjeJyXeBVoE/SutKuLe36XFvN01epHpSk7YCzgeEWvfHMbIGZPRS3D4w9qimS/i6pUSyfI+mqWD5dUodYfrA2Bjd9Q1JjBW7QxvBFJ8e6/SVNkDQWmKUQ5ugGSZMkTZP0g1hPkm6W9Jak54FdKnPOWRxB8Dz73ML8q+dI2TweS3folPlmNiUuLwNmE7wwE8cCaQ3DlFptkG59rq14ktJX2SG+rwMfmdkXZTcozEm6HDjMzHoArwM/yaryWSy/Fbgoll0EnGshusBBwErgREIIo26EsEg3SGoR6/cALjCzvYHvAUvNrBfQCzhbUjuCO3t7oCMh3XsmXXtlaQl8nLU+l5TcYGsaktoC+xKeylJBfOCZCnxKeBBxbQWSZn2urXiS0FeVThJ9CEbhpXhSQ4Dds7Zn5uZMBtrG5ZeA30s6H2hqZmuBbwCjLYTFWQD8m2CAAF4zsw/i8kDgjHisVwlRKfYiDHNl9v8EeKHUJ+oUT+xVPwxcmOtBJyksxWGY0qwN0q3PtRVPEvoqa6DeBdooRIUoiwhWtnv862hm2ZO7MhN01xGD1prZdcBZQEOCYetQzvG/LHO84VnHa2dmzxZzUgUyD2idtd4qljkFohAp5GFglJmlcTIxluIwTGnWBunW59qKpzr1VcpAmdkKQgihPykEJc2EsPk28ArQV9LXY/n2kvbeUnuS9jSz6WZ2PSG6QAdgAnBy7F42J/SIXsux+zPAOfGmh6S9FSIpvJi1fwvgkMqcc5njDZS0o6QdCT24Z0rU9laPJBG+O7MtRPhIDUpxGKY0a4N063NtxZOUvlLEpLsc+BXBUWEVoVfzSzNbKGkoMFohdFGm7ttbaOtCSYcQci/NBJ4CviLkiXqT8FLuEjP7X47e1Z2EocIp8ea3EDiBkKTvUGAWIefUxEqdbcTMPpd0DcGQAlxtZqlySFAIT9IfaKYQrukKM7srWVUb6EtI9jg9DssCXGpmTyYnaQMtgHsVIuvXAh6ylIRhIt3aIN36XFvxJKLPQx05juM4qcQjSTiO4zipxA2U4ziOk0rcQDmO4zipJBEDJelKSfOyokZMzXiIlPAYl5ayvTzHSG0sPqU8theApKaSxkj6r6TZkg5IWhOApPZlvptfSLowaV2Qbm2Qbn2urXiS0peIk4SkK4HlZnZjFR5juZk1qsL2dyJEx+hJ8C6cDOwXwx4lTvRk3N7MlkfX+/8Qom68krC0DUi6F5hgZnfGaQrbxTkWqSF6Lc0D9jezD5PWk02atUG69bm24qlOfaka4pP0iqROWevjJfVUmEN1d+wRvCHp+Lh9qKRHJD0dezG/jeXXAQ2jpR8V938i9iZmKMbzqySpjsWX9thekpoQ5rTdBWBmX6XNOEUGAO+l8UZBurVBuvW5tuKpNn1JGqgfZ3UXx8WyB4GTABQm1bYws9eBywjp4XsTJtreoDAJF0KcvpOBLoQJua3NbASwMkaUOJ1gOD4xs25m1hl4ugT6Ux+LT+mO7dWOMFftnvjQcWfWZ5omTgFGJy0iD2nWBunW59qKp9r0JWmg/pAVligT3eEhIJOr6SRgTFweCIyIN9vxQAOgTdz2LzNbamarCJNxs+P9ZZgOHC7pekkHmdnS0p9O+kh5bK86hGC/t5rZvoQJ3iOSlbQpcdhxEPD3pLWUJc3aIN36XFvxVLe+VA3xmdk8YJGkroRe0YNxk4BvZRm0NmY2O25bndXEhrh+Zdp9m3AznA78StIvSyC3xsTiS2lsr7nA3Kxe3RjCZ5QmjgKmxCDFaSPN2iDd+lxb8VSrvlQZqMiDwCVAEzObFsueAYbHF/9I2reAdtZoY1y+3YAVZvY34AZKcyNMdSw+pTy2l5n9D/hYUvtYNIDQA04Tp5LeoZY0a4N063NtxVOt+koRi69YfizpO1nrJ5jZHMKT9J+Aa7K2XQP8EZgmqRbwAXBsOe3fHutPAe4jvLdaD6wBzqms+BoQiy/tsb0AhgOj4rDB+8CZCevZQHwfdjjwg6S1lCXN2iDd+lxb8SShz2PxOY7jOKkkjUN8juM4juMGynEcx0knbqAcx3GcVOIGynEcx0klbqAcx3GcVOIGynEcx0klbqAcx3GcVOIGynEcx0klbqAcx3GcVOIGynEcx0klbqAcx3GcVOIGynEcx0klbqAcx3GcVOIGynEcx0kl5RooSSbpd1nrF0m6spx9TpDUsQT6crV9Z3ltF3p8ST+UdEaJdI2UNLj8mo7jOE4hFNKDWg2cKKlZBdo9AagSA2VmZ5lZeZlXCzq+mf3FzO4riTDHcRynpBRioNYSstP+uOwGSW0lvSBpmqR/SWoj6UBgECGD7VRJe5bZZ6SkWyW9Iul9Sf0l3S1ptqSRWfVulfS6pJmSrsoqHy+pZ1xeLulaSW/G9nbNdXxJZ0uaFOs9LGm7uP+Vki7Kavd6Sa9JelvSQbG8tqQb4v7TJP0glkvSzZLekvQ8sEsFrrvjOI5TDoW+g7oFOF1SkzLlfwbuNbOuwCjgJjN7GRgLXGxm3c3svRzt7QgcQDB6Y4E/AJ2ALpK6xzqXmVlPoCtwsKSuOdrZHnjFzLoBLwJn5zn+I2bWK9abDXwvz3nWMbPewIXAFbHse8BSM+sF9ALOltQO+CbQntBTOwM4ME+bjuM4ThEUZKDM7AvgPuD8MpsOAO6Py38FvlHgcf9pIdf8dGCBmU03s/XATKBtrHOSpCnAGwTjlWvI7ivg8bg8OWvfsnSWNEHSdOD02F4uHsnR1kDgDElTgVeBnYG9gH7AaDNbZ2afAC9s6YQdx3GcilGnAnX/CEwB7inBcVfH/+uzljPrdWIP5SKgl5ktjkN/DXK0syYaOoB15D+fkcAJZvampKFA/3J0ZbclYLiZPZNdUdLRedpwHMdxSkDBbuZm9jnwEJsOj70MnBKXTwcmxOVlQONK6NoB+BJYKmlX4KgK7l/2+I2B+ZLqRp0V4RngnLgvkvaWtD1hSPHk+I6qBXBIBdt1HMdxtkBF50H9Dsj25hsOnClpGvBd4IJY/gBwsaQ3yjpJFIKZvUkY2vsvYQjxpQo2Ufb4vyAMz70U26wIdwKzgCmSZgC3EXpXjwLvxG33ARMr2K7jOI6zBbRxhMxxHMdx0oNHknAcx3FSiRsox3EcJ5W4gXIcx3FSSeIGSlIdSQslXVem/NIKtLHF+HzZ0SeK0NdP0hRJa7Nj7UnqLmlijHQxTdLJxbRflUhqKmmMpP/GSB0HJK0pg6TWksZJmhWv4QXl71V9SDoyRgl5V9KIpPVkk2ZtkG59rq14EtFnZon+EVzIXwLeIzptxPLlBe5fu4A644GeReprS4hmcR8wOKt8b2CvuLwbMB9omvT1LKP9XuCsuFwvTfqAFkCPuNwYeBvomLSuzHcqfh/3iNftTddW8/W5tpqnL/EeFHAq8CfgI0JkCmJvqmGMpTeq7A4xBt/vJL0JHJDpIcU5SSMlzZA0XdKPy+xXK27/VaHizGyOmU0jTCLOLn/bzN6Jy58AnwLNK3bqVUcMS9UPuAvAzL4ysyWJisrCzOab2ZS4vIwQgqplsqo20Bt418zeN7OvCNMWjk9YU4Y0a4N063NtxZOIvkQNlKQGwGHAP4HRBGOFmY0AVlqIpZdrYu32wKtm1s3M/pNV3h1oaWadzawLm0a9qEOIF/iOmV1e4vPoTXiqyBV3MCnaAQuBe+J8sDvjBOPUIaktsC9hrloaaAl8nLU+l/QYzzRrg3Trc23Fk4i+pHtQxwLjzGwl8DBwgqTaBey3LtYvy/vAHpL+LOlI4IusbbcBM8zs2sqKziZGkfgrcKaFeIJpoQ7QA7jVzPYlROZI47h2I8JneaGFmI+O4zhA8gbqVOAwSXMIAVp3Bg4tYL9VZraubKGZLQa6Ed45/ZAQBSLDy8AhsddWEiTtADxBiLz+SqnaLRFzgblmlumVjCEYrNQQw0c9DIwys0fKq1+NzANaZ623imVpIM3aIN36XFvxJKIvMQMVb+4HAW3MrK2ZtQXOJQ7zAWsy8e8q0GYzoJaZPQxczqY35LuAJ4GHJFUkSG6+Y9UjhDu6z8zGVLa9UmNm/wM+ltQ+Fg0ghGVKBZJE+Exmm9nvk9ZThknAXpLaxc/5FEIKlzSQZm2Qbn2urXgS0VfpG3Ul+CbwgpllRzP/B/BbSfUJSRKnSZqS5z1ULloS3rlkDO/Pszea2e+j88BfJZ1eyJCcpF4EQ7QjcJykq8ysE3ASwQlhZ4UI6QBDzWxqgVqrg+HAqPiFeh84M2E92fQlxG+crpDKBOBSM3syOUkBM1sr6TxCoODawN1mNjNhWUC6tUG69bm24klKn8ficxzHcVJJ0u+gHMdxHCcnbqAcx3GcVFJpAyWprqTrJL0TQwJNlFTRBIPVRpyoO7j8mgW1NSSe9zuShpSizVIi6W5JnyrksUoVrq140qzPtRVPmvUlpa0UPahrCGFrOptZD+AEKpdNt8KUwiuviGPuBFwB7E+YZX2FpB2rW0c5jASOTFpEHkbi2oplJOnVNxLXViwjSa++kSSgrVIGStJ2wNnA8Iw3npktMLOH4vaBsUc1RdLf46RMJM2RdFUsny6pQyw/OIY3mhqjHzRW4Ias8EUnx7r9JU2QNBaYpRDm6AZJkxSCt/4g1pOkmxWCHD4P7FKZc87iCOA5M/s8zr96jpR9uczsReDzpHXkwrUVT5r1ubbiSbO+pLRVtgf1deCjXBEA4pyky4HDYs/qdeAnWVU+i+W3AhfFsouAc82sO2GO1ErgREIIo26EsEg3KERvgDDP6QIz2xv4HrDUzHoBvYCzJbUjuLO3BzoCZwAHVvKcM6Q9NInjOE6NpiqHxvoQjMJLYU4m9YCJWdszkQMmE4wQhKjmv1cIEPuImc2V9A1gdIwcsUDSvwkG6AvgNTP7IO47EOia9X6pCbAXYa5SZv9PJL1QBefqOI7jlJjKGqh3gTaSdsjRixJhCOzUHPsBZCborsvoMLPrJD0BHE0wbEeUc/wvyxxvuJk9s4kI6egCzqMY5gH9s9ZbEUIsOY7jOCWgUkN8ZraCEK7mTzFaAZKaS/o28ArQV9LXY/n2kvbeUnuS9jSz6WZ2PSG0RgdgAnByfMfUnNAjei3H7s8A5yiGR5K0t0L07hez9m8BHFKZcy5zvIGSdozOEQNjmeM4jlMCSuHFdzkhrcOs6IL4OPCFmS0EhgKjJU0jDO91KKetC6MzxDRgDfAUIczQNEKCrBeAS2KcubLcSYg1NyXquI3QM3sUeCduu49NhxmLxsw+J3gwTop/V8ey1CBpNOF820uaK+l7SWvK4NqKJ836XFvxpFlfUto81JHjOI6TSjyShOM4jpNK3EA5juM4qcQNlOM4jpNKEjFQkq6UNC8rasRUSU1LfIxLS9lenmOkPRbfnBh9Y6qk15PWk42kC6JDzExJFyakYbP4YpK+HTWtl9QzCV1p15Z2fa5t69GXZA/qD2bWPetvSYnbr1IDpZoRiw/gkHh9E/3iZyOpMyFEVm9ChJBjM9MRqpmRbB6eagZh4viL1a5mU0aSXm2Qbn0jcW3FMpIU6UvVEJ+kVyR1ylofL6lnnEN1t6TXFGL0HR+3D5X0iKSnYy/mt7H8OqBh7DmMivs/IenN+NR+cgnkpj4WX4rZB3jVzFaY2Vrg32yMJlJt5IovZmazzeyt6tZSljRrg3Trc23FkzZ9SRqoH2cN742LZQ8SUqkTJ9W2MLPXgcsI6eF7Eyba3hAn4UKI03cy0IUwIbe1mY0AVsaew+kEw/GJmXUzs87A0yXQXxNi8RnwrKTJkr6ftJgsZgAHSdpZIeDw0UDrhDU5jpMy0jLEl4nu8BCQiaV3EjAmLg8ERkiaSggn1ABoE7f9y8yWmtkqwmTc3XMcazpwuKTrJR1kZktLfzqp5BsxIO9RwLmS+iUtCMITGXA98CzhYWEqIeSV4zjOBlI1xGdm84BFkroSekUPxk0CvpVl0NrEmxxsjOkHWXH9yrT7NiHy+XTgV5J+WQK589j0qb9VLEsN8XpiZp8SImr0TlbRRszsLjPbz8z6AYuBt5PW5DhOukiVgYo8CFwCNDGzabHsGWC4FMKiS9q3gHbWaGNcvt2AFWb2N+AGgrGqLKmOxRffuzXOLBP0pSZTp6Rd4v82hPdP9yeryHGc1GFm1f4HXEnobUzN+msbt+0KrAWuyKrfkBBbbzowE3g8lg8Fbs6q9zjQPy5fD8wGRhEcGqbF40wCepboPIYRIrq/C5yZxLXcgrY9CPEL34zX7LKkNZXRN4EwJPsmMCAhDaOB+YS4j3MJOcW+GZdXAwuAZ1xbzdLn2rYefR6Lz3Ecx0klaRzicxzHcRw3UI7jOE46cQPlOI7jpBI3UI7jOE4qcQPlOI7jpBI3UI7jOE4qcQPlOI7jpBI3UI7jOE4qcQPlOI7jpBI3UI7jOE4qcQPlOI7jpBI3UI7jOE4qcQPlOI7jpBI3UI7jOE4qcQPlOI7jpBI3UI7jOE4qcQPlOI7jpBI3UI7jOE4qcQPlOI7jpJL/B6EFjgWNgYzQAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD+CAYAAABvEpGeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7LklEQVR4nO3dedyVc/7H8de7naSo0IIYlPZSRCT7MoTBYBhiLGPIMCOTZcY25sdkxhhm7EkmYazN2LdGiAppRZZQkkQJlZbP74/v9+TqdO7uc3ef+5zr1Of5eNyP+zrX8r0+57rPfT7X+vnKzHDOOefSplapA3DOOedy8QTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBudST1F/SS6WOA0DSFEl9C9DODEn7Vj+i9Zskk7RdBdNGSTq12DGVWpr+X6rLE1TKxC+uRZK+kTRH0lBJG61lW9X6B5XURtILkr6T9HZav1Al/SLGtzBus8clNaqJdZlZBzMbVRNtr2/K6YtU0q2S3pG0QlL/HNO3lfTf+Bn8QtKfSxDmOscTVDodamYbAd2BHsAlVVlYQSH+tiOAN4GmwMXAA5KaF6DdgpG0J/An4DgzawTsCNy3lm3VKWRs6wJJtUsdQ0q8BfwKeCN7gqR6wDPA88AWQGvgX0WNbh3lCSrFzGwW8ATQUdImcQ9trqSv4nDrzLzxaOkqSS8D3wF3A3sAN8ajsRsl/UPSX5LrkDRS0nnZ65a0AyFBXmpmi8zsQWAScGQ+sUv6saQ3JX0t6RNJlyWmtYmnZk6O076S9EtJPSVNlDRf0o15bqaewBgzezNusy/N7C4zW5jYLiuPIrP32mMcZ0maDkyXdJOka7Pey6OSfhOHZ0jaV1LLeKS7aWK+bnHvua6kH0l6XtK8OG64pCZ5vqc1qqxtSd3jtl8o6d+S7pP0x8T0CyTNlvSppFOTp8niEftN8Sj0W2Cv+F4fjJ+9DyWdk2hrA0l3xb/htNj2zMT0QZLej7FMlXREHL8jcDOwa/x8zo/j60u6VtLHCkfDN0vaINHewETsp+SxuX4kaWz8HD6a+XtJekzSgKztOjETXzYz+4eZPQcszjG5P/Cpmf3VzL41s8VmNjFXO5Iul3RDHK4r6VtJgxPbcnEixl6SXon/D28pcWpZUmNJd8RtMUvSH1XBzoSkwZJeistsJ+l/khbEz85a7cwVjZn5T4p+gBnAvnF4S2AKcCXhKOZIYEOgEfBv4JHEcqOAj4EOQB2gbhx3amKenYFPgVrxdTNCMts8RxxHANOyxt0I3BCHdwfmr+F99AU6EXaCOgNzgMPjtDaAEb6gGgD7E/7xHwE2A1oBnwN7xvn7Ay9VsJ49gEXA5UBvoH7W9OxtsEpbMY5ngE2BDYA+wCeA4vRNYvstc/x9ngdOS7Q1GLg5Dm8H7AfUB5oDLwJ/y/V3XovPSIVtA/WAj4Bfx8/AT4DvgT/G6QcCn8XPyYaEPX0DtovThwIL4rasFed5HfhDbHtb4APggDj/1cD/4nZqDUwEZiZiPRpoGds6BvgWaFHR3xW4DhgZ/x6NgP8A/5eIfQ7QEWgI3JOMPcd2GgXMSsz/IPCvOO2nwGuJebsA84B6lWz7l4D+WeOGEHYInwC+iOvtVMHyewOT4vBuwPuZOOK0t+JwqxjPwXHb7RdfN4/THwZuie9rM2AscEZyu8blbgOeAjaM00YQzobUIvzv7V7q77w1bu9SB+A/WX+Q8MX1DTCf8EXzT2CDHPN1Bb5KvB4FXJE1zygSX85x3DRgvzh8NvB4BXH8HHg1a9xVwNC1fF9/A66Lw23iF0urxPR5wDGJ1w8C58bh/lSQoOL0gwhfZPPjtvsrUDvXNshuK8axd+K1CIm+T3x9GvB81t8nk6BOzUyLy32SWS5HjIcDb+ZqpwCfmZVtExLsLGKCjeNe4ocENYT4hR9fb8fqCWpYYvouwMdZ67sQuDMOr0xWiW0ycw2xTgAOq+BvIUIC+1Fi3K7Ah4nYr05M24HKE1Ry/vaEZF2b8OX8FbB9nHYt8M88tnWuBPU0sDR+DusBA+N2WS3ZEXaCFhN2OAcBFwEzgY0IO1l/j/P9Drg7a9mngJOAzYElJL4XgOOAFxLb9TXCqe4Hk3EAw4BbgdaF+OzV9I+f4kunw82siZltbWa/MrNFkjaUdIukjyR9TdhrbpJ1WP9JHm3fBZwQh08g7Pnl8g2wcda4jYGF+bwBSbso3GAxV9IC4JeEI7akOYnhRTle53VziJk9YWaHEva6DyP8g1bl5pCV283Cf/G9hH94gJ8BwytY7kHCKaoWhMSwAhgNIGlzSffG0y9fE45Ust//aiRtFU95fSPpmwrmWVPbLYFZ8X2s9v7i9E8qmJZr3NZAy3iaaX48FXcR4Uuy0vYknShpQmLZjlS8HZoTj9gS8z8Zx+da10cVtFPRe/mIcFTZzMwWE77AT1C4XnscFf8vVGYRIdE+YWbfE5JdU8L10FWY2SJgPLAn4TPzP+AVwhHrnvE1hO1+dNZ23x1oEafVBWYnpt1COJLK2I7wv3B5jCnjAsKOwFiFO1LzOU1aMp6gysdvgbbALma2MeHDDeHDlpFdmj5Xqfp/AYdJ6kL4B3qkgvVNAbbVqnfDdYnj83EP4VTNlmbWmHA6T2tepHrMbIWF6wTPE74IIeyRb5iYbYtci2a9HgEcJWlrwhHEgxWs7yvC3vMxhER2byIx/Cm22yn+vU4gj/dvZh+b2UaZnwpmW1Pbs4FWkpLr2jIxPJtwKi7XtJVhJIY/IRzBNEn8NDKzgytrL26/2whH6k3NrAkwORFr9nb/gvBl3yGxrsaJ7TA7K96tcsSeLXv+pXE9EHbWjgf2Ab4zszF5tJfLRHL/r1Xkf4TTed2AcfH1AYRT8C/GeT4hHEElt3tDM7s6TltCSLSZaRubWYfEOqYBJwNPSGqbGWlmn5nZaWbWEjgD+KcquE0/DTxBlY9GhH/e+fEi6qV5LDOHcM1gJTObSfinuBt4MO7RrcbM3iWcjrlUUoN48bgzFXxZVxDvl2a2WNLOhC/wgpN0mKRjFW4iUVzXnsCrcZYJwE/iEeh2wC8qa9PCDRdfALcDT5nZ/DXMfg9wInBUHM5oRDgKXSCpFeG0T6Gsqe0xwHLgbEl1JB1G+OLLuB84WdKOkjYEfl/JusYCCyX9Ll7Ery2po6SeifYujNu/FSEZZTQkfHHPBZB0Mj/sOED4fLZWuAsOM1tBSGjXSdosLtNK0gGJdfWX1D7Gns//wAmJ+a8AHjCz5XF9YwhHvX+hkqMnSfUkNSAk17rxfyLz/fkvoJfCzTO1gXMJn59pFTT3P8JnZmo8uhlFOOL/0MzmJto8VNIBcZs3kNRXUmszm03YMfqLpI0l1VK4cWbP5ErMbAThaPdZST+K7+No/XBz1VeEv8+KNW7BEvIEVT7+Rjh//QXhy/fJPJa5nnAk8JWkvyfG30W4gaGyUxrHEm5z/4pwMfyozD+QpD0qOgUV/Qq4QtJCwgX2+/OId218RbhONB3InO4abGaZ03LXEa47zCG874pO12W7B9iXVZNOLiOB7YHPzOytxPjLCXdBLgAeAx7Kc735qLDt+IX3E0Iink84uvovYY8bM3sC+DvwAvAePyTyJblWFL/MDyFc8/yQHxJ34zjLFYRrKB8CzwIPJNY1lfDlP4aw/TsBLyeaf55wRP6ZpMxRze8yccXTl88SzhxkYv9bXO69+LsydxOuq31GuO50Ttb0YTGuym4Lf5qwg7gb4RrOIuJZDDN7h7CdbyZ8Hg8D+mWdWkt6hfC/nDlamkq4LpV5jZl9Etu5iJDgPyHsiGS+s08kXO+aGtf5AOH03yrM7C7C3+h5SW0Id72+Fv93RwK/NrMPKnnvJZO5U8mtRyT1IfxDbm3+AVjnSXqNcHfhnTmm7Ug47VbfzJYVYF1nAsea2Z6VzpwCkk4ETjez3Usdi1udH0GtZyTVJdyCfLsnp3WTpD0lbRFP8Z1EODX7ZGL6EQrPG20CXAP8Z22Tk6QWknrH00xtCddKHy7E+6hp8bTfrwhHRC6FPEGtR+Le8nzCqYC/lTQYV5PaEiofzCckjKPidYuMMwjPmb1PuF51ZjXWVY9wB9lCwim3RwmPRqRavK41l3DqsbLTuK5E/BSfc865VPIjKOecc6nkCco551wqlXX1ZklDCLfAfm5mHXNMF+FW64MJNef6m9lq1YizNWvWzNq0aVPgaJ1zzmV7/fXXvzCznL0klHWCIjzfcCPhWYZcDiI8o7I9oSLATfH3GrVp04bx48cXKETnnHMVkVRhyaqyTlBm9mJ8+KwihxEKXxrhwb8mklpk3dFUMK/+8zQaza/o4fGa8832R7DL0b8t+nqdc64mlXWCykMrVi0WOTOOWy1BSTodOB1gq63yKfGV24pl1X7WsUq2Xv4RH037N+FuYuecW3es6wkqb2Z2K/GBvR49eqzVvfe9fnVbQWPKx5Q/7Q5FTorOOVcM63qCmsWq1Yxbx3HOuTKzdOlSZs6cyeLFuTq1dWnXoEEDWrduTd26dfNeZl1PUCMJVZ3vJdwcsaCmrj8552rWzJkzadSoEW3atGHV3kRc2pkZ8+bNY+bMmWyzzTZ5L1fWCUrSCELX4s0kzSSU368LYGY3A48TbjF/j3Cb+cmlibRmGTByQnEPDBttUJe92m5W+YzOFcjixYs9OZUpSTRt2pS5c+dWPnNCWScoMzuukukGnFWkcErHjOaNGhR1lXMX+mkWV3yenMrX2vztvJKEc87lSRInnHDCytfLli2jefPmHHLIIWtcbvz48ZxzTnZXVKuaP38+//xnfnV2d9ttt7zmq8yMGTPo2HG1Ggep4QnKOefy1LBhQyZPnsyiRaEj6meeeYZWrVpVulyPHj34+9//vsZ5qpKgXnnllbzmK3eeoJxzrgoOPvhgHnvsMQBGjBjBccf9cKVh7Nix7LrrrnTr1o3ddtuNd955B4BRo0atPMq67LLLOOWUU+jbty/bbrvtysQ1aNAg3n//fbp27crAgQP55ptv2GeffejevTudOnXi0UcfXbmejTbaaGW7ffv25aijjqJdu3Ycf/zxZHqoeP3119lzzz3ZaaedOOCAA5g9e/bK8V26dKFLly784x//qOGtVT1lfQ3KBUvNz8u79cvl/5nC1E+/Lmib7VtuzKWHdqh0vmOPPZYrrriCQw45hIkTJ3LKKacwevRoANq1a8fo0aOpU6cOzz77LBdddBEPPvjgam28/fbbvPDCCyxcuJC2bdty5plncvXVVzN58mQmTJgAhNOHDz/8MBtvvDFffPEFvXr1ol+/fqtdy3nzzTeZMmUKLVu2pHfv3rz88svssssuDBgwgEcffZTmzZtz3333cfHFFzNkyBBOPvlkbrzxRvr06cPAgQOrv+FqkCeodcDSFaWOwLn1R+fOnZkxYwYjRozg4IMPXmXaggULOOmkk5g+fTqSWLp0ac42fvzjH1O/fn3q16/PZpttxpw5c1abx8y46KKLePHFF6lVqxazZs1izpw5bLHFFqvMt/POO9O6dWsAunbtyowZM2jSpAmTJ09mv/32A2D58uW0aNGC+fPnM3/+fPr06QPAz3/+c5544olqb5Oa4gnKOVd28jnSqUn9+vXj/PPPZ9SoUcybN2/l+N///vfstddePPzww8yYMYO+ffvmXL5+/forh2vXrs2yHNVghg8fzty5c3n99depW7cubdq0yfmQcq62zIwOHTowZsyYVeadP39+Fd9pafk1KOecq6JTTjmFSy+9lE6dOq0yfsGCBStvmhg6dGiV2mzUqBELFy5cpa3NNtuMunXr8sILL/DRRxUW/V5N27ZtmTt37soEtXTpUqZMmUKTJk1o0qQJL730EhCSYJp5gnLOuSpq3bp1ztvGL7jgAi688EK6deuW86hoTZo2bUrv3r3p2LEjAwcO5Pjjj2f8+PF06tSJYcOG0a5du7zbqlevHg888AC/+93v6NKlC127dl1559+dd97JWWedRdeuXVfeUJFWSnuApdCjRw8rl/6gpvxpdxYuXsblTf9c1PV2btWYa47qUtR1uvXbtGnT2HHHHUsdhquGXH9DSa+bWY9c8/s1qDJy3TPvcv1z01cZd2+9sJc2bfbC1eZvtlE9mjeqv9r46vpo3ncsXeZ3ZjjnapYnqDJy3n47cN5+O6w68s5/8uqH83h636rVuKqO88c3Zun3S4q2Pufc+skT1Dpi6QbNi7Yuqw1a6gnKOVez/CYJ55xzqeRHUG6tGOZdfDjnapQnKLdW6taq5V18OOdqlJ/iWwc0r7+81CE4t16QxG9/+9uVr6+99louu+yyNS7zyCOPMHXq1BqJ59RTT6207XzXf/PNNzNs2LCCxNW/f38eeOCBarfjCWodsFn9qj0Q6JxbO/Xr1+ehhx7iiy++yHuZmkxQt99+O+3bty/I+n/5y19y4oknFiq0gvAE5ZxzeapTpw6nn34611133WrTZsyYwd57703nzp3ZZ599+Pjjj3nllVcYOXIkAwcOpGvXrrz//vurLNO/f3/OPPNMevXqxbbbbsuoUaM45ZRT2HHHHenfv//K+c4880x69OhBhw4duPTSS1eO79u3L5miAhtttBEXX3wxXbp0oVevXsyZMyfn+m+77TZ69uxJly5dOPLII/nuu++A0A3Itddeu7Ld3/3ud+y8887ssMMOK6u1L1++nIEDB9KzZ086d+7MLbfcAoTCtmeffTZt27Zl33335fPPPy/M9i5IK845V0xPDILPJhW2zS06wUFXVzrbWWedRefOnbngggtWGT9gwABOOukkTjrpJIYMGcI555zDI488Qr9+/TjkkEM46qijcrb31VdfMWbMGEaOHEm/fv14+eWXuf322+nZsycTJkyga9euXHXVVWy66aYsX76cffbZh4kTJ9K5c+dV2vn222/p1asXV111FRdccAG33XYbl1xyyWrrb9KkCaeddhoAl1xyCXfccQcDBgxYLa5ly5YxduxYHn/8cS6//HKeffZZ7rjjDho3bsy4ceNYsmQJvXv3Zv/99+fNN9/knXfeYerUqcyZM4f27dtzyimn5LXZ18SPoJxzrgo23nhjTjzxxNV6yB0zZgw/+9nPgNCNRaYga2UOPfRQJNGpUyc233xzOnXqRK1atejQoQMzZswA4P7776d79+5069aNKVOm5DxlV69evZWdIu60004rl802efJk9thjDzp16sTw4cOZMmVKzvl+8pOfrNbW008/zbBhw+jatSu77LIL8+bNY/r06bz44oscd9xx1K5dm5YtW7L33nvn9d4r40dQzrnyk8eRTk0699xz6d69OyeffHK128p0l1GrVq1Vus6oVasWy5Yt48MPP+Taa69l3LhxbLLJJvTv3z9ntxt169Zd2ZlhRV14QDit+Mgjj9ClSxeGDh3KqFGj1hhXsi0z44YbbuCAAw5YZd7HH3+8am86T34EtQ6o9/1XpQ7BufXKpptuyk9/+lPuuOOOleN222037r33XiB0Y7HHHnsAq3ejUVVff/01DRs2pHHjxsyZM6fKHQxmr3/hwoW0aNGCpUuXVrm7jQMOOICbbrppZUeM7777Lt9++y19+vThvvvuY/ny5cyePZsXXnihSu1WxBPUOqC+Jyjniu63v/3tKnfz3XDDDdx555107tyZu+++m+uvvx4IXcQPHjyYbt26rXaTRD66dOlCt27daNeuHT/72c/o3bt3lZbPXv+VV17JLrvsQu/evavUhQeE29rbt29P9+7d6dixI2eccQbLli3jiCOOYPvtt6d9+/aceOKJ7LrrrlVqtyLe3UYO5dTdBnf+GD56iSn73VO0VQ56BbT8ey48uk/R1gnhQd1+XVsVdZ0uPby7jfJX1e42/AjKOedcKnmCcs45l0p+F986os34K4u2rj8tghdr7QIU9xSfc2794gmqnLzwf/C/3LfXNvxq2mrjvm/QrEb6idpmxUd8v8L4vuAtO7dmZrbyVmpXXtbmfgdPUOVkrwvDT7bLGhf1JokVz13JUi//54qsQYMGzJs3j6ZNm3qSKjNmxrx582jQoGo9IHiCcmut6Yf/Ker6ln5fl5H0Leo6XXrUsto0/XouH836DPC7j9OilkSDurUrna9Bgwa0bt26Sm2XdYKSdCBwPVAbuN3Mrs6a3h8YDGR61rvRzG4vapDrsGJ2Mw+wBXOpW+Q+qFzaNPRTyykzd+Fi+u1YM49/lG2CklQb+AewHzATGCdppJllF6m6z8zOLnqAzjnnqqVsExSwM/CemX0AIOle4DCgZjpecasZ9Epx17d38wb02qa463TOlU45J6hWwCeJ1zOBXXLMd6SkPsC7wHlm9kmOeZB0OnA6wFZbbVXgUGvWpy32q7G2h78D97y76rh764Xfk+atPv9mG8DmGxY+jg8WgJbXp1fhm3bOpVQ5J6h8/AcYYWZLJJ0B3AXkrANvZrcCt0IodVS8EKvv05YHsEENtX182/CT1GZ8SE6PHVpDK81h0CuA92zv3HqlnCtJzAK2TLxuzQ83QwBgZvPMbEl8eTuwU5Fic845V03lnKDGAdtL2kZSPeBYYGRyBkktEi/7Aas/zeqccy6VyvYUn5ktk3Q28BThNvMhZjZF0hXAeDMbCZwjqR+wDPgS6F+ygJ1zzlVJ2SYoADN7HHg8a9wfEsMXAjlKL6xbNqhXm/nfFe/pkGXLV1DeB9/OuXJQ1gnKBR1abgwbNSveCt+rS7PFy1hS+ZzOObfWfDfYrZXmDcrqRkfnXBnyBOWccy6VPEE555xLJU9QzjnnUskTlHPOuVTyBOWccy6VPEE555xLJU9Qrmx8trjyXjudc8X1+KTZNda2JyhXNj73BOVc6jw5ZU6Nte2VJNxaafjtx7QZf2XR1venRXBP7d3Ae4Rybr3hCcpV3bZ9+XbR0qJ+eLZZ8RGH1QavX+Hc+sMT1LqgQWP4puYOs1fTshsfLNqCes23LdoqVzx3Jawo2uqccyngCWpdsH3NdflekeUzbiv6OgGu+O+UkqzXOVd8nqBc6gx/B+55d9Vx99YLv6fNXrja/M02qkfzRvWLEJlz66+5C5fwxTe5u/VpM+ix1cb9ep/tOW+/Haq1Tk9QLnWObxt+ktqMh0nzYMRpfpOEc2ly3G2vMuPqH9dI236buXPOuVTyBOWccy6VPEE555xLJb8G5dbKBvVqM/+73BdMa8Ky5Svw/Snn1i+eoNxa6dByY9ioWfFW+F5dmi/6ns+Lt0bnXB4O7LB5jbXtu6SubGxWf1mpQ3DOZTm4U4saa9sTlHPOuVTyBOWccy6VPEE555xLJU9QzjnnUskTlHPOuVTyBOWccy6V/Dkot3aK3QfV8qWYvMt359YnnqDc2il2H1Tj78S+n1vcdTrnSqqsE5SkA4HrgdrA7WZ2ddb0+sAwYCdgHnCMmc0odpyuMGrVEnMXLi51GM65hEYb1K2xtss2QUmqDfwD2A+YCYyTNNLMpiZm+wXwlZltJ+lY4BrgmOJH6wqhUf069OvaqtRhOOeKpJxvktgZeM/MPjCz74F7gcOy5jkMuCsOPwDsI0lFjNE559xaKucE1Qr4JPF6ZhyXcx4zWwYsAJrmakzS6ZLGSxo/d65f60idLTrBptuWOgrnXBGV7Sm+QjOzW4FbAXr06GElDsdlO+jqyudxzq1TyvkIahawZeJ16zgu5zyS6gCNCTdLOOecS7lyTlDjgO0lbSOpHnAsMDJrnpHASXH4KOB5M/OjI+ecKwMq5+9rSQcDfyPcZj7EzK6SdAUw3sxGSmoA3A10A74EjjWzD/Jody7w0VqG1Qz4Yi2XLRWPuTg85uLwmIujUDFvbWbNc00o6wSVRpLGm1mPUsdRFR5zcXjMxeExF0cxYi7nU3zOOefWYZ6gnHPOpZInqMK7tdQBrAWPuTg85uLwmIujxmP2a1DOOedSyY+gnHPOpZInKOecc6nkCaqAJB0o6R1J70kaVOp4cpG0paQXJE2VNEXSr+P4yyTNkjQh/hxc6liTJM2QNCnGNj6O21TSM5Kmx9+blDrODEltE9tygqSvJZ2btu0saYikzyVNTozLuV0V/D1+vidK6p6imAdLejvG9bCkJnF8G0mLEtv75hTFXOFnQdKFcTu/I+mAFMV8XyLeGZImxPE1s53NzH8K8EN4WPh9YFugHvAW0L7UceWIswXQPQ43At4F2gOXAeeXOr41xD0DaJY17s/AoDg8CLim1HGu4bPxGbB12rYz0AfoDkyubLsCBwNPAAJ6Aa+lKOb9gTpx+JpEzG2S86VsO+f8LMT/x7eA+sA28Xuldhpizpr+F+APNbmd/QiqcPLp/qPkzGy2mb0RhxcC01i9Cny5SHanchdweOlCWaN9gPfNbG2rk9QYM3uRUGUlqaLtehgwzIJXgSaSWhQl0IRcMZvZ0xZ6LAB4lVCbMzUq2M4VOQy418yWmNmHwHuE75eiWlPMsduinwIjajIGT1CFk0/3H6kiqQ2hDNRrcdTZ8RTJkDSdLosMeFrS65JOj+M2N7PZcfgzYPPShFapY1n1HznN2xkq3q7l8hk/hXCkl7GNpDcl/U/SHqUKqgK5PgvlsJ33AOaY2fTEuIJvZ09Q6ylJGwEPAuea2dfATcCPgK7AbMLhe5rsbmbdgYOAsyT1SU60cJ4hdc9MxELG/YB/x1Fp386rSOt2rYiki4FlwPA4ajawlZl1A34D3CNp41LFl6WsPgtZjmPVna4a2c6eoAonn+4/UkFSXUJyGm5mDwGY2RwzW25mK4DbKMEphTUxs1nx9+fAw4T45mROMcXfn5cuwgodBLxhZnMg/ds5qmi7pvozLqk/cAhwfEysxNNk8+Lw64TrOTuULMiENXwW0r6d6wA/Ae7LjKup7ewJqnDy6f6j5OK54zuAaWb218T45LWEI4DJ2cuWiqSGkhplhgkXxCezancqJwGPlibCNVplTzPN2zmhou06Ejgx3s3XC1iQOBVYUpIOBC4A+pnZd4nxzSXVjsPbAtsDlfZoUAxr+CyMBI6VVF/SNoSYxxY7vjXYF3jbzGZmRtTYdi72nSHr8g/hLqd3CXsPF5c6ngpi3J1wymYiMCH+HEzolmRSHD8SaFHqWBMxb0u4q+ktYEpm2wJNgeeA6cCzwKaljjUr7oaEDjIbJ8alajsTkudsYCnhWscvKtquhLv3/hE/35OAHimK+T3CdZvMZ/rmOO+R8TMzAXgDODRFMVf4WQAujtv5HeCgtMQcxw8Ffpk1b41sZy915JxzLpX8FJ9zzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOVcJSU0k/SrxuqWkBwrU9mWSzo/DV0jat0DttpD030K0VUH731Rh3mdTWnPQpZwnKOcq1wRYmaDM7FMzO6rQKzGzP5jZswVq7jeE8jlpcDeJ7edcvkqWoCSZpH8lXteRNLeyvT5JPST9vZJ5VtnjrWTeV/KLuNJ22iQ79nLrlKuBH8WO2AYn/9aS+kt6RKFjvxmSzpb0m1jV+VVJm8b5fiTpyViNfbSkdtkrkTRU0lFxeIakyyW9odBRY7s4vmGsfD02rqOiLl2OBJ6MyzwmqXMcflPSH+LwFZJOi8MDJY2LlbUvT8R0QlzXBEm3ZMrZJKY3kzRG0o/jUduLcd7JiYrWIwkln5yrklIeQX0LdJS0QXy9H3kURDSz8WZ2TiWzNSHPPTYz2y2f+dx6bRChP6euZjYwx/SOhOKZPYGrgO8sVHUeA5wY57kVGGBmOwHnA//MY71fWKjgflNcBkIJnOfNbGdgL2BwrE+4Uqzf9pWZLYmjRgN7SGpMqPTdO47fA3hR0v6E2mk7Eypr7ySpj6QdgWOA3mbWFVgOHJ9Yz+bAY4RO6x4DfgY8FeftQih7g5l9BdSX1DSP9+zcSqU+xfc48OM4nF1Uc+e4Z/ampFcktY3j+2aOsuL5+yGSRkn6QFImcWXv8W4k6bnE3uhhifV8k2h3lKQHFLqOHi5JcdpOCn2cvC7pKf1Q6XknSW9Jegs4q2Y3lUuxF8xsoZnNBRYA/4njJwFtFLo22Q34t0IX2bcQejauzEPx9+uEHkshFModFNsZBTQAtspargUwN/F6NKF31N6EhLKRpA2Bbczsndjm/sCbhDpq7QgJax9gJ2BcXN8+hLqIAHUJ9fouMLNn4rhxwMmSLgM6WegQM+NzoGUe79m5leqUeP33An+ICaczMISwVwfwNrCHmS2LF47/RDhtka0dYU+yEfCOpJsIe7wd455cpjz8EWb2taRmwKuSRtrqhQi7AR2AT4GXgd6SXgNuAA4zs7mSjiHsJZ8C3AmcbWYvShpciA3iytKSxPCKxOsVhP+xWsD8zOdxLdpdzg//qwKOjImlIosIiStjHNCDUF36GaAZcBoh8WXa/D8zuyXZiKQBwF1mdmGOdSyLyx8A/A9CD6wK/XT9GBgq6a9mNizO3yDG5VzeSnoEZWYTCXuGxxGOppIaE/Y4JwPXERJHLo9Z6IvkC8JeWq5eVQX8SdJEQnXmVhXMN9bMZlron2VCjK0t4RTOM3Ev8hKgtaQmQBML3SJDuBDs1k0LCTtAa8VCh5AfSjoaQpcnkrqsZXNPAQMSR/fdcszzLj8ccWFm3xMqfR9NOO04mnDKMPPZfQo4JR7pIamVpM0IR0hHxWEkbSpp60yzhJ20dpJ+F6dvTehl9TbgdqB75v0CWwAz1vI9u/VUqY+gIFxAvRboSyjzn3El4dTJEQpdk4+qYPnk3mtyTzPpeKA5sJOZLZU0g1X3MNfUloApZrZrcsaYoNx6wMzmSXo57iw9QehyoqqOB26SdAnh9Ni9hO5DqupK4G/AREm1gA8JnfQl4/1W0vuStjOz9+Lo0cA+ZrZI0mhCJ3ij4/xPx+tNY2Le+wY4wcymxnifjutaSjiV/VFcbrmk44CRkhYSrisPlLQ0tpG5/rYT8KqZLVuL9+vWY2lIUEMIpz8mSeqbGN+YH26a6F/FNrP3eBsDn8fktBewde7FcnoHaC5pVzMbo9Ab7Q5mNkXSfEm7m9lLJC4eu3WPmf0sa1THOH4ooX+czHxtEsMrp5nZh8CBOdq9LDHcv4J2xhN24DCzRcAZeYR8I+H/5pK43O+B38fhTwk7Xsk4rgeuzxHffSR6Tk2M3yj+XkI4zZdxV45Yfk5+N4U4t4pS3yRBPKWW67bxPwP/J+lNqphILXQ9/HK81XUwMBzoIWkSYa/u7Sq09T1wFHBNvBliAuGCN8DJwD/iqT/lbMC5EjCzh0nPKbXJZvZcqYNw5cc7LHTOOZdKJT+Ccs4553LxBOWccy6VPEE555xLpZInKP1Qg+/qrPEXVaGN2yW1X8P0UZJ6rGV8fWIFimWKddIS056Md/LVWNXomiDpPElT4k0kIyTluuU+NWK1kM9VRrUOPebi8JiLo1QxlzxBEWrwvQscnXn4MMorQUmqbWanmtnUGokOPibcrntPjmmDCbfQlg1JrYBzgB5m1hGoDRxb2qgqNZQct2in3FA85mIYisdcDEMpQcxpSFDHEZ6/+BjYFSAeTW2gUEtvePYCkr6R9Jd42/eumSMkSbUVKkJPVqi5d17WcrXi9D/mG5yZzYgVL1bkmPYc4ZmrclOHsH3rABsSSjulVqzW8WWp46gKj7k4PObiKFXMJX1QN55a2pfw4GETQrJ6xcwGSTp7DbXLGgKvmdlvYzuZ8V2BVvHIILvaQx3C81CTzeyqgr6RMmJmsyRdS9ghWAQ8bWZPlzgs55xbTamPoA4hlDNaBDwIHK6s/mYqsDzOn+0DYFtJN0g6EPg6Me0W1vPkBKDQs+lhwDaE6tINJZ1Q2qicc251pU5QxwH7xtp4rxNq8e2dx3KLzWx59sjY70wXQt2+XxIKVma8AuyV9hsCimBf4EMzm2tmSwldOnifWM651Cllj7obE7rW2MrM2sTaY2fxQ8+bS2Pdu6q02QyoZWYPEmqQdU9MvoNQMf3+eO1lffUx0EvShvGmlH2AaSWOyTnnVlPKI6gjCD2DJiuIPwocKqk+oQfSiblukliDVsCoWBvvX8Aq/diY2V8JnbLdHaszV0pST0kzCV0V3CJpSmLaaODfwD6SZko6oKJ20sLMXgMeIHRMN4nwGbi1pEFVQtIIQjcRbeN2/kWpY6qMx1wcHnNxlCpmr8XnnHMulUp9Dco555zLyROUc865VKp2gpJUV9LVkqbHkkBjJB1UiOBqQnxQ96jK58yrrZPi+54u6aRCtFnTyrTMSgNJYyW9FUs0XV7qmCrjMReHx1wcpYq5EHezXQm0ADqa2RJJmwN7FqDdvEmqU+zupCVtClwK9AAMeF3SyHire5oNJfS2OqzEcVTFEmBvM/sm3tn5kqQnzOzVUge2Bh5zcXjMxVGSmKt1BCVpQ+A0YEDmbjwzm2Nm98fp+8cjqjck/VvSRnH8DEmXx/GTJLWL4/dUKG80QdKbkhopGKwfyhcdE+ftK2m0pJHAVIUyR4MljZM0UdIZcT5JulHSO5KeBTarzntOOAB4xsy+jEnpGcqgvlaZllkxM/smvqwbf1J9d4/HXBwec3GUKubqnuLbDvjYzL7OnqDwTNIlwL5m1h0YD/wmMcsXcfxNwPlx3PnAWbHE0R6EUjw/IZQw6kJ4yHSwpBZx/u7Ar81sB+AXwAIz6wn0BE6TtA3hdva2QHtCd++Feii1FfBJ4vXMOM7VgLgDMgH4nLBj8FqJQ6qUx1wcHnNxlCLmmrxJohchKbwc39RJwNaJ6Q/F368DbeLwy8BfJZ0DNImn7XYHRpjZcjObA/yPkIAAxprZh3F4f+DEuK7XCFUptgf6JJb/FHi+0G/U1bz49+sKtAZ2ltSxxCFVymMuDo+5OEoRc3UT1HvAVgpVIbKJkGW7xp/2ZpZ8uCvzgO5y4rUwM7saOBXYgJDY2lWy/m+z1jcgsb5targI6ixgy8Tr1nGcq0FmNh94gTI4nZrhMReHx1wcxYy5WgnKzL4jlBC6XlI9AEnNJR0NvAr0lrRdHN9Q0g5rak/Sj8xskpldA4wD2gGjgWPi4WVzwhHR2ByLPwWcGS/gIWkHSQ2BFxPLtwD2qs57zlrf/pI2USjAun8c5wosfqaaxOENCH2IvV3SoCrhMReHx1wcpYq5EHfxXQL8kXCjwmLCUc0fzGyupP7ACIXSRZl5311DW+dK2ovQ99IU4Ange0I/UW8RLspdYGaf5Ti6up1wqvANSQLmAocDDxMK0E4l1KEbU613G5nZl5KuJCRSgCvMLPU3HyiULOkLNFMo4XSpmd1R2qgq1QK4S6HSfS3gfjNLey/GHnNxeMzFUZKYvdSRc865VPJKEs4551LJE5RzzrlU8gTlnHMulUqSoCRdJmlWomrEhMwdIgVcx0WFbK+CdZRjLb4tJb0gaapCTa1flzqmykhqm/VZ+VrSuaWOa0085uLwmIujVDGX5CYJSZcB35jZtTW4jm/MbKMabH9TQnWMlbX4gJ3SXosv3mrfwszekNSIEPfhZja1xKHlJd5FNAvYxcw+KnU8+fCYi8NjLo5ixpyqU3ySXpXUIfF6lKQe8RmqIQrVdN+UdFic3l/SQ5KejEcxf47jrwY2iJl+eFz+MYVKvJMV6/lVU7nW4pttZm/E4YWE7t7LqUTTPsD75fLPHHnMxeExF0fRYi5lgjovcbj4Qhx3H/BTWGVPfzxwMaF7+J0JD9oOVngIF0KdvmOAToQHcrc0s0HAolhR4nhC4vjUzLqYWUfgyQLEX/a1+CS1AboRSkOVi2OBEaUOooo85uLwmIujaDGXMkFdlyhLlKnucD+Q6avpp8ADcXh/YJBCnb1RQANgqzjtOTNbYGaLCQ/jJuv9ZUwC9pN0jaQ9zGxB4d9OeVGoLP8gcG6uYr9ppFCtpB/w71LHki+PuTg85uIodsypOsVnZrOAeZI6E46K7ouTBByZSGhbmdm0OG1JoomVdf2y2n2XUPl8EvBHSX8oQLhlW4tPoRzUg8BwM3uosvlT5CDgjVg0uFx4zMXhMRdHUWNOVYKK7gMuABqb2cQ47ilgQCxhhKRuebSzVD/U5WsJfGdm/wIGE5JVdZVlLb64De8AppnZX0sdTxUdR/mdDvGYi8NjLo6ixpyWa1AT4vUQCKf1jiWc7su4ktBB1kRJU+Lrytwa5x9OuD41Np4ivJRQO7BaYt29TC2+cZRJLT6gN/BzYO/Etj+41EFVJl5z3I8fumlJPY+5ODzm4ihFzF6LzznnXCql8RSfc8455wnKOedcOnmCcs45l0qeoJxzzqWSJyjnnHOp5AnKOedcKnmCcs45l0qeoJxzzqWSJyjnnHOp5AnKOedcKnmCcs45l0qeoJxzzqWSJyjnnHOp5AnKOedcKlWaoCSZpL8kXp8v6bJKljlcUvsCxJer7dsrazvf9Uv6paQTCxTXUElHVT6nc865fORzBLUE+ImkZlVo93CgRhKUmZ1qZlMLsX4zu9nMhhUkMOeccwWVT4JaRuid9rzsCZLaSHpe0kRJz0naStJuQD9gcOyt9UdZywyVdJOkVyV9IKmvpCGSpkkampjvJknjJU2RdHli/ChJPeLwN5KukvRWbG/zXOuXdJqkcXG+ByVtGJe/TNL5iXavkTRW0ruS9ojja0saHJefKOmMOF6SbpT0jqRngc2qsN2dc85VIt9rUP8AjpfUOGv8DcBdZtYZGA783cxeAUYCA82sq5m9n6O9TYBdCUlvJHAd0AHoJKlrnOdiM+sBdAb2lNQ5RzsNgVfNrAvwInBaBet/yMx6xvmmAb+o4H3WMbOdgXMJXcMT511gZj2BnsBpkrYBjgDaEo7UTgR2q6BN55xzayGvBGVmXwPDgHOyJu0K3BOH7wZ2z3O9/7HQ1/wkYI6ZTTKzFcAUoE2c56eS3gDeJCSvXKfsvgf+G4dfTyybraOk0ZImAcfH9nJ5KEdb+wMnSpoAvAY0BbYH+gAjzGy5mX0KPL+mN+ycc65q6lRh3r8BbwB3FmC9S+LvFYnhzOs68QjlfKCnmX0VT/01yNHO0pjoAJZT8fsZChxuZm9J6g/0rSSuZFsCBpjZU8kZJR1cQRvOOecKIO/bzM3sS+B+Vj099gpwbBw+HhgdhxcCjaoR18bAt8ACSZsDB1Vx+ez1NwJmS6ob46yKp4Az47JI2kFSQ8IpxWPiNaoWwF5VbNc559waVPU5qL8Aybv5BgAnS5oI/Bz4dRx/LzBQ0pvZN0nkw8zeIpzae5twCvHlKjaRvf7fE07PvRzbrIrbganAG5ImA7cQjq4eBqbHacOAMVVs1znn3BrohzNkzjnnXHp4JQnnnHOp5AnKOedcKnmCcs45l0olT1CS6kiaK+nqrPEXVaGNNdbnS1afWIv4+kh6Q9KyZK09SV0ljYmVLiZKOmZt2i8FSU0kPSDp7VjBY9dSx1QZSQfGqh3vSRpU6njy4TEXh8dcHCWJ2cxK+kO4hfxl4H3iTRtx/Dd5Ll87j3lGAT3WMr42hGoWw4CjEuN3ALaPwy2B2UCTUm/PPN/TXcCpcbhe2uMGasfPx7Yx3reA9qWOy2Mu/Y/HvG7HXPIjKOA44HrgY0JlCuLR1Aaxlt7w7AViDb6/SHoL2DVzhBSfSRoqabKkSZLOy1quVpz+x3yDM7MZZjaR8BBxcvy7ZjY9Dn8KfA40r9pbL75YrqoPcAeAmX1vZvNLGlTldgbeM7MPzOx7wmMEh5U4psp4zMXhMRdHSWIuaYKS1ADYF/gPMIKQrDCzQcAiC7X0cj1Y2xB4zcy6mNlLifFdgVZm1tHMOrFq1Ys6hHqB083skgK/j50JexW56g6mzTbAXODO+JzY7fHB4zRrBXySeD0zjkszj7k4PObiKEnMpT6COgR4wcwWAQ8Ch0uqncdyy+P82T4AtpV0g6QDga8T024BJpvZVdUNOilWkbgbONlCPcG0qwN0B24ys26Eih1lcQ7cObd+KXWCOg7YV9IMQoHWpsDeeSy32MyWZ480s6+ALoRrTr8kVIHIeAXYKx61FYSkjYHHCJXXXy1UuzVsJjDTzF6Lrx8gJKw0mwVsmXjdOo5LM4+5ODzm4ihJzCVLUPHLfQ9gKzNrY2ZtgLOIp/mApZn6d1VosxlQy8weBC5h1S/eO4DHgfslVaVIbkXrqkcodzTMzB6obnvFYmafAZ9IahtH7UMo15Rm44DtJW0Tt/uxhC5V0sxjLg6PuThKEnO1v6ir4QjgeTNLVjN/FPizpPqEThInSnqjgutQubQiXFvJJN4LkxPN7K/xJoG7JR2fzyk5ST0JiWgT4FBJl5tZB+CnhJsNmsYK6QD9zWxCnrGW0gBgePygfQCcXOJ41sjMlkk6m1C4tzYwxMymlDisNfKYi8NjLo5Sxey1+JxzzqVSqa9BOeecczl5gnLOOZdK1U5QkupKulrS9FgSaIykqnYwWDTxQd2jKp8zr7ZOiu97uqSTCtFmTZM0RNLnCn1blQWPuTg85uLwmPNXiCOoK4EWQEcz6w4cTvV6062yQtyVtxbr3BS4FNiF8JT1pZI2KXYca2EocGCpg6iioXjMxTAUj7kYhuIx56VaCUrShsBpwIDM3XhmNsfM7o/T949HVG9I+rekjeL4GZIuj+MnSWoXx+8ZyxtNiFUOGikYnChfdEyct6+k0ZJGAlNjmaPBksYpFG89I84nSTcqFDl8FtisOu854QDgGTP7Mj5/9Qxl8KEzsxeBL0sdR1V4zMXhMReHx5y/6h5BbQd8bGZfZ0+IzyRdAuwbj6zGA79JzPJFHH8TcH4cdz5wlpl1JTwjtQj4CaGEURdCWaTBsXoDhOecfm1mOwC/ABaYWU+gJ3CapG0It7O3BdoDJwK7VfM9Z5RjuRLnnCsbNXlqrBchKbwsCUKtujGJ6Q/F368TkhCEquZ/VSgQ+5CZzZS0OzAiVo6YI+l/hAT0NTDWzD6My+4PdE5cX2oMbE94Vimz/KeSnq+B9+qcc67Aqpug3gO2krRxjqMoEU6BHZdjOYDMA7rLM3GY2dWSHgMOJiS2AypZ/7dZ6xtgZk+tEoR0cB7vY23MAvomXrcmlFhyzjlXANU6xWdm3xFKCF0fqxIgqbmko4FXgd6StovjG0raYU3tSfqRmU0ys2sIpTXaAaOBY+I1puaEI6KxORZ/CjhTsTySpB0UqnS/mFi+BbBXdd5z1vr2l7RJvDli/zjOOedcARTiLr5LCN03TI23IP4X+NrM5gL9gRGSJhJO77WrpK1z480QE4GlwBOEMkMTCR1kPQ9cEOvJZbudUFPujRjHLYQjs4eB6XHaMFY9zbjWzOxLwh2M4+LPFXFcqkkaQdgGbSXNlPSLUsdUGY+5ODzm4vCYq7BeL3XknHMujbyShHPOuVTyBOWccy6VPEE555xLpZIkKEmXSZqVqBoxQVKTAq/jokK2V8E6yq4WH6ys5DEpbvfxpY4nI1e9L0lHS5oiaYWkHqWMLxePuTg85uJIW8ylPIK6zsy6Jn7mF7j9Gk1QKt9afBl7xe2epn+SoaxeLmoy4UHuF4seTX6G4jEXw1A85mIYSopiTtUpPkmvSuqQeD1KUo/4DNUQSWMVavQdFqf3l/SQpCfjUcyf4/irgQ3iEcLwuPxjkt6Kt7EfU4Bwy7IWX5rlqvdlZtPM7J0ShVQpj7k4PObiSFvMpUxQ5yVO770Qx91H6Eqd+FBtCzMbD1xM6B5+Z8KDtoPjQ7gQ6vQdA3QiPJC7pZkNAhbFI4TjCYnjUzPrYmYdgScLEH851+Iz4GlJr0s6vdTBOOdcLmk5xZep7nA/kKml91PggTi8PzBI0gRCOaEGwFZx2nNmtsDMFhMext06x7omAftJukbSHma2oPBvp6zsHgv1HgScJalPqQNyzrlsqTrFZ2azgHmSOhOOiu6LkwQcmUhoW5nZtDhtSaKJlXX9stp9l1D5fBLwR0l/KEC4s4AtE69bx3GpF7czZvY5odLGzqWNyDnnVpeqBBXdB1wANDaziXHcU8AAKZRFl9Qtj3aW6oe6fC2B78zsX8BgQrKqrrKsxRevxzXKDBPiLpuePZ1z6xEzK/oPcBnhaGNC4qdNnLY5sAy4NDH/BoTaepOAKcB/4/j+wI2J+f4L9I3D1wDTgOGEGxomxvWMA3oU6H2cQqjo/h5wcim25VrEvC2hruFbcVteXOqYErGNAGYT6jDOJPTxdUQcXgLMAZ4qdZwes8fsMRcnHq/F55xzLpXSeIrPOeec8wTlnHMunTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecSyVPUM4551LJE5RzzrlU8gTlnHMulTxBOeecS6X/B99Op9sAZLGEAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Chi2=3.396389, p=0.065339 for all events secure, exploiting aggregates\n", "Chi2=3.396391, p=0.065339 for all 161 time moments secure\n" ] } ], "source": [ "import sys\n", "from kmsurvival import main\n", "sys.argv[1:] = ['-i2', '--plot-curves']\n", "await main()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Complete Run: 5 logrank tests + survival tables\n", "To run the demo with three parties on localhost, for instance, we add `-M3` as command line option and run [kmsurival.py](kmsurvival.py) outside this notebook using a shell command. The plots are not shown this way, so instead we print the survival tables:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using secure fixed-point numbers: SecFxp64:32\n", "Dataset: aml, with 3-party split, time 1 to 41 (stride 4) months\n", "2024-04-10 18:16:13,332 Logrank test on all events in the clear.\n", "Chi2=2.675902, p=0.101878 for all events in the clear\n", "2024-04-10 18:16:13,348 Start MPyC runtime v0.9.11\n", "2024-04-10 18:16:14,999 All 3 parties connected.\n", "2024-04-10 18:16:15,024 Logrank test on own events in the clear.\n", "Chi2=0.510517, p=0.474915 for own events in the clear\n", " removed observed censored entrance at_risk\n", "event_at \n", "0.0 0 0 0 4 4\n", "3.0 1 1 0 0 4\n", "5.0 1 1 0 0 3\n", "8.0 1 1 0 0 2\n", "12.0 1 1 0 0 1\n", " removed observed censored entrance at_risk\n", "event_at \n", "0.0 0 0 0 4 4\n", "2.0 1 1 0 0 4\n", "3.0 1 1 0 0 3\n", "7.0 1 1 0 0 2\n", "11.0 1 1 0 0 1\n", "2024-04-10 18:16:15,070 Logrank test on aggregated events in the clear.\n", "Chi2=2.685357, p=0.101275 for aggregated events in the clear\n", " removed observed censored entrance at_risk\n", "event_at \n", "0.0 0 0 0 11 11\n", "4.0 3 2 1 0 11\n", "8.0 4 3 1 0 8\n", "12.0 3 2 1 0 4\n", "44.0 1 0 1 0 1\n", " removed observed censored entrance at_risk\n", "event_at \n", "0.0 0 0 0 12 12\n", "4.0 6 5 1 0 12\n", "8.0 3 3 0 0 6\n", "12.0 3 3 0 0 3\n", "2024-04-10 18:16:15,105 Optimized secure logrank test on all individual events.\n", "2024-04-10 18:16:15,105 Interval 1 (time 1 to 4) # observed events = 7\n", "2024-04-10 18:16:15,106 Interval 2 (time 5 to 8) # observed events = 6\n", "2024-04-10 18:16:15,106 Interval 3 (time 9 to 12) # observed events = 5\n", "2024-04-10 18:16:15,107 Interval 4 (time 13 to 16) # observed events = 0\n", "2024-04-10 18:16:15,107 Interval 5 (time 17 to 20) # observed events = 0\n", "2024-04-10 18:16:15,107 Interval 6 (time 21 to 24) # observed events = 0\n", "2024-04-10 18:16:15,107 Interval 7 (time 25 to 28) # observed events = 0\n", "2024-04-10 18:16:15,107 Interval 8 (time 29 to 32) # observed events = 0\n", "2024-04-10 18:16:15,107 Interval 9 (time 33 to 36) # observed events = 0\n", "2024-04-10 18:16:15,107 Interval 10 (time 37 to 40) # observed events = 0\n", "2024-04-10 18:16:15,107 Interval 11 (time 41 to 41) # observed events = 0\n", "Chi2=2.675902, p=0.101878 for all events secure, exploiting aggregates\n", "2024-04-10 18:16:15,685 Secure logrank test for all 41 time moments.\n", "Chi2=2.675903, p=0.101878 for all 41 time moments secure\n", "2024-04-10 18:16:17,651 Stop MPyC -- elapsed time: 0:00:02.651|bytes sent: 1465252\n" ] } ], "source": [ "!python kmsurvival.py -M3 -i2 --print-tables --collapse" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To try out other runs of the demo for yourself, remember to consult MPyC's help message, using the `-H` option:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "usage: kmsurvival.py [-V] [-H] [-h] [-C ini] [-P addr] [-M m] [-I i] [-T t]\n", " [-B b] [--ssl] [-W w] [-L l] [-K k] [--log-level ll]\n", " [--no-log] [--no-async] [--no-barrier] [--no-gmpy2]\n", " [--no-numpy] [--no-uvloop] [--no-prss] [--mix32-64bit]\n", " [--output-windows] [--output-file] [-f F]\n", "\n", "MPyC help:\n", " -V, --VERSION print MPyC version number and exit\n", " -H, --HELP print this help message for MPyC and exit\n", " -h, --help print help message for this MPyC program (if any)\n", "\n", "MPyC configuration:\n", " -C ini, --config ini use ini file, defining all m parties\n", " -P addr use addr=host:port per party (repeat m times)\n", " -M m use m local parties (and run all m, if i is not set)\n", " -I i, --index i set index of this local party to i, 0<=i