select_on_test.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Description**: Simulation check for math in [my answer on stats.stackexchange.com](https://stats.stackexchange.com/a/570680) that there's a way to select a model w/ the best test set performance while also guaranteeing *something* about the selection procedure's performance. Also gives us a sense of how loose the theoretical upper bound is, at least for a particular linear problem.\n",
                "\n",
                "**Estimated runtime**: 5 min.\n",
                "\n",
                "**Required packages**:\n",
                "```\n",
                "matplotlib==3.4.3\n",
                "pandas==1.3.4\n",
                "scikit-learn==0.24.2\n",
                "tqdm==4.62.3\n",
                "```"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "metadata": {},
            "outputs": [],
            "source": [
                "from __future__ import annotations # in case you're on <= python3.8\n",
                "import warnings\n",
                "\n",
                "import matplotlib.pyplot as plt\n",
                "import numpy as np\n",
                "import pandas as pd\n",
                "from sklearn.datasets import make_classification\n",
                "from sklearn.linear_model import LogisticRegression\n",
                "from sklearn.model_selection import train_test_split\n",
                "from tqdm.notebook import tqdm"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "metadata": {},
            "outputs": [],
            "source": [
                "# globally configure figure and font sizes to be kinda big\n",
                "plt.rcParams['figure.figsize'] = (10, int(4*(10/6)))\n",
                "plt.rcParams['font.size'] = 12"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "metadata": {},
            "outputs": [],
            "source": [
                "def hoeffding_bound(num_models: int,            # m in math\n",
                "                    max_acceptable_diff: float, # eps in math\n",
                "                    test_size: int              # n in math\n",
                "                    ) -> float:\n",
                "    return 2*num_models * np.exp(-2 * (max_acceptable_diff**2) * test_size)\n",
                "\n",
                "\n",
                "def simulate(test_size: int,\n",
                "             num_models: int=2,\n",
                "             num_simulations: int=20,\n",
                "             train_size: int=1_000,\n",
                "             test_size_huge: int=50_000,\n",
                "             max_acceptable_diff: float=0.01,\n",
                "             noise: float=0.1,\n",
                "             random_state: int=None):\n",
                "    '''\n",
                "    Simulates the following model selection process `num_simulations` times:\n",
                "      1. `num_models` models are trained and tuned on a training set with `train_size` examples\n",
                "      2. Each of these models is evaluated on the test set with `test_size` examples\n",
                "      3. The model with the highest test set performance is selected\n",
                "\n",
                "    and returns:\n",
                "      * observed fraction of `num_simulations` in which the selected model's\n",
                "        test set error is outside the interval\\n\n",
                "            `(true_error - max_acceptable_diff, true_error + max_acceptable_diff)`\n",
                "        where `true_error` is the selected model's \"true\" out-of-sample error.\n",
                "      * theoretical upper bound on the above fraction.\n",
                "\n",
                "    Implementation notes:\n",
                "      * The \"true\" out-of-sample error is actually just an estimate on a test\n",
                "        set with `test_size_huge` examples.\n",
                "      * the classification data is generated via a sparse linear model with some\n",
                "        flipped labels to allow for greater variance in model errors.\n",
                "    '''\n",
                "    # input checks\n",
                "    if num_models < 1:\n",
                "        raise ValueError(f'num_models must be at least 1. Got {num_models}.')\n",
                "    if max_acceptable_diff < 0 or max_acceptable_diff > 1:\n",
                "        raise ValueError('max_acceptable_diff must be in the interval [0,1]. '\n",
                "                         f'Got {max_acceptable_diff}.')\n",
                "    # before looking into any data, compute the upper bound\n",
                "    frac_oob_ub = hoeffding_bound(num_models, max_acceptable_diff, test_size)\n",
                "    # simulation\n",
                "    num_oob_observed: int = 0\n",
                "    total_size = train_size + test_size + test_size_huge\n",
                "    for i in range(num_simulations):\n",
                "        _random_state_data = None if random_state is None else i+random_state\n",
                "        X, y = make_classification(total_size, flip_y=noise,\n",
                "                                   random_state=_random_state_data)\n",
                "\n",
                "        X_tr, X_te_huge, y_tr, y_te_huge = train_test_split(X, y,\n",
                "                                                            test_size=test_size_huge)\n",
                "        X_tr, X_te, y_tr, y_te = train_test_split(X_tr, y_tr,\n",
                "                                                  test_size=test_size)\n",
                "\n",
                "        # [10^-4, 10^4] is a fine range for make_classification data\n",
                "        Cs = np.logspace(start=-4, stop=4, num=num_models)\n",
                "\n",
                "        best_accuracy_obs = -np.inf\n",
                "        best_model = None\n",
                "        for j, C in enumerate(Cs):\n",
                "            with warnings.catch_warnings():\n",
                "                warnings.filterwarnings('ignore',\n",
                "                                        message='The max_iter was reached which '\n",
                "                                                'means the coef_ did not converge')\n",
                "                _random_state_model = None if random_state is None else j+i+random_state\n",
                "                model: LogisticRegression = (LogisticRegression(C=C, penalty='l1',\n",
                "                                                                solver='liblinear',\n",
                "                                                                random_state=_random_state_model)\n",
                "                                            .fit(X_tr, y_tr))\n",
                "            accuracy = model.score(X_te, y_te)\n",
                "            if accuracy > best_accuracy_obs:\n",
                "                best_accuracy_obs = accuracy\n",
                "                best_model = model\n",
                "        best_accuracy_almost_true = best_model.score(X_te_huge, y_te_huge)\n",
                "        is_oob = np.abs(best_accuracy_obs - best_accuracy_almost_true) > max_acceptable_diff\n",
                "        num_oob_observed += is_oob\n",
                "    return num_oob_observed/num_simulations, frac_oob_ub"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 4,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "fc674801d8a948f3b7df189e6f2b1aef",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Simulating:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "202902851141475cb25df5f7a57c6695",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=2000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "55385a3f91a949819ea426237cd3f80b",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=4000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "aea47ec800d34ed489ad6595daa62c01",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=6000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "9be9d494b0f9439ab73589769c636282",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=8000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "a8e9d926af1d42bba1d885135be30a04",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=10000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "503b4c6e98a34b3682fc4413497ee676",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=12000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "5f0c976d253e46199bb4c89a98167d04",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=14000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "b33508301d214277b7179a9e50ecf19e",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=16000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "6cfe6a3ea59349888de3aa5e4d99f4f9",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=18000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                },
                {
                    "data": {
                        "application/vnd.jupyter.widget-view+json": {
                            "model_id": "2fbfae517f484d06b0a20a0c2f1cbc76",
                            "version_major": 2,
                            "version_minor": 0
                        },
                        "text/plain": [
                            "Replicating n=20000, m=2:   0%|          | 0/10 [00:00<?, ?it/s]"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "num_models = 2\n",
                "\n",
                "random_state = 43045\n",
                "num_replications = 10\n",
                "test_sizes = np.linspace(2_000, 20_000, num=10).astype(int)\n",
                "_records = []\n",
                "for test_size in tqdm(test_sizes, desc='Simulating', position=0):\n",
                "    for _ in tqdm(range(num_replications), leave=False, position=0,\n",
                "                  desc=f'Replicating n={test_size}, m={num_models}'):\n",
                "        frac_oob_observed, frac_oob_ub = simulate(test_size, num_models,\n",
                "                                                  random_state=random_state)\n",
                "        _records.append({'test set size': test_size,\n",
                "                         'frac oob (observed)': frac_oob_observed,\n",
                "                         'frac oob (upper bound)': frac_oob_ub})\n",
                "sim_results_df = pd.DataFrame(_records)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 5,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAGKCAYAAABJgObQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABNaUlEQVR4nO3dd3xUVd7H8c8vISEhoSdSQwCRqhSJilhAFOviY1kboGvX3VXXtmuDVextdd3y7K6uPlasqyKKoLgiIoICikoV6U1DJ4QQkpznjzsJQ0yZwMzcmcn3/Xrd18ztvzMJmR/n3HOOOecQERERkdiQ5HcAIiIiIrKHkjMRERGRGKLkTERERCSGKDkTERERiSFKzkRERERiiJIzERERkRii5ExEJMDMOpqZM7Oj63DO4MA57SMZm4jUH0rORCRqzOxHMzss8P5TMxvud0yJwMyOMbP/mNlqM9tpZt+b2V1m1tDv2ESk7hr4HYCI1A9m1gXIAL4ys1QgD/jM36gSxlHAD8ATwCrgUOCfQCvg1z7GJSL7QDVnIhItRwEznXMlwGHARufcitpOCjQZXmtmr5rZDjNbaWa/NLOmZvaSmW03s6Vmdnal87qZ2XtmVhBYxgcSxOBjzjWzJWZWZGbTgd5V3L9LoFZqi5ltNrMPzOyQGuJNMbPHArVYu8xsnZm9EuqHtC+ccw865/7gnJvqnFvmnPsP8CBwbiTvKyKRoeRMRCIqkNRswavJOTrw/gOgddC+2twBTAD6AO8CzwOvAB8C/YD3gOfNrGXgnumBe6QBgwJLJjAxUGuHmfULXOP1wHUfxat5Co69FTAN+Ak4BhgALAKmmFl2NbFei5cUjQQOAk4HZtRUODN7PyiJrG45JoTPKVhTYEMdzxGRGGCaW1NEIsnMOgbezgB+A8zBS4peBsYBOOeW13C+A55wzl0fWM/GS5b+5py7NrCtObAJGOace9fMLgP+CnRwzm0IHNMKWA5c5Zx73sxeBDo75wYG3euawHnHOOemmdldwMnOuQFBxxiwBPirc+7PZjYY+BjIcc6tNrMngEOA412If2DNrB2QXstha5xzO0O8Xg9gJnC7c+5voZwjIrFDz5yJSEQ555abWW8gBS8Za4RXU3W6c+6nEC8zN+h6+WZWCnwTtG2zmRUDBwQ29QLmlydmgWN+NLNFgX0APYGPKt1nWqX1w4D+ZlZQaXs6Xq1YVf4Pr0ZviZl9GHg/3jlXXF3hnHNrqttXV2Z2EF6t4StKzETik5IzEYkYM5sH5OL9rUkBtuI9TpEGLPUqoejpnFtZy6V2h7DNsfejGlXVWlnQdqvmmGBJeAncNVXs21rVCc65r82sEzAUOA6vqfQeMxvgnNtW1Tlm9j5es2lNTnHOfVrTAWZ2MF4yOA51BBCJW0rORCSSTsVLyp4BJgKvAX8EivEeWAdYG4H7zgOuNrOsSs2aXfGeLSs/5qhK51VenwVcTB2aFAGccwXAW8BbZnY/sA7vubfx1ZxyOSE0a9a0MzBEyUTgReD6UJtURST2KDkTkYhxzq0ws2S8XpBXOeeWBGp37nbOLYngrcfiJYGvmtnv8WrJHsVLcF4NHPM48KWZ3Qc8h9fceVOl6/wNuAx428zuxRumoj1wCvCec2565RsH7rcW+BooBC4ASoHF1QW7v82aZnYsXkeJN4AHgFaBWkmcc+v359oiEn3qrSkikdYPKHbOLTCzpngPy0+N5A0DtVwnArsC9/oE2IH3cH9x4JjZwHDgfOBb4FbghkrX+RE4Eq/X45t4PTVfwmuqXVfN7bcBNwKfB657JnC2c25R+Er4M5cCjYFLAnEFLyISZ9RbU0RERCSGqOZMREREJIYoORMRERGJIUrORERERGKIkjMRERGRGKLkTERERCSGJMw4Z1lZWa5jx45+hyEiIiJSq9mzZ29wzmVXtS9hkrOOHTsya9Ysv8MQERERqZWZrahun5o1RURERGKIkjMRERGRGKLkTERERCSGKDkTERERiSEJ0yFARERkf+3evZvVq1dTVFTkdyiSANLS0mjfvj0pKSl1Ok/JmYiISMDq1atp3LgxHTt2xMz8DkfimHOOjRs3snr1ajp16lSnc9WsKSIiElBUVETLli2VmMl+MzNatmy5T7WwSs5ERESCKDGTcNnX3yUlZyIiIiIxRMmZiIhIHFq+fDlmRklJid+hVCkzM5OlS5fu1zUuvvhiRo0aVeW+Z599lqOPPnq/rh8uZsaSJUvCdj0lZ6Eaaz9fREREoqRjx45MnjzZ7zCqNHjwYP7973/vta2goIDOnTv7FFF8U3IWiuoSMSVoIiJSbtlL8HZHGJvkvS57ye+IwiJWa+YSmZIzERGR/bXsJfjiSihcATjv9Ysrw5agXXjhhaxcuZJhw4aRmZnJww8/XLHvpZdeokOHDmRlZXHfffdVbC8rK+PBBx/kwAMPpGXLlpx77rls2rSpYv8777xDr169aNasGYMHD2bBggUV+zp27MhDDz1E7969ycjIoKSkhBkzZjBw4ECaNWtGnz59mDJlCgB33HEHn376Kddccw2ZmZlcc801wN5NfTt37uSmm24iNzeXpk2bcvTRR7Nz504AzjnnHFq3bk3Tpk059thjmTdvXsifi3OOa6+9lqZNm9K9e3c++uijin1r167l9NNPp0WLFnTp0oWnnnqqYl/l5tIpU6bQvn37vcr/6KOP0rt3b5o2bcp55523V6/LRx55hDZt2tC2bVueeeaZkOMNlcY5ExERqc7s62Hz17Uft2EGlO3ae1tpIcy8DH54qupzyjXvC/3/XOMhL7zwAp9++in//ve/OeGEEwDvmTOAadOmsWjRIhYvXszhhx/OWWedRY8ePfjLX/7C22+/zSeffEJ2djbXXXcdv/3tb3n55ZdZvHgxF1xwAW+//TaDBw/m8ccfZ9iwYcyfP5/U1FQAXn75Zd577z2ysrL48ccfOe2003jhhRc4+eST+eijjzj77LNZuHAh9913H5999hkjR47k8ssvrzL+m2++mXnz5jF9+nRat27NzJkzSUry6odOOeUUnnnmGVJTU7nlllsYMWIEX3/9dc2fWcDMmTP55S9/yYYNG3jzzTc566yzWLZsGS1atOCCCy6gV69erF27loULFzJ06FA6d+7M8ccfH9K1X3vtNSZOnEhaWhpHHXUUzz77LFdffTUTJ07k0Ucf5aOPPqJTp05cccUVIV2vLlRzJiIisr8qJ2a1bQ+jO++8k/T0dPr06UOfPn2YO3cuAP/617+47777aN++PQ0bNuSuu+7ijTfeoKSkhFdffZXTTjuNoUOHkpKSws0338zOnTuZPn16xXWvu+46cnJySE9P58UXX+TUU0/l1FNPJSkpiaFDh5KXl8eECRNqja+srIxnnnmGJ554gnbt2pGcnMzAgQNp2LAhAJdeeimNGzeuiHHu3Lls3bo1pLIfcMABXH/99aSkpHDeeefRrVs33nvvPVatWsW0adN46KGHSEtLo2/fvlx++eW88MILIX+u1113HW3btqVFixYMGzasImF87bXXuOSSSzj44IPJyMjgrrvuCvmaoVLNmYiISHVqqdGq8HbHQJNmJY1y4YQpYQzo51q3br3ndo0aUVBQAMCKFSs488wzK2qoAJKTk/nxxx9Zu3Ytubm5FduTkpLIyclhzZo1FdtycnIq3q9YsYLXX3+d8ePHV2zbvXs3xx13XK3xbdiwgaKiIg488MCf7SstLeWOO+7g9ddfJz8/vyLWDRs20LRp01qv3a5du73GEsvNzWXt2rWsXbuWFi1a0Lhx4732zZo1q9Zrlqv8ua5duxbwmkv79++/13XDTTVnoRju6rZdRETqlz73QXKjvbclN/K2h0ldBzTNycnh/fffZ8uWLRVLUVER7dq1o23btqxYsSeZdM6xatUq2rVrV+X9cnJyuPDCC/e61o4dO7j11ltrjS0rK4u0tDR++OGHn+0bO3Ys48aNY/LkyWzdurWiqda50L5f16xZs9exK1eupG3btrRt25ZNmzaxffv2vfaVly8jI4PCwsKKfevXrw/pfgBt2rRh1apVe1033JSchWq427P0vMXbtjH0DFxERBJYpxFw+JNeTRnmvR7+pLc9TFq1alWnccOuvvpq7rjjjookLD8/n3HjxgFw7rnn8t577/HRRx+xe/du/vSnP9GwYUMGDhxY5bVGjhzJ+PHjmTRpEqWlpRQVFTFlyhRWr15da2xJSUlceuml3Hjjjaxdu5bS0lI+//xzdu3axfbt22nYsCEtW7aksLCQ22+/vS4fCT/99BN/+ctf2L17N6+//joLFizg1FNPJScnh4EDB3LbbbdRVFTEN998w9NPP82IEd7Po2/fvkyYMIFNmzaxfv16/vznP4d8z3PPPZdnn32W+fPnU1hYyJgxY+oUcyiUnO2LXrdD2gEw53oIMbsXEZEE12kEnLEchpd5r2FMzABuu+027r33Xpo1a8ajjz5a6/G/+93vOP300znxxBNp3LgxAwYMYObMmQB069aNF198kWuvvZasrCzGjx/P+PHjKzoDVJaTk8O4ceO4//77yc7OJicnh0ceeYSysrKKe73xxhs0b96c66677mfnP/rooxxyyCEcdthhtGjRgltuuYWysjIuuugicnNzadeuHT179mTAgAF1+kyOOOIIvv/+e7Kysrjjjjt44403aNmyJeB1aFi+fDlt27blzDPPZMyYMQwdOhTwer/26dOHjh07cuKJJ3LeeeeFfM9TTjmF66+/niFDhtClSxeGDBlSp5hDYaFWHca6vLw8V5e25P225Cmvm/RRr0LuudG7r4iIRMyCBQvo0aOH32FIAqnud8rMZjvn8qo6RzVn+6rzpdCsN3z9Byit+4zzIiIiIlVRcravkpK9Xjw7VsDCx/2ORkRERBKEkrP90eo4aH8GzLsfdq7zOxoRERFJAErO9le/R7xBBueOqv1YERERkVooOdtfjbtAt9/B0v+DTXP8jkZERETinJKzcOg1ChpmwZwbNLSGiIiI7BclZ+GQ2hR63wM/TYVVb/odjYiIiMQxJWfhcuBl0OwQ+Or3GlpDRERE9pmSs3BJagCHPg47lsGiJ/yORkRE6pGLL76YUaPiq2PaXXfdxciRI/0OIyZFJTkzs4Zm9rSZrTCz7Wb2lZmdUs2xF5tZqZkVBC2DoxHnfmt9PLQbBt/dBztDn0RVREREpFy0as4aAKuAQUBTYDTwmpl1rOb4z51zmUHLlOiEGQb9HoXSnfDNaL8jERERCZuSkhK/Q6g3opKcOed2OOfucs4td86VOefeBZYB/aNx/6hq0hW6Xgs/PA2bv/Y7GhERSSALFixg8ODBNGvWjF69evHOO+9U7NuwYQNDhw6lcePGDBo0iBUrVgDgnOOGG27ggAMOoGnTpvTu3ZvvvvsOgF27dnHzzTfToUMHWrVqxdVXX83OnTsBmDJlCu3bt+ehhx6idevWXHLJJfTo0YN333234p4lJSVkZWUxZ443lNSMGTMYOHAgzZo1o0+fPkyZMqXi2GXLljFo0CAaN27M0KFD2bBhQ6Q/rrjlyzNnZtYK6ArMq+aQfma2wcwWm9loM2sQxfD23yGjoWELmK2hNURE6pOyMkf+9l24CPzt3717N8OGDePEE0/kp59+4q9//SsjRoxg0aJFALz00kuMHj2aDRs20LdvX0aMGAHABx98wNSpU1m8eDFbtmzh1VdfpWXLlgDccsstLF68mK+//polS5awZs0a7r777op7rl+/nk2bNrFixQqefPJJLrjgAl5++eWK/ZMmTSIrK4tDDz2UNWvWcNpppzFq1Cg2bdrEo48+ytlnn01+fj4Aw4cPp3///mzYsIHRo0fz3HPPhf0zShjOuaguQAowGfhXNfs7A53wEsdDgPnAbdUceyUwC5jVoUMHF1MW/d25l3Bu5Vt+RyIiIiGaP3/+Pp9bWlrmzv3ndHfgbe+5c/853ZWWloUxMuemTp3qWrVq5UpLSyu2nX/++e7OO+90v/rVr9x5551XsX379u0uKSnJrVy50n300UfuoIMOcp9//vle55aVlblGjRq5JUuWVGybPn2669ixo3POuY8//tilpKS4nTt3Vuz//vvvXWZmptuxY4dzzrnhw4e7MWPGOOece/DBB93IkSP3ivnEE090zz77rFuxYoVLTk52BQUFFfsuuOACN2LEiHB8NDGtut8pYJarJleKas2ZmSUBLwDFwDVVHeOcW+qcW+a85s9vgbuBX1Zz7JPOuTznXF52dnbE4t4nXa6Epj3hq5uhdJff0YiISIRt3FHM7BWbKSlzzF6xmY07isN6/bVr15KTk0NS0p6v7tzcXNasWQNATk5OxfbMzExatGjB2rVrGTJkCNdccw2//e1vadWqFVdeeSXbtm0jPz+fwsJC+vfvT7NmzWjWrBknn3xyRU0XQHZ2NmlpaRXrXbp0oUePHowfP57CwkLeeecdhg8fDsCKFSt4/fXXK67VrFkzpk2bxrp161i7di3NmzcnIyNjr9ilalFLzszMgKeBVsDZzrndIZ7qAItYYJFSPrRGwQ+w+K9+RyMiIhGWlZlK/9zmNEgy+uc2JyszNazXb9u2LatWraKsrKxi28qVK2nXrh0Aq1atqtheUFDApk2baNu2LQDXXXcds2fPZt68eSxevJhHHnmErKws0tPTmTdvHlu2bGHLli1s3bqVgoKCiut4X917K2/aHDduHD179qRLly6AlxxeeOGFFdfasmULO3bs4NZbb6VNmzZs3ryZHTt27BW7VC2aNWf/AHoAw5xzO6s7yMxOCTyThpl1x+vZOS46IYZZmxOh7Wnw3T1Q9JPf0YiISASZGS9fMYDPbzueV64cUGVisz+OOOIIMjIyePjhh9m9ezdTpkxh/PjxnH/++QBMmDCBadOmUVxczOjRozniiCPIycnhyy+/ZObMmezevZuMjAzS0tJITk4mKSmJK664ghtuuIGffvK+o9asWcOkSZNqjOP888/ngw8+4B//+EdFrRnAyJEjGT9+PJMmTaK0tJSioiKmTJnC6tWryc3NJS8vjzvvvJPi4mKmTZvG+PHjw/r5JJJojXOWC1wF9AXWB41fNsLMOgTedwgcfjzwjZntACYAbwL3RyPOiOj3KJQUwjd/9DsSERGJsKQkI7txw7AnZgCpqam88847vP/++2RlZfGb3/yG559/nu7duwPeA/djxoyhRYsWzJ49m5deegmAbdu2ccUVV9C8eXNyc3Np2bIlN998MwAPPfQQXbp0YcCAATRp0oQTTjihooNBddq0acORRx7J9OnTOe+88yq25+TkMG7cOO6//36ys7PJycnhkUceqajpGzt2LDNnzqRFixaMGTOGiy66KOyfUaIwlyC9CfPy8tysWbP8DqNqs6/3mjZP/gqa9/Y7GhERqcaCBQvo0aOH32FIAqnud8rMZjvn8qo6R9M3RcPBf4SUZjDnRg2tISIiIjVSchYNDVvAIXfBjx/BGrWxi4iISPWUnEXLQVdDk+4w5yYoDW/3ahEREUkcSs6iJSkFDn0MCpbA93/3OxoRERGJUUrOoqntKdDmZPh2DBRpTjERkViUKB3lxH/7+ruk5CzaDv0TlBTAt3f6HYmIiFSSlpbGxo0blaDJfnPOsXHjxr1mWAhVfE0ongia9oSDfg3f/6/32uxgvyMSEZGA9u3bs3r16r2mMBLZV2lpabRv377O5yk588Mhd8GyF72hNY6bBBEYrFBEROouJSWFTp06+R2G1HNq1vRDw5Zegrb+Q1g7we9oREREJIYoOfNL199Ak25e7VlZqHPAi4iISKJTcuaXpBTo9yfYvhgW/6/f0YiIiEiMUHLmp7anQusT4du7YNdGv6MRERGRGKDkzE9m3sC0Jdu8BE1ERETqPSVnfmvWC7pcBd//A7Yu8DsaERER8ZmSs1hwyBhokOnNuykiIiL1mpKzWJCWDQf/Eda9D2vf9zsaERER8ZGSs1jR9RrI7OLVnmloDRERkXpLyVmsSE715t3ctgC+/5ff0YiIiIhPlJzFknbDoNXx3qTouzb5HY2IiIj4QMlZLCkfWmP3Fvjubr+jERERER8oOYs1zXvDgVfA4r/D1oV+RyMiIiJRpuQsFvW+Gxo0gq9u9jsSERERiTIlZ7Eo7QA4eDSsfQ/WfeB3NCIiIhJFSs5iVddrIfNAmHMjlJX4HY2IiIhEiZKzWJXcEPo9AlvnwQ9P+R2NiIiIRImSs1jW/gw4YDB8MxqKt/gcjIiIiESDkrNYZgb9H/fGPPvuHr+jERERkShQchbrmveFAy+DRX+BbYv9jkZEREQiTMlZPOh9LySnw1e/9zsSERERiTAlZ/EgvRUcfAeseQfWT/Y7GhEREYkgJWfxotvvIKMTzL5BQ2uIiIgkMCVn8SI5LTC0xnfww9N+RyMiIiIRouQsnuScBQccGxhaY6vf0YiIiEgEKDmLJ2Zw6OOwawPMu8/vaERERCQClJzFmxaHQudLYNGfYfsSv6MRERGRMFNyFo/63AtJqfDVH/yORERERMJMyVk8Sm8DvW6H1W/Bjx/7HY2IiIiEkZKzeNXtBsjIDQytUep3NCIiIhImSs7iVYN06PswbJkLS//P72hEREQkTJScxbMO50D2UfDNHbB7m9/RiIiISBgoOYtnZnDon6HoJ5h3v9/RiIiISBgoOYt3LfOg069g4eNQsNTvaERERGQ/KTlLBH3uB2ugoTVEREQSgJKzRNCoLfS6DVb9B378xO9oREREZD8oOUsU3W+CRjkwR0NriIiIxDMlZ4mifGiNzV/Bsuf9jkZERET2UVSSMzNraGZPm9kKM9tuZl+Z2Sk1HH+Dma03s61m9oyZNYxGnHEv9zzIOhLm3g67t/sdjYiIiOyDaNWcNQBWAYOApsBo4DUz61j5QDM7CbgVOB7oCHQGxkQpzvhmBoc+DkXrYf6DfkcjIiIi+yAqyZlzbodz7i7n3HLnXJlz7l1gGdC/isN/BTztnJvnnNsM3ANcHI04E0LWEdBxJCz4ExQs9zsaERERqSNfnjkzs1ZAV2BeFbt7AXOD1ucCrcysZTRiSwh9HwBLgq9v8TsSERERqaOoJ2dmlgK8BDznnFtYxSGZwNag9fL3jau41pVmNsvMZuXn54c/2HjVqD30vAVWvgY/TfM7GhEREamDqCZnZpYEvAAUA9dUc1gB0CRovfz9z55wd8496ZzLc87lZWdnhzXWuNfj916SNud6cGV+RyMiIiIhilpyZmYGPA20As52zu2u5tB5QJ+g9T7Aj865jREOMbE0aAR9HoRNs2HZC35HIyIiIiGKZs3ZP4AewDDn3M4ajnseuMzMeppZc2AU8GwU4ks8HS+AlkfA3Ntgd4Hf0YiIiEgIojXOWS5wFdAXWG9mBYFlhJl1CLzvAOCcmwg8DHwMrAgsd0YjzoRjSdD/z7BzHSx42O9oREREJAQNonET59wKwGo4JLPS8Y8Bj0U0qPoiawDkDocFj8CBl0NGB78jEhERkRrsU82ZmaWbWWq4g5EI6fsgYPD1rX5HIiIiIrUIKTkzs0fN7PDA+9OATcAWMxsWyeAkTDJyvN6bK16G/Ol+RyMiIiI1CLXmbATwXeD9H4GRwOnA/ZEISiKg5x8gvS3MuUFDa4iIiMSwUJOzRs65wsAo/Z2dc/9xzk0GciMYm4RTgwzo8wBs/AKWj/U7GhEREalGqMnZYjMbgTdw7IcAZpYF1DQkhsSaTiOhRZ737FnJDr+jERERkSqEmpz9BvgtcBwwOrDtJOCDSAQlEVIxtMYamP+I39GIiIhIFUJKzpxzXzrnBjrnBjvnfghse8k5d2Fkw5Owyz4KOpznjXu2Y5Xf0YiIiEgl1Y5zZmZDQrmAc+6/4QtHoqLfQ7D6bW/mgIEv+h2NiIiIBKlpENqnK623AxywEWiJN6jsaqBzZEKTiMnIhR43w7z7oOu1kHWE3xGJiIhIQLXNms65TuUL8BTwV6C5c64t0Bz4S2C7xKOet0Jaa5h9PTjndzQiIiISEGqHgBuAW51zhQCB19uAGyMVmERYSib0fQA2zoAVr/gdjYiIiASEmpztAA6vtO0woDC84UhUdboImh8KX98CJfpRioiIxIJQk7PRwEQzG2tmD5nZWGAiMCpyoUnElQ+tUbgKFvzJ72hERESE0IfSeAE4AlgANAEWAgMC2yWeHXAM5PwS5j8IhWv8jkZERKTeq6m35l6cc/OB+RGMRfzS72FY9Qa83X7v7cPVUUBERCTaQkrOzKwFcDPQF8gM3uecOzb8YUlUvVPNaChjTQmaiIhIlIVaczYWaAi8hjoBiIiIiERMqMnZQCDbObcrksGIiIiI1Heh9tb8Bmhf61EiIiIisl9CrTn7L95QGv8HrA/e4Zx7JuxRiYiIiNRToSZnx+DNozm00nYHKDmLd8Od9/B/ZZYEa96Fdr+IfkwiIiL1VEjJmXPuuEgHIj6r3CuzZAdMHgSfnQ9DP4PmffyJS0REpJ4J6ZkzM0uqbol0gOKTBhlw7DuQ2hw++QXsXOd3RCIiIvVCqMlVCbC7mkUSVaO2MGg8FG+GT4Z5tWkiIiISUaEmZ52AzkHLUcB44MoIxSWxonlfGPgybJoDn18ErszviERERBJaqHNrrqi0zAB+BdwS2fAkJrQfBoc+BqvehLm3+x2NiIhIQgt5bs0qNAGywxWIxLhuv4Pti2H+Q9D4IDjwMr8jEhERSUihzq35At6wGeUaAccCL0YiKIlBZtD/L7D9B/jiasjoBK2H+B2ViIhIwgn1mbMlwA9BywxguHPu2kgFJjEoqQEc/Ro06Qqfng3bFvkdkYiISMIJdZyzMZEOROJEalMY9B58cARMOQ1OnAFpWX5HJSIikjBCHqfMzC4xs/+a2aLA6yWRDExiWGZHOOZtKFwNn54Jpbv8jkhERCRhhDoI7R3ArcArwHWB1z8Etkt9lH0kHPkc5E+DmVeAc7WfIyIiIrUKtbfm5cBg59yK8g1mNgmYCtwXicAkDuSeB9u/h29Ge8+hHTzK74hERETiXqjJWQaQX2nbRiA9vOFI3Ol1B2xb7CVojQ/yEjYRERHZZ6E+czYReMnMuplZupl1B54DJkUuNIkLZnDEU5B9NHz+K8j/3O+IRERE4lqoydk1wHZgLrAj6FVDaQgkN4Rj3oJG7WHq/0DBMr8jEhERiVuhTt+0zTl3Ed7gs62BdOfcRc65LZEMTuJIWhYMfg/KdsMnv4DirX5HJCIiEpfqMpTGQcBtwJ3AbYF1kT2adINj/uM9gzbtXCgr8TsiERGRuBPqUBrDga+A3njNmYcAcwLbRfZoPQQO/yes/wBmXashNkREROoo1N6a9wKnOuemlm8ws2OAF4CxkQhM4tiBl3m1Zwse9mrTul/vd0QiIiJxI9TkrDFQuRveDLwhNkR+ru8DULAE5twImQdC+2F+RyQiIhIXQn3m7DHgfjNLAzCzdLzBZx+LVGAS5ywJjnwBWvSH6RfA5q/9jkhERCQumKvmmSAzWwWU7zS8XpoO2Aw0D2xb55zrEIU4a5WXl+dmzZrldxhS2c51MOlw79mzk76ARm39jkhERMR3ZjbbOZdX1b6amjVHRigeqU/S28Cgd+HDo2Hq6XDCJ9BAreEiIiLVqTY5c859Es1AJIE17wNHveIlZ9NHesNtWMijuIiIiNQrdf6GNLNtkQhEEly706DfY7D6bfj6Vr+jERERiVn7Un1h+3IjM7vGzGaZ2S4ze7aG4y42s1IzKwhaBu/LPSXGdLsODvoNLHgElvzb72hERERiUrXNmmbW3Dm3OYz3Wos3XtpJQHotx37unDs6jPeWWGAG/Z+AgqXw5a8hsxO0Pt7vqERERGJKTTVnK8rfmNnkoO099+VGzrk3nXNvAxv35XxJEEkNvOfPmnSDT8+GrQv9jkhERCSm1JScFZrZwWaWDBxuniRgjZkllS8RiqufmW0ws8VmNtrMQh0sV+JBalOvB2dyQ/jkNCja4HdEIiIiMaOm5GoM8AVQjDcTQAmwO2gpXw+3qcDBwAHA2cAFwO+rOtDMrgw8xzYrPz8/AqFIxGR2hGPHwc618OkZULrL74hERERiQrXJmXPuH0ATIBfYCXQGDgy8dgY6BV7Dyjm31Dm3zDlX5pz7Frgb+GU1xz7pnMtzzuVlZ2eHOxSJtKwBMOA5yP8MZl6uSdJFRESoZW5N51wJsNrM+jnnVtR0bAQ59rGHqMSB3HNh+/fwzSho3BUOGe13RCIiIr4K9Zmx5WY2xsyWmVmRmS0NrKeGeiMzaxCYmzMZSDaztKqeJTOzU8ysVeB9d2A0MC7U+0gc6nU7dLoIvv0jLH/Z72hERER8FWpy9jBwAnAV0Ae4GhgCPFSHe43Cax69FW9qqJ3AKDPrEBjLrHyOzuOBb8xsBzABeBO4vw73kXhjBoc/CdnHwIxLIH+63xGJiIj4ptqJz/c6yGw10Mc5tzFoWxYw1znXLoLxhUwTnyeAXRth0gDYvRVOmumNgyYiIpKAapr4PNSas+qe+dKzYBI+DVvC4HehbDdMOQ2Kt/gdkYiISNSFmpy9Dow3s5PMrIeZnQy8DbwWscikfmrSDY590+skMO1cL1ETERGpR0JNzv4ATAb+DswG/gp8DNwSobikPmt1nPcM2voPYda1GmJDRETqlZBG3nfOFQN/DCwikXfgJbB9Ecx/yKtN636D3xGJiIhERZ2nXzKzWyMRiMjP9Lkfcs6COTfB6nf8jkZERCQq9mVuzNvDHoVIVSwJjnwBWvSH6cNh01d+RyQiIhJx+5KcqYemRE+DRjDoHUhtAZ8Mg8I1fkckIiISUfuSnL0Y9ihEapLeBgaN98Y/+2QYlOzwOyIREZGICSk5M7Oby987534dtP3GSAQl8jPN+8BRr8CWuTB9BJSV+h2RiIhIRIRac1ZdL81R4QpEpFbtToNDH4fV42Cu+qWIiEhiqnEoDTMbEnibbGbHsffzZp2B7ZEKTKRKXa+FbYthwaPQuCt0ucLviERERMKqtnHOng68pgHPBG13wHrg2kgEJVItM+j/Zyj4Ab78jTf/ZusT/I5KREQkbGpMzpxznQDM7Hnn3EXRCUmkFkkN4OhX4YOj4NNfwomfQ9MefkclIiISFiE9c6bETGJOShNvkvTkht4k6UX5fkckIiISFqH21lxlZiurWiIdoEi1MnLh2HegaB18eiaUFvkdkYiIyH4LaW5NYGSl9TbA74BXwhuOSB1lHQFHPg/TzoUZl8HAF73n0kREROJUqBOff1J5m5lNASYCT4Q5JpG66XAO9LkP5t4BTbrCIXf6HZGIiMg+C7XmrCq7gE7hCkRkv/S8zRti49u7oPFB0HG43xGJiIjsk5CSMzO7u9KmRsCpwPthj0hkX5jB4U/CjmUw4xLvebTso/yOSkREpM5CnSEgp9KSBjwG/CpCcYnUXXIqHPMmNOoAU8+AgqV+RyQiIlJnoT5zdkmkAxEJi4YtYfB78MEAmPILOHE6pDbzOyoREZGQhVpzhpkdZ2bPmNmkwOuQ2s8S8UGTrl4NWsESmHYOlO32OyIREZGQhTrO2eXAq3hTNr0JrAPGmpkmNpTY1GowHPYvWD8ZZl0DzvkdkYiISEhC7a35B2Coc25u+QYzexX4D/BUJAIT2W8HXgLbv4f5D0DjbtDjRr8jEhERqVWoyVlLYH6lbYuAFuENRyTM+tzrJWhf3eQtwYarNk1ERGJPqM+cTQMeM7NGAGaWATwCTI9UYCJhYUmw6o2q943VTAIiIhJ7Qk3OrgZ6A1vN7EdgC9AHuCpCcYmIiIjUS6EOpbEOGGRm7YG2wFrn3OqIRiYiIiJSD9Vp+qZAQqakTBJH6S5Ibuh3FCIiIhVCHudMJCFNPhZ2rPI7ChERkQpKziTxVdcr85g3YesCmNgffvw4ujGJiIhUo07NmiJxq7oErUkP+PRM+O9Q6PsQdL/Rm0RdRETEJyEnZ2bWFOgGZAZvd879N9xBiURN0+5w0hcw4xL46mbY+AUc8TSkZNZ+roiISASElJyZ2cXA34ECoDBolwM6hz8skShKaQxHvw4LHoa5t8PWeXDMW9DkIL8jExGReijUZ87uA37pnGvlnOsUtCgxk8RgBj1vgeMmQdF6mJQHq8f7HZWIiNRDoSZnDYAPIhmISExofQKcPBsyu8DU0+GbO8GV+R2ViIjUI6EmZw8Bo8xMvTsl8WXkwtBp0Pli+O5u+GQYFG/2OyoREaknQk22bgBGAdvNbGXwEsHYRPzTIB2OeAYO+wes/xAm5sHmuX5HJSIi9UCovTVHRjQKkVhkBgddDc37wqdnwwdHwuFPQacRfkcmIiIJLNS5NT+JdCAiMStrgPcc2mfnwecjYdOX0O8RSErxOzIREUlAITVrmlmKmY0xs6VmVhR4HWNmqZEOUCQmpLeGIZOh2/Ww6An46HjYud7vqEREJAGF+szZw8AJwNVAn8DrELyOAiL1Q1IK9H8cBo6FTbO8aZ/yP/c7KhERSTChJmfnAKc75z5wzi1yzn0AnAmcG7nQRGJUxwvgxBmQnAYfDYLv/wGumumhRERE6ijU5Ky6yQY1CaHUT817w8mzoPVQ+PI3MPNSKNnpd1QiIpIAQk3OXgfGm9lJZtbDzE4G3gZei1hkIrEutTkMGg8H3wlLn4UPj4aC5X5HJSIicS7U5OwPwGS8+TVnA38DPgZuiVBcIvHBkqD3XV6SVvCD9xzaug/9jkpEROJYSMmZc67YOfdH51wX51yjwOto59yuUG9kZteY2Swz22Vmz9Zy7A1mtt7MtprZM2bWMNT7iPii3S/gpC8hvS1MORnmPajn0EREZJ9UO86ZmR3rnJsaeD+kuuOcc/8N8V5rgXuBk4D0Gu57EnArXm/QtcBbwJjANpHY1eQgOGkGzLwc5t7mjYc24P8gpYnfkYmISBypaRDa/wUODrx/uppjHNA5lBs5594EMLM8oH0Nh/4KeNo5Ny9w/D3ASyg5k3jQIMMbaqPl4fDV72HSEXDMW9C0u9+RiYhInKi2WdM5d3DQ+07VLCElZnXUCwiexHAu0MrMWkbgXiLhZwbdb/AGrd21ESYdBqve9DsqERGJE6HOEDCumu2R+MbJBLYGrZe/b1zF/a8MPMc2Kz8/PwKhiOyHVoPhlDnQtKc3N+fXt0FZqd9RiYhIjAu1t+Zx1WwfHKY4ghUAwQ/plL/fXvlA59yTzrk851xednZ2BEIR2U+N2sMJU6HLVTD/QZhyChRt8DsqERGJYTVOfG5mdwfepga9L9cZWBGBmObhTRFVPoZaH+BH59zGCNxLJPKSG8Lh/4SWh3kD1k7Kg2PehBaH+h2ZiIjEoNpqznICS1LQ+xy8B/pX4U3rFBIza2BmaUAykGxmaWZWVXL4PHCZmfU0s+bAKODZUO8jErMOvAyGTgNXBh8M9AauFRERqaTGmjPn3CUAZjbdOffUft5rFHBn0PpIYIyZPQPMB3o651Y65yaa2cN4g9ymA/+pdJ5I/Gp5GJw8Gz47H2ZcAhu/gEP/DMmpfkcmIiIxwlwIA2Wa2UXA1865b4K29QF6O+deiGB8IcvLy3OzZs3yOwyR0JSVwNzbYcEjkHUkHP06NGrnd1QiIhIlZjbbOZdX1b5QOwTcg9eMGWwV3qCyIlJXSQ2g38Nw9Guw5Rtv2qefpvodlYiIxIBQk7MmwLZK27YCzcIajUh90+EcOHEmpDSFj4bAwic07ZOISD0XanI2Hzi70rYzgQXhDUekHmrWC076wpufc871MH0klBT6HZWIiPikxg4BQW4BJpjZecAPQBfgeODUSAUmUq+kNvWG15j3AHwzGrZ+5603PtDvyEREJMpCqjlzzk3Dm2fzSyAD+AI42Dn3WQRjE6lfLAkOvgMGT4DCVTAxD9ZM8DsqERGJslCbNQkMc/Ggc+63gdfKHQREJBzangwnz4LMjvDJL+Dbe7yx0UREpF4ItVkTMzsdGARkAVa+3Tl3UQTiEqnfMjvD0M/gi6vg2z/Cpi/hyOchtZnfkYmISISFOvH5ncC/AsefA2wETgK2RCwykfquQSMvIev/V1j7Pkw8DLZ853dUIiISYaE2a14KDHXO3QAUB16HAR0jFZiIAGbQ7Ro4/mMoKYBJR8CKV/2OSkREIijU5KyZc678v+zFZpbinPsCr5lTRCLtgKO9aZ+a9/WmfppzszfLgIiIJJxQk7MfzKxX4P13wK/N7EJgc2TCEpGfadTWq0E76Lew8E/w36FQ9JPfUYmISJiFmpyNAloG3t8KXAc8AtwYiaBEpBrJqXDY32DAc7Bxhjft04aZfkclIiJhVGtvTTNLAoqAGQCB5swuEY5LRGrS+SJodgh8ehZMPhby/gZfXPnz44ZrKigRkXhTa82Zc64MGOecK45CPCISqhb9vPHQDhhcdWIGMNaq3i4iIjEr1GbNqWY2IKKRiEjdNWzpzSggIiIJI9RBaFcA75vZOGAVUNFW4pz7YyQCE5EQJSX7HYGIiIRRqMlZOvB24H37oO16oEUk1hXlQ1q231GIiEiIqk3OzOwa59zfAqv3OeeWRCkmEQmndzpD95ugx42Q0sTvaEREpBY1PXN2X9D7OZEORET2Q3W9Mk+bD21Ogu/GeEnagsegtCi6sYmISJ3U1Kz5g5n9CZgHpJjZpVUd5Jx7JiKRiUjdVJegHfMGbPwS5t4OX90Eix6HQ+6CTr+CpFCfbBARkWgx56r+g25mXYE/ALnAccCnVRzmnHNDIhde6PLy8tysWbP8DkMktq3/L8y9DTZ+AU26Qe97IOdssFA7bouISDiY2WznXF5V+6r9b7NzbjFweeACHznnjo9QfCISLa2HQKsZsHocfHMHTDsXWvSHPvdD66HeROsiIuKrkP67rMRMJIGYQc4ZcMo3MOBZ2LUBPj4JPhoCG2b4HZ2ISL2ntgyR+iopGTr/Cn6xCPr/BbbOgw+OhKlnwJbv/I5ORKTeUnImUt8lN4Ru18LpS71n0H78GCb0hukXQcEyv6MTEal3lJyJiCclEw4e5SVpPW6CVa/Du91g1rWw80e/oxMRqTeUnInI3hq2hH6PwLAl0PkS+P4fMP5AmDsKirf4HZ2ISMJTciYiVWvUDg7/lzeQbdtfwLz7vIFs5z8MJYV+RycikrCUnIlIzZp0haNfgZPnQMsB8PUtMP4g+P5fULbb7+hERBKOkjMRCU2LfnDcBDjhE8joCF9eDe/2hOUvgyvzOzoRkYSh5ExE6uaAY2HoNBg0Hhqkw/Th8P6hsGYCVDPjiIiIhE7JmYjUnRm0+wWc8jUc+SKUbIdPToPJx8JP0/yOTkQkrik5E5F9Z0nQaQSctgDy/g7bl8DkY2DKL2DzXL+jExGJS0rORGT/JadC19/A6UugzwOQ/xm83w8+G+4lbCIiEjIlZyISPg0yoNet8D9LoectsPpteLcHfPFrKFzrd3QiInFByZmIhF9qc+j7AJz+A3S5En74N4zvAl/fCsWb/Y5ORCSmKTkTkchJbwOH/R1+sRByzvIGsB3XCebdDyU7/I5ORCQmKTkTkchrfCAMfBFOnesNxTH3DnjnQFj8dygt9js6EZGYouRMRKKn2SEw6B0Y+hk06QazroF3u8OyF6Gs1O/oRERigpIzEYm+7IFw/BQY/D6kNoPPL4T3+8LqdzSQrYjUe0rORMQfZtD2ZDh5Fhz1CpTtgqn/Ax8eBT9+4nd0IiK+UXImIv6yJMg9D06bB4c/CTtWwkeD4eOTYdMcv6MTEYm6Bn4HICICQFIKdLkCOo6E7/8O8x6Aif2hwznQ+x7v2bTKhqsJVEQSj2rORCS2NEiHHjfD6Uvh4NGwdkLViRnAWItubCIiUaDkTERiU2pT6H23l6SJiNQjSs5EJLalHVDz/q3zoxOHiEiURC05M7MWZvaWme0wsxVmNrya4y42s1IzKwhaBkcrThGJM+/18iZZX/AoFK72OxoRkf0WzZqzvwPFQCtgBPAPM+tVzbGfO+cyg5Yp0QpSROJM/ycgKRW++j283QEmHwdL/q05PEUkbkUlOTOzDOBsYLRzrsA5Nw14B7gwGvcXkThXXa/M4Q66XQcnzYRfLIZD7oKda+GLK+DN1jD1TFj5BpQWRTVcEZH9Ea2hNLoCpc65xUHb5gKDqjm+n5ltADYBLwAPOOdKIhyjiMSy2obNaHIQHPJHr4fnptmwfCysfAVWvw0pTbyJ1zuOgAOOg6TkqIQsIrIvopWcZQJbK23bCjSu4tipwMHACqAX8CpQAjxQ+UAzuxK4EqBDhw5hDFdE4pYZtMzzln6PwE8fe4naqv/A0mchrTXknu8lai36e8eLiMQQc1GYx87M+gGfOecaBW27CRjsnBtWy7nnA793zvWv6bi8vDw3a9assMQrIgmoZCesfc9L1Na+B2XF0LgrdBzuJWqNu/gdoYjUI2Y22zmXV9W+aHUIWAw0MLODgrb1AeaFcK4D9F9bEdk/DdKhwy/h2DfhrPVw+FOQ3ha+HQPjD4JJR8DCJ2Dner8jFZF6LirJmXNuB/AmcLeZZZjZUcD/4D1PthczO8XMWgXedwdGA+OiEaeI1BOpzaHL5XDCx3DGSq/5s6wY5lwPb7eD/54ES5+D3dv8jlRE6qFoDqXxGyAd+Al4Gfi1c26emXUIjGVW/tDY8cA3ZrYDmICX1N0fxThFpD5p1N6bLuqUr7zJ13veBtu/hxkXw5utYNp5sHoclBb7HamI1BNReeYsGvTMmYiEjXOwYQYsfwlWvgq7Nni1bR3O8Z5Pyz4aTBOsiMi+q+mZMyVnIiI1KdsN6yd7idrqt6FkBzTKgdwLvESt2SHq8SkidVZTchatoTREROJTUgq0PcVbSnZ4TZzLx8LCx2DBw9C0l5ek5V4AmR39jlZEEoBqzkRE9kVRPqx8HVaMhfzPvG3ZR3mJWs45kJblb3wiEtPUrCkiEkkFy2DFy17T59b5YA2gzUleotb+dGiQ4XeEIhJjlJyJiESDc7DlGy9JW/EyFK72ErP2Z3iJWusTvGZSEan39MyZiEg0mEHzPt7S90H46dNAj8/XvdeG2dDhXC9RyxqgjgQiUiXVnImIRFrpLlg30UvQ1oyH0iLI7Ay5w73po5r2gLFVJGq1TfYuInFLzZoiIrFi9zZY9ZaXqP34Ebiymo9XgiaSkGJhbk0REQFIaQKdfwVDPoAzVsOhj/sdkYjEGCVnIiJ+SW8D3a+v+ZhJA2DOzd74ars2RiUsEfGXOgSIiMSypAaw+K+w8E/eetOekH2MtxxwNGTk+hufiISdkjMRkVg2dJrXgWDjF5A/zesBuuJlWPIvb3+jnECiFkjYmvbQvJ8icU7JmYiI34a7mntrJqfBAcd6Sy+grNQbTy1/GuR/Cj/+15upACC1hTdTQXnC1vxQSE6NWlFEZP+pt6aISLxzDgp+8GrVyhO27d97+5LToeURe2rWso6ElEx/4xURDUIbLiUlZSzJL6Brq0ySktRsICIxwgwad/GWAy/xtu1cv6cZNH8azLvPG7bDkqF5P8g+OpCwHQ1pB/gbv4jsRTVnISopKaPfvR+yvaiExmkN+GrUUBo0UIImInFi9zbI/9yrVcufBhtnes+yATTpFuhkEEjYMjpp9gKRCFPNWRgsyS9ge1EJANuLSliSX0D3Nk18jkpEJEQpTaDtSd4C3qwFm2Z7ydpPn8LKN+CHf3v70tsGdTI4Gpodok4GIlGk5CxEXVtl0jitQUXNWddWemZDROJYckPIHugtPW/xmjy3zgs0g5YnbK96x6Y03buTQYs873wRiQg1a9aBnjkTkXrDOdixYk+ilv8pbFvo7UtqCC0P39PJIHugVzMnIiHT3JoiIrL/ivIh/7M9CdvmOeBKvSbPZr33Hm8tvfXPz9fk7iIVlJyFSXFxKbNXbubwTs1JTk6O6L1ERGLe7gKvY0F5zdqGGVBa6O3L7OLNYFA+m8G7Xau/jhI0qYfUISAMiotL6XHnREodJBssGHMyqalK0ESkHkvJhNbHewtA2W7YNGfPWGtrxsPSZ30NUSQeKTkL0eyVmykN/Oeu1HnrR3bJ8jcoEZFYkpQCWUd4S4+bvE4G2xZ6ydoXV1V/3swroEkPb+qpJt29+ULVO1TqMSVnITq8U3OSjYqas8M7Nfc7JBGR2GZJ3kTtTXvWnJytfht2/XvPenK6N/Zak0Cy1rSH977xQeolKvWCkrMQJScns2DMyXrmTEQk3M7Oh6INXi3btgXe69YFsOFzWPEKEGi2sCTI6LwnWSuvaWvSA1Kb+loEkXBSclYHSUlG84xUTCNni4jUTW2Tu6dlQdrRXieCYCWFsH2xl6xtWxB4XQjrJkFZ8Z7j0tvsqWmrSNx6eNv1N1vijJKzEGn6JhGR/bQvvTIbNILmfb0lWFkJFCzbU9NWnrgtf9GbqqpcSpNKCVvgfWZnSNJXoMQm/WaGSNM3iYjEkKQG0OQgb+H0Pdudg6L1e2rayptI138Iy54LOj/Ve4atck1bk67QICPqxREJpuQsRJq+SUQkDph5TZnpbaD1kL33FW+FbYsCSVugpm3Lt16HBFe657iM3ECi1gOadt/zPq2WHvoaZFfCRIPQ1oGmbxIRSUClu2D7kr1r2srfl+7cc1zDrCpq2rpDRgd4uYZOYkrQpAoahDZMGjRIUlOmiEiiSW4IzXp5SzBXBoWrft5Euvot+OGpoPMbRTdeSXhKzuog0WvOEr18IiJ1YkleE2dGLrQ9ee99wUN/bF0Aix6v/jr/yYZG7SG9HTRqB+ntA6+B9UbtIaWpepVKBSVnIUr03pqJXj4RkbCqPPRHTclZztlQuAZ2roGNX8Cu/J8fk9xoT6JWZRLXHtJaQZLG2KwPlJyFKNF7ayZ6+UREfHP4P/deL90FO9d5yVrh6j2JW+Fq7zV/mvdatnvv8ywZ0lrXksS184Yfkbim5CxEid5bM9HLJyISUbUNshssuSFkdvSW6rgy2LXBS9zKk7bCNbAzkMxtWwjrJ+89plu51OZ7atuCk7bgJK5hy7o3o6o3atSot2YdJPozWYlePhGRhLO7oPoauPLEruhHKqbAKpfUcO8m0yqTuDbeZPZQdWJWTgnaPlFvzTBJ9N6aiV6+sjLHxh3FZGVqCi4RSRApmZDSzZsovjplu2Hn+qpr4HaugY0zYdUaKNtV6UTznnNr1K7mGIryvdo6zbgQNvokpV4oK3Nc8NQMZq/YTP/c5rx8xQCSkpSgiUg9kJQCGTneUh3noHhT9c2oNXnzAO81pSmktoCGLSC1ZeC1hdeEmtqi6n1K6qqkT0TqhY07ipm9YjMlZY7ZKzazcUcx2Y0b+h2WiEhsMPOSqIYtoXnvn++vqVmz/1+8xG7XJijeGHjdBAVLvdfizfysWTVYRVLXMoTkrvx9s8gkdTHyXJ2SM6kXsjJT6Z/bvKLmLCsz1e+QREQSQ7dra97vymD3VtgVlLjt2hhI3ILelyd3ISd1zYJq4EJM7lKbVz8cSXUJ6FiLeoKm5EzqBTPj5SsG6JkzEZF9UZfeqJVZkpcUpTaHxnW4Z1mpl9RVVStXVXJXsNQ7pngLoSV1lRK4GKLkTOqNpCRTU6aIyL6KdvNeUrKXPDVssR9J3cbak7vtS7zXGKLkrA4SvbdfwY5iXpm1iouOzCE1NfGa/YqKSvh4cT5De2bToEFi/uoXFu5m/LfrOOvQNqSkpPgdTtgVF5cye+VmDu/UnORkjZQuEm1x8T24V1LXJfTzanquLsoS8xsqAhK9t1/BjmIOvudDAO59fyHfjR5KZkbiJGhFRSV0v2tSxfrCu04iLS2xfv0LC3fT8+4PALjlrW+Z/8cTadQocRK04uJSetw5kVIHyQYLxpxMaqoSNJFoSfTvwViikUZDVFVvv0TyyqxVNa7Hu48X59e4ngjGf7uuxvV4N3vlZkoDrSqlzlsXkehJ9O/BapttfeitqeQsROW9/RokWUL29rvoyJwa1+Pd0J7ZNa4ngrMObVPjerw7vFNzkgP/SU82b11EoifRvwcBLxGrvPggatM3mVkL4GngRGADcJtzbmw1x94A3AKkA/8Bfu2cqzx08V6iMX1TXLS17wc9cxb/9MyZiERSon8PRlNN0zdFMzl7Ga+m7jKgL/AeMNA5N6/ScScBzwNDgLXAW8AM59ytNV0/GsmZiIiISDjUlJxFpVnTzDKAs4HRzrkC59w04B3gwioO/xXwtHNunnNuM3APcHE04qzvSkrKWLhuG2VlZX6HEhGFhbt5deZKdu/e7XcoEZPoZdxWsIs/f7CIXbtqrEiPWxu2FnLjq19RWFjodygRk+hlTPTybdq2kzve/IadO3f6HUrExMLPMCo1Z2bWD5junEsP2nYzMMg5N6zSsXOB+51zrwbWs4B8IMs5t7G6e6jmbP+UlJTR794P2V5UQuO0Bnw1aigNGiTOI4nBPRmBhOvJCIlfxm0Fu+h97+SK9W9GnUCTzMQZt27D1kLyHvi4Yn3WbceR1bSRjxGFX6KXMdHLt2nbTg69/78V63NuH0KLJuk1nBF/ovkz9L3mDMgEtlbatpWqh5WrfGz5+58da2ZXmtksM5uVn594ve+iaUl+AduLSgDYXlTCkvwCnyMKr0TvyQiJX8Znpi+vcT3e3T9xUY3riSDRy5jo5fvT5O9rXE8EsfIzjFZyVgA0qbStCbA9hGPL3//sWOfck865POdcXnZ24vW+i6aurTJpHBj3q3FaA7q2yvQ5ovBK9J6MkPhl/PWgjjWux7t7h3WrcT0RJHoZE718o045qMb1RBArP8NoNWtmAJuBXs657wPbngfWVn7Q38zGAsucc3cE1ocAY51zrWu6h5o1919JSRlL8gvo2iqTpKTEadIsl+g9GSHxy7itYBfPTF/Orwd1pGHDxGnSLLdhayH3T1zEvcO60ahR4jSHBUv0MiZ6+TZt28mfJn/PqFMOIj09sZo0y0XrZxgrvTVfwZuJ9HK83poTqLq35snAs3i9NdfhDaXxhXprioiISKKIhWfOAH6DN27ZT8DLeGOXzTOzDmZWYGYdAJxzE4GHgY+BFYHlzijGKSIiIuKbqI3E6ZzbBJxRxfaVeJ0Agrc9BjwWnchEREREYkfiPVgkIiIiEseUnImIiIjEECVnIiIiIjFEyZmIiIhIDFFyJiIiIhJDojbOWaSZWT7esBuRlgVsiMJ9/KLyxb9EL6PKF/8SvYwqX/yLRhlznXNVTm+UMMlZtJjZrOoGjUsEKl/8S/QyqnzxL9HLqPLFP7/LqGZNERERkRii5ExEREQkhig5q7sn/Q4gwlS++JfoZVT54l+il1Hli3++llHPnImIiIjEENWciYiIiMQQJWciIiIiMaReJmdm1tDMnjazFWa23cy+MrNTgvYfb2YLzazQzD42s9ygfWZmD5nZxsDysJlZ0P6OgXMKA9c4IdrlC2ZmB5lZkZm9GLQtkcp3vpktMLMdZvaDmR0T2B73ZQzEMcHMNpvZejP7m5k1iNfymdk1ZjbLzHaZ2bOV9kWsPGY2PPBvfYeZvW1mLaJdRjMbYGYfmtkmM8s3s9fNrE28lbGmn2HQMXeamQuOMRHKZ2aNzOx/zWyDmW01s6nxVr4QyniueX9Pt5vZfDM7I97KaD5+v4e1jM65ercAGcBdQEe8BPUXwPbAehawFTgHSAMeAWYEnXsVsAhoD7QD5gNXB+3/HHgMSAfOBrYA2T6W9QPgU+DFwHrClA8Yijfw8IDAz7FdYEmIMgITgGcDZWgNfAtcF6/lA84CzgD+ATwbtD1i5QF64f3bPhbIBMYCr/hQxlMC5WsCNAKeASbGWxmrK1/Q/gMDv6drgRMSqXzAi8ArQDaQDPSPt/LV8jvaDigO/K4acBpQCBwQT2XEp+/3cJcxIj/8eFyAbwIf9pXA9Eo/6J1A98D6dODKoP2Xlf9wga7ALqBx0P5Pg3+4US7T+cBrgV/U8uQskco3Hbisiu0JUUZgAXBq0PojwL/ivXzAvez9pRCx8gD3A2OD9h2I9wXUOJxlqq2MVew/FNhe6Xc5bspYXfmA94FTgeXsnZzFdfmAbsA2oEk1x8dV+aop4xHAT5WOyQeOjNcyBt0v4t/v4S5jvWzWrMzMWuF98PPwst+55fucczuAHwLbqbw/8D5431Ln3PZq9keNmTUB7gZuqrQrUcqXDOQB2Wa2xMxWm9fsl06ClBF4Ajg/0JzSDu9/tBNJnPKVi2R5Kl/7B7w/mF3DGP++OBbv7025uC+jmZ0DFDvnJlSxO97LdwReLf2YQLPmt2Z2dtD+eC8fwCxggZmdbmbJgSbNXXiJzc/iJE7KGMXv97CWsd4nZ2aWArwEPOecW4hXHbm10mFbgcaB95X3bwUyA+3StZ0bTfcATzvnVlXanijlawWkAL8EjgH6Av2AUSROGT/B+we/DViN98fzbRKnfOUiWZ6YK6+Z9Qb+CPw+aHNcl9HMMvFqDq6v5pC4Lh9eM9fBgfu2Ba4BnjOzHoH98V4+nHOlwPN4zXG7Aq9XBRIYiMMyRvn7PaxlrNfJmZklAS/gZbfXBDYX4D0XEqwJXltyVfubAAXOq8es7dyoMLO+wAnA41XsjvvyBewMvP7VObfOObcB71mAU0mAMgZ+NycBb+JVvWcBzYGHQogx5stXSSTLE1PlNbMueE1/v3POfRq0K97LOAZ4wTm3rJr98V6+ncBu4F7nXLFz7hPgY+DEwP54Lx+Bh9sfBgYDqcAg4N+B7xOIszL68P0e1jLW2+QskAk/jVcDc7Zzbndg1zygT9BxGXhtx/Oq2h94H7yvs5k1rmZ/tAzGe/hxpZmtB24GzjazOSRG+XDObcarTXJV7E6EMrYAcoC/Oed2Oec2Av+Hl3wmQvmCRbI8la/dGWgILA5j/CEJ9AqbDNzjnHuh0u54L+PxwHXm9Spej/e7+5qZ3VJVjMRf+b6pZX+8lw+81oepzrlZzrky59yXwEy8/+j/LE5iuIw+fb+Ht4zReBgvFhfgn8AMILPS9my8qsiz8XpzPMTevTmuxntQux1e9fY89u7NMQN4NHDumfjQ0w+vN1jroOVR4I1A2eK+fEGx3A18CRyAV6v0KV5zbkKUEVgK3Ao0AJoBb+FV0cdl+QLlSAMewPsfbVpgW8TKw55m4WPwaiBfJLI94aorYzu8Z1t+X815cVHGGsrXkr3/5qzC6xGXmSDlSwGWAKMD60fh1Yh0j6fy1VLGQcAGoG/guH7ARuDEOCxj1L/fw13GiHwwsb4AuXg1LkV4VZHly4jA/hOAhXhV2VOAjkHnGl7V76bA8jCBabAC+zsGztmJ1yX3hGiWrZry3kWgt2YilQ/vD+b/Bv6BrAf+AqQlShnx/ic7BdiM90fzdfZ0a4+78gV+D12l5a5IlwcYDqwEdgDjgBbRLiNwZ+B98N+bgngrY00/w0rHLWfv3ppxXz68L9/PAzHMB86Mt/KFUMZr8JLQ7Xj/Obwp3sqIj9/v4Syj5tYUERERiSH19pkzERERkVik5ExEREQkhig5ExEREYkhSs5EREREYoiSMxEREZEYouRMREREJIYoORMRiTIzKwiMIC4i8jNKzkQk5pjZ8sBcf/t7nYvNbFo4Yqri2i4wV2adOecynXNLwx2TiCQGJWciIiIiMUTJmYjEFDN7AegAjA80//0hsH2AmU03sy1mNtfMBgedc7GZLTWz7Wa2zMxGmFkPvDn2jgxcZ0s19/vZuUH7LjWzBWa22cwmBSYvx8ymBg6ZG7j2eVVct4uZfWJmW81sg5m9GrTPBfa3DZxfvhSamQs6rsr7i0hi0/RNIhJzzGw5cLlzbnJgvR3wDXAhMBE4HngF6A4UAuuAw5xzi8ysDd6cdvPM7OLAdY6u5j4ZNZx7Bt4kx8OA7/EmoT/VOTcwcK4DDnLOLanm2i8D3+FNMJ0K5DnnptV0rpm9BCQ55y6o7f4ikrhUcyYi8WAkMME5N8E5V+ac+xCYBZwa2F8GHGxm6c65dc65eXW4dnXnXgU84Jxb4JwrAe4H+tah9mo33iTMbZ1zReWJWXXM7Ba8ZPPSMN1fROKUkjMRiQe5wDmBJs0tgSbKo4E2zrkdwHnA1cA6M3vPzLqHctFazs0Fngi63ybAgHYhxvyHwPFfmNk8M7u0ugPN7BTgd8AZzrmdYbq/iMQpJWciEosqP2+xCnjBOdcsaMlwzj0I4Jyb5JwbCrQBFgJPVXOdn9+o+nNXAVdVume6c256SAVwbr1z7grnXFu8WrD/rap3p5l1A54DznXOrapU5n2+v4jELyVnIhKLfgSCxwF7ERhmZieZWbKZpZnZYDNrb2atzOz0wPNju4ACoDToOu3NLLWqm9Ry7j+B28ysV+DYpmZ2Tg0xVr72OWbWPrC6GS9RLK10TBNgHDCqimbP2u4vIglKyZmIxKIHgFGBJr2bAzVK/wPcDuTj1Sr9Hu9vWBJwE7AWr+lvEPCbwHX+C8wD1pvZhiruU+25zrm3gIeAV8xsG97D/acEnXsX8FwgxnOruPZhwEwzKwDeAX7nnFtW6ZhDgW7AY8G9NkO8v4gkKPXWFBEREYkhqjkTERERiSFKzkRERERiiJIzERERkRii5ExEREQkhig5ExEREYkhSs5EREREYoiSMxEREZEYouRMREREJIYoORMRERGJIf8PAN04+KVVo38AAAAASUVORK5CYII=",
                        "text/plain": [
                            "<Figure size 720x432 with 1 Axes>"
                        ]
                    },
                    "metadata": {
                        "needs_background": "light"
                    },
                    "output_type": "display_data"
                }
            ],
            "source": [
                "plt.scatter(\n",
                "    sim_results_df['test set size'],\n",
                "    sim_results_df['frac oob (observed)'],\n",
                "    5,\n",
                "    label='observed')\n",
                "plt.plot(\n",
                "    sim_results_df['test set size'],\n",
                "    sim_results_df['frac oob (upper bound)'],\n",
                "    label='theoretical bound',\n",
                "    color='orange',\n",
                "    marker='o')\n",
                "plt.ylim(-0.1) # give room for those 0s\n",
                "plt.legend()\n",
                "plt.xticks(test_sizes)\n",
                "plt.title(f'# models = {num_models}')\n",
                "plt.xlabel('test set size')\n",
                "plt.ylabel('fraction out-of-bounds');"
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3.9.7 ('base')",
            "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.9.7"
        },
        "orig_nbformat": 4,
        "vscode": {
            "interpreter": {
                "hash": "18983e9e827802b63835b8dbc18d4f8a4ebd5dd283f63362764c2c714492b89c"
            }
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}