From c43453a4619f77bacc33fcc278523bdcadb3e78d Mon Sep 17 00:00:00 2001 From: Tijana Zrnic Date: Tue, 4 Jun 2024 08:01:20 -0700 Subject: [PATCH] ppboot updates --- examples/census_healthcare_ppboot.ipynb | 382 ++++++++++++++++++++++++ ppi_py/baselines.py | 99 +++++- ppi_py/ppi.py | 6 + 3 files changed, 486 insertions(+), 1 deletion(-) create mode 100644 examples/census_healthcare_ppboot.ipynb diff --git a/examples/census_healthcare_ppboot.ipynb b/examples/census_healthcare_ppboot.ipynb new file mode 100644 index 0000000..1f979ba --- /dev/null +++ b/examples/census_healthcare_ppboot.ipynb @@ -0,0 +1,382 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "23c114b7-6751-4192-9939-86d40967caba", + "metadata": {}, + "source": [ + "# Correlation between income and private health insurance\n", + "\n", + "The goal is to investigate the correlation between income and the procurement of private health insurance using US census data. The target of inference is the Pearson correlation coefficient when regressing the binary indicator of health insurance on income. The data from California in the year 2019 is downloaded through the Folktables interface (1). Predictions of health insurance are made by training a gradient boosting tree via XGBoost (2) on the previous year’s data.\n", + "\n", + "Since the basic PPI method is not applicable to this estimation problem, we use PPBoot for prediction-powered inference. We also use the bootstrap for classical inference.\n", + "\n", + "1. F. Ding, M. Hardt, J. Miller, L. Schmidt, “Retiring adult: New datasets for fair machine learning” in Advances in Neural Information Processing Systems 34 (2021), pp. 6478–6490.\n", + "2. T. Chen, C. Guestrin, “XGBoost: A scalable tree boosting system” in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2016), pp. 785–794." + ] + }, + { + "cell_type": "markdown", + "id": "fa0b89de-40f4-4225-ba6f-f35428d8f648", + "metadata": {}, + "source": [ + "### Import necessary packages" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "bec4524b-d6bd-4d3c-ac59-2d6b77ac8a21", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "import os, sys\n", + "\n", + "sys.path.append(os.path.abspath(os.path.join(os.getcwd(), os.pardir)))\n", + "import numpy as np\n", + "import pandas as pd\n", + "from ppi_py.datasets import load_dataset\n", + "from ppi_py import ppboot, classical_bootstrap_ci\n", + "from sklearn.linear_model import LogisticRegression\n", + "from tqdm import tqdm\n", + "from scipy.optimize import brentq\n", + "from utils import *" + ] + }, + { + "cell_type": "markdown", + "id": "5cf90ae6", + "metadata": {}, + "source": [ + "### Import the census healthcare data set\n", + "\n", + "Load the data. The data set contains reported indicators of health insurance (```Y```), predicted indicators of health insurance (```Yhat```), and reported income (```X```)." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a6da3138", + "metadata": {}, + "outputs": [], + "source": [ + "dataset_folder = \"./data/\"\n", + "data = load_dataset(dataset_folder, \"census_healthcare\")\n", + "Y_total = data[\"Y\"]\n", + "Yhat_total = data[\"Yhat\"]\n", + "X_total = data[\"X\"][:,0] # first coordinate is income; second is constant term" + ] + }, + { + "cell_type": "markdown", + "id": "8969f9db", + "metadata": {}, + "source": [ + "### Problem setup\n", + "\n", + "Specify the error level (```alpha```), range of values for the labeled data set size (```ns```), and number of trials (```num_trials```).\n", + "\n", + "Compute the ground-truth value of the estimand." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "5b3c8f29", + "metadata": {}, + "outputs": [], + "source": [ + "alpha = 0.1\n", + "n_total = Y_total.shape[0] # Total number of labeled examples\n", + "ns = np.array([500, 750, 1000, 1500, 2000]).astype(\n", + " int\n", + ") # Test for different numbers of labeled ballots\n", + "num_trials = 100\n", + "# define Pearson correlation coefficient\n", + "def pearson(X, Y):\n", + " return np.corrcoef(X, Y)[0,1]\n", + " \n", + "# Compute ground truth\n", + "true_theta = pearson(X_total, Y_total)" + ] + }, + { + "cell_type": "markdown", + "id": "83ce18be", + "metadata": {}, + "source": [ + "### Construct intervals\n", + "\n", + "Form confidence intervals for all methods and problem parameters. A dataframe with the following columns is formed:\n", + "1. ```method``` (one of ```PPI```, ```Classical```, and ```Imputation```)\n", + "2. ```n``` (labeled data set size, takes values in ```ns```)\n", + "3. ```lower``` (lower endpoint of the confidence interval)\n", + "4. ```upper``` (upper endpoint of the confidence interval)\n", + "5. ```trial``` (index of trial, goes from ```0``` to ```num_trials-1```)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "812f8fd5", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|█████████████████████████████████████████| 100/100 [09:29<00:00, 5.70s/it]\n", + "100%|█████████████████████████████████████████| 100/100 [09:27<00:00, 5.68s/it]\n", + "100%|█████████████████████████████████████████| 100/100 [09:26<00:00, 5.66s/it]\n", + "100%|█████████████████████████████████████████| 100/100 [09:29<00:00, 5.69s/it]\n", + "100%|█████████████████████████████████████████| 100/100 [09:28<00:00, 5.69s/it]\n" + ] + } + ], + "source": [ + "# Run prediction-powered inference and classical inference for many values of n\n", + "results = []\n", + "for i in range(ns.shape[0]):\n", + " for j in tqdm(range(num_trials)):\n", + " # Prediction-Powered Inference\n", + " n = ns[i]\n", + " rand_idx = np.random.permutation(n_total)\n", + " _X, _X_unlabeled = X_total[rand_idx[:n]], X_total[rand_idx[n:]]\n", + " _Y, _Y_unlabeled = Y_total[rand_idx[:n]], Y_total[rand_idx[n:]]\n", + " _Yhat, _Yhat_unlabeled = (\n", + " Yhat_total[rand_idx[:n]],\n", + " Yhat_total[rand_idx[n:]],\n", + " )\n", + "\n", + " ppi_ci = ppboot(\n", + " pearson,\n", + " _Y,\n", + " _Yhat,\n", + " _Yhat_unlabeled,\n", + " _X,\n", + " _X_unlabeled,\n", + " alpha=alpha\n", + " )\n", + "\n", + " # Classical interval\n", + " classical_ci = classical_bootstrap_ci(pearson, _X, _Y, alpha=alpha)\n", + "\n", + " # Append results\n", + " results += [\n", + " pd.DataFrame(\n", + " [\n", + " {\n", + " \"method\": \"PPI\",\n", + " \"n\": n,\n", + " \"lower\": ppi_ci[0],\n", + " \"upper\": ppi_ci[1],\n", + " \"trial\": j,\n", + " }\n", + " ]\n", + " )\n", + " ]\n", + " results += [\n", + " pd.DataFrame(\n", + " [\n", + " {\n", + " \"method\": \"Classical\",\n", + " \"n\": n,\n", + " \"lower\": classical_ci[0],\n", + " \"upper\": classical_ci[1],\n", + " \"trial\": j,\n", + " }\n", + " ]\n", + " )\n", + " ]\n", + "\n", + "# Imputed CI\n", + "imputed_ci = classical_bootstrap_ci(pearson, X_total, (Yhat_total > 0.5).astype(int), alpha=alpha)\n", + "results += [\n", + " pd.DataFrame(\n", + " [\n", + " {\n", + " \"method\": \"Imputation\",\n", + " \"n\": np.nan,\n", + " \"lower\": imputed_ci[0],\n", + " \"upper\": imputed_ci[1],\n", + " \"trial\": 0,\n", + " }\n", + " ]\n", + " )\n", + "]\n", + "\n", + "df = pd.concat(results, axis=0, ignore_index=True)\n", + "df[\"width\"] = df[\"upper\"] - df[\"lower\"]" + ] + }, + { + "cell_type": "markdown", + "id": "d15ba288", + "metadata": {}, + "source": [ + "### Plot results\n", + "\n", + "Plot:\n", + "1. Five randomly chosen intervals from the dataframe for PPI and the classical method, and the imputed interval;\n", + "2. The average interval width for PPI and the classical method, together with a scatterplot of the widths from the five random draws." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "6077b2c4", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3IAAAEcCAYAAACCtobiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABh5ElEQVR4nO3deXxU1fk/8M+dubNPMpNkkskewr5vSqGAFEVRXGvLT9u6QGut0talKopa168t2rr1+7UttFq0VWltqbZU0dYFFGQRBAwkICSBhCSTkHWSmcx6z++PkGmGSchCkplJPu/XKy+Yc+7c+8xlQubJOec5khBCgIiIiIiIiOKGKtoBEBERERERUe8wkSMiIiIiIoozTOSIiIiIiIjiDBM5IiIiIiKiOMNEjoiIiIiIKM4wkSMiIiIiIoozTOSIiIiIiIjiDBM5IiIiIiKiOMNEjoiIiIiIKM4wkSMiIiLqI4/Hg5tuuglWqxUZGRl45plnun3O1q1bMXLkyIh2q9UKSZLCvlpaWvp8HSIa2uRoB0BEREQUr1auXIndu3fjww8/xPHjx7Fs2TLk5eVh6dKlnR5fUFCApUuXQq/Xh7VXVFSgqakJxcXFMBqNoXaTydSn6xDR0CcJIUS0gyAiIiKKNy6XCzabDZs2bcLChQsBAE888QTef/99bN68OeL4tWvX4p577sHIkSPR1NSEY8eOhfref/993HjjjaisrDzr63RGCIHm5mYkJCRAkqRevlIiikWcWklERETUB/v374ff78fcuXNDbfPnz8fOnTuhKErE8Zs2bcIrr7yCn/zkJxF9hYWFGDt2bL9cpzPNzc2wWCxobm7u0fFEFPuYyBERERH1QVVVFWw2G7RabajNbrfD4/Ggrq4u4vi33noL3/jGNzo9V1FREdxuNxYuXIiMjAxceuml+PLLL/t0HQDwer1wOp1hX0Q0tDCRIyIiIuoDt9sNnU4X1tb+2Ov19upchw4dQn19PX7605/iH//4BwwGAxYtWoTm5uY+XWf16tWwWCyhr5ycnF7FQ0Sxj8VOiIiIiPpAr9dHJFLtjzsWLOmJd999F36/H2azGQDw2muvIScnBxs3buzTde6//37cddddocdOp5PJHNEQw0SOiIiIqA+ysrJQW1uLQCAAWW77SOVwOGAwGGC1Wnt1Lp1OFzbqptfrkZ+fj4qKCsybN6/X1zn9fD0hhIDH40EwEIBalqHX61kYhSiGcWolERERUR9Mnz4dGo0GO3bsCLVt3boVs2bNgkrV849YQgiMGjUKL7/8cqjN5XLhyJEjGD9+fL9d50xcLS0oKy3FsaNHcay4GMeOHkVZaSlcp/axI6LYwxE5IiKKa8Lrgvf9/w1r0114OySdKUoR0XBhNBqxbNky3HrrrVi3bh0qKirw9NNPY926dQDaRs0sFgsMBsMZzyNJEi677DI88sgjGDFiBFJTU/HQQw8hOzsbl156KdRq9Rmvc7ZcLS04UVYGv88HnV4PnV6PYDCIZqcTHo8H2bm5MJ2a8klEsYOJHBERxT+fO9oR0DD17LPPYsWKFTj//PNhsVjw2GOPhSpTZmRkYN26dVi+fHm35/nFL34BjUaD73znO2hqasIFF1yAd955B2q1utvrnA0hBGprauD3+WA0mUJTKWVZhtpkgtvlQm1NTVgfEcUGbghORERxTXhd8L7zZFib7tJVHJEj6sDpdMJisaCpqQmJiYmh9tbWVhw7ehSyRhNaf9dRIBBAwO/HiNGjux1ZJKLBxTVyRERERMNUMBBAMBgMjfydTq1WQ1EUBAOBQY6MiLrDqZVERBTfNDpo5n8voo2IuqeWZajVagSDwU5H5ILBIFQqFdSd9BFRdPG7koiI4pqkkqFOzY92GERxSa/Xw2gyodnphPq0dXBCCHg9HiQkJkKv10cxSiLqDKdWEhEREQ1TkiTBlpYGjVYLt8uFQCAAIQQCgQDcLhe0Wi1saWksdEIUg5jIEREREQ1jJrMZ2bm5SEhMRMDvR6vbjYDfj4TERGRx6wGimMWplURERETDnMlshtFkgsfjQTAQgFqWodfrORJHFMOYyBERERERJEniFgNEcYRTK4mIiIiIiOIMR+SIiCiuCa8b3s2/DWvTLVwBSWeMUkREREQDj4kcERHFOQG4GyPbiIiIhjBOrSQiIiIiIoozTOSIiIiIiIjiDKdWEhFRfJN10Hz1+og2IiKioYyJHBERxTVJLUOdPi7aYRAREQ0qTq0kIiIiIiKKM0zkiIiIiIiI4gwTOSIiIiIiojjDRI6IiIiIiCjOsNgJERFF3Y767dEOIS7MSf5qtEMgIqIYwUSOiIiiLkVrg01ri3YYMc8T9ECv1kc7DCIiigFM5IiIKOokACqJs/2JiIh6ij81iYiIiIiI4gwTOSIiIiIiojjDRI6IiIiIiCjOMJEjIiIiIiKKMyx2QkREMUERSrRDICIiihtM5IiIKOpGm8dEOwQiIqK4wqmVREREREREcYYjckREREQEIQQ8Hg+CgQDUsgy9Xg9JkqIdFhF1gYkcERER0TDnamlBbU0N3C4XgsEg1Go1jCYTbGlpMJnN0Q6PiDrBRI6IiOKa8LXCt+2VsDbtvGWQtIYoRUQUX1wtLThRVga/zwedXg+dXo9gMIhmpxMejwfZublM5ohiEBM5IqIYtrN+J9xBd7TDiHECmDgyvKlpO4DYnBJ2fur50Q6B+pHH48GPfvQjbNiwAQaDAffccw/uvvvuMz5n69atuPHGG1FSUhJqE0LgF7/4BdasWYO6ujrMmjUL//d//4eJEycCAPbu3YuZM2eGneecc87B7t27zyp+IQRqa2rg9/lgNJlCUyllWYbaZILb5UJtTU1YHxHFBiZyREQxzCybMdo0OtphUD9qDbbCoOZo4VCxcuVK7N69Gx9++CGOHz+OZcuWIS8vD0uXLu30+IKCAixduhR6vT6sfe3atXj66aexbt06jB07Fr/4xS+wZMkSFBUVwWg0orCwENOnT8emTZtCz9FoNGcdv8fjgdvlgq6T9XCSJEGn18PtcsHj8cBgGJ7vW64dpFjFRI6IKIZJkKCSWGB4SBHRDoD6i8vlwosvvohNmzZh5syZmDlzJg4ePIgXXnih00Ru7dq1uOeeezBy5Eg0NTWF9b388su45557cPnllwMAfvvb3yIpKQnbtm3DRRddhKKiIkyYMAHp6en9+hqCgQCCwSB0pyWW7dRqNXxeL4KBQL9eN15w7SDFMn46ICIiIuqD/fv3w+/3Y+7cuaG2+fPnY+fOnVCUyA3uN23ahFdeeQU/+clPIvqefvppXHfddaHHkiRBCBFK+AoLCzF27Nh+fw1qWYZarUYwGOy0PxgMQqVSQS0Pv9/9t68dbHY6IWs0MJpMkDUaNDudOFFWBldLS7RDpGGOiRwRERFRH1RVVcFms0Gr1Yba7HY7PB4P6urqIo5/66238I1vfKPTc82fPx/Z2dmhxy+++CICgQDmz58PACgqKsK+ffswZcoU5Obm4pZbboHT6ewyNq/XC6fTGfbVGb1eD6PJBK/HAyHCh4uFEPB6PDCaTBFTQYe609cOyrIMSZIgyzKMJhP8Ph9qa2oi7hnRYGIiR0RERNQHbrcbOp0urK39sdfr7fN5d+7cibvvvhsrV65Eeno6/H4/iouL4fP5sG7dOrz00kvYtm0bbrjhhi7PsXr1algsltBXTk5Op8dJkgRbWho0Wi3cLhcCgQCEEAgEAnC7XNBqtbClpQ27NWG9WTtIFC3Db5yciIiIqB/o9fqIhK39sdFo7NM5t2/fjiVLlmDJkiV4/PHHAbQVNamtrYXBYAgVOHnllVdw7rnnorKyEpmZmRHnuf/++3HXXXeFHjudzi6TOZPZjOzc3NBaMJ/XC5VKhYTExGG7FoxrBykeMJEjIophGkkDRUSutaE4NrwGNoa0rKws1NbWIhAIQD61hszhcMBgMMBqtfb6fJs3b8bll1+OxYsXY/369VCp/jtxKjExMezYCRMmAAAqKio6TeR0Ol3EaOGZmMxmGE0mVmc8pePaQbmT9YHDee0gxQ6++4iIYtiYhDHRDoGIujB9+nRoNBrs2LEjtJZt69atmDVrVlgS1hMHDhzAlVdeiSVLlmD9+vVhyUNhYSFmz56NL774Avn5+QCAffv2QZZljB7df9uTSJI0bLcYOF372sFmpxPq0/bQa187mJCYOOzWDlJsYSJHRERE1AdGoxHLli3DrbfeinXr1qGioiK0FxzQNjpnsVh6lBzdcsstyMnJwbPPPova2tpQu8Viwfjx4zF69GjcfPPNeP7559HY2IhbbrkFN998M5KSkgbs9Q1n7WsHO66Vax+h83o8w3btIMUWJnJERBTXhK8Vvp2vh7VpZ38HkpYjCzTwnn32WaxYsQLnn38+LBYLHnvssVBlyoyMDKxbtw7Lly8/4zkcDgc+/fRTAEBubm5YX/vz//nPf+KOO+7AeeedB5VKheuuuw6//OUvB+Q1URuuHaRYJwnWTSUiojgmvC5433kyrE136SpIOlOUIiKKPU6nExaLBU1NTRHr7ejMhBBcO0gxiSNyREQ0ZH148gO0BLhp72Axy2ZckLoo2mEQ9SuuHaRYxUSOiIiGrETZgkkJk6MdxrDiCrhgkjkaSkQ00JjIERFRfJO1kKddHtEGtP0mXSX1rnogERFRPGAiR0REcU1SayCPnB3tMIiIiAYVf01JREREREQUZ5jIERERERERxRkmckRERERERHGGa+SIiGjI0klaKEKJdhhERET9jokcERENWZMtU6IdAhER0YBgIkdERHFN+D3wf/bXsDbNrP8HSaOPUkREREQDj4kcERHFNyUIpfrLiDYiIqKhjMVOiIiIiIiI4gwTOSIiIiIiojgjCSFEtIMgIooXfyx7BQ3+hmiHQR0JAKdXppRUgNTzU9wx6s7+jIgo5jidTlgsFjQ1NSExMTHa4RBRP+AaOSKiXqj11WKCeQLGJ0yIdijUj5r8TbBoLNEOg4iIqMeYyBER9ZIkSVBL6miHQf2Ik1OIiCjecI0cERERERFRnGEiR0REREREFGeYyBEREREREcUZJnJERERERERxhsVOiIh6waa1AZAQFMFoh0L9SJJ6sVcBERFRDOA+ckREFNeE1wXvO0+GtekuXQVJZ4pSRESxh/vIEQ09HJEjIiIiIgghUONwwO/zwZ6ZCY1GE+2QiOgMmMgREREREUqOHEH5sWMAgLJjxzD93HNhNHFkmyhWMZEjIqL4ptZAPf78iDYi6p26kydDf/d5vdi/Zw9mzJoFvcEQxaiIqCtM5IiIKK5JshaaCRdEOwyiuOZqaYFaDv9Y6PV4sPezzzBz9mzodLooRUZEXWEiR0QUp+49cE+0QxhSUnWpWDnmvmiHQTToXC0tOFFWBo1GA41GA7/fH+rzejzYdyqZ45o5otjCRI6IKE5Ns0zHOdZzox3GkFLnq0OKNiXaYRANGiEEamtq4Pf5YEowQ6VXo6G6FkH/f7dYaXW78cWePZh27rmQZX50JIoV/G4kIopTEiSoJXW0wxhauCEPDTMejwdulwtBOYjj7mNoCbogzAoSmkxQKarQcc1OJwr27sXUmTOhVvP/HaJYoOr+ECIiIiIaioKBANx+F074KuAMOKGVNDBqjPBbA1BUStixTQ0NOLBvHxRF6eJsRDSYmMgRERER9ZHH48FNN90Eq9WKjIwMPPPMM90+Z+vWrRg5cmRE+/r16zFq1CgYjUZcffXVqK2tDfUJIbBq1SqkpqYiOTkZ9957b78kVCq1Ck1BJwIBP0xqM2SVDEmSoNbICFoBRQq/RkNdHQq/+ILJHFEMYCJHRERE1EcrV67E7t278eGHH+I3v/kNHnvsMfztb3/r8viCggIsXbo0IhHatWsXbrrpJjzyyCPYsWMHGhoasHz58lD/s88+i9dffx1vvvkmNmzYgNdeew3PPvvsWcffqm5Fs9wCg2KAdFqfJANBYxCKFD7nuLamBocPHoQQnItMFE2S4HchEVFcer38NcxK+kq0wxhSkjXJSNGx2An1jMvlgs1mw6ZNm7Bw4UIAwBNPPIH3338fmzdvjjh+7dq1uOeeezBy5Eg0NTXh2KnNtwHgxhtvhEqlwssvvwwAKC8vR15eHoqLi5Gfn4/c3Fw8/vjjoeTu1VdfxU9/+tOwc5yJ0+mExWJBU1MTEhMTQ+2Vnkp8VPYB0pw2SAEJ0KLt1/wKAB8gZIEGXRPsNTaIYHjymZmdjTETJkCSTk8BiWgwsNgJEVGckiAhKILdH0g9x8+j1Av79++H3+/H3LlzQ23z58/Hz372MyiKApUqfOLTpk2b8Morr8DpdOLRRx8N69uxYwdWrVoVepyTk4Pc3Fzs2LEDOp0O5eXlWLBgQdh1jh8/jqqqKmRkZPT5NehVOsAgoVX2wtCih8qDtqI/EqAYgVazB35NAPn20Tj2xdGwkcTKEyeglmWMHDOGyRxRFDCRIyKKU9/O+U60Q4gJwuuC950nw9p0l66CpDNFKSIaLqqqqmCz2aDVakNtdrsdHo8HdXV1SE1NDTv+rbfeAoDQqNvp58rMzAxrs9vtOHHiBKqqqgAgrN9utwMATpw4cVaJnFWTBLvOjhNKOeyp6VAFJCAIQA0oskC9rwHZuhxkp+bCNM2EA/v2hU2pLD92DGq1GiNGjepzDETUN0zkiIgovqllqEfPi2gjGmhutxs6nS6srf2x1+vtl3N5vV643e6wc/fkOl6vN6zP6XR2epxKUmFK4hQ0+BtQ7XPAqrFCq9HBp3jR6GuESTZjSuIUqCQVUlJTMXHqVBzcvz/sHMeKi6GWZeTk5fXqNRPR2eFPOiIiimuSrINmyiXRDoOGIb1eH5FItT82Go39ci6j0Qi9Xh963PHvZ7rO6tWr8dhjj/Xo2nZ9OhakLECBswDV3mo0+Z2QVTKyDTmYkjgFdn166NhUux3jJ03CoYMHw85RfPgw1Go1MrOze/aCieisMZEjIiIi6oOsrCzU1tYiEAhAlts+UjkcDhgMBlit1l6fy+FwhLU5HA5kZGQgKysr9HjEiBGhvwPoclrl/fffj7vuuiv02Ol0Iicnp8vr2/XpSNWlodHfAI/ihV6lg1WTBJUUWeA8PSsLwWAQRw4dCmv/srAQarUa9rOY6tlOCAGPx4NgIAC1LEOv13MdHtFpuP0AERERUR9Mnz4dGo0GO3bsCLVt3boVs2bNiih00p05c+Zg69atocfl5eUoLy/HnDlzkJmZidzc3LD+rVu3Ijc3t8tETqfTITExMeyrOypJhWRtCjL1mUjWpnSaxLXLys3FyDFjItqLDhxAbU1Nt9c6E1dLC8pKS3Hs6FEcKy7GsaNHUVZaCldLy1mdl2io4YgcERERDTv//ve/sX//fng8noj90B5++OEencNoNGLZsmW49dZbsW7dOlRUVODpp5/GunXrALSNmlksFhgMhm7PtWLFCixcuBBf/epXMWvWLNxxxx24/PLLkZ+fH+q/7777kH1q6uKqVatw99139+Yl97vc/HwEg0EcLyn5b6MQOLh/P6bMnInklN5v5eFqacGJsjL4fT7o9Hro9HoEg0E0O53weDzIzs2FyWzux1dBFL+YyBEREdGwcu+99+KZZ57BtGnTYLFYwvokSepxIge0bdS9YsUKnH/++bBYLHjsscfwjW98A0DbtMd169aFbezdla9+9atYu3YtHn74YdTX12Px4sX4/e9/H+pfuXIlampqcPXVV0OWZdx00034yU9+0uM4B8qIUaMQCARQUVYWahNC4MDevZh6zjmwJiX1+FxCCNTW1MDv88FoMoWmUsqyDLXJBLfLhdqamrA+ouGMG4ITERHRsJKUlIS1a9fimmuuiXYog6arDcH7gxAChwsL4aioCGtXyzKmn3suEnp4vdbWVhw7ehSyRhNac9hRIBBAwO/HiNGjezTKSTTUcUSOiIjimvB7ETj477A2edJiSBpdF8+g4U6WZcycOTPaYQwZkiRh3MSJUIJB1HQo2BIMBLB/zx7MmDWrR9Mhg4EAgsEgdKcqc55OrVbD5/UiGAj0W+xE8YzFToiIKL4pAQRLd4V9QeEHPeraj3/8YzzyyCPweDzRDmXIkCQJ4ydPRsppm6AH/H7s37MHraf2wjsTtSxDrVYjGAx22h8MBqFSqaDuZLSOaDjidwIRERENefn5+aF1VYqioKysDBs2bIDdbodarQ47tqRj8Q7qMZVKhYlTp6Jg71401teH2n1eL/bt3o0ZX/lKaB+8zuj1ehhNJjQ7nVCftg5OCAGvx4OExMQznoNoOGEiR0REREPeo48+Gu0QhgW1Wo3J06fji88/h7OxMdTu9Xiwf/duzJg1C1pd59OeJUmCLS0NHo8HbpcLOr0+NELn9Xig1WphS0tjoROiU5jIERFRfFPJUOd/JaKNqKNly5aF/v7444/jnnvugdFoDDvG6XTiscceG+zQhhxZljFlxgzs370bLc3NofZWtxv79+zB9FmzoNFoOn2uyWxGdm4uamtq4Ha54PN6oVKpkJCYCFtaGrceIOqAVSuJiIhoyDt8+DCqq6sBAOeffz7+/ve/I+m00vgHDhzAypUr4XK5ohHigBrIqpVd8fl82PfZZ3Cfdj8TEhMx7dxzO61M2U4IAY/Hg2AgALUsQ6/XcySO6DT8lSURERENeZWVlVi0aFHo8dVXXx1xjMlkwp133jmIUQ1tWq0W0845B3s/+wye1tZQe7PTiYK9ezF15syI9YntJEniFgNE3eCIHBEREQ0r+fn5+Oyzz2Cz2aIdyqCJxohcu1a3G3s/+ww+rzesPdlmw+Tp06FSsYg6UV/wO4eIiIiGldLS0mGVxEWbwWjEtHPPjVgXV19bi8IvvoCiKFGKjCi+cUSOiIiIhryO2w90ZyhuPxDNEbl2Lc3N2PfZZwictqG3PSMD4ydP5ho4ol7iGjkiIiIa8jpuP1BcXIznn38eK1aswKxZs6DVavH555/jhRde4Bq5AWROSMCUmTOxf88eKB02/a6uqoJaljFm/Hgmc0S9wBE5IiKKayLgRaDoo7A2ecL5kOTO96oiOvfcc3Hffffh//2//xfW/o9//AMPPvggDhw4EKXIBk4sjMi1a6ivR8Hnn0dMqcwZMQIjx4xhMkfUQxyRIyKi+BYMIHh0W1iTPPY8gIkcdeHQoUOYMmVKRPvIkSNx/PjxKEQ0vCQlJ2PStGk4sG8fOo4nlB87BlmWkTdyZBSjI4ofLHZCREREw8p5552HO++8ExUVFaG2kpIS3Hbbbbj44oujGNnwkZKaigmdJNOlR4/iBJPpuCCEQGtrK1qam9Ha2gpO8ht8HJEjIiKiYeUPf/gDvvnNbyI3NxfJyckQQqChoQEXXHABfv/730c7vKhRhIJGfwM8ihd6lQ5WTRJU0sD9zj8tPR3BYBCHDx4Maz96+DDUajUysrMH7Np0dlwtLaitqYHb5UIwGIRarYbRZIItLQ0mszna4Q0bTOSIiCi+qdRQ5c6IaKuprkbA749OTNStzCh+SM/IyMCnn36KgwcPoqioCAAwefJkjB8/PmoxRVu1x4ECZwGqvdUIKH7IKg3sOjumJE6BXZ8+YNfNyMpCMBjE0UOHwtoPFxZCpVbDnpExYNemvnG1tOBEWRn8Ph90ej10ej2CwSCanU54PB5k5+YymRskLHZCRERDUuWJEzAYjTAYjdEOhTohyzJkefB+n1xWVoacnBxIkoSysrIzHpubmztIUQ2eMxU7qfY48HHdx3AFWmDVJEGr0sKn+NDob4BJNmNByoIBTeYA4HhJCUqPHg1rkyQJk6ZNgy0tbUCvTT0nhEBZaSmanU4YTaawwjRCCLhdLiQkJiK3F9t9UN9xRI6IiIY0fpggABgxYgQcDgfS0tIwYsQISJIUtqan/bEkSQh2KI0/1ClCQYGzAK5AC+y69ND3i16th12Vjmpv20hdqi5tQKdZ5o0ciWAwiLLS0lCbEAIHv/gC4ydNQqrdDpWKpR2izePxwO1yQafXR/zfKkkSdHo93C4XPB4PDAZDlKIcPpjIERER0ZBXWlqK1NRUAEBWVhZWrFiBSy65BMnJycM62W/0N6DaWw2rJqnTD+ZWjRXV3mo0+huQrE0Z0FjyR49GMBhERYcRU6EoKCoowNFDh5CWkQF7RgYSEhOH9b9ZNAUDAQSDQej0+k771Wo1fF4vgqdt+k4Dg4kcERERDXl5eXmhvz/88MP4z3/+g+eeew5qtRqLFy/GJZdcgosvvhgpKQObrMQaj+JFQPFDq9F22q9V6dDkd8KjeAc8FkmSMHrcOAQDATgqK8P6/H4/KsrKUFFWBqPJhPTMTNgzMrpMKGhgqGUZarUawWCw06nRwWAQKpUK6kGcNj2ccYyaiIiIhpWbb74Zb7zxBmpqavDOO+9gypQpWLduHTIyMjB79uxohzeo9CodZJUGPsXXab9P8UJWydCrBmdfRkmSMG7SJKSld70mz+1yoeTIEWz/+GPs37MH1VVVw2o6bDTp9XoYTSZ4PZ6I7QaEEPB6PDCaTNAzwR4UTJeJiIho2AkGg/j888+xdetW7NixAwUFBdDr9bBardEObVBZNUmw6+w40VoOuyo9onhFo78R2YYcWDVJgxaTJEmYMGUKkm02VJ44AWdjY5fHNtTVoaGuDmq1Gqnp6UjPzITFauXUywEiSRJsaWlha+XaR+i8Hg+0Wi1saWm8/4OEVSuJiCiuiYAPgSNbw9rkMfNxsq4BGo2GVStj1GBXrexo4cKF2L17N5KSkjBnzhzMmzcP5513HmbMmDFkC2r0vGqlFVqVDj7Fi0Z/46BVrTwTt9uN6spKOCor4fV4uj1ebzDAnpGB9MxMfv8PkI77yCmKApVKxX3kooAjckREFN+CfgQPfRTWJI+cjTS7PUoBUazTaDRQqVSw2WzIzMxEVlYWsrKyhmwS1x27Ph0LUhaE9pFr8jshq2RkG3IGfB+5njAajcgfPRojRo1CU0MDHJWVOFld3eV0Sk9rK46XlOB4SQksVivSMzORardD1mgGOfKhy2Q2w2gywePxIBgIQC3L0HdSyZIGFkfkiIgorgmvC953ngxr0126CpLOFKWIKB4EAgHs2bMHH3/8MT7++GPs2LEDFosF5513HtatWxft8PrdmUbk2ilCQaO/AR7FC71KB6smaUC3HDgbwUAAJ2tqUF1VhYa6um6PV6lUsKWlIT0zE0kpKUw4aEjgiBwRERENO7IsY/bs2UhISIDRaIRer8fGjRvx4YcfRju0qFFJqgHfYqC/qGUZ6ZmZSM/MhMfjQXVVFaorK+F2uTo9XlEU1DgcqHE4oNXpQlMvOQ2Q4hlH5IiIKK4Jvwf+z98Ka9PM/DokDaumUed+85vfYPPmzfj444/hcrmwYMECLF68GIsXL8aECROiHd6A6MmIXLwTQqDZ6YSjshI1DgcCfn+3zzEnJiI9MxNp6enQajvfgoEoVjGRIyIiomFl+vTpuPjii7F48WKcd955w+ID/HBI5DpSFAV1J0/CUVmJ+traiFL5p5MkCck2G9IzM5GSmjps10tSfGEiR0REUVF27Fi0Q6BuaDQaZGRlRTsM6gc9SeSEECh1l6A50IyJCZOgUQ2N4iA+nw81VVVwVFWhxens9nhZo4E9PR32zEwkJCZyPR3FLCZyREQUFbUnT3J9SoxztbTAlpoa7TCoH/Qkkfvg5Pt42/EvAECyJgU3j7g56hUr+1tLc3PberqqKvi83m6PN5pMSM/MhD0jAzpuck0xhokcERFFRV1tLRO5GOdqaUGKzRbtMKgf9CSR++WRp1DlqQo91qv0+G7e9zDGPHawwhw0iqKgsb4ejspK1NbUQFGUbp+TlJKC9MxM2NLSoFarByFKojPjBGAiIiKiPvJ4PLjppptgtVqRkZGBZ555pstj9+7di9mzZ8NoNGLWrFnYs2dPqE+SpE6//vjHPwIA3nzzzYi+pUuX9utrydKHT6P1KB6sLV2DnfU7+vU6sUClUiHZZsPEqVMx92tfw9iJE2GxWs/4nIa6OhQVFODTLVtwuLAQTY2N3a69IxpI3H6AiIiIqI9WrlyJ3bt348MPP8Tx48exbNky5OXlRSRZLpcLl156Ka677jq8/PLLWLNmDS677DIUFxfDZDKhqqoq7PjnnnsOf/nLX3DVVVcBAAoLC3HFFVfgd7/7XegYfT9P9bsy4+uo8dagrLUs1KZAwV8q/oxaXy2W2C+N2X3lzoas0SAzOxuZ2dlodbvhqKxEdVUVPK2tnR4fDARQdeIEqk6cgMFoRHpWFtI59ZKigFMriYgoKji1MvZxauWZuVwu2Gw2bNq0CQsXLgQAPPHEE3j//fexefPmsGP/8Ic/4IknnkBxcTEkSYIQAmPHjsWDDz6I5cuXhx1bWlqKiRMnYuPGjbjwwgsBANdffz1yc3Px85//vE+x9rRqpU/x4fXy1/CFc39E33TLdHwr+zvQqoZ+lU8hBJoaG3Hi+HHUnTzZo5G35JQUpGdlISU1lVMvaVAMvV+rEBEREQ2C/fv3w+/3Y+7cuaG2+fPnY+fOnRFrrnbs2IH58+eHKiBKkoR58+Zh+/btEed9+OGHsWjRolASB7SNyI0dO/Br1bQqLW7MXYbzbRdE9O1r2oc1pb9BS6BlwOOINkmSoNFooNFqkWyzISklpdsR0Pq6OhR+8QW2b9mCLwsL4Wxq4tRLGlBM5IiIKCqEEPyK8S/upXVmVVVVsNlsYfvQ2e12eDwe1NXVRRybmZkZ1ma323HixImwtrKyMrz++ut46KGHQm1CCBw+fBjvvfcexo4di1GjRmHVqlXw+Xxdxub1euF0OsO+ekolqXBFxpX4f1nXQnXaR8Vj7mN4/uhzqPY4eny+eCSEQG1NDfw+H4xmEzRmLYxpCUjJSoMlyQrVGUbcAoEAKk+cwOc7d+KzTz9FWWkpvD2okEnUW1wjR0REUdFfZe2F1wXvO0+GtekuXQVJZ+qX8w9nBoMh2iHENLfbDZ1OF9bW/vj0D+5dHXv6cS+99BLOPfdczJ49O9RWVlYWev4bb7yB0tJS3H777WhtbcWvfvWrTmNbvXo1HnvssT6/NgD4avJXkaxJxitl6+BRPKH2en8d/rf4V1ie990hWdESaCti43a5EJSDOO4+hpagC4oIQiWpYdaakJyWDMknQdZoUH/yJILBYKfncbtcKDlyBCVHjoQ2HLelpfGXJNQvmMgREVF8U6mhso+NaCMaaHq9PiIRa39sNBp7dOzpx/3tb3/DrbfeGtaWl5eHuro6JCUlQZIkTJ8+HYqi4Prrr8ezzz7b6Xqs+++/H3fddVfosdPpRE5OTq9f47iEcbht1B148djv0OBvCLW3Kq1YW7oG12Rdi68kzz7DGeJTMBCA2++CAzXwCx/0Kj3UKj2CIghnwIlWyYN0yY5R+WMxfuJEnKypgaOyEo319V2es762FvW1tZBlGWkZGUjnhuN0lpjIERFRXJM0emjn3hDtMGgYysrKQm1tLQKBAGS57SOVw+GAwWCA9bRS9llZWXA4wqcjOhwOZGRkhB6Xl5ejsLAwVKmyo+Tk5LDHEyZMgMfjQX19PVI7Gd3W6XQRI4B9laHPwB2jfoKXjr+I8tMqWv65Yj1qfbW4xL5kSFW0VKlVaAo6ERB+mLRmtOdasiTDJJnh9rnQJDVBpVZBLctIz8xEemYmPK2tcFRWwlFZ2WXVy0AggMryclSWl8NkNiM9MxNpGRn99u9Fw8fQ+Y4jIiIiGkTTp0+HRqPBjh3/3Wdt69atmDVrVsTUuTlz5uDTTz8NFb8QQmDbtm2YM2dO6JidO3ciJycHubm5Yc997733kJKSArfbHWrbt28fUlJSOk3iBkKiJhE/GvljTEmcGtH3/sn/4NXyP8Gv+PvteopQUO+rQ6WnEvW+Oiii+w27+1OruhXNcgsMigGnj5dJAAyKHs1yC1rV4cma3mDAiFGjMHv+fEyfNQvpmZlnXE/namlB8ZdfYvvHH6Ng716crK7u0ebkRABH5IiIqAeKDx+Odgg0BIwaNy7aIfQro9GIZcuW4dZbb8W6detQUVGBp59+GuvWrQPQNuJmsVhgMBiwdOlSrFq1CnfeeSduueUWrF27Fi6XC9dcc03ofAcOHMDEiRMjrjN37lwYDAZ8//vfxyOPPIKSkhKsXLkS995776C9VqCtouWy3OV42/EvfFT7YVjfvqa9aPQ34Ht534dZPrttRao9DhQ0fYHalloowSBUajVsZhumWKbCrk8/q3P3lFf40GJywRg0AB4AWrQNfygAfIAkq9BicsMrOi84I0kSrElJsCYlYfT48aitrm6betnQ0OnxEAJ1J0+i7uRJyBoN7KemXpoTEjj1krrEfeSIiKhbjspKJFgs0Q6D4pxWq4VGo4l2GP3K7XZjxYoV2LBhAywWC1auXIk777wTQNuH+XXr1oX2idu1axduvfVWFBUVYerUqVizZg1mzJgROteKFSvQ2NiI9evXR1zn4MGDuPPOO7Fjxw4kJCTglltuwcMPP9zjD/k93Ueupz6t24a/V26AgvDRoxRtCr6f9wPY9fY+nbfa48C2yq2QGhWY/CaohBqKFIRL44KwqjAvc/6gJHP1vjpsqt6EBL8ZhhY9VB4AAoAEKHqg1dyKZo0LS+xLkKxN6fF5W91uVFdVwVFRAY/H0+3xJrMZ6VlZsKenQ8upl3QaJnJERNSt6qoqJnJ01mRZDivVT4OnvxM5ADjUXIRXyl6GVwkv4mJQGbA873sYYx7Tq/MpQsGH5e/DXdUCozBGjIK5JTeMGWZckHPhgK/HU4SCD09+gBOt5bBr06EKSEAQgBpQZIFqnwPZhhxckLqoT7EIIdDY0ABHZSVOOhzdTqeUJKmt6mVWFlJsNla9JABM5IiIqAeYyFF/YCIXPQORyAFAlacqoqIlAKigwjXZ1+IrST2vaFnnrcXWwo9h8Oqh0qsQtjhNAIoniFadF/MnLkCKztZPr6Br1R4HPq77GK5AC6waK7QqHXyKF43+RphkMxakLOiX0cFAIICTp6ZeNnU19bIDTfvUy6wsmBMSzvr6FL+YyBERUbeYyFF/YCIXPQOVyAGA0+/ES8d/j/LW8oi+C1Mv6nFFy+ONx1BQtA86rQ6SHDllVAQEvD4vpkyYjjzriP4IvVvVHgcKnAWo9lYjoAQgq2TYdXZMSZwyIFM8W93uUNVLbw+mXpoTEkJVL/m9NfwwkSMiom4xkaP+wEQuegYykQMAn+LDa+WvosD5RUTfDMsMfCv7O9Cozrw+sqL+BPYf/hwqnQqyOrIeXyAYgOJVMG3cTGQlZ/db7N1RhIJGfwM8ihd6lQ5WTdKAT+0UQqCxvr5t6mUPKllKkoSU1FSkZ2YimVMvhw1WrSQiIiKis9Je0fJfjo3YXPtRWN/epr1o8Dfie3k3nbGiZZI+CXpZj+ZAM9QqGR3ruAgBeANeJMgJSNInDdTLiBmSJCEpJQVJKSkYM2ECahwOOCor4Wxs7PR4IQRqa2pQW1MDjVYbVvWShi6OyBERUbeqq6pgHoDf4tPwMhSrVsaLgR6R6+hMFS1vHvEDpOk6r2gphEDhkS9QXlcOj8YDvVoPtaRGUAThCXqg9+uRk5KDiWOmDlpJ/vCplX7IKs2ATq3sjtvlgqOyEtWVlfB6vd0eb05MRHpmJuwZGfzeG4KYyBERUVwTvlb4dr4e1qad/R1IWkOUIiKKPT1J5PpzCmGXFS3VRnw397sY3UVFS1dLC46UHkK9ux4ulRtBKQi1UMOkGJFsTMaY/PEwmc9un7qeCi92kgStSguf4kOjv6Ffi530hRACDfX1cFRUoLampkdTL21paUjPzERSSgqnXg4RTOSIiIiIhrjuErmBGHmq9FTixWO/Q6O/MaxdLalxTda1mJX0lU6f52ppwcnqajS2NCAQDEBWy7AmJCE1zT5oSVzY9gO69LARQCEEqr1nt/1Af/L7/W1VLysq4Gxq6vZ4rVYLe2Ym0jMzB+1+0sBgIkdEREQ0xJ0pkRvIkSenvwkvHn8RJzqpaHlR2mJckrak02mSQgh4PB4EAwGoZRl6vX7QplMC/90Q3KQ2Qa/WR/R7gq1wBd293hB8oLlcLlSfqnrp68HUy4TERKRnZSEtPb3XUy+jUQSGwjGRIyIiIhriukrkBmPkyat48Vr5qzjgLIjom2GZiW9lf7vbipaDrdJTif9Uv4dUXVqnr1sRCk56T+Ii+2Jk6jOjEOGZCSFQX1cHR2UlamtqILqbeqlShaZeJqekdJs0x9raweGKVSuJiIiIhqlGfwOqvdWwapIiPrxLkgSrxopqbzUa/Q19HnnSqXRYnvvdLipafo5GfwO+201Fy8GmV+kgqzTwKb5OR+R8iheySoZepYtCdN2TJAkpNhtSbDb4/f62qpcVFWh2Ojs9XigKTjocOOlwQKvTtRVIycyEyWSKODZiBFfTNoJ7orUcDf6GqK4dHG6YyBERDRGFX3wBT2trtMOgGKM3GDBx6tRoh0ExyqN4EVD80Go6399Pq9Khye+ER+l+mt6ZqCQVrsy4CinaFLxZ+fewipal7lL8qvh53Dzi5i4rWg42qyYJdp29baRSSocqIAFBAGpAkQUa/Y3INuTAqon9rRA0Gg2ycnKQlZMDV0tLqOqlz+fr9Hif14uy0lKUlZYi0WJBelYWUu12aDQaKEJBgbMArkBL2AiuXq2HXZWOam/bSF1XI5nUv5jIERENEa1uN5JtNiTZbNEOhWKM1+uFThebIwcUXYM98jQvZT6StSn442kVLet8tfhV8a/w3dzvYbR5dL9c62yoJBWmJE6Bs6UJnqoWmPwmqIQaihSER+NCojURUxKnxF2yYjKbMWrsWOSPHo2GjlMvu1hp5WxqgrOpCUcPHYItLQ2mtARU+wZ2BJd6jokcEdFQIkksK00RuByeuhI28qSKXCM3ECNPExIm4LaRd+DF4+EVLVuDbqw99tszVrQcTOaAGSNdI1Dvb98KwQO1UMPityLZlQxzIHamgvaWSqVCSmoqUlJT4ff5UONwoKqyEi1dTL1UFAU1DgfgcMAuJ0MkS1CSAXFaft9fI7jUM0zkiIiIiIap9pGnBn8Dqr0OWDVWaFU6+BQvGv2NMMnmARl5yjRk4s5RP4moaBkUQaw/8TrqfHW4OO2SQa1U2ZEQArU1NZCDMvKSR8CreBEQAciSDJ1Kh1aXG7U1NTCaTFGLsb9otFpk5eYiKzcXLc3NbVMvq6rg72LqpRyQgRoANYBiFAgmCwQTBaCJ/bWDQw2rVhIRDRF7duxA8qnfsBJ1pNVqoddHTpuj4aN3+8gFIKvkQalC2FbR8k844DwQ0TcpYTIuSlsMs2wa9NL2ra2tOHb0KGSNBrIcOe4RCAQQ8PsxYvRoGAyGQYtrsCiKgvraWjgqK1F38mSPRvUVrYDL4IbRYsKc3HkwGoxxn+TGOiZyRERDBBM56goTOeoukQOity+YIhRsrPoHttRtiegzqo0YZRqNHEPOoJa2b2luxrHi4i5H3IQQaHW7kTdyJMwJCYMSU7T4fD7UVFXBUVmJlubmHj9Pp9PBkpQEa1ISLElJQ2L0MtYwkSMiGiKYyFFXmMhRTxK5aNtWtxV/r9wAgciPpgaVASnaFCxKvRBTLFMHPMEc7iNyXWlpboajogJVVRUI+oO9eq5Go4HlVFJnTUqCOSGBid1ZYiJHRDREFH7xBQwmE5JTWCmMOpAk6PV6Vq0c5uIhkVOEgj+fWI/PG/eEbU9wukQ5ERMTJmFS4iSMMY+FVtX51glnQwiBstJSNDudESNJQgi4XS4kJCYiNz9/WCYjiqKg9mQNTlSUoaXBCSV45g3HO6OWZVis1lBil5CYyGJdvcREjoiIiGiIi4dErt5Xh03VmyCEwMHmg/Aqnm6fI0sajDWPwcSESZiYOAlWjbXf4nG1tOBEWRn8Ph90ej3UajWCwSC8Hg+0pwqEmMzxW7myvwgh0NLcjKaGBjQ2NKCpoQF+v7/X51GpVEi0WEKJXaLVCrVaPQARDx1M5IiIiIiGuHhI5Co9lfhP9XtI1aXBr/hx1HWkrfiKCPT4HFn6LExMnIRJCZORbcg+6ymYrpYW1NbUwO1yQVEUqFQqGE2mtj3VhnkSJ4SAx+NBMBCAWpah1+shSVLbiKXb/d/Err4eXm/vtyOQJAkJiYlhiZ1GoxmAVxK/mMgRERERDXHxkMi1j8iZ1KbQ5uRtBVgacdJbgxpvDVqV1h6fL0FOxMSEiaEpmLo+lsTvKmEZzjomuMFgEGq1ussEt/3+dRyxa3W7+3Rdc0LCfwuoWK3QDvMp40zkiIiIiIa4eEjkFKHgw5MftG1OrovcnLxtn7skJGosKGouRKmr5Ixr6TqSJRmjTWMwKXESJiZMQpK2/zY4H27ONOVUo9UiuwdTTr1eb1hi52pp6VMsRpMprDLmcCvqxESOiIiIaIiLh0QOaNvP7uO6j+EKtHS6OfmClAWhLQjcQTcONReh0FmIopYitAZ7PsqTqc8MFUzJMeQO6h518WygisD4/X40nUrqGhsa0NzcDPQhRdHr9WGJncE4tPeyYyJHRNQHn336Kdx9nBpC/U1E/sCXJACD88PbaDRi1ty5g3Itor6Kl0QO6Nvm5EERRKmrFIXNB1HYfBA13poeX88smzExYSImJkzGOPM46NTDe7remQzWtgyBQADOpqZQYudsaoJQel8ZU6vVhm15YDKbBzyxG8ypuEzkiIj6oLCgAPaMjGiHQTHCaDINq72kKP7EUyIHnP3m5Ce9J1HYfBAHnQdR4iru8RRMtaQOm4KZrE3u60sYkqK1UXowGESz0/nf6ZiNjVCCvdvHDgBkWY7Yy64/tzzozdrB/sBEjoioD4oOHEB6Zma0w6AYodfrYTAaox0GUZfiLZHrT61BNw41H8LB5oMoau7dFMwMfQYmJUzGxISJyDXmDfspmLGyUbqiKBFbHgQCPa9u2k6lVv93LzurFQkWS5+3POiPtYO9xUSOiKgPmMhRR0zkKNYN50Suo6AI4rj7GA4626ZgVnure/xcs9qMUaZRGHnqK0OfMewSu1jdKF0IAVdLS1ti19iIpvp6+Hy+Xp9HkqSIvew6S1g7u3407gsTOSKiPmAiRx0xkaNYx0Suc7XeWhxsPoBCZyGKXUd7PAUTAPQqPUaaRiLfNBKjjKOQbciBrOr+Q3+8i4eN0oUQaG1tDRux87T2fOuKEElCwmlbHmi02ojDojVSyUSOiKgPmMhRR0zkKNbFWyIXjb3bWoOtONxyGIXOAyhqLoIr6OrV8zWSBnnGERhpGolRplHIM46AVhX5oX8oiMeN0k/fy87t6t2/bzuT2RyW2On0+rC1g0Db9grta+R0p/a6G4i1g0zkiIj6gIkcdcREjmJdPCVyg10wojOKUNqmYDYfRKHzIBxeR6/PoYIKOYac0FTMfFM+jOqh8/+EoihwNjXB5/VCq9Mh0WLp18IhA83n9aKpsTGU2LU0N/fpPAajEeaEBHg9HqhUKrjdbgQDAQghIEkSZI0GZrMZskbDETkiolhw9PBhJKWkRDsMigESAAOrVlKMi5dELhoFI3qizleH4pajKHYXo8RVgjpfba/PIUFCuj6jbZ2dcSRGmkYiUWMZgGgHXiwk2/0t4PejqbExlNw1NzWhr2mSJElQqVSQJAmKogCSBFtaGiZMnsw1ckRERESxwOPx4Ec/+hE2bNgAg8GAe+65B3fffXenx+7duxe33norCgoKMGnSJKxZswbnnHNOqN9qtaKpqSnsOc3NzTCbzb26TmfiIZGL1UIanWn0N6LUVYISVwmKXcVweKv6dJ5UbWrbGrtTo3bJmuSov7buxGqy3d+CwWD4XnaNjW1JWS+1j1KaExIwfdasfh21HPorMomIiIgGyMqVK7F79258+OGHOH78OJYtW4a8vDwsXbo07DiXy4VLL70U1113HV5++WWsWbMGl112GYqLi2EymVBRUYGmpiYUFxfD2GGarunUmpueXieeeTweuF0u6DpZDydJEnR6PdwuFzweT9RHwK0aK2ZYZ2KGdSYAwBVwodRdihJXMUpcxTjReqJHhVNO+k7ipO8kdjXsBABYZMupqZhtyV2azh5TlTGFEKitqYHf5wtLtmVZhtpkgtvlQm1NTZf7zMUTtVqNpORkJCW37SWoKErEXnbBHmx5oFKrIWs08Hq9cDY1wZqU1G8xckSOiIiIqA9cLhdsNhs2bdqEhQsXAgCeeOIJvP/++9i8eXPYsX/4wx/wxBNPoLi4GJIkQQiBsWPH4sEHH8Ty5cvx/vvv48Ybb0RlZeVZXacr8TAiF63NpgeCN+jFMfcxlLjbErvj7jIEhL/X5zGpTcg35WOkcRRGmUYh05AFtdS3fc76Q6zsIxcLhBChveyqq6rQ7HR2epzp1KbjXo8Ho8eNQ1p6er/FwBE5IiKKayIYgHKyOKxNlToKkpo/4mhg7d+/H36/H3Pnzg21zZ8/Hz/72c9Clfza7dixA/Pnzw8lKJIkYd68edi+fTuWL1+OwsJCjB079qyvE8/UshyaptdZkhAMBqFSqaDuwb5e0aZT6zAuYRzGJYwDAASUAMpby1DsaltjV+ougVfxdnseV9CFA84DOOA80HZelQ4jjCNCBVRyDbnQqDQD+lo6CgYCCAaD0On1nfar1Wr4vN4ejVTFO0mSkJCYiITERJgTE1H4xRdtm4kLgUAgAAFAq9VCVqsRCAQgSRK0pypY9pfY/04gIiI6k4AX/u2vhjXpLl0FMJGjAVZVVQWbzQZth32l7HY7PB4P6urqkJqaGnbspEmTwp5vt9tx4EDbB/SioiK43W4sXLgQhw8fxowZM/D8889j7NixvbpOO6/XC6/3v4mCs4vRglii1+thNJnQ7HRC3ckaOa/Hg4TEROi7SCJimaySkX9qzzmgrSpmpacitMau1FWClmBLt+fxKl4cbjmMwy2HAQBqSY1cQx5GmUa1nd+YD7164O7PUEq2+1OixQLjqamlOr0+LGFThIDf74fRZEKipX+L2wyvu0xERETUT9xud2iPqHbtjzsmUWc6tv24Q4cOob6+Hj//+c+RmJiIp556CosWLUJhYWGvrtNu9erVeOyxx/r+4qJAOlXZr+NaudM3m7alpcX92isAUEkqZBtykG3IwQLb1yCEQI235tRUzBKUuIrR4G/o9jxBEUSpu22EDyfbKmNmGbJPVcVsW2tnlvuv8MhQTrbPhkqlQnZuLkqOHGkr+qLRQKVWQwkG4ff7IcsysnNz+330nIkcERERUR/o9fqIRKr9sfG0fQW7Orb9uHfffRd+vx/mU9X+XnvtNeTk5GDjxo29uk67+++/H3fddVfosdPpRE5OTm9f4qAzmc3Izs0Nlbb3eb1QqVRISEyM69L23ZEkCXa9HXa9HV9NbptCW++rDyV1Je5i1Hhruj2PgMCJ1nKcaC3Hx3VbAAB2nb0tqTuV3CVp+15sYzgl271lS0sDAJwoK0Or2w3h90OSJBhNJmTn5ob6+xMTOSIiinMSYLRGthENsKysLNTW1iIQCISmmTkcDhgMBlit1ohjHY7wTaUdDgcyMjIAtI2wdRx10+v1yM/PR0VFBebNm9fj67Q7/XzxxGQ2w2gywePxIBgIQC3L0HdSyXKoS9YmI1mbjHOTzgUANAeaUXpqKmaJqwSVngoIdF+zsNpbjWpvNbbXf9p2Xk0yRpr+O2KXqu1d4tWebJ+srkZjSwMCwQBktQxrYhJS0+xDNtnuCVtaGpJttkHbKJ2JHBERxTVJZ4T+4p7vp0XUX6ZPnw6NRhMqZAIAW7duxaxO9oqaM2cOnnzySQghQlUrt23bhgcffBBCCIwePRoPPfQQli9fDqCtUuWRI0cwfvz4Xl1nqBAQaFW74ZG80Kt00EEHaZj/giZBTsBUyzRMtUwDALQGW3HMXRpaY1fWWoagCHZ7nnp/Peob67G7cTcAwCybMco4KrSfXYY+s9stD1rkFhw1F6MWtVCCQajUathMNhhkE0wYvokc0DbNsj+3GDgTbj9ARERE1Ee33nortm7dinXr1qGiogLLli3DunXr8I1vfAMOhwMWiwUGgwFOpxOjR4/Gt7/9bdxyyy1Yu3Yt3njjDRw9ehQmkwm33347/vGPf+CVV15BamoqHnroIRw5cgT79u2DWq0+43V6Ih62H2hX7XGgwFmAam81AoofskoDu86OKYlTYNf3X+n2ocan+FDmLgtNxTzmPgaf4uv1efQqPfJNI9tG7YyjkGPIgaz679hPtceBj+s+hivQAqsmCVqVFj7Fh0Z/A0yyGQtSFvDfaZAwkSMiIiLqI7fbjRUrVmDDhg2wWCxYuXIl7rzzTgBt64nWrVsXGmXbtWsXbr31VhQVFWHq1KlYs2YNZsyYAaBtM+wHH3wQ69evR1NTEy644AL85je/Ca1rO9N1eiJeEjkmCf0nKII40Xri1CblJShxl6A16O71eTSSBrnGU5UxjfkodZfC4amCXZceUeyk2utAtiEHF6QuiqmNzIcqJnJEREREQ1w8JHKKUPDhyQ9worWcScIAUISCam81SlzFp9bZFcMZ6Mu2FBISZDOStSlI0iQhSWOFRtW2NYYn2ApX0I0l9iVI1qb07wuIE4pQ0OhvgEdpmxZs1SQN2PuVa+SIiIiIKOoa/Q2o9lbDqkmKKL4hSRKsGiuqvdVo9DcM2yThbKgkFTL0GcjQZ2BeynwIIVDnqzu15UExil0lqPPV9uBMAs2BZjQHmnEcxwAAWpUWRrUJRrUBihD4oukLjDWPg01ng1alPfPphpDBnhbMRI6IiIiIos6jeBFQ/NBqOv/gr1Xp0OR3wqN0vnce9Y4kSbDpbLDpbPhK0mwAQJO/qcNUzGJUeap6dC6f4gtNgQWAfzr+EeqzyBbYdKlI1aae+tOGVF0qUrQ2aFSa/n9hURIxLVjTNi34RGs5GvwNAzItmIkcERHFNaEEoNSVh7WpUnIgqfgjjiie6FU6yCoNfIoPenXkhtI+xQtZJUOvis9tFeKBRWPBDOtMzLDOBAC4Aq5QZcwSVzFOtJ6AAqVX52wKNKEp0IRi19Gwdglto6w2bSpSdTbYOiR6KVpbWIGVWKcIBQXOArgCLWHTgvVqPeyqdFR720bqUnVp/TrNMn7uEBERUWf8Xvi3/iGsSXfpKkDHH3FE8cSqSYJdZ29bI6eKXCPX6G9EtiEHVs3glHZvN5hrnmKNSTZhUuJkTEqcDADwBr3Y17QX2+q2osHfAHfQ3aO97DojINDgb0CDvwFHXF+G9UmQkKRJ6jCSZwuN6KVoU6CW1Gf92vpTtKYF86ccEREREUWdSlJhSuIUNPgbUO11wKqxQqvSwad40ehvhEk2Y0rilEFNorgVQjidWofZyXNgls34pPYTnPTWwC/8AACDyoBkXTICIoha70k0BZr6fB0B0bbfnb8eX+JwWJ8KKiRpk5GqtZ02ZTMVSdqkqCR50ZoWzESOiIiIiGKCXZ+OBSkLQslTk98JWSUj25Az6MlTNNY8xYNqjwOFzYXQqbQYYx4LlSRBEQKtQTfMmoTQffEqXtR563DSdxK13pMd/qxFc5+qZbZRoKDOV9tWmKXlUFifCiqkaFNCa/BSOyR6SQM4khqtacFM5IiIKP5pjdGOgIj6iV2fjlRdWlSnM0ZrzVOs63hf0vUZYdMILcISdl90Kh0yDZnINGRGnMcT9KDWV/vfBK/D31sCLX2PDwpO+trOczq1pG5L8jqM4LX/adFYzurfMVrTgpnIERFRXJN0Jugvuz/aYQwNhWsBb92Zj9GlABNvGZx4aNhSSaqobjHArRA611/3Ra/WI9uQjWxDdkRfa7D1v4md9yRqfe1/1sIVdPU59qAIosZbgxpvDdAc3idLGti0KRHVNW26VFhkS8RrPV3HacEOTxUMagNUkgqKUNAabIVZkzAg04KZyBEREVEbbSJgm9H9ca4qwJQx8PEQRQm3QujcYNwXg9qAHEMOcgw5EX3uoLstqQtN1axtS/R8tWgNuvt8zYDww+F1wOF1RPRpJS1SdDbYTk3VbKuy2ZboJciJoSTPrk/HxISJ+KT2E1R4KhAUQaglNWzaVHwlYSL3kSMiIqKBJAE9KhTQtyp1RPGCWyF0Ltr3xag2Is+YhzxjXkSfK+DqdD1erfckPIqnz9f0CR+qPJWo8lRG9OlUOthOjdwZVQbU+E5CCIE8w4jQHnmtSisKmwth09q4jxwRERER0UCK1a0Qoi2W74tJNsEkmzDCOCKsXQiBlmBL2Chex9E871mMHnoVLyo8FajwVET0aSQt8k0jkGcYwX3kiIiIiIgGQyxuhRAL4vG+SJKEBDkBCXIC8k0jw/qEEGgONIfW4XVM9Oq8tfAJX5+v6xc+fNnyJZI0ydxHjoiIiIhosMTSVgixZCjdF0mSkKhJRKImESNNo8L6hBBoCjSFTdFsn7JZ66tD4NT+ed0JCD8SVAncR46IiOh0QglCNFWFtUmWDEiqwd8UloiGlljYCiEWDYf70l6F06qxYjTGhPUpQkGTvyk0kneitRxFzYfgFz60BlshTq0jTtGkwKpJ4j5yREREnfJ74Nu8NqxJd+kqQGeKUkBxTDYCQunBgWcuxU00lER7K4RYNZzvi0pSIUmbhCRtEsaYx0IRCj48+QFOtJYjTWuHV7SNvOlVbQVhuI8cERERDawRV0Y7AiKiuHP62sH2feS8QQ/3kSMiIiIiIopV3EeOiIioL9SaaEdARETDWLXHgcLmQujUOow2jYEkSRBCcB85IiKirkg6E/RXPhztMIiIaJhShIICZwFcgRak68L317MIy4DtIzd0SssQERERERENskZ/A6q91bBqksKSOOC/1S/b95HrT0zkiIiIiIiI+sijeBFQ/NCqtJ32a1U6BJRAv+8jx0SOiIiIiIioj/QqHWSVBj7F12n/QO0jx0SOiIiIiIioj6yaJNh1djT6GyCECOsTQqDR3wi7zt7v+8gxkSMiIiIiIuqj9n3kTLIZ1V4HPMFWKEKBJ9iKaq8DJtk8IPvISeL0tJGIiCiOCCUI0VIb1iaZbZBU6ihFRBR7nE4nLBYLmpqakJiYGO1wiIakak9bdcpqbzUCSgCySoZdZ8eUxCncR46IiCiC3wPfBy+ENekuXQXoTFEKiIiIhiO7Ph2pujQ0+hvgUbzQq3SwapL6fSSuHRM5IiIiIiKifqCSVEjWpgzOtQblKkRERERERNRvOCJHRMOeEALNzc3RDoP6SHhd8LrD9+bROZ2QdMEoRUT9LSEhIWKTXeqd9pIITqczypEQUWf68v8ci50Q0bDXXgSAiGITC3ScvRMnTiAnJyfaYRBRF/ry/xwTOSIa9uJ9RM7pdCInJwfl5eVx+WGX8UdXPMTPEbmzpygKKisru72X8fB+iAbel87xvnSuL/elL//PcWolEQ17kiQNiR9AiYmJcf06GH90xXv8dGYqlQrZ2dk9Pp7vh87xvnSO96VzA31fWOyEiIiIiIgozjCRIyIiIiIiijNM5IiI4pxOp8MjjzwCnU4X7VD6hPFHV7zHT/2L74fO8b50jvelc4N1X1jshIiIiIiIKM5wRI6IiIiIiCjOMJEjIiIiIiKKM0zkiIiIiIiI4gwTOSKiGOTxeHDTTTfBarUiIyMDzzzzTJfHvv3225g+fTrMZjOmTp2Kf/7zn6E+IQSeeuop5OfnIzExEYsWLUJhYWHcxN/RX//610HbFLo/4//b3/6GsWPHwmQyYfHixTh+/PhAh9+v759HH30U2dnZSEpKwrXXXouTJ08OePw0MN58801IkhT2tXTpUgDA3r17MXv2bBiNRsyaNQt79uwJe+769esxatQoGI1GXH311aitrY3GS+hXXq8XkydPxubNm0NtpaWluPDCC2EymTBx4kT8+9//DnvO+++/j8mTJ8NoNOKCCy5ASUlJWP/zzz+PrKwsJCQk4KabboLb7R6Ml9KvOrsvd9xxR8R754UXXgj1n+n9IYTAqlWrkJqaiuTkZNx7771QFGUwX9JZqaiowNKlS5GcnIysrCzcdddd8Hg8AGLg/SKIiCjm/PjHPxZTp04Ve/bsEX//+99FQkKC+Otf/xpx3P79+4VWqxW/+tWvxJEjR8QLL7wgNBqN2LdvnxBCiN/+9rfCZrOJjRs3isOHD4ubbrpJ5ObmCpfLFRfxt2toaBDp6elisH5s9Vf827ZtE7IsizVr1ohDhw6Jyy+/XMyZMydu4l+zZo3Izs4WmzdvFgUFBWL+/PniyiuvHPD4aWA88cQT4oorrhBVVVWhr4aGBtHS0iLS09PF3XffLQoLC8Xtt98u7Ha7aGlpEUIIsXPnTmEwGMQrr7wi9u/fL772ta+Jyy67LMqv5uy0traKq6++WgAQH330kRBCCEVRxNSpU8V1110nCgsLxc9//nNhNBrF8ePHhRBCHD9+XJhMJvH000+LAwcOiGuuuUZMmTJFKIoihBDib3/7m7BYLGLjxo1i165dYuLEieJHP/pRtF5in3R2X4QQ4sILLxSrV68Oe++0/xzp7v3x9NNPi5ycHPHJJ5+IDz/8UGRmZopf/vKXg/3S+kRRFDFnzhyxZMkSceDAAfHxxx+L0aNHi3vuuScm3i9M5IiIYkxLS4vQ6/VhP0T/53/+R3zta1+LOPa+++4Tl1xySVjb4sWLxQMPPCCEEGL27NniySefDPX5fD5hMpnEv//97wGJXYj+jb/d97//fTFv3rxBSeT6M/6rr75aLF++PNRXUlIi8vLyxMmTJwckdiH6N/4rr7xS3H333aG+f/7zn8JkMg1I3DTwrrvuOnH//fdHtL/00ksiPz8/9AFTURQxevRosW7dOiGEEDfccINYtmxZ6PiysjIhSZIoKSkZjLD73cGDB8W0adPE1KlTwxKWDz74QJhMplACK4QQixYtEo888ogQQoiHHnoo7PvI5XKJhISE0PPPO++80LFCCPHJJ58Ig8Ew4L846y9d3RchhMjKyhLvvfdep8/r7v2Rk5MTei8JIcSf/vQnkZeXNwCvoP8VFRUJAMLhcITaXn/9dZGZmRkT7xdOrSQiijH79++H3+/H3LlzQ23z58/Hzp07I6ajLFu2DE8++WTEOZqamgAATz/9NK677rpQuyRJEEKE+gdCf8YPAFu2bMHmzZvx4IMPDljMHfVn/Js3b8Y3vvGNUHt+fj6OHTsGm802QNH3b/wpKSl4++23UVFRgdbWVqxfvx4zZswYsNhpYBUWFmLs2LER7Tt27MD8+fNDU5clScK8efOwffv2UP+CBQtCx+fk5CA3Nxc7duwYnMD72ZYtW3D++eeHXl+7HTt2YObMmTCZTKG2+fPnd3kfjEYjZs6cie3btyMYDOKzzz4L658zZw58Ph/2798/wK+of3R1X5xOJyoqKjp97wBnfn9UVlaivLw8rH/+/Pk4fvw4qqqqBuaF9KP09HS8++67sNvtYe1NTU0x8X5hIkdEFGOqqqpgs9mg1WpDbXa7HR6PB3V1dWHHTpgwAdOmTQs9PnjwID744AMsWrQIQNsPlezs7FD/iy++iEAggPnz58dF/F6vFz/4wQ/w61//GgaDYcBiHoj4Gxsb0dDQgEAggIsvvhjp6em46qqrUFFRERfxA8DDDz8MWZaRnZ2NhIQEfPLJJ1i/fv2Axk8DQwiBw4cP47333sPYsWMxatQorFq1Cj6fD1VVVcjMzAw73m6348SJEwDQbX+8WbFiBZ577jkYjcaw9rO5D42NjfB4PGH9siwjJSUlbu5TV/elqKgIkiThZz/7GbKzszFt2jS88sorof4z3Zf2ZK1jf3tSFA/3xWq14uKLLw49VhQFL7zwAhYtWhQT7xcmckREMcbtdkOn04W1tT/2er1dPq+2thbf/OY3MW/ePFx11VUR/Tt37sTdd9+NlStXIj09vX+D7qA/4/+f//kfzJw5E4sXLx6weE/XX/G3tLQAAG6//XZcf/312LhxI7xeLy6//PIBXejfn/f/2LFjMBqN2LhxI7Zs2YLs7Gx873vfG7DYaeCUlZWF3htvvPEGnn76abz22mtYuXJll++Z9vdLd/1Dxdnch/YiFUPxPh06dAiSJGH8+PF455138P3vfx8/+MEP8OabbwLo/X3pyf9Hseree+/F559/jp/97Gcx8X6R+/IiiIho4Oj1+oj/yNsfn/6b0nbV1dW46KKLoCgK/va3v0GlCv893fbt27FkyRIsWbIEjz/++MAEfkp/xX/gwAH87ne/Q0FBwYDGe7r+il+W237Efv/738cNN9wAAHjttddgt9uxY8eOsKmPsRi/EAI33ngjfvnLX+Lyyy8HALzxxhvIy8vDzp07MXv27AGJnwZGXl4e6urqkJSUBEmSMH36dCiKguuvvx4LFy7s9D3T/n7p6j3V1fspXun1+ohR657cB6vVCr1eH3rc1fPj1Y033ogrrrgCycnJAICpU6fiyy+/xG9/+1tcffXVZ3x/dLwvp9+jeLsv9913H55//nn85S9/weTJk2Pi/cIROSKiGJOVlYXa2loEAoFQm8PhgMFggNVqjTi+oqICCxYsgNfrxebNm5GamhrWv3nzZlx00UW44IILsH79+ogkL1bj37BhA+rr6zFq1CiYzWYsWbIEAGA2m/Haa6/FfPw2mw0ajQbjx48PHZuSkoKUlBSUl5fHfPwnT55EeXl52NTLnJwc2Gy2QdlCgfpfcnJy2BYeEyZMgMfjQXp6OhwOR9ixDocDGRkZANreU2fqHyq6e51n6k9JSYFerw/rDwQCqKuri/v7JElSKIlrN2HChNA08TPdl6ysrNDjjn0A4uq+3HbbbXjmmWfw6quv4pvf/CaA2Hi/MJEjIoox06dPh0ajCSsksHXrVsyaNSsiCXO5XLjkkkugUqmwZcuWiPn4Bw4cwJVXXoklS5bgjTfegEajiZv4b7vtNhw6dAj79u3Dvn378OKLLwIA9u3bhyuvvDLm45dlGeecc07YwvXa2lrU1tZixIgRMR9/cnIydDpd2L6DtbW1qKurQ35+/oDFTwPjvffeQ0pKStg+Vfv27UNKSgrOO+88fPrppxBCAGhbT7dt2zbMmTMHQFsRhq1bt4aeV15ejvLy8lD/UDFnzhx8/vnnaG1tDbVt3bq1y/vgdruxd+9ezJkzByqVCrNmzQrr3759OzQaTdgvQ+LRww8/jAsvvDCsbd++faFfUp3p/ZGZmYnc3Nyw/q1btyI3NzduErnHHnsMa9aswZ///Gd861vfCrXHxPult2U4iYho4N1yyy1i0qRJYteuXeLNN98UiYmJYsOGDUIIIaqqqoTb7RZCCPHAAw8Ig8Egdu7cGba/T2NjoxBCiLlz54qJEyeKsrKysP7258d6/B199NFHg7aPXH/F/9e//lWYTCbxxhtviMLCQnH55ZeLmTNnhsq8x3r8t956q8jPzxdbtmwRBQUF4uKLLxZz584d8Pip/zmdTpGVlSW+/e1vi0OHDol33nlHZGZmiqeeeko0NTWJ1NRUcfvtt4uDBw+K22+/XaSnp4fKqn/66adCq9WKF198Uezfv18sXLhQXHHFFVF+Rf0DHcrsBwIBMXHiRHHttdeKAwcOiNWrVwuz2RzaF6y0tFTo9XqxevXq0L5gU6dODX0/rF+/XiQmJoo333xT7Nq1S0yaNEncdttt0XppZ6Xjfdm1a5eQZVn88pe/FEePHhW/+c1vhE6nE59++qkQovv3x+rVq0VmZqb46KOPxEcffSQyMzPFM888E42X1WuFhYVCrVaLn/70p2H/R1ZVVcXE+4WJHBFRDHK5XOLGG28UJpNJZGZmiueeey7UByC0J8+4ceMEgIivZcuWiaqqqk77Oj4/luM/3WAmcv0Z/+9+9zuRl5cnDAaDWLJkiSgvL4+b+FtbW8Xdd98tsrKyRHJysrj22mtFTU3NgMdPA+PAgQPiwgsvFGazWWRkZIhHH3009KFy586dYsaMGUKv14uvfOUr4vPPPw977rp160ROTo4wmUzi6quvFrW1tdF4Cf2uY8IihBBHjhwRCxYsEDqdTkyaNEn85z//CTv+nXfeEWPHjhUGg0EsWrQoYi+91atXi7S0NGGxWMT3vvc90draOhgvo9+dfl/eeustMXXqVKHX68X48eNDvxhqd6b3RyAQED/5yU+E1WoVNptN3HfffXHzy6DVq1d3+XNUiOi/XyQhTo2jExERERERUVzgGjkiIiIiIqI4w0SOiIiIiIgozjCRIyIiIiIiijNM5IiIiIiIiOIMEzkiIiIiIqI4w0SOiIiIiIgozjCRIyIiIiIiijNM5IiIiIiIiOIMEzkiIiIiIqI4w0SOiIiIiIgozjCRIyIiIiIiijNM5IiIiIiIiOIMEzkiIiIiIqI4w0SOiIiIiIgozjCRIyIiIiIiijNM5IiIiIiIiOIMEzkiIgrT0NCAu+++G/n5+TAajZgwYQKef/55KIoyaDFs3rwZkiT1+Pi//vWvqKmpAQA8+uijWLhw4QBFNjA8Hg++/vWvw2AwYOHChTh69CimT58OvV6Phx56KNrhERFRDJKjHQAREcWOuro6zJkzB5mZmXjppZeQn5+PXbt24bbbbkNxcTH+7//+L9ohRjh+/DiuueYalJaWAgDuuece3H777VGOqnfeffddvPvuu9i2bRsyMzPx1FNPAQAOHjyI5OTkKEdHRESxiIkcERGFrFq1CjqdDu+99x70ej0AhEbmrrrqKtx2220YO3ZslKMMJ4QIe2w2m6MUSd81NTXBbrfjnHPOCT2eNm0aRo0aFeXIiIgoVnFqJRERAQC8Xi/+/Oc/48c//nEoiWt3+eWX44MPPkBeXh6AtumXP/jBD2C322GxWHDDDTegoaEBQNu0yBEjRmDFihWwWCx46qmnsHz5cixfvhzTpk1DWloajhw5gsbGRtxwww1ITExEZmYmbrvtNrS2tnYa27Zt2zB//nwYjUaYTCZceumlqKqqAtCWaLb/+fLLL0dMrdy+fTvmz58Pk8mE/Px8rFmzJtS3fPly3HXXXbj22mthNBqRk5ODP/3pT13eo5qaGlx77bVITExEeno6HnjggVAieeLECVxzzTVITk6GzWbD7bffDq/XG3ruJ598gnPPPRcGgwFTpkzBhg0bAAAvv/wyli9fjrKyMkiShBEjRuDll1/GH//4R0iShGPHjvXkn4+Ihrhjx45BkiT8/e9/x6hRo6DX63H55Zejvr4+2qFRlDCRIyIiAEBxcTFaWlowa9asiD5JknD++edDp9MBAK6++mrs27cP//rXv/Cf//wHRUVFWL58eej448ePw+PxYM+ePfj2t78NAPjTn/6EJ554Am+//TbGjBmDm266CU1NTdi2bRveeustfPbZZ/jxj38cce2mpiZcdtllWLx4MQ4ePIh///vfOHr0KFavXg0A2LVrV+jPa6+9Nuy5RUVFuOCCC7BgwQJ8/vnnePTRR3H33XfjzTffDB3zwgsv4JxzzsGBAwfwzW9+E7fccguampo6vUdf//rXUVVVhS1btuCNN97AunXr8Otf/xo+nw8XXHABXC5XqO/tt9/GvffeCwBwOBy4/PLLsXz5chQUFOC+++7D8uXL8cknn+Daa6/F888/j+zsbFRVVWH//v245pprcM0116Cqqgo5OTk9/SckomHg5z//OdavX48tW7bgs88+wzPPPBPtkChKOLWSiIgAAI2NjQAAi8VyxuO++OILbNmyBYcPHw5Ns3z11VcxYcIEHD58OHTcfffdh9GjR4cez5o1C1dccQWAtqTxrbfeQn19feh6v//97zF9+nQ8++yzYddrbW3FQw89hLvuuguSJCE/Px/f/OY3Qwlcampq6E+DwRD23N///veYMWMGfv7znwMAxo0bh6KiIvziF7/A1VdfDQCYNm1aKOF6/PHH8atf/QoHDx7E3LlzI1739u3bUVJSEhoFXLNmDVpaWvDuu++ioqICO3fuRFJSEgDg17/+Na644gr87Gc/w69//WtceOGFoUR19OjR2Lt3L55//nls2LABFosFarUa6enpABB6He2PiYjaPfbYY/jKV74CALjuuuvw2WefRTkiihYmckREBABISUkBgNAUya4UFRXBarWGrZUbP348kpKSQn0AMGLEiLDndXxcVFQERVGQlZUVdoyiKDh69GhYW3p6OpYtW4bnnnsO+/btQ2FhIfbv34958+Z1+5qKioowe/bssLa5c+eGTa8cM2ZM6O+JiYkAAL/fH3Guw4cPIzk5OZTEAcBVV10FAHjqqacwduzYUBLXfp1AIICjR4+iqKgIGzduDFu/5/f7Y269IRHFvtP/z+rs/ysaHpjIERERAGDUqFGwWCzYs2dPp9Mr24udnL5+rl0wGEQwGAw9Pv24jo8DgQAsFgt2794dcZ6srCzs3Lkz9LiiogLnnnsuzjnnHFx00UW4+eab8fbbb2PHjh3dvqbOYj09Tq1WG3HM6QVUAECj0fT6Ou1/BgIBXH/99XjggQd6fE4ios6c/n9WZ/9f0fDANXJERAQAkGUZ3/rWt/DCCy/A5/OF9W3cuBH//Oc/kZmZiXHjxqGxsTFsGmVhYSGcTifGjRvXo2uNGzcOTU1NkCQJo0ePxujRo9Ha2oqVK1eGFQgBgDfffBPJycn417/+hTvuuAPnnXceSkpKQh9ezrTf3Lhx4yISvu3bt/c4zo7GjBmD+vp6lJeXh9r+93//F1//+tcxbtw4fPnll2FFB7Zv3w5ZljFq1CiMGzcOR44cCb3W0aNH4x//+Adee+21XsdBREQEMJEjIqIOHn30UTidTlx88cXYsmULiouL8dJLL2HZsmW44447MHHiRIwfPx5LlizBjTfeiM8++wy7du3CsmXLsGDBAkyePLlH15kwYQIuueSS0PqOzz//HMuXL0dLS0toama7lJQUlJWV4YMPPkBJSQmeeuopbNiwIZTwmUwmAMD+/fvR0tIS9twf/vCH2LdvHx544AF8+eWXeOWVV/DrX/8aP/rRj3p9byZNmoQLLrgAN910EwoKCrB582Y8+eSTuOiii3DRRRdh5MiRuOGGG1BQUICPPvoIt912G77zne/AarXihz/8IXbv3o2f/vSnOHLkCF5//XU88MADoSqgREREvcVEjoiIQtLT07Ft2zaMHDkS1113HSZPnoznnnsOjz/+eFhltD/+8Y8YOXIkFi1ahIsvvhiTJk3CW2+91atr/elPf0J+fj4WLVqECy+8EOPGjcOf//zniOOuueYaXH/99Vi6dCnOPfdcfPjhh3jmmWdQVFQEr9cLm82G66+/Htdccw1efPHFsOfm5ubiX//6F959911MmTIFTzzxBJ599ll897vf7dP9efXVV2EymTBnzhx85zvfwQ9+8AP88Ic/hFqtxj//+U8AwOzZs/Gtb30LV111FdauXQsAyMvLw8aNG7Fp0yZMnjwZP/3pT/HMM8/guuuu61McREREkuDEWiIiIiIiorjCETkiIiIiIqI4w0SOiIiIiIgozjCRIyIiIiIiijNM5IiIiIiIiOIMEzkiIiIiIqI4w0SOiIiIiIgozjCRIyIiIiIiijNM5IiIiIiIiOIMEzkiIiIiIqI4w0SOiIiIiIgozjCRIyIiIiIiijNM5IiIiIiIiOLM/wfY6jBaditclQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "make_plots(\n", + " df,\n", + " \"./plots/census_corr.pdf\",\n", + " intervals_xlabel=\"Correlation coeff\",\n", + " n_idx=-1,\n", + " true_theta=true_theta,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "2b114272-e293-4375-a5b8-0bf870857408", + "metadata": {}, + "source": [ + "### Power experiment\n", + "\n", + "For PPI and the classical approach, find the smallest value of ```n``` such that the method has power 80% against the null $H_0: \\theta^* < 0.5 \\cdot 10^{-5}$." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "c4fd41f6", + "metadata": {}, + "outputs": [], + "source": [ + "# Find n such that we reject H0: Corr coeff < 0.15 with probability 80% using a test at level alpha\n", + "num_experiments = 100\n", + "list_rand_idx = [\n", + " np.random.permutation(n_total) for i in range(num_experiments)\n", + "]\n", + "\n", + "\n", + "def _to_invert_ppi(n):\n", + " # print(f\"PPI: {n}\")\n", + " n = int(n)\n", + " nulls_rejected = 0\n", + " # Data setup\n", + " for i in range(num_experiments):\n", + " # print(f\"PPI: {n}, {i}\")\n", + " rand_idx = list_rand_idx[i]\n", + " _X, _X_unlabeled = X_total[rand_idx[:n]], X_total[rand_idx[n:]]\n", + " _Y, _Y_unlabeled = Y_total[rand_idx[:n]], Y_total[rand_idx[n:]]\n", + " _Yhat, _Yhat_unlabeled = (\n", + " Yhat_total[rand_idx[:n]],\n", + " Yhat_total[rand_idx[n:]],\n", + " )\n", + "\n", + " ppi_ci = ppboot(\n", + " pearson,\n", + " _Y,\n", + " _Yhat,\n", + " _Yhat_unlabeled,\n", + " _X,\n", + " _X_unlabeled,\n", + " alpha=alpha\n", + " )\n", + " if ppi_ci[0] > 0.15:\n", + " nulls_rejected += 1\n", + " return nulls_rejected / num_experiments - 0.8\n", + "\n", + "\n", + "def _to_invert_classical(n):\n", + " # print(f\"Classical: {n}\")\n", + " n = int(n)\n", + " nulls_rejected = 0\n", + " # Data setup\n", + " for i in range(num_experiments):\n", + " rand_idx = list_rand_idx[i]\n", + " _X, _X_unlabeled = X_total[rand_idx[:n]], X_total[rand_idx[n:]]\n", + " _Y, _Y_unlabeled = Y_total[rand_idx[:n]], Y_total[rand_idx[n:]]\n", + " _Yhat, _Yhat_unlabeled = (\n", + " Yhat_total[rand_idx[:n]],\n", + " Yhat_total[rand_idx[n:]],\n", + " )\n", + "\n", + " classical_ci = classical_bootstrap_ci(pearson, _X, _Y, alpha=alpha)\n", + " if classical_ci[0] > 0.15:\n", + " nulls_rejected += 1\n", + " return nulls_rejected / num_experiments - 0.8" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "8ca727f2", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The PPBoot test requires n=264 labeled data points to reject the null.\n", + "The classical test requires n=314 labeled data points to reject the null.\n" + ] + } + ], + "source": [ + "n_ppi = int(brentq(_to_invert_ppi, 50, 800, xtol=50))\n", + "n_classical = int(brentq(_to_invert_classical, 50, 1000, xtol=50))\n", + "print(\n", + " f\"The PPBoot test requires n={n_ppi} labeled data points to reject the null.\"\n", + ")\n", + "print(\n", + " f\"The classical test requires n={n_classical} labeled data points to reject the null.\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a3676a6-fcb2-44e5-a0ff-85935df5f0c0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/ppi_py/baselines.py b/ppi_py/baselines.py index f6c4e2e..36e93ca 100644 --- a/ppi_py/baselines.py +++ b/ppi_py/baselines.py @@ -5,7 +5,7 @@ from statsmodels.stats.weightstats import _zconfint_generic, _zstat_generic from sklearn.linear_model import LogisticRegression, LinearRegression from sklearn.isotonic import IsotonicRegression -from .utils import dataframe_decorator +from .utils import dataframe_decorator, bootstrap from .ppi import _ols, _wls import pdb @@ -309,3 +309,100 @@ def classical_logistic_ci(X, Y, alpha=0.1, alternative="two-sided"): return _zconfint_generic( pointest, np.sqrt(np.diag(cov_mat) / n), alpha, alternative ) + +""" + BOOTSTRAP CI + +""" + +def classical_bootstrap_ci( + estimator, + Y, + X=None, + n_resamples=1000, + alpha=0.1, + alternative="two-sided", + method="percentile", +): + """Classical bootstrap confidence interval for the estimator. + + Args: + estimator (callable): Estimator function. Takes in (X,Y) or (Y) and returns a point estimate. + Y (ndarray): Gold-standard labels. + X (ndarray, optional): Covariates corresponding to the gold-standard labels. Defaults to `None`. If `None`, the estimator is assumed to only take in `Y`. + n_resamples (int, optional): Number of bootstrap resamples. Defaults to `1000`. + alpha (float, optional): Error level; the confidence interval will target a coverage of 1 - alpha. Must be in (0, 1). Defaults to `0.1`. + alternative (str, optional): Alternative hypothesis, either 'two-sided', 'larger' or 'smaller'. Defaults to `'two-sided'`. + method (str, optional): Method to compute the confidence interval, either 'percentile' or 'basic'. Defaults to `'percentile'`. + + Returns: + float or ndarray: Lower and upper bounds of the bootstrap confidence interval for the estimator. + """ + + if X is None: + + pointest = estimator(Y) + + bootstrap_distribution = np.array( + bootstrap( + [Y], + estimator, + n_resamples=n_resamples + ) + ) + + else: + + pointest = estimator(X, Y) + + bootstrap_distribution = np.array( + bootstrap( + [X, Y], + estimator, + n_resamples=n_resamples, + paired=[[0, 1]] + ) + ) + + # Deal with the different types of alternative hypotheses + if alternative == "two-sided": + alpha_lower = alpha / 2 + alpha_upper = alpha / 2 + elif alternative == "larger": + alpha_lower = alpha + alpha_upper = 0 + elif alternative == "smaller": + alpha_lower = 0 + alpha_upper = alpha + + # Compute the lower and upper bounds depending on the method + if method == "percentile": + lower_bound = np.quantile( + bootstrap_distribution, alpha_lower, axis=0 + ) + upper_bound = np.quantile( + bootstrap_distribution, 1 - alpha_upper, axis=0 + ) + elif method == "basic": + lower_bound = 2 * pointest - np.quantile( + bootstrap_distribution, 1 - alpha_lower, axis=0 + ) + upper_bound = 2 * pointest - np.quantile( + bootstrap_distribution, alpha_upper, axis=0 + ) + else: + raise ValueError( + "Method must be either 'percentile' or 'basic'. The others are not implemented yet... want to contribute? ;)" + ) + + if alternative == "two-sided": + return lower_bound, upper_bound + elif alternative == "larger": + return -np.inf, upper_bound + elif alternative == "smaller": + return lower_bound, np.inf + else: + raise ValueError( + "Alternative must be either 'two-sided', 'larger' or 'smaller'." + ) + \ No newline at end of file diff --git a/ppi_py/ppi.py b/ppi_py/ppi.py index d499f6a..e9aa10c 100644 --- a/ppi_py/ppi.py +++ b/ppi_py/ppi.py @@ -1071,6 +1071,12 @@ def ppi_logistic_ci( ) +""" + PPBOOT + +""" + + def ppboot( estimator, Y,