From 779123f0e96030aa1b2bf6bc76517d724017f923 Mon Sep 17 00:00:00 2001 From: seanohagan Date: Tue, 5 Aug 2025 14:10:21 -0500 Subject: [PATCH 1/2] add ai powered bayesian inference via ppi_py/ai_priors.py, with examples/galaxies.aip.ipynb and tests/test_ai_priors.py --- examples/aip_utils.py | 249 ++++++++++++ examples/galaxies_aip.ipynb | 547 ++++++++++++++++++++++++++ ppi_py/__init__.py | 1 + ppi_py/ai_priors.py | 750 ++++++++++++++++++++++++++++++++++++ tests/test_ai_priors.py | 253 ++++++++++++ 5 files changed, 1800 insertions(+) create mode 100644 examples/aip_utils.py create mode 100644 examples/galaxies_aip.ipynb create mode 100644 ppi_py/ai_priors.py create mode 100644 tests/test_ai_priors.py diff --git a/examples/aip_utils.py b/examples/aip_utils.py new file mode 100644 index 0000000..2a72776 --- /dev/null +++ b/examples/aip_utils.py @@ -0,0 +1,249 @@ +import os +import numpy as np +import matplotlib.patheffects as pe +import matplotlib.pyplot as plt +import seaborn as sns +import pdb + + +def plot_interval( + ax, + lower, + upper, + height, + color_face, + color_stroke, + linewidth=5, + linewidth_modifier=1.1, + offset=0.25, + label=None, +): + label = label if label is None else " " + label + ax.plot( + [lower, upper], + [height, height], + linewidth=linewidth, + color=color_face, + path_effects=[ + pe.Stroke( + linewidth=linewidth * linewidth_modifier, + offset=(-offset, 0), + foreground=color_stroke, + ), + pe.Stroke( + linewidth=linewidth * linewidth_modifier, + offset=(offset, 0), + foreground=color_stroke, + ), + pe.Normal(), + ], + label=label, + solid_capstyle="butt", + ) + + +def make_plots( + df, + plot_savename, + n_idx=-1, + true_theta=None, + true_label=r"$\theta^*$", + intervals_xlabel="x", + plot_classical=True, + aip_facecolor="#E6D7FF", + aip_strokecolor="#8B5CF6", + classical_facecolor="#EEEDED", + classical_strokecolor="#BFB9B9", + imputation_facecolor="#FFEACC", + imputation_strokecolor="#FFCD82", + empty_panel=True, +): + # Make plot + num_intervals = 5 + num_scatter = 3 + ns = df.n.unique() + ns = ns[~np.isnan(ns)].astype(int) + n = ns[n_idx] + num_trials = len(df[(df.n == n) * (df.method == "AIP")]) + + aip_intervals = df[(df.n == n) & (df.method == "AIP")].sample( + n=num_intervals, replace=False + ) + if plot_classical: + classical_intervals = df[ + (df.n == n) & (df.method == "Classical") + ].sample(n=num_intervals, replace=False) + imputation_interval = df[df.method == "Imputation"] + + xlim = [None, None] + ylim = [0, 1.15] + + if empty_panel: + fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(9, 3)) + else: + fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(6, 3)) + sns.set_theme(style="white", font_scale=1, font="DejaVu Sans") + if true_theta is not None: + axs[-2].axvline( + true_theta, + ymin=0.0, + ymax=1, + linestyle="dotted", + linewidth=3, + label=true_label, + color="#F7AE7C", + ) + + for i in range(num_intervals): + aip_interval = aip_intervals.iloc[i] + if plot_classical: + classical_interval = classical_intervals.iloc[i] + + if i == 0: + plot_interval( + axs[-2], + aip_interval.lower, + aip_interval.upper, + 0.7, + aip_facecolor, + aip_strokecolor, + label="AI Posterior", + ) + if plot_classical: + plot_interval( + axs[-2], + classical_interval.lower, + classical_interval.upper, + 0.25, + classical_facecolor, + classical_strokecolor, + label="classical", + ) + plot_interval( + axs[-2], + imputation_interval.lower, + imputation_interval.upper, + 0.1, + imputation_facecolor, + imputation_strokecolor, + label="imputation", + ) + else: + lighten_factor = 0.8 / np.sqrt(num_intervals - i) + yshift = (num_intervals - i) * 0.07 + plot_interval( + axs[-2], + aip_interval.lower, + aip_interval.upper, + 0.7 + yshift, + lighten_color(aip_facecolor, lighten_factor), + lighten_color(aip_strokecolor, lighten_factor), + ) + if plot_classical: + plot_interval( + axs[-2], + classical_interval.lower, + classical_interval.upper, + 0.25 + yshift, + lighten_color(classical_facecolor, lighten_factor), + lighten_color(classical_strokecolor, lighten_factor), + ) + + axs[-2].set_xlabel(intervals_xlabel, labelpad=10) + axs[-2].set_yticks([]) + axs[-2].set_yticklabels([]) + axs[-2].set_ylim(ylim) + axs[-2].set_xlim(xlim) + + sns.despine(ax=axs[-2], top=True, right=True, left=True) + + aip_widths = [ + df[(df.n == _n) & (df.method == "AIP")].width.mean() for _n in ns + ] + if plot_classical: + classical_widths = [ + df[(df.n == _n) & (df.method == "Classical")].width.mean() + for _n in ns + ] + + axs[-1].plot( + ns, + aip_widths, + label="AI Posterior", + color=aip_strokecolor, + linewidth=3, + ) + if plot_classical: + axs[-1].plot( + ns, + classical_widths, + label="classical", + color=classical_strokecolor, + linewidth=3, + ) + + n_list = [] + aip_width_list = [] + if plot_classical: + classical_width_list = [] + for _n in ns: + trials = np.random.choice( + num_trials, size=num_scatter, replace=False + ).astype(int) + aip_width_list += df[ + (df.n == _n) & (df.method == "AIP") & df.trial.isin(trials) + ].width.to_list() + if plot_classical: + classical_width_list += df[ + (df.n == _n) + & (df.method == "Classical") + & df.trial.isin(trials) + ].width.to_list() + n_list += [_n] * num_scatter + + axs[-1].scatter(n_list, aip_width_list, color=aip_strokecolor, alpha=0.5) + + if plot_classical: + axs[-1].scatter( + n_list, + classical_width_list, + color=classical_strokecolor, + alpha=0.5, + ) + + axs[-1].locator_params(axis="y", tight=None, nbins=6) + axs[-1].set_ylabel("width") + axs[-1].set_xlabel("n", labelpad=10) + sns.despine(ax=axs[-1], top=True, right=True) + + if empty_panel: + sns.despine(ax=axs[0], top=True, right=True, left=True, bottom=True) + axs[0].set_xticks([]) + axs[0].set_yticks([]) + axs[0].set_xticklabels([]) + axs[0].set_yticklabels([]) + + plt.tight_layout() + os.makedirs("/".join(plot_savename.split("/")[:-1]), exist_ok=True) + plt.savefig(plot_savename) + + +def lighten_color(color, amount=0.5): + """ + Lightens the given color by multiplying (1-luminosity) by the given amount. + Input can be matplotlib color string, hex string, or RGB tuple. + + Examples: + >> lighten_color('g', 0.3) + >> lighten_color('#F034A3', 0.6) + >> lighten_color((.3,.55,.1), 0.5) + """ + import matplotlib.colors as mc + import colorsys + + try: + c = mc.cnames[color] + except: + c = color + c = colorsys.rgb_to_hls(*mc.to_rgb(c)) + return colorsys.hls_to_rgb(c[0], 1 - amount * (1 - c[1]), c[2]) diff --git a/examples/galaxies_aip.ipynb b/examples/galaxies_aip.ipynb new file mode 100644 index 0000000..514a597 --- /dev/null +++ b/examples/galaxies_aip.ipynb @@ -0,0 +1,547 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5db36929-b470-46a0-a291-facca0846b1c", + "metadata": {}, + "source": [ + "# Galaxy classification\n", + "\n", + "The goal is to determine the demographics of galaxies with spiral arms, which are correlated with star formation in the discs of low-redshift galaxies, and therefore, contribute to the understanding of star formation in the Local Universe. A large citizen science initiative called Galaxy Zoo 2 (1) has collected human annotations of roughly 300000 images of galaxies from the Sloan Digital Sky Survey (2) with the goal of measuring these demographics. The target of inference is the fraction of galaxies with spiral arms. This notebook shows that prediction-powered inference allows for a decrease in the requisite number of human-annotated galaxies by imputing labels via computer vision.\n", + "\n", + "1. K. W. Willett, C. J. Lintott, S. P. Bamford, K. L. Masters, B. D. Simmons, K. R. V. Casteels, E. M. Edmondson, L. F. Fortson, S. Kaviraj, W. C. Keel, T. Melvin, R. C. Nichol, M. J. Raddick, K. Schawinski, R. J. Simpson, R. A. Skibba, A. M. Smith, D. Thomas, Galaxy Zoo 2: detailed morphological classifications for 304 122 galaxies from the Sloan Digital Sky Survey. Monthly Notices of the Royal Astronomical Society 435(4), 2835–2860 (2013).\n", + "2. D. G. York, J. Adelman, J. E. Anderson Jr, S. F. Anderson, J. Annis, N. A. Bahcall, …, N. Yasuda, The Sloan digital sky survey: Technical summary. The Astronomical Journal 120(3), 1579 (2000)." + ] + }, + { + "cell_type": "markdown", + "id": "c8248171-8a76-461a-a31a-2d4604f02b10", + "metadata": {}, + "source": [ + "### Import necessary packages" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "806c9b57-90d0-41ca-a2e0-6b88d171225d", + "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 (\n", + " classical_mean_ci,\n", + " sample_ai_posterior,\n", + " calibrate_dp_alpha_empirical_coverage_estimate,\n", + ")\n", + "\n", + "from tqdm.auto import tqdm\n", + "from aip_utils import *" + ] + }, + { + "cell_type": "markdown", + "id": "5cf90ae6", + "metadata": {}, + "source": [ + "### Import the galaxies data set\n", + "\n", + "Load the data. The data set contains human-annotated labels indicating whether the galaxy has spiral arms (```Y```) and corresponding predicted labels based on computer vision (```Yhat```)." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a6da3138", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Dataset galaxies not found at location ./data/; downloading now...\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Downloading...\n", + "From: https://drive.google.com/uc?id=1pDLQesPhbH5fSZW1m4aWC-wnJWnp1rGV\n", + "To: /Users/seanohagan/projects/for-tijana-temp/ppi_py/examples/data/galaxies.npz\n", + "100%|██████████| 268k/268k [00:00<00:00, 2.96MB/s]\n" + ] + } + ], + "source": [ + "dataset_folder = \"./data/\"\n", + "data = load_dataset(dataset_folder, \"galaxies\")\n", + "Y_total = data[\"Y\"]\n", + "Yhat_total = data[\"Yhat\"]" + ] + }, + { + "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.linspace(50, 1000, 10).astype(\n", + " int\n", + ") # Test for different numbers of labeled examples\n", + "num_trials = 5\n", + "\n", + "# True mean (computed on all labels)\n", + "true_theta = Y_total.mean()" + ] + }, + { + "cell_type": "markdown", + "id": "f0d0322d", + "metadata": {}, + "source": [ + "### Define loss function\n", + "\n", + "We define the loss `squared_error` which we use to do inference on the mean." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "0e65423b", + "metadata": {}, + "outputs": [], + "source": [ + "def squared_loss(y, theta):\n", + " \"\"\"Squared loss for mean estimation\"\"\"\n", + " return (y - theta) ** 2" + ] + }, + { + "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": 5, + "id": "812f8fd5", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "a941aba2a9d24c778a5d38015d40acee", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/10 [00:00 0.5, 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": 6, + "id": "6077b2c4", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3IAAAEcCAYAAACCtobiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABW1UlEQVR4nO3deXhU5d0+8Puc2TOTmclCAmRICBDCKkJYLEUFWzfqhkux7pdatb62peJatVpt61bqT2x9q75udSlaF9S6lgq1LggKAoYthCWQkD2ZLZnlnPP8/phkZEgC2WfJ/bmuuWCeOXPmOZP1zrN8JSGEABERERERESUNOd4dICIiIiIiop5hkCMiIiIiIkoyDHJERERERERJhkGOiIiIiIgoyTDIERERERERJRkGOSIiIiIioiTDIEdERERERJRkGOSIiIiIiIiSDIMcERERUQoQQsDj8UAIEe+uENEgYJAjIiIiSgFerxcOhwNerzfeXSGiQcAgR0RERERElGT08e4AERFRV0TQj+CHy2LaTKcuhWSyxqlHREREiYFBjoiIEpsajncPiIiIEg6nVhIRERERESUZBjkiIiIiIqIkw6mVRESUuAxmGOdf26GNiIhoqGOQIyKihCXJOkgZrnh3gyilaJqG+tpahEMhWNPTYXc4IMucpEWUbBjkiIiIiIaI+tpalG3fjmAgAADQ6fVIt9vhys9Hdk5OnHtHRD3BP78QERERDQH1tbUo37kzGuIAQFUU+H0+7C4rQ31tbRx7R0Q9xSBHRERElOI0TcOBigooitLhMb1eD0VRcKCiApqmxaF3RNQbDHJEREREfVBbW4tzzjkHTqcT2dnZWLJkSaeBCQD++te/ori4GOnp6Rg/fjwef/zxmMcfeughuFwuWK1WzJ8/Hzt27OiXPnrcbrS2tMBoNEKn08U8pqkqDAYDWlta4HG7++X1iGjgMcgRERER9cHixYths9lQVVWFdevWYdWqVXjkkUc6HLdy5UrcfvvteP755+HxePD888/jjjvuwOuvvw4AeP7557F8+XJ8+OGHaGhoQElJCc477zwIIfrcx1AwCKFpkHU66PWxWyQoigJZp4MQAqFgsM+vRUSDQxL98d2BiIhoAIigH8FVy2PaTD/8BSSTNU49Ioq1a9cuFBUVobKyEiNHjgQAvPLKK7jllluwb9++mGMff/xxeL1e3HrrrdG2c889F6NGjcKjjz6KefPmYeHChfj1r38NAAiHw8jKysJbb72FBQsWHLUvHo8HDocDbrcbdrs95rHmpiZs27IFer0eQgi0+P0xj6dZrVBVFROnToUzI6NX7wURDS6OyBERUWILtcTeiBJIaWkpMjMzoyEOACZNmoSKigo0NzfHHHv99dfHhLja2lp88sknKCkpiZ5r6tSp0ccNBgOKioqwadOmPvfT7nDAkpaGcDgM+bCplQAQCoVgSUuD3eHo82sR0eBgkCMiIiLqJa/XC6s1doQ4LS0NAODz+bp8XnV1NU4//XSUlJTgoosuOuK5ujpPMBiEx+OJuXVFlmW48vOh1+sRCgY7rRvnys9nPTmiJMKvViIiIqJeslqtaGmJHSluv5+ent7pc9auXYtZs2ahuLgYb7/9dnTNWlfn6uo8999/PxwOR/Q2atSoI/Y1OycHY4qKkGa1QpKkmMeMRiPryBElGQY5IiJKXAYTDPOujLnBYIp3r4iipkyZgoaGBtTU1ETbtm7dCpfLBUcn0xSfeeYZ/OAHP8CSJUvw8ssvw2QyxZyrtLQ0ej8cDqOsrAxTpkzp9LVvv/12uN3u6G3//v1H7W92Tg6OmTEDBWPHxrS3trRA7WKnTSJKTAxyRESUsCRZD92wwpibJOuP/kSiQVJUVIR58+ZhyZIl8Hq92LNnD+677z5cddVVHY59/fXX8bOf/QxvvPEGli5d2uHxK6+8Eo899hg2bdqEQCCA2267Dbm5uTjhhBM6fW2TyQS73R5z6w5ZljHS5QIOGZUTQsDN0gNESYVBjoiIiKgPXnvtNSiKgsLCQsyZMwennXYa7rrrLgCAzWbDSy+9BAD47W9/C0VRcN5558Fms0Vv1113HYBIkPvVr36FRYsWYdiwYdi4cSPeffddGAyGfu2vEALhcLjDerzmxsZ+fR0iGlgsP0BERESUAo5UfqCd3+dDfW0tWvz+aJHwdnanEzNmzx6s7hJRH3FEjoiIiGgI8Pt8OFBRAa/HA73BANthm6h43W6ukyNKIgxyRERERClOCIH62lqEQyGkWa3Q6/UwWywdjjm89h0RJS4GOSIiIqIUFwgE0OL3w2Q2R0sPyLIMoyl2F9iGurp4dI+IeoFbfxERUcISwRYE1/xvTJtp/s8gmdLi1COi5KQqClRVhclsjmk3WywIBYPR+26OyBElDQY5IiJKYAJoae7YRkQ9otProdPpoKpqtAA5AJjNZngOOa7F54OqKNDp+SsiUaLj1EoiIiKiFGc2m5FmtSIYCODQDcsPH6FjPTmi5MEgR0RERJTiJElCdk4ODEYjWvx+KIoCIQQ0TYsZoQNYT44oWXDcnIiIOti5TkGw5ejHDTwTkHFXbNMXMoD+3SLdlAaMn80fiZTarDYbXPn50TpyoWAQsizDarPFrI1rbmqKXyeJqNv4U4uIiDqw2CTkjUuUSRu6QXmVcAAwmI9+HFEys9psSLNaEQgEomvhWnw+bNm4MXpMez05rpMjSmz8CiUiog4kSYIkS/HuxqAS3EOFhghJkmA5pIacwWAAJCn6RdC+Ti4zKyteXSSibkiUP7cSERERURzo9XrY7faYNq6TI0p8DHJEREREQ5wzMzPmPtfJESU+BjkiIiKiIc6ZkRFzv32dHBElLgY5IiIioiHO7nRG1sm1YT05osTHzU6IiKgDnUFAaENr9w9JGlqbuxAdSq/XI91uh/eQ8Nbc2MgNT4gSGIMcERF1MGLc4Gz5fzQi6EfwvQdi2kwLb4NkssapR0Spy5mRERvkuE6OKKFxaiURERERddjwhOvkiBIbgxwRERERwcF1ckRJhVMriYgocelNMMxe3KGNiPof18kRJRcGOSIiSliSTg9d3pR4d4NoyOA6OaLkwSBHRERHtP79obFGZtbp/JFI5MzMxP69e6P329fJ6fT8+iBKNPyqJCKiI8oaISHblfpLqkOtgNES714QxVd0nZyIlB9pXyfH6ZVEiYdBjoiIjkiSJMhy6tdYE0OrbB5RpzpbJ+duamKQI0pAqf8nViIiIiLqNmdGRsz95sbGOPWEiI6EQY6IiIiIog6vJ+dxu6Gqapx6Q0RdYZAjIiIioqjO6sl5mpvj1h8i6hyDHBERERFFta+TOxTLEBAlHm52QkRERyGgaam/E4gkpf6GLkTd1aGeHNfJESUcBjkiIjqiwmm6uL22CPoRfO+BmDbTwtsgmaxx6hHR0HB4Pbn2dXI6Xfy+HxBRLAY5IiJKXHoj9DMWdWgjooHVWT05T3MzMliGgChhMMgREVHCknQG6AtmxLsbREOOXq9Heno6vB5PtK25qYlBjiiBcLMTIiIiIurg8DIEXCdHlFgY5IiIiIiog8MLg7OeHFFi4dRKIiIadF/+M4wWz9GPS1VpdmDOGYZ4d4PoiBwZGVwnR5TAGOSIiGjQ+ZuBEYUyRowZmhNDDu7W4t0FoqPiOjmixMYgR0REcSFJgCwPzdptLFlHycKZmRkb5LhOjihhDM0/hRIRERHRUXGdHFHiYpAjIiIiok5F18m1aV8nR0TxxyBHRERERJ1qXyd3qOampjj1hogOxSBHRERERF1iPTmixMTNToiIaNBZnYAAoGki3l2JC4Mp3j0g6j5nRgb2790bvd++Tk6n08WvU0TEIEdERINvqNdQc47gL8CUPFhPjigxcWolEREREXWJ6+SIEhODHBEREREdEdfJESUeBjkiIiKiPqitrcU555wDp9OJ7OxsLFmyBIqiHPE5r7/+OsaMGRPTpmkabDYbrFYrbDZb9Ob3+wey+93CenJEiYdBjoiIiKgPFi9eDJvNhqqqKqxbtw6rVq3CI4880umx4XAYDz30EC688EJomhbz2NatWxEOh9HU1ASfzxe9Wa3WwbiMI2I9OaLEw81OiIhoULzzlxA8DUNnl0p7loQz/8cY727QANu1axfWrFmDyspKpKWlYcyYMbjrrrtwyy234Oabb+5w/CmnnAKz2YzbbrsNL7zwQsxj69evxzHHHAOjMfE+b9rXyXk9nmhbc1MTNzwhiiMGOSIiGhRZw2XMOW3o7Na46xtOOxsKSktLkZmZiZEjR0bbJk2ahIqKCjQ3N8PpdMYc/8ILL8DlcuG5557rcK7169ejtbUVs2bNwt69ezFx4kQ88MADmDt3bqevHQwGEQwGo/c9h4SsgeDMzIwNclwnRxRXnFpJRESDQpIBWScNmZvEn7BDgtfr7TD1MS0tDQDg8/k6HO9yubo8l8ViwZw5c7By5UpUVFTgrLPOwqmnnoo9e/Z0evz9998Ph8MRvY0aNaoPV3J0XCdHlFj4Y4aIiIiol6xWK1paWmLa2u8fvmX/0SxbtgxPP/008vLyYLFYcNNNNyE/Px/vvvtup8fffvvtcLvd0dv+/ft7dxHd5DgsyHGdHFF8McgRERER9dKUKVPQ0NCAmpqaaNvWrVvhcrngcDh6dK477rgDGzdujGkLBoOwWCydHm8ymWC322NuA0mv1yP9sNdgPTmi+GGQIyIiIuqloqIizJs3D0uWLIHX68WePXtw33334aqrrurxub799lv88pe/RHV1NYLBIO699154PB4sWrRoAHreO6wnR5Q4GOSIiIiI+uC1116DoigoLCzEnDlzcNppp+Guu+4CANhsNrz00kvdOs+zzz6LsWPHYtq0acjKysKaNWuwatUqZB4WnuKJ6+SIEockhBg6e0ETEVHcbPhQgWvc0Nm18sAuFTNO5ebQNHg8Hg8cDgfcbveATbNUFAWffvxxTNu0khKWISCKA/6EISKiQTHUQk3O2KF1vTQ0tK+TYz05ovjj1EoiIiIi6rYO6+S44QlRXDDIEREREVG3dVgn19zMdXJEccAgR0RERETd1mk9Obc7Tr0hGro4gZ+IiBKWCAcQ3rAyps0w4xxIBnN8OkREna+Ta2xERgLtrkk0FDDIERFR4tJUaFWlsW3HnhmfvhBRlDMzs8OGJ0Q0uBjkiIjoqJ77dQBN1XGoViMkQLsltu1tCZBaB/RlM4ZLuOIPHPUj6oozIwP79+6N3m9fJ6fTDZ0SI0TxxiBHRERHNbxAhx9eNHR+ZHz7mRLvLhAltK7WyXF6JdHgGTo/lYmIqNckCdDppXh3Y9BIQ+dSiXqF6+SI4o+7VhIRERFRj7GeHFF8McgRERERUY+xnhxRfDHIEREREVGPOZzOmPusJ0c0uBjkiIiIiKjH9AYD0u32mLbmxsY49YZo6OFmJ0REdFSWdEBV4lB+IE4s6fHuAVFycGZksJ4cUZwwyBER0VGd8BNDvLswqFxThtb1EvWWMzMT+/fti95nPTmiwcMgR0RECUuEAwhvfi+mzXDMQkgGFusmSgRdrZNjGQKigccgR0REiUtToVVsjG2bcmp8+kJEHbSvk2M9OaLBx81OiIiIiKjXDi9DwHVyRIODQY6IiIiIeu3wwuAet5v15IgGAadWEhFR4tLpoRv3/Q5tRP1B0zT84x//QHl5ORRFiXnsN7/5TZx6lXw6rJPTNK6TIxoE/GlIREQJS9KbYJh6Wry7QSnquuuuw8svv4xp06bBaDRG2yVJYpDrAa6TI4oPBjkiIiIakl599VWsX78eEydOjHdXkh7ryRENPq6RIyIioiHJ4XCguLg43t1ICVwnRzT4GOSIiIhoSLrooouwbNmyeHcjJXS2Tq7m4EG0trZCCBGfThGlOEnwq4uIiIiGkMLCQkiShHA4jMrKSjgcDmQctoX+7t2749S73vN4PHA4HHC73bDb7YP++l+vXRszvdJkNsOWno50ux3DcnNhtdkGvU9EqYxr5IiIiGhIueeee+LdhZRktdliglwoGIQfQGtLC/w+H0aPHcswR9SPGOSIiChhiXAQSulHMW36yadAMpji1CNKBZdffjkA4I9//CNuuummDo/fddddg92lpCeEQCgc7tAWDocBRYGiKKg6cADjioshSVKcekmUWhjkiIgocWkK1D3rYpr0E08CwCBHvVNfX4+tW7cCAO6++27MmTMnZg2X2+3GI488gvvuuy9eXUxKgdZWtPh8HdolSYIAoCgKGuvrEcjPhyUtbfA7SJSCGOSIiIhoyDAajTj//PNRX18PADjxxBNjHjeZTLj22mvj0bWk1uL3IxwKdWgXQkCn10MJhxEKBtHi9zPIEfUTBjkiIiIaMux2O2prawEAEyZMwPbt2+Pco9QQCoUghIAkyxCaFm3XNA06ISDrdNBUFaFOwh4R9Q6DHBERJS5ZD13h7A5tRP2BIa7/GAwGQJIgCYFDt0NvXycnyTIkSYocR0T9gj8NiYgoYUkGEwzHnhnvblCKaS8/cCTJWH4gntJsNphMJgQDgQ6PCSEATYPRbEYad60k6jcMckRERDSktJcf2LBhA1auXImlS5di7Nix2L9/P5YtW4Zzzjknrv1LRhaLBZnZ2airroaqaVAVJebx9rVyZrM5Tj0kSj0sCE5ERERD0rRp07BixQpMnDgx2rZr1y4sXLgQO3fujGPPeifeBcH9Ph/2lJfD7/VCUZRONz8pHDcOBWPGDHrfiFIRgxwRERENSenp6WhsbIxZt9Xa2oqcnBx4vd449qx34h3kgEiYq6upgc/jQaC1Fb5OShKMnzQJI12uOPSOKLXI8e4AERERUTzMnDkTN910E4LBIACgpaUFP//5zzFv3rw49yx5WW02FIwZg7ETJmDC1KkYW1zc4ZidW7eirqYmDr0jSi0ckSMiIqIhaceOHfjRj36E/fv3Izs7G3V1dSguLsZ7772HUaNGxbt7PZYII3KdObBvH3bt2BHTJkkSjpkxAxlZWXHqFVHyY5AjIqKEUnXgQGyD0GLvS5xMkmriOc1OURR89tlnqKqqwqhRozB37lzIcnJ+jiVqkAOAPbt2Yd9hO4HqdDpMmzkTdocjTr0iSm4MckRElFCaGhthSUuLdzdoEOn1euj1g7eR9oEDB+ByuVBRUdHlMfn5+YPWn/6SyEFOCIGd27bh4GF/qDEYDJg+ezbSrNY49YwoeTHIERFRQmlqbOQvdUOMTqcb1CBnt9vh8XggtxWpBiJBQ5Kk6L+qqg5af/pLIgc5IPIeb928ucP6OJPZjOmzZ7M0AVEPsY4cERERDSmlpaUAgFGjRmHp0qWYMWNGUq6JSzaSJGHi1KlQwmE0NTZG24OBADZ//TWmz5oFg9EYxx4SJZfknARORERE1Evtoe2WW27Bu+++ix/+8If48Y9/jBUrViAUCqGgoCDOPUxdsixj8rHHIv2wEcMWvx9bNm7sUEiciLrGqZVERJRQOLVy6BnsqZWH83q9eP/99/HOO+9g5cqVKCwsxObNm+PWn97q7tRKoQl4GgRCAcBoBuxZEiRZGsSeAqFQCBvXrUNrS0tMe0ZWFqZOn560G84QDSZ+lRAREdGQ5fF48O9//xv/+c9/sG7dOiiKguHDh/foHLW1tTjnnHPgdDqRnZ2NJUuWQDnKyNLrr7+OMWPGdGh/6KGH4HK5YLVaMX/+fOw4bNv+vmqo1LD2nTBWvxTCmpdDWP1SCGvfCaOhUjv6k/uR0WjEtJISmEymmPamhgZs//ZbcJyB6OgY5IiIKOEIIXgbQrd4Oe6445CTk4N7770XZrMZy5cvR1NTEz766KMenWfx4sWw2WyoqqrCunXrsGrVKjzyyCOdHhsOh/HQQw/hwgsvhKbFhqfnn38ey5cvx4cffoiGhgaUlJTgvPPO67f3qKFSw+dvhrH9CxV1BwTcDQINBwW2f6Hi8zcHP8yZLRYcU1ICvcEQ015bXY1d27czzBEdBadWEhER0ZA0f/58bNiwAbNnz8Ypp5yCU045Bccee2yPzrFr1y4UFRWhsrISI0eOBAC88soruOWWW7Bv374Oxy9YsABmsxkzZ87ECy+8gL1790YfmzdvHhYuXIhf//rXACKhLysrC2+99RYWLFhw1L4caWql0ARWvxxG+QYVoSDQ6om024cBegOghoFxM3SYf5Fh0KdZepqb8c3XX0M7bKfQ0WPHYvTYsYPaF6JkwhE5IiIiGpLWrFmDgwcP4pe//CX27duHCy+8EMOHD8ell17a7XOUlpYiMzMzGuIAYNKkSaioqEBzc3OH41944QW8//77GNtJQCktLcXUqVOj9w0GA4qKirBp06aeXVgn3HUa9m9TAQlo9X7X7qkDDCZAkoD921W46wZ3VA4A7E4npkybFi0F0W5veTkq9+8f9P4QJQsGOSIiIhqyrFYrTj75ZPzoRz/CySefDE3TsHr16m4/3+v1wnrY5jxpbQXtfT5fh+NdLlePz9XZeQAgGAzC4/HE3LrSVAME/IDO0PGxxipAZwRafZHj4iEzOxsTDwmx7cq2bUNtdXUcekSU+BjkiIgoYQklhPC2j2NuQgnFu1uUIpYvX47TTz8dmZmZuPHGG2G1WvH222/jwIED3T6H1WpFy2E7L7bfT09P71F/ujpXV+e5//774XA4orcj1sITApoCtLgB3WEbhIaDgLcR0NTIcfGSM3w4iiZM6NC+bcsWNNbXR+8LIdDa2gqf14vW1laupaMhiwXBiYgocalhqNtjR0f0Y+YAehYNpr574YUXsGjRIixbtgyTJk3q1TmmTJmChoYG1NTUIDc3FwCwdetWuFwuOByOHp+rtLQUZ5xxBoDIGrmysjJMmTKl0+Nvv/123HjjjdH7Ho+nyzDnzJUgACghwGQFWppjHw8HImvlnLmDuz7ucHn5+QiHw9hbXh5tE0Lg202bcOzMmdDpdKivrUWL3w9VVaHT6ZBmtSI7JwdWmy2OPScafAxyRERENCStX7++z+coKirCvHnzsGTJEjz55JOor6/Hfffdh6uuuqrH57ryyitx991347TTTkNxcTHuuOMO5Obm4oQTTuj0eJPJ1GH7/q5IEmBKiwQ5NRxZFxcOfve40CJBTopvjgMAFIwZg3AoFLM+TlNVbP76azizsgAhYDKbYTKboaoqvB4PAoEAXPn5DHM0pDDIERHRgDlYWYlwONy3k0y7/LCT1gGo69s5aUAYDAaMyMuLdzcG3WuvvYYbbrgBhYWFkGUZl112Ge666y4AgM1mwxNPPIGLL774qOe58sor0dzcjEWLFqGurg6zZs3Cu+++C4Ohk4VtPRQOSbBnS5AkgVYvIPQAgrHHKKHIcfEmSRLGTZiAcDgcsz5OURQ01tVheF5etIC8Xq+HzmpFi9+P+tpapFmtHTZNIUpVLD9AREQDpmLvXqRZrfwr+RDh9/mQPWxYvLsxZB2p/IC7TsNnr4chJMDXIOBrFmhxA8HYJXlYdKMBY6cnxt/5NU3Dlo0b0dTQENOuNxiQnZMDSZIgyzL0ej0URYESDmP0uHGwWCxx6jHR4EqMr1QiIkpZkiTxL+RDBD/OicueJSFzpIyavRpGjpcQapURahXY/oUGVfnuuE0fqwkT5GRZxpRp0/DN11/D63ZH25VwGLUHD8JoMkGn08FgMMCSlgZN06AqyhHOSJRauGslERERUYqTZAlFJTqkpUtorIqshbNlSMgtjA3fu7/RULVr8GvJdUWn1+OY6dNhPmyUTdM0hMNhSJKEYDAId3MzNFWFTp8YIZRoMDDIEREREQ0BWXkyZpyiR+5oGQGfQFO1gD2r4yawn73Wx3Wt/cxgNGJaSQkkOfbXVlVR0NrSAlmWoSgKVE3r9uYvRKmAf7YgIiIiGiKy8mRkjpDgaRAIBQCjGcgtVPGfv383JXFfqYb921SMmqiLY08PI0lw2O1obm6OaVZVFX6fDwajETpZRjAY5Bo5GjI4IkdEREQ0hEiyBMcwGcNGyXAMk3HsD/WwOmOP+ex1JaEKbauKAshylyNu4VAIXq8XoWCw08eJUhGDHBERDRiDwQAhBG9D5CbL/LUiGRmMEo47K7bEwYEdGvZ9mzhr5WSdDuFQCJIsIy0trdONdULBILZt2QKvxxOHHhINPk6tJCKiAdPXmmJCCUHdvS6mTTdmNqTDF/VQQuCUtuQ1db4O695V4G34bhTu09fCKJgiJ8RupNEeCAG9yYR0vR6tra0d6lQGWlux4csvMaaoCK6CgoToO9FA4Z/OiIgocalhKKUfxtygJtZGDESpQG+QMPec2L/vV+8WKN+YGKNyqqrCaDRCp9cjHApBADCnpcFsNnc4VgiB8p07sXnDBgQ51ZJSGIMcEREREWHSPB2cubEjWJ+9FobQ4r9WTqfXw2gywWazwWQyRWvG6fR6pNvt0HdSdqCpoQFfff456uvq4tBjooHHIEdERERE0OklzF0UG4jq9gvsXB//UTmz2Yw0qxWqpsHudMKZmQlHRgacmZnIyMqCIyMDDqezw/PC4TC+3bgRO7dtg6qqg99xogHEIEdERIlL1kHOHR9zg5xAW6ITpZgJ39MhK++wUbk3wtDiPConSRKyc3JgMBrR4vcDiGymBAAtfj9MJhPGT5qEaTNnwtjJzpZV+/fj67Vr4fN6B7XfRANJEom0tywRUYqo2LMH4VAo3t2gfmQwGpFfWBjvbhB1yePxwOFwwO12w2639/o8O9apeOex2O9fp19rwOR58d8jz+/zob62Fi1+PzRNgyzLSLNakZ2TA6vNBiBSimDH1q2or63t8HxJljF2/HjkjRrFjVAo6THIERENgOqqKqQ7HPHuBvUzo9EYHQUgSjT9FeSEJvDCb4Ko3ffdr4iOHAlXPmiCTh//8COEQCAQiK6RM5vNHUKZEAIHKyuxa/t2aFrHqaGZ2dmYMHlyp6N3RMmCQY6IaADUHDzIIJeC9Ho9jEaWPqDE1F9BDgDKN6p480+xo3KnXGnAMQviPyrXE36/H9s2b+50SqXBaMSEKVOQlZ0dh54R9R3XyBERERFRjDHHyhgxNnaU64u3FCjh5Pr7v9VqxYw5c+AqKOjwWDgUwpYNG7Br+/aYjVCEEGhtbYXP60Vrays45kGJikGOiIiIiGJIkoR558dOI/Y2CGxek3w7P8qyjHHFxThmxoxOR9QPVFRgw5dfwu/zwe/zoWLPHuzdtQt7y8uxd9cuVOzZA7/PF4eeEx0ZgxwRERERdZA/WYZrQuyvil++FUY4mJwjVJnZ2Zg5d26nUyn9Ph++WrsWO7Zuhcftht5gQJrVCr3BAK/HgwMVFQxzlHAY5IiIiIiog8ioXOyaOL8b+GaVEqce9Z3RaMSU6dNRNGECZDn212ChafA0N6PF74ckSZAkCXq9HmlWK8KhEOpraznNkhJKcq1YJSJKErJOxx/4KYjbldNQ4yrWYfRUGXu3fLfz47p/Kph2kh5GS3J+PUiShLz8fDgyMrBty5YOI22tLS2oOnAAzowMWNLSoNfrYTKb0eL3IxAIwGKxxKnnRLEY5IiIBsCwnJx4dyEliKAfwfceiGkzLbwNkskapx4RDT3fP8+AvVuC0futPuDrDxV875zkLsVhS0/HjDlzsLusDJUVFTGPaaqKxvp6yLIMs9kMu9MJTdOgKsk7Gkmph1MriYiIiKhLI8bKGDsj9lfGr95XEPAn/6wDnU6HogkTUDx5cqcj7pqmoaWlBdVVVfB5PFAY5CiBMMgREVHikmRI2aNjbpD4o4tosH3/3NjRt2BLJMylitwRI5Bmsx3xmNbWVnyzfj22bNyIhro6Tp+nuOPUSiIiSliS0QLT8VfFuxtEQ15OgYziOTrs+PK78gNff6hgxil6pNmTc63coQKBADRVhSzLEIhsfNKVhro6NNTVwWyxYKTLheF5eZ2WNSAaaPyzJhEREREd1dxFehw6+zAcANa9mxqjcq1+P1RVhdligdFohE6vj+xqeYQNjgKtrdhdVoYv/vMfbN28Gc1NTRylo0HFETkiom7Y8OWX8e4C9aMZc+bEuwtESScrT8bEuTps/ey7Ublv/qVg5ul62JxJPionSZAQKR6uNxhg1DQIISBJEjQhEAoEEA6HO32qEAK11dWora6G1WbDSJcLuSNHQq/nr9k0sPgZRkTUDRlZWcjopIgsJadgMAiTyRTvbhAlnbmL9Nj2hQrRNvNQCQNfvh3GDy5L7qmFaWlp0BsMUMJhGHW6mBpzEiIlZdKMRuQMH4662lr4vd5Oz+P3+VC2fTvKy8qQO2IE8kaNgi09PeYYTdPgcbsRCgZhNJlgdzg61LQj6g4GOSKi7pAk/qBNIZz+RNQ7zlwZU0/QYfOa70blNq9WMWuhBnt28n6PNFsscGZmor6uDqFgEAaDAZIsQ2gawuEwJElCRlYWCsaMQcGYMfC43ajavx+1NTWdrqfTVBUHDxzAwQMHYHc4MHLUKAzLzUVTQwP279uHVr8fmqZBlmVYrFaMKihANsvWUA8xyBERERFRtx13th6ln6pQ25bHqQqw9i0Fp1yVvKNykiRhpMuFUDAIv88XU2ZAluXolMn2EgUOpxMOpxPjiotRXVWFqgMH0NrS0um5PW43PG43yrZtAyQpEuAkCZAkqKoKr8eDXTt2AADDHPVI8v7phIiIiIgGnT1bxjELdDFtWz5R0VTT9U6PycBqs2H02LHIGT4caVYrzGYz0qxW5IwYgdFjx8LaSXkCg9GIUaNHY/b3v49jSkoiQayLDVJUVYWqKBCaBk3TomvyJElCKBjEvj17oB1ht0yiw3FEjoiIiIh6ZM5ZBmxZo0Jp2/9DaMAXbypYeF3yjsoBkTCXZrUiEAhAVRTo9HqYzeZOi4UfSpIkZGZlITMrC8FAAFUHDuBgZSVCwWCnxwshoqN+7ef2e71obmpCZlZW/14UpSyOyBERERFRj9icEo49OXY8YOvnKhoqk39ESZIkWCwW2NLTYbFYjhriDmcym1E4bhyOO/54TJ42DRmZmUc8XggBIQQ0TUPppk3YXlqKupqaLnfJJGrHETkiom46UoFYSi49/cWMiDqa/SM9Nn2sIBxoaxDA528qOPOG5B6V6y+yLGNYbi6G5eZi3+7d2LNr11GfoyoKqisrUV1ZCUgSHA4HMrOzkZmdDVt6Or93UQxJcOsuIiJKUCLUitBnz8e0Gb9/OSSjJU49IkpcHo8HDocDbrcbdrt9UF7z09fCWPtWbFHwy35nQk4BJ30dqsXvx8b16xEOh6HT6aCpao93zzUYjZHpm9nZyMjKgtHYeWAWQvR4aiglJwY5IiIiohQQjyAX8As8dWMAwUM2bBw7XcZJlxoQCgBGM2DPkiDJQztICCGwfcsW1NfVQRMCep0usoOlqkJV1aOfoBPpdnt0tC7dbocsy/D7fKitroanuRmKokCv18PudCJn+PBON2uh5MYgR0RERJQC4hHkAGDtW2F8+lrsqJxrggSdHjCYJAwvlFA0U4+svKE9Suf3+bBrxw74PB6oqgqBSLFxnU4Ha3o6XPn5aG1tRWN9PZqbmno0nV+v1yPd4UAgEEA4GIQmBCAEIEnQyTKsdjuKiosZ5lIMgxwRERFRCohXkAu1RkblWn3ftemNQHrb5ouyTkJ2noS5iwwMcz4f6mpq4D5kxMzhdGJYbm5MyFJVFc2NjWhsaEBjfX2XNeqOpL20gSYEZElCVk4OJk6ZwmmWKYRBjoiIiCgFxCvIAcC6f4bxySuxo3LO4YDBCIRDgKYC42bosOAiA6dZ9mINW2tLCxrr69HY0ICmxkZovZiOKcsyxk6YgNzhw6HXc7/DVMAgR0SUItZ//jlaevFXW0oOaWlpmDV3bry7QQksnkGufr+Cv90ZxqGzASUZMFkBoylSZ85sA87+hRHOXF3XJ6Kj0jQN7qamaLDz+3xHf9KhJAnpdjucGRlwZmbC4XQy2CUpftSIiFKENT0dY8aPj3c3aAC1trbCYuGOnZR4vI2R6ZShwHdtQovc2ts0D9BULeDMjU8fU4Usy8jIykJGVhbGAggEAthdVobagwe7dwIh4HW74XW7sX/vXkCSYLfb4czMhDMjA/ajBDvuipk4GOSIiFKEJEmQ5aG9/iTlcRINJSghAAFEdu845NM06I8EPEmOTK/kp3D/M5vNGJmXh4baWmiaBp1eD6Fp0DSteyUOhIDH7YbH7UbFnj2Q2kfsOgl27Wv8vG0btuh0OqTb7R3W+NHgYJAjIiIioj4xmAFZBnR6QA3HPqaEIv/KOiAcYpIbCHanE1abDV6vF5qmQa/TRQOdoqoQmhapX9eNcCe6CHZWmw0tfj9CoVAkKLZtotLa0gK/z4fRY8f2S5jjiF/3McgRERERUZ+YLBLM6UCLJzKdUutkLw5NBd79XwUHdwFzztTDks5fzvuLLMvILyzErh07EA6FYmrTyZIEg8WCccXFyMjKgru5Gc2NjWhuaoLX7e5RsDuUJEnRgBUKhaDT6zF+4sQ+zQxpH/HzeTxQVBV6nQ42jvh1iUGOiIiIiPrElCYhe6SMBmho9QNq6LuRuENpCvDV+wo2r1Ewa6EeJafpYTQz0PWH7JwcAMCBffvQ0tICoWmQZBlpVitc+fnRxzOzspCZFakNoSoK3G53j4JdOyFE9FhN01BdWYmaqiqYTCaYzObozXzI/01mMwwGQ6cjbH6fD3vKy9Hi80XPGwTQ0tICn8+Hwn4a8Usl3LWSiChFbPv2WwwfOTLe3aABZDabYUlLi3c3KEHFc9dKoQl8+U8F+7epUBSBFnckyIWDkXVyXf22mWYHjjvHgGkLdNDpGej6g6Zp8LjdCAWDMJpMsDsc3R4liwl2jY3wejzdDnbdJctyJNQdFvga6+vhdbshyTJ0Ol007KlqZHFlVk4OxhUXc5rlIRjkiIhSxK4dO5DR9ldWSj0SAIvVyl0rE1BtbS2uueYarFmzBnq9Hpdccgn++Mc/drrz33vvvYdbb70Vu3fvRn5+Ph5++GGcccYZACK/gNvtdgghYn5ZrampgdVqPWo/4hnkAKChUsOGjxS0eARMaZENToQG+N0Cnnqg8aCA0Dp/rj1bwrzz9ZjwPR3kIV5nLpEoigJPczMOVFSgsb5+0F+/fRMvSZZhNJkw5dhjkcY/ZkUxyBERERH1wYIFC5CXl4cnn3wS1dXVOOuss3D55Zfj5ptvjjmurKwMxxxzDP7+97/jjDPOwBtvvIErrrgCZWVlyMvLw7fffouSkhJ4vV4YjcYe9yPeQQ6IhLmyr1U0VmlQwoDeAGSNlDGuRAedAfjsdQXbvlBjdrY8VPYoCcdfYMCYY2WOvCSQ+tpabN2yBUBkzV37tMpDp1cOBoPBALvTCVt6Omzp6bCmp8NisQzZzxUGOSIiIqJe2rVrF4qKilBZWYmRbVObX3nlFdxyyy3Yt29fzLF33nkn1q1bh48++ijadvrpp2P27Nn47W9/i2effRaPP/441q9f36u+JEKQAyLTLD0NAqEAYDQD9iwJ0iGjbLX7NHz6Whi7v+lieA5A3ngZx/9YD1dxx+LhRzs/9b+WlhZs2bABoWCwQ6kbra3UgdFoxJjiYkiI1LYLHnILBAIIBYMD0jedTgdrW7CLBjybDTpd54XnVVVFfV0dgq2tMFksyB42rMtjEx03OyEiIiLqpdLSUmRmZkZDHABMmjQJFRUVaG5uhtPpjDl26tSpMc+fNGkSNm3aBABYv349WltbMWvWLOzduxcTJ07EAw88gLlz53b62sFgEMFDfjn2eDz9eGW9J8kSHMO6DlY5BTLOXWrCgR0q/vuqgsqdHQNd5U4NK34XQuE0GcdfYEBOQSQ4dDbilzlSRlGJDll5rKM5UCwWCzKzs1FXXQ1NCGjaIR8zSYJer0dWTg6G5eR0OTqmaRpCoVBMuGv/v9/nQ2tLS6/6pqoqPM3N8DQ3x7SnWa3fBTybDTa7HbXV1Tiwbx/CoVB0CvNuoxGjRo+GKz+/V69/uMEsn8AgR0RERNRLXq+3w/q19jU8Pp8vJsh1dazP5wMQ+WV5zpw5uPfee5GZmYm//OUvOPXUU7F582YUFhZ2eO37778fv/3tb/v5igaPq1iHC++UsXuThv++Gkb9/o6TxPZs0rBncxATj9Nh0jwddn2twu/RYEqTYDJGShpU71HhbRCYcYqeYW6ASJKEkS4XgsEg/G216trJsgxbejpGulxHDCyyLMPctovl4bweD0o3bUKwfcRPkiCAaGHz3mjx+9Hi96OuurrT65FlGUIIBAMB7CkrgxACowoKevVa7Qa7YDqDHBEREVEvWa1WtBw2ktB+Pz09vVvHth+3bNmymMduuukmPPvss3j33Xdxww03dHjt22+/HTfeeGP0vsfjwahRo3p/MXEgSRLGHqtD4TEytn+h4rPXFbjrDgt0Atj2hYrta1VYM4D0TKChUkBVIgXIrQ4gHNCw62sVmSM4zXKgWG02FI4dG63z1h5U+qPOm95giP6RQ1EUQAhIACRZhqzTQafTwWw2Izs3F6FAAD6fDz6PB+Fw+Mgn7oIQIqbWnqqqKN+xA/v37IHRZIreTCYTjEZj7H2TqdNdQP0+H/aWl8Pv80FT1QErmH4oBjkiIkpYQlWg1ZXHtMnDxkLS8ccXJYYpU6agoaEBNTU1yM3NBQBs3boVLpcLDoejw7EbNmyIadu6dStmzpwJALjjjjtw/vnnY/r06dHHg8FglzuVmtp+sUwFsixh0vf1KJ6jw+bVKr5YGUbLYTNFhQB8jZGbzgDojYAkAa0+wJwmUGVQ4WnQHXFaJ/WN1WZDmtXa71MHzWYz7E4n1LYRuHAoFK2DZzAaIcsyMjIz4crPj76WEAKhUAg+jycS7Lxe+DyeXk/RBCKFzUOhEOD1HvE4vV4fE+4MRiO8bjd8Pl+k32191CQJQlGgKAqqDhzo9/IJ3OyEiIgSlgj6EXzvgZg208LbIJmOvhU70WA5/vjj4XK58OSTT6K+vh5nnnkmzj//fNxzzz0xx23fvh3Tp0/H888/j3PPPRdvvPEGLr/8cmzatAnjx4/H2WefjaamJrz66qvIyMjAgw8+iD//+c/Yvn07MjMzj9qPRNnspD+EAgIbPlSw7l0FodYuDpIia+QkOfL/NDtw9s+NyClIzo0rhjq/z4cDFRUIh0LQ6fWQ2nbHVBUFRqMRefn53RrRUhUFfp8vOmrn8/l6VOi8P7XXw9M0DSazGcfMmNGvtUA5kZiIiIioD1577TUoioLCwkLMmTMHp512Gu666y4AgM1mw0svvQQAmDBhAlauXIk//OEPyMjIwL333ovXX38d48ePBwA8++yzGDt2LKZNm4asrCysWbMGq1at6laISzVGs4Tjzjbgp8vMmLlQD7mzbCbaio4HgHAr4G0ADuxUITSOUSQjq80GV34+0u12oC3AQQik2+3dDnEAoNPrYXc6MdLlwvhJkzBj9mwUT5kSXRen0+kitekGoWSBqqqQZBmQJISCQfj9/n49P0fkiIgoYXFEjqj7UmlE7nDffhLGR88o0NSjH5ueKWH8bBkT5ugxfKw0ZGuMJauB2PVRVVWs++wzhIJB6A0GyIecT9U0KOEwDCYTxk+YgHA4jFAwGLmFQgi17Q4bCgZ7NapnNBojO32qKoomTMAIl6tP13IoLjIgIqIEJgFpzo5tRDSkWB0SLDYgHAaUII4Y6LyNAl9/oOLrD1TYsyUUz9FhwnE65BQw1CUDSZK6XBfaWzqdDqNGj8beXbughMORKY+yDKFpUFUVer0eBYWFGNa2zrUzQggoihINecFDwp7X7Yb7sPIHQGQtHYDoujlDP69p5YgcERERUQpI5RG55hoVby0PIeADIAOhAKCGAFUBRDd3p3fmSpgwR4fi43QYNoqri4aiAxUV2L93b0wdOYPJhFEFBX2qI9fa0oLN7QXTZTn650ZJkqJ190wmU7+vkWOQIyIiIkoBqRzkhCaw+qUwdm1UodN9t2OlEECwJXKTZESCXjdk5X03Upc5IjbUCU3A0yAQCgBGM2DPYkmDVKKqKurr6hBsbYXJYkH2sGHQ6fq2QY4QArt27EBtTQ00VY2ZNyIQGREclpvLXSuJiIiIqKNUDnIA0FCp4fM3w6ivFNDU7359lXUSsvMkzF1kgBDA9i9V7Firoqm6e7/iDsv/bqRODQM7v1JQs0dDOAgYTEBuoYzxM1lsnI6svY6cr4uC6QNRR45BjoiIiCgFpHqQAyJhruwrBdV7BMJBAYNJwvBCCUWHBS0hBOoqBLavVbHjS7VjkfEumK2ApIvUqZPbyhro9cCwUTLmLjIwzNER+X0+1NXUwOvxQFNVyDod0h0ODMvJ6fcQBzDIEREREaWEoRDkgJ5PfRRCoHqPwI61KnasU+Ft6N6vvpIMyLrIzWAEJhwnY8HFRk6zpCMaiF03u8IgR0RERJQChkqQ6wuhCVTt0rB9rYqd61T43d1/riQD+RMl5E/RYeRYGbmFMoxmhjqKHwY5IiJKWEJToDXsj2mTs0ZBklk9h+hwDHI9o2kCB3Zo2LFWxbYvVIRae/Z8SQaGjZIwYpyMkeNkjBgrI2M4SxzQ4GGQIyKihMWC4ETdxyDXe5vXhLDqeRWSFClpgF7+dmy2ASPGRoLdyHEyho+RYUqLDXZqSEXZRgFfo4AtU0LRdAk6Y992TaShiX/SJCIiIqIhzeaUoDcAkABLeiTMqeHIv5oSKXPQHQEfsGeThj2b2nYtlICskRJGFskYOVZGc52GbWtVtLgj9e8kGfivE5h5mh7Tf2gYqMujFMUgR0RERERDWuYIGVYn4GsClFBk10q9HtBE5L4aBowWoGCyjMYqgdoK0b1C5AJoqBRoqFSxZY0abZbbdsYEAG898N9XFQDolzDHOnhDB4McEdFg2/oEEGzo2zmm/7p/+pIMjGnx7gERpTjHMBmFU3XYuV6FEm4bjWt7TJIiZQnGz9ZhwUUGSLKEcFCgeo+Gg7s0HCzXUFWm9WjjFE2N3NqFWoGP/6Zg66cqMkfIyBghIWO4jMzhEpzDJRiM3QtiDZUayr5W0VilQQkDegOQOVJGUYmOpRNSENfIERENtl1/B+xj+3YO6yjAOqJ/+kNEKYFr5PomWnD8QCQEtU991BuAbNeR68gJERkFO7hLQ9WuSMCr2StiwlpfpGdJyBwuIWO4hIwRMjJHRP5vz5Ygt422NVRq2PCRAp9bhayTIuv8JEBTBWwOHWac0j9FzTnilzg4IkdENOikSMXZPuHf4IiI+lNWXiSsdafg+OEkSYIjW4IjW8aE4yJtSkigdp+Gql0CpZ8qqKvo/fdtb4OAt0FgXynw3VghoNMDjpxIqAv4BXxNAkoQCIUiUz9lHWCxAq0OFbu+lpA5om+hq6FSw86vFNTs0RAOAgYTkFsoY/xR3p+eYFDsPgY5IiIiIiJEwlzmCEO/BAm9UcLIIh1GFgG2TOCDp8LQ6SNTNTUFUFVAqIDWnbV2XVAVoLFKoLGq85AoSUC4FfA2As11Ctz1Ahm5Eix2CWnRG5CWLkF/lOmbh45Yaoe8XFO1ivr94ogjlt3FqaE9wyBHRERERNRGkiU4hvXvCFDRdAn/dQK+xkiJAoPpu8dUNbLbpdkGzDpVh+Y6oKlaoKlawNvYt9kX7QuohAq0uIHNq7ue62m0oEO4a79vSQd2rldxsFyLbARjAnS6SAhVw0D1bg2bVyuY37aGsDfap4a2eARMaYDJEJneWrNHg7dB9NvU0FTCIEdERERENIB0Rh1mnqbHp/9QEPABBnMkCKkqEA4ABiPwvbM7liAIBQSaawQaqzU0HRRtAU9D40GBYEv/9jHUCoRaI6/XXZIMyDIACfj2UxWtfgHHMBkWqwSzDTDbpO/+b5VgsiK6pu9QQhMo+1pFc60GVRGor4yMWsp6wOoAQgEJu75W+zw1NNUwyBERERERDbD2kPbVBwr8zUC4bTMVW2bXdeSMZgk5BRJyCmJHooQQaPUiEuqqBcq/VlD+jQBE92ve9QehAWrb1FA1DGz9VANw5LmiprRIwDNbAYtNgtkqQdZFNoppbYmcB4hMC5UkoNULmK0CVQYNngbR59HSVFqDxyBHRDTY9GnoXgGiI0nOHzo9JTQVwn0wpk1yjIAk93WzGCKiwTf9hwYcc4KMso0CvkYBW6aEoukSdMaefU+TpLbpj3Yd8sYDGbkSqveEoIQiW2HJushPCVWNrKMTWmTkbGRRJBC2eARaPJEw2F87a3ZXsAUItghEqjV0L3UG/ICvScPql8LIHS3Dni3BnhXZtdOWIUGn7355hh3rw6jcoUWDXF6xjOJZfV/fFw8sP0BERAlLBP0IvvdATJtp4W2QTNY49YgocbH8wNClKRpeXxZCbYWAzhAZ1Tq0fIISBnLyJZy31AhZ/11gEZpAoAVo9Yq2cAe0uEU06LW3eeo1ePpY/nSgSBJgy2wLdm3hzp4tIb39/1kSjGYJDZUaVr8cQs1eASWMaHkGvQHIHS1hwUXGfglzakjtc1DvLo7IERERERElMVkvY8YpenyyQkGgRcBii5QmUJXICFiaTcKMU/QxIQ6IbOxisUWmOGYeoTRp7T4Vby8PwecBdHKkgo6Ets1OlMionsEIjJooQ5IiI2gBn0CrTyDgB5TQwF27EN+VZ6js4hizLXJcqDUS/HSGyKhle9uBbQKrXwph3gUGGIwS9EZAb4j8qzOgbbfRo4/6bVwVjk6dbQ/S/3V2PXW2rxjkiIiIiIiS3NjpkV/rN3ykoKlGINQa2Swk2xUJce2P94YpTYqMVskagv5IABKIhCKDMbLuLWukjAUXG+AY1nFUKxwSCPrRFuwEAj60/RsJerUVGiq2ahBa2xq/fl7rF/B993+BzqeT7t0isHdLF4mzbeQuEvAiIU/Xdt9gjPzf746UgRBa5H3X6SOB19cIfPoPBQD6PcwxyBERUWLT9f9fMYmIUtHY6XoUTpVRVR6ZFplmlzByrNRhJK6n7FkSRoyToYQFVKeA3x27q6ROL2HkOBn2rM5HrQxGCQYjYMvo/PHmGhVvLQ8h4G+rfRdqC1tSW5kDJfKjYNQEGYEWwFMv4KkXAzrSF0NERhW/ez0R++Bh1HDkpjMAJisQ9Ec2uTnmBLlfp1kyyBERUcKSTFaYz/pNvLtBRJQ0ZL0MV3H/nlOSJRSV6OBtiAREZ05k2qDQ2qZu2iWMK9H1evdHxzAZoybqUL5BhawHjGmALAGaiAQiTQHGzdDF1KkTIjKy52kQ0WDnadDa/o3cb/X257vQc2o4ssmMwQz4m4GyjQIT5vTf+RnkiIiIiIjoiLLyIuvwyr5W0VilQQlEphkOL5QxrkTXp41CJFnCtAV6+JsF6g9okfCGth04ZSBnjIxjFuhjgqIkRQqVW9Il5I7u/LzhYCTUlW9U8PkbKrS23Tvbp3BGC6Zr7f3oh02lDxFZWxcZVQxrgK+PBd4PxyBHRERERERHlZUnI3OENCB12LLyZMxdZMDOrxTU7NEQDgIGUyQoFs3U9yooGkwSskZKyMjRo+wrDXVtu3rqjd+N+CmhyMhZToGEC39tBCQJSjiy06cSElDDbdMq2+4r4fbnCITDgBoCqnap2L5Wg6RrK5AuIpupGEyRfqhqe83A/i0dxCBHRERERETdIslSn4tydyUrT8ZxIwz9HhRlvYzjztLj4xfDaPECShCRbTdFJMxZncCcM7/b1dOoi7x2d2u2Tj1BQuWuEHyNgNHSFubaaBoQDkQKvxdNZ5AjIiIiIqIUNFBBMbqr54cKGg4KqEpk6mPWyL7v6qkz6jDzND0+/YeCgC+yJk6ni4zEhdumoM48Td/v9eQY5IiIiIiIKOUN1K6ewHelBdrryIW19umUA1dHThKiP6s0EBER9R+hqRC++pg2yZYNSe7fv2oSpQKPxwOHwwG32w273R7v7hANSWpIRdlGAV+jgC1TQtF0qd9H4tpxRI6IiBJXOIDQv/8c02RaeFukMA8REVGC0Rl1/Vpi4Ej6Po5IREREREREg4pBjoiIiIiIKMlwaiURDUlCCHi93nh3g45CBP0ItgRj2kweDySTGqceUSJIT0+HJA3M9ufJrH3bA4/HE+eeEFFP9eb7Gjc7IaIhqX1TACJKPtzMo3MHDhzAqFGj4t0NIuqF3nxfY5AjoiGpNyNyHo8Ho0aNwv79+/lLZC/w/es9vnexOCLXOU3TUFVVlbLvz1D4OuA1pobeXGNvvm45tZKIhiRJknr9A8Rut6fsD5/BwPev9/je0ZHIsgyXyxXvbgy4ofB1wGtMDQN9jdzshIiIiIiIKMkwyBERERERESUZBjkiom4ymUy4++67YTKZ4t2VpMT3r/f43hENja8DXmNqGKxr5GYnRERERERESYYjckREREREREmGQY6IiIiIiCjJMMgRERERERElGQY5IhrSamtrcc4558DpdCI7OxtLliyBoiidHvvXv/4VxcXFSE9Px/jx4/H444/HPP7QQw/B5XLBarVi/vz52LFjx2BcQlz11/unaRpsNhusVitsNlv05vf7B+tSBl133ztN03DPPfdg1KhRsNlsmDp1Kl599dWYY4bi5x6llk2bNuHkk09GZmYmhg8fjssuuwz19fUAgJ/97GcwmUwx3xuefPLJ6HOff/55jBs3DlarFTNnzsQXX3wRr8s4oldeeQV6vT7mOi699FIAwJdffok5c+bAZrOhsLAQTz/9dMxzk+EaX3rppZhrs9lsMBqN0Q0/kv3jWFdXh3HjxmHNmjXRtr583FRVxc0334zc3Fykp6fj7LPPxsGDB3vWKUFENITNnz9fXHzxxcLv94vy8nIxefJk8dBDD3U47s033xROp1N88cUXQtM08fnnnwun0ylee+01IYQQzz33nMjLyxPffvutaG1tFTfeeKOYPHmy0DRtsC9pUPXX+7dlyxZhNBpFMBgc7EuIm+6+d8uXLxeFhYVi165dQggh3nnnHSHLcvT+UP3co9TR0tIiRowYIX7zm9+IYDAo6uvrxcKFC8UZZ5whhBCipKREPPfcc50+d/Xq1SI9PV18+umnIhQKiT/96U8iOztb+P3+wbyEblm6dKm44oorOrQ3NjaKzMxM8ec//1mEw2Hx73//W6Snp4svv/xSCJFc13ioAwcOiBEjRogXXnhBCJHcH8dPP/1UjB07VgAQq1evFkL0/eN2zz33iGOOOUZUVFQIt9stFi9eLBYuXNijfjHIEdGQVVZWJgCIysrKaNuKFStEfn5+h2P/8pe/iAceeCCmbdGiReIXv/iFEEKI73//++L3v/999LFQKCTS09PFxx9/PEC9j7/+fP+eeeYZMXPmzIHtcALpyXunqqrw+XxCCCECgYB45plnRHp6uqiqqhJCDM3PPUot27dvF6eddppQFCXa9tZbbwm73S4CgYAwGo3i22+/7fS5F198sfjpT38a0zZhwgTxzDPPDGife+OEE04Qf/7znzu0P/XUU6KoqCim7brrrhOXXXaZECK5rrGdpmliwYIF4uqrrxZCiKT+OD733HMiPz9frFixIibI9fXj5nK5xEsvvRR9rLq6WkiSJMrLy7vdN06tJKIhq7S0FJmZmRg5cmS0bdKkSaioqEBzc3PMsddffz1uvfXW6P3a2lp88sknKCkpiZ5r6tSp0ccNBgOKioqwadOmgb2IOOrP92/9+vVobW3FrFmzMGzYMJxwwgn4/PPPB+U64qEn750sy7Barfjoo4+QlpaGq666Cvfddx9GjBgRPddQ+9yj1FJcXIz3338fOp0u2vbaa6+hpKQEmzZtQjgcxm9+8xvk5uZi/PjxePDBB6FpGoCOn/9A5Gsp0T7/NU3Dhg0b8O6776KgoAAulwvXXHMNmpqajnoNyXKNh3rxxRdRWlqKP/3pTwCQ1B/HU089FeXl5Vi8eHFMe18+bm63GwcOHIh5PDc3FxkZGdi8eXO3+8YgR0RDltfrhdVqjWlLS0sDAPh8vi6fV11djdNPPx0lJSW46KKLjniuI50n2fXn+2exWDBnzhysXLkSFRUVOOuss3Dqqadiz549A3cBcdSb9+7EE09EMBjEv/71L9x555145ZVXjniuVP7co9QlhMCdd96Jd955B48++ijcbjfmz5+PX/ziFzhw4ABefPFFLF++HMuWLQOQPJ//dXV1mD59Os4//3xs27YNn3/+OcrKynDJJZcc9RqS5RrbaZqG++67D3fccQfS09MBIKk/jsOHD4der+/Q3pePm9frBYA+XzODHBENWVarFS0tLTFt7ffbf/gcbu3atZg1axaKi4vx9ttvR7+5d3Wurs6TCvrz/Vu2bBmefvpp5OXlwWKx4KabbkJ+fj7efffdgb2IOOnNe2cymaDX6/GDH/wAl156KV5++eUjniuVP/coNXk8Hpx//vl48cUX8cknn2Dq1Kk4+eST8fHHH+PEE0+EwWDA7NmzsWTJkugfMpLl8z83NxeffPIJrrzySqSlpSE/Px8PPfQQ3n//fQghjngNyXKN7VavXo2DBw/iqquuiralysfxUEfr85Eebw9wfb1mBjkiGrKmTJmChoYG1NTURNu2bt0Kl8sFh8PR4fhnnnkGP/jBD7BkyRK8/PLL0Z242s9VWloavR8Oh1FWVoYpU6YM7EXEUX++f3fccQc2btwYc3wwGITFYhm4C4ijnrx3S5cuxdKlS2PagsEgMjMzo+caap97lHrKy8sxa9YseDwefPXVV9EpZytXrsQTTzwRc+yh3xsO//wHIl9Lifb5v3nzZtx2220QQkTbgsEgZFnG7Nmzj3gNyXKN7V5//XUsWrQoZrQpVT6Ohzpan4/0eEZGBvLy8mIer66uRmNjY8+uuccr/oiIUsi8efPEhRdeKDwej9i9e7eYPHmyuPvuuzsc99prrwmj0Sg++OCDTs/zf//3fyIvL09888030Z0Dx40bJ0Kh0ABfQXz11/t31llnieOPP14cPHhQBAIB8dvf/lYMGzZMNDQ0DPAVxE9337uVK1eKtLQ08Z///EeoqirefvttkZaWJj777DMhxND93KPU0djYKPLz88UVV1whVFWNeeyNN94QFotFrFq1KrrjbXZ2dnQnxFWrVkU39wmFQuKRRx4RGRkZCfe9Y//+/cJqtYoHH3xQhMNhsW/fPnHccceJq666StTX1wun0ykeeeQREQqFxMcffxyzYVGyXGO7qVOniqeeeiqmLVU+jjhks5O+ftzuvPNOMWXKFLF7927h8XjE4sWLxYknntiz/vTnxRERJZvq6mpx/vnni6ysLDFs2DCxdOnS6M5pVqtVvPjii0KIyA8mWZaF1WqNuV177bVCiMgOXX/84x9FYWGhsNlsYsGCBWLHjh1xu67B0l/vX0NDg7jiiitETk6OsFqtYsGCBWLTpk1xu67B0N33Tgghnn76aVFUVCTsdruYOXNmTCAeqp97lDqWLVsmAIi0tLQO3yOEEOKvf/2rGD9+vEhLSxNjxowRf/nLX2Ke/8ILL4ji4mJhtVrF7Nmzxdq1a+NxGUe1Zs0a8b3vfU+kp6eLYcOGiZ///OeitbVVCCHE+vXrxdy5c0V6eroYM2aMePbZZ2OemyzXKETk+9d7773XoT0VPo6HBjkh+vZxC4VC4tZbbxV5eXnCbreLs88+W9TU1PSoP1Jbp4iIiIiIiChJcI0cERERERFRkmGQIyIiIiIiSjIMckREREREREmGQY6IiIiIiCjJMMgRERERERElGQY5IiIiIiKiJMMgR0RERERElGQY5IiIiIiIiJIMgxwREREREVGSYZAjIiIiIiJKMgxyRERERERESYZBjoiIiIiIKMkwyBERERERESUZBjkiIiIiIqIkwyBHRERERESUZBjkiIgoqqysLN5dICIiom5gkCMiSmCjR4+G2WyGzWaLuZ1yyin9/lobN27E5MmTo/evu+46XHfddf3+Ot1x8803Iz09HVlZWWhsbOy389psNvz3v//t1XPvuecezJ8/v9/6cihJkrBmzZo+naOiogI2mw0VFRX90ykiIkpo+nh3gIiIjuyvf/0rrrjiigF/HbfbjXA4HPO68fL//t//wyuvvIJzzz23X8/r8/n69XyJJD8/P6Wvj4iIYnFEjogoic2fPx9XXHEFCgoKkJ+fD6/Xi3feeQdz585FTk4O0tLScOKJJ8ZMmXz55ZcxZcoU2Gw2TJw4Ea+++ip2796N008/HUBk1OqLL77AFVdcERMg/+///g+TJ0+G3W7H1KlT8dJLL8X04/bbb8cJJ5wQc96ubNmyBQsXLkRmZiZcLheuv/56uN1uNDQ0wGazQVEUXHzxxZ0G2KqqKpx++unR55577rk4ePAgAOCKK67AVVddhZNOOglWqxUTJ07EypUro889dORr9OjRuO666zB8+HBMnz4dmqbhmWeeQUlJCbKyspCeno4zzjgDdXV13fpYLF++HAUFBcjKysKFF16I8847D/fcc0+0z4sXL0ZhYSHS0tIwZswYPPPMM52eZ9u2bTjjjDOQn58Pi8WCSZMm4Z///CcAYMWKFTAajdi0aROAyCiqxWLBBx98gL1790KSJOzduxcAUFNTg0suuQTDhw/HyJEjcd1118Hr9QIAFEXB9ddfj+HDhyM7OxvHH388Pvvss25dJxENnvav66effhqjR4+Gw+HAySefjAMHDsS7a5QAGOSIiJLcqlWr8Pnnn2Pz5s1wu9244IILcPvtt6O2thb79++HEAL33nsvAGDNmjW48sor8dBDD8Hj8eCRRx7BJZdcgkAggPfffx9AZNTqe9/7XsxrPPfcc1i6dCkee+wxNDU14dFHH8X111+PN998M3rMk08+iUcffRSNjY0477zzcM011yAQCHTob0NDA+bPn49JkyahsrISX331FXbs2IHLLrsMWVlZ0VGl999/H88991yH599+++1wuVyoqanBtm3b4PP58MADD8T09dprr4Xb7cbtt9+OH//4x9i+fXun792XX36J7du3Y/Xq1fjqq6/w85//HP/7v/+LhoYGbNu2DWVlZVi+fPlRPwYrVqzAPffcg7///e+orq7GCSecgDfeeCP6+NVXXw2j0YjS0lJ4vV7ccMMNuOGGGzodQTvvvPMwdepUlJeXw+1249RTT8XPfvYzAMCFF16ISy65BJdddhmampqwePFi/OpXv8Jpp50Wcw5N03D22WdDlmWUlZVhy5YtqKysxDXXXAMAeOGFF/D5559j+/btqKmpwQknnIDrr7/+qNdJRPHxz3/+E9988w127tyJmpoa/O53v4t3lygBMMgRESW466+/Hk6nM+bm9/ujj59++unIy8uD0+lETk4OSktLceaZZ8Lr9WL//v3Izs5GZWUlAOD555/Hueeei4ULF0KWZZx22mn47LPPkJeXd8Q+PPPMM7j22mtx0kknQafT4aSTTsK1116LJ554InrMBRdcgOnTp8NoNOLyyy+H2+1GbW1th3O99dZbMBqNePDBB2GxWDB8+HA89thjePvtt1FdXX3U98NiseDTTz/FihUr4PV68cEHH+DRRx+NPn7GGWdg8eLF0Ov1uOyyyzBz5kysWLGi03Odf/750fd06tSpKC0txezZs9HU1ISqqioMGzYs+t4dydNPP41rr70Wc+fOhcFgwPXXX49Zs2ZFH3/qqafw+OOPw2g0oqKiAunp6Whtbe10/d+7776Le+65B5qmYe/evcjIyIjpw2OPPYZgMIjp06djxIgRuO+++zqc46uvvsLXX3+Nxx9/PLrWcNmyZVixYgUaGhpgsViwZ88ePP3009ixYwfuu+++6CgfESWeW2+9FU6nE7m5uTjzzDOxc+fOeHeJEgCDHBFRgnv88cfR3Nwcc7NardHHR44cGf2/wWDA3//+d7hcLkyaNAm//vWvUVtbC03TAAAHDx5EQUFBzPlnzZoFh8NxxD7U1NRgzJgxMW2FhYXRaXwAMHz48Jh+AIi+7uHnKigogE6nizkXgJjzdWX58uVYvHgxHn74YbhcLpSUlMRsYFJUVBRzfH5+fnTq5eEOfe90Oh0effRR5OTkoKSkBH/4wx/g8Xg6vYbD7d+/H6NHj45pO/T92r17N04++WTk5OTgJz/5CT755BMAnb8/33zzDWbOnAmXy4Wrr74amzZtghAi+rjVasWVV16Jffv24Yorroh5H9vt3bsXqqrC5XJFg+rs2bNhMpmwe/duXHjhhXjsscfw1ltvYfr06SgoKIjrmkgiOrLDv7925/sSpT4GOSKiJCdJUvT/r776Kh577DGsWbMG+/fvx3vvvYfp06dHHx81alSHXQ2XLVuGzz///IivMXr0aJSXl8e0lZeXY8SIET3u7+jRo7Fv3z6oqhpzLgDdOt+GDRtw7bXXYvPmzaipqcG8efNiNkU5fARtz549yM/P7/Rch753jzzyCD766CNs2bIFu3fvxsqVKzuE3q4UFBRg3759MW3t98PhMM444wxceumlaGhowNq1a7FkyZJOz1NVVYULLrgAf/jDH1BXV4dPPvkEF110Ucwx5eXl+N3vfoerr74aN910U6drZVwuFywWCxoaGqLhv6amBt988w2mT5+OnTt3oqSkBJ988gmam5vx+9//Hj/72c9QWlrareslIqL4Y5AjIkohbrcbOp0OFosFQgh88MEH+Nvf/oZQKAQgshnIG2+8gY8++giapuHDDz/E3XffDYfDAbPZHD3H4a6++mo88cQT+Pjjj6GqKlavXo0nn3wSV155ZY/7uHDhQkiShFtvvRWtra2orq7GL3/5S5x00kndCk6///3vccMNN8Dj8SAjIwNWqxXZ2dnRx998802sWrUKiqLgmWeewZYtWzqEoc643W4YDAYYjUYoioIXX3wRH3zwQfS9O5JrrrkGTz31FNavXw9FUfDss89i7dq1AIBQKISWlhakpaVBkiRUVFTglltuiT52KK/XC1VVoyOuW7duja5vDIVCCIfD+MlPfoILL7wQTz31FE488URceumlHf46P3v2bBQVFWHp0qXw+XxobW3Fr371K/zgBz+Aoih45513sGjRIuzduxcWiwVZWVnQ6/VHHZklIqLEwfIDREQp5PLLL8enn36KyZMnQ6/XY8KECViyZAn+/Oc/IxQK4fvf/z7+9re/4aabbsLevXtRUFCAFStWYPLkyfD7/Zg3bx5GjhyJf/zjHzHnveCCC+DxePDzn/8c+/btg8vlwsMPP4xLL720x310OBz417/+haVLl8LlckGSJJx99tl4+OGHu/X8J598Etdffz0KCwsRCoUwc+bMmP4ef/zxePDBB3HuueeiqKgI7733XnTq5pHcdNNN2LJlCwoKCmA2mzFjxgz8z//8D/79738f9bnnnXceysvLcfbZZyMQCGDhwoWYOXMmjEYjrFYrnn32Wdx11134xS9+gZycHPz0pz9FaWkptmzZgvHjx0fPU1xcjIcffhgXX3wxWlpa4HK5cM011+Dmm2/Gli1b8I9//AP19fX405/+BAB44oknMHnyZNx///24+OKLo+fR6/X45z//iZtuugnjxo1DIBDA7Nmz8a9//Qtmsxm//OUvUVlZiblz58LtdmP06NF45ZVX4HK5uvUxICKi+JPEoRPviYiIklh7uYLOdrscSJs2bYLT6YwZUSwpKcF1112Hn/70p4PaFyIiGho4tZKIiKiPPv74Y5x55pmorq6GEAKvvPIKtm7dih/+8Ifx7hoREaUoTq0kIiLqo/Ypp9OnT4fP58OECRPw9ttvd2tKJxERUW9waiUREREREVGS4dRKIiIiIiKiJMMgR0RERERElGQY5IiIiIiIiJIMgxwREREREVGSYZAjIiIiIiJKMgxyRERERERESYZBjoiIiIiIKMkwyBERERERESUZBjkiIiIiIqIk8/8BlHsUJudxVkoAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "make_plots(\n", + " df,\n", + " \"./plots/galaxies.pdf\",\n", + " n_idx=3,\n", + " intervals_xlabel=\"Fraction of spiral galaxies\",\n", + " true_theta=true_theta,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "affb506c-56ce-4d44-b8bc-62e64dbeca99", + "metadata": {}, + "source": [ + "### Calibrating `dp_alpha`\n", + "\n", + "For $n=1000$, find the largest value of the concentration parameter `dp_alpha` such that the AI posterior has frequentist coverage of 90%" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "dc025663", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Calibrating dp_alpha using 10 candidates: [ 10. 20. 40. 80. 160. 320. 640. 1280. 2560. 5120.]\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "f93e8b256aca4d8baa7a8a4e8f68a409", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Testing alphas: 0%| | 0/10 [00:00 ndarray\n" + "where y is the data array and theta is the parameter(s) to optimize.\n" + "Example: lambda y, theta: (y - theta) ** 2 # for squared loss" + ) + + # When using base_measure, X_unlabeled is not required since base_measure generates it + if Yhat_unlabeled is not None and (X is None) != (X_unlabeled is None): + raise ValueError( + "When using Yhat_unlabeled, X and X_unlabeled must both be provided or both be None." + ) + + Y = np.asarray(Y) + if Yhat_unlabeled is not None: + Yhat_unlabeled = np.asarray(Yhat_unlabeled) + if X is not None: + X = np.asarray(X) + X_unlabeled = np.asarray(X_unlabeled) + + n = len(Y) + N = len(Yhat_unlabeled) if Yhat_unlabeled is not None else None + + if X is not None and len(X) != n: + raise ValueError( + f"X must have same length as Y. Got len(X)={len(X)}, len(Y)={n}" + ) + + if ( + X_unlabeled is not None + and Yhat_unlabeled is not None + and len(X_unlabeled) != N + ): + raise ValueError( + f"X_unlabeled must have same length as Yhat_unlabeled. Got len(X_unlabeled)={len(X_unlabeled)}, len(Yhat_unlabeled)={N}" + ) + + if dp_alpha <= 0: + raise ValueError(f"dp_alpha must be positive. Got {dp_alpha}") + + if n_samples <= 0: + raise ValueError(f"n_samples must be positive. Got {n_samples}") + + if Yhat_unlabeled is not None: + if m is None: + m = N + elif m > N: + warnings.warn( + f"Truncation size m={m} exceeds unlabeled data size N={N}. Using m=N." + ) + m = N + elif m <= 0: + raise ValueError(f"m must be positive. Got {m}") + else: + # With base_measure, m is required + if m is None: + raise ValueError("m must be specified when using base_measure") + elif m <= 0: + raise ValueError(f"m must be positive. Got {m}") + + if n_jobs == -1: + import os + + n_jobs = os.cpu_count() + elif n_jobs <= 0: + raise ValueError(f"n_jobs must be positive or -1. Got {n_jobs}") + + if theta0 is None: + raise ValueError( + "theta0 (initial parameter guess) is required. " + "Please provide an appropriate initial value for optimization." + ) + + theta0 = np.asarray(theta0) + + # Default optimizer kwargs + if optimizer_kwargs is None: + optimizer_kwargs = { + "method": "Nelder-Mead", + "options": {"maxiter": 1000, "xatol": 1e-8, "fatol": 1e-8}, + } + + return ( + Y, + Yhat_unlabeled, + X, + X_unlabeled, + loss, + dp_alpha, + n_samples, + m, + n_jobs, + optimizer_kwargs, + theta0, + ) + + +def _validate_loss(loss, has_covariates=False): + """ + Validate loss function signature. + + Args: + loss (Callable): Loss function to validate + has_covariates (bool): Whether covariates are provided + + Raises: + ValueError: If loss function signature is incorrect + """ + try: + sig = inspect.signature(loss) + params = [ + p + for p in sig.parameters.values() + if p.kind in (p.POSITIONAL_ONLY, p.POSITIONAL_OR_KEYWORD) + ] + + actual_count = len(params) + + if has_covariates and actual_count == 2: + raise ValueError( + "Loss function must accept 3 parameters (y, theta, X) when covariates are provided. " + f"Got {actual_count} parameters." + ) + elif not has_covariates and actual_count == 3: + raise ValueError( + "Loss function accepts 3 parameters (y, theta, X) but no covariates were provided. " + "Either provide covariates or use a loss function with 2 parameters (y, theta)." + ) + elif actual_count not in [2, 3]: + raise ValueError( + f"Loss function must accept either 2 parameters (y, theta) or 3 parameters (y, theta, X). " + f"Got {actual_count} parameters." + ) + + except (ValueError, TypeError) as e: + if isinstance(e, ValueError): + raise + else: + warnings.warn( + "Could not validate loss function signature. " + "Potential signature errors will be caught at runtime." + ) + + +def _validate_base_measure(base_measure, Y, X=None): + """ + Validate that base_measure returns outputs compatible with provided data. + + Args: + base_measure (Callable): Function to validate + Y (ndarray): Labeled data to check compatibility against + X (ndarray, optional): Covariates to check compatibility against + + Raises: + ValueError: If base_measure returns incompatible types/shapes + """ + if base_measure is None: + return + + has_covariates = X is not None + + try: + # Test call with size=2 + test_output = base_measure(size=2) + except Exception as e: + raise ValueError(f"base_measure failed with size=2: {e}") + + if has_covariates: + # Must return (X, Y) tuple + if not isinstance(test_output, tuple) or len(test_output) != 2: + raise ValueError( + "base_measure must return (X, Y) tuple when covariates are provided. " + f"Got {type(test_output)}" + ) + + X_test, Y_test = test_output + + X_test = np.asarray(X_test) + Y_test = np.asarray(Y_test) + + if X_test.shape[0] != 2: + raise ValueError( + f"base_measure X output has wrong first dimension. Expected 2, got {X_test.shape[0]}" + ) + + if X.ndim == 1 and X_test.ndim == 2 and X_test.shape[1] != 1: + raise ValueError( + f"base_measure X output has wrong shape. Expected (2,) or (2,1), got {X_test.shape}" + ) + elif X.ndim == 2 and X_test.shape != (2, X.shape[1]): + raise ValueError( + f"base_measure X output has wrong shape. Expected (2, {X.shape[1]}), got {X_test.shape}" + ) + + else: + # Must return Y only (not a tuple) + if isinstance(test_output, tuple): + raise ValueError( + "base_measure must return Y array only when no covariates are provided, not a tuple" + ) + + Y_test = np.asarray(test_output) + + if Y_test.shape[0] != 2: + raise ValueError( + f"base_measure Y output has wrong first dimension. Expected 2, got {Y_test.shape[0]}" + ) + + if Y.ndim == 1 and Y_test.ndim == 2 and Y_test.shape[1] != 1: + raise ValueError( + f"base_measure Y output has wrong shape. Expected (2,) or (2,1), got {Y_test.shape}" + ) + elif Y.ndim == 2 and Y_test.shape != (2, Y.shape[1]): + raise ValueError( + f"base_measure Y output has wrong shape. Expected (2, {Y.shape[1]}), got {Y_test.shape}" + ) + + # Check for numeric types + if not np.issubdtype(Y_test.dtype, np.number): + raise ValueError( + f"base_measure Y output must be numeric. Got dtype {Y_test.dtype}" + ) + + if has_covariates and not np.issubdtype(X_test.dtype, np.number): + raise ValueError( + f"base_measure X output must be numeric. Got dtype {X_test.dtype}" + ) + + +def sample_ai_posterior( + Y, + loss, + dp_alpha, + theta0, + Yhat_unlabeled=None, + n_samples=1000, + m=None, + n_jobs=1, + verbose=False, + optimizer_kwargs=None, + X=None, + X_unlabeled=None, + base_measure=None, +): + """ + Samples from the posterior on the risk minimizer with an AI prior. + + The AI prior is DP(dp_alpha, F_AI) where the base measure F_AI + is the empirical distribution of Yhat_unlabeled. + + Args: + Y (ndarray): Gold-standard labels of shape (n,) or (n, d). + loss (Callable): Loss function with signature loss(y, theta) -> ndarray. + Must accept vectorized inputs where y is an array and theta is a scalar or array. + dp_alpha (float): Concentration parameter for Dirichlet Process. + theta0 (ndarray): Initial guess for optimization. Must be provided. + Yhat_unlabeled (ndarray, optional): AI predictions on unlabeled data of shape (N,) or (N, d). + Either Yhat_unlabeled or base_measure must be provided. Defaults to None. + n_samples (int): Number of posterior samples to draw. Defaults to 1000. + m (int, optional): Truncation size for F_AI (not used by default). If None, uses all unlabeled data. + n_jobs (int): Number of parallel jobs. Use -1 for all CPUs. Defaults to 1. + verbose (bool): Whether to show progress bar. Defaults to True. + optimizer_kwargs (dict, optional): Arguments for scipy.optimize.minimize. + Defaults to {"method": "Nelder-Mead"}. + X (ndarray, optional): Covariates for labeled data of shape (n,) or (n, p). Defaults to None. + X_unlabeled (ndarray, optional): Covariates for unlabeled data of shape (N,) or (N, p). Defaults to None. + base_measure (Callable, optional): Custom base measure for DP prior. If None, uses empirical distribution of Yhat_unlabeled. + When X is None: should have signature base_measure(size) -> Y array of shape (size,) or (size, d). + When X is provided: should have signature base_measure(size) -> (X, Y) tuple. Defaults to None. + + Returns: + ndarray: Array of shape (n_samples,) or (n_samples, p) containing posterior samples. + + Notes: + This implements the posterior bootstrap algorithm from the AI-priors paper, + where we sample from F|D_n ~ DP(alpha+n, G_n) by: + 1. Using all points from F_AI (no truncation by default) + 2. Sampling weights from Dirichlet distribution + 3. Minimizing weighted loss + """ + + ( + Y, + Yhat_unlabeled, + X, + X_unlabeled, + loss, + dp_alpha, + n_samples, + m, + n_jobs, + optimizer_kwargs, + theta0, + ) = _validate_args( + Y, + Yhat_unlabeled, + X, + X_unlabeled, + loss, + dp_alpha, + n_samples, + m, + n_jobs, + optimizer_kwargs, + base_measure, + theta0, + ) + + _validate_loss(loss, has_covariates=(X is not None)) + + _validate_base_measure(base_measure, Y, X) + + n = len(Y) + N = len(Yhat_unlabeled) if Yhat_unlabeled is not None else None + + if n_jobs == 1: + # sequential execution + samples = [] + for i in tqdm( + range(n_samples), desc="Posterior sampling", disable=not verbose + ): + samples.append( + _sample_one_helper( + ( + i, + Y, + Yhat_unlabeled, + n, + m, + dp_alpha, + loss, + optimizer_kwargs, + X, + X_unlabeled, + theta0, + base_measure, + ) + ) + ) + else: + # parallel execution + samples = list( + tqdm( + Parallel(return_as="generator", n_jobs=n_jobs)( + delayed(_sample_one_helper)( + ( + i, + Y, + Yhat_unlabeled, + n, + m, + dp_alpha, + loss, + optimizer_kwargs, + X, + X_unlabeled, + theta0, + base_measure, + ) + ) + for i in range(n_samples) + ), + total=n_samples, + desc="Posterior samples", + disable=not verbose, + ) + ) + + return np.array(samples) + + +def _sample_weights(n, m, alpha): + """ + Sample weights from Dirichlet distribution for n labeled and m unlabeled points. + + Args: + n (int): Number of labeled observations. + m (int): Number of unlabeled observations. + alpha (float): Concentration parameter. + + Returns: + tuple: (w_y, w_ystar) weights for labeled and unlabeled data. + """ + alpha_vec = np.array([1] * n + [alpha / m] * m) + all_weights = np.random.dirichlet(alpha_vec, size=1).flatten() + w_y = all_weights[:n] + w_ystar = all_weights[n:] + return w_y, w_ystar + + +def _sample_one_helper(args): + """Helper function for parallel sampling.""" + ( + seed, + Y, + Yhat_unlabeled, + n, + m, + dp_alpha, + loss, + optimizer_kwargs, + X, + X_unlabeled, + theta0, + base_measure, + ) = args + + if seed is not None: + np.random.seed(seed) + + # sample from base measure if provided, else use Yhat_unlabeled + if base_measure is None: + ystar = Yhat_unlabeled + xstar = X_unlabeled + else: + if X is None: + ystar = base_measure(size=m) + ystar = np.asarray(ystar) + xstar = None + else: + xstar, ystar = base_measure(size=m) + xstar = np.asarray(xstar) + ystar = np.asarray(ystar) + + w_Y, w_Y_star = _sample_weights(n, m, dp_alpha) + return _minimize_weighted_loss( + Y, ystar, w_Y, w_Y_star, loss, optimizer_kwargs, X, xstar, theta0 + ) + + +def _minimize_weighted_loss( + y, + y_star, + w_y, + w_y_star, + loss, + optimizer_kwargs, + x=None, + x_star=None, + theta0=None, +): + """ + Internal function to minimize the weighted loss criterion. + + Args: + y (ndarray): Labeled data of shape (n,) or (n, d). + y_star (ndarray): Unlabeled data of shape (m,) or (m, d). + w_y (ndarray): Weights for labeled data of shape (n,). + w_y_star (ndarray): Weights for unlabeled data of shape (m,). + loss (Callable): Loss function with signature loss(y, theta) -> ndarray. + optimizer_kwargs (dict): Arguments for scipy.optimize.minimize. + x (ndarray, optional): Covariates for labeled data of shape (n,) or (n, p). + x_star (ndarray, optional): Covariates for unlabeled data of shape (m,) or (m, p). + theta0 (ndarray): Initial guess for optimization. + + Returns: + ndarray: Optimal parameter value. + """ + + def criterion(theta): + + # call user-provided loss function with appropriate signature + if x is None: + loss_y = loss(y, theta) + loss_y_star = loss(y_star, theta) + else: + loss_y = loss(y, theta, x) + loss_y_star = loss(y_star, theta, x_star) + + weighted_loss_y = w_y * loss_y + weighted_loss_y_star = w_y_star * loss_y_star + + return np.sum(weighted_loss_y) + np.sum(weighted_loss_y_star) + + result = minimize(criterion, theta0, **optimizer_kwargs) + + if not result.success: + warnings.warn( + f"Optimization procedure did not succeed: {result.message}" + ) + + return result.x + + +def calibrate_dp_alpha_empirical_coverage_estimate( + Y, + loss, + theta0, + target_coverage=0.9, + Yhat_unlabeled=None, + m=None, + n_jobs=1, + verbose=False, + optimizer_kwargs=None, + X=None, + X_unlabeled=None, + calibration_kwargs=None, + base_measure=None, + coordinate=None, +): + """ + Calibrate dp_alpha via empirical coverage using bootstrap. + + This function tests multiple alpha values and selects the largest alpha + that achieves coverage within tolerance of the target coverage. + + Args: + Y (ndarray): Gold-standard labels of shape (n,) or (n, d). + loss (Callable): Loss function with signature loss(y, theta) -> ndarray. + theta0 (ndarray): Initial guess for optimization. Must be provided. + target_coverage (float): Target coverage level for credible intervals. Defaults to 0.9. + Yhat_unlabeled (ndarray, optional): AI predictions on unlabeled data. + m (int, optional): Truncation size for F_AI. + n_jobs (int): Number of parallel jobs. Use -1 for all CPUs. Defaults to 1. + verbose (bool): Whether to show progress. Defaults to False. + optimizer_kwargs (dict, optional): Arguments for scipy.optimize.minimize. + X (ndarray, optional): Covariates for labeled data. + X_unlabeled (ndarray, optional): Covariates for unlabeled data. + calibration_kwargs (dict, optional): Calibration parameters with keys: + - num_bootstrap_samples: Number of bootstrap samples. Defaults to 20. + - num_posterior_samples: Number of posterior samples per bootstrap. Defaults to 200. + - tolerance: Acceptable deviation from target coverage. Defaults to 0.02. + - alpha_grid: Grid of alpha values to test. Defaults to n-scaled grid. + base_measure (Callable, optional): Custom base measure for DP prior. + coordinate (int): if theta is a vector, specify which coordinate of theta to empirical validate coverage + + Returns: + float: Calibrated dp_alpha value. + """ + + if calibration_kwargs is None: + calibration_kwargs = {} + + num_bootstrap_samples = calibration_kwargs.get( + "num_bootstrap_samples", 100 + ) + num_posterior_samples = calibration_kwargs.get( + "num_posterior_samples", 100 + ) + tolerance = calibration_kwargs.get("tolerance", 0.02) + + n = len(Y) + default_grid = 0.01 * n * (2 ** np.arange(10)) + alpha_grid = calibration_kwargs.get("alpha_grid", default_grid) + + ( + Y, + Yhat_unlabeled, + X, + X_unlabeled, + loss, + _, + num_posterior_samples, + m, + n_jobs, + optimizer_kwargs, + theta0, + ) = _validate_args( + Y, + Yhat_unlabeled, + X, + X_unlabeled, + loss, + 1.0, + num_posterior_samples, + m, + n_jobs, + optimizer_kwargs, + base_measure, + theta0, + ) + + _validate_loss(loss, has_covariates=(X is not None)) + + theta0 = np.asarray(theta0) + if theta0.ndim > 0: + if coordinate is None: + raise ValueError( + "coordinate must be specified when theta0 is not scalar. " + f"theta0 has shape {theta0.shape}, please specify which coordinate to validate coverage for." + ) + if ( + not isinstance(coordinate, int) + or coordinate < 0 + or coordinate >= len(theta0) + ): + raise ValueError( + f"coordinate must be a valid index for theta0. " + f"Got coordinate={coordinate} but theta0 has length {len(theta0)}" + ) + + if verbose: + tqdm.write( + f"Calibrating dp_alpha using {len(alpha_grid)} candidates: {alpha_grid}" + ) + + best_alpha = None + alpha_pbar = tqdm( + alpha_grid, desc="Testing alphas", position=0, disable=not verbose + ) + + coverage_results = {} + + for alpha in alpha_pbar: + coverage_count = 0 + + alpha_pbar.set_description(f"Testing alpha={alpha:.2f}") + + bootstrap_pbar = tqdm( + range(num_bootstrap_samples), + desc=f"Bootstrap (alpha={alpha:.2f})", + position=1, + leave=False, + disable=not verbose, + ) + + for b in bootstrap_pbar: + # bootstrap sample from Y (and X if provided) + boot_idx = np.random.choice(n, size=n, replace=True) + Y_boot = Y[boot_idx] + X_boot = X[boot_idx] if X is not None else None + + # compute "true" parameter on bootstrap sample + def bootstrap_loss(theta): + if X_boot is None: + return np.mean(loss(Y_boot, theta)) + else: + return np.mean(loss(Y_boot, theta, X_boot)) + + result = minimize(bootstrap_loss, theta0, **optimizer_kwargs) + theta_boot = result.x + + posterior_samples = sample_ai_posterior( + Y_boot, + loss, + alpha, + theta0, + Yhat_unlabeled=Yhat_unlabeled, + n_samples=num_posterior_samples, + m=m, + n_jobs=n_jobs, + verbose=False, + optimizer_kwargs=optimizer_kwargs, + X=X, + X_unlabeled=X_unlabeled, + base_measure=base_measure, + ) + + lower_q = (1 - target_coverage) / 2 * 100 + upper_q = (1 + target_coverage) / 2 * 100 + + # extract the relevant coordinate if theta is vector + if theta0.ndim > 0: + posterior_samples_coord = posterior_samples[:, coordinate] + ci_lower = np.percentile(posterior_samples_coord, lower_q) + ci_upper = np.percentile(posterior_samples_coord, upper_q) + theta_boot_coord = theta_boot[coordinate] + else: + ci_lower = np.percentile(posterior_samples, lower_q) + ci_upper = np.percentile(posterior_samples, upper_q) + theta_boot_coord = theta_boot + + if ci_lower <= theta_boot_coord <= ci_upper: + coverage_count += 1 + + bootstrap_pbar.close() + + empirical_coverage = coverage_count / num_bootstrap_samples + + coverage_results[alpha] = empirical_coverage + + if empirical_coverage >= target_coverage - tolerance: + best_alpha = alpha + + alpha_pbar.close() + + if verbose: + print("Coverage results (dp_alpha: estimated coverage):\n") + pprint(coverage_results) + + if best_alpha is None: + if verbose: + tqdm.write( + "Warning: No alpha achieved target coverage. Using smallest alpha." + ) + best_alpha = alpha_grid[0] + + if verbose: + tqdm.write(f"\nSelected alpha={best_alpha:.2f}") + + return best_alpha + + +# Example losses + + +def squared_loss(y, theta): + """Squared loss for mean estimation""" + return (y - theta) ** 2 + + +def absolute_loss(y, theta): + """Absolute loss for median estimation""" + return np.abs(y - theta) + + +def quantile_loss(y, theta, tau=0.5): + """Check loss for quantile estimation""" + diff = y - theta + return diff * (tau - (diff < 0).astype(float)) diff --git a/tests/test_ai_priors.py b/tests/test_ai_priors.py new file mode 100644 index 0000000..ae1788d --- /dev/null +++ b/tests/test_ai_priors.py @@ -0,0 +1,253 @@ +import numpy as np +from ai_priors import ( + sample_ai_posterior, + squared_loss, + absolute_loss, + quantile_loss, +) + + +def test_posterior_consistency(): + """Test that larger vs smaller dp_alpha shifts the mean in the expected direction""" + trials = 10 + n = 100 + N = 1000 + true_mean_Y = 0.0 + true_mean_Yhat = -0.5 + dp_alphas = np.array([0.01, 1000.0]) + + posterior_means = np.zeros((trials, len(dp_alphas))) + + for i in range(trials): + Y = np.random.normal(true_mean_Y, 1, n) + Yhat_unlabeled = np.random.normal(true_mean_Yhat, 1, N) + + theta0 = np.mean(Y) + + for j, dp_alpha in enumerate(dp_alphas): + samples = sample_ai_posterior( + Y, + squared_loss, + dp_alpha, + theta0, + Yhat_unlabeled=Yhat_unlabeled, + n_samples=1000, + n_jobs=1, + verbose=False, + ) + posterior_means[i, j] = np.mean(samples) + + mean_posterior_low_alpha = np.mean(posterior_means[:, 0]) + mean_posterior_high_alpha = np.mean(posterior_means[:, 1]) + + distance_low_to_Y = np.abs(mean_posterior_low_alpha - true_mean_Y) + distance_low_to_Yhat = np.abs(mean_posterior_low_alpha - true_mean_Yhat) + distance_high_to_Y = np.abs(mean_posterior_high_alpha - true_mean_Y) + distance_high_to_Yhat = np.abs(mean_posterior_high_alpha - true_mean_Yhat) + + assert distance_low_to_Y < distance_low_to_Yhat + assert distance_high_to_Yhat < distance_high_to_Y + + +def test_different_loss_functions(): + """Test that different loss functions produce appropriate posteriors.""" + np.random.seed(42) + n = 200 + N = 2000 + + # Generate data with known median and mean with a slightly skewed distribution + Y = np.concatenate( + [np.random.normal(0, 1, n // 2), np.random.normal(2, 0.5, n // 2)] + ) + Yhat_unlabeled = np.concatenate( + [np.random.normal(0, 1.2, N // 2), np.random.normal(2, 0.6, N // 2)] + ) + + true_mean = np.mean(Y) + true_median = np.median(Y) + + # squared loss (mean estimation) + samples_mean = sample_ai_posterior( + Y, + squared_loss, + dp_alpha=10.0, + theta0=true_mean, + Yhat_unlabeled=Yhat_unlabeled, + n_samples=500, + n_jobs=1, + verbose=False, + ) + posterior_mean_squared = np.mean(samples_mean) + + # absolute loss (median estimation) + samples_median = sample_ai_posterior( + Y, + absolute_loss, + dp_alpha=10.0, + theta0=true_median, + Yhat_unlabeled=Yhat_unlabeled, + n_samples=500, + n_jobs=1, + verbose=False, + ) + posterior_median_absolute = np.median(samples_median) + + # quantile loss (75th percentile) + tau = 0.75 + + def quantile_loss_75(y, theta): + return quantile_loss(y, theta, tau=tau) + + true_quantile = np.percentile(Y, tau * 100) + samples_quantile = sample_ai_posterior( + Y, + quantile_loss_75, + dp_alpha=10.0, + theta0=true_quantile, + Yhat_unlabeled=Yhat_unlabeled, + n_samples=500, + n_jobs=1, + verbose=False, + ) + posterior_quantile = np.median(samples_quantile) + + # assert that each loss function estimates the appropriate quantity + # mean estimator should be closer to true mean than to true median + assert np.abs(posterior_mean_squared - true_mean) < np.abs( + posterior_mean_squared - true_median + ) + + # median estimator should be closer to true median than to true mean + assert np.abs(posterior_median_absolute - true_median) < np.abs( + posterior_median_absolute - true_mean + ) + + # quantile estimator should be reasonably close to true quantile + assert np.abs(posterior_quantile - true_quantile) < 0.5 + + +def test_covariates(): + """Test that covariate handling works correctly.""" + np.random.seed(123) + n = 150 + N = 1500 + p = 3 # number of covariates + + true_beta = np.array([1.0, -0.5, 2.0]) + + X = np.random.normal(0, 1, (n, p)) + X_unlabeled = np.random.normal(0, 1, (N, p)) + + noise_std = 0.5 + Y = X @ true_beta + np.random.normal(0, noise_std, n) + + ai_beta = true_beta + np.array([0.2, -0.1, 0.3]) + Yhat_unlabeled = X_unlabeled @ ai_beta + np.random.normal( + 0, noise_std * 1.2, N + ) + + # define regression loss + def regression_loss(y, theta, x): + """Squared loss for linear regression.""" + predictions = x @ theta + return (y - predictions) ** 2 + + theta0 = np.linalg.lstsq(X, Y, rcond=None)[0] + + samples = sample_ai_posterior( + Y, + regression_loss, + dp_alpha=20.0, + theta0=theta0, + Yhat_unlabeled=Yhat_unlabeled, + n_samples=500, + n_jobs=1, + verbose=False, + X=X, + X_unlabeled=X_unlabeled, + ) + + assert samples.shape == (500, p) + + posterior_mean = np.mean(samples, axis=0) + + # should be closer to true coefficients than AI coefficients for small dp_alpha + for i in range(p): + distance_to_true = np.abs(posterior_mean[i] - true_beta[i]) + distance_to_ai = np.abs(posterior_mean[i] - ai_beta[i]) + # with moderate dp_alpha, should be between true and AI + assert distance_to_true < 1.0 # reasonable bound + + +def test_base_measure(): + """Test using custom base measure instead of Yhat_unlabeled.""" + np.random.seed(456) + n = 100 + m = 1000 + + Y = np.random.normal(1.0, 1.0, n) + + def example_base_measure(size): + # 70% from N(0, 1), 30% from N(3, 0.5) + n1 = int(0.7 * size) + n2 = size - n1 + samples = np.concatenate( + [np.random.normal(0, 1, n1), np.random.normal(3, 0.5, n2)] + ) + np.random.shuffle(samples) + return samples + + theta0_simple = np.mean(Y) + + samples_base = sample_ai_posterior( + Y, + squared_loss, + dp_alpha=10.0, + theta0=theta0_simple, + base_measure=example_base_measure, + m=m, + n_samples=500, + n_jobs=1, + verbose=False, + ) + + assert samples_base.shape == (500,) or samples_base.shape == (500, 1) + assert np.all(np.isfinite(samples_base)) + + # test with covariates + p = 2 + X = np.random.normal(0, 1, (n, p)) + true_beta = np.array([1.5, -1.0]) + Y_cov = X @ true_beta + np.random.normal(0, 0.5, n) + + # base measure that returns (X, Y) pairs + def example_base_measure_covariates(size): + X_sample = np.random.normal(0, 1, (size, p)) + + base_beta = np.array([1.0, -0.5]) + Y_sample = X_sample @ base_beta + np.random.normal(0, 0.7, size) + return X_sample, Y_sample + + def regression_loss(y, theta, x): + return (y - x @ theta) ** 2 + + theta0 = np.linalg.lstsq(X, Y_cov, rcond=None)[0] + + samples_base_cov = sample_ai_posterior( + Y_cov, + regression_loss, + dp_alpha=15.0, + theta0=theta0, + base_measure=example_base_measure_covariates, + m=m, + n_samples=500, + n_jobs=1, + verbose=False, + X=X, + ) + + assert samples_base_cov.shape == (500, p) + assert np.all(np.isfinite(samples_base_cov)) + + posterior_mean = np.mean(samples_base_cov, axis=0) + assert np.all(np.abs(posterior_mean - true_beta) < 1.0) From 198ea318bfdf02e5d763533098a7aa911c9a027f Mon Sep 17 00:00:00 2001 From: seanohagan Date: Tue, 5 Aug 2025 14:36:51 -0500 Subject: [PATCH 2/2] fix import in test_ai_priors.py --- tests/test_ai_priors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_ai_priors.py b/tests/test_ai_priors.py index ae1788d..a4dd17c 100644 --- a/tests/test_ai_priors.py +++ b/tests/test_ai_priors.py @@ -1,5 +1,5 @@ import numpy as np -from ai_priors import ( +from ppi_py import ( sample_ai_posterior, squared_loss, absolute_loss,