{ "cells": [ { "cell_type": "markdown", "id": "2d3d3b21-222b-44ed-b58b-bede7e548ee1", "metadata": {}, "source": [ "# Ash & Hughes-Oliver (2022): Confidence Bands for Hit Enrichment Curves\n", "\n", "The following is a replication of the applied analysis of Ash & Hughes-Oliver (2022), which uses data from Empereur-Mot et al. (2016). The authors discuss estimation of the 'hit enrichment curve', which is a function commonly used to summarize effectiveness of a virtual drug screening campaign. In the example, the authors are focused on estimating the hit enrichment curve for the protein regulating gene peroxisome proliferator-activated receptor gamma (PPARg). This was done using three different docking methods: Surflex-dock, ICM, and Vina. To compare the hit enrichment curves from these different docking methods, the authors discuss the use of confidence bands. The original materials for the paper can be found [HERE](https://github.com/jrash/Enrichment-Inference-Supplemental-Materials). \n", "\n", "Here, we will show how `delicatessen` can be used to construct confidence bands in the context of hit enrichment curves.\n", "\n", "## Setup" ] }, { "cell_type": "code", "execution_count": 1, "id": "40f0959b-ac89-438b-abc1-bff678931110", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Versions\n", "NumPy: 2.3.5\n", "SciPy: 1.16.3\n", "pandas: 2.3.3\n", "Matplotlib: 3.10.8\n", "Delicatessen: 4.1\n" ] } ], "source": [ "import numpy as np\n", "import scipy as sp\n", "import pandas as pd\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "\n", "import delicatessen as deli\n", "from delicatessen import MEstimator\n", "from delicatessen.estimating_equations import ee_regression\n", "from delicatessen.utilities import inverse_logit\n", "\n", "print(\"Versions\")\n", "print(\"NumPy: \", np.__version__)\n", "print(\"SciPy: \", sp.__version__)\n", "print(\"pandas: \", pd.__version__)\n", "print(\"Matplotlib: \", mpl.__version__)\n", "print(\"Delicatessen:\", deli.__version__)" ] }, { "cell_type": "markdown", "id": "0e5df441-6814-4886-84c7-555bc5c63c20", "metadata": {}, "source": [ "Data is available from the GitHub page linked above." ] }, { "cell_type": "code", "execution_count": 2, "id": "30175ca8-2923-4e4e-b59f-6e7243ab347b", "metadata": {}, "outputs": [], "source": [ "d = pd.read_csv(\"data/pparg.csv\")\n", "\n", "# Extracting Relevant data for hit enrichment curve estimation\n", "active = np.asarray(d['surf_actives']) # Outcome variable for hit enrichment\n", "x_mxz = np.asarray(d['maxz_scores']) # Max Z-score docking\n", "x_srf = np.asarray(d['surf_scores']) # Surflex-docking\n", "x_icm = np.asarray(d['icm_scores']) # ICM docking" ] }, { "cell_type": "markdown", "id": "33af2ac6-8856-4f8a-a348-ba150152f4dd", "metadata": {}, "source": [ "## Empirical Distribution Function\n", "\n", "As in Ash & Hughes-Oliver, we are interested in estimating the hit enrichment curve. They define the hit enrichment curve at the point $s$ as the recall curve divided by the percent-captured curve. Here, $s$ is a set of scores the correspond to the belief that a ligand is active. For a given threshold $s$, the estimator of the hit enrichment curve is \n", "$$ \\hat{H}(s) = \\frac{\\hat{\\beta}(s)}{\\hat{\\rho}(s)} $$\n", "where \n", "$$ \\hat{\\beta}(s) = \\frac{\\sum_{i=1}^n A_i I(S_i \\ge s)}{\\sum_{i=1}^n A_i}, $$ \n", "$$ \\hat{\\rho}(s) = \\frac{1}{n} \\sum_{i=1}^n I(S_i \\ge s), $$ \n", "$A_i$ denotes a 'hit' or the activity value, and $S_i$ denotes the active ligand score. To plot the estimated hit enrichment curve, we will plot the ordered pairs $(\\rho(s), \\beta(s))$ for $s$.\n", "\n", "The estimating function for these parameters are\n", "$$ \\psi(O_i; \\beta(s)) = A_i \\left( I(S_i \\ge s) - \\beta(s) \\right) $$\n", "$$ \\psi(O_i; \\beta(s)) = I(S_i \\ge s) - \\rho(s) $$\n", "respectively. However, we are interested in plotting the complete curve (not only one point). To estimate this function, let $s_k \\in \\{s_1, ..., s_K\\}$ denote the set of unique observed values of $S_i$ where $A_i = 1$ in the data. Here, the stacked estimating functions expand to the number $K$.\n", "Note the the number of estimating functions (and parameters) is now a function of $n$ (via $K$). This seemingly subtle change has important implications for inference. Specifically, the theory from Stefanski & Boos 2002 is no longer adequate. Instead, we need an extension. See Kosorok 2008 for those details.\n", "\n", "We will start by estimating the cumulative distribution function for Max Z-score docking." ] }, { "cell_type": "code", "execution_count": 3, "id": "7bfd0658-22b2-41bd-bb4d-894fe50ef476", "metadata": {}, "outputs": [], "source": [ "unique_mxz = np.unique(x_mxz[active == 1])[1:] # All unique values for A=1\n", "n_unique_mxz = len(unique_mxz) # Number of unique values, K" ] }, { "cell_type": "code", "execution_count": 4, "id": "45d00a5d-21be-48e6-b71e-85afb9db5383", "metadata": {}, "outputs": [], "source": [ "def psi_edf(theta):\n", " # Separating parameters into subsets\n", " beta = theta[: n_unique_mxz]\n", " rho = theta[n_unique_mxz: n_unique_mxz*2]\n", "\n", " # A fast way to compute the indicator terms without a for-loop\n", " indicator_matrix = (x_mxz >= unique_mxz[:, None]).astype(int)\n", "\n", " # Recall fractions at all ranks\n", " ef_hedf = active*(indicator_matrix.T - beta).T\n", " # Fraction tests at all ranks\n", " ef_fedf = indicator_matrix.T - rho\n", "\n", " # Stacking the estimating functions\n", " return np.vstack([ef_hedf, ef_fedf.T])" ] }, { "cell_type": "code", "execution_count": 5, "id": "2c8c66b8-5b25-451a-91de-8e832be892d1", "metadata": {}, "outputs": [], "source": [ "init_vals = list(np.linspace(0.99, 0.01, n_unique_mxz)) + list(np.linspace(0.99, 0.01, n_unique_mxz))\n", "estr = MEstimator(psi_edf, init=init_vals)\n", "estr.estimate(solver='lm')" ] }, { "cell_type": "markdown", "id": "17e80bca-7225-43c0-a459-e9e55176e704", "metadata": {}, "source": [ "We now have all the pieces estimated that we need to create our figure. \n", "\n", "While we could present confidence intervals in our figure, these can be misleading when presenting functions. Confidence intervals claim to cover the parameter at the advertised rate. As such, coverage of a *function* will be below that level. Therefore, when presenting functions we should instead consider presenting confidence bands. Confidence bands modify the critical value so that a set of modified confidence intervals have coverage at the advertised rate. The most familiar method is the Bonferroni correction. However, the Bonferroni correction can be overly conservative (and is not appropriate to apply in this setting). Instead, we use the sup-t method as described in the referenced paper. See Zivich et al. (2025) for further details.\n", "\n", "The following function inside the `MEstimator` class allows for computation of the confidence bands." ] }, { "cell_type": "code", "execution_count": 6, "id": "b48a6254-69d7-4b35-a2e1-61c8b250d20e", "metadata": {}, "outputs": [], "source": [ "cbands = estr.confidence_bands(method='supt', seed=2101101)[:n_unique_mxz, :]" ] }, { "cell_type": "markdown", "id": "7aa84252-5aa9-40ba-9926-4afc6f0e71e2", "metadata": {}, "source": [ "To illustrate the difference between the two methods to compute confidence regions, we will also compute the confidence intervals here." ] }, { "cell_type": "code", "execution_count": 7, "id": "77b9567c-6a13-4456-ad2d-bd430f705903", "metadata": {}, "outputs": [], "source": [ "cints = estr.confidence_intervals()[:n_unique_mxz, :]" ] }, { "cell_type": "markdown", "id": "08328d84-2ed0-4a4c-a197-6adab2ecf1f7", "metadata": {}, "source": [ "Now we will recreate the plot in Figure 5 of Ash & Hughes-Oliver." ] }, { "cell_type": "code", "execution_count": 8, "id": "3db16cae-8ff8-46f2-8716-b83b1c8bcb23", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHWCAYAAAARl3+JAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAS29JREFUeJzt3Ql8VPW5//EnCWYhECBNBaSItuKKoAWhaF1LS1uvV+tyqVrABfty19JFaRGLG7VeFatWrnDVarVSvVfbf2utFVFcoFbEWr1F0KJYKyDFErYQQvJ/Pb944slkZnLOmbOfz/v1mhfMZGYyCWN88vs9z/dX1tbW1iYAAABItPKoXwAAAABKR1EHAACQAhR1AAAAKUBRBwAAkAIUdQAAAClAUQcAAJACFHUAAAApQFEHAACQAj0kY1pbW+Uf//iH9O7dW8rKyqJ+OQAAAKJnQWzatEl22203KS/3tuaWuaJOC7rBgwdH/TIAAAC6ePfdd+VTn/qUeJG5ok5X6KxvWl1dXdQvBwAAQBobG82ik1WneJG5os7actWCjqIOAADESSmtYQxKAAAApABFHQAAQApQ1AEAAKQARR0AAEAKUNQBAACkAEUdAABAClDUAQAApABFHQAAQApQ1AEAAKQARR0AAEAKUNQBAACkAEUdAABAClDUAQAApABFHQAAQAr0iPoFAAAAlGr7dpGdO0v/PlZUiFRVJfPfg6IOAAAkupjbtk3k1VdFtm4t/fl69hQZPlyksjJ5xR5FHQAASOQKXHPzx8Xc1q0i9fUiu+zi/XPt2CGyfr3IkiX5i70xY+Jd2FHUAQCA2BZ0f/xj8RU4q5jr06fr6ppb1dXtRVtra9fXoZ/Hj+3dIFHUAQCAUPrV3Gpqai+mtFgrtELmRzFnV+i5dFUw7ijqAABASatlQdLP63fhllYUdQAApIjfq2pOVsuCREGXkKJu0aJFcsMNN8jSpUvl/fffl0ceeUROOOGEoo95+umnZerUqfL666/L4MGDZfr06XLGGWeE9poBAMjaqlomVsuam0VaC1TDTfrxCt2clTiLtKjbsmWLjBgxQs466yw58cQTu73/qlWr5Nhjj5Vzzz1X7r//flmwYIFMmTJFBg4cKOPHjw/lNQMAEFca7RHEqlomCrrXXmtflsz78XKR1lqRw/cW6Rnf8ddIi7qvfOUr5uLUnDlzZM8995Qbb7zRXN9vv/3kueeek5tvvpmiDgAgWV+ls+I9Ul+EuV1p68725vaCTvNQ8mWitO0UaWyK/fhronrqFi9eLOPGjet0m67QXXrppQUfs337dnOxNDY2BvoaAQCIonfO6n3TeI9MFnSvFVlpc0If26tX/qJOiz6J//hrooq6NWvWSP/+/Tvdpte1UNu2bZvU1NR0ecysWbNk5syZIb5KAACi6Z2zVukyV9DpvnNTkZU2JwoVdAmSqKLOi2nTppnBCosWgDpgAQBAUqdZC02kZm7b1b5C11RkpS0jElXUDRgwQNauXdvpNr1eV1eXd5VOVVVVmQsAAGmaZs1M71yxXjl7L1yvbBd0iSvqxo4dK4899lin2/7whz+Y2wEASBtdoSs0zZqZgq67XjkfV+ja2tqkafu2rp9i+w5pam42H4+zSIu6zZs3y5tvvtkpsuSVV16R+vp62X333c3W6XvvvSf33nuv+bhGmdx2223yve99z8SgPPXUU/LLX/5Sfvvb30b4VQAAEEwosFXLaEGn55Jmjq7Qddcr52NBd/bV35BXVy4reJ91Z6+TWqmVuIq0qHvppZfk6KOP7rhu9b5NnjxZ7rnnHhNIvHr16o6Pa5yJFnDf+ta35JZbbpFPfepTMm/ePOJMAACpDQWOxfBDKXEhpTBTp9JetAXcStW0fVvRgi4JIi3qjjrqqKJLmVrY5XvMsmXJ/qYDALK34uY1FDjybVY/4kJKYW2vhuiJ256VmqqPe/WbNu2QzR82S8+ePSXOEtVTBwBA0lfcEtcH52QLNEghDEC0tbXJNlsvnRZ0NdUfF3Blzc3SUlkhZWVlEmcUdQAA+DC40J3QCjq/t0pD3AKNQpuDXrqkoKgDAMDhFmvsBxeC2iqNYAs0LE05vXQj9v6sVNu2XpOEog4AABdbrLEYXAh7qzQjGXBP3Pas9Kurj/02ayEUdQAAuNhijXVf3I6WVG+Vus2XcyK3l66spUWktbXznVoimPz1gKIOABDZVGgcxX6LtdjW68qVqd4qDbwnrmWHyIcbu/7D7ygXqaoVqaiQOKOoAwBEOhUaR7HeYi00EGEdmVVXl4mtUr/z5UZoL10PreRbRQ44QKTKthyrhf6Wis63xRBFHQBkkN+raqXksMVRrLdYiw1EZGyVrrt8OTd0OEKjSwwt3qpznmeHxB5FHQBkTFCraonNYUuiQgMRGRlocJovlzUUdQCQMaVkrRVDQRdi7lzKs+Oyli/nF4o6AMjYKl1iBwGyqFjuXMa3WoPKl2veUSatTV3/u0kCijoAyOC2ayIGAVA8dy5jW61h5Ms17yiT9z+skp6by7r00OmxrzEffqWoA4AsbruyVRoBL8d3ZXybtVj+XJd8ubISA4N37JDW5hbpWdUqo0e2SnV95w9rQRf3fwJW6gAgIxOu9m1XhhkSdHxXRrdZQ+2Z27FDZP16kfKeJo+uurbCrMwlDUUdAGRowpVt1wQe35XRbVan+XO+9NK1trY3mH5mP5GWmtjn0RVCUQcAGZpwZds1YhndRg0yf87ky/l1VmtVpUh5Mgs6RVEHACnDhGuMe+Mydq5qKcifc4+iDgBShAnXgGWsN44suGShqAOAFGHCNWAZ643z61zVUjjumduxo703zgt9bApQ1AFAgjHhGvIWa4YjRko5V7UUjnrmrOnVau9p2s0VNbK9JeZBdN2gqAOAhGLCNaIt1hhvo/rd/5aYvjZrevWAAzxNruo/+fvrekjPtl0SETJcCEUdACQUE64RbbHGdBuV/jdpL+iq3a8m6qZtzz4io0eL1NYmdxGWog4AEkpXF2J1hquXqdA4SugWa5D9b36dqxp31dWJ+ifvgqIOABK69frqqzEKEy5lKjSOYrzFGkX/m69ZcAgMRR0AJHjrtb4+Jkd+lTIVGkcx3WItpk3aktH/5jcdkkjJ9GqpKOoAIME81x1+b5UmdMsyLbSfbso1EyVz7FOv1dUi5QmdcPAJRR0AJLifLlZbpQnfskwy7adb8c5y8/e9h+ybif63LlOvNTUxWbaODkUdAGStny6ordIEblmm0bzp92Wv/62qUpqlUlqbvP83lQYUdQCQ1X46tkoTKzePzp4nVyYZK+isnLkPxWTMeZXkfDoLRR0AJJSnRTH9v1+CD5cHeXSFdmF79mzPmfMa76MFXdLbQSnqACAr7L109L8lVrE8ulTmyRU70zVn6rW6urTVuqSjqAOArLD30tH/lgq5eXSpy5NzcqYrU68dKOoAIK2Tr7mxJfbYEQYaEiHfWa6JOY81rDNdNcaktVIkJcMOpaCoA4A0Tr4Wii1h2zUxOMvVxZmuKTnIpFQUdQCQxsnXQrElbLum5izXVPbPeZ39YZXOoKgDgARyvHtKbEkq5DvLNXX9c16jTN5vH47omYJIklJR1AFA2vrpiC1JRG9cd2LVO1dsAjXoz+swyqS2NvmRJKWiqAOANPXTEVsSO4nvjXMygRokB9OtepeqjBd0iqIOANLUT0dsSeJ647oTee+ckwnUIGlBl/EzXZ2iqAOABKzQaUFnDbJ26qcjtiTxvXHdiU3vXHcTqIgcRR0AxLyg++Mf21foVKetV2JLQpf43jikGkUdACRgy1V3n7RnSAu6jp0oYktClZjeOL+HGroZVgiaWYwu8OUQZdIZRR0AJGDLVQu6gn3qxJaEIhG9cUENNUR0FJc9sqQQokw+RlEHAEnacrXb0RL2S0Pce+OCGmqIaFjBHllSqE7VbDomX9tR1AFAkrZc7UsYK1dy7FeIfXOJ6o1L2VCDFnTFVuvQjqIOAGKq6Jar1U9XV+fieAmkpm8OyIOiDgBi0j9nsfroirK2XnvwYzzsvrnIc+OKaWFLPsv4aQAAMeqfs3R7agRbr5H1zcUmNy7fkERjo0jfvpEMNSB6FHUAEKP+OUveProYb716yW+Lo0T1zRUakhg61PehhmKxIkEissQdijoAiIj+D8tRZEkhMdl6pQ8tZnbpEXqsSJCILHEuHj8RACDD265Ft1ozkN8WR7Hum4thrEiQiCxxjqIOACLedi261ZqB/LY4im3fXISIFYk/ijoAiJBuu4Zd0Pnd/5boPrQkyz0OLOLjvBA9ijoAiID2KXlWwkkS9L+lRKHjwCI6zgvxQFEHABH00736qsdeuhLjTILsf6MPLUSFjgMrr5BmqZRWJ1mHDjGBmhwUdQAQUT9dfb2HrVcf40z87n+jDy3648CCmlRlAjUZKOoAICIl1WQe4kx065X+t3QLalKVCdRkoKgDgJCPAnN0DJjP/XT00qVoIMK6rQgmVbOpPOoXcPvtt8see+wh1dXVMmbMGHnxxReL3n/27Nmyzz77SE1NjQwePFi+9a1vSVNJPyEBIJxMukWL2i/6Y063X8vLS+inc/ng3F46+t8SNhBhBRpaF72doQjEaaVu/vz5MnXqVJkzZ44p6LRgGz9+vLzxxhuy6667drn/Aw88IJdffrncddddcuihh8qKFSvkjDPOMFlCN910UyRfAwB4OQrMUzadT/102kvXr66eHLYkD0QonXJNS8Ahkr9Sp4XYOeecI2eeeabsv//+prjr2bOnKdryeeGFF+Swww6T0047zazufelLX5JTTz2129U9AIhDfIl1FJheSvp/sYN+OtM/17T140tOlhzBugkdiLBf8ryJiKrLtshW6pqbm2Xp0qUybdq0jtvKy8tl3Lhxsnjx4ryP0dW5n//856aIGz16tPztb3+Txx57TCZOnFjw82zfvt1cLI2NjT5/JQAQUHxJp9PUd4psdxZuR/9cNunbZMMGkYaG9sEGZE9kRd369etl586d0r9//0636/Xly5fnfYyu0OnjPv/5z5sfWi0tLXLuuefK97///YKfZ9asWTJz5kzfXz8ABB5fYv2f+rXXPp6ucJBPVyyLjl66mA08OHmMy8nX4cM/3uZHtiRq+vXpp5+W6667Tn7605+aHrw333xTLrnkErn66qvliiuuyPsYXQnUvj37Sp0OWABAmDy3wFl9dPoEetGCzsWT5WbRkSUXsxMgnHA5EEGbXXZFVtQ1NDRIRUWFrF27ttPten3AgAF5H6OFm261TpkyxVw/8MADZcuWLfLNb35TfvCDH5jt21xVVVXmAgBRRJh4Hs7P3XLVQs7DzzLOYk3AwEN3GIhA3Iu6yspKGTlypCxYsEBOOOEEc1tra6u5fuGFF+Z9zNatW7sUbloYKt2OBYA4RZjotqty3U/nYcu1UMAw4n0CBJCa7VfdFp08ebKMGjXKDD5opImuvOk0rJo0aZIMGjTI9MWp4447zkzMHnzwwR3br7p6p7dbxR0AxC3CxHV8icctVwYkssUs5tpa9DijFZEWdRMmTJAPPvhAZsyYIWvWrJGDDjpIHn/88Y7hidWrV3damZs+fboZw9c/33vvPfnkJz9pCrprr702wq8CAPKzIkw8c7nlSsBwOgYenCh0xitntGZbWVvG9i11UKJPnz6yceNGqdMATwDw4egvO11k0/jM3r0dFnVW/1zHEzeLvP56+/+hXRR1mkd3+DmjzN8JGI7pwMOwYb5MMuh7bNOmrme8ckZrtuuTRE2/AkBc++ZyOe6jy+2fszjoo9PfyXV1zkLAcPYGHjjjFXYUdQDg09Ffdo776HL75yzd9NHRPxcxBh4QQxR1AOBhu9VaWCu5b86n/jk7AoaBbKKoAwCP260lHf1lbb06PPqrGAKGQxyGiOBw1dwpV8WkK/KhqAMAj9utrqNKCvXSucihs7TJxzNuBAyHPAzh8oSHIKZcFZOuyEVRBwAu+Lbdau+lc3n0l/bTTblmog8vAp6GIUI84cE6zzV3ylUx6YpcFHUA4HDFxFc7Wtr/zB2QcNhPt+Kd5ebvew/Z15zninQPQzDlCie6HpYKAOjSv/Tqq+3br3mOmPZWIa5c2b5SV+ITzpt+nwllBwBW6gDAYT9dfb1Pu27W1qsGjLpcpcvtpyuTsmSftBBHEQxDAH6gqAMAhzzUX938BHb/Izi2/XSlnLQQRwEMQ+SbYu0OU65wg6IOABxm0pXEfhRYCTEmse2nK+WkhTjyeRii2BRrd5hyhVMUdQAQdCZdvqPAPMSYJKKfLkbDBXFSbIq1O0y5wimKOgAIOpMu31FgLmNM8impn87v/jf60BxhihVBoqgDgKAz6TweBZa4/rcQQ3kBdEVRBwBBszLp4iKo/rcQQ3mjxMAD4oqiDgCCZM+kK7GHznf0v7nGwAPijKIOAIJUYiZdsYw6hI+BB8QZRR0AhPLTtvQft75k1Gk/HUMNJWPgAXFEUQcAOXEmvubSlZBJ53tGnX1AgqEGIHUo6gAgTz6dr7l0AfTTecqosw9I1NRkYqgByBKKOgDIk0/nay6dD5l0vp75qgMSFHSeJlw5tgtxRlEHADk0Ss6XmsfHXLrYnvmawQlXju1CXFHUAYDXXjr7ea52PvbRldRPl3tqBAMSvky4cmwX4oqiDkDmeeqly3eeq12AuXSO+ukKnRrBgIRjTLgiaSjqAGSep166fOe52vnVR9fWZlbptm3f5q6frtCpERk59QHIIoo6ACilly7A81y1oDv76m/IqyuXJeLUCC/HZ8URwxBIKoo6APDSQxdA31wuXaHLLehG7P1Z9/l0MT8+K44YhkASUdQBgHjsoQvxPNcnbntWaqpqTEHnqJ8u5KGIUo7PiiOGIZBEFHUAIB576Hzqm3PSR6cFXU11z9ifGsFwARAdijoAcGJHS+A9dL700XFqBJBZFHUAMk93Vru9w8qVoWy3+tZH182pEX4PNTBcAESPog5Apmkx8uqr3eTTWVuvdXWBbLeW1EfnIWA4qKEGhguAaFHUAcg0K6Ouvt5BnEmPHoH1z1lc9dF5DBgOaqiB4QIgWhR1APBRq1zY8SW+9s95CBhmqAFIF4o6AJlV8LzXkOJL8vXPeeqjCzFgGEB8UdQByKSi572GGF+S2z9ncZRHBwA2FHUAMsnRea8BHwHmKYfOruWjmBUXdBGSSVUgnSjqAGSap/Ne43Cmqw5JNDaK9O3rOGDYPvXKpCqQPhR1ADJFV6l0lS5vL11IcnvpPOXQWUMSQ4c6rkrtU6+1tYFmKAOIAEUdgEz20ami2XQh0V66fnX13vvndnH/Y1xrQQo6IH0o6gBkso9Oi5qCvXQBa5O2Tr10ngo6D/10ANKNog5A5mhBVzB010wSBJNLZ/XTTblmYmlP4rCfLvcoMAYkgHSjqAOAfPl0AZ3zqv10K95Zbv6+95B93ffSOeynK3QUGAMSQHpR1AFAvny6gHPp1Lzp95WWRVekn67QUWAc5QWkF0UdAOSyhw77fL6rPZuuTIIPF+YoMCA7KOoAIMBeOl8y6QDAgXJx6aqrrpKtVh6AzbZt28zHACAudDDAOgZMLwWz6axeutdfb79Tuesfja7Pd/WUTecQp0YA2VTWpr9GulBRUSHvv/++7Lrrrp1u/+c//2lu26mZATHW2Ngoffr0kY0bN0pdXV3ULwdASJl0Fr0+cGDOfEHTNpGXl318LJiPW6/bmrbK4eeM6nK+a0lnu1rV6mcPFqmuKXpqxJgxZNIBSeBHfeJ6+1VrwHw/iP785z9LfX29pxcBAEFn0lmKZtP53EuXL5PO9fmuLnFqBJBdjou6fv36mWJOL3vvvXenwk5X5zZv3iznnntuUK8TAPzPpAuYL5l0HoOHOTUCyB7HRd3s2bPND6izzjpLZs6caZYILZWVlbLHHnvI2LFjg3qdAODoTFdLlGe7+ppJ5zF4GED2OC7qJk+ebP7cc8895dBDD5VdAs5vAgA/+ueiPtvVt0y6AsHDnBoBwHNP3ZFHHimtra2yYsUKWbdunfm73RFHHOH2KQEg/P65AFm5dIFm0u3Sg1MjAJRW1C1ZskROO+00eeedd8wPLjv9LTTu068A0stz/9yO7nvU4phLx6kRAEoq6nQYYtSoUfLb3/5WBg4c6M92AgBE1T+n+5crV/p21mu+XDpfM+nyDElwagQAT0XdypUr5eGHH5a99tqL7yCA5PfPWee9ai6Uz73CVi5dSZl0xYYkOne/AMg410XdmDFj5M0336SoA5Cu/rkePTyf52pn76PzPZcuZ0hixyb/nhpABou6iy66SL797W/LmjVr5MADD+wyBTt8+HA/Xx8ABJc/5/G818jPc/1oSGLDBpGGBj3pJ5qXASBeXB9weNJJJ8lf//pXk1d3yCGHyEEHHSQHH3xwx59u3X777Sbjrrq62qwCvvjii0Xv/69//UsuuOAC089XVVVlgpAfe+wx158XQMaVcN5rofNcA+uj021X3XfWP3OGJPT3aPuKJYDscr1St2rVKt8++fz582Xq1KkyZ84cU9BpwPH48ePljTfe6HK2rGpubpYvfvGL5mPa1zdo0CAzhdtX+0sAZKYW84XVS6e7DTog4bGfzn6eq52vfXTr13+8LKl/2vrpoohsAZCSom7IkCG+ffKbbrpJzjnnHDnzzDPNdS3udKr2rrvukssvv7zL/fX2DRs2yAsvvNCx7aurfACyQRerXn3V51DhEs97Dfw8V6uP7oADRKoq2ws6reRicGIGgIRvv6r77rtPDjvsMNltt93MSpnSVbZf/epXjp9DV92WLl0q48aN+/jFlJeb64sXL877mF//+tfmKDLdfu3fv78MGzZMrrvuOrLxgIwNSdTXl7BCpUt9Tdtc99JpH922pq3tlzwDEoHTgq66puMLt+3EAoC3lbo77rhDZsyYIZdeeqlce+21HQWVboFqYXf88cc7ep7169ebx2pxZqfXly9vPysx19/+9jd56qmn5PTTTzd9dDqFe/7558uOHTvkyiuvzPuY7du3m4ulUeMAACSa54U1q4/OCrZzmE0X+WBEDoYkAPiyUnfrrbfK3Llz5Qc/+IFU2EauNJD4L3/5iwRJjyTTfro777xTRo4cKRMmTDCvQ7dtC5k1a5b06dOn4zJ48OBAXyOAYOjvZiWFDOf20emUgY6OOqgQCw1G+DoM4WA4wsKQBADfBiXyTbnqJOqWLVscP09DQ4MpCteuXdvpdr0+YMCAvI/RiVftpbMXk/vtt5+JV9Ht3Mo8+zHTpk0zwxj2lToKOyC5gcO+9NNpIedxZNQ+GOHbMITT4YgcDEkAKGmlbs8995RXXnmly+2PP/64KbCc0gJMV9sWLFjQaSVOr2vfXD7ax6dbrno/y4oVK0yxl6+gs4rNurq6ThcAyQ0cHjgw/GKmTdq6DEboJbBjEu3DEZ89WGTYMCo4AP6v1Omqlw4qNDU1mT4TzZX7xS9+YbY5582b5/q5Jk+ebLZuR48ebXrydLXPmoadNGmSiS3R51bnnXee3HbbbXLJJZeYEGQ9skwHJS6++GK3XwaABNLFtdALurY2mXLNRImENRwBAEEUdVOmTJGamhqZPn26bN26VU477TQzBXvLLbfI17/+dVfPpT1xH3zwgRm80C1UDTDWFT9reGL16tVmItai26a///3v5Vvf+pY5uUILPi3wLrvsMrdfBoAsZtN5oP10K95pH97ae8i+wfTQ6XarbQeC0VYAXpS16a+hHmlRt3nz5rxBwXGlPXU6MLFx40a2YoEE9dNpi1nJW68aZfLysvYhCYc9dRphcvg5o8zfF839k/SsrpVA++cser3AtqvOemzaJHLEEe1fCoDk86M+cb1SZ9ezZ09zAYBYZ9O5pL/r6gqdsmfSlUlZ8OHCFitkGACCKur++c9/mu3ShQsXyrp16zoNLSg98QEAYpFNp/u2GmNicRA4HFkmHf1zAMIu6iZOnGgmUM8++2zT+xbY9BeATG61fpRn3sFzNl1u0LD9CYsEDkeSSecSp0kA8KWoe/bZZ+W5556TESNGuH0oADjKosvlKZvOHjRsX+bTgs7hsl/gmXRanbms0DhNAoBvRd2+++4r27ZFcO4hgMxk0eXOMGhB57m9rISgYSuTLhD2AYkC4cL5cJoEAN+Kup/+9Kdy+eWXm766YcOGmRMe7Aj3BVAKrb9yB0GD6p+LlH1AoqbGddXKDAWAkou6vn37mrHbY445pktzsW5N7MxtiAGAoHnsn4sFHZBwUaHRTwfAt6Lu9NNPN6tzDzzwAIMSAOLBh/650LW0uH4I/XQAfC3qXnvtNVm2bJnss88+bh8KAMEqoX/O2nGw59IFRpfbGht168NxL52inw6Ar0WdntP67rvvUtQBiE4A/XOh5tNZ/XRDh3pqjqOfDoAvRd1FF11kzlv97ne/KwceeGCXQQk9kxUAktY/l5tPF0ou3S4lHeoDAJ24/okyYcIE8+dZZ53VcZsOSDAoAaCUsGHHIcMh9M9pPl2/unp/c+l0y9U6gcdDNp0+VL9fAOBbUbdq1Sq3DwEAR2HDrkKGS+yfK9ZLp/l0vhd0ViadxWE2nRZ077+vZ223X9c/K5y34QHIENdF3ZAhQ4J5JQAk62HDjkKGd7ifGo28l86eSacRJkoLOgfNcdZwxOjR7U+hBZ1PtSyAlPHU0LFy5UpZuHChrFu3Tlqt7YSPaCgxAAQSNqzLVitX+po/F2ovnRZ01d6eW79H1modAPhS1M2dO1fOO+88aWhokAEDBnTaotC/U9QBCIzVT1dXF0j+XCC9dAAQ16LummuukWuvvVYuu+yyYF4RgEzQRTfPegQzNep7L12JQcMMRwBwo9zVvUXkww8/lFNOOcXtwwCg05DEq6+299SVl7uocuJ+nmuxoGGHgxH24YhNm9r/znAEgECKOi3onnjiCbcPA4AuQxL19Q6DdK1sutdfb99+dVwJxoCHoGH7cMQRR4iMGcNwBIDuOdrD+MlPftLx97322kuuuOIKWbJkSd7w4YsvvtjJUwKA87Y4ezZdnM9z9TlomOEIAG44+ilz8803d7req1cveeaZZ8zFTntRKOoAdLf16jhoOFdu4LAP2qRNAmGFDTsIGrb65yyEDAMIrKgjcBiA34HDroKGfc6ms2fUTblmov9PnBs2XKSfLjdc2EIfHQC3OHgQQCSBw46ChgPKprNn1K14Z7n5+95D9vUvny43bLhI0HBuuLCFkGEAbrnuNj7ppJPk+uuv73L7j3/8Y6ZiATgOHHY4MxB4Np1l3vT7/I8zscKGHXyxVv+cdeHUCACBF3WLFi2Sr371q11u/8pXvmI+BgCB8DmbLve81zIp87bNqnvKuRcHfXQA4DfXPyU3b94slXl+69Qp2EbNYgKAj2h9o1uuFs8DEj7z5bzX3L65XA5y6Uz03nbvLwEASirqNMZk/vz5XY4De/DBB2X//fd3+3QAMjAUYVd0QMKMge7MeSL/A4d9Oe81t28uV5E+utwBCYYiAERS1GlG3YknnihvvfWWHHPMMea2BQsWyC9+8Qt56KGHfHlRANI1FGHvDys4IGEFDOdbzgtgSMK3816tvjmX7AMStbX00AGIoKg77rjj5NFHH5XrrrtOHn74YampqZHhw4fLk08+KUceeaQPLwlAmmhBV2iHsmDAcO5AhE+Bw7rtqqt09l46V+e9Wtlz1t99oN8bhiIA+MFT5/Gxxx5rLgBQqIeupIDhAKqckvvo8vXQuTjPNd/TAYCfyKkDEFgPnauA4YDl9tG57qXL10PXTd9cIbrTvGGDSENDex4dAIRW1NXX18uKFSukoaFB+vXrV3SrYoP+pAKQOfl66BwHDIdM++h021ULOte9dB576PL10w0fztYrgAjOfu3du7f5++zZs3389AAy20MXIS3oaqpzzuWKQBwLXgApL+omT55s/mxpaTG/1Y4fP1769+8f9GsDgPhoCeb8WQCI5ESJHj16yLnnnitNcUkQBZBs2lzWtK39EkAenW90qkHD1UsYjDBfalP7hcBhALEYlBg9erQsW7ZMhgwZEsgLApAR+XLpAsyjK4k1JDF0qOfBCCto2ELgMIDIi7rzzz9fvv3tb8vf//53GTlypNRqaqaNZtYBgKdcOp/y6OyZdBZ7Np1nu3gLDLAHDVv9hjr1Sj4dAD+5/gn19a9/3fx58cUXd9ymfXb6A1T/3Gk/6BEAIsil8+1s1wCChu2rdQAQaVG3atUqX18AgOTTHjHXrbY7WkLNpHOVTUfQMIAsFHX00gEoFDrsOGxYm8xWrgylh87KpLM4yqYjaBhAAnlqEFm5cqUsXLhQ1q1bJ63W9sRHZsyY4ddrA5Cw0GHHYcNWP11dna9nulpyz3b1nElH0DCANBd1c+fOlfPOO8+cLjFgwIBOv/Hq3ynqgGzStjjXg6E9esSjfy4EBA0DCJrrn6jXXHONXHvttXLZZZcF84oAIMz+OQDIalH34YcfyimnnBLMqwGQ/gEJ7acLKGjYU/9cQPTL1O4UgoYBxPJECaUF3RNPPBHMqwGQuAGJF19s76krL3cROPz66+3VoKMHdd1u3da09eNLnv456xJlQadhw5s2tf+doGEAsVyp22uvveSKK66QJUuWyIEHHii75DQ52/PrAKRXSQMS+nPDQ9BwKP1zGmdSYi5dbtgwQcMAYlnU3XnnndKrVy955plnzMVOfyumqAOyxdOAhP0UiTj1z9nz6Uo459VC2DCAMBE+DCDgxjLbKTM+9tIF0j9nz6erqem2WrX65nLRRwcgCqXnCQDIHEcDElb/XO4dSwgcbpM2f/LnnOTTOSjotG+u0LFf9NEBiG1Rt//++8tzzz0n9fX15vr5558vV111lcmrUxpEvMcee8hWbbIBkFqOT5Cw98/Zt1o99NJZ/XRTrpkocemly+2by0UfHYCwOR49W758ubS0fHxW489//nNpbGzsnOjuOtsAQJIHJAYOdNBPpwWcNt5ZF48nSGg/3Yp3lpu/7z1kX3/z56xeOv3TZS+d1TeXe9EvFQASsf2qRVyuqOIDAMRkQMLeQxdQFp2aN/0+f3/euOylA4A4oqcOgD/y9dCV0D/XJZvOlkdXJmWR9NLZByMYhgCQ2KJOfyvO/c2YlTkgO7SI0a3Xgl0W+XroPPbPxfFs13yDEQxDAEhkUac/WL/whS9Ij48O4N62bZscd9xxUvnRb7X2fjsA6R2OUEUHJKweOp/kZtNFdZ5rvsEIhiEAJLKou/LKKztdP/7447vc56STTvLnVQGI7XCE1muOT5DwmWbT9aurj3SXgEBhAKkr6gBkjxZ0+eI7wqLZdFEVdCWeHAYAgWNQAkDe3jm7KNOK7IHDUfbTbdggorGcuuUKAInOqQvS7bffboKLq6urZcyYMfLiiy86etyDDz5ofms/4YQTAn+NQJZ65xYt6nzR/yR1+7U85J8YgQcOW7rpCbb66YYPJ38OQHxFXtTNnz9fpk6darZ3X375ZRkxYoSMHz/enFBRzNtvvy3f+c535PDDDw/ttQJZ6p3r3bvzxVHQsM8CDRy276tqkLqD0GHi6wDEWeRF3U033STnnHOOnHnmmeYosjlz5kjPnj3lrrvuKviYnTt3yumnny4zZ86UT3/606G+XiBLvXP2S9GCRvcnAwwbDiRwODd4eOjQTl+kfkm67awXMukApLKou/fee2V7np9wzc3N5mNu6GOWLl0q48aN+/gFlZeb64sXLy74OD1zdtddd5Wzzz7b5asHEFjo8Ouvt1dAHvdoTcBw09bOlzAChy279OiSSbdpU/tFr5NJByB1gxK6ovblL3/ZFFV2mzZtMh+bNGmS4+dav369WXXr379/p9v1up41m89zzz0n//3f/y2vvPKKo8+hBai9CLWfVwvAYbBwMfbQYY9hw5EFDOvWa56xVjLpAGSiqNMfvvm2QP7+979Ln4JppP7QwnHixIkyd+5cadAxNAdmzZpltmkB+BAsXIz9JIkSA4ZzBRI4rMXc+vUf7y/n6acjkw5AKou6gw8+uOOoMPvJEkpX21atWmVW8NzQwqyiokLWrl3b6Xa9PmDAgC73f+utt8yAhJ5kYWn96CBGfT1vvPGGfOYzn+n0mGnTpplBDPtK3eDBg129TiDN4hIsbA8Y1jw6Oy3ofO+ns3rpDjhApKaGKQgA2SnqrNgQ3fbU6dRetkO69agwjSRxe6KEPm7kyJGyYMGCjufXIk2vX3jhhV3uv++++8pf/vKXTrdNnz7drODdcssteYu1qqoqcwEQQLBwAAMSWtDVVNsOWA1aVSUFHYBsniihxduECRNMppwfdBVt8uTJMmrUKBk9erTMnj1btmzZYvrzlPboDRo0yGyj6uccNmxYp8f37dvX/Jl7O4CQBiSsEVHbL3qxptuuukrHEREAst5TpwWYn7RA/OCDD2TGjBmyZs0aOeigg+Txxx/vGJ5YvXq1mYgFEEw/nefTInwYkAidvY9OOcimA4BUFXX19fWyYsUK0wPXr1+/or0tG/QsHZd0qzXfdqt6+umniz72nnvucf35AHQekPA8HFHigETo7H10uu2qBR2JwgCyVNTdfPPN0lsj5T/6e1QHagMIZkDC9XCED710OkmvU6/KnkcXCi3oqmsKf2ldozgBIB1FnX3L9Ywzzgjy9QCIYEDCdUFXYi9dZLl03bBChzVomLBhAKntqXMa2ltXV1fK6wEQdz700hXKpQskj84Fe+hwbW17wQsAqSvqdMq02LarFUqsmXUAUnp6RAC9dPZcukDy6DzQtjsKOgCpLeoWLlzYqYD76le/KvPmzTNxIwAydnpEidqkLbpcugLopQOQmaLuyCOP7HRdT4L43Oc+J5/+9KeDeF0A4nZ6hFY9uvXqw4DElGsmSqiZdPbredBLByCTOXUAMnh6hH04QpUQNqz9dCveWW7+vveQfYProcvNpLPkyaajlw5AGlDUAXA3HOFj2PC86fcF10OXm0lnKZJNRy8dgMwWdXFoaAZQfCDCzpfhiBImCHTr1Z5JVyZlkWbSAUAmi7oTTzyx0/WmpiY599xzpVbn/m3+93//179XB8CXgQi7bocjrN65Tk9YWh9dKNl0DvvnAECyXtT1yfm/wDe+8Y0gXg+AAAYi7IoOR+T2ztmV0EeXL5vO10w6F/1zACBZL+ruvvvuYF8JgOgHInJ75+x86qOzsun61dX718LhoX8OANKGQQkAvvfO5Z7tquy9dJpNF0hPrsf+OXZqAaQBRR2Qwn66kgciMtQ/pzvOGzaINDRo/qY/Lw8AokBRB6R0QMLVaRE+BQt3d7arL710PvfPWRl1w4dzNBiAZKOoA1I6IOH4tAgfg4W7O9vVl/NdA+qfo/UOQNJR1AEppO1wjosUn4KF7T10uf1zgZzt6kP+HOe9AkgTijoAJQ9HBN5DZ9fS4svTcN4rgLShqANSdHKE6wEJs1TVHFgPna9ZdFY/XWOjSN++JefPcd4rgLShqANSdnKE4wEJey+dj3109h66kvvnCvXTDR3qWxMc570CSAuKOiBlJ0c4HpCw99L5GCwcWA+d3S7ef3SZQd/W9oIYANKEog7I4skRdvlOj4gj3XotMSXY3ken9E+y6QCkBUUdgPizZ9OVcJ6rvY9On0YLuhIPzgCA2KCoA7LG56DhUNiz6WpqSu6n06eyVusAIC0o6oAETrlaPE27+hw0rHEm9ly6QGk2HSnBAJAXRR2Q0ClXi6vjwHwKGg49n86nbDoASDOKOiChU64Wx9OuPgUNF8un8z2XzudsOgBIM4o6ICtTrgHTfLp+dfX+5tIFlE0HAGlEUQekfSDCzsfhiNxeOs2n872gKzGbzsqks5BNByDNKOqABAxHlDwQYefTcEQgvXS61WqvwqzbfMiks5BNByCtKOqAhAxHlDQQYefD6RGB9NLZs+hyecimy82ks5BNByCtKOqAhAxHRDkQEUovnT2LTqNL7LSg89hPRyYdgKygqAPSNhyh+44BBgsH3kunBV21+1U/+ucAZB1FHZAm9l46H3rnIsulc4n+OQCgqANixX5qhOvhiNxeOh965yLJpfOA/jkAoKgDYn1qhKvhCLt8AxJJyaUrAf1zALKM7VcgxqdGeBqOCIBuu+oqne+9dPYIExfRJfTPAUBXFHVAGgYjrJDhAAYkAs2ky40wcRBdQv8cAORHUQckWb6QYZ8HJHL76HzrpcsXYeIguoT+OQDIj6IOSLJ8IcMBDEjY++h021ULOt966TxGmNA/BwCdUdQBEU235vI07epzyLDVO2eX20dXU51z7pZXLS2eHubx1DAASD2KOiDC6dZcnqddk5ZBp5VZY6NI376ujv/S3eYNG0QaGtqP+wIAfIyiDohwujVX0WlX+0CExcfBiHy9c3a+ZtJZ/XRDh7oa77X66YYPD/z0MwBIHIo6IO7TrYUGIiwBnBxh9c7Z+dpHZ9nF24+gOMS8AEDcUNQBSR2IsJQ4GFEog8633rl8W68eGuPMkbbbA3lFAJAKFHVAnI/9CmggIrKzXO3ZdA4y6fJl0+mFfjoA6IqiDkjCsV8BCSyDzkk2XU2N431UezZdbS39dACQD0UdEBBfjv2yhiMCOCkilAy6Ytl0HhrjtB5kQAIA8qOoA+I4GJFvOMLngQjdeg28j85+tqt13ePTAACKo6gDkjIc4eNJEaH00uU721W56KVTZNMBgDMUdUDc+Twcka+XLpA+unxnuzo83zX3acimA4DuUdQBPh/3ZfE07WoPGA6hj87qpetXV+9PH519u9XaM/V4tmsusukAoDiKOiCA4748TbvmCxgOIFg4l/bS+VbQ5W63ethqtbfgKbLpAMAZijoggOO+LK6mXfMFDPvYR1coaDjQ7VYXW632LLpcZNMBQPco6oCwp1oj6KELNWjY43arPYsu9/uqYcNEmQBAcRR1QEaEHjTskRZ0+VbrAADFUdQBYbIPQuQKaTAi9KBhAEAoKOqAMKZaCw1C5HtyHwYjrN45u8CDhgEAkYpFUXf77bfLDTfcIGvWrJERI0bIrbfeKqO1sSaPuXPnyr333iuv6f8cRWTkyJFy3XXXFbw/EPlUa7FBiFw+DEaE1jsHAIiV8qhfwPz582Xq1Kly5ZVXyssvv2yKuvHjx8u6devy3v/pp5+WU089VRYuXCiLFy+WwYMHy5e+9CV57733Qn/tSO9Ua+/exS8DB5aQm2YNQuS7+DDpmq93LvA+Oo0z8XCWly5eap2rF6JLAKA0ZW36a32ExowZI4cccojcdttt5npra6sp1C666CK5/PLLu338zp07pV+/fubxkyZN6vb+jY2N0qdPH9m4caPU1dX58jUgHbSgW7SovWgLZKq1aZvIy8vapwACHOXc1rRVDj9nVKfeOTvf++js+XR6GTbMUcWbL8JE/z5mDJOuALKn0Yf6JNLt1+bmZlm6dKlMmzat47by8nIZN26cWYVzYuvWrbJjxw6pr68P8JUCHoV8QoT+jhZ675w9n66mxvESZr4IE6JLAMC7SIu69evXm5W2/v37d7pdry9fvtzRc1x22WWy2267mUIwn+3bt5uLvRIGQhHyCRGR99JpPp2HPWkiTAAgRYMSXv3oRz+SBx980PTZVRfYL5s1a5bMnDkz9NeG5E25ep5qjfCEiGK9dKFk0HnspbMeCgBISVHX0NAgFRUVsnbt2k636/UBAwYUfex//ud/mqLuySeflOHDhxe8n27t6iCGfaVOe/aQXcWmXD1NtUZ0QkQx2kvXr64+2Ay63F46l2e8btigPwPat1wBAAmffq2srDSRJAsWLOi4TQcl9PrYsWMLPu7HP/6xXH311fL444/LqFHtDeGFVFVVmYZD+wXZVmzKtaSp1giYHrqmre2XnF66wEOF7b10Docjcvvp9Pcxjv8CgJRsv+oq2uTJk01xpllzs2fPli1btsiZZ55pPq4TrYMGDTLbqOr666+XGTNmyAMPPCB77LGHybZTvXr1MhcglLNbi50MYQl4MCLSHjr7tqvHXjqVpAIaAOIu8qJuwoQJ8sEHH5hCTQu0gw46yKzAWcMTq1evNhOxljvuuMNMzZ588smdnkdz7n74wx+G/vqRQU5OhghhMKJQHl3gvXQlbLta3z4y6QAghUWduvDCC80lHx2CsHv77bdDelVACSdDWAIcjLCz59EFfp6rxwiT3Gw6vdBPBwApK+qAIM9qzeXblGuIAxCxPMvVw7arPZuutpZ+OgDwE0UdMnFWayhTrgGJPH/Ovu2qVZmLLBLTetj68XVr21UX+hiQAAB/UdQhFVOsbgsELeiS0qQfyVmuxfrolINeunzHgCm2XQEgGBR1yPYUawynWosJ5SzX7vrodNtVC7puquJ8x4ApjgIDgGBQ1AExm2rN7aMLvXfO2mbNvU1pQVftblWQY8AAIBwUdUAMp1oj66PL3Wa1c7jlatWDxJYAQLgo6pC46dXAzmqN8FgvJ310ofTO5W6z2nWz5Zqvh47+OQAID0UdEjm9msQp1lL76ELpnbN42GbN10NH/xwAhIeiDomcXi15ijXfQESEAxC5WXSR9dG5iCsphB46AIgGRR3SNb1a6kBECAMQsTzH1WVcCQAgfijqkD3FBiJCOtYrNue4eowrAQDED0UdsisGAxGRn+NaYh8dACA+KOqAkOQ7vzXyc1x169WHPjrd0SbCBACiRVGHyKJKfIkkcXICRK4IBiJic35roV66Evro7FEmRJgAQHQo6hBpVElJkSRuToCIeCCiu/NbQ+2hy9dLV1PjuY/OHmVSWxvLHW0AyASKOkQaVeI5ksTtCRC5IhqIKHR+qyWUHrrc+BLtpfP4j2DfdtX6kIIOAKJDUYfkR5XEdOAh0ty5EOJL2HYFgHihqAOy0kPnc3wJ264AEC8UdUDacudCji9h2xUA4oGiDtHyMr0ak2O9Yp87Z/XO5d7mAyJMACB+KOpSGBUSNcfDqKVMr0Z4rFexvDlL5D10ub1zdiUeA0YvHQDEE0VdSqNCouYoqqSU6dWIplhj2yvXXe+cXYnHgNFLBwDxRFGX0qiQqLmKKonp9GopeXOR9tDZT4nwsXfO7JS3EmECAHFFUZdAoUWFwHPenCXUHjofT4kotuWqODkCAOKHog7RDUDEbNAhEb1yIZ0SUezUCH36iorELK4CQGZQ1MF/bgYgIhh0SHSvXHdTrj6dEpE7NGs/NcJarQMAxAtFHfznZgAiwuO6Etcr53TK1adTInKx5QoA8UZRl6DokVKSPyKRoAGI2PfKuZly9emUiNy+TbZcASDeKOoSFj3iKCoErnvnYt8r150Sp1zzTbayzQoAyUJRl7DoEVdRIUhv75z9thIx2QoA6UBRl7ToEf0/cFPMj5SI2VSrm9652PTKhXhCBJOtAJAOFHVJ4sexWmGJyVSr29652PTKhXhChIUtVwBINoq6JPHjWK2whDDV6iRXLlfieuesLdcAToiwfk+w+ugAAMlGUZdECZ0q9VOie+O8brn6eEJEbi8dcSUAkHwUdUgkt7lyiemdK7bl6tM2a774ktrazP+eAACJR1GHxHOSK5crtr1zhU6H8HnLNTe+JOMLvwCQChR1QQliSjXGU6Vh9cQltjcugtMh8iG+BADSi6IuCLoE8sYKkfItIpU5uWIZmCp1IhM9cRGfDuHkxAhOiQCA9KCoCyp9eHuTSN0uIj39W2WJ01mpUffEJao3rrsAYfvHfNxutbZZ7TgxAgDSi6Iu0O9uj66ZYvClJy72vXFuA4R93m7N3Wa1Y9IVANKJog6R9MqlsieulABhi0/brbnbrHZsuQJAOlHUwVf0yjng8zRrsUBhTokAgOygqEOkvXKJ64nzI6IkQAQKA0B2UdQh0l65RPXExTCiJBeBwgCQXRR1Gcl1C0sme+UijCjJnXIlUBgAsouiLiHoVUtQNIn9PgH30OWbcmW6FQCyiaIuY7luYUl1r5yTaJIItlutl8N0KwBkE0VdxnLdwpLqXjkn0SQWH7dbCRMGABRDURdU71vzNumxfae0NbX48pz0qmVrWzUXYcIAgO5Q1AVQ0I07dpws+dMSv58aGdpWzUWYMACgOxR1Ptu6dWugBV2qe9Uyuq3qZMuVMGEAQHco6gL0/370lPT9RJ2vz5nqXrWot1kj2FZ1s+XKVCsAoBiKugBVV2Y8py2J26whb6t2lzlnn2xlqhUAUAxFHdI/2JDvMYW2WUPaVnWTOVdbK1JVFclLAgAkCEUdsjHYkEsfU1MTWQGXD5lzAIBSUNQhG4MNuSJakcuXNWdhGAIAUAqKOp/17NlT1r29Tp7/779KdaWHFaSs8bJ9an9sjAYbSsmaszAMAQDwiqLOZzqZWltba4YkmFINcPs0JoMNfmXNWRiGAAAkuqi7/fbb5YYbbpA1a9bIiBEj5NZbb5XR+n++Ah566CG54oor5O2335ahQ4fK9ddfL1/96ldDfc2IePs0BoMNXqZZlX7JxVbrAADwolwiNn/+fJk6dapceeWV8vLLL5uibvz48bJu3bq893/hhRfk1FNPlbPPPluWLVsmJ5xwgrm89tprob92dLMKp1VMsUvu9qmXS8wLOt1u3bSp/aLX2V4FAASlrE3PtYrQmDFj5JBDDpHbbrvNXG9tbZXBgwfLRRddJJdffnmX+0+YMEG2bNkiv/nNbzpu+9znPicHHXSQzJkzp9vP19jYKH369JGNGzdKXZ2/wcCWreu3yqI5r0vv+kqprotv0RGb47aGDYt1ceZ24MGidavez77dyvYqACCo+iTS7dfm5mZZunSpTJs2reO28vJyGTdunCxevDjvY/R2Xdmz05W9Rx99VDKhlMGCsBTLgUvY9qnXgQcLOXMAgLBEWtStX79edu7cKf379+90u15fvnx53sdo312+++vt+Wzfvt1c7JVwpgcLwhLDHDg/V+C6G3iwsDIHAMjUoESQZs2aJTNnzoywr6wtmhWwqLECBwBAdoq6hoYGqaiokLVr13a6Xa8PGDAg72P0djf3161d+3atrtRpz16gdCmnqlqkZYvI1mZ/nzuhK2BJWJFjBQ4AkGSRFnWVlZUycuRIWbBggZlgtQYl9PqFF16Y9zFjx441H7/00ks7bvvDH/5gbs+nqqrKXEK1c2d78TVoV5H+Pg9jJHAFLEk9cfTAAQCSKvLtV11Fmzx5sowaNcpk082ePdtMt5555pnm45MmTZJBgwaZbVR1ySWXyJFHHik33nijHHvssfLggw/KSy+9JHfeeafEhlVE9uqViJMOksxJD1yuYity9MABAJIq8qJOI0o++OADmTFjhhl20GiSxx9/vGMYYvXq1WYi1nLooYfKAw88INOnT5fvf//7JnxYJ1+HaSxGXPTuLbJ7pUhtyCuEGeNmCjUXK3IAgLSJPKcubKHk1G0VWbSovbZLwqBq0lbaiuXAOcWKHAAgThKfU4dsK2WlzcKKGwAA7SjqkIi8t0JYcQMAoB1FXYa2K8OicXobNnDiAgAAYaKoC5DtIIvQi6WoNTSIDB/effoKK20AAPiDoi4AWqho4aUDE7qyFkWxFDWKNQAAwkVRF1BM3Zgx7RnEfqNYAgAA+VDUBSTsQywAAEC2fZzqCwAAgMSiqAMAAEgBijoAAIAUoKgDAABIAYo6AACAFKCoAwAASAGKOgAAgBSgqAMAAEgBijoAAIAUoKgDAABIAYo6AACAFKCoAwAASAGKOgAAgBSgqAMAAEiBHpIxbW1t5s/GxsaoXwoAAECnusSqU7zIXFG3adMm8+fgwYOjfikAAABd6pQ+ffqIF2VtpZSECdTa2ir/+Mc/pHfv3lJWVmZuO+SQQ+RPf/pTt491cr9C99EKXAvJd999V+rq6iQNnH7fkvJ5/XheL8/h9jF+vV+LfZz3a/rfr14f7+ZxYfxsVbxf4/1+jepna9Ler1qOaUG32267SXm5t+64zK3U6TfqU5/6VKfbKioqHBVaTu7X3X30Y2kp6px+35Lyef14Xi/P4fYxfr1fnTwP79f0vl+9Pt7N48L82ap4v8bz/RrVz9Ykvl+9rtBZGJQQkQsuuMC3+zl9rjSI6msN6vP68bxensPtY/x6v2bpvap4v/rz/XDzOH62epem92tUP1uz+H7N3PZrVHS5VSvwjRs3pmalDunF+xVJwvsVSdIYYD3ASl1Iqqqq5MorrzR/AnHH+xVJwvsVSVIVYD3ASh0AAEAKsFIHAACQAhR1AAAAKUBRBwAAkAIUdQAAAClAURdD//rXv2TUqFFy0EEHybBhw2Tu3LlRvySgIE1FP+qoo2T//feX4cOHy0MPPcR3C7H1ta99Tfr16ycnn3xy1C8F6OI3v/mN7LPPPjJ06FCZN2+euMX0awzt3LlTtm/fLj179pQtW7aYwu6ll16ST3ziE1G/NKCL999/X9auXWt+CVmzZo2MHDlSVqxYIbW1tXy3EDtPP/20OYrpZz/7mTz88MNRvxygQ0tLi/nleOHChSbHTn+WvvDCC67+389KXQzp8SJa0Ckt7jQfmoxoxNXAgQNNQacGDBggDQ0NsmHDhqhfFpCXrirr2d9A3Lz44otywAEHyKBBg6RXr17yla98RZ544glXz0FR58GiRYvkuOOOM4fulpWVyaOPPtrlPrfffrvsscceUl1dLWPGjDH/WG63YEeMGGHOqf3ud79r/kcJxPX9alm6dKlZadbDqoE4v1eBuL1///GPf5iCzqJ/f++991y9Boo6D3RLVAsu/cfJZ/78+TJ16lSTGP3yyy+b+44fP17WrVvXcR+rXy73ov+oqm/fvvLnP/9ZVq1aJQ888IDZ3gLi+n5Vujo3adIkufPOO/mHQqzfq0Bc378l07Nf4Z1+Cx955JFOt40ePbrtggsu6Li+c+fOtt12261t1qxZnj7Heeed1/bQQw/xz4TYvl+bmpraDj/88LZ7772XfyXE/mfrwoUL20466ST+pRCr9+/zzz/fdsIJJ3R8/JJLLmm7//77XX1eVup81tzcbLagxo0b13FbeXm5ub548WJHz6GrctrIq/TAX13S1WkYII7vV/35dcYZZ8gxxxwjEydO5B8JsX2vAnF+/44ePVpee+01s+W6efNm+d3vfmdW8tzo4fsrz7j169ebnqL+/ft3ul2vL1++3NFzvPPOO/LNb36zY0DioosukgMPPDCgV4ws8+P9+vzzz5ttBY0zsXpI7rvvPt6ziN17Ven/RLW1RbfKtGdZI3jGjh3LvxYif//26NFDbrzxRjn66KOltbVVvve977lOvaCoiyGt1l955ZWoXwbgyOc//3nzAwhIgieffDLqlwAU9O///u/m4hXbrz7TKVWNJMkdbNDrGvcAxAnvVyQF71UkWUNItQFFnc8qKytNYOCCBQs6btNVDL3OEj/ihvcrkoL3KpKsMqTagO1XD7SB8c033+y4rrEjul1aX18vu+++uxlZnjx5sjnqS7dSZ8+ebfo3zjzzTN/+4QDer0gbfrYiyTbHoTYoeW43g3QcXr91uZfJkyd33OfWW29t23333dsqKyvNGPOSJUsifc3ILt6vSAreq0iyhTGoDTj7FQAAIAXoqQMAAEgBijoAAIAUoKgDAABIAYo6AACAFKCoAwAASAGKOgAAgBSgqAMAAEgBijoAAIAUoKgDAABIAYo6AJn19NNPS1lZmfzrX/+SrNPvw6OPPhr1ywBQAoo6AKE744wzTBGRe7Efhu23o446Si699NJOtx166KHy/vvvS58+fQL7vPm+Tvvlhz/8YUnPTSEGwNKj428AEKIvf/nLcvfdd3e67ZOf/GSX+zU3N0tlZWUgr0Gfd8CAARIkLRot8+fPlxkzZsgbb7zRcVuvXr0C/fwAsoOVOgCRqKqqMgWV/VJRUWFW1C688EKzqtbQ0CDjx48397/pppvkwAMPlNraWhk8eLCcf/75snnz5k7P+fzzz5vH9+zZU/r162ce++GHH5qVwWeeeUZuueWWjhWyt99+O+/26//8z//IAQccYF7fHnvsITfeeGOnz6G3XXfddXLWWWdJ7969Zffdd5c777yz4Ndp//p0RVA/n/22Bx98UPbbbz+prq6WfffdV3760592Kmj1ezFw4EDz8SFDhsisWbM6Xof62te+Zp7Tuq5+9atfyWc/+1nzmE9/+tMyc+ZMaWlp6fj4ypUr5YgjjjAf33///eUPf/hDCf+SAOKCog5A7PzsZz8zq2hapM2ZM8fcVl5eLj/5yU/k9ddfNx9/6qmn5Hvf+17HY1555RX5whe+YIqUxYsXy3PPPSfHHXec7Ny50xRzY8eOlXPOOcesnOlFC8NcS5culf/4j/+Qr3/96/KXv/zFbI1eccUVcs8993S6nxZ6o0aNkmXLlpni8rzzzuu0+ubU/fffb1burr32WvnrX/9qikX9fPr1Kf16f/3rX8svf/lL8/x6f6t4+9Of/mT+1NVO/Xqs688++6xMmjRJLrnkEvm///s/+a//+i/z+vVzqNbWVjnxxBPN9/ePf/yj+f5edtllrl87gBhqA4CQTZ48ua2ioqKttra243LyySebjx155JFtBx98cLfP8dBDD7V94hOf6Lh+6qmnth122GEF76/Pe8kll3S6beHChW36Y/DDDz8010877bS2L37xi53u893vfrdt//3377g+ZMiQtm984xsd11tbW9t23XXXtjvuuKPb13z33Xe39enTp+P6Zz7zmbYHHnig032uvvrqtrFjx5q/X3TRRW3HHHOM+Rz56Gt/5JFHOt32hS98oe26667rdNt9993XNnDgQPP33//+9209evRoe++99zo+/rvf/S7vcwFIFnrqAETi6KOPljvuuKPjum6rWkaOHNnl/k8++aTZely+fLk0Njaa7cSmpibZunWr2W7VlbpTTjmlpNekq2XHH398p9sOO+wwmT17tlnx0+1hNXz48I6PW9up69atc/W5tmzZIm+99ZacffbZZgXRol+XNbih28Zf/OIXZZ999jE9iP/2b/8mX/rSl4o+75///GezwmmtzCl97db3Sr9GXaXcbbfdOj6uq5gAko+iDkAktIjba6+9Cn7MTvvftKDRbU4tVurr6832qhZE2nemRV1NTU1Ir1xkl1126XRdCzvd1nTD6gecO3eujBkzptPHrOJR++JWrVolv/vd70xRq1vD48aNk4cffrjo82oPnW6x5tIeOgDpRVEHIPa0102LJu1l0946pX1mdrp6tmDBAlPQ5KM9ZLpiVYwOLOgql51e33vvvTsKLb/079/frJb97W9/k9NPP73g/erq6mTChAnmcvLJJ5sVuw0bNpjCVovL3K9JC0HtvytUMOvX+O6775o+PB3AUEuWLPH1awMQDYo6ALGnBcqOHTvk1ltvNcMP9gEKy7Rp08x0rA4unHvuuaaIW7hwodmS1SlaHTDQwQBd9dMYES2Kcn3729+WQw45RK6++mpTROnAxW233dZpItVPWoBefPHFZrtVi7Xt27fLSy+9ZCZ2p06daiZ+tfA6+OCDTTH70EMPma3evn37msfr16SFrG4R67SuTvzq4IWuaupUrhaB+jjdkn3ttdfkmmuuMSt9WqROnjxZbrjhBrOV/YMf/CCQrw9AuJh+BRB7I0aMMAXO9ddfL8OGDTNToFa0h0ULlSeeeMIUMKNHjzZ9Yhrt0aNH+++u3/nOd8xqm07Hah7e6tWru3weXeXSFUCNGdHPowXSVVddZXrbgjBlyhSZN2+emWDVgvTII480k6p77rmn+bhGpvz4xz82k7ZabGpB+thjj3WsVurKpcaRaI+cFn5KY1x+85vfmO+FPuZzn/uc3HzzzSYOReljH3nkEdm2bZv5PulrsPffAUiuMp2WiPpFAAAAoDSs1AEAAKQARR0AAEAKUNQBAACkAEUdAABAClDUAQAApABFHQAAQApQ1AEAAKQARR0AAEAKUNQBAACkAEUdAABAClDUAQAApABFHQAAgCTf/wciT6Qp1KN4gAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x_vals = estr.theta[n_unique_mxz: ]\n", "plt.step(x_vals, estr.theta[:n_unique_mxz], where='post', color='k')\n", "plt.fill_between(x_vals, cints[:, 0], cints[:, 1], alpha=0.2, color='red', step='post')\n", "plt.fill_between(x_vals, cbands[:, 0], cbands[:, 1], alpha=0.2, color='blue', step='post')\n", "plt.xscale('log')\n", "plt.xlabel(\"Fraction Tested\")\n", "plt.ylabel(\"Hit Enrichment\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "351462c2-a216-43b4-b907-c8181ddeb0f0", "metadata": {}, "source": [ "The estimated hit enrichment curve is comparable to the one reported in Figure 5 of the paper. There might be slight differences due to minor differences in estimators (we ignore ties) and the sup-t confidence bands rely on a resampling process (so there can be Monte-Carlo errors).\n", "\n", "Importantly, we see that the confidence intervals are much narrower than the confidence bands. This illustrates how the confidence *intervals* can be misleading when making inferences for functions. Hence why focusing on the confidence *bands* in settings like this is preferred for statistical inference.\n", "\n", "Now let's compute and plot the functions for the their methods described in the paper and their corresponding confidence bands." ] }, { "cell_type": "code", "execution_count": 9, "id": "9c25560b-1fba-4a8a-b691-8808a872f1de", "metadata": {}, "outputs": [], "source": [ "# Storing outputs from max Z-scores\n", "beta_mxz = estr.theta[:n_unique_mxz]\n", "rho_mxz = estr.theta[n_unique_mxz:]\n", "cb_mxz = cbands" ] }, { "cell_type": "code", "execution_count": 10, "id": "7802f301-a84c-4e18-a393-f410ee90192b", "metadata": {}, "outputs": [], "source": [ "# Computing for Surflex method\n", "unique_srf = np.unique(x_srf[active == 1])[1:] # All unique values for A=1\n", "n_unique_srf = len(unique_srf) # Number of unique values, K\n", "\n", "def psi_edf(theta):\n", " # Separating parameters into subsets\n", " beta = theta[: n_unique_srf]\n", " rho = theta[n_unique_srf: n_unique_srf*2]\n", "\n", " # A fast way to compute the indicator terms without a for-loop\n", " indicator_matrix = (x_srf >= unique_srf[:, None]).astype(int)\n", "\n", " # Recall fractions at all ranks\n", " ef_hedf = active*(indicator_matrix.T - beta).T\n", " # Fraction tests at all ranks\n", " ef_fedf = indicator_matrix.T - rho\n", "\n", " # Stacking the estimating functions\n", " return np.vstack([ef_hedf, ef_fedf.T])\n", "\n", "\n", "# Estimation\n", "init_vals = list(np.linspace(0.99, 0.01, n_unique_srf)) + list(np.linspace(0.99, 0.01, n_unique_srf))\n", "estr = MEstimator(psi_edf, init=init_vals)\n", "estr.estimate(solver='lm')\n", "\n", "# Confidence bands computation\n", "cbands = estr.confidence_bands(method='supt', seed=2101101)[:n_unique_srf, :]\n", "\n", "# Storing outputs\n", "beta_srf = estr.theta[:n_unique_srf]\n", "rho_srf = estr.theta[n_unique_srf:]\n", "cb_srf = cbands" ] }, { "cell_type": "code", "execution_count": 11, "id": "42b20020-4679-46b4-b21c-d705b4dfecfd", "metadata": {}, "outputs": [], "source": [ "# Computing for ICM method\n", "unique_icm = np.unique(x_icm[active == 1])[1:] # All unique values for A=1\n", "n_unique_icm = len(unique_icm) # Number of unique values, K\n", "\n", "def psi_edf(theta):\n", " # Separating parameters into subsets\n", " beta = theta[: n_unique_icm]\n", " rho = theta[n_unique_icm: n_unique_icm*2]\n", "\n", " # A fast way to compute the indicator terms without a for-loop\n", " indicator_matrix = (x_icm >= unique_icm[:, None]).astype(int)\n", "\n", " # Recall fractions at all ranks\n", " ef_hedf = active*(indicator_matrix.T - beta).T\n", " # Fraction tests at all ranks\n", " ef_fedf = indicator_matrix.T - rho\n", "\n", " # Stacking the estimating functions\n", " return np.vstack([ef_hedf, ef_fedf.T])\n", "\n", "\n", "# Estimation\n", "init_vals = list(np.linspace(0.99, 0.01, n_unique_icm)) + list(np.linspace(0.99, 0.01, n_unique_icm))\n", "estr = MEstimator(psi_edf, init=init_vals)\n", "estr.estimate(solver='lm')\n", "\n", "# Confidence bands computation\n", "cbands = estr.confidence_bands(method='supt', seed=2101101)[:n_unique_icm, :]\n", "\n", "# Storing outputs\n", "beta_icm = estr.theta[:n_unique_icm]\n", "rho_icm = estr.theta[n_unique_icm:]\n", "cb_icm = cbands" ] }, { "cell_type": "code", "execution_count": 12, "id": "f40fabfa-8301-45f9-af8c-a0f9caf8dade", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHWCAYAAAARl3+JAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcB5JREFUeJzt/QmYXGWZN/5/z6nu2rrTnYQYAhgIW0ACCRgIAiouKC7DK6MobhDW+QOCOLggI+IgmyMDxiEILzC4jQwIrzr+FRFFEFEIEAIYRpKACIhsMSSddHXt53fdT1V1n6o+VXWq6uz1/VxXXemqruqurnS67zzP974fzTAMA0REREQUarrfT4CIiIiIeseijoiIiCgCWNQRERERRQCLOiIiIqIIYFFHREREFAEs6oiIiIgigEUdERERUQSwqCMiIiKKgAH0mXK5jL/97W+YMWMGNE3z++kQERERtSTnRGzduhU77rgjdL35elzfFXVS0M2fP9/vp0FERETUkeeffx6vf/3rm76/74o6WaGrvTAjIyN+Px0iIiKilsbGxtSCVK2GaabvirralqsUdCzqiIiIKCzaxcbYKEFEREQUASzqiIiIiCKARR0RERFRBLCoIyIiIooAFnVEREREEcCijoiIiCgCWNQRERERRQCLOiIiIqIIYFFHREREFAEs6oiIiIgigEUdERERUQSwqCMiIiKKABZ1RERERBHAoo6IiIgoAgb8fgJEREREoVPOAShVr8QAPeHzE2JRR0REREGTywGlUu0Px8RiQCLhwB1LW4HsGsDIVK5raSB1sO+FHVfqiIiIKDikklu1CvnNGaxfB2RlQcwhyQSwYFdgcHDqNl0D4vGGO6bTwOLFFu+QAi4DlP8XKG8E9NmyZFct7hysPrvEoo6IiIjcZ7HsZrkSl81Cfy2DvBHHuJ7A4Agw4EC1UioBfx8Dxp+tv10W5BYtMtVvhQKwcSPwwAMWHyUPjD4L7DUXSO4A6ENAOQsYeQQBizoiIqI+5/Q25zRbt0J7bA20iep2ZbV2+ssz1itxsVwGuVmjmCjHMWdG/cpaLwaGgXIZdc9hogCUpKBLVm9MJiuV3uQdpWCrvp3PArktAPasFHQBw6KOiIioj1V3O5GZqrccpeVzGF23BvGxjSgMz4YxMFWhZfPAyGglwlZHH0VqMI4h3bmCTlh9LCnsplQLuMld1wKgPyXPtHJ1IF+pQjX/myKssKgjIiLq41W1bLZS0Mn2Y2NvgBRkWrm3T6YbWQzrGZTnzoaenlH3vrTDRVtv8oC+dqqAmyTXRyolkxEDShqgWWTt+r2ou/fee3H55Zdj9erVePHFF/HjH/8YRx99dMvH3HPPPTjnnHPwxBNPYP78+Tj//PNxwgknePaciYiIoraqJh9vdLS+L0AKuvRTqxDL9v7JtHIGxeFRmBbpAqhcLeDkSZqf6HDDdf8bIgJZ1I2Pj2PJkiU46aST8MEPfrDt/Z955hm8//3vx2mnnYYf/OAHuOuuu3DKKadghx12wJFHHunJcyYiIvLLxETzVbVeSEGXQA7IThUsei6rCrryYBzGYI+fbHgUxmAwV7emkwIumNurgS7q3vve96qLXddeey123XVXXHHFFer6G97wBtx33334xje+waKOiIgiv0r3+OPWq2q9Uqtya6evymnZDIxQFWT9LVSZuvvvvx9HHHFE3W2yQveZz3zGt+dERETkRXauln2bPbu7gk4KNzTJxzVdleuDgk7Tat2tDg7E80moirqXXnoJ22+/fd1tcn1sbAwTExNIpVLTHpPL5dSlRu5LREQUxuxcbZXOqZW4fl+V07Q8hkfXQhvIAjqqmTrJ0IVTqIq6blx22WW48MIL/X4aREREHa/KNWbnZqZyiMtqW2ODZhu28nF9VtAJXc9Bj5mbIxqbIsIlVEXdvHnz8PLLL9fdJtdHRkYsV+nEeeedp7plzSt10jVLREQUhlW52larWm1b3303aj+uxLVbpUsObUCxJCdChLuYC2VRd8ghh+D222+vu+1Xv/qVur2ZRCKhLkREREEkK3RWq3LTmiHKpd66Ufu8oDMMYCKn9lgVTZciqIiJidkwjMH6OzY7cDaXlw9SuU/tvuWJ6lFh1dv6tajbtm0bnnpKJjVPjSx59NFHMXv2bOy8885qle2FF17A9773PfV+GWWycuVKfOELX1BjUH7zm9/ghz/8IX7+85/7+FUQERF1PxhYtlqFFHRyQlVdY0PDiBEhBZ2RMN2R2jIM4JMX7YU1GxrzcpVFoY2/+W+kkqXKHU+6CHhsQ+sP+NdfAbJDaGSBv76rctserwAY6t+i7uGHH8bb3/72yeu1bdLly5fjO9/5jhpI/Nxzz02+X8aZSAH3z//8z/jmN7+J17/+9bjhhhs4zoSIiEI9GLixAaLViBFZcaPOTOR0i4LOgqzQtSvoAszXou5tb3sbjBbLlVLYWT1mzZo1Lj8zIiLqZ92uuHU7GLixAaKfR4y47XcrH0MqUYYe24pY/AlkMmmkk42HzwL49UoglZi+/Sp/yelgrpSGKlNHREQUhhW3TubINWuAYGODTYYBrbo13UiWjbI5HbmcjjTG1W1DGEdayyKZ/F8Y2ALN0KBl5VxXtaQ39WAp6FINxZuuAUZJuixMnyE4WNQRERHZaFywo6uTHpo1QHBVrj3DwC4XnYz0hsfb3rVS0gE407nPjZdOQZCwqCMior7UbIu1WeNCL9qd5iDYANHF65rL2iroOrb/QiDZpqKXJon8+srbg7sBWhp+Y1FHRER9p90Wa7cnN/RymgMbIHqzfuWdKCemZtZmC3kcdsbB6u3fXXMzUokSUoni5M6pAR353ID6XtizsYaTK5NbrDa87t87u79LWNQREVHfabfF2tU2ajN25stxq7X5Fmcuq3Jxk+e0anJOa4Wem5h8O6PHUNanGh6y2iAy1REjiZE0ksmSSsCZU3CGVjn1FVILJjp8XsbU5wb8L+gEizoiInK9MzRonNhibbWlasbtVffzcuKwMw6eLOLcZQCbzgAKf0TQsKgjIiJPOkODppctVjtbqnX35/aqq3m5+3AYMrDOtC3d+xW19eoYWS00F3SJ/QAtGCdXsagjIiJPOkODpqct1k6P7OL2ak/m4mX8dMWjmDXrf1EupeuP9QIwK5HAWu1my8dKQeda3G3+nZW9W2MbgoBFHRFRxDm5XepGZ2gQtdta5Zaq+3PnNFNebhxDmLnD00inCigWB6VVuOFBUzk7T2kSxpPZdQgEFnVERBHmxnapk52hQWR3a5Vbqt7m6IzyIIrFORYFnZcMYM4KBBWLOiKiCK/O9XJ0VTOOdoYGcFWu6RFdjbil6kmOTvJy++y9Gcm41nNBVywC5XL9dUuFQv0dJ28fBwZfqLwdXwhoScAwnULhMxZ1REQRX53r5uiqftFsVY5HdAUjRyfbrndevRa37fQjGGXJ0nX/8YpFYPPm6f8O4gkgVpmYMlXQbdxonS/QClNvz7uhMpsuIFuvgkUdEVHEmxlY0LXAI7p8O5+1di6rmS7/I6mSgk5GlKSSZUdqp3K58m9it92BQdOCnxR05uvqjlLQLVo0fXm7nAFeCdZsOjMWdUREEWuE6JdmBieaH9jwELx5c4008+qYAwYH258AZvkPSJYJXzwJQcaijogogo0QUW9mcLL5gQ0PwTuftTZ37o0LxzBzdB00LQsDw/CVYTrrtZanCxgWdUREEZwbxy3XDubKseHB8/NZJ3I63nzmEnXbr674I1KJ+qaEVDyJh7VHkU5lEBvIolwa8bnrtUEtTxcwLOqIiEIsn6/82c9brXaO6+I2a3OyqyhFllv0rFaXkysjhQnopnNZE0gmGztNK6e01uomw/CpXDHkjNdqniGAZ702YlFHRBTirdfHH+/vrdZOjuviNqt1zfLJi/bCmg1ubW0aeARvnLwmq3PenM/q0Ivz0slArrssoB9Y1BERhXzrdfbsPh5X0slxXdxmnUZW6Nwr6IA0MjgAj6q312D/aeezvnHhtmlbr4FhZK0LusSSQObpBIs6IqKQqxvH0KekoDMSwfxFGxa/W/mY4wWWLkd9nVm9ctVKPJysFHg18vl0PQ9o1p9X02wO9i02GRZcnU8ni27lIqDLsOFqZKEpmVNndcarOhJMnlSyPk9nONud2wsWdUREIT2/tTa6pF+0GknSL5zOv5k/lhRYaXO2zWK2XKd0TOXQUkkDRkN2TtPySKbXQo81/zxtO1+LBQxs3ggjnrQs6LZtBeLV//gk40BMnpLlSRKGfJNV7yiFm6mglIJOT1k8JA+UNwH6HJl4B7+xqCMiCvH5rf2Sp+NIEi/yb87NlrNNK6uCTs51NZp0t6qCrlXna7msCrqJ3Reh3LAFX8gD2XFgt32mTo4YjDd5cTedART+OHXbyzaev1EGtDSQXAzoDp3D1wMWdUREIR1b0lejSziSxNX8W2O2rZvZcq1kFi5puT0uBZ1h9FYUlaWga/gc8hUZRSA+KiuFrR48UV/QdZqh04Lxj5BFHRFRSPTD2JJm40k4ksTd/Jt8rGZj12qz5XqhCroAznWzZM7PWWXoAoxFHRERhWI8Sb+PJDGffTot/9bzBzegmUKaqsGhSgo6I9lbURcqWpP8XAiwqCMionCMJ+njkSQS+Tru4r18yc9JMwN0d1aqbHe3ki0s6oiIQnJqRNTpkmrneJKmebonn63MeNt7l4yjW6+t8nOZhYuRmPUU9AH3iq9AnOsaESzqiIgCrF9OjZCt1+SGx/t+i9WO758vB9y78/cwLT+XAtIDa1p2p/aqbXdrsYBivqxWK61oxYKaQVeQ6SI2Rs45RsaZGMFaaWRRR0QUYH1zakR167U0Mrtvt1hbzaMzv91VPddi5pxVfk5tuWrSPJFzrDu1KzJU+O8bsS2bnJw1Z6UcT2J8Qledro2kuUhGmbQ827XuXFebBV3xxco4E7kEYEadYFFHRBQCUT01otbtOtndOhDRL9TPeXQdzpxrHAjs6/ZouYzyQBKFPRZht8UJNWvOkq43/V9PrNm7ejnbtTafLrUM0IYCMaNOsKgjIqJAdLv2e3ernXl03ZyVanfmXG2WnKzOmQcCt90e9UBpIIH4aLL1rDknznbt9FxXuW9ACjrBoo6IiILR7drH3a1259G1midnR6uZc42z5HzbcvVzNp0Wnpl0VljUEREFWJQ6XxsHC3OgcPv8XN08ui7PYu105pymBeeAencZkZhNZ8aijogooKLU+dpssHC/b7nazs95dBar5OkSyQ3RHzNiSJ7uFEQNizoiooCKVOdrs8HCfb7lajc/58RZrO3OX618orLK05VLI/7m6KTrVZoksgUULDpaHcnT5ddX3o4v7CxHF2As6oiIAi5Kna9S0LUtLPqU3fxct2exdnL+qmH4WB4UCxjYvBEFPYmxrcDAUBLJtN58LEmv5t3QeY7OCOYWNYs6IqKAbr2ajuIMdXbOnJ/rp2ycHXbzc311Fmu5DEPmzs1fhGw+gX3205Eaibu4Wq11PqOuvAnQ5wRmPl0NizoiogAWdKtWVbZew5ana5adi3p+ztHZch7l50KxqqslER8OWPzAqM6oSy4O1DgTwaKOiCigWTr5RSYFXaB+oXWbnYt4fq5VNs4Ou/k5W7m4HpokaidIkA1a8L6XWdQREQV02zWRCFlBJ9uEcgBnn2fnGrNxdtjNz3WSi+uE+RSJyHe+RhiLOiKigAjztmtt6zW54fHQb7M6mo2z+Ql9z89Vu17lFImS4f8JEp7OqIsQFnVERAER6m1X09ZraWR2aLdZHT93NWT5OTlFwrWCrjqmpO6mYuU1N9OKBZSL1cHbbnS8GtGcUSdY1BERBUwYt13NjIHwrvI4mY2zw+v8nGyzyqrc9Ntznowpka7WyZuKwLatQNzi26UcT6JQ0pFMw/lRJkaPM+oCOs5EsKgjIgrAtqus0rk9wsRq1IiToja2xMlsnB1u5+fMuTnr97uYpauOKZnYfRHK1QYaiV9mx4Hd9gHijU2kuq7+ZyMFnav/wZnX4Yy6AI8zESzqiIgCkqMTbmXpWo0acfTzhCRP1yw311M2rsUna3Zmq6f5OVNuTm2zNj5NuJ+lUwVddfVRXlmjCMRHgVTSpb9kw+J1NyZ6mFEX3HEmgkUdEVFAcnSy7epalq7VqBEnhWBsiae5OR8yc+22WKWgM4zgFSTO5+ZOBnKP9804E8GijojIh63WGvP4kmTSva3V2tZoP48a6SQ312k2rhm7Z7Y6lZ/zdYs1SIxs+4IusSQyZ77WsKgjIvJpq7Wm1y1Xu1urYdkaDUJurpdsXDOtzmx1LD8XgC3WwJl/J6BZvO6aOzP//MSijojIp63Wmp63XO1urTqwNdrNDLegcSQ31yInZ+bXma2dbrFKJ2rDtBFH6dkiyoVKc0Tt0xQKDmXkLO87MfW2FHR6D6+7NEdIlk69HewTN1jUERF5rNet1mbc3lr1fIZbUPk8W65ZZq7yvs6LDinoNm92r8tU5s4Nbh1DcWgmxid01RxRI/8ObI0scTsj16qgK75YaY6oUW8Hr/NVsKgjIvL4+C+n1HJ0Xo0S6XWGW9B0m5uzm5PzIzPXTW5OVuikoNttd2DQjZ3ZfBlaJonyoj2BGfWVo+2RJXYycm7k5mrdrqllpo8TC2Tnq2BRR0QUwuO/GnN0XuflupnhFjRO5OZa5eT8yMz1kpuTgi7pVq0iq3Mz5BO4mJGz4lRuTj6OblqtCygWdUREYTz+qzFH5/EoEcdmuIVcNzm5Vtun7R/bR2NJmuk1IxdhLOqIiFymzrC0efyX3VMf7I4ocbKxIewNEk41QZibH9zYPm3/MfpkLEldg0T3r3k/YVFHROTy1uvjj9vbdu301Id2W65sbAhgE4SN7dN2isVhlErOhd+kUaIpaVHttS224zbXADRIhBSLOiIiD7ZeZ8+2se3a6akPbbZc3WpscGowbxB10gTRS/NDt9unbnWqytmr07pQpRjbuNGZVm35GHKea6caGyTcHhhsmMaXhGCESSMWdUREHuikq9CN0SRONja4MZg3iNo1QTjW/NABtzpVpaCb9vHkk0kxtmhR/WDFbkhB12slKg0S+iz3XnPDYnxJwEeYNGJRR0Tk4lFg7caYmDN0To0mqeXoXDmcPuL5Ob+GBQeqU9WrwYrdNEi4WUQbVuNLgj3CJHBF3dVXX43LL78cL730EpYsWYKrrroKy5Yta3r/FStW4JprrsFzzz2HOXPm4JhjjsFll12GZBC+4Yio71kdBdYsT2eVoet1NAlzdOEZIkwBpYVjfEngirpbbrkF55xzDq699locfPDBqmA78sgjsW7dOsydO3fa/W+66SZ88YtfxI033ohDDz0U69evxwknnABN03DllVf68jUQEbU7CqzpGBOrDF2Po0mscnRRzsC5lZ9zalhw3efT8l2d+EAUiqJOCrFTTz0VJ554orouxd3Pf/5zVbRJ8dboD3/4Aw477DB8/OMfV9cXLFiAj33sY1gl/y0mIgrY+BLZQJjcXs12P5ak1xxdv2TgnMzPOZ2XM48y6WUkSW1b3xWNna69dK1SfxV1+Xweq1evxnnnnTd5m67rOOKII3D//fdbPkZW5/7rv/4LDz74oNqi/fOf/4zbb78dxx13nIfPnIjI3vgSOyNKnNhuNWfnmKOzP3vO0/ycaZRJyejuxAfpfB0fB4Zn2DwvtRPNOl277Vp1jOHj5w4f34q6jRs3olQqYfvtt6+7Xa4/+eSTlo+RFTp53Jvf/GYYhoFisYjTTjsN//Iv/9L08+RyOXWpGRsbc/CrICJqMb4ka2NESQ/brczPhS87p2bTdTmfrtb5uvN8F85obdbp6kTXak8z6k7x53OHVKjGg99zzz249NJL8a1vfQuPPPIIfvSjH6nt2osuuqjpY6SJYnR0dPIyf/58T58zEfWfxl+4te1Vy4vD+bka5ujsz55zIz/XSNOc28qMubkcU8sN1C5+FXS1GXX59ZW34wvdnU8XEb6t1EnnaiwWw8svv1x3u1yfN2+e5WO+/OUvq63WU06pVO777bcfxsfH8U//9E/40pe+pLZvG8n2rjRjmFfqWNgRkdNyYzlkx0rQs4A+CGgOjijpZg4dc3T2Z8+5PW9O8nSJ5Ib+O97LSfNu8HwmYBj5VtTF43EsXboUd911F44++mh1W7lcVtfPPPNMy8dkMplphZsUhkK2Y60kEgl1ISJys6D703dXobAlg3QOGJ01tVrXa2auFfNPvb6dQ9fmnNZAzJ6r5unKpZGut16JBV3gu19lBW358uU48MADVeODjDSRlbdaN+zxxx+PnXbaSW2hiqOOOkp1zB5wwAFqBMpTTz2lVu/k9lpxR0TktVK+pAq6WCqO0TkJtVI32aTY44iSZuT/scddvBf6WoCycnYYxkDHjRHmZtSWZ7Ta0eocV3a6RoKvRd2xxx6LV199FRdccIEaPrz//vvjjjvumGyekAHD5pW5888/X82kkz9feOEFvO51r1MF3SWXXOLjV0FEVBFLJTAwnPSkX0/ydE8+WxmQuvcumb6cQ9fJOa2dZudky1RW2JzQzWy6Zme8Wp7Raoedc1x973T1iVE97zVk57xa0Yxm+5YRJZk6aZjYsmULRkZG/H46RBTCo78aZTdl8PR37kV89gzER7wJc2eyOg489QD19kPXr8FQH269atkJ7H3qW2yd09pJds48U86x56plUSzOsb39KvMO5Yi5xjNeLc9otfsNLK3Zrc5x9bPT1Up5Aniu8veLnX8H6Cn3z3vV5JiwgwN3LJjd2sX3Y8KIiMJ49JeZNEdIli7p0iJH4yw6Yb7eF2kjL+fMmWbKqREkDlANEl18LMfPeA3KOa7mb27pcrV834QHn7/ccN5reM55tcKijoioi6O/zCRDJ80R8qfTWx+cRedfdk4KOsMI7y/4wFNz6E4GcgHIRGrhPe/VjEUdEVGPixxadUXFjROcWs2i65d5dF7OmeP5rB6SFTo7BV1iCWfU2cSijoioTY5Osk2t6IXqga8ez6Lrx3l0bs6Zc+p81m7Y7nRt1cHaStC7W+ffCWip5qto/fRN3gMWdURENnJ0tfNcG8n5rskNj7s6j65vZtE1mTnn2Zw5B85n7YbtTlc7HaytBLm7VQo6Nxoh2jECXux2iEUdEZGNHJ0UdJaNgeXK+a6lkdmuzKPrGwGYOVc7yquX81l7OdO1badrs/NZ7Qpad6vfjDxQ3gTocyoNEhHAoo6IqMlICTvNgrWtV2OAJwW4PXPOzTNag3CUl+1O16B1sIaVUe18TS4OdcerGYs6IiKLrdfHH2++5erH1ms/aTZzztUzWnmUV//SorN6yaKOiKjJ1uvs2W12qzzYeu2f6fCGv+ezdnmUV+S1miPX88f2YA5dn+F3LxFRE3Yn97u19do357saBhZcfIrrn6bV0V/dHOVlp4vV7mNsdbp63cEapDlyZAuLOiKihmPA2o0wqW296hadmk7ql/NdJU+XfHa9eju7y0JXcnN2jv7qNU/XrIvVDtudrl52sNqdI9crzqFzDIs6Iup7VseAtcrTSUGXXrtKbb16laf7/vnr+mJU11/Ov8Gd3JyNo7+6PcqrXRerHbY7Xf3qYG01R65XnEPnGBZ1RNT3rI4BazrCxJSlKw/GYQyPejLKRIvwXDrzHDq3v1Ivjv5y9LxWPzpdazk6c+bNrzly1BEWdUREXf7+NAYTnE0X4Ll05gydU5m5yGOOLtRY1BERUSDm0jl9hmtjhq6TzJyjDQ92+N0U0SpHx8xbaLCoIyKiQMylc3QOnUWGzm5mztGGBzuC0BTRKkcX1cybEa0jwgSLOiKigDIifK6rV+e5dpOhc7Thwe4nDExThNEfOTojekeECRZ1REQBFKkZdS7m55rNnnMiQ+dowwNCcPyXytO5Py8wEIzoHREmWNQREfp9nInduXTS9Srcnk8XtRl1rc517SVH1272nJ/nuIaS5OnylXmBiC+sbLtGnRadI8IEizoi6lvm+XR259JN3ubhea9RmlHXeK5rTzm6NrPnep0719fmuTQvkFzFoo6I+pZ5Pp3tuXSD1a0aF+fTyS6YrNTVaCHOznmRn+t29lyzDteeuljtdrc2vi9wQvVdR1Us6oio70mUyU4eXc2lc+EIq7rPYQCfvGgvrNkQwm3DEM2ea9fh2lUXa6fdrUHqdI0qI1/Jz027PZpzC1nUEVFfapWlM+fnvMrQ1cgKnbmge+PCbaHJ07XKzjkxh67X2XOddLh21cXaaXermV/Hf0W9oCu+WGmIsKJuj07nq2BRR0R9p1WWzio/53WGruZ3Kx/D7JFiKKNNjdk50fMcuh5mzwWiw9Xv7tZ+Y1Q7XFPLmjR9xCLV+SpY1BFR32mZpbPKzwkPznhtzNLJCl24CjojsLPnqO2LWn/Wa5RosrXdZLUuYljUEVHfapWl8yI/F9YsnWXm3zCw8KtTM87yMttV62y2XDu6nkNePnex8nr1wpFmiFbND433CzKe9xoZLOqIiAIgLFm6Zpl/yR2mn6/MONu200JsKyaBqVhiXUE3PNp8tlz7z59FITfcc1HXczOEneaHsDRCNJ73GuSzXps1PljeN4d+w6KOiChggpylM2f+pSiaZNq5i337BhyQbvbky9AGpKAb7DILJ4WvMx0MPTVD2Gl+CGMjhJz3qs8K5oy6do0PfdIM0QqLOiIiH9VydK5n6Szmx3VDz8uqHBAvAynzapkpj5VKaUCrhR71pUo1FYFcXJiaH1RuzuJ7wJylk/Neg1jQ2Wp8QF80Q7TCoo6IyCee5eg8mB9nLa9W5ur135ZYIEQpN9dHjQ+dYlFHRH1HQvy18SVaQ/DLz5l0bmXp2s2Pc9SSJdWVK1nSWwvA6vWU24LfEBIpjbk5K0HO0pEtLOqIqO9m1D3+OJDdksNOf1+FZKl+Hp2fM+mkmHNj67Vg6vRce+X0+XEdfaxC5bJkfyBl9ftfCjr1BUhhWs3OTcvAOZeLc027ztagd7S2y83JNmsjKegCu/UqDRJc5W2HRR0R9eWMutmjJSRfsZhH5+JMusY5dI05unSy7HheTgq6zS9OZaa2lVIol3qYH6cDyVEgNiTto3YeEMLsnN3O1iB2tNrNzenuzBB0vUGizxofOsWijoj6Uq3r0at5dK7l5zrMy8kKG1K9d422b+QM8UqW3c7WoHW0Rik317RBYqivGh86xaKOiPryvFetUA3W+ZifcyJH11FebskSpGbKFhtcJnm6DeHPzoWps7UfcnOqQYIFXSss6oio7857lTzd0N8f9yU7Z87P1TiVo7M6b1WaQuRrVhm4mV5lpmp5upHgZ+eiKoy5OeoZizoi6rvzXmfNreTpSiOzXT/PtdmZrl3n50wfVFbp9NyE5XmrtZx/QZfbq1uunv8uD8mvmMamiCA3QTTLzIU9N0f99C+OiMi5HbXa2pExMBjOLF2bHF1jzl/+7Po4rK5m0uXC3xQR1CaIKGbmyDEs6ogounI55DIltUonWTpdLoOArvk3i86JOXRWObrMwiWTDR+NR3nZa2zohdVMupDk6Zo1RQStCcJuZi7suTnqCYs6IoqmXA75+1Zh/SMZZKsLR+kcMDqr0vnq1yw6p890reXoVEHX8IGloLOcJec4q5l0IZhFF+amiGaZOcHcXN9iUUdE0VQqobwtg4lSHLGRBAYGgKReWakruTiLrpH5eNS2DRE2z2dtlqPzXwhn0gUpE9f2sczMUWss6ogo0sqDCcSHkmrVqrHI8uL393EX7+X6+az+5fxrObqQZOisTonwsimCmThyGYs6IiIX83RPPls5eHzvXTIts3TdnM8qObq8nmya83e3OaIxRxfwDF2rUyK8aoqwm4lrp98yc0aAu5EDhkUdEZEHvn/+OttZOqt5c1YkR1fOa3VNETXuN0c05ugCnqFrdUqEH00RrTJx7fRTZk6OCCtvAvQ5PB7MBhZ1REQuzaQzz6bT2mToes3JedcUEfIcXVAaIjhHrrMjwpKLeZqEDSzqiIj8mEnXQ4bOXyHL0VE0aAEbLxNQHYcIvvrVryIjY9kbTExMqPcREfklN5ZDZmOmctmURd7juqPZ+a5Ws+maZejM8+bsKBbhoWqOTn8E0J+obr/qwcjLyTlwrS5BPiWCyK+VugsvvBCnnXYa0ulK+LdGCj153wUXXODUcyMi6qig+9N3V6GwZeo/nbFcBhPpUaR9qDvM57u2G2ViztBZzZtrRuqUsTFg5kwvTowIaI6uVQNEoyCeEkHkZ1FnSAbE4gfOY489htmzZzv1vIiIOlLKl1RBF0vFEUtVM176KLZLxNWwYS80zqRLJ0qVzFzOnVlztez/nnt6nfMPUI6uVQNEIz9PiVDz6Uxz5oj8LOpmzZqlijm5LFy4sK6wK5VK2LZtm1rBIyLykxR08RHvg/DTZtJ5mJlz/wjbEOTogtIAYYXz6ShoRd2KFSvUKt1JJ52ktllHR6eO14nH41iwYAEOOeQQt54nEVEdiUnJma61K7kt3p3namcmXVrL2CroOs3QeS9k8+iCqHE+Xb/NmaPgFXXLly9Xf+6666449NBDMejVfgYRkUVBt2qVZHkBLZ/D6LpVKj9XHMsgOdfb81ybzqSzOXfObobOv8MQApSj8/tECKfm0+mz+mfOHAU7U3f44YejXC5j/fr1eOWVV9TbZm9961udfH5ERNPICp0UdBKPSg2WMEPPoDwjDswdxUDa+8yUUTZUo0Ya4+p6LDcBHc6dz9ruMARvmiR8ztEF4UQIJ9KWMp+OBV3Dy5OvzKOzfOkCvOUfhaLugQcewMc//nE8++yzajvWTHJ2kq8jIvIqRiVlhhR3pXQCRsKfgi7//zsTh2ZXVUs6AGe61wtgPjWi85Mjatm4TgTkl2rQToToKE93it/PItgFXfHFyoDhZtT7Yl4+q/4p6qQZ4sADD8TPf/5z7LDDDpadsEREruXnZDMwW9l21bMl6Jq/Wbrc1jwOyK7yJDPX26kRjdm4TgQoRxfkhohmebr8+srb8YXM0jU7MSK1rMVrE+NpEm4VdRs2bMBtt92GPfbYo9OHEhH1lJ+rkYIu+fgqjKYzalyJls0Aw/5n6VZd9huMzBmsy9N1MnfO22xcJwIwjy4K5t0QkO+FAJKCTm+xWkfuFHUHH3wwnnrqKRZ1ROR5fq626yYrdFLQxdJxlAYTqqAzBv3ffovPkFWkRN28uuAJ0Iw5u80QtdtDjwUdBayoO+uss/DZz34WL730Evbbb79pXbCLFy928vkREU3bdZNfjfKjRwo6L8aBSCxKRpZYyTW5nVw4HcKrhgg1KNihbX0OHKYgF3Uf+tCH1J8yr65GcnW1kyY6bZS4+uqrcfnll6siccmSJbjqqquwbNmypvffvHkzvvSlL+FHP/oRNm3ahF122UXN0Hvf+97X6ZdCRGTr9/snL9rL8kxXIR2vkw0Snuuk8SEgDQ+9nA7hRUMEBwVTPxV1zzzzjGOf/JZbbsE555yDa6+9Vm3rSnF25JFHYt26dZg7d+60++fzebzrXe9S75Nc30477aS6cGfKwYdEFEn5/FSOTkPlP426HL3lEVmha1bQNUpWz3r1RjeNDwFqeAhqM0TjoGCncOAwBbGok5Uxp1x55ZU49dRTceKJJ6rrUtxJV+2NN96IL37xi9PuL7fL6twf/vCHyW1fOcmCiKLbJPH440B2Sw47/X0VkqWpbgk/miN+t/IxdabrtDNcz5zaFjYC3fjAhoeOBwXLXDmnGgHYJEEu6yqc8P3vfx+HHXYYdtxxR7VSJmSV7X/+539sfwxZdVu9ejWOOOKIqSej6+r6/fffb/mYn/70p+oosk996lPYfvvtse++++LSSy/lbDyiiDdJzB4tqYKuPBhHKT1DXYpzdnC1OUJ24TJZvS5LJwVdOll/SSXcn81ZLNppfLBzCXgHa+sv1KMs3dTQaFXQ6Q5dWNBREFfqrrnmGlxwwQX4zGc+g0suuWSyoJItUCnsPvCBD9j6OBs3blSPleLMTK4/+eSTlo/585//jN/85jf4xCc+gdtvv1114Z5xxhkoFAr4yle+YvmYnJwJKf/drxobG+vgqyWiIKj1YxkeNka0ytGZ77jg4lNc7x3Yui2PmTPLDSdHhCAj1+kXKj+fJU7jx+kQzNJ1d+JDz697xL6Pw1bUSSPD9ddfj6OPPhpf+9rXJm+XgcSf+9zn4CY5kkzydNdddx1isRiWLl2KF154QTVaNCvqLrvsMlx44YWuPi8icp78X2xyyLBTnYg95OjeuHDbtK1XLZdF8tnKYNnsLgtdKTgNI485c9di992ziCdDmpHrpElizz39OR2iMUvHDJz9Ex96xRMj/G2UOOCAA6bdnkgkMD5uvwdszpw5qjB7+eWX626X6/PmzbN8jJxgIVk6eVzNG97wBtU5K9u5cYsfBOedd55qxjCv1M2fP9/28yQi/wYOT2yuDBmekc5AK/kzYLiWo5NLqx20v5zv0mBZrQw9loUes8rORTAj1zAmy7csnT6LW6a2T3zoFU+McErHa9y77rorHn300Wm333HHHarAsksKMFlpu+uuu+pW4uS65OasSI5PtlzlfjXr169XxZ5VQVcrNkdGRuouRBSOLF1ioITXVYcMu52hMzM3O9RydNb1muHdYFnDKjsXgAIoMoz6LB0zcNYnPrhyCfAw7Kiv1MmqlzQqZLNZNZvuwQcfxH//93+rbc4bbrih44+1fPlytXUrs+kkkyerfbVu2OOPP16NLZGPLU4//XSsXLkSZ599thqCLEeWSaPEpz/96U6/DCJCOCZbTA4Z9qqgM4DjLt7L1h3dzNNpWr6ySqdHPHNUO0HCzxMjVJ7O3WwkUSCLulNOOQWpVArnn38+MpkMPv7xj6su2G9+85v46Ec/2tHHOvbYY/Hqq6+qxgvZQt1///3Vil+teeK5555THbE1sm36y1/+Ev/8z/+sTq6Qgk8KvHPPPbfTL4OIAiy/NaeOAvM6S1fL0z35bCU/tPcumWk5Oi/ydFLQJdNr1bZrviD1ToSyc61OkPDqxIhG8n2Wr/xdIr7QxW3GEDZHsJEhVDRDltu6JEXdtm3bLAcFB5Vk6kZHR7FlyxZuxRIFUG4shz99dxW2vZLB7FlAvJTxdOtVxpgceGolN/zQ9WswlGxS1GUnsPepb1FvP3n9vTCSzgXJNT2L9PAjMMqDyOUGkc3qWLjXoBwrG73wpOyz106Q8OLECCvlCeC5yt8ldr6XB8s3NkeoTN3B3CYNQe3S8UqdWTqdVhciIqeU8iUUtmQwNCsOfTSBoj7qekFnPtvVPJeuaUrOMCpDh9vfs8fnNQjDSKjnF2lenCDR6jzXuvNZXc5GhrI5go0MYdFxUff3v/9dbZfefffdeOWVV+qaFoSc+EBE1Cs9GbCZdKYH7HLRyUhvcOEoKXIHZ9D11hxB0S3qjjvuONWBevLJJ6vsm8YOISLqYfetOr988obcFv9n0jWbS1fL0pkLuszCJY4Un7XGiMrbU80Rda9PVJoizNeDdJ4rZ9NVXy8fm1bI26Lud7/7He677z4sWbKkt89MRH2tNotOIlVChgyPrluFWC6D4lgGybn+zaQT7ebSifUr70RppPd5ZubGiKnbsigWhyHjP4dnoOE0iQg0RdR43RzR6jxXns9aaZAobwL0OZVtV4p2Ubf33ntjYsKcPyAi6n4WneTiVUY+W8IMPYPyjDgwdxQDae8D87WZdNMYhlqhE+YsXTnh0Dyz6oBhaYyQHJ36lBhGqTSoXp+d5wdjJq8jJ0bUmiJq3GqOMGforM5zpdZ5uuRiNkb0Q1H3rW99C1/84hdVrm7fffdVJzyYcbgvEXWTkZfSSH63l9KSpfOhAzIAGbpaY0SjWE8tbX3aFPHSyfa2XMmaFqB/g2Rbxz8qZs6cqVpr3/GOd9TdLpNRJF9XilQAhIj6XWOGzuksnfoc2vQMU7FYuUQiO1e7zSvNMnTMzFHEdVzUfeITn1CrczfddBMbJYior0iGTm25St0gBZ0DW6+Sp0skN6gMnWy5CinmNm+urFzGEyHL0zXLzvk1XNicoWNmjiKu46Ju7dq1WLNmDfbay8YxOkREAWaeT2eVn7PK0BlJh/NY1TxduTRSPd+1ssglBd1uu0MNHA5Vnq5Zds6N/Fyz2XPM0FGf6riok3Nan3/+eRZ1RBRqlvPpfJxBZxjTfxxLMReqgs7L7Bxzc0S9F3VnnXWWOm/185//PPbbb79pjRJyJisRUdA1zqeTuXRpLdO0oHMyQ9cuT0cOzZ5jho76TMdF3bHHHqv+POmkkyZvkwYJNkoQUScDh7MNu2Z6Ie/bCyjz6WaPFGGa+VuXn3MyQ2cnTxfaBgnhx5NvNnuOGTrqMx0Xdc8884w7z4SI+m7gsPw5OloZPJzc8Di0bAYYHvUmS5fVkMa4uj6Ecei5svv5uTZ5ulA3SNSaJMbGZEyCOw0RnD3n/uBhw/Q/G4p+UbfLLru480yIqO8GDktBp3Lz2RJi2QxKI7NhDMbdz9J9dSGueuo9GMcfKjee6eqntPGcBsLfIGFukthzT+cHCjND535BV3yxMnhYLjxNIpS6Gmm5YcMG3H333XjllVdQbphFJEOJiYi6ydAbA4OeZOnWPaXjsFpB52F+zq5QN0gIN548Z895c5JEahmgDfE0iX4p6q6//nqcfvrpmDNnDubNm6fydDXyNos6IgqTdSvvhGHKzrmZn5Pt1vrbuNXVFc6ec4/kEPXpp5pQRIu6iy++GJdccgnOPfdcd54REUVSbiyH7FgJehbQByvHgtXoprlwXtBgTL5teJCdk4IumV6r8nPT35dFsSjnvHrQY2B10oPTH7/b2XJtH8fzW4kcL+pee+01fPjDH+70YUTU5wXdn767CoUtGaRzwOis6Tt0XjVJSFHxO7wFnqo2RBhlOdu1/guXgu611wYnI2iuNUi0OunBSa1OjWAuLiDNEBaFPRsk+rOok4LuzjvvxGmnnebOMyKiyCnlS6qgi6XiGJ2TUCt1006JHh51vUlC6PksDsCj6u3M/L08zc5JQWcY9VtbskJXa46QQlcKOlfydK1OenBSq1Mj7MyWa4ez55xphrDCBon+KOr+4z/+Y/LtPfbYA1/+8pfxwAMPWA4f/vSnP+38sySiSIilEhgYTpo2P/217gv/iaRD57c25uU6zc7Jj1LpeA39SQ+9zpZrh7PnHGqGsPoeiDFP1w9F3Te+8Y2668PDw/jtb3+rLmbSKMGijogm5XLIZUrIbfE2M2d1hmuNJsPyJq9orubl6u9XGTAsuTlzrM2xHF27vJydvJuTrLJzzMUFqBmiyWodRb+o48BhIupYLof8fauw/pEMsjmgOJZBcq4HmTmvz3BtkZere0pS0BUGJ4cLm/Wco7Obl2uVd3MSs3P+Ym6ub3U1p46IqK1SCeVtGUyU4oiNJDA8dxQDafczc0JW6NoVdKsThyA1I+5qXq6RebiwObnSc47Obl6uVd7NSe2yc8zFufjaMzfXzzou6j70oQ9h2bJl00aafP3rX8dDDz2EW2+91cnnR0QhVx5MID6UxIBPo6/MZ7jK4OE3n7lEvf3bq9ZB0/1J97mWnwtKXq5ddo65OPcwN9fXOl6Hv/fee/G+971v2u3vfe971fuIiHwhGbrshLpYneEqf45jCJnqxTw4vZc8HQcIW/5lmF6kFKA3XBwe7EwtcnPTLhwsHGUdr9Rt27YNcYvle+mCHZODnImor+fRyfgSJZtF0asDE9pk6NR5rxfthTUbhh37lOYGiVoTRCvSEOH4cOFac4TXTRBt83Sn+P0siPpSx0WdjDG55ZZbph0HdvPNN2OfffZx8rkRUUgHDNfEchlMpEeR1v3J0NXOcJVtV3NB98aF25BKlB1rkCgZw0CLJgkp5moNEo4NF25sjvCqCcJOni6/vvJ2fGGT0RnkygBh9T4ePdfPOi7qZEbdBz/4QTz99NN4xzveoW6766678N///d/M0xH1MfOAYZlHp+ij2C4R9/RwenOGzuoM19+tfAyzR4qO7QCqjtcWBV1jg4Rk6Rx5PRqbI7xqgujEvBu41ep1I4TgEOG+1XFRd9RRR+EnP/kJLr30Utx2221IpVJYvHgxfv3rX+Pwww9351kSUWhIQRcf8XB1xjAsM3TVd2Eiq6uVuhpZoeuloKsNGu4mSyfFnOMFbhCaI8wz6cyz6OpO+CVvBggLDhHuV12NNHn/+9+vLkTU38wZOr8GDDfL0rmdo6tcb5+lqx0F5tgQYfP9goAz6fzBAcJkgXPqiMixDJ2nA4YtsnS1DJ1ozNE5kqVrGDSsCro2W6+SpxsfB4Zn2MjS2R0iXBOEHF2zmXScRdddJq7t683MHPVY1M2ePRvr16/HnDlzMGvWrJajADZt2mTnQxJRBDN0SQ8HDFtl6UojsywzXJKjk2Ku163XTgYNN+bpdp5vY+vV7hDhmqDl6Mwz6TiLrvtMXDvMzFGvZ7/OmDFDvb1ixQo7DyGiPuF5hs50rmtjlq5ZxSbFXDrZY7drj2IDIcvJ2c3R8TxXFzJx7TAzR9Zs/ZhZvny5+rNYlI4xDUceeSS23357Ow8lIgrvua4ONEdEEnN0vWMmjvzO1A0MDOC0007Dn/70JzeeCxFRVzPpzFk6p3XTHCE5ulqvQ9uBw+bGiKA0P3STo2OGjih8jRJy7uuaNWuwyy67uPOMiIg6nElnNY/OMR02R5gHDdc0HThs1RgRhOaHbnJ0zNARha+oO+OMM/DZz34Wf/3rX7F06VIMDQ3VvV9m1hERuZGhEy1n0pnm0ZnfduYp2GuOMA8arjVGSEFn2SRh1RgRtOYHM+boiKJV1H30ox9Vf37605+evE1ydob80NU0lDoayERE1HuGzo2ZdL2SIk5Oj4hUY8RLJ1uPLyGicBZ1zzzzjDvPhIhCNaPOq2HD7c51bTaTrnE2Xa3ZoePP32FzxLT/17YaJhyWDJ1gjo4oekUds3RE/c08dNjrYcPtznU1z6Srkbd1vb7ZoVN2T46YNmjYzjDhsGXoBHN0PQwdZgc1BexEiQ0bNuDuu+/GK6+8gnLD/0AvuOACp54bEQV86LDXw4ZrGbpW+TnLmXQNzQ6dsnNyhOWg4ZyNYcJBztAJ5uicHzrM4cEUlKLu+uuvx+mnn65Ol5g3b17d6RLyNos6ov4ZOuzH6RG95Oc6OQnC0UHDYcjMWWGOzqWhw0OA7v73IfWfjou6iy++GJdccgnOPfdcd54REVEL7fJz6WQWml6/UsehwV1ijs6locMs6CggRd1rr72GD3/4w+48GyIKNC8bJOxozM9JQZcass7O2c3FdauYLaCcK0OXYcP5EDZCtMMcHVH0ijop6O688051sgQR9Q9fGiTUfLqJulW6CTTPz8kKXbPsnN1cXLcF3cTzGxEbSiIZB2LylIsBa4SoZeM6eszUa68GDOuVJhXq5vWPUIFP0Snq9thjD3z5y1/GAw88gP322w+DDRM1zfPriCg6PG+QsJhP9+YzlyCDocBk5yY/X6msCrr571mE5EgCg/GANUIwG+fz658HypsAfY4kLn1+MhRlHRd11113HYaHh/Hb3/5WXcykUYJFHVG0edUg0Tif7j4chgzS0+bPBcngcALxkWQ4snGd4LmuzjRJJBczT0eu4vBhIgq8uXgZP175F6xOPTB5W2WgMDxriJAZdM1mCJdrW61hysZ1gue6OkML8Oga6t85dUQUcblc3dEIctWzBonqOa/mM17HMYRZs5/A8NC2tg93oyFCCrrNm5vvokpzhMrS6QHNyDEb1+WgYIdWgjlwmIJW1O2zzz647777MHv2bHX9jDPOwFe/+lU1r07IIOIFCxYgk8m492yJyH1Swa1aBVT/LefzwPp1QDYH9xskWpzzand4sBsNEbWhwrvtXh0q3ChfaY6oy9K5iRk57wYFO4UDhylIRd2TTz6Jovx3teq//uu/8LnPfW6yqDMMA9lscEYdEFGXZIVOCjqpYhIJlLLAuA4MjgDDLjdIWJ3zKlm6XXYuI5UowiinPW2AaCQFXbLZpy+GICPHbFwXg4KdykjGmKej4G6/ShHXyHy6BBGFnOkUhLJsLaaBAQ/rqceu+DUO/eybVHPEw195QOXnpv/UoY4ycszGdTEo2MHVOiKXMVNHRIEkObra+BJd82fGV605wrRJ0fyOXmbnmJEjol6KOlmFa1yJ48ocUQQbI0wxCsnT5d1rKG3pXZ/db/LtRPJp10+EaNccEU80aYSQUyPGxoCZM50fMszsnE8NEj590xN5VdTJdus73/lODAxUHjIxMYGjjjoK8epPPHPejojC2xihZDLIp0axdn2lxpPLsHf1VJ03LhzDUHobyqUR106EsNMcIQWdZZOE3FG2qffc0/khw+2yc8zIudcgwcYGinJR95WvfKXu+gc+8IFp9/nQhz7kzLMiIt8aI5TRUZTKcVXMSSEjBZ1lQePB2a7bzdym8nTl8kDwmiMa7+h1do4ZORcbJIbY2ED9U9QRUXQbI5TsVK3iSUFnGHWz6ZoNGPaKaUxfa27sUtRydMzO+dgg4V+XNVG32ChB1O85umxWZedkdImZp1m6JvPpND3r6ikRreq08XFgeEabgcJu5OmYo/OP4U9DDpFTvJx/3tTVV1+tBhcnk0kcfPDBePDBB2097uabb1bNGkcffbTrz5Eokjm6e+9F/r4HseHxDB59TMeaRzB5eeKJSpbO6ey/3XNeFywYx3bbPYzU0BPVJgnd8zzdzvPbrFK6kaezytExO+dNnq68iVk6CjXfV+puueUWnHPOObj22mtVQbdixQoceeSRWLduHebOndv0cX/5y1/U8OO3vOUtnj5foqjl6EqDCYylK0OFGwsYP7J0cs7rq3gd1l7yvcqwYZdOibAjZvcnpFsvUi1Hx+ycd3m65GJuvVJo+b5Sd+WVV+LUU0/FiSeeqI4ik+IunU7jxhtvbPqYUqmET3ziE7jwwgux2267efp8iaKYozMGKwWdXDVf/GiOkPl0gATp4ur0CHWChAcFnWy5qhEueZsxOdl6lYubpKDTpajjYHfPaF6d9UYUgKLue9/7HnKyddMgn8+r93VCHrN69WocccQRU09I19X1+++/v+nj5MxZWcU7+eSTO3z2RBQEmpYHtCwm8nl1yfo1DK9hJl1tfItcbzqXTkgxt3Fj5U/Zfu11j1pydOWJysXcHEEubLFmrS+cTUf9uP0qK2rvec97pm2Nbt26Vb3v+OOPt/2xNm7cqFbdtt9++7rb5bqcNWvlvvvuw3/+53/i0UcftfU5pAA1F6FjEmom6mfy78HUHJHXvC/oEqm1+PCX3o/VT+6oNlcfwRvhp8aZdKLpXLraA6SYW7QISKV6y9OxMcL7GXTNcDYd9VtRJ0OIrU6S+Otf/4rR0VG4SQrH4447Dtdffz3mzJlj6zGXXXaZ2qYloqkGifzmjGqOkCyd7Gx6OlxYKyNXLFYLOiCNDA5A5T9pa7C/Out16d6vIJWwO1PEh5l0NbJH3WuDRLMBw2yOcHEGnWl0T50Y83TUH0XdAQccMHlUmPlkCSGrbc8884xaweuEFGaxWAwvv/xy3e1yfd68edPu//TTT6sGCTnJoqYs/2OWL2RgQDVX7L777nWPOe+881Qjhnmlbv78+R09T6KoNUiUBuN1zRF+DRcW917zI+D0ytvJ//w01iZvVgVdX8bIzAOG2Rzh4gy6Fqt1RP1Q1NXGhsi2p3SnDpv+Wy9HhclIkk5PlJDHLV26FHfdddfkx5ciTa6feeaZ0+6/9957449//GPdbeeff75awfvmN79pWawlEgl1ISLzP75EXXOE62SFP5ednD2HgRzSGFfXh4zKnyKVLMFIer9C5y/p721ojCAHz3Gt/MefmTnqBx2fKCHF27HHHqtmyjlBVtGWL1+OAw88EMuWLVMjTcbHx1U+T0hGb6eddlLbqPI5991337rHz5Shn8C024koGDTksPNFpyG9of4/ZOP458obZ6B/qTzdKX4/i/7J0DEzRxHXcaZOCjAnSYH46quv4oILLsBLL72E/fffH3fcccdk88Rzzz2nOmKJyLnTI/LVxQsvmiJS+iPTCjor2b0Xwui3VXXJ0+XXV96OL2yR9SJnMnTMzFG02SrqZs+ejfXr16sM3KxZsywbJWo2bdrU8ZOQrVar7VZxzz33tHzsd77znY4/H1Ffnh6RyaiCrtYgkS150ByhlaHHps4fe+7GK2EkExjPxnHQSR9Ttz30nz9EWrZdpaDryyBd1bwb+vvrdwszdNRHbBV13/jGNzBjxozJt1sVdUQUjtMj5vgwXLgUn4GMMYQMBpBRQ4YBQ4Yf+5ijk5l0tocNVxuznBs6bMrTqYHLREQuF3XmLdcTTjihh09HRL5RW5tyekSlmPOj2/W4C9+J32/YFUFRGzosU0lsDRs2Z4l7HTrMPJ3zzRB1t/s70Joo0Jk6u0N7R0ZGenk+RBSZEyPK0LT6X6xrNtQPLfdrHp3V0GGZT2dr2HAt9ycFXU9Dh5mnc32gMBsjqM/YLuqky7TVtmttKLHMrCMinxohrHh8eoQUdMn02sksnWZxrOCDN96KdKIYmHl0tlcuq2flOo55OpcGCrMxgvqL7aLu7rvvrivg3ve+9+GGG25Q40aIKBiNEFbMzRGenB5RbY4wSgNA1kApO/2TSUEnzRF+apmlM+fnatd7IVutsjJXd5v5jNcAVLZhxmYIos6KusMPP7zuupwE8aY3vQm77bab3Q9BRC43QlhNEpYVOs9PjzAMbH/eFUiuexpB1DJLZ5Wf6yVDx7NdncWBwkTOzakjooBqsTVYa47wagyclstPK+juw2G+nutqO0tnlZ/rJUPX7GzXGp7x2sFryYHCRK2wqCOijhsgWt+nPkP35LeuwdIzjlMF3YM33obtRnKByNG1zdK5kZ8zn+1awzNe7eNAYSL3ijrOqyMKxikRZqopwlR35XPuNEBYn++aV2/qua2TN49rQ5Mz6SRLF5SCzjFWeblm2Tme7eoMZuiIeivqPvjBD9Zdz2azOO200zA0VPlhXfOjH/3I7ockIqeaI+TP0VFV0K1dO63Oc6Y5otYAUR6EIfu5ZoaBef/yNcsM3VtPr//ZESnMyzk7W67t4zh7jsiRom50dLTu+ic/+Um7DyUit5sj5N+nnBiRrRRwjduKTjZHSEFnGPXhPC2btSzoajk6EYQsnePa5eXMmJ2zN1uuHc6eI+q9qPv2t79t965E5JUmuS8vmyIaPXvjSnWOayYbw1tOlhNoNDWXLkhZOldY5eXMmJ2zOVuuHc6eI2qGjRJEZKsxorEBoplxIw0DSWQ0+fGiBTJLNzmnWc2fa9gGtDuTTmXpmJfrGnNxRI5jUUcUoeYIp5oimjVGaFoWBoYt65uag07+yGRjRFBn1I2PAzNSBcRe2wgMW6wWtZtJxyxd99k55uKIXMOijihizRGOnRhh0RihCrrGJgkAE7lY0w8TtCxdbUbd/J3KGITFPDo7M+kas3TMy3WWnWMujsgVLOqIItYc4fSJEVaNEa3ce82PkByZ+tESlPNdG8XkKRYdmEcnWTp9lixjOvn0Ip6dYy6OyA0s6ojCJiDNEUbZQHascnhqbuvUSlwqXkQqGcwCp+6816YHv9qYSdeYpWNB1xyzc0SeYVFHRF0VdFuPuwpLsg+F5tUzn/eaiBUQ2zYGzJlp/zxX5ug6Y9hsOCEix7CoIwpDY4TwoDnC3CTRqttVVuisCrpHU8swYtp6DZK6815RxmAhCey5p/3zXK1m0jFL1+S1ygPlTYA+p7LVSkSeCOZPXyKa3hjhRXNEQ9drs25XszUrrkdiRuUXtxR0mh7MrdcalT2UN2QhqdvwYW0mHWfPtc7TJRcDuk8DE4n6EIs6orA0RggvmiNMXa8lw6Lb1TCgS8FZJQVdaqaDnRlOKhZQzJfVzmm5COgSo5PjabVC7zk6vcWgYaq+TjZXQYnIESzqiIKuRXemm80RaoyJRUG3w5cuQnLdBgResQD8fSO2ZZOIV7+MZByISW02aGMWnWCOrruZdJxFR+QLFnVEZJuWy9UVdHK26yxVVXZxOLvbymWUB5Io7LEIuy1OIJ4AYjowGLc5i04wR9f9TDrOoiPyHIs6IurqOLC5eBmv4nVYq90c6FewNJBAfDSJVA+j6BTm6DqcScdZdEReY1FHFKQOV7NsttLhWt/0qvhyHJhhqJW6mnF1FFiwmyImqVBdtovHMUfXEc6kI/IVizqioHW4VklBt+HxDMbSo1Ync3l7HFiYsnTTGIhvOhkoNIwjIedydMzQEQUCizqioHW41t6VhSroBtJxy85WL48Ds8rSZZAO3LmuVnQti1ivBR3n0bXP0TFDR+Q7FnVEAe5wlcUyr4//spule/DG27DdSC5cJ2TVcnGd4jw6Gzk6ZuiI/MaijqgPmh9a3y/XMjun7mORpUsnioEr6IpZ2Q6sZOe0YgGlckGt1E3ifDnnMUdHFBgs6oj8bI5oOPrLqeO/rJofWt+/2hgR4uycFHQ77vJVDG33rN9PJdqYoyMKLBZ1RH43R1SP/hKOHf9l0fzQSq0xQstlWxZ0gc7SGdnWBR1zcQ68xszREQUZizoiv5sjqkd/qXeVnT3+q1XzQzvP3rgSRjXMl8nGcNDJH1EFXSiydFbZOebiesccHVGgsagj6rPjv6YxZejM2Tkp6IxkUo14G88lkFFZOviUpTOmZf+KRaBsVN7WikUYhW1T72R2zl3M0REFEos6on7WJkMnBd1HvvRurF431/OnZnoW2GHBRUimw5fziwzm6IhCgUUdkdenRjQ0R9QaI9TbOWe6W+0e9dU4f64mu/dCtVI3kYvVFXR+ZOnka+mooGN2zlnM0RGFBos6Ij9Ojag2RzQ2RthtjrDT3Wp51JfNDJ36s2GP9cEbb/U9S/fsupUwygnkC0AuC+y5EEgiD0xkgCX7V7azmZ1zFnN0RKHBoo7Ij1Mjqs0RcmqEuTFC2GqOsNHdOu2orzZqGbrJ6waQyU39iAjCXLpcLgGjnFR5OkMqu7wMwNXlyQN6CtCtM4rkAOboiAKPRR1RiBsjeulubf1xg5ClU88E8xZcPHlNVufKUscVCxjKbkRsIglI3Sqvq677+kxDlY3r6DE9DEwkIk+xqCOiaYKQpZvM06Uqs+dKhV2w58JqAZsvq4JucP9FlWpYCrrqWBiymY3r6C9CHhPjS0sUcCzqiNxQltWNahFUzgDqqCrrrVCpR7TqxS67jRBOCEKWThTGz0cyYXoSxdarn9QuG9cJnutKFAYs6oic7GytbVflHgLKE5XrhRwK+nqUyimgWL+aVDKAGbMqtUm8wyHDnTZCWM2mM8+lkwHDRnU1xp8s3fRZdJpuul6QKq56vVDw4glFE7NxRJHFoo7I6c5WWZVLrAUM+ec1qOqP55/SMJ4cgjEwvXIr5IHEIFDucHez00aIdrPpKidGVAYMB3EWXWzz34GhkakbmKPrLCPHbBxR5LGoI3K6s1W2WQdSgCHHVCVQigHb4jpi8UEMWPyLSyWBWKzSnOAFq9l0tTNdG3mVpWs3i66U2wODQzOARdUMnej3HF03GTlm44gijUUdkRMas126FBtSfFTnvg1AFXRBq0Ge/NY1WHrGcZNnuspWq5kUdF5n6Wqz6MTkPLoFkuqaYIau54wcs3FEUcaijqifmM55FePaUN2ZrulkyffsnBR0hlGZRSdb02oeXSEH+NyoEVjMyBFRFYs6ol6aItRREHl5w3wnlemv3S0wmX6LLN1bT/9gIM9xlYJu82YgoVfn0RWTwIw+nkVnlZ1jRo6IGrCoI+qlKUIKOuNxIJ4G9ErTghR0zz2fxfjW4cmcnBwHlu5iPJibWTpzji4o2blsZqEapiwDhmWretfXl5EsVefRpVLB27/2OzvHjBwRmbCoI+qpKSJbKegG5Bfu4OTdpKCLxaYaI6Sgs2qS8MtcvIxX8brJHJ3f2bmayukYU09kMA4M1ubR9WNB1zY7x4wcEU0J0K8ZorA2RUgxJ5epAkVW6ALTGGExk25c5eg0D3J09fk5q+yclXKuAL1QruxdM0tXffFk+9nn5V4iCjQWdURR1mQmXZDyc42K2QIKL21EalYSMVml69csXS1Hx+wcEdnEoo6oI9ObIoKs1Uw6t3N0rfJzteycFaNURmwoiXnvXITB7fr0XNfGHB2zc0RkA4s6IrudrhNbUNDXoVTW6477Usd1laeO6wpMt2ubmXRenufamJ9rzM5ZGRjq43Ndp+XomJ0jovZY1BHZ7HTNF7L420tjGB+ci5I522RIl2v9cV2edrs2zJ4zazWTzr2CrpKjs5Ofk9El0uk6+XyLBZSyAa2K/cAcHRF1gEUdkc1O11JuEJnELGgDaSR0663DGs+6XTvIzHkzk85+jq42i662syoF3eDWjSgPJjEwnERsQO/fWXTM0RFRF1jUEXXQ6WoMDGIwBgzWL8wFKjNnxauZdFY5umb5udosut12r76e+TK0TBLlNyxCbCiF+HC8v2fRMUdHRB1iUUcUEc/euBLG5Py8ikw2hoNO/kjd2a5ezaSr5eja5eekoEvWnrZ0u45KId1HBV3TWXTM0RFRZ1jUEdntdG04rzRos+ekoCsnkpjIxSZvy2gDDp7tOv3M1kZ259BR40tbzREyQ0dEPWBRR/3d0drs5oktiGnrYJg6XUtlma4hRd5Up6svmuToZODxR770bqxeNzcwM+fIzkubB8qbAH1OZXWOiCjMRd3VV1+Nyy+/HC+99BKWLFmCq666CsuWLbO87/XXX4/vfe97WLt2rbq+dOlSXHrppU3vT9T67NZKp+r6dUDWtAilx7IYGR1DPl3f6ZrPDWN0dDBwObrs3gvVFmuzgq7XHF27M1s7mUNHTbZek4uBNg04RESBLupuueUWnHPOObj22mtx8MEHY8WKFTjyyCOxbt06zJ07/RfUPffcg4997GM49NBDkUwm8W//9m9497vfjSeeeAI77bSTL18DhfnsVqCUBcZ1YHBkqmNV0wehz5qFQaQxYCpOUslgneFay9GpLF1uKrf24I23qu3WGidzdFZntjZqlaOrjTGRP6e9o59pfZYjJCLH+T4z4Morr8Spp56KE088Efvss48q7tLpNG688UbL+//gBz/AGWecgf333x977703brjhBpTLZdx1112eP3cKeUer6VKOJ9UpBvGRqctgalDVf+aLrwWd5OiyWcscXSY3oC41tfxc7dJdQScZumzlYpGVa3VpVdDJGJNstvJ2PAHE9OrE5rGxyt9Hvx0JVsvTERH1yNc1h3w+j9WrV+O8886bvE3XdRxxxBG4//77bX2MTCaDQqGA2bNnW74/l8upS82Y/OIgChvPc3TuZOgax5hIQafGmeTKlYJuzz3760gw5umIyEG+/pd448aNKJVK2H777etul+uSr7Pj3HPPxY477qgKQSuXXXYZRkdHJy/z58935LlT2Dta6y/SAKGZL0HqdO0wR+fEHLpmGTqnsnK1MSbT5v0FZQCgV5inIyIHBSgd1Lmvfe1ruPnmm1XOTvJ1VmQVUDJ75pU6FnZ92uXa5OzWkgHMmFXZlY2bagp1pqvfna5d5OicnkNnztDZObO1mWK2gHKuDL1Yra3NgnpgrleYpyOisBd1c+bMQSwWw8svv1x3u1yfN29ey8f++7//uyrqfv3rX2Px4sVN75dIJNSF+kizLtdmZ7dKTZEHEoNA2VQHqoKu4UzXIFAFncV/YnqfQ9fk8zkwb04KuonnN6rcoswVjk1UBw2bRT1PZz4GbPK2YK0IE1G4+VrUxeNxNZJEmhyOPvpodVut6eHMM89s+rivf/3ruOSSS/DLX/4SBx54oIfPmELd5dri7Fbpao3FKhm1oA0YFubmiOq71JBhc3OEg5+4rjHCkY9YKquCbv57FiE5ksCgVWxOCrqo5umsjgGr4XFgRBSV7VfZGl2+fLkqzmTWnIw0GR8fV92w4vjjj1ejSiQbJ2SEyQUXXICbbroJCxYsmMzeDQ8PqwtRs3Nbg3h2a6eNEaEbMlwsoJgvo5StbK8ODidUZ3HfsTwGrIbHgRFRRIq6Y489Fq+++qoq1KRAk1Eld9xxx2TzxHPPPac6YmuuueYa1TV7zDHH1H2cr3zlK/jXf/1Xz58/kVeNEW43R1g1SPTUGFEsAH/fiG3ZpMoqDgwnERvQ+2drte791ZVPHgNGRFEu6oRstTbbbpUmCLO//OUvHj0rCryy/KK0KGTKGUCTrtaG5biAdbR22hghvGqOUJ933UqUSyNdN0bI/JLyQBKFPRZht8UJpIZ0xIfj/bW1asZtViLqh6KOqJNzWhUjh/zWh1AuSeK+QSEHTV8Po5yq73INytmtTXJzZubbxo00DFS37HKwHDLc45OYHOHSOGS464LOpDSQQHw0iXiyH7dWzbjNSkTuYlFHoepgNXeyvrBxLXK5AZTL0wNyekFDYcaQytDVPS4AZ7fayc2ZHXTyR5DBkFtPwtEMnXS5SlOE0IoFlPrp5C9urRKRz1jUUag6WCffnRtEdlMKWiqFAatD0DUd8YaCLmhntzbLzZndh8NUhs5K0IYMm8eW1Mjxa8m0XjkKLLJZunBu6xNR9ATk1xuR/Q7WmvJgHIOxBAYHwz+H0JybE5lsrLpCl8aDN96mtlkbBW3IsHlsiXS5KrqOWDIezUkl5iwd83JEFAAs6oh8ytGZc3O1gcKT8+e0gckt1+5yc1M5uVYaM3S9DhkWUtClZkY1QNcsSzcEWK0YExF5iEUdUeTmzzk8a85mhq42iy7yauNL6saUsKAjIv+xqCNf5MZyKOVbrD5ls9DzeRg56VadrpAPZ47JKkcns+dkpU5W6JyYP9csJ9eKUxm6SM+isxpfwm1XIgoQFnXkS0H3p++uQmGLdWer0PQ8Zgw/jtLM9LQO1ppSKYvEYIDGk3SZo1NZuoZwnFPz58w5uVacytBJQRfJWXRNx5dwTAkRBQeLOvKcrNBJQRdLxRFLWRccmp7FwKw0dD0Nw2JkiZCCLhYLyHiSLs5wlWKunEiqFTrR2fy55pk5N3JyrfRNhs6M40uIKIBY1JFvpKBrdg6opgODqUGUS4PdH1MV2Qydd5k5apKjIyIKIBZ1RAE6w9VOjs5uZq6nc1tpCnN0RBQSLOqIAnSGq+gkR9cqM9dNTo6sXkjm6IgoHFjUkS+kEUJyc7LNavl+GzPWgnt2q6G+PqEVckCqcmtGi8HQK/k5FIqQCSDp1Li6OpTahrTNLlevM3PNOl9DN8KktoXa8ePMo0usT/cgIgoCFnXk+biS3OYtGJq1DvFRHYNDzTslNS0LA8Phy8v9AMD+puuPVP54A06f9iHGHzkeYWMeZRKaESaNW6id4ugSIgoBFnXk+bgSPZZFKjUG6HNRLjX/JasKOmMwXHk5WZV7ozef36/MnHmUSXJmKhwjTKZtoXaKo0uIKPhY1JEP40oGkdpuFrSBdGSC/LW8nGyN7oIz1W3PPX4FxjNpvPWMD6qmiHuv+TFSVme4xrubRedHZs687SqjTEJR0JlxC5WIIoxFHfkyriSWkHElIXvxG/Jzk2+nqhepsUw7kdvKszBuDGPjRKXDNTk4gFTcqgiTsS0IvNBsu1pl5ziKhIj6AIs6ol7zc2/E5Oqc2UEnfQSZiaHIvL6h2HZtlZ1jLo6IIo5FHZEL+bn7Vh+GzES6pzNcfVUsAOWp1a5iEeHYdm2ZnWMujoiijUUdOS4y40o6yM/V5sVlcrHqCl0aD954W1ez53xXLGBg80YY8eRkQbdtKxAfhL/brnZGknD8CBH1MRZ1ERsX4rcojStplp+rzZvTazPnZDTJhJxRm1Tnt9a2XNuf3xpQ5bIq6CZ2X4TyYAKFPJAdB3bbB0gN6f6s0nUykoTbrETUp1jURWxciN8iMa6kTX7Oat5c1PJzQgo6JJKQtTGjCMRHgeriXcBHknCblYj6E4u6iI0L8VsUxpX0mp8Le4ZOk7erCoXKJTA4koSIqCkWdREbF+K30I4raZOfm8gXJ1fo1jx0PRKDU1uvs7QE1v7g5rrHhTlDJ0qJJIolHRs3Aclk5RLzc4KJEaTKkogomFjUUf+wPKdV3YyJ3FSRpueK0/Jzpnergi4lXQOTujhPNMAZOkXXUSrHVTG3aBGQSgFxvxpeJU9X3gTocypbq0REZIlFHfWHZjm5Zlrk5yKjuuUq263S4ZozEiijulpXntp2jSd8LOjMebrkYkAP55Y+EZEXWNRRX7DMyTXTIj/36J+WYWQwAv9sTFuuUtC9lk0iP6Grhggzz7ddW50GoQV0Nh4RUUBE4LdT/4wfyW3J9vwxNC0PaO5tF4ZhBl0tJycy2RgOOvkj6u17r/kRUvGimj9XW6F78vFrUC5PrQ5JQaeFJixnb8tVVuikoNtn/7halTOTgs6zVTqeBkFE1BMWdSEbP1IcyyA5d7Trgi6ZXqvGjrjJ9xl0Ftk5qzlzYiI2UFmZk3fNLCItDQ76VCEu57UaQR290q1ioW7LNacl1XQZKehSfvbf8DQIIqKesKgL2fgRKegG0l0unWhlVdAZ5UFXCxVfZ9C1ys41ycmNP3I8+kZ127WgJ6e2XAeDts2aBHQbQ4aJiKgOi7o+HD8iBV1YZ8h1nZ1rkZNrJptZGL3XqbrtOj5/EfLF1OSWK7dZiYjCj0UdRZY5O9dqzpxIxafPlasUdBHIzzXpdDUG4/5suXKblYjIFSzqKFQz5abebSBbmN6Eostjqhk5c3au9Zw5IauXiCaLTtdCSUcy7dGWa+NWK7dZiYhcwaKOwjlTTvJxb2ryvkcqf0R6xpwDna6eDBRu1tGqrnOQMBGRk1jUUfhmynWRj4vcnLkuyGkRMlxYZtF5NlC46VZrjIOEiYgc1p+/3Tyi6XloeladhxoEYZgh1yoXV2POxz32+5WIx6av+BiJuHzB026PzJy5Jpk5K7UcXSEPFLz60mtbrtxqJSLyDIs6txg5DM1ah/iojsGh4EzC932GnFq8MZAdazi6wMTq7FUzcz5ucCiN5LR8XJ8xZeYs320+McKL8SWNW67caiUi8gSLOteUoA/kYBgzUC4FZ+aWrzPkqgXd1uOuwpLsQ63v2A9nr7qQmZMt1kayQmc+McL18SXTtly51UpE5AUWdS6L8ky4bsgKXduCzmZmrp/zcY1jSoQq6BJJFBp2YmXL1bETI6yGBk+7D4cIExH5oY9/I5Lf1qy4HokZ0/Nwrc5e7Yt8XBdbrqWEnMKgq4Ju48bKFquZI1uurc5mbcQtVyIiz7GoI0/zcrmt1dlyKSA5WkTK4hjbyJ+96saWq64Dg3GUc5UCbtGiyspcjSNbri2HBjfilisRkddY1JFjmJfzecs1B/WncOykCPN2K7dViYgCjUUdBTIvF9mzVz3YcnWsu9Vqu5XbqkREgcWijnzLyz27biWMJnm5yJ296uGWq2PdrZbbrdxWJSIKKhZ1Ecqr+a3TvJwUdIbh9Wny4R0q7NmW6+S2K7tYiYjChEVd1PJqfuN8OVeHCru+5dq47crtViKi0GBRF6W8mt+Yl3N/qLDbW67Ttl2HeEYrEVFIsKiLUF7Nb8zLOXtua+3M1pyRQBnVJTm5m6ddrmxUISIKCxZ1rubfDOiFEmK56m9gJ/JqgCroUjO7md1mQNOqv7BdwLycs+e21p3ZahGlZJcrERGZsahzqaB7eo8PYMm2+xEcBnZYcBGS6Q1+PxFqmDfXuMUq7zIMIJ8HCqWpM1sbscuViIjMWNS5ILMxg8UuFXSPppZhZKTzvzZZofOqoON8uc7nzZUTKZWVUw0QW6sNEDqQlGhbqsu8nJ1zWtX92OVKRBQFLOpc9vDXViD9uhHHPp4UdJre2+y2dvPhesX5ct3Nm6u+y5kGiE7OaRXsciUiCj0WdW5S89oKSM+aysP1rtR1E0MN58MFc96ccKwBoqNzWgWHChMRhR2LOjdIIKo6r21ffN6VT0HRnDfXUQNEq+1VbqkSEfUdFnWuyNg+39RrzLsFe96csLXlamd7lVuqRER9hUWdC9LbpYHXKm8/veYyxJJzEBTMu3k0Z04rQFND5QCtXESpXEDeMKq3VMNzuSxKBUAiklLEpcyFXLv+BlmJa7u9yi1VIqJ+wqLOBeZGBgM837Tv5sxpBQwkXkO52oxSKgJbSwnkc+MwytNnFg7Jdqu8YaNRtY46xosnPhARUQWLOqIeVuWs5sxpyMIobcO23FJ1EkQuB+SKMSx9YwJJq3lzsamt185wJY6IiKY4dQR4T66++mosWLAAyWQSBx98MB588MGW97/11lux9957q/vvt99+uP322z17rtTHqqtysWxm8qIXCyglYzASGrQE1KWka/j7tgReK8zGa/k5yGhzkBidhaHhNNIWl0QqDejdXHiEFxERBWil7pZbbsE555yDa6+9VhV0K1aswJFHHol169Zh7ty50+7/hz/8AR/72Mdw2WWX4R/+4R9w00034eijj8YjjzyCfffdF4GgpZGZ/Rf89f//n9CSTp2yTkFseNCQx2DsFcQGZCZJZS5JyQAG4mkcuDiGZGpqNS7BGoyIiFykGYbM3/CPFHIHHXQQVq5cqa6Xy2XMnz8fZ511Fr74xS9Ou/+xxx6L8fFx/OxnP5u87U1vehP2339/VRi2MzY2htHRUWzZsgUjI84NBW6UeXUjnv/JtdDTszGQcu/zkLdbreVcBuN7LkG5mqHTtRx05LGttAxlo3JbLg/kcjG85a0JpG3O/iUiIuq1dvF1+zWfz2P16tU44ogjpp6Qrqvr999vfcyW3G6+v5CVvWb3z+Vy6sUwX4i62WqVgi6jZZDLjSOf3aou2Yk8tmbS2Dw2hC1b0+qSzaWRSifU6hwREVFfbL9u3LgRpVIJ22+/fd3tcv3JJ5+0fMxLL71keX+53Yps01544YUOPmuK3PgRublYmRldtyqnJZGZvwhGdau1WMiiXB7H7ge8GQnztGAtVgnTmXC7lYiI+i5T57bzzjtPZfZqZKVOtnepD1mNH6kWdNu2AvHB+rsbiRgyJQ1GdT1b1zWk0nGkh5KV5gYiIqIA8bWomzNnDmKxGF5++eW62+X6vHnzLB8jt3dy/0QioS7UZyxW5GT1raAnMW5afRP5fOXM1d32MY0WMfLQ8AowONUAIfSBNOIJ7qsSEVHw+Jqpi8fjWLp0Ke66667J26RRQq4fcsghlo+R2833F7/61a+a3p/6kEUeTmXisgVszlVW3zIlTF6KMSA5CqTkMlK7GEjOGEVy1jIkZ7918hIfOZijRIiIKJB8336VrdHly5fjwAMPxLJly9RIE+luPfHEE9X7jz/+eOy0004qGyfOPvtsHH744bjiiivw/ve/HzfffDMefvhhXHfddT5/JeRV/q1ZDq5VHk4U8nno2it4wz45xNUK3BRdBgDLvwbzp+OJDUREFCK+F3UyouTVV1/FBRdcoJodZDTJHXfcMdkM8dxzz6mO2JpDDz1UzaY7//zz8S//8i/Yc8898ZOf/CQ4M+rI1fzb5LuLwPh4AfFB64k85XiiLg8n9JiBeHIUidFl9Y0OTfHEBiIiCg/f59R5jXPqgqGYLcAoNV+Fk9U2vZTBxG57wxiYnoksSgiu/Bp23zM9rcFBiWkyAXjazSoTxy1UIiKKYO3i+0od9V/BVs4Xkd84hpicZN+EphUQG8kgWyrCMK3UTooBg+mdkJi7GImOTu3g6hsREUUTizpyvKCbeH5jy4JNJOcNYYe37IaBlNUyW7X7VC+gPPRmQLP+WLGBGBJJdjYTEREJFnUuK00UYJTrQ/lRVsoWMDAcw+vfvTsGh5sUXEYBA7EtGBySK4UmH0kDtJlAaojdpkRERDawqHNJfmsOsaEEimMTKOaL6BeaXkRquwKSo0XE0y0m5mg7AMnFgNZq65RbpURERHaxqHNJLJHAQGom5iyZgeTMFPqJHh9GfPZSFmxEREQeYlHnkhk7zcHg2z6C5GiTzFikcYWNiIjIayzqXJScNermhyciIiIKxjFhREREROQMFnVEREREEcCijoiIiCgCWNQRERERRQCLOiIiIqIIYFFHREREFAEs6oiIiIgigEUdERERUQSwqCMiIiKKABZ1RERERBHAoo6IiIgoAljUEREREUUAizoiIiKiCGBRR0RERBQBA+gzhmGoP8fGxvx+KkRERERt1WqWWg3TTN8VdVu3blV/zp8/3++nQkRERNRRDTM6Otr0/ZrRruyLmHK5jL/97W+YMWMGNE2bvP2ggw7CQw891PKxvd5HKm0pJp9//nmMjIwgzOy8FmH5vL1+zG4f38nj7N633f365fvTr+9Rtz6nH9+jfnx/9tP3KH+G9v5a9NPPUMMwVEG34447QtebJ+f6bqVOXozXv/71026PxWJt/xKcuo+8P+w/kOx8nWH5vL1+zG4f38nj7N633f365fvTr+9Rtz6nH9+jfnx/9tP3KH+G9v5a9NvP0NEWK3Q1bJSo+tSnPuXZfaLAr6/Tjc/b68fs9vGdPM7ufdvdr1++P/36Wt36nH58j/rx/dnp5w0z/gzt/bXgz9Dp+m771U+yNCuV9pYtW0L/v0yKHn5/UtDxe5SCbCwAv+O5UuehRCKBr3zlK+pPoqDh9ycFHb9HKcgSAfgdz5U6IiIiogjgSh0RERFRBLCoIyIiIooAFnVEREREEcCijoiIiCgCWNQF0ObNm3HggQdi//33x7777ovrr7/e76dEVEcmpr/tbW/DPvvsg8WLF+PWW2/lK0SB8o//+I+YNWsWjjnmGL+fChF+9rOfYa+99sKee+6JG264wbVXhN2vAVQqlZDL5ZBOpzE+Pq4Ku4cffhjbbbed30+NSHnxxRfx8ssvq/94vPTSS1i6dCnWr1+PoaEhvkIUCPfcc486Vum73/0ubrvtNr+fDvWxYrGo/gN89913qzl28vPyD3/4gyu/07lSF0ByDIkUdEKKO5kPzRnRFCQ77LCDKujEvHnzMGfOHGzatMnvp0U0SVaS5YxvIr89+OCDWLRoEXbaaScMDw/jve99L+68805XPheLui7ce++9OOqoo9TBupqm4Sc/+cm0+1x99dVYsGABkskkDj74YPWX2ukW7JIlS9Q5tZ///OfVL02iIH2P1qxevVqtLstB1kRB+/4k8vv79W9/+5sq6Grk7RdeeAFuYFHXBdkSlYJL/hKt3HLLLTjnnHPUZOlHHnlE3ffII4/EK6+8MnmfWl6u8SJ/+WLmzJl47LHH8Mwzz+Cmm25SW11EQfoeFbI6d/zxx+O6667jXw4F7vuTKCjfr56Rs1+pe/IS/vjHP667bdmyZcanPvWpyeulUsnYcccdjcsuu6yrz3H66acbt956K/+aKFDfo9ls1njLW95ifO973+PfDAXyZ+jdd99tfOhDH+LfDvn6/fr73//eOProoyfff/bZZxs/+MEPDDdwpc5h+XxebUcdccQRk7fpuq6u33///bY+hqzKScBXyMHAsvQrXTNEQfkelZ9tJ5xwAt7xjnfguOOO418MBer7kyhI36/Lli3D2rVr1Zbrtm3b8Itf/EKt5LlhwJWP2sc2btyo8kXbb7993e1y/cknn7T1MZ599ln80z/902SDxFlnnYX99tvPpWdM/caJ79Hf//73astBxpnU8iXf//73+X1Kgfj+FPJLVSIssnUm2WQZu3PIIYfwb4g8/34dGBjAFVdcgbe//e0ol8v4whe+4No0CxZ1ASRV/aOPPur30yBq6s1vfrP64UQUVL/+9a/9fgpEk/7P//k/6uI2br86TLpUZSRJY2ODXJfRD0R+4/coBRm/PylM5gTsdz6LOofF43E1WPCuu+6avE1WNOQ6l/4pCPg9SkHG708Kk3jAfudz+7ULEnR86qmnJq/L2BHZLp09ezZ23nln1dq8fPlyddSXbKWuWLFC5TpOPPFEJ//uiPg9SqHEn6EUJtvC9DvflZ7aiJM2eXnpGi/Lly+fvM9VV11l7LzzzkY8Hlftzg888ICvz5n6C79HKcj4/UlhcneIfufz7FciIiKiCGCmjoiIiCgCWNQRERERRQCLOiIiIqIIYFFHREREFAEs6oiIiIgigEUdERERUQSwqCMiIiKKABZ1RERERBHAoo6IiIgoAljUEREBuOeee6BpGjZv3tz3r4e8Dj/5yU/6/nUgChsWdUTkqxNOOEEVEY0X8wHaTnvb296Gz3zmM3W3HXrooXjxxRcxOjrq2ue1+jrNl3/913/t6WOzECPqbwN+PwEiove85z349re/XfdCvO51r5v2wuTzecTjcVdeMPm48+bNc/UvQ4rGmltuuQUXXHAB1q1bN3nb8PCwq5+fiKKNK3VE5LtEIqEKKvMlFoupFbUzzzxTrarNmTMHRx55pLr/lVdeif322w9DQ0OYP38+zjjjDGzbtq3uY/7+979Xj0+n05g1a5Z67GuvvaZWBn/729/im9/85uQK2V/+8hfL7df/9//+HxYtWqSe34IFC3DFFVfUfQ657dJLL8VJJ52EGTNmYOedd8Z1113X9Os0f32yIiifz3zbzTffjDe84Q1IJpPYe++98a1vfauuoJXXYocddlDv32WXXXDZZZdNPg/xj//4j+pj1q6L//mf/8Eb3/hG9ZjddtsNF154IYrF4uT7N2zYgLe+9a3q/fvssw9+9atf9fA3SUR+YlFHRIH23e9+V62iSZF27bXXqtt0Xcd//Md/4IknnlDv/81vfoMvfOELk4959NFH8c53vlMVKffffz/uu+8+HHXUUSiVSqqYO+SQQ3DqqaeqlTO5SGHYaPXq1fjIRz6Cj370o/jjH/+otka//OUv4zvf+U7d/aTQO/DAA7FmzRpVXJ5++ul1q292/eAHP1Ard5dccgn+9Kc/qWJRPp98fUK+3p/+9Kf44Q9/qD6+3L9WvD300EPqT1ntlK+ndv13v/sdjj/+eJx99tn43//9X/zf//t/1fOXzyHK5TI++MEPqtd31apV6vU999xzO37uRBQQBhGRj5YvX27EYjFjaGho8nLMMceo9x1++OHGAQcc0PZj3HrrrcZ22203ef1jH/uYcdhhhzW9v3zcs88+u+62u+++25Afia+99pq6/vGPf9x417veVXefz3/+88Y+++wzeX2XXXYxPvnJT05eL5fLxty5c41rrrmm7XP+9re/bYyOjk5e33333Y2bbrqp7j4XXXSRccghh6i3zzrrLOMd73iH+hxW5Ln/+Mc/rrvtne98p3HppZfW3fb973/f2GGHHdTbv/zlL42BgQHjhRdemHz/L37xC8uPRUTBx0wdEfnu7W9/O6655prJ67KtWrN06dJp9//1r3+tth6ffPJJjI2Nqe3EbDaLTCajtltlpe7DH/5wT89JVss+8IEP1N122GGHYcWKFWrFT7aHxeLFiyffX9tOfeWVVzr6XOPj43j66adx8sknqxXEGvm6ao0bsm38rne9C3vttZfKIP7DP/wD3v3ud7f8uI899pha4aytzAl57rXXSr5GWaXccccdJ98vq5hEFE4s6ojId1LE7bHHHk3fZyb5NyloZJtTipXZs2er7VUpiCR3JkVdKpXy6JkDg4ODddelsJNtzU7U8oDXX389Dj744Lr31YpHycU988wz+MUvfqGKWtkaPuKII3Dbbbe1/LiSoZMt1kaSoSOiaGFRR0ShIlk3KZokyybZOiE5MzNZPbvrrrtUQWNFMmSyYtWKNCzIKpeZXF+4cOFkoeWU7bffXq2W/fnPf8YnPvGJpvcbGRnBscceqy7HHHOMWrHbtGmTKmyluGz8mqQQlPxds4JZvsbnn39e5fCkAUM88MADjn5tROQdFnVEFCpSoBQKBVx11VWq+cHcQFFz3nnnqe5YaVw47bTTVBF39913qy1Z6aKVBgNpDJBVPxkjIkVRo89+9rM46KCDcNFFF6kiShouVq5cWdeR6iQpQD/96U+r7VYp1nK5HB5++GHVsXvOOeeojl8pvA444ABVzN56661qq3fmzJnq8fI1SSErW8TSrSsdv9J4Iaua0pUrRaA8TrZk165di4svvlit9EmRunz5clx++eVqK/tLX/qSK18fEbmP3a9EFCpLlixRBc6//du/Yd9991VdoLXRHjVSqNx5552qgFm2bJnKicloj4GByv9jP/e5z6nVNumOlXl4zz333LTPI6tcsgIoY0bk80iB9NWvflVl29xwyimn4IYbblAdrFKQHn744apTddddd1Xvl5EpX//611WnrRSbUpDefvvtk6uVsnIp40gkIyeFn5AxLj/72c/UayGPedOb3oRvfOMbahyKkMf++Mc/xsTEhHqd5DmY83dEFC6adEv4/SSIiIiIqDdcqSMiIiKKABZ1RERERBHAoo6IiIgoAljUEREREUUAizoiIiKiCGBRR0RERBQBLOqIiIiIIoBFHREREVEEsKgjIiIiigAWdUREREQRwKKOiIiIKAJY1BEREREh/P4/BvMwmuyeWLUAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Max Z-score Plot\n", "plt.step(rho_mxz, beta_mxz, where='post', color='blue')\n", "plt.fill_between(rho_mxz, cb_mxz[:, 0], cb_mxz[:, 1], alpha=0.2, color='blue', step='post')\n", "# Surflex Plot\n", "plt.step(rho_srf, beta_srf, where='post', color='red')\n", "plt.fill_between(rho_srf, cb_mxz[:, 0], cb_mxz[:, 1], alpha=0.2, color='red', step='post')\n", "# ICM Plot\n", "plt.step(rho_icm, beta_icm, where='post', color='gold')\n", "plt.fill_between(rho_icm, cb_icm[:, 0], cb_icm[:, 1], alpha=0.2, color='gold', step='post')\n", "# Formatting plot\n", "plt.xscale('log')\n", "plt.xlabel(\"Fraction Tested\")\n", "plt.ylabel(\"Hit Enrichment\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "50e7b044-1987-4662-b7f6-431ae1d31247", "metadata": {}, "source": [ "The plot closely aligns with the complete Figure 5.\n", "\n", "## References\n", "\n", "Ash JR, & Hughes-Oliver JM. (2022). \"Confidence bands and hypothesis tests for hit enrichment curves\". *Journal of Cheminformatics*, 14(1), 50.\n", "\n", "Empereur-Mot C, Zagury JF, & Montes M. (2016). \"Screening explorer–An interactive tool for the analysis of screening results\". *Journal of Chemical Information and Modeling*, 56(12), 2281-2286.\n", "\n", "Kosorok MR. (2008). Introduction to empirical processes and semiparametric inference. New York, NY: Springer New York.\n", "\n", "Zivich PN, Cole SR, Greifer N, Montoya LM, Kosorok MR, & Edwards JK. (2025). \"Confidence Regions for Multiple Outcomes, Effect Modifiers, and Other Multiple Comparisons\". *arXiv:2510.07076*" ] } ], "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.13.7" } }, "nbformat": 4, "nbformat_minor": 5 }