{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MaldiAMRKit - Evaluation & Splitting\n",
    "\n",
    "This notebook covers AMR-specific evaluation metrics, stratified splitting\n",
    "utilities, and integration with scikit-learn's cross-validation tools."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset\n",
    "\n",
    "These notebooks use the **MALDI-Kleb-AI** dataset (Rocchi *et al.*, 2026; [Zenodo DOI 10.5281/zenodo.17405072](https://zenodo.org/records/17405072)), a curated archive of MALDI-TOF spectra of *Klebsiella* clinical isolates from three Italian centres (Rome, Milan, Catania) with Amikacin / Meropenem resistance annotations. For simplicity we restrict the demo to the **Rome sub-cohort** (single site, no batch correction needed). The helper in [`notebooks/_demo.py`](_demo.py) caches the 370 MB tarball under `~/.cache/maldiamrkit/` (or `$MALDIAMRKIT_CACHE_DIR`) on first use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:44.037448Z",
     "iopub.status.busy": "2026-05-13T10:10:44.037349Z",
     "iopub.status.idle": "2026-05-13T10:10:44.686051Z",
     "shell.execute_reply": "2026-05-13T10:10:44.685228Z"
    }
   },
   "outputs": [],
   "source": "import numpy as np\nimport pandas as pd\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.model_selection import StratifiedKFold, cross_val_score\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler\n\nfrom maldiamrkit.evaluation import (\n    CaseGroupedKFold,\n    SpeciesDrugStratifiedKFold,\n    amr_classification_report,\n    case_based_split,\n    categorical_agreement,\n    major_error_rate,\n    me_scorer,\n    sensitivity_score,\n    specificity_score,\n    stratified_species_drug_split,\n    very_major_error_rate,\n    vme_me_curve,\n    vme_scorer,\n)\nfrom maldiamrkit.susceptibility import LabelEncoder"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## AMR Evaluation Metrics\n\nIn clinical antimicrobial resistance testing, standard ML metrics (accuracy,\nF1) do not capture the clinical severity of errors. MaldiAMRKit provides\ndomain-specific metrics following EUCAST conventions:\n\n- **VME (Very Major Error)**: resistant isolate called susceptible - the most\n  dangerous error (could lead to treatment failure)\n- **ME (Major Error)**: susceptible isolate called resistant - wasteful but\n  not immediately dangerous\n- **Sensitivity**: TP / (TP + FN) = 1 - VME\n- **Specificity**: TN / (TN + FP) = 1 - ME\n- **Categorical Agreement**: overall accuracy (TP + TN) / N\n\nWe illustrate the metric definitions with a small hand-crafted vector so\nthe arithmetic is easy to follow, then move on to a full run on the real\ndataset."
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:44.703705Z",
     "iopub.status.busy": "2026-05-13T10:10:44.703528Z",
     "iopub.status.idle": "2026-05-13T10:10:44.708357Z",
     "shell.execute_reply": "2026-05-13T10:10:44.707951Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VME (Very Major Error):   25.00%\n",
      "ME (Major Error):         16.67%\n",
      "Sensitivity:              75.00%\n",
      "Specificity:              83.33%\n",
      "Categorical Agreement:    80.00%\n"
     ]
    }
   ],
   "source": [
    "# Toy 10-sample example for illustrating the metric definitions\n",
    "y_true = np.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0])\n",
    "y_pred = np.array([1, 1, 0, 1, 0, 0, 1, 0, 0, 0])\n",
    "\n",
    "# Very Major Error: 1 out of 4 resistant samples was missed\n",
    "vme = very_major_error_rate(y_true, y_pred)\n",
    "print(f\"VME (Very Major Error):   {vme:.2%}\")\n",
    "\n",
    "# Major Error: 1 out of 6 susceptible samples was falsely called resistant\n",
    "me = major_error_rate(y_true, y_pred)\n",
    "print(f\"ME (Major Error):         {me:.2%}\")\n",
    "\n",
    "# Sensitivity and specificity\n",
    "print(f\"Sensitivity:              {sensitivity_score(y_true, y_pred):.2%}\")\n",
    "print(f\"Specificity:              {specificity_score(y_true, y_pred):.2%}\")\n",
    "\n",
    "# Categorical agreement\n",
    "print(f\"Categorical Agreement:    {categorical_agreement(y_true, y_pred):.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AMR Classification Report\n",
    "\n",
    "Get all metrics in a single dictionary with `amr_classification_report`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:44.709692Z",
     "iopub.status.busy": "2026-05-13T10:10:44.709601Z",
     "iopub.status.idle": "2026-05-13T10:10:44.712698Z",
     "shell.execute_reply": "2026-05-13T10:10:44.712461Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                        vme: 0.2500\n",
      "                         me: 0.1667\n",
      "                sensitivity: 0.7500\n",
      "                specificity: 0.8333\n",
      "      categorical_agreement: 0.8000\n",
      "                n_resistant: 4\n",
      "              n_susceptible: 6\n",
      "                    n_total: 10\n"
     ]
    }
   ],
   "source": [
    "report = amr_classification_report(y_true, y_pred)\n",
    "for key, value in report.items():\n",
    "    if isinstance(value, float):\n",
    "        print(f\"  {key:>25s}: {value:.4f}\")\n",
    "    else:\n",
    "        print(f\"  {key:>25s}: {value}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the Dataset\n",
    "\n",
    "Everything below runs on the real **MALDI-Kleb-AI** Rome cohort. The first\n",
    "call downloads 370 MB; later calls are instantaneous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:44.713880Z",
     "iopub.status.busy": "2026-05-13T10:10:44.713791Z",
     "iopub.status.idle": "2026-05-13T10:10:53.883293Z",
     "shell.execute_reply": "2026-05-13T10:10:53.882955Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "Processing spectra:   0%|          | 0/472 [00:00<?, ?spectrum/s]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "Processing spectra: 100%|██████████| 472/472 [00:00<00:00, 4855.36spectrum/s]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'source': 'Zenodo MALDI-Kleb-AI',\n",
       " 'doi': '10.5281/zenodo.17405072',\n",
       " 'record_id': '17405072',\n",
       " 'md5_tar': 'c14b6c6b4210553962faa7f1dc27d275',\n",
       " 'n_samples': 470,\n",
       " 'n_bins': 6000,\n",
       " 'bin_width': 3,\n",
       " 'antibiotic': 'Amikacin',\n",
       " 'city': 'Rome',\n",
       " 'cache_dir': '/home/ettore/.cache/maldiamrkit/maldi-kleb-ai'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pathlib\n",
    "import sys\n",
    "\n",
    "sys.path.insert(0, str(pathlib.Path.cwd()))  # _demo.py sits next to this notebook\n",
    "from _demo import load_maldi_kleb_ai\n",
    "\n",
    "ds = load_maldi_kleb_ai(antibiotic=\"Amikacin\", verbose=True)\n",
    "ds.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:53.884870Z",
     "iopub.status.busy": "2026-05-13T10:10:53.884759Z",
     "iopub.status.idle": "2026-05-13T10:10:53.890438Z",
     "shell.execute_reply": "2026-05-13T10:10:53.889955Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X shape: (469, 6000)\n",
      "Class counts: {1: 272, 0: 197}\n",
      "Species: {'Klebsiella pneumoniae': 469}\n"
     ]
    }
   ],
   "source": [
    "# Drop species with fewer than 2 samples - stratified splits need >=2\n",
    "species_counts = ds.meta[\"Species\"].value_counts()\n",
    "keep = ds.meta[\"Species\"].isin(species_counts[species_counts >= 2].index)\n",
    "\n",
    "X = ds.X.loc[keep]\n",
    "species = ds.meta.loc[keep, \"Species\"].values\n",
    "\n",
    "enc = LabelEncoder(intermediate=\"susceptible\")  # binary R vs rest\n",
    "y = enc.fit_transform(ds.meta.loc[keep, \"Amikacin\"].values)\n",
    "\n",
    "print(f\"X shape: {X.shape}\")\n",
    "print(f\"Class counts: {pd.Series(y).value_counts().to_dict()}\")\n",
    "print(f\"Species: {pd.Series(species).value_counts().to_dict()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### VME / ME Trade-off Curve on Real Predictions\n",
    "\n",
    "The `vme_me_curve` function computes VME and ME rates at varying decision\n",
    "thresholds. This is analogous to an ROC curve but uses clinically meaningful\n",
    "error types. We fit a quick logistic-regression baseline and read its\n",
    "out-of-fold probability scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:53.891554Z",
     "iopub.status.busy": "2026-05-13T10:10:53.891452Z",
     "iopub.status.idle": "2026-05-13T10:10:54.545855Z",
     "shell.execute_reply": "2026-05-13T10:10:54.545134Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAGGCAYAAADmRxfNAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjcJJREFUeJzs3XdYU2cbBvA7YW9QtqK490RFHHXhrIMOwVF3HVXr+tyLYl3Vaq11tVpnVRyt2ta9sK6qteJEXDhaleFgiMyc74+3CUaGhHUC3L/rOleSN+ecPAkvkOe8SyFJkgQiIiIiIqJcUModABERERERFX5MLIiIiIiIKNeYWBARERERUa4xsSAiIiIiolxjYkFERERERLnGxIKIiIiIiHKNiQUREREREeUaEwsiIiIiIso1JhZERERERJRrTCyIiGTwxRdfQKFQyPb6Fy5cQJMmTWBhYQGFQoHg4GAAwIEDB1C3bl2YmppCoVDg5cuXssVI2tavXw+FQoH79+9rytzd3dG5c+cCjUOhUOCLL77I8/N26tQJgwcPzvPzyqVx48aYOHGi3GEQFSgmFkSFTNeuXWFubo7Y2NhM9+nduzeMjY3x7NkzAOKLgEKhwKeffprh/tOmTdPsExUVpSnv37+/pvztzdTUNNsx/+9//0P16tUzff7+/fua886ePTvT96RQKGBpaalV3rJly0xjrFq1aqavmdVxb2758QVKbsnJyejevTueP3+Ob775Bps2bULZsmXx7Nkz+Pr6wszMDMuXL8emTZtgYWGR5bnu3r2LoUOHonz58jA1NYW1tTWaNm2Kb7/9Fq9fvy6gdyS/ly9fapKxkJAQucMpdE6fPo1Dhw5h0qRJmrKgoCCt30UDAwM4Ojri448/LhSf8aRJk7B8+XI8ffpU7lCICoyh3AEQkW569+6N3377Dbt27ULfvn3TPR8fH489e/agQ4cOKFmypKbc1NQUP//8M1asWAFjY2OtY7Zu3QpTU1MkJCSkO5+JiQnWrFmTrtzAwCDbMe/duxddunR5536mpqbYunUrpk+frlX+6tUr7NmzJ9NkpnTp0pg3b166chsbm0xfa9q0aVqJ1oULF7B06VJMnToV1apV05TXrl37nXEXNnfv3sWDBw+wevVqrc/gwIEDiI2NxZdffglvb+93nmfv3r3o3r07TExM0LdvX9SsWRNJSUk4deoUJkyYgOvXr+OHH37Iz7eiN3bs2AGFQgFnZ2ds3rw50wQ5N/r06YMePXrAxMQkz8+ti9evX8PQMG+/PixcuBBt2rRBxYoV0z03atQoNGzYEMnJybhy5QpWrVqFoKAgXLt2Dc7OznkaR17q1q0brK2tsWLFCsyaNUvucIgKhkREhUp8fLxkZWUltW/fPsPnt2zZIgGQAgMDNWUAJB8fH0mpVEq7d+/W2v/06dMSAOmjjz6SAEiRkZGa5/r16ydZWFjkKt67d+9KAKTjx49nuk9YWJgEQPrwww8lAFJwcLDW85s3b5aMjIykLl26pIunRYsWUo0aNXIVoyRJ0o4dO94ZpyRJUlxcXK5fS5Ikyd/fX5LrT/CJEyckANKOHTu0yjds2CABkC5cuPDOc9y7d0+ytLSUqlatKj1+/Djd87dv35aWLFmSJ/Hm1Ween9577z3pww8/lMaOHSuVK1euwF63bNmy0vvvv19gr5cfwsPDJUNDQ2nNmjVa5cePH8+wnq5cuVICIH311VcFGWaOjBw5UipbtqykUqnkDoWoQLArFFEhY2Zmhg8//BBHjx5FREREuue3bNkCKysrdO3aVau8VKlSeO+997Blyxat8s2bN6NWrVqoWbNmvsS7d+9e2NjYoFmzZu/c18vLC+XKlcswxg4dOqBEiRL5EmNm1OMgbty4gV69esHOzk7zPq5cuYL+/ftrugA5Oztj4MCBmu5nbzp16hQaNmwIU1NTVKhQAd9//32mr/nTTz/Bw8MDZmZmKFGiBHr06IFHjx5lO+Zjx46hefPmsLCwgK2tLbp166bVbaR///5o0aIFAKB79+5QKBRo2bIlWrZsiX79+gEAGjZsCIVCgf79+2f6OgsWLEBcXBx+/PFHuLi4pHu+YsWKGD16NIC0rm7r169Pt9/b3c0y+8y//vprKBQKPHjwIN05pkyZAmNjY7x48UJTdu7cOXTo0AE2NjYwNzdHixYtcPr06Sw/u5x6+PAhTp48iR49eqBHjx4ICwvDmTNn0u3XsmVL1KxZE1euXEGLFi1gbm6OihUrYufOnQCAEydOwNPTE2ZmZqhSpQqOHDmidXxGYywysmHDBhgaGmLChAkAgOfPn2P8+PGoVasWLC0tYW1tjY4dO+Ly5cvpjk1ISMAXX3yBypUrw9TUFC4uLvjwww9x9+5dzT6Z/czu3LmD/v37w9bWFjY2NhgwYADi4+Pf+fnt3bsXKSkp2WolA4DmzZsDgFZMAHDp0iV07NgR1tbWsLS0RJs2bfDnn39q7aP+DE+dOoVRo0bBwcEBtra2GDp0KJKSkvDy5Uv07dsXdnZ2sLOzw8SJEyFJktY5VCoVlixZgho1asDU1BROTk4YOnSoVv1Ta9u2LR48eKAZw0RU1DGxICqEevfujZSUFGzfvl2r/Pnz5zh48CA++OADmJmZpTuuV69e+O233xAXFwcASElJwY4dO9CrV68sXy8qKirdFhMTk61Y9+3bh7Zt22a760TPnj0RGBio+WceFRWFQ4cOZRljampqhjG+evUqW6/5Lt27d0d8fDzmzp2rGVx6+PBh3Lt3DwMGDMB3332HHj16IDAwEJ06ddL6InL16lW0a9cOERER+OKLLzBgwAD4+/tj165d6V5nzpw56Nu3LypVqoTFixdjzJgxOHr0KN57771sDaI+cuQI2rdvr3mtcePG4cyZM2jatKnmy+jQoUMxdepUAKKLyaZNmzBt2jRMmzYNQ4YMAQDMmjULmzZtwtChQzN9rd9++w3ly5dHkyZNsvsx6uTtz9zX1xcKhSJdnQeA7du3o127drCzswMgkqv33nsPMTEx8Pf3x9y5c/Hy5Uu0bt0a58+fz/NYt27dCgsLC3Tu3BmNGjVChQoVsHnz5gz3ffHiBTp37gxPT08sWLAAJiYm6NGjB7Zt24YePXqgU6dOmD9/Pl69eoWPP/44y7FUGfnhhx8wYMAATJ48GQsXLgQA3Lt3D7t370bnzp2xePFiTJgwAVevXkWLFi3w+PFjzbGpqano3LkzAgIC4OHhgUWLFmH06NGIjo7GtWvX3vnavr6+iI2Nxbx58+Dr64v169cjICDgncedOXMGJUuWRNmyZbP1HtV1Wf3zBoDr16+jefPmuHz5MiZOnIgZM2YgLCwMLVu2xLlz59Kd4/PPP8ft27cREBCArl274ocffsCMGTPQpUsXpKamYu7cuWjWrBkWLlyITZs2aR07dOhQTJgwQTOWaMCAAdi8eTPat2+P5ORkrX09PDwAIN+SWiK9I3OLCRHlQEpKiuTi4iJ5eXlpla9atUoCIB08eFCrHIA0YsQI6fnz55KxsbG0adMmSZIkae/evZJCoZDu37+v6ZrzdlcoABlumXXFetOrV68kU1NTad26dVnup+4KtXDhQunatWsSAOnkyZOSJEnS8uXLJUtLS+nVq1cZds1q0aJFpjEOHTr0nTGqZdQVSv2Z9OzZM93+8fHx6cq2bt0qAZD++OMPTZmPj49kamoqPXjwQFN248YNycDAQKsr1P379yUDAwNpzpw5Wue8evWqZGhomK48I3Xr1pUcHR2lZ8+eacouX74sKZVKqW/fvpqyzLqYrFu3LltdoaKjoyUAUrdu3d4ZkySl/XwzqgcAJH9/f83jrD5zLy8vycPDQ6vs/PnzEgBp48aNkiRJkkqlkipVqiS1b99eq/tJfHy8VK5cOalt27bZilkXtWrVknr37q15PHXqVMne3l5KTk7W2k9dV7ds2aIpu3nzpgRAUiqV0p9//qkpP3jwYLrPTP3zCQsL05S92RXq22+/lRQKhfTll19qvW5CQoKUmpqqVRYWFiaZmJhIs2bN0pStXbtWAiAtXrw43Xt887PM7Gc2cOBArWM++OADqWTJkunO9bZmzZql+7lKUlo9Xbt2rRQZGSk9fvxYOnDggFSxYkVJoVBI58+f1+zr4+MjGRsbS3fv3tWUPX78WLKyspLee+89TZn6M3y7fnh5eUkKhUIaNmyYpiwlJUUqXbq01KJFC03ZyZMnJQDS5s2btWI9cOBAhuWSJEnGxsbSZ5999s7PgagoYIsFUSFkYGCAHj164OzZs1rdIrZs2QInJye0adMmw+Ps7OzQoUMHbN26VbN/kyZNsrxSaGpqisOHD6fb5s+f/844jx07hsTERHTs2DHb761GjRqoXbu2VozdunWDubl5pse4u7tnGOOYMWOy/bpZGTZsWLqyN1uEEhISEBUVhcaNGwMA/v77bwDiCvDBgwfh4+ODMmXKaPavVq0a2rdvr3W+X375BSqVCr6+vlqtLs7OzqhUqRKOHz+eZYxPnjxBcHAw+vfvr9VlrHbt2mjbti327dun+xvPhLq1ysrKKs/O+baMPnM/Pz9cvHhRqwvMtm3bYGJigm7dugEAgoODcfv2bfTq1QvPnj3Tar1q06YN/vjjD6hUqjyL88qVK7h69Sp69uypKevZsyeioqJw8ODBdPtbWlqiR48emsdVqlSBra0tqlWrBk9PT025+v69e/eyFceCBQswevRofPXVV+kmPzAxMYFSKf7dp6am4tmzZ7C0tESVKlU0dRUAfv75Z9jb2+Pzzz9Pd/7sTI389s+sefPmePbs2TtbN589e6bV+vC2gQMHwsHBAa6urujQoQOio6OxadMmNGzYUPOeDh06BB8fH5QvX15znIuLC3r16oVTp06li2HQoEFa78nT0xOSJGHQoEGaMgMDAzRo0EDrZ7Bjxw7Y2Nigbdu2Wr+nHh4esLS0zPD31M7OTmu2PaKijLNCERVSvXv3xjfffIMtW7Zg6tSp+Oeff3Dy5EmMGjUqyxmbevXqhT59+uDhw4fYvXs3FixYkOXrGBgYZLvv89v27t2LBg0awMnJSafjevXqhUWLFmHs2LE4c+aMputOZiwsLHIcY3aUK1cuXdnz588REBCAwMDAdGNdoqOjAQCRkZF4/fo1KlWqlO74KlWqaH3Zv337NiRJynBfADAyMgIAxMXFabqyAeLn4+DgoBl7UKVKlXTHVqtWDQcPHsSrV6/eOX1sdlhbWwOAzt10dJHRZ969e3eMGzcO27Ztw9SpUyFJEnbs2KHpVw+IzxGAZrxIRqKjozP9Ivv21KA2NjYZditU++mnn2BhYYHy5cvjzp07AEQy7u7ujs2bN+P999/X2r906dLpvqTb2NjAzc0tXRmADPvtv+3EiRPYu3cvJk2apBlX8SaVSoVvv/0WK1asQFhYGFJTUzXPvTlz3N27d1GlSpUcz/j0ZvIMpHVVevHihebnkxnprXEMb5o5cyaaN2+OuLg47Nq1C4GBgZpECRC/Z/Hx8ZnWfZVKhUePHqFGjRqZxqr+vDP6Obz5M7h9+zaio6Ph6OiYYawZjXuTJEnWNWuIChITC6JCysPDA1WrVsXWrVsxdepUbN26FZIkoXfv3lke17VrV5iYmKBfv35ITEyEr69vvsW4b98+DBgwQOfjevbsiSlTpmDw4MEoWbIk2rVrlw/RZV9GXyx9fX1x5swZTJgwAXXr1oWlpSVUKhU6dOiQoyviKpUKCoUC+/fvzzAxVK/f8fXXX2v1Wy9btuw7B/PmNWtra7i6umar3z2Q+dXuN7/gvi2jz9zV1RXNmzfH9u3bMXXqVPz55594+PAhvvrqK80+6s9+4cKFqFu3bobnfnstlDe9PRB93bp1mQ5ilyQJW7duxatXrzJcpyUiIgJxcXFar5dZ0p9ZeVZfuNVq1KiBly9fasbFvJ2UzZ07FzNmzMDAgQPx5ZdfokSJElAqlRgzZkyett7k9D2ULFkyywSqVq1amgsHPj4+iI+Px+DBg9GsWbN0iUBuY82o/M34VSoVHB0dMx1D4+DgkK7s5cuXsLe3z1GcRIUNEwuiQqx3796YMWMGrly5gi1btqBSpUqa7gGZMTMzg4+PD3766Sd07Ngx3/7hXbt2DQ8fPkx3xTY7ypQpg6ZNmyIoKAifffZZns+Zn1svXrzA0aNHERAQgJkzZ2rK1VfL1RwcHGBmZpauHABCQ0O1HleoUAGSJKFcuXKoXLlypq/dt29frRm21F/A1d3Z3j4vANy8eRP29vZ50lqh1rlzZ/zwww84e/YsvLy8stxXfeX67QHoGc3w9C5+fn4YPnw4QkNDsW3bNpibm2utkVKhQgUAIvnJSSvW4cOHtR6/eZX7bSdOnMA///yDWbNmaa19Aog6MmTIEOzevRuffPKJznHowt7eHjt37kSzZs3Qpk0bnDp1Cq6urprnd+7ciVatWuHHH3/UOu7tL7wVKlTAuXPnkJycrGkhKwhVq1bFzz//nO3958+fj127dmHOnDlYtWoVHBwcYG5unmndVyqVOU5A3lahQgUcOXIETZs2zbIlS+3ff/9FUlJSuvpBVFRxjAVRIaZunZg5cyaCg4Pf2VqhNn78ePj7+2PGjBn5Ftu+ffvg5OSEBg0a5Oj42bNnw9/fP8P+3nJTX9V8+0rskiVL0u3Xvn177N69Gw8fPtSUh4SEpOt//+GHH8LAwAABAQHpzitJkmYa2/Lly8Pb21uzNW3aFIC40l63bl1s2LBB6wv8tWvXcOjQIXTq1ClX7/ltEydOhIWFBT799FOEh4ene/7u3bv49ttvAYgv+fb29vjjjz+09lmxYoXOr/vRRx/BwMAAW7duxY4dO9C5c2ethMnDwwMVKlTA119/rdVlTC0yMjLL87/52Xp7e2c4la6auhvUhAkT8PHHH2ttgwcPRqVKlTK9sp3XSpcujSNHjuD169do27at1rTHBgYG6erUjh078O+//2qVffTRR4iKisKyZcvSnT87LSc55eXlhRcvXmR7PEmFChXw0UcfYf369Xj69CkMDAzQrl077NmzR6v1Ljw8HFu2bEGzZs3e2RUru3x9fZGamoovv/wy3XMpKSnpkueLFy8CQL7Nnkakb/TrMiAR6aRcuXJo0qQJ9uzZAwDZTizq1KmDOnXqZGvflJQU/PTTTxk+98EHH2R6FXzv3r3o2LFjjvsWt2jRQrPewrtER0dnGmN+XC22trbGe++9hwULFiA5ORmlSpXCoUOHEBYWlm7fgIAAHDhwAM2bN8fw4cORkpKC7777DjVq1MCVK1c0+1WoUAGzZ8/GlClTcP/+ffj4+MDKygphYWHYtWsXhgwZgvHjx2cZ18KFC9GxY0d4eXlh0KBBeP36Nb777jvY2NhorTuQFypUqIAtW7bAz88P1apV01p5+8yZM9ixY4dWF6JPP/0U8+fPx6effooGDRrgjz/+wK1bt3R+XUdHR7Rq1QqLFy9GbGws/Pz8tJ5XKpVYs2YNOnbsiBo1amDAgAEoVaoU/v33Xxw/fhzW1tb47bffcvv2kZiYiJ9//hlt27bNdEX4rl274ttvv0VERESmffLzUsWKFXHo0CG0bNkS7du3x7Fjx2BtbY3OnTtj1qxZGDBgAJo0aYKrV69i8+bNWgOdAdEatnHjRowbNw7nz59H8+bN8erVKxw5cgTDhw/XDJDPa++//z4MDQ1x5MgRzZTH7zJhwgRs374dS5Yswfz58zF79mwcPnwYzZo1w/Dhw2FoaIjvv/8eiYmJ7xxHposWLVpg6NChmDdvHoKDg9GuXTsYGRnh9u3b2LFjB7799lt8/PHHmv0PHz6MMmXKoF69enkWA5E+Y2JBVMj17t0bZ86cQaNGjVCxYsU8P39iYiL69OmT4XNhYWEZJhbR0dE4c+YMRo4cmefxZOSff/7JNMb86oayZcsWfP7551i+fDkkSUK7du2wf/9+rS4ogJiV6eDBgxg3bhxmzpyJ0qVLIyAgAE+ePNFKLABg8uTJqFy5Mr755hvNOAo3Nze0a9cu3YKHGfH29saBAwfg7++PmTNnwsjICC1atMBXX32V4WDo3OratSuuXLmChQsXYs+ePVi5ciVMTExQu3ZtLFq0SLPmByBa1SIjI7Fz505s374dHTt2xP79+3P0hdvPzw9HjhyBlZVVhi0xLVu2xNmzZ/Hll19i2bJliIuLg7OzMzw9PbNcm0MXe/fuxcuXL7W6Yb2tS5cuWLRoEQIDAzFq1Kg8ed13qVWrFvbv3w9vb2906dIFBw4cwNSpU/Hq1Sts2bIF27ZtQ/369bF3715MnjxZ61gDAwPs27cPc+bMwZYtW/Dzzz+jZMmSaNasGWrVqpVvMTs5OaFTp07Yvn17thOLBg0aoGXLlli5ciWmTJmCGjVq4OTJk5gyZQrmzZsHlUoFT09P/PTTT1qzbeWFVatWwcPDA99//z2mTp0KQ0NDuLu745NPPtG0IAJiPMbPP/+cbgYqoqJMIeVn+yYRFUvbt29H7969ERUVpZlthYgoMydPnkTLli1x8+bNTGdGK2x2796NXr164e7du1l2qSMqSjjGgojynK2tLZYuXcqkgoiypXnz5mjXrl2edluS21dffYWRI0cyqaBihS0WRERERESUa2yxICIiIiKiXGNiQUREREREucbEgoiIiIiIco2JBRERERER5VqxW8dCpVLh8ePHsLKy4rzSRERERERZkCQJsbGxcHV1hVKZdZtEsUssHj9+DDc3N7nDICIiIiIqNB49eoTSpUtnuU+xSyysrKwAiA/H2tpalhhUKhUiIyPh4ODwzsyPijbWBQJYDygN6wKpsS4QoB/1ICYmBm5ubprv0FkpdomFuvuTtbW1rIlFQkICrK2t+ceimGNdIID1gNKwLpAa6wIB+lUPsjOEgDWViIiIiIhyjYkFERERERHlGhMLIiIiIiLKNSYWRERERESUa0wsiIiIiIgo15hYEBERERFRrjGxICIiIiKiXJM1sfjjjz/QpUsXuLq6QqFQYPfu3e88JigoCPXr14eJiQkqVqyI9evX53ucRERERESUNVkTi1evXqFOnTpYvnx5tvYPCwvD+++/j1atWiE4OBhjxozBp59+ioMHD+ZzpERERERElBVZV97u2LEjOnbsmO39V61ahXLlymHRokUAgGrVquHUqVP45ptv0L59+/wKk4iIiIiI3kHWxEJXZ8+ehbe3t1ZZ+/btMWbMmEyPSUxMRGJiouZxTEwMALFEukqlypc4s/L0KdCsmQKpqfa4e7fgX5/0i0qlgiRJstRF0h+sB6TGukBqrAsE6Ec90OW1C1Vi8fTpUzg5OWmVOTk5ISYmBq9fv4aZmVm6Y+bNm4eAgIB05ZGRkUhISMi3WDPz7JkSYWGOAAwRHv4YBgYcP1+cqVQqREdHQ5IkKJWsC8UV6wGpsS6QGusCAYDV2LEoeeYMYmbORNL778sSQ2xsbLb3LVSJRU5MmTIF48aN0zyOiYmBm5sbHBwcYG1tXeDxGL7xiZcs6QhjY/6xKM5UKhUUCgUcHBz4j6MYYz0gNdYFUmNdIADAixdQPnwIG6USCkdHWUIwNTXN9r6FKrFwdnZGeHi4Vll4eDisra0zbK0AABMTE5iYmKQrVyqVsvyivhlKaqo8MZB+USgUstVH0h+sB6TGukBqrAsk/Xerrgty0OV1C1VN9fLywtGjR7XKDh8+DC8vL5ki0p2RUdr95GT54iAiIiIiykuyJhZxcXEIDg5GcHAwADGdbHBwMB4+fAhAdGPq27evZv9hw4bh3r17mDhxIm7evIkVK1Zg+/btGDt2rBzh5wgTCyIiIiIqimRNLP766y/Uq1cP9erVAwCMGzcO9erVw8yZMwEAT5480SQZAFCuXDns3bsXhw8fRp06dbBo0SKsWbOmUE01q1QCBgaiYSssTOZgiIiIiIjyiKxjLFq2bAlJkjJ9PqNVtVu2bIlLly7lY1T5S6EASpUCHj4E1qxRoFEjuSMiIiIiIsq9QjXGoqjw8RG3bLEgIiIioqKCiYUMmjQRrTQyLKNBRERERJQvmFjIQL3G34MHABfUJCIiIqIMZTFkQB8xsZBB7dqAhYUKjx4p8OOPckdDRERERHpJnVgUkrVMCkeURYytLTB27CsAwKZN8sZCRERERHpK3bVFoZA3jmxiYiGTatXEIhZxcTIHQkRERET6jYkFZcXERNwmJsobBxERERHpKXVXKCYWlBUrK9G0FRUlcyBEREREpJ+YWFB2lC+fCgCIiABevpQ3FiIiIiLSQxy8TdlhaSnB1VVUlps3ZQ6GiIiIiPQPB29TdpUpI243bpQ3DiIiIiLSQ+wKRdllbi5uV66UNw4iIiIi0kNMLCi7ZsxIW00xPFzGQIiIiIhI/3CMBWVX8+Zp93ftki8OIiIiItJDHGNB2aVQAHPnivtHj8obCxERERHpGXaFIl14eYnbixfljYOIiIiI9Exysrg1MpI3jmxiYiGz8uXFbVgYsGKFvLEQERERkR5JShK3xsbyxpFNTCxkVqpU2n2OsyAiIiIiDbZYkC4MDIDdu8X9Fy9kDYWIiIiI9AkTC9KVvb24ff5c3jiIiIiISI8wsSBdlSsnbh88ABIS5I2FiIiIiPRESoq4NTSUN45sYmKhB1xcACsrMVXxrVtyR0NEREREekGdWLDFgrJLoUgb7L9+vayhEBEREZG+YIsF5YR6nAW7QhERERERgLQxFkwsSBfTp4vbH38ELlyQNxYiIiIi0gNssaCc6N0b+OADsQ7K//4ndzREREREJCuVCorXr8V9c3N5Y8kmJhZ6QqEAli0T61qcPAncuSN3REREREQkG3VSAQAWFvLFoQMmFnrE1RWoU0fcDw2VNxYiIiIiklFcXNp9MzP54tABEws94+wsbp8+lTcOIiIiIpLRq1cAAJWZGaAsHF/ZC0eUxYiTk7h98kTeOIiIiIhIRv8lFlIh6QYFMLHQO5UqidubN+WNg4iIiIhk9F9XKKmQDNwGmFjonRo1xO3mzcCzZ/LGQkREREQyUbdYMLGgnKpbN+3+11/LFgYRERERyYmJBeVWmTJAx47i/tmz8sZCRERERDLhGAvKC198IW5PnAAiImQNhYiIiIjkwDEWlBcqVky7v3OnfHEQERERkUzYFYryQokSQNOm4v6DB/LGQkREREQyYIsF5ZXBg8Xt3r2AJMkbCxEREREVsNhYAIBkaSlzINnHxEJPdekCmJsD168DGzfKHQ0RERERFSh1YsHB25RbJUoAM2aI+xMmACqVvPEQERERUQF68gQAoLK3lzmQ7GNiocd69xa3kZFAhw7yxkJEREREBej+fQBASpky8sahAyYWeqx06bQF806dkjUUIiIiIipIkZEAAJWDg8yBZB8TCz2mUAB//CHuv37NNS2IiIiIio2YGACAZG0tcyDZx8RCz1lZAbVri/vqJIOIiIiIirDUVM10syrOCpV9y5cvh7u7O0xNTeHp6Ynz589nuf+SJUtQpUoVmJmZwc3NDWPHjkVCQkIBRSuP+vXFbUiIvHEQERERUQH4L6kAAMnKSsZAdCNrYrFt2zaMGzcO/v7++Pvvv1GnTh20b98eEZn0+dmyZQsmT54Mf39/hISE4Mcff8S2bdswderUAo68YFWuLG5v3ZI3DiIiIiIqADdvAgCkkiUBExOZg8k+WROLxYsXY/DgwRgwYACqV6+OVatWwdzcHGvXrs1w/zNnzqBp06bo1asX3N3d0a5dO/Ts2fOdrRyFHRMLIiIiomIkLEzc1qwpbxw6ki2xSEpKwsWLF+Ht7Z0WjFIJb29vnD17NsNjmjRpgosXL2oSiXv37mHfvn3o1KlTgcQsF3ViERrK9SyIiIiIiryXL8WtnZ2sYejKUK4XjoqKQmpqKpycnLTKnZyccPO/5p+39erVC1FRUWjWrBkkSUJKSgqGDRuWZVeoxMREJCYmah7H/DfCXqVSQSXTt3SVSgVJkrL9+hUrApaWCkRHK3DunAqenvkcIBUYXesCFU2sB6TGukBqrAvF3IsXUELMCCV3PdDltWVLLHIiKCgIc+fOxYoVK+Dp6Yk7d+5g9OjR+PLLLzFDvUz1W+bNm4eAgIB05ZGRkbIN+lapVIiOjoYkSVAqs9do1KaNDfbsMcOkSckIDHyBbB5Gei4ndYGKHtYDUmNdIDXWheLN+b+L5vHGxnj58qWs9SA2Njbb+8qWWNjb28PAwADh4eFa5eHh4XB2ds7wmBkzZqBPnz749NNPAQC1atXCq1evMGTIEEybNi3DD3zKlCkYN26c5nFMTAzc3Nzg4OAAa5nmBVapVFAoFHBwcMh2JfHzA/bsAU6eNMHBg47o1y+fg6QCkZO6QEUP6wGpsS6QGutCMfZf7xoAMGvUCLa2trLWA1NT02zvK1tiYWxsDA8PDxw9ehQ+Pj4AxC/R0aNHMXLkyAyPiY+PT/ehGhgYAAAkScrwGBMTE5hkMJpeqVTK+ouqUCh0iuGDD9LuX7+uZItFEaJrXaCiifWA1FgXSI11oZh6+FBzVzFoEBQREbLWA11eV9aaOm7cOKxevRobNmxASEgIPvvsM7x69QoDBgwAAPTt2xdTpkzR7N+lSxesXLkSgYGBCAsLw+HDhzFjxgx06dJFk2AUVaamwLffivuhofLGQkRERET55PlzcVu1qrxx5ICsYyz8/PwQGRmJmTNn4unTp6hbty4OHDigGdD98OFDrSxp+vTpUCgUmD59Ov799184ODigS5cumDNnjlxvoUB5eIjb4GBZwyAiIiKi/LJzp7jNZGiAPpN98PbIkSMz7foUFBSk9djQ0BD+/v7w9/cvgMj0j6OjuH2j6x0RERERFRWbNwPLl4v75cvLG0sOsNNeIWJuLm7j4+WNg4iIiIjywZIlafenT5ctjJxiYlGIqBOLlBTg8GEgKkreeIiIiIgoj8THA3/9Je5v3AiUKydvPDnAxKIQUScWANCuHdCoEZDJZFhEREREVJgcPZp2389PvjhygYlFIWJiAkyaBNSsCSiVQFgY8OiR3FERERERUa6p13arUwcwNpY3lhxiYlHIzJ8PXL0q6hwA7NrFVgsiIiKiQi8lRdwWwi5QakwsCqnu3cXtmDFA/frAihVAdLSsIRERERFRTr16JW7NzOSNIxeYWBRSY8YA/fqJlrLgYGDECKBsWeDcObkjIyIiIiKdbd0qbsuWlTeOXGBiUUiZmQHr1wOPH4uZyapUES0WXbsC9+/LHBwRERERZZ8kAQ8fivvsCkVyKVkSGD1azE5Wty4QEQF4ewNPnsgdGRERERFly8uXQGSkuN+zp6yh5AYTiyLC0hL47TfAxQW4e1e0XBARERFRIZCQIG6VSsDKSt5YcoGJRRFSujTwww/i/l9/ATt3csYoIiIiIr0XGytuLS3ljSOXmFgUMR06AM7O4n737kD79sDt2/LGRERERERZ2L5d3JYsKW8cucTEoogxNBSzRE2eLBbUO3xYjL1Yt46tF0RERER6ads2cVuzprxx5BITiyLIyQmYNw+4fh1o1QqIjwcGDgQGDQKSkuSOjoiIiIi0PHsmbj//XN44comJRRFWoYJosZg7V4wFWrcOaN5cJBnjx3PmKCIiIiK9oB5jUYinmgUAQ7kDoPxlYABMmSK6Q/n6AufPiw0ADh4ETp0CbGxkDZGIiIio+Dp3DoiLE/etreWNJZfYYlFMdOwI/P03sGCBaMFwcQGuXQNmzJA7MiIiIqJirF+/tPuF/GovWyyKkUqVgAkTxP1SpUQ9vnxZ3piIiIiIiq2EBOD+fXF/yRIx804hxsSimFJPOnD6NLB4MWBhAdjaAj4+hb5OExERERUOhw8DiYmAqyswapTc0eQau0IVU/XrAx9/DKSmAv/7HzBsGNCjB7B+vdyRERERERUTJ0+K2zJlAIVC3ljyAFssirEVK8Sq8S9fAiEhwM2bYopaIiIiIioAKSni1sFB3jjyCBOLYszBAVi7VtzfvBn45BPgwAHRIsfuUERERET5LCxM3LZtK28ceYRdoQgA0LWrGGdx+zbQoIGYQYqIiIiI8okkpc2iU7GivLHkESYWBEB0idq5U7RiXLsGNGoE+PtzpW4iIiKifHHlimixMDcXKxgXAUwsSKNDBzHGQj2oe9YskWAcOiSSaiIiIiLKIxcuiNsmTQBLS3ljySNMLEiLgwOwYwcQGAiUKCFa6Nq3B957Dzh7Vu7oiIiIiIqIhw/Fbfny8saRh5hYUIb8/MRMUWPGiIHcp04B3t5AcrLckREREREVAfv2idvq1eWNIw8xsaBMOToC33wjBnQDQHw8YGwsEg0HB655QURERJQjMTHAxYvivouLvLHkISYW9E5ubmKshVpSEhAVBQwYAHz/vXxxERERERVK6quzxsaiz3kRwXUsKFvOnAGePEl7vHixaM0YPRpo2BAwNc34OFtbsUo9EREREf1H/aXq008BGxt5Y8lDTCwoWwwMgNKl0x4vWiS6BoaGAh4eWR87YQIwb544BxEREVGxd/WquK1SRd448hi7QlGOKBTAjBmiNcLePuOtZEmx78KFQK9enLKWiIiICJIEnDsn7nt6yhtLHmNiQTnWuzfw779AZGTGW1QUsGWLaKnYvh04f17uiImIiIhkFhIiviQZGQF168odTZ5iYkH5qmdPoFIlcb9xY+DgQXnjISIiIpLV8uXitkMHMdVmEcLEgvLd3LliZilA/A5ZWoqpbH//Xd64iIiIiArU8ePAmjXi/tix8saSD5hYUL774APg1i1g6FDx+NUr0VXql1/kjYuIiIiowJw+DbRuLebtf+89oGVLuSPKc0wsqECYmgKrVgHh4cDUqaLsr78AlUreuIiIiIgKxOHDafd/+knMhFPEMLGgAuXoCIwbJ7pDXb0KHDokd0RERERE+SwlJW1sxerVaX3EixgmFlTgSpYE2rQR9+/flzUUIiIiovwXESFmgjIwAPr3lzuafMPEgmTh4iJuv/8eeP1a3liIiIiI8pV6xhoHB8Cw6K5PzcSCZDFtmvjdCg4GunYFYmPljoiIiIgon1y+LG6L+OBSJhYki9KlgZ07xViLI0eAJk2AkyfljoqIiIgoj4WEpE0xu2KFvLHkMyYWJJv33gOCgsSA7mvXxOMPPxStGXv3ihXviYiIiAq1ESPEFLOdOokvOkVYjhKLu3fvYvr06ejZsyciIiIAAPv378f169fzNDgq+jw8xOxQQ4aIWdd27RIL6nXuDNSrJ1o1inirIRERERVV166JRfEMDUVrRRGcYvZNOicWJ06cQK1atXDu3Dn88ssviIuLAwBcvnwZ/v7+OgewfPlyuLu7w9TUFJ6enjh//nyW+798+RIjRoyAi4sLTExMULlyZezbt0/n1yX94egoBnH//TcwcSIwaJDoInX5MtC9e5FcmJKIiIiKuitXgFq1xH0vL6BsWXnjKQA6JxaTJ0/G7NmzcfjwYRgbG2vKW7dujT///FOnc23btg3jxo2Dv78//v77b9SpUwft27fXtIK8LSkpCW3btsX9+/exc+dOhIaGYvXq1ShVqpSub4P0UN26wFdfiW6I9+8DU6aI8pUrgY4d07avvxbTQRMRERHpjZSUtO3HH4H69dOeK+JdoNR0nu/q6tWr2LJlS7pyR0dHREVF6XSuxYsXY/DgwRgwYAAAYNWqVdi7dy/Wrl2LyZMnp9t/7dq1eP78Oc6cOQMjIyMAgLu7u65vgQqBkiVFl6gTJ4AzZ4ADB9KeO3BAdJFavx6oWlW2EImIiIiE4cPFldC3ffABsGQJUKZMgYckB51bLGxtbfHkyZN05ZcuXdKp5SApKQkXL16Et7d3WjBKJby9vXH27NkMj/n111/h5eWFESNGwMnJCTVr1sTcuXORmpqq69ugQuKXX4BNm0QSsX49sHAhYGMDnDsHVK8OfPyxuE9EREQki4sXgVWrtMsMDIDJk4Gffy42SQWQgxaLHj16YNKkSdixYwcUCgVUKhVOnz6N8ePHo2/fvtk+T1RUFFJTU+Hk5KRV7uTkhJs3b2Z4zL1793Ds2DH07t0b+/btw507dzB8+HAkJydnOr4jMTERiYmJmscxMTEAAJVKBZVMo4JVKhUkSZLt9QsTBwegVy/tsu7dgZEjFfj9dwV+/ln8zjo4SFAqRdIxd66EDz6QJ15dsS4QwHpAaVgXSI11oZDYswfK/7o5SR06QPrpJ1FubAxYWIgpLnMxzaU+1ANdXlvnxGLu3LkYMWIE3NzckJqaiurVqyM1NRW9evXC9OnTdT2dTlQqFRwdHfHDDz/AwMAAHh4e+Pfff7Fw4cJME4t58+YhICAgXXlkZCQSEhLyNd7MqFQqREdHQ5IkKJWc8VdXJibA6tXAzZuGWLnSArt2mSIyUsyyEB4O+PoCs2fHol+/eOj7x8u6QADrAaVhXSA11gX9Z3jjBkr6+QEAkitXxvNvvoGUnCyeTE4GXr3K9WvoQz2I1WEVY50TC2NjY6xevRozZ87E1atXERcXh3r16qFSpUo6ncfe3h4GBgYIDw/XKg8PD4ezs3OGx7i4uMDIyAgGBgaasmrVquHp06dISkrSGkyuNmXKFIwbN07zOCYmBm5ubnBwcIC1tbVOMecVlUoFhUIBBwcH/rHIBUdHsfZFVJSEx4/F1YDvvlNg7VoFpk61RmCgFWbNktCwoe7ntrcXM8PlN9YFAlgPKA3rAqmxLhQC69dDkZwMqU4dGJw7B4f/xv/mJX2oB6amptneV+evTrNmzcL48ePh5uYGNzc3Tfnr16+xcOFCzJw5M1vnMTY2hoeHB44ePQofHx8A4sM7evQoRo4cmeExTZs2xZYtW6BSqTQf7q1bt+Di4pJhUgEAJiYmMDExSVeuVCpl/UVVKBSyx1BUODqKDRAzSlWrBsyeDVy5ooCPT87mi3Z2BkaPBoYNA2xt8y7WjLAuEMB6QGlYF0iNdUHPXbgAAFB07w5FBt8184rc9UCX19U5woCAAM3aFW+Kj4/PsMtRVsaNG4fVq1djw4YNCAkJwWeffYZXr15pZonq27cvpqjnHAXw2Wef4fnz5xg9ejRu3bqFvXv3arpmEQFi3Znx44F794BJk8SYC6VSt02hAJ4+FdPdurkBAweKiR7Onwdk6j1HRERE+iQiQswwAwBvTERU3OncYiFJEhQZrBp4+fJllChRQqdz+fn5ITIyEjNnzsTTp09Rt25dHDhwQDOg++HDh1pZkpubGw4ePIixY8eidu3aKFWqFEaPHo1Jkybp+jaoiCtRApg/X2y6SkoCAgPFDFTXrgHr1okNEJM81KghVgz38ADatgUqV87b2ImIiEjPbdqUdr9BA/ni0DPZTizs7OygUCigUChQuXJlreQiNTUVcXFxGDZsmM4BjBw5MtOuT0FBQenKvLy8dF6Ij0gXxsZA375Anz7AkSNAUJCYSe7iRSAqSiykeeWKSDYMDUVrxqefyh01ERERFZgzZ8TtrFniqiMB0CGxWLJkCSRJwsCBAxEQEAAbGxvNc8bGxnB3d4eXl1e+BEkkB4VCtEi0bSseSxLwzz9pSUZQEHDqFDB4sPj7Mm8e8NbsyURERFQUqZdGaNRI3jj0TLYTi379+gEAypUrhyZNmmhWviYqLhQKMebCzQ3w8RGJxpdfAv7+ovVi505g+nRgxAgxdTUREREVQQ8fArdvi/tVq8obi57RefB2ixYtNElFQkICYmJitDai4kKhAGbOFK0WHh5AbKwYMO7uDsyZAzx+DDx/DvDXgoiIqIi4eBGoXVusU1G+vLjaSBo6Jxbx8fEYOXIkHB0dYWFhATs7O62NqLhp2lTMGLV+vfgbExUlWi5KlQJKlhQzU3l5ialwdVhjhoiIiPTJn3+KgdrR0UClSsChQ9D7lXgLmM6fxoQJE3Ds2DGsXLkSJiYmWLNmDQICAuDq6oqNGzfmR4xEek+pBPr1A0JDgZ9+AmrV0n7+zz/FWAwXF6BrV6B7d7H5+irw+ec2uHNHnriJiIgoC0lJwJYtgKenuEqotno1UKGCfHHpKZ2nm/3tt9+wceNGtGzZEgMGDEDz5s1RsWJFlC1bFps3b0bv3r3zI06iQsHQEOjdW2ypqaIsPFzMSrd2LXDrFvDbb28eoQBgBkmSsHOnDAETERFRxv74A+jZU/RtBsS0kb16idVzPT3ljU1P6ZxYPH/+HOXLlwcAWFtb4/nz5wCAZs2a4bPPPsvb6IgKMfXsc66uYuzFxInA6dNiqlq1Fy9UmD5diV9/BW7cAKpXlydWIiIiesOzZ6KrwePHgJ0dMHYsMHQo4Ogod2R6TefEonz58ggLC0OZMmVQtWpVbN++HY0aNcJvv/0GW1vbfAiRqGhQKIBmzcSmplIBx44l4tgxEzRpAtjb637e1q2BH37IuziJiIiKrZs3xYra330HPH0KWFuL7gY5+QddDOmcWAwYMACXL19GixYtMHnyZHTp0gXLli1DcnIyFi9enB8xEhVpS5ZEo317Bzx5okB0tO7H370LzJjBiSmIiIhy7KefgLlzgZCQtLJq1YDNm5lU6EDnxGLs2LGa+97e3rh58yYuXryIihUronbt2nkaHFFx4OCgwuXLEm7fVrx757c0bSpuq1QBPvoI6N9frNVjbAyYmORtnEREREXWxInAkyfifocO4p9q796AmZm8cRUyOicWbytbtizKli0LANi5cyc+/vjjXAdFVNyULAk4OOh+3LZtYoG+mzfFxZaffhLlRkZiEgv+OhIREb3DwYNpSUVICBe9ywWdEouUlBTcvHkTxsbGqFy5sqZ8z549mDlzJm7evMnEgqgA+fqKaWvPnwc2bAC2bgVevhTr9kydCuzalfXxpUoBAwfybygRERVTkgQEBIj7Q4fyH2IuZTuxuHbtGjp37oxHjx4BALp164aVK1fC19cX165dw+DBg7F37958C5SIMqZQiFnvPD2BZcvEmLPu3YHbt8X2LgsXAi1bAsOHi5ZfrvVDRERF3r17wJw5wJ07wNmzomzQIHljKgKynVhMmjQJFStWxLJly7B161Zs3boVISEhGDRoEA4cOAAz9kEjkp1SCXzwgRhrFhGR9b6SBJw4IdbVCAoS24QJwIIFBREpERGRTCQJGDMmbWEpIyPgq6+Ahg1lDasoUEiSJGVnR0dHRxw6dAh169ZFdHQ07OzssGHDBvTp0ye/Y8xTMTExsLGxQXR0NKytrWWJQaVSISIiAo6OjlDy8nCxpg914eFDYPlykVAYGACXLwM1asgSSrGlD/WA9APrAqmxLuQTSRItE+vWicd9+ojFpvT0H58+1ANdvjtnO8KoqCi4uroCAGxsbGBhYYHGjRvnLlIikl2ZMuJCTefOYrXwHTvkjoiIiCiPSZLo/jRpUlpSsXgxsHGj3iYVhVG2EwuFQoHY2FjExMQgOjoaCoUCr1+/RkxMjNZGRIVTp07idscOIClJ3liIiIjyxK5dQJcugJMTUKGCGFgIAB4eYjVtylPZHmMhSZLWTFCSJKFevXpajxUKBVJTU/M2QiIqEL6+wMyZwI0bwNKlwPjxckdERESUAyoVcO6cuFL2zTdp5UZGQL16YraToUPli68Iy3Zicfz48fyMg4hkVrIkMG8eMHgwsGIF4OfH1byJiKiQkCTgwgWxwNOOHcB/s5hqnD0L1K0LmJrKEl5xke3EokWLFvkZBxHpgZ49xcxQYWFA+fKiFWPCBPG3mIiISO+8fi2mjd28Gbh/P63cykp0gSpfHvDxEV2fKN/leuVtIio6LCyAAweAyZPF9LNbtgCBgWKmqJo15Y6OiIiKrYQE4PRp4MwZIDExrXzLFnE1DBD/xLp0EU3uHTqwdUIGTCyISIunJ3D8OPD330DfvsD168CaNWK8m5GR3NEREVGxEB8P7NwpVno9fx44eVK0TmSmcmXg0iXA3LzgYqR0mFgQUYbq1weGDAFGjwa+/RbYvx9YtEhMS0tERJSnwsKAX34RLRMREcCmTcCLF9r7uLgArVqJQYFvMjUFhg1jUqEHmFgQUaZGjBCL5gUEALduiRbmvXvTpqYlIiLKkeRkMdD68WPgr7+An38Wszm9yd0daNsWqF497VahkCVcyh6dEovk5GSYmZkhODgYNdnhmqjIMzAQyUWfPsDIkeIC0vvvA19+CUyfLnd0RERUqEgSsGePWKjuxAng11+1n2/TBihXDjA0FFewOnUS/4io0NApsTAyMkKZMmW4VgVRMWNtDaxeLVqqT50SLRjduwNVqsgdGRER6ZXERDG167VrIpF408mTYipYNQMDMR2hvT0wYABQu3bBxkp5TueuUNOmTcPUqVOxadMmlChRIj9iIiI9ZGIC/PEH0Lq1mDHqyy+B//0PUCqBihXFZBxERFSMxMWJVofoaOD5c/HP4dQpMU4iMwYGwIcfinERvXsD7dsXWLiU/3ROLJYtW4Y7d+7A1dUVZcuWhcVb3yb+/vvvPAuOiPSLQpE2Fe3mzWIDRHJRsybQqJGYVcrTU3SFZQs2EVERFBUFfP018P33wMuX6Z93cgIaN04/3auxMTBoEMC10YosnRMLHx+ffAiDiAqL1q3FIG71NYTERPE/5soVsa1ZI8otLIDmzYEvvhCJBhERFQHHj4uWhidPxOOKFYE6dUSztpeX+CdRrRoHWRdTOicW/v7++REHERUSRkbpx9v9+6+YZvz8eeDcOeDCBdFCfuCA2Lp3F91ou3QRY/KIiKgQevRIzODx+rVIHubPF3OQK5VyR0Z6Isf/4i9evIiQkBAAQI0aNVCvXr08C4qICpdSpYAPPhAbAKSmAiEhwOLFwPr1Yqzejh1Ajx5ikVReyCIiKkQePBBdn9atE0mFpydw7BjXjaB0dE4xIyIi0Lp1azRs2BCjRo3CqFGj4OHhgTZt2iAyMjI/YiSiQsbAQIy5WLtWTE+u7kEZGChaPHTZjI2B8eNlfTtERMXXH38A9eoBy5YBr16JP+4//sikgjKkc2Lx+eefIzY2FtevX8fz58/x/PlzXLt2DTExMRg1alR+xEhEhVj9+sCuXeL/kImJaM1IScn+lpwMrFghulYREVEBevxYrHT94gXQoAFw+LAYTFejhtyRkZ7SuSvUgQMHcOTIEVSrVk1TVr16dSxfvhzt2rXL0+CIqOgYOBD46CPdEgRJEjNNPXkiLpItWQJUrSrGCnKsBhFRPrhxA/jlF+DQIbHuBACULStaLszM5I2N9J7O/5pVKhWMjIzSlRsZGUH19lLsRERvsLERmy4CA8XK3w8epI3hqFJFtMp7e+d9jERExVJwMPDbb2IF1DcXQq5QQQyYY1JB2aBzYtG6dWuMHj0aW7duhaurKwDg33//xdixY9GmTZs8D5CIirf33gOuXwf8/UWS8eIFEBoKtG0rulm963+dkxPg7p5+s7LK99CJiPKHJAE3bwJPn+bNuX74Adi2La2sTRsxnV/btkD58rl/DSo2crRAXteuXeHu7g43NzcAwKNHj1CzZk389NNPeR4gEZGlJbBokdhevgRmzgSWL09bSyMnSpbUTjTKlRP/PytUEI+NjfMkdCIi3SQliVkvkpNh9OIFYGeXNp2rJAFnzgA//SSuuOQlAwOgVi2xKva0aZxClnJE58TCzc0Nf//9N44cOYKbN28CAKpVqwZv9kkgogJgawssXQp8/jlw7VrW+6pUYnzG/ftAWJi4vX8feP4cePZMbBcvpj9OqRRdiitWTNsqVBC35cuzRwAR5QNJAu7cEV/sr12DEkDJrPY3MRF/mPJi/m4nJ7EmRcOGuT8XFWs6JRbJyckwMzNDcHAw2rZti7Zt2+ZXXEREWapUSWw5EROTlmSok46wMODePfF//fXrtLLDh9Mfn91ZFg0NxQVALy+gSRNx6+ycs5iJqAiJixMDo69fF4Ol1bfq2S2srSE5OyM1NRUGBgbQSh3c3MSKox99JK60EOkRnRILIyMjlClTBqlvDuohIipkrK2B2rXF9jZJEt2W79xJ2+7eFbe3b4ukJD4++691+rTY1MqVE0nGyJFA48a5fy9EVAhIEnDrFpCQACQmitVCw8LS72doKFoNNm6EVL48oiIi4OjoCAW7JVEhoXNXqGnTpmHq1KnYtGkTSpQokR8xERHJRqEAXFzE1ry59nOSJLpPxcZm71zx8aKr9JkzwNmzouuWuiXk1CnRQsLvC0TFwOzZYnDYmxwdgZYtgerVxboQ1auL/pbqAV6caZMKoRwN3r5z5w5cXV1RtmxZWFhYaD3/d25GUxIR6TGFArC3F1t21agB9Osn7kdHA+fOAb6+Yvrc/fuB99/Pn1iJSGZRUWL61t27gV9/TSt3cQHKlAHWrhXJBFERonNi4ePjkw9hEBEVfTY2QLt2wKefihmuBg5UL2CrQFKSHYyN82AQZjYZGABjxjCxIcoTjx+Lbk6AuIKwahWwYYPo9qRWrZpoqmRvDyrCdEosUlJSoFAoMHDgQJQuXTq/YiIiKtJmzAA2bxZjOSIiAEABwKTA47hzB3j1SrvMxESsHWJnV+DhEBVOixcD//tfxs/VqSNmeerWTQzqyosZnIj0mE6JhaGhIRYuXIi+ffvmVzxEREWejQ3w559iAwCVSoWYmBhYW1tDWQCDLlQqYNAgMSOWn1/65w0NgdatxaQzPj6iKzgRvSU+Hti0STupsLQUyUOLFsCkSUCzZvLFRySDHK28feLECbi7u+dZEMuXL8fChQvx9OlT1KlTB9999x0aNWr0zuMCAwPRs2dPdOvWDbt3786zeIiI8lvZsmIDxBf9iIgEODpaF9hg7tevxRpbkqRdHh4OhIQAhw6J7bPPxCxaah06iF4eNjYFEydRgUtNFYvdZEaSxJiJadPUTY7il+TECaBu3QIJkUhf6ZxYdOzYEZMnT8bVq1fh4eGRbvB2165ddTrftm3bMG7cOKxatQqenp5YsmQJ2rdvj9DQUDhmcZns/v37GD9+PJq/PW0LERG908CBYsvIrVvAzz+L7eJFsdq5WmAgEBwsWjLepFSKlg31jFrq7a1/EUT6R6US60gcPw4cOyYShDcrfVbKlQNGjQIGDGC2TQRAIUlvX6/KWlbN9AqFQuc1Ljw9PdGwYUMsW7YMgOgS4Obmhs8//xyTJ0/O8JjU1FS89957GDhwIE6ePImXL19mu8UiJiYGNjY2iI6OhvWbl+EKkEqlQsR/c1MXRLcH0l+sCwTodz14+lSMRQWAR4+A/v2Bf//N/vHW1umTDVdX7cemplmfw9BQ7GdgkOO3UWjoc10oEiRJzPf8119p28WLYoEaXdjbA1OnAiNGpE0Pm8dYFwjQj3qgy3dnnVssVHk4r3JSUhIuXryIKVOmaMqUSiW8vb1x9uzZTI+bNWsWHB0dMWjQIJw8eTLL10hMTETiG7MyxPz3x0OlUuXpe9GFSqWCJEmyvT7pD9YFAvS7Hjg6po2xqFQJOH8e+OEHIDpaexBqcrLoRvX0KfDkiZgk5/VrBWJixHe20NDcxWFsLKFiRRFD5cpApUoSKlcW9x0di86YWH2uC4WOSgUkJYkFY06cgOLECXGr7r70BsncHGjaFFKrVkCrVkC9eiKjzYq60uXTz4p1gQD9qAe6vLbOiUVeioqKQmpqKpycnLTKnZyccPPmzQyPOXXqFH788UcEBwdn6zXmzZuHgICAdOWRkZFIUE8NV8BUKhWio6MhSRKvQhRzrAsEFK56oFQCw4a9ez9JAmJjFQgPVyIiwgDh4cr/NgNERGjfT07OOitISQGSkhS4cQO4cUNdmnaMtbUK5cunoG3bRAwf/uqdLSD6rDDVBb2SlATDW7dgdO0aDK9ehdHVqzC8fh3K+Ph0u0rGxkiuXh0pdeoguU4dJNeujZTKlQEjo7SdXrwowOAzxrpAgH7Ug9jsrgoLHRKLTp06YevWrbD5rw/h/PnzMWzYMNja2gIAnj17hubNm+NG2l/9PBcbG4s+ffpg9erVsM/mClVTpkzBuHHjNI9jYmLg5uYGBwcHWbtCKRQKODg48I9FMce6QEDRrQdOTmIh4XfLukduairw6JGE0FDg9m3g1i0Fbt0S9x88AGJilAgONkZwsDF+/90SU6dKMDd/96saGIjZr/RpHEhRrQv54tkz4Pffodi9Gzh8GIrXrzPcTTIzA5o0gdSihZhLuVEjGJqYwBCAPuegrAsE6Ec9MNXhak22E4uDBw9qdSmaO3cufH19NYlFSkoKQnVs67a3t4eBgQHCw8O1ysPDw+Hs7Jxu/7t37+L+/fvo0qWLpkzdPGNoaIjQ0FBUqFBB6xgTExOYmKSfH16pVMr6i6pQKGSPgfQD6wIBrAdZUSqB8uXF1rGj9nMJCcDdu2JF86lTgZAQBfr0yX6/KHX3rv/+lekF1oUMSJLIIs+cEdvZs2IWgTe7aNjYAPXri61ePXHr6gqFuTlgZITC2FuOdYEA+euBLq+b7cTi7THeOo75zpCxsTE8PDxw9OhRzYreKpUKR48exciRI9PtX7VqVVy9elWrbPr06YiNjcW3334LNze3XMdERESFh6mpWL28Rg0xU9XMmcClS9k7Vt0CUqYMYGaW+X5GRkCtWkCjRkDDhmJ7qwcv5VZCAvDPP2KGgMy2jGZqqlNH/OC7dRNTvRaVwTZEhZSsYywAYNy4cejXrx8aNGiARo0aYcmSJXj16hUGDBgAAOjbty9KlSqFefPmwdTUFDVr1tQ6Xt1i8nY5EREVLyVKAP9NMJgtly6JrlAvXwLv6kL877/AgQNpj8uUEYlGkyZAly7Z7fJVTC1YAMyeLfq0ZUSSxMIq72JoKFoimjYVH7yXF1C6dN7GSkS5ku3EQqFQQPHWlYC3H+eEn58fIiMjMXPmTDx9+hR169bFgQMHNAO6Hz58yCZAIiLKc/Xqid41Dx5kvV9cnEhCzp8HLlwQCwg+fCi2nTuBceOAmjXFhXMfH3Fe/tsC8OqVWAxl0qTs7W9uDri5aW+lS6fdL1cO2Ro8Q0SyyfY6FkqlEh07dtSMV/jtt9/QunVrzQJ5iYmJOHDggM7rWBQ0rmNB+oR1gQDWg8ImJkYsfXD+PHD4MBAUpH0x3sICqF1b9MxRbzVrZu87cZGoCxERwKJFwMqV2k1Bx48D7u4ZH2NtDdjZsSvTG4pEXaBc04d6kC/rWPTr10/r8SeffJJun759+2b3dERERIWStbVY6qBVK3Ex/sULYO9eYNcu4OBBcaH+7FmxqSmVQJUqQJs2wPz5+jUTVZ6RJGDFCmDiREA9zWv58mJV6n79RKsDERVp2U4s1q1bl59xEBERFUp2dsAnn4gtJUUMCA8OTtsuXQIiI0UXqpAQsdjzypVicb9C37Pn1i2xYuLu3SKjevpUlDdoIEbSv/8++4URFSOyD94mIiIqKgwNgWrVxNazpyiTJPF9+9QpYOhQ4M8/xTgMQFzEV68gXrmyGAReooQB7OyADGZK1x+pqUCfPsDWrdrlSiXw1VfA//7Hbk1ExRATCyIionykUAAuLkD37kD16sCoUaIV48WLtJlUjx5V760E4ABDQwlOTtrfzZVKwNtbdL+qXFmGN/KmefNEUqFUisVFPv1UTJXl4iI2IiqWmFgQEREVkBo10pKIZ89ETyLtTcKtW0BCggL//pv++LVrgfXrAV9fMVajbFkdA7h1C9i8Gfj99+xN8ZqRpCSxKiEggunTJ2fnIaIih4kFERGRDEqWFEsxeHmllalUEp4+jUBSkiOePdMem/D8ObB0qcgJAgOB06eBffsAZ+csXuTff8UcuRcuiNHkV68AAMwRD3PkMLFQa9dODCwhIvoPEwsiIiI9olSKXkUZzczatq0YEN6zJ3DzplgRPGul/tt8tEpNjFJx6rtgNKgal/MgPTw4joKItDCxICIi0meSBBw6BHzzDfDiBeoCOGZaEj4WC3D+Vc0cnTIx2QB7/vFAg6F5GikRFXNMLIiIiPRRQgLw66/AqlVicbk3uAD4E/sh4b8WAzNzMaWUeqtQQSyc4eGRbtGMNWvE7FTLlwOtW4tdnZz0fBYqIioUmFgQEREVlKgo4LPPxMIWGVAAKJGUBIWREXDlCvDypXjC2BgYMUJkAm/sq7CxEYmEs3O2uyX16ycGgZ87p3U62NiIBMPZWdxmtZmZ5eztE1HRxsSCiIhILSZGzHh0507a9ugRoFLlzfnT5pXNkAKA8ZsFpUuLTODTTzMedJEDJibAL78AQ4YAly8D4eFAcjIQHS22W7fefQ5ra6BEibS178zMgPHjgf798yREIiqkmFgQEVHxJUnAihXAli0iiYiIKJjXLVsWWLgwXbFKpUJ0TAxsrK2hdHEBmjYFDAzy/OVdXcXsUoD4CF68EAnGu7anT8VsszExYnvTgAFiEcAqVdK/npmZGGhepw5ga5vnb4eI9AQTCyIiKn6ePQNmzAAuXgTOn9d+zsFBe7xC2bKAkVHevbapKdC+fbqxDwAAlQqJERGAo2Nac0A+UyhE60OJEmLF8KxIkmjVCA8XyYgkifJ9+4DZs4Eff3z365UtC9Stq72VLcsJpoiKAiYWRERUfPzzD3DwILBoERASIsqMjMS3Ym9vMZLZxkbeGPWYQiFaHN5udfDyAho1AnbtyrjX2IsXotvVgwdp2549ac8PGQJ8/31+Rk5EBYGJBRERFU2RkcCNG8D162L74w/g2rW0511dAX9/4L33gKpV5YuziOjSRWxZefFCjEkPDhbbX3+JH8nmzUBsbN7H1KED0Ldv3p+XiDLGxIKIiPTPs2eim9Jff4mWhdTU7B0nSWIgwPXrGc+8pFAAnp7iG+fgwSK5oAJjZwe0aCE2QCQazs7Aq1fA1q15/3pbt4phM40a6XachQVQvTpnvyLSFRMLIiKSlySJmZgOHgSCgkQycf9+3py7XDnxDbFGDaBePbF0dcmSeXNuyjU7O+DIEeDvv/P+3JcvA+vWARMm5Ox4pRKoXFmMAalTJ21zceF4EKLMMLEgIqKciY8HLl0Sg591aVV4U2IicOYMEBaW/rlKlYAGDYDatcWA5+yysxOJRLVqGQ+QJr3SvLnY8ppKJcbAvzmWI7uiosR286bYAgPTnrO3F1PrTpqUd7ESFRVMLIiIKGuSJFoR9uwR37YSEkQH+WvXcpZMZMTISEyt2rYt0LgxUL8+5yWlXFEqgfnzxaYrdY+6y5dFVb98WWyhoeJXICBArAeibrlwdgY+/jhfZgYmKlSYWBARkRiPcOiQaH14M1lISAAOHwYePsz4OGdn0YG9Tp2cdUhXKICaNYGWLQFLyxyFTpTXFArR5cnFRQzHUXv9WjSEPXgAjB2rfcy+fUDHjgUbJ5G+YWJBRFRcvXoFLF0q5gj966+0RQkyYmEBdO4sxisYGopvVw0bAqVKscM5FRtmZmKtjnXr0qbVVQ86X7gQsLISDW/8laDiiokFEVFxc/26+Da0YYNY10GtTh2gVSvx7UhNoRCjVzt04BQ5RADatBGbmpcXMGoUcPy4GCtSp474lVEqgT59xK8UUXHBxIKIqLgIDAQWLwYuXEgrK1sWmD4d6NSJU68S5cDnn4uEYvlysR6HejwGIFo2zp4V962tRUMfWzOoKGNiQURU1D19Kro8zZsnHhsais7gbdqItRzMzeWNj6iQq1sXWL0aWLAA+Pln4M8/RZcpQLRoqPn6AmvXcrIyKrqYWBARFSUxMWLq1nv3gNu3gf37gRMn0sZP/O9/wMSJYh5OIspTdnbAp58CAweKlcT/+ivtuYcPge3bxfS1x45xORUqmphYEBHpO5VKJAzPngHPn4tb9f0nT9ISiXv3xFyYGWnYEBgxAujXr2BjJyqGlEpg2zbtslOnxJS0V64AH30E+PllfryLC9ClS/7GSJQfmFgQEcktNBQWGzdCkZAAvHihnTg8eybKdFkvwt4eKF9ebA0aiG8zZcvmX/xE9E7NmomZmxs3Fo2IJ05kvf+BA2IgOFFhwsSCiEhOz55B0akTrO7ff/e+5uai/0SJEuK2ZEnAwSEtiShfHihXTowSJSK9U6uWSBhWrACSkzPe5/Jl4M4dYOhQBZo0sUaLFmKcRo0aXICP9B8TCyKi/PDwIXDwIJCYmPHzkiSmi/n5ZyiSkkTRmDFQODtrJw9v3pqaFuAbIKL80Ly52DLz+++Ajw/w4IECDx6Ya9bJsLQUDZAVK6bNLFW7tpjS1sYm38MmyhYmFkREeUWSgJ07ge+/F6Mzs1pw7s3D6tbFc39/2HXtCoVSmc9BEpE+69wZePwYOHNGhePH43HtmgUuXFAgNhYIChLbmyZPFsnFqFFiOlsiOTGxICLKC6mpwJQpYvldtebNxSjMzDg7A337QqpbF8mRkfkfIxEVCo6OQNeuQOPGcXB0NIckKRASApw7J2aPBoCkJDG17fXrwKpVYs2MM2eA+vXljZ2KNyYWRES5kZoK7N0LTJ0q/sMDwIQJwPDhgLt79s6hUuVbeERU+BkYADVriu1NX3whWjCmTRM9KydNEj0w2fBJcmHVIyLKifBwseBcxYpAt24iqbCzA9avF6tkZTepICLKIYUCaNUK2LgRMDICjhwRs0qfOgX8849uk8kR5QW2WBBR0SJJQGSkGDz9zz+iv0BeSk4Gfv0V2LUrbVoXW1tg6FBxudDOLm9fj4joHSpWFCt99+0rukWtWiXKjYyAMmXEbNPu7tpb2bJAqVKcaYryFhMLIio8JAk4fVokDpIEREeLBEK9PXgAPHoEJCQUTDyNGwPDhgG+voCZWcG8JhFRBvr0Ebdr16b9KUxOBu7eFVtGLC2BPXuA1q0LLk4q2phYEJF+u3ZNrCwNABs2iNGK2eHiAri55c8X/urVgSFDgLp18/7cREQ51KdPWoKRkiJml3rwALh/P21TP37wAIiLE2MymFhQXmFiQUT6JTkZiIoSrRK//y5GJb7JyAho2FB0Lra0FO35ZcqIzc1N3JYqBZiYyBM/EZEeMDRM+9OY0boZX30lpqr9/nvAw0Ms3le1atoaGUQ5wcSCiPJfUpKYsiQhAXj+XCQNERHi9u37L16kP75OHZEo2NkBAQGAp2fBvwcioiJkyBDgt99E71I/P1HWsaMYQmbIb4eUQ6w6RJRzKhXw99/Aq1dZ77dsmVg4LruUSsDeHnBwAPr3B/73P15GIyLKQ3Z2Yh3PiROB7dvFdZ39+0UPT11W8jYwEN2vPv2Uf6aJiQUR5YQkAdu2AXPmiDEQ2VW5suim5OCQtjk6pr9vZ8epSoiI8pmxMbBkidh+/hn4+OO05Xh0cfKkaJRu1SqtzMwMaNo06zVCqehhYkFEugkNBfz9RWIBiHEObm5ZH6NUAmPGiEtaRESkdz76CAgJAW7e1O244GBg1iyx8ve6demfr1EDaNMG8PYGWrQArK3zJFzSU0wsiEjb69fiv0tG9uwR/0HUJk0So/9sbQskNCIiyj9Vq4pNFz4+gJcXsGIFkJiYVh4RIZKO69fFtnSpuMZkaiqSjFWr2JpRFDGxICqu4uKAe/fSl3t6Zm8diCFDxMrT7FRLRFSstW8vtrc9ewYcPy5WBD9yRKynER8vBoifPg306gUMHMiZu4sSJhZExU1CArB8OTB7NvDyZeb7GRkBTk7pyytVEi0XVlb5FiIRERV+JUuKcRsffyweh4eLNTSGDgUuXwa++w5YvRpYswZo0EC0ZpQpw+tVhZlS7gAAYPny5XB3d4epqSk8PT1x/vz5TPddvXo1mjdvDjs7O9jZ2cHb2zvL/YkIYrB1aKiY/sPNDRg/XiQVtrYieXh78/MTbdqPHqXfjh1jUkFERDpzchKN4ufOickCGzYU17o++UR0wXJ3B95/H3jyRO5IKadkb7HYtm0bxo0bh1WrVsHT0xNLlixB+/btERoaCkdHx3T7BwUFoWfPnmjSpAlMTU3x1VdfoV27drh+/TpKlSolwzsg0kNPnoi/3H/9JbaLF8Wic2qlS4uxEn37cvYlIiIqUCYmwIgRouVi0iRgwwYgNRWIjRVT3taoIZYvAsTifgEBbMUoLBSSJElyBuDp6YmGDRti2bJlAACVSgU3Nzd8/vnnmDx58juPT01NhZ2dHZYtW4a+ffu+c/+YmBjY2NggOjoa1jJNTaBSqRAREQFHR0colXrRaEQyyXVdSE3VvrQTFwcsXCim5nj7V9vQEGjbVvwlf/99roCkR/g3gdRYF0itONaF69fF9a6//9YunzMHqFBB3Dc2Fv/KLC0LPj456EM90OW7s6zfLJKSknDx4kVMmTJFU6ZUKuHt7Y2zZ89m6xzx8fFITk5GiRIl8itMIv2TmChGv02YADx4kPE+tWuLdmYPD9F5tVYt0YGViIhID9WoAfz5pxjorW69WL8emDZNe7+hQ8WsUqR/ZE0soqKikJqaCqe3Bog6OTnhZjYnUp40aRJcXV3h7e2d4fOJiYlIfGP+s5iYGAAiA1SpVDmMPHdUKhUkSZLt9Ul/qFQqIDoaKoVCTPMaGSm2iAggKgoK9eO3NkVsrOYckoGBdncmT09I8+cDjRtn9IIF8K5IV/ybQGqsC6RWXOuCgUHaDFMdOwIqlQKPHonHr18Df/6pwPbtEhYulGBhIV+cBUUf6oEur12o+0LMnz8fgYGBCAoKgmkmV2LnzZuHgICAdOWRkZFIyM6UmvlApVIhOjoakiQVm+bNYi0+HoqUFM1DgydPYHTuHIzPn4fR+fNwVv/F1JHK0hLxQ4cibvhwwNw8/Q4RETmNmAoY/yaQGusCqbEuCF99lXY/NRVo3NgB//xjgP/9Lx6zZsVmfmARoQ/1IDY2+5+zrImFvb09DAwMEB4erlUeHh4OZ2fnLI/9+uuvMX/+fBw5cgS1a9fOdL8pU6Zg3LhxmscxMTFwc3ODg4ODrGMsFAoFHBwcivUfi0JBksQlkpwcd/w4FCtWQHHwYPYOMTICHBwAR0fA3l7cd3CA5OCQ9tjRUVMOW1uYKxTIIKWgQoZ/E0iNdYHUWBcy9v33Ypjg2rXmCAgwy3BW9KJEH+pBZhfvMyJrYmFsbAwPDw8cPXoUPj4+AMQHePToUYwcOTLT4xYsWIA5c+bg4MGDaNCgQZavYWJiAhMTk3TlSqVS1l9UhUIhewyUCUkCkpOBlBQxL961a3l7fjMz0U2peXOomjRBZPnycChXDkoDgwynveBEGMUD/yaQGusCqbEupNepE9CoEXD+vAJff63AokVyR5T/5K4Huryu7F2hxo0bh379+qFBgwZo1KgRlixZglevXmHAgAEAgL59+6JUqVKYN28eAOCrr77CzJkzsWXLFri7u+Pp06cAAEtLS1gWlykCKP9cuwZ8+CFw+3buz2VnJ5YUHTpUrPijZmQEqH9JVSpIERHiMefSIyIieqdx44AePYDFiwEXF7E0E+kH2RMLPz8/REZGYubMmXj69Cnq1q2LAwcOaAZ0P3z4UCtTWrlyJZKSkvCxehnH//j7++OLL74oyNCpsIuKAq5cAa5eFduVK2J7Y7A/AMDbG9i9W/fzm5pyjQgiIqI85ucnJkScNAmYPBlo3RqoX1/uqAjQg8QCAEaOHJlp16egoCCtx/fv38//gKhoi4oSf5WOHcv4+WbNgK1b0ybJtrFhawIREZEemThRrP+6YwfQv7+4b2wsd1SkF4kFUYFJSgJatABu3BCPK1QQ6zvUri1ua9UCKldmIkFERKTnli0Djh8XnQ6cnERP4zcZGYnlnsaMkSW8YomJBRUPsbHAwYPAkiUiqTA2BoKDgWrV5I6MiIiIcsDREVi5EujeHXj5MuN9xo4F7t0DbG3TP2dmBgwZApQsmZ9RFi9MLKjoW7wYmDJFtFao9evHpIKIiKiQ+/hj4OFD4L/1j7V8/z3w3Xdiy8zNm8CGDfkXX3HDxIKKti1bgP/9T9yvVAno2lVszZrJGxcRERHlCTe3jMsXLwaqVBHJw9sSEoA1a8TXhFmzgLJl8zfG4oKJBRUNf/8NrFoFHDgg1qBQi4wUt//7H/D11/LERkRERAXO0BAYMSLz5+/dE/O4dOkiEouqVcVK31w2JOeYWJD+S0wETp0C9u8HjhwBnj/Xfj4lBXjyJPPjfX3FXwoiIiKi/0ydKhIL9azzv/8uWj/Kl0+/r5ER8N57YlwGZY6JBemX4GBg1Cjgn3/SyiIigFevsj7OyEh0tBwwQIzmUjM15SxPRERElE6bNsDhw2KMxqJFYm6X0aMz33/gQODHHwsuvsKIiQUVLJVKdGj89tuMR1rdv689yFrNyQno0AHo2BGoWDH98+7unNaBiIiIdOLtLW6dnYE5c7R7U6s9fw7cvQuEhRVsbIUREwsqOImJIjE4fjzr/Tp3FrM4qTs5WlmJGZzY6ZGIiIjyQadOYsvIrl3Ahx+K5GL2bKBOHTEug9JjYkEFZ//+tKRizhzRWfFtVlZisTp2XSIZqVQqJGXUcpaPr5ecnIyEhAQomUAXa3LVBSMjIxgYGBTY6xEVJuo1MB4+BGbMEF9RLl8GXF3ZWeJtTCwo/4WHi3UjDh4Uj4cOFSOmiPRQUlISwsLCoFKpCuw1JUmCSqVCbGwsFEyqizU564KtrS2cnZ1ZB4ne0ry5aKl4+FC0XkRGimugADBsmFikjwQmFpR/9u4F5s0DTp/WLm/VSp54iN5BkiQ8efIEBgYGcHNzK7ArxpIkISUlBYaGhvxSV8zJURckSUJ8fDwiIiIAAC4uLgXyukSFhaEhMG2auF+unLivvva0aZNY9XvCBKB+fdlC1BtMLCh/hIYCgwaJ1gq1SZMAHx/A01O2sIiykpKSgvj4eLi6usLc3LzAXpeJBanJVRfM/ptDMyIiAo6OjuwWRZSJyZPF15noaMDFRUxaGRgopqo9dw6oXl3uCOXFzryU9zZuFG2E4eFiutcdO4ALF4D584HGjTl+gvRWamoqAMDY2FjmSIgKnjqZTs5oWhwi0lAoxLiLoCDA31+UxcUBNWpot2YUR2yxoLwTHw/s2yfGUwBAy5bAhg1AmTKyhkWkK7YaUHHEek+kG09PoGFD4NEj4Px54No1YO5c4NYtcY21OC6mx8SCci8pCRgyRLRMxMeLshIlgEOHxMJ1REREREWQUpm2aN6GDcDgwcDOncCDB2LdXgAoVQro1at4dNhgVyjKvaNHxW9TfDxQtizwv/8Bp04xqSAqpn788Ue0a9cuV+e4f/8+FAoFgoOD8yaofPbFF1+gbt26coeRoaioKDg6OuKff/6ROxSiIq1fP+DIEXFt9cIFMRZj0iTgk09EsvHqldwR5j8mFpQ7qanAnj3ifo8eYlnKr78WC9oRUb7r0qULOnTokOFzJ0+ehEKhwJUrVzRf1A0MDPDvv/9q7ffkyRPNYOH79+8DSPtin9H2559/ZhpPQkICZsyYAX91x+MccnNzw5MnT1CzZs1cnaegjB8/HkePHs3WvvmZhPTv3x8+Pj5aZfb29ujbt2+ufyZE9G7vvSe6RX3+OdC/P9CokSj39RVLdX33nazh5TsmFqS7mzfFjE9+fmKk0vffi/JOnYpHOx+RHhk0aBAOHz6c4dXodevWoUGDBqitnnAdQKlSpbBx40at/TZs2IBSpUpleP4jR47gyZMnWpuHh0em8ezcuRPW1tZo2rRpDt+RYGBgAGdnZxga5n+P3bxYDNHS0hIl9XilrAEDBmDz5s14/vy53KEQFXkVKgBLlwLr1gFbtqQtoidJwKhRIvEoqphYkG7CwwFvb2DtWmD7djGtrJ2dmPHpk0/kjo4ob0mSaLuWY5OkbIXYuXNnODg4YP369VrlcXFx2LFjBwYNGqRV3q9fP6xbt06rbN26deinnnThLSVLloSzs7PWZpRFN8fAwEB06dJFq0x9FX3u3LlwcnKCra0tZs2ahZSUFEyYMAElSpRA6dKlteJ6uyvUixcv0Lt3bzg4OMDMzAyVKlXS7J+UlISRI0fCxcUFpqamKFu2LObNm5dpjOp45syZA1dXV1SpUgUA8OjRI/j6+sLW1hYlSpRAt27dNC04ABAUFIRGjRrBwsICtra2aNq0KR48eAAgfStEZvuuX78eAQEBuHz5sqYFSP2zW7x4MWrXrg1bW1uUKVMGw4cPR1xcnOac69evh62tLQ4ePIhq1arB0tISHTp0wJMnTzQxbNiwAXv27NGcOygoCABQo0YNuLq6YteuXZl+LkSU9ypUEF+dbt9OK1OviVEUMbEg3Xz7LfDvv0ClSiIdX7tWdH+aNImtFVT0xMcDlpb5vimsrGBkZweFlVVauXoihHcwNDRE3759sX79ekhvJCM7duxAamoqevbsqbV/165d8eLFC5w6dQoAcOrUKbx48SJdMpBTp06dQoMGDdKVHzt2DI8fP8Yff/yBxYsXw9/fH507d4adnR3OnTuHYcOGYejQoZmOA5gxYwZu3LiB/fv3IyQkBCtXroS9vT0AYOnSpfj111+xfft2hIaGYvPmzXB3d88yzqNHjyI0NBSHDx/G77//juTkZLRv3x5WVlY4efIkTp8+rfninpSUhJSUFPj4+KBFixa4cuUKzp49iyFDhmQ4k1JW+/r5+eF///sfatSooWkB8vPzAwAolUp8++23CA4Oxvr163Hs2DFMnDhR69zx8fH4+uuvsWnTJvzxxx94+PAhxo8fD0B0x/L19dUkG0+ePEGTJk00xzZq1AgnT57M8nMhorxnYABUrAjMmiUev5lkFDWcFYqyb8sWYOFCcX/SJNEdiohkN3DgQCxcuBAnTpxAy5YtAYhWiI8++gg2NjZa+xoZGeGTTz7B2rVr0axZM6xduxaffPJJpq0QTZo0SbcC+ZtX0d/08uVLREdHw9XVNd1zJUqUwNKlS6FUKlGlShUsWLAA8fHxmDp1KgBgypQpmD9/Pk6dOoUePXqkO/7hw4eoV6+eJml5M3F4+PAhKlWqhGbNmkGhUKBs2bIZf1BvsLCwwJo1azRrlvz0009QqVRYs2aNJllYt24dbG1tERQUhAYNGiA6OhqdO3dGhQoVAADVMhlLFhMTk+W+lpaWMDQ0hLOzs9ZxY8aM0SyQV7FiRcyePRvDhg3DihUrNPskJydj1apVmvOOHDkSs/77tmJpaQkzMzMkJiamOzcAuLq64tKlS+/8bIgof3z6KTBzppieNjm5aM5xw8SChOfPRdLw8mXGz8fEiMQCALp3B/r0KbDQiGRjbi5WPcpnGa62rMPK31WrVkWTJk2wdu1atGzZEnfu3MHJkyc1XzjfNnDgQDRp0gRz587Fjh07cPbsWaSkpGS477Zt2zL9Av22169fAwBMTU3TPVejRg2tBMXJyUlrYLaBgQFKliyJiIiIDM/92Wef4aOPPsLff/+Ndu3awcfHR3M1vn///mjbti2qVKmCDh06oHPnzu+clapWrVpaCyFevnwZd+7cgZWVldZ+CQkJuHv3Ltq1a4f+/fujffv2aNu2Lby9veHr6wsXF5d05y5RokS2933TkSNHMG/ePNy8eRMxMTFISUlBQkIC4uPjNYvXmZuba5IKAHBxccn0M3ubmZkZ4rPZEkZEec/ZWawbnJAAzJkjFtmrXFkMUS0qmFiQsGaNGCfxLmPGAIsWiYmbiYo6hQKwsMj/15EkICUFMDTMcZfCQYMG4fPPP8fy5cuxbt06VKhQAS1atMhw31q1aqFq1aro2bMnqlWrhpo1a2Y6raubmxsqVqyYrRhKliwJhUKBFy9epHvu7RYRhUKRYZkqkyVrO3bsiAcPHmDfvn04fPgw2rRpgxEjRuDrr79G/fr1ERYWhv379+PIkSPw9fWFt7c3du7cmWmsFm/9XOPi4uDh4YHNmzen29fBwQGAaMEYNWoUDhw4gG3btmH69Ok4fPgwGjdunO4YXfYFxJiSzp07Y9iwYQgICICDgwNOnz6NQYMGISkpSZNYZPSZSdkcj/P8+XPNeyGigqdQiJ7kV68CAQFp5VevAoVkArx3YmJBYgD20qXivpFR5qOK6tYFunUrsLCIKPt8fX0xevRobNmyBRs3bsRnn32W5UrKAwcOxPDhw7Fy5co8i8HY2BjVq1fHjRs3cr2ORUYcHBzQr18/9OvXD82bN8eECRPw9ddfAwCsra3h5+cHPz8/fPzxx+jQoQOeP3+OEiVKZOvc9evXx7Zt2+Do6Ahra+tM96tXrx7q1auHKVOmwMvLC1u2bMk0WchsX2NjY6Smpmrte/HiRahUKixatAgqlQqGhobYsWNHNj+ZNBmdW+3atWuarnJEJI9vvxXDU1UqseZFRARQq5ZYY1i9oF5hxsSiOLpwQaw9kZwMPH4MbN0q1qOoUAE4exbgFS2iQsfS0hJ+fn6YMmUKYmJi0L9//yz3Hzx4MLp37w5bW9ss93v27BmePn2qVWZra5thdycAaN++PU6dOoUxY8boEP27zZw5Ex4eHqhRowYSExPx+++/a7poLV68GC4uLqhXrx6USiV27NgBZ2fnd763N/Xu3RsLFy5Et27dMGvWLJQuXRoPHjzAL7/8gokTJyI5ORk//PADunbtCldXV4SGhuL27dvo27dvunOFhYVlua+7uzvCwsIQHByM0qVLw8rKChUrVkRycjK+++47dOzYEefOncOqVat0/pzc3d1x8OBBhIaGomTJkrCxsYGRkRHi4+Nx8eJFzJ07V+dzElHeadVKbIDoKDJlirjfvbv4KlbYO4QwsSiK7t0DNm4UicPb/vwTOHYsfXmXLmLVFiYVRIXWoEGD8OOPP6JTp04ZDqB+k6GhoWZWpax4e3unK9u6dWuGA6zVMagHOr89cDw3jI2NMWXKFNy/fx9mZmZo3rw5AgMDAQBWVlZYsGABbt++DQMDAzRs2BD79u1LN+g8K+bm5vjjjz8wadIkfPjhh4iNjUWpUqXQpk0bWFtb4/Xr17h58yY2bNiAZ8+ewcXFBSNGjMDQoUMzPFdW+3700Uf45Zdf0KpVK7x8+RLr1q1D//79sXjxYixYsABTp07Fe++9h3nz5mWYuGRl8ODBmsHmcXFxOH78OFq2bIk9e/agTJkyaN68uU7nI6L8M3ky4OiYNhfOyZNAJj1YCw2FlN3OmUVETEwMbGxsEB0dnWVzd35SqVSIiIiAo6OjTv/4ss3VFfhvXvMMGRqK9rZSpcQcaF27ArlczIpyJt/rAukkISEBYWFhKFeuXKZX5PNDhoO3C7Hu3bujfv36mKK+FEfZll91oXHjxhg1ahR69eqV6T5y1X/KGP8/FA9JSYCJibj/1VfAWzNM60U90OW7M1ssCqu9e4GMpg28ciUtqfjgA6BMGe3n7ezEGvPZmI6RiCgnFi5ciN9++03uMOg/UVFR+PDDD9OtaUJE8jM2BlauBD77TExFO26cuP5bWBXi0Iuxf/8VrQyZzJ4CALC3B375peBiIiL6j7u7Oz7//HO5w6D/2Nvbp1toj4j0x4ABIrFITASuXRNz5RRWTCwKowsXRFLh7CwSjLcZG4saSkRERER6zcREjLWIiBBbYcbEorBRqcSaEwDQsSPw/ffyxkNEREREudKqFbBtG9C+PRAeLhKNwoijgQqbEyfE+ApjY2D0aLmjISIiIqJcenNeha++ki+O3GJiUZgkJqZNF9C1K1CnjrzxEBEREVGude0KtG4t7n/3nbyx5AYTi8LgyhVg3TqgZ0/gr79E2YAB8sZERERERHnmiy/EbXIycPSorKHkGMdY6KvXr4HffgPCwoCpU7VngBozRoyvICIiIqIioVkz0Rnl8mXA2xu4fx9wc5M7Kt0wsdBHCQnAe++ltU4AQMOGYiTPhx8CAwfKFxsRERER5TmFQqy+rV6DbtQoYNcueWPSFbtC6SM/v7SkonZtwN8fOHsW+P13JhVEJJukpCRUrFgRZ86c0ZTdvHkTjRs3hqmpKepmc/L1L774Itv76gt9jjkqKgqOjo74559/5A6FiHLJyiptxYBffwXu3ZM3Hl0xsdAXycnAoUMiifj1V1E2dqxoD/viC8DAQNbwiEh/9e/fHwqFAsOGDUv33IgRI6BQKNC/f/90+7+9dejQIcvXWbVqFcqVK4cmTZpoyvz9/WFhYYHQ0FAcLaydgrNh/Pjx2X5/+ZmE9O/fHz4+Plpl9vb26Nu3L/z9/fPlNYmoYE2dmnZ/yxb54sgJJhb6okcPMXnxrFlpZYsXyxcPERUqbm5uCAwMxOvXrzVlCQkJ2LJlC8qUKZNu/w4dOuDJkyda29atWzM9vyRJWLZsGQYNGqRVfvfuXTRr1gxly5ZFyZIl8+4N5aGkpKRcn8PS0lJv3x8ADBgwAJs3b8bz58/lDoWIcql0aWDRInH/xg2FvMHoiImF3KKigGHDgF9+SSt77z1g3z75YiIiAIAkAa9eybNJkm6x1q9fH25ubvjljb8lv/zyC8qUKYN69eql29/ExATOzs5am52dXabnv3jxIu7evYv3339fU6ZQKHDx4kXMmjULCoUCX/w3pcmkSZNQuXJlmJubo3z58pgxYwaSk5MzPXdQUBAaNWoECwsL2NraomnTpnjw4IHm+T179qB+/fowNTVF+fLlERAQgJSUlEzPp76qP2fOHLi6uqJKlSoAgEePHsHX1xe2trYoUaIEunXrhvv372crjrdbITLbd/369QgICMDly5c1LUHr168HACxevBi1atWChYUF3NzcMHz4cMTFxWnOuX79etjZ2eHQoUOoXr06LC0tNQmgOoYNGzZgz549mnMHBQUBAGrUqAFXV1fsKmwdsokoQ7a24raw9XDk4G25NW8O3Lwp7jdsKBbAMzOTNyYiAgDExwOWlgXxSgoARlolcXGAhYVuZxk4cCDWrVuH3r17AwDWrl2LAQMGaL585sbJkydRuXJlWFlZacqePHkCb29vdOjQAePHj4flfx+WlZUV1q9fD1dXV1y9ehWDBw+GlZUVJqrX4XlDSkoKfHx8MHjwYGzduhVJSUk4f/48FAqF5nX79u2LpUuXonnz5rh79y6GDBkCAFl2/Tl69Cisra1x+PBhAEBycjLat28PLy8vnDx5EoaGhpg9ezY6dOiAK1euQKlUZhlHdmP28/PDtWvXcODAARw5cgQAYGNjAwBQKpVYunQpypUrh3v37mH48OGYOHEiVqxYoTl3fHw8vvnmG2zcuBEGBgb45JNPMH78eGzevBnjx49HSEgIYmJisG7dOgBAiRIlNMc2atQIJ0+eTNeqRESFj/p/z+nTCsTEKArNStxMLOR082ZaUtG3r1irQslGJCLKmU8++QRTpkzRXGU/ffo0AgMDM0wsfv/9d00ioDZ16lRMfbNz7xsePHgAV1dXrTJnZ2cYGhrC0tISzs7OmvLp06dr7ru7u2P8+PEIDAzMMLGIiYlBdHQ0OnfujAoVKgAAqlWrpnk+ICAAkydPRr9+/QAA5cuXx5dffomJEydmmVhYWFhgzZo1MDY2BgD89NNPUKlUWLNmjSZZWLduHWxtbREUFIQGDRpkGYcuMVtaWsLQ0FDrMwGAMWPGaH0us2fPxrBhw7QSi+TkZCxbtgxVqlSBQqHAyJEjMeu/LrKWlpYwMzNDYmJiunMDgKurKy5dupTpZ0JEhUerVmn3g4ONULGifLHogomFjBTqNds7dgQ2bJA3GCJKx9xctBzkN0mSkJKSAkNDQ82XXnNz3c/j4OCA999/H+vXr4ckSXj//fdhb2+f4b6tWrXCypUrtcrevPr9ttevX8PU1DRbcWzbtg1Lly7F3bt3ERcXh5SUFFir5098S4kSJdC/f3+0b98ebdu2hbe3N3x9feHi4gIAuHz5Mk6fPo05c+ZojklNTUVCQgLi4+NhnskHVatWLU1SoT7PnTt3tFpcADEO5e7du2jXrl2WcegSc2aOHDmCefPm4ebNm4iJiUFKSkq692Fubq5JVgDAxcUFERERWZ5XzczMDPHx8dnal4j0m4MDYGoqViA4e9YYH38sd0TZoxeXx5cvXw53d3eYmprC09MT58+fz3L/HTt2oGrVqjA1NUWtWrWwrxCOR1D+8w8QHCwesNmaSC8pFKI7khxbBj1wsmXgwIFYv349NmzYgIFZTE9tYWGBihUram1ZJRb29vZ48eLFO1//7Nmz6N27Nzp16oTff/8dly5dwrRp07IcQL1u3TqcPXsWTZo0wbZt21C5cmX8+eefAIC4uDgEBAQgODhYs129ehW3b9/OMtGxeKsfWVxcHDw8PLTOExwcjFu3bqFXr17vjEOXmDNy//59dO7cGbVr18bPP/+MixcvYvny5QC0B5cbGWl3iVMoFJCyOeDm+fPncHBwyNa+RKT/WrcWt5s2FZ4u8rInFtu2bcO4cePg7++Pv//+G3Xq1EH79u0zvUJz5swZ9OzZE4MGDcKlS5fg4+MDHx8fXLt2rYAjz4W//4Zjw4ZQXLkiHvMfARHlkQ4dOiApKUkzpiCv1KtXDzdv3nznl9wzZ86gbNmymDZtGho0aIBKlSppDcTO6vxTpkzBmTNnULNmTWz5b47F+vXrIzQ0NF0SVLFiRSh16Dpav3593L59G46OjunOox4DkVUcusRsbGyM1NRUrX0vXrwIlUqFRYsWoXHjxqhcuTIeP36c7fjVMjq32rVr1zIcqE9EhdNHH4nbZ88MdJ7QQy6yJxaLFy/G4MGDMWDAAFSvXh2rVq2Cubk51q5dm+H+3377LTp06IAJEyagWrVq+PLLL1G/fn0sW7asgCPPhQsXAACSmZlYs71RI5kDIqKiwsDAACEhIbhx4wYMslj/JjExEU+fPtXaoqKiMt2/VatWiIuLw/Xr17N8/UqVKuHhw4cIDAzE3bt3sXTp0ixnKgoLC8OUKVNw9uxZPHjwAIcOHcLt27c1YxZmzpyJjRs3IiAgANevX0dISAgCAwO1xnFkR+/evWFvb49u3brh5MmTCAsLQ1BQEEaNGoV//vnnnXHoErO7uzvCwsIQHByMqKgoJCYmomLFikhOTsZ3332He/fuYdOmTVi1apVO70F97itXriA0NBRRUVGa2bbi4+Nx8eJFtGvXTudzEpF++vDDtPv//itfHLqQNbFISkrCxYsX4e3trSlTKpXw9vbG2bNnMzzm7NmzWvsDQPv27TPdX++kpEA5fLi4//HHwOHDohMdEVEesba2znRMg9qBAwfg4uKitTVr1izT/UuWLIkPPvgAmzdvzvK8Xbt2xdixYzFy5EjUrVsXZ86cwYwZMzLd39zcHDdv3sRHH32EypUrY8iQIRgxYgSGDh0KQPx9//3333Ho0CE0bNgQjRs3xjfffIOyZctmGUdGr/PHH3+gTJky+PDDD1GtWjUMGjQICQkJsLa2fmccusT80UcfoUOHDmjVqhUcHBywdetW1KlTB4sXL8ZXX32FmjVrYvPmzZg3b55O7wEABg8ejCpVqqBBgwZwcHDA6dOnAYgpecuUKYPmzZvrfE4i0k/qKWcB4I8/ZAtDJwopu50388Hjx49RqlQpnDlzBl5eXpryiRMn4sSJEzh37ly6Y4yNjbFhwwb07NlTU7ZixQoEBAQgPDw83f6JiYlITEzUPI6JiYGbmxtevHjxzn+8+eLpUyhLlQIApC5bBoV63XYqllQqFSIjI+Hg4KBTtw7KHwkJCbh//z7KlSuX7YHKeSU5OTld/3p9c+XKFbRr1w537txJN6MU5Z2c1AUvLy98/vnnmvEiOZGQkICwsDDNmEeSF/8/EACULKnAy5cKrFuXir595VksLyYmBnZ2doiOjn7nd+ciPyvUvHnzEBAQkK48MjISCQkJBR6P8tkz2NvaIr5mTcR06wZlNmf7oKJJpVIhOjoakiTxH4ceSE5OhkqlQkpKSpYLsOU1SZI0/eYzWjdBX1SvXh1z587F7du3UatWLbnDKZJyUheioqLQrVs3dO/ePVf1NiUlBSqVCs+ePdP7JLc44P8HAoBmzWxw8qQR4uNjERGR+O4D8kFsbGy295U1sbC3t4eBgUG6lobw8PAM5+gGxLzpuuw/ZcoUjBs3TvNY3WLh4OAgT4uFoyNUkZGIjYyEI69CFHsqlQoKhYJXpPREQkICYmNjYWhoCEPDgv/zWBi+zGU10xTlHV3qgrOzMyZPnpzr1zQ0NIRSqUTJkiXZYqEH+P+BAGDXrjdbrmzefUA+0OXvgayJhbGxMTw8PHD06FH4+PgAEL9IR48exciRIzM8xsvLC0ePHtVaaOjw4cNaXaneZGJiAhMTk3TlSqVS1l9UhUIhewykH1gX9IdSqYRCodBsBUWSJM3r6XOLBeU/OeuCut7z75H+4M+DAPnrgS6vK3tXqHHjxqFfv35o0KABGjVqhCVLluDVq1cYMGAAAKBv374oVaqUZpDb6NGj0aJFCyxatAjvv/8+AgMD8ddff+GHH36Q820QERERERVrsicWfn5+iIyMxMyZM/H06VPUrVsXBw4cgJOTEwDg4cOHWplSkyZNsGXLFkyfPh1Tp05FpUqVsHv3btSsWVOut0BERYyMc1oQyYb1nohyS/bEAgBGjhyZadenoKCgdGXdu3dH9+7d8zkqIipu1Os+JCUlwcys8Kx0SpQX4uPjARSOsT5EpJ/0IrEgItIHhoaGMDc3R2RkJIyMjAqsP6skSUhJSYGhoSHHWBRzctQFSZIQHx+PiIgI2NraZrmwIhFRVphYEBH9R6FQwMXFBWFhYXjw4EGBva4kSVCpVJrB41R8yVkXbG1tM51hkYgoO5hYEBG9wdjYGJUqVUJSUlKBvaZ67YCSJUty9pdiTq66YGRkxJYKIso1JhZERG9RKpUFOo+/SqWCkZERTE1NmVgUc6wLRFSY8a8WERERERHlGhMLIiIiIiLKNSYWRERERESUa8VujIV6AaCYmBjZYlCpVIiNjWUfWmJdIACsB5SGdYHUWBcI0I96oP7OnJ1FNItdYhEbGwsAcHNzkzkSIiIiIqLCITY2FjY2Nlnuo5Cyk34UISqVCo8fP4aVlZVs88XHxMTAzc0Njx49grW1tSwxkH5gXSCA9YDSsC6QGusCAfpRDyRJQmxsLFxdXd/ZalLsWiyUSiVKly4tdxgAAGtra/6xIACsCySwHpAa6wKpsS4QIH89eFdLhRo77RERERERUa4xsSAiIiIiolxjYiEDExMT+Pv7w8TERO5QSGasCwSwHlAa1gVSY10goPDVg2I3eJuIiIiIiPIeWyyIiIiIiCjXmFgQEREREVGuMbEgIiIiIqJcY2KRT5YvXw53d3eYmprC09MT58+fz3L/HTt2oGrVqjA1NUWtWrWwb9++AoqU8psudWH16tVo3rw57OzsYGdnB29v73fWHSocdP2boBYYGAiFQgEfH5/8DZAKjK514eXLlxgxYgRcXFxgYmKCypUr839EEaFrXViyZAmqVKkCMzMzuLm5YezYsUhISCigaCk//PHHH+jSpQtcXV2hUCiwe/fudx4TFBSE+vXrw8TEBBUrVsT69evzPc5skyjPBQYGSsbGxtLatWul69evS4MHD5ZsbW2l8PDwDPc/ffq0ZGBgIC1YsEC6ceOGNH36dMnIyEi6evVqAUdOeU3XutCrVy9p+fLl0qVLl6SQkBCpf//+ko2NjfTPP/8UcOSUl3StB2phYWFSqVKlpObNm0vdunUrmGApX+laFxITE6UGDRpInTp1kk6dOiWFhYVJQUFBUnBwcAFHTnlN17qwefNmycTERNq8ebMUFhYmHTx4UHJxcZHGjh1bwJFTXtq3b580bdo06ZdffpEASLt27cpy/3v37knm5ubSuHHjpBs3bkjfffedZGBgIB04cKBgAn4HJhb5oFGjRtKIESM0j1NTUyVXV1dp3rx5Ge7v6+srvf/++1plnp6e0tChQ/M1Tsp/utaFt6WkpEhWVlbShg0b8itEKgA5qQcpKSlSkyZNpDVr1kj9+vVjYlFE6FoXVq5cKZUvX15KSkoqqBCpgOhaF0aMGCG1bt1aq2zcuHFS06ZN8zVOKjjZSSwmTpwo1ahRQ6vMz89Pat++fT5Gln3sCpXHkpKScPHiRXh7e2vKlEolvL29cfbs2QyPOXv2rNb+ANC+fftM96fCISd14W3x8fFITk5GiRIl8itMymc5rQezZs2Co6MjBg0aVBBhUgHISV349ddf4eXlhREjRsDJyQk1a9bE3LlzkZqaWlBhUz7ISV1o0qQJLl68qOkude/ePezbtw+dOnUqkJhJP+j7d0ZDuQMoaqKiopCamgonJyetcicnJ9y8eTPDY54+fZrh/k+fPs23OCn/5aQuvG3SpElwdXVN90eECo+c1INTp07hxx9/RHBwcAFESAUlJ3Xh3r17OHbsGHr37o19+/bhzp07GD58OJKTk+Hv718QYVM+yEld6NWrF6KiotCsWTNIkoSUlBQMGzYMU6dOLYiQSU9k9p0xJiYGr1+/hpmZmUyRCWyxINJT8+fPR2BgIHbt2gVTU1O5w6ECEhsbiz59+mD16tWwt7eXOxySmUqlgqOjI3744Qd4eHjAz88P06ZNw6pVq+QOjQpYUFAQ5s6dixUrVuDvv//GL7/8gr179+LLL7+UOzQiDbZY5DF7e3sYGBggPDxcqzw8PBzOzs4ZHuPs7KzT/lQ45KQuqH399deYP38+jhw5gtq1a+dnmJTPdK0Hd+/exf3799GlSxdNmUqlAgAYGhoiNDQUFSpUyN+gKV/k5G+Ci4sLjIyMYGBgoCmrVq0anj59iqSkJBgbG+drzJQ/clIXZsyYgT59+uDTTz8FANSqVQuvXr3CkCFDMG3aNCiVvFZcHGT2ndHa2lr21gqALRZ5ztjYGB4eHjh69KimTKVS4ejRo/Dy8srwGC8vL639AeDw4cOZ7k+FQ07qAgAsWLAAX375JQ4cOIAGDRoURKiUj3StB1WrVsXVq1cRHBys2bp27YpWrVohODgYbm5uBRk+5aGc/E1o2rQp7ty5o0kuAeDWrVtwcXFhUlGI5aQuxMfHp0se1AmnJEn5FyzpFb3/zij36PGiKDAwUDIxMZHWr18v3bhxQxoyZIhka2srPX36VJIkSerTp480efJkzf6nT5+WDA0Npa+//loKCQmR/P39Od1sEaFrXZg/f75kbGws7dy5U3ry5Ilmi42NlestUB7QtR68jbNCFR261oWHDx9KVlZW0siRI6XQ0FDp999/lxwdHaXZs2fL9RYoj+haF/z9/SUrKytp69at0r1796RDhw5JFSpUkHx9feV6C5QHYmNjpUuXLkmXLl2SAEiLFy+WLl26JD148ECSJEmaPHmy1KdPH83+6ulmJ0yYIIWEhEjLly/ndLPFwXfffSeVKVNGMjY2lho1aiT9+eefmudatGgh9evXT2v/7du3S5UrV5aMjY2lGjVqSHv37i3giCm/6FIXypYtKwFIt/n7+xd84JSndP2b8CYmFkWLrnXhzJkzkqenp2RiYiKVL19emjNnjpSSklLAUVN+0KUuJCcnS1988YVUoUIFydTUVHJzc5OGDx8uvXjxouADpzxz/PjxDP/vq3/2/fr1k1q0aJHumLp160rGxsZS+fLlpXXr1hV43JlRSBLbz4iIiIiIKHc4xoKIiIiIiHKNiQUREREREeUaEwsiIiIiIso1JhZERERERJRrTCyIiIiIiCjXmFgQEREREVGuMbEgIiIiIqJcY2JBRERERES5xsSCiKgYcnd3x5IlS/J83/wkVxwKhQK7d+/O1TlatmyJMWPGZLmPvnzOREQ5xcSCiEhP9O/fHwqFAgqFAkZGRnByckLbtm2xdu1aqFSqPH2tCxcuYMiQIXm+b068+b4z2tzd3fPttYmIKO8wsSAi0iMdOnTAkydPcP/+fezfvx+tWrXC6NGj0blzZ6SkpOTZ6zg4OMDc3DzP982Jb7/9Fk+ePNFsALBu3TrN4wsXLuT43MnJyXkVJhERvQMTCyIiPWJiYgJnZ2eUKlUK9evXx9SpU7Fnzx7s378f69ev1+z38uVLfPrpp3BwcIC1tTVat26Ny5cva53rt99+Q8OGDWFqagp7e3t88MEHmufe7HYjSRK++OILlClTBiYmJnB1dcWoUaMy3BcAHj58iG7dusHS0hLW1tbw9fVFeHi45vkvvvgCdevWxaZNm+Du7g4bGxv06NEDsbGxGb5nGxsbODs7azYAsLW11Tx2cHDQ7BsfH4+BAwfCysoKZcqUwQ8//KB57v79+1AoFNi2bRtatGgBU1NTbN68GQCwZs0aVKtWDaampqhatSpWrFihOS4pKQkjR46Ei4sLTE1NUbZsWcybN08rxqioKHzwwQcwNzdHpUqV8Ouvv2o9f+LECTRq1AgmJiZwcXHB5MmTs0wEIyIi0KVLF5iZmaFcuXKaOImICjMmFkREeq5169aoU6cOfvnlF01Z9+7dERERgf379+PixYuoX78+2rRpg+fPnwMA9u7diw8++ACdOnXCpUuXcPToUTRq1CjD8//888/45ptv8P333+P27dvYvXs3atWqleG+KpUK3bp1w/Pnz3HixAkcPnwY9+7dg5+fn9Z+d+/exe7du/H777/j999/x4kTJzB//vxcfxaLFi1CgwYNcOnSJQwfPhyfffYZQkNDtfaZPHkyRo8ejZCQELRv3x6bN2/GzJkzMWfOHISEhGDu3LmYMWMGNmzYAABYunQpfv31V2zfvh2hoaHYvHlzuu5XAQEB8PX1xZUrV9CpUyf07t1b81n/+++/6NSpExo2bIjLly9j5cqV+PHHHzF79uxM30f//v3x6NEjHD9+HDt37sSKFSsQERGR68+HiEhWEhER6YV+/fpJ3bp1y/A5Pz8/qVq1apIkSdLJkycla2trKSEhQWufChUqSN9//70kSZLk5eUl9e7dO9PXKlu2rPTNN99IkiRJixYtkipXriwlJSW9c99Dhw5JBgYG0sOHDzXPX79+XQIgnT9/XpIkSfL395fMzc2lmJgYzT4TJkyQPD09M3/zbwAg7dq1K8M4PvnkE81jlUolOTo6SitXrpQkSZLCwsIkANKSJUu0jqtQoYK0ZcsWrbIvv/xS8vLykiRJkj7//HOpdevWkkqlyjSe6dOnax7HxcVJAKT9+/dLkiRJU6dOlapUqaJ1/PLlyyVLS0spNTVVkiRJatGihTR69GhJkiQpNDRU6/OSJEkKCQmRAGg+ZyKiwogtFkREhYAkSVAoFACAy5cvIy4uDiVLloSlpaVmCwsLw927dwEAwcHBaNOmTbbO3b17d7x+/Rrly5fH4MGDsWvXrky78YSEhMDNzQ1ubm6asurVq8PW1hYhISGaMnd3d1hZWWkeu7i45MkV+dq1a2vuKxQKODs7pztvgwYNNPdfvXqFu3fvYtCgQVqf1ezZszWfVf/+/REcHIwqVapg1KhROHToUJava2FhAWtra83rhoSEwMvLS/PzAYCmTZsiLi4O//zzT7pzhYSEwNDQEB4eHpqyqlWrwtbWVsdPg4hIvxjKHQAREb1bSEgIypUrBwCIi4uDi4sLgoKC0u2n/nJqZmaW7XO7ubkhNDQUR44cweHDhzF8+HAsXLgQJ06cgJGRUY7iffs4hUKRJzNbZee8FhYWmvtxcXEAgNWrV8PT01NrPwMDAwBA/fr1ERYWhv379+PIkSPw9fWFt7c3du7cme/vh4ioKGGLBRGRnjt27BiuXr2Kjz76CID4Ivz06VMYGhqiYsWKWpu9vT0AcYX96NGj2X4NMzMzdOnSBUuXLkVQUBDOnj2Lq1evptuvWrVqePToER49eqQpu3HjBl6+fInq1avn8p3mPScnJ7i6uuLevXvpPit1ogYA1tbW8PPzw+rVq7Ft2zb8/PPPmjEU71KtWjWcPXsWkiRpyk6fPg0rKyuULl063f5Vq1ZFSkoKLl68qCkLDQ3Fy5cvc/5GiYj0AFssiIj0SGJiIp4+fYrU1FSEh4fjwIEDmDdvHjp37oy+ffsCALy9veHl5QUfHx8sWLAAlStXxuPHjzUDths0aAB/f3+0adMGFSpUQI8ePZDy//bumCW1MI7j+O+Ii0GoCJU5JA45aQmKQ+ABcThzU+AUhYZYhEuLEK4NIS4NBjq4hJMKh4bAXXoDbvUKmgXl3uFC3IteuPc+93Idvp/5OTycs32f/wNnPpfrurq5uVnas9vtarFYKJvNamNjQ71eTz6fT3t7e0trC4WCEomEisWims2m5vO5KpWKbNv+4QrSOmk0Grq6upLf75fjOJrNZnp9fdXHx4dqtZru7+8VDoeVSqXk8XjU7/e1s7Pzy1eTKpWKms2mLi8vVa1WNZ1OdXt7q1qtJo9n+fwuHo/LcRyVy2U9PDzI6/Xq+vr6t6ZMALCOmFgAwBp5fn5WOBxWNBqV4zgaj8dqtVoaDAafV3csy5Lrusrlcjo9PdX+/r5OTk70/v6u7e1tSd/+9Nzv9zUcDnV4eKh8Pq/JZLJyz0AgoHa7raOjIyWTSb28vGg0GikUCi2ttSxLg8FAwWBQuVxOhUJBsVhMT09P/+6jGDo/P9fj46M6nY4SiYRs21a32/2cWGxuburu7k7pdFqZTEZvb29yXXdlFKwSiUTkuq4mk4kODg50cXGhs7Mz1ev1nz7T6XS0u7sr27Z1fHysUqmkra2tv/K+APC/WF++n90CAAAAwB9gYgEAAADAGGEBAAAAwBhhAQAAAMAYYQEAAADAGGEBAAAAwBhhAQAAAMAYYQEAAADAGGEBAAAAwBhhAQAAAMAYYQEAAADAGGEBAAAAwBhhAQAAAMDYV7DbsCA4gW0EAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "\n",
    "pipe = Pipeline(\n",
    "    [\n",
    "        (\"scaler\", StandardScaler()),\n",
    "        (\"clf\", LogisticRegression(max_iter=2000, class_weight=\"balanced\")),\n",
    "    ]\n",
    ")\n",
    "cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "y_scores = cross_val_predict(pipe, X.values, y, cv=cv, method=\"predict_proba\")[:, 1]\n",
    "\n",
    "vme_rates, me_rates, thresholds = vme_me_curve(y, y_scores)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.plot(thresholds, vme_rates, \"r-\", label=\"VME (miss resistant)\")\n",
    "ax.plot(thresholds, me_rates, \"b-\", label=\"ME (false resistant)\")\n",
    "ax.set_xlabel(\"Decision Threshold\")\n",
    "ax.set_ylabel(\"Error Rate\")\n",
    "ax.set_title(\"VME / ME Trade-off Curve - Amikacin (Rome)\")\n",
    "ax.legend()\n",
    "ax.set_ylim(-0.05, 1.05)\n",
    "ax.grid(alpha=0.3)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Splitting Utilities\n",
    "\n",
    "Standard random splitting can cause data leakage or imbalanced species\n",
    "representation. MaldiAMRKit provides:\n",
    "\n",
    "- **`stratified_species_drug_split`**: stratified train/test split preserving\n",
    "  species-drug label distributions\n",
    "- **`case_based_split`**: keeps all samples from the same patient in the same\n",
    "  split to prevent leakage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified Species-Drug Split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:54.547226Z",
     "iopub.status.busy": "2026-05-13T10:10:54.547024Z",
     "iopub.status.idle": "2026-05-13T10:10:54.557297Z",
     "shell.execute_reply": "2026-05-13T10:10:54.556534Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train: 375 samples, Test: 94 samples\n",
      "Train R-S: 217R - 158S\n",
      "Test  R-S: 55R - 39S\n"
     ]
    }
   ],
   "source": [
    "X_train, X_test, y_train, y_test = stratified_species_drug_split(\n",
    "    X, y, species, test_size=0.2, random_state=42\n",
    ")\n",
    "\n",
    "print(f\"Train: {len(X_train)} samples, Test: {len(X_test)} samples\")\n",
    "print(f\"Train R-S: {y_train.sum():.0f}R - {len(y_train) - y_train.sum():.0f}S\")\n",
    "print(f\"Test  R-S: {y_test.sum():.0f}R - {len(y_test) - y_test.sum():.0f}S\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case-Based Split\n",
    "\n",
    "In clinical data, multiple spectra may come from the same patient.\n",
    "`case_based_split` ensures no patient appears in both train and test sets.\n",
    "\n",
    "The Zenodo metadata does not expose patient IDs, so we simulate them\n",
    "here by grouping every three consecutive spectra into one synthetic\n",
    "patient. In your own datasets, pass the real patient ID column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:54.558532Z",
     "iopub.status.busy": "2026-05-13T10:10:54.558414Z",
     "iopub.status.idle": "2026-05-13T10:10:54.569295Z",
     "shell.execute_reply": "2026-05-13T10:10:54.568723Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train: 325 samples (109 cases)\n",
      "Test:  144 samples (48 cases)\n",
      "Case overlap: 0 (should be 0)\n"
     ]
    }
   ],
   "source": [
    "case_ids = np.array([f\"patient_{i // 3}\" for i in range(len(X))])\n",
    "\n",
    "X_train_c, X_test_c, y_train_c, y_test_c = case_based_split(\n",
    "    X, y, case_ids, test_size=0.3, random_state=42\n",
    ")\n",
    "\n",
    "train_pos = X.index.get_indexer(X_train_c.index)\n",
    "test_pos = X.index.get_indexer(X_test_c.index)\n",
    "train_cases = set(case_ids[train_pos])\n",
    "test_cases = set(case_ids[test_pos])\n",
    "print(f\"Train: {len(X_train_c)} samples ({len(train_cases)} cases)\")\n",
    "print(f\"Test:  {len(X_test_c)} samples ({len(test_cases)} cases)\")\n",
    "print(f\"Case overlap: {len(train_cases & test_cases)} (should be 0)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-Validation Splitters\n",
    "\n",
    "MaldiAMRKit provides sklearn-compatible CV splitters for use with\n",
    "`cross_val_score`, `GridSearchCV`, etc.\n",
    "\n",
    "- **`SpeciesDrugStratifiedKFold`**: stratified by species-drug combinations\n",
    "- **`CaseGroupedKFold`**: grouped by patient/case ID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:54.570516Z",
     "iopub.status.busy": "2026-05-13T10:10:54.570371Z",
     "iopub.status.idle": "2026-05-13T10:10:54.574131Z",
     "shell.execute_reply": "2026-05-13T10:10:54.573743Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SpeciesDrugStratifiedKFold folds:\n",
      "  Fold 1: train=375, test=94, test R%=58.5%\n",
      "  Fold 2: train=375, test=94, test R%=58.5%\n",
      "  Fold 3: train=375, test=94, test R%=57.4%\n",
      "  Fold 4: train=375, test=94, test R%=57.4%\n",
      "  Fold 5: train=376, test=93, test R%=58.1%\n"
     ]
    }
   ],
   "source": [
    "cv_stratified = SpeciesDrugStratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "\n",
    "print(\"SpeciesDrugStratifiedKFold folds:\")\n",
    "for i, (train_idx, test_idx) in enumerate(\n",
    "    cv_stratified.split(X.values, y, species=species)\n",
    "):\n",
    "    print(\n",
    "        f\"  Fold {i + 1}: train={len(train_idx)}, test={len(test_idx)}, \"\n",
    "        f\"test R%={y[test_idx].mean():.1%}\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:54.575344Z",
     "iopub.status.busy": "2026-05-13T10:10:54.575241Z",
     "iopub.status.idle": "2026-05-13T10:10:54.593823Z",
     "shell.execute_reply": "2026-05-13T10:10:54.593333Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CaseGroupedKFold folds:\n",
      "  Fold 1: train=375, test=94, case overlap=0\n",
      "  Fold 2: train=373, test=96, case overlap=0\n",
      "  Fold 3: train=376, test=93, case overlap=0\n",
      "  Fold 4: train=376, test=93, case overlap=0\n",
      "  Fold 5: train=376, test=93, case overlap=0\n"
     ]
    }
   ],
   "source": [
    "cv_grouped = CaseGroupedKFold(n_splits=5)\n",
    "\n",
    "print(\"CaseGroupedKFold folds:\")\n",
    "for i, (train_idx, test_idx) in enumerate(\n",
    "    cv_grouped.split(X.values, y, groups=case_ids)\n",
    "):\n",
    "    overlap = set(case_ids[train_idx]) & set(case_ids[test_idx])\n",
    "    print(\n",
    "        f\"  Fold {i + 1}: train={len(train_idx)}, test={len(test_idx)}, \"\n",
    "        f\"case overlap={len(overlap)}\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## sklearn Scorers\n",
    "\n",
    "Use `vme_scorer` and `me_scorer` directly with `cross_val_score` or\n",
    "`GridSearchCV` to optimise models for clinical error rates.\n",
    "\n",
    "The scorers return **negative** values because sklearn maximises scores,\n",
    "and we want to **minimise** VME and ME."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:54.594983Z",
     "iopub.status.busy": "2026-05-13T10:10:54.594849Z",
     "iopub.status.idle": "2026-05-13T10:10:56.386636Z",
     "shell.execute_reply": "2026-05-13T10:10:56.385958Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ROC AUC:  0.714 +/- 0.029\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VME:      -0.312 +/- 0.024  (negative = better)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ME:       -0.402 +/- 0.066  (negative = better)\n"
     ]
    }
   ],
   "source": [
    "cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "\n",
    "auc_scores = cross_val_score(pipe, X, y, cv=cv, scoring=\"roc_auc\")\n",
    "print(f\"ROC AUC:  {auc_scores.mean():.3f} +/- {auc_scores.std():.3f}\")\n",
    "\n",
    "vme_scores = cross_val_score(pipe, X, y, cv=cv, scoring=vme_scorer)\n",
    "print(\n",
    "    f\"VME:      {vme_scores.mean():.3f} +/- {vme_scores.std():.3f}  (negative = better)\"\n",
    ")\n",
    "\n",
    "me_scores = cross_val_score(pipe, X, y, cv=cv, scoring=me_scorer)\n",
    "print(\n",
    "    f\"ME:       {me_scores.mean():.3f} +/- {me_scores.std():.3f}  (negative = better)\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified CV + AMR Scorers\n",
    "\n",
    "Combine the domain-specific splitters with the AMR scorers for a complete\n",
    "clinical evaluation workflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-13T10:10:56.388275Z",
     "iopub.status.busy": "2026-05-13T10:10:56.388153Z",
     "iopub.status.idle": "2026-05-13T10:10:56.990067Z",
     "shell.execute_reply": "2026-05-13T10:10:56.989362Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stratified CV with AMR scorers:\n",
      "  VME: -0.312 +/- 0.024\n",
      "  ME:  -0.402 +/- 0.066\n"
     ]
    }
   ],
   "source": [
    "cv_clinical = SpeciesDrugStratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "\n",
    "vme_scores = cross_val_score(\n",
    "    pipe,\n",
    "    X.values,\n",
    "    y,\n",
    "    cv=cv_clinical.split(X.values, y, species=species),\n",
    "    scoring=vme_scorer,\n",
    ")\n",
    "me_scores = cross_val_score(\n",
    "    pipe,\n",
    "    X.values,\n",
    "    y,\n",
    "    cv=cv_clinical.split(X.values, y, species=species),\n",
    "    scoring=me_scorer,\n",
    ")\n",
    "\n",
    "print(\"Stratified CV with AMR scorers:\")\n",
    "print(f\"  VME: {vme_scores.mean():.3f} +/- {vme_scores.std():.3f}\")\n",
    "print(f\"  ME:  {me_scores.mean():.3f} +/- {me_scores.std():.3f}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "maldiamrkit",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}