notebooks/700_metrics/701b_aupimo_advanced_i.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AUPIMO\n",
    "\n",
    "Advance use cases of the metric AUPIMO (pronounced \"a-u-pee-mo\").\n",
    "\n",
    "> For basic usage, please check the notebook [701a_aupimo.ipynb](./701a_aupimo.ipynb).\n",
    "\n",
    "Includes:\n",
    "- selection of test representative samples for qualitative analysis\n",
    "- visualization of the AUPIMO metric with heatmaps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# What is AUPIMO?\n",
    "\n",
    "The `Area Under the Per-Image Overlap [curve]` (AUPIMO) is a metric of recall (higher is better) designed for visual anomaly detection.\n",
    "\n",
    "Inspired by the [ROC](https://en.wikipedia.org/wiki/Receiver_operating_characteristic) and [PRO](https://link.springer.com/article/10.1007/s11263-020-01400-4) curves, \n",
    "\n",
    "> AUPIMO is the area under a curve of True Positive Rate (TPR or _recall_) as a function of False Positive Rate (FPR) restricted to a fixed range. \n",
    "\n",
    "But:\n",
    "- the TPR (Y-axis) is *per-image* (1 image = 1 curve/score);\n",
    "- the FPR (X-axis) considers the (average of) **normal** images only; \n",
    "- the FPR (X-axis) is in log scale and its range is [1e-5, 1e-4]\\* (harder detection task!).\n",
    "\n",
    "\\* The score (the area under the curve) is normalized to be in [0, 1].\n",
    "\n",
    "AUPIMO can be interpreted as\n",
    "\n",
    "> average segmentation recall in an image given that the model (nearly) does not yield false positives in normal images.\n",
    "\n",
    "References in the last cell.\n",
    "\n",
    "![AUROC vs. AUPRO vs. AUPIMO](./roc_pro_pimo.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install `anomalib` using `pip`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO(jpcbertoldo): replace by `pip install anomalib` when AUPIMO is released  # noqa: TD003\n",
    "%pip install ../.."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change the directory to have access to the datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "# NOTE: Provide the path to the dataset root directory.\n",
    "#   If the datasets is not downloaded, it will be downloaded\n",
    "#   to this directory.\n",
    "dataset_root = Path.cwd().parent.parent / \"datasets\" / \"MVTec\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import torch\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.ticker import PercentFormatter\n",
    "from scipy import stats\n",
    "\n",
    "from anomalib import TaskType\n",
    "from anomalib.data import MVTec\n",
    "from anomalib.data.utils import read_image\n",
    "from anomalib.engine import Engine\n",
    "from anomalib.metrics import AUPIMO\n",
    "from anomalib.models import Padim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option(\"display.float_format\", \"{:.2f}\".format)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basics\n",
    "\n",
    "This part was covered in the notebook [701a_aupimo.ipynb](./701a_aupimo.ipynb), so we'll not discuss it here.\n",
    "\n",
    "It will train a model and evaluate it using AUPIMO.\n",
    "We will use dataset Leather from MVTec AD with `PaDiM` (performance is not the best, but it is fast to train).\n",
    "\n",
    "> See the notebooks below for more details on:\n",
    "> - datamodules: [100_datamodules](https://github.com/openvinotoolkit/anomalib/tree/main/notebooks/100_datamodules);\n",
    "> - models: [200_models](https://github.com/openvinotoolkit/anomalib/tree/main/notebooks/200_models)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train the model\n",
    "task = TaskType.SEGMENTATION\n",
    "datamodule = MVTec(\n",
    "    root=dataset_root,\n",
    "    category=\"leather\",\n",
    "    image_size=256,\n",
    "    train_batch_size=32,\n",
    "    eval_batch_size=32,\n",
    "    num_workers=8,\n",
    "    task=task,\n",
    ")\n",
    "model = Padim(\n",
    "    # only use one layer to speed it up\n",
    "    layers=[\"layer1\"],\n",
    "    n_features=64,\n",
    "    backbone=\"resnet18\",\n",
    "    pre_trained=True,\n",
    ")\n",
    "engine = Engine(\n",
    "    pixel_metrics=\"AUPIMO\",  # others can be added\n",
    "    accelerator=\"auto\",  # \\<\"cpu\", \"gpu\", \"tpu\", \"ipu\", \"hpu\", \"auto\">,\n",
    "    devices=1,\n",
    "    logger=False,\n",
    ")\n",
    "engine.fit(datamodule=datamodule, model=model)\n",
    "# infer\n",
    "predictions = engine.predict(dataloaders=datamodule.test_dataloader(), model=model, return_predictions=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute AUPIMO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Metric `AUPIMO` will save all targets and predictions in buffer. For large datasets this may lead to large memory footprint.\n"
     ]
    }
   ],
   "source": [
    "aupimo = AUPIMO(\n",
    "    # with `False` all the values are returned in a dataclass\n",
    "    return_average=False,\n",
    ")\n",
    "\n",
    "anomaly_maps = []\n",
    "masks = []\n",
    "labels = []\n",
    "image_paths = []\n",
    "for batch in predictions:\n",
    "    anomaly_maps.append(batch_anomaly_maps := batch[\"anomaly_maps\"].squeeze(dim=1))\n",
    "    masks.append(batch_masks := batch[\"mask\"])\n",
    "    labels.append(batch[\"label\"])\n",
    "    image_paths.append(batch[\"image_path\"])\n",
    "    aupimo.update(anomaly_maps=batch_anomaly_maps, masks=batch_masks)\n",
    "\n",
    "# list[list[str]] -> list[str]\n",
    "image_paths = [item for sublist in image_paths for item in sublist]\n",
    "anomaly_maps = torch.cat(anomaly_maps, dim=0)\n",
    "masks = torch.cat(masks, dim=0)\n",
    "labels = torch.cat(labels, dim=0)\n",
    "\n",
    "# `pimo_result` has the PIMO curves of each image\n",
    "# `aupimo_result` has the AUPIMO values\n",
    "#     i.e. their Area Under the Curve (AUC)\n",
    "pimo_result, aupimo_result = aupimo.compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statistics and score distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MEAN\n",
      "aupimo_result.aupimos[labels == 1].mean().item()=0.742841961578308\n",
      "OTHER STATISTICS\n",
      "DescribeResult(nobs=92, minmax=(0.0, 1.0), mean=0.742841961578308, variance=0.08757792704451817, skewness=-0.9285678601866055, kurtosis=-0.3299211772047075)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAHHCAYAAABa2ZeMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHeUlEQVR4nO3deVyVZeL+8eugrCKgiKAJgktuqZVLWo0raoqW6aiZlFtuqVk21tdpUcuyZbLFXBvDmkDLMkctNXMrjVwwTVNJc8FccAVUEBHu3x8ez29OuHAQzgH8vF8vXuO5n+c8z3W4G7t6tmMxxhgBAABAbq4OAAAAUFRQjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYEUxAgAAsKIYAQAAWFGMABQba9askcVi0Zo1a2xj/fr1U3h4uFP2Hx4ern79+tlez5kzRxaLRZs3b3bK/lu1aqVWrVo5ZV/ArYpiBBQB06ZNk8Vi0T333HPV5QcOHJDFYtG//vWvqy7/17/+JYvFogMHDtjGWrVqJYvFYvspX768mjRpoo8//lg5OTm29fr16ydfX1+77V15b82aNa+6vxUrVti2++WXX+Za/ttvvyk6Olq33XabPD09VblyZfXp00e//fbbjX4VTrFz506NHz/e7vdVVBTlbMCtoLSrAwCQYmNjFR4ero0bN2rv3r2qUaNGgWy3SpUqmjRpkiTpxIkT+vTTTzVw4ED9/vvveuONN677Xi8vL+3du1cbN25U06ZNc+X18vLShQsXcr1vwYIF6t27t8qXL6+BAwcqIiJCBw4c0OzZs/Xll19q3rx5evjhhwvk80nSRx99ZFf08mLnzp2aMGGCWrVq5dDRpsTERLm5Fe5/T14v23fffVeo+wbAESPA5fbv36+ffvpJkydPVlBQkGJjYwts2/7+/oqOjlZ0dLSeeeYZrV+/XlWqVNGHH36orKys6763evXqqlWrlubOnWs3fuHCBX399deKiorK9Z4//vhDjz32mKpVq6Zff/1VEydO1MCBA/Xqq6/q119/VbVq1fTYY49p3759BfYZ3d3d5enpWWDb+ytjjDIyMiRJnp6ecnd3L7R93YiHh4c8PDxctn/gVkAxAlwsNjZW5cqVU1RUlP7+978XaDH6Kx8fHzVr1kznz5/XiRMnbrh+79699fnnn9sdkVm8eLHS09PVs2fPXOu//fbbSk9P16xZsxQUFGS3rEKFCpo5c6bOnz+vt95664b7/vPPP9W1a1eVKVNGFStW1DPPPKPMzMxc613tGqN58+apUaNGKlu2rPz8/FS/fn29//77ki5fF9SjRw9JUuvWrW2nBK9ctxQeHq7OnTtr+fLlaty4sby9vTVz5kzbsv+9xuiK9PR0DRkyRIGBgfLz89Pjjz+uM2fO2K1jsVg0fvz4XO/9323eKNvVrjE6fvy4Bg4cqODgYHl5ealhw4b65JNP7Nb531Oxs2bNUvXq1eXp6akmTZpo06ZNuTIBtzJOpQEuFhsbq27dusnDw0O9e/fW9OnTtWnTJjVp0qRQ9rdv3z6VKlVKAQEBN1z30Ucf1fjx47VmzRq1adNGkhQXF6e2bduqYsWKudZfvHixwsPD9be//e2q22vRooXCw8P1zTffXHe/GRkZatu2rZKSkvTUU0+pcuXK+s9//qNVq1bdMPOKFSvUu3dvtW3bVm+++aYkadeuXVq/fr1GjRqlFi1a6KmnntIHH3ygf/7zn6pTp44k2f5XunzKrHfv3hoyZIgGDRqkWrVqXXefI0aMUEBAgMaPH6/ExERNnz5dBw8etF0snld5yfa/MjIy1KpVK+3du1cjRoxQRESE5s+fr379+iklJUWjRo2yWz8uLk5nz57VkCFDZLFY9NZbb6lbt27at2+fS4+EAUUJxQhwoYSEBO3evVtTpkyRJN1///2qUqWKYmNjC6QYZWdn6+TJk5KkkydPavr06dqyZYu6dOkiHx+fG76/Zs2aaty4seLi4tSmTRulpKTo22+/1UcffZRr3dTUVB05ckQPPfTQdbfZoEEDLVq0SGfPnlXZsmWvus6sWbP0+++/64svvrAdQRk0aJAaNmx4w8zffPON/Pz8tHz5cpUqVSrX8mrVqulvf/ubPvjgA7Vr1+6qd3nt3btXy5YtU4cOHW64P+nyKa6VK1faykXVqlX13HPPafHixXrwwQfztI28Zvtfs2bN0q5du/TZZ5+pT58+kqShQ4eqZcuWevHFFzVgwAC733FSUpL27NmjcuXKSZJq1aqlhx56SMuXL1fnzp3znBMoyTiVBrhQbGysgoOD1bp1a0mXT7f06tVL8+bNU3Z29k1vf/fu3QoKClJQUJDq1KmjKVOmKCoqSh9//HGet/Hoo49qwYIFunjxor788kuVKlXqqhdPnz17VpKuWXauuLI8LS3tmut8++23qlSpkv7+97/bxnx8fDR48OAb5g0ICND58+e1YsWKG657LREREXkuRZI0ePBguyMuw4YNU+nSpfXtt9/mO0NefPvttwoJCVHv3r1tY+7u7nrqqad07tw5rV271m79Xr162UqRJNuRvYK85gso7ihGgItkZ2dr3rx5at26tfbv36+9e/dq7969uueee5ScnKyVK1c6vM2/nrYJDw/XihUr9P3332vdunU6duyYlixZogoVKuR5m4888ohSU1O1dOlSxcbGqnPnzlctP1fGrhSka8lLgTp48KBq1KiR6/Pc6JSWJD355JO6/fbb1bFjR1WpUkUDBgzQsmXLbvi+/xUREeHQ+n99rIGvr68qVapU6LfcHzx4UDVr1sx1p9yVU28HDx60Gw8LC7N7faUk/fV6KOBWxqk0wEVWrVqlo0ePat68eZo3b16u5bGxsWrfvr2ky7fOS7LdHfVX6enpdutdUaZMGUVGRt5UzkqVKqlVq1Z65513tH79en311VdXXc/f31+VKlXSr7/+et3t/frrr7rtttvk5+d3U7mupWLFitq6dauWL1+upUuXaunSpYqJidHjjz+e66Lka/H29i6UbFdTEEcG8+pqpxaly3feAbiMI0aAi8TGxqpixYqaP39+rp/evXvr66+/thWhoKAg+fj4KDEx8arbSkxMlI+Pj0NHghzx6KOP6scff5Sfn586dep0zfU6d+6s/fv3a926dVdd/uOPP+rAgQM3vJ6latWq+uOPP3L9C/tan/+vPDw81KVLF02bNk1//PGHhgwZok8//VR79+6VlPvI2s3as2eP3etz587p6NGjdnfLlStXTikpKXbrXbx4UUePHrUbcyRb1apVtWfPnlzPcdq9e7dtOQDHUIwAF8jIyNCCBQvUuXNn/f3vf8/1M2LECJ09e1aLFi2SdPm/9Nu3b6/FixcrKSnJbltJSUlavHix2rdvf80jAjfr73//u8aNG6dp06Zd9zk6Y8aMkbe3t4YMGaJTp07ZLTt9+rSGDh0qHx8fjRkz5rr769Spk44cOWL3VO0rjwG4kb/u183NTQ0aNJAk2+3+ZcqUkaRcRSW/Zs2aZfdcqOnTp+vSpUvq2LGjbax69er64Ycfcr3vr0eMHMnWqVMnHTt2TJ9//rlt7NKlS5oyZYp8fX3VsmXL/Hwc4JbGqTTABa7clXWtO5aaNWtme9hjr169JEmvv/66mjVrprvvvluDBw9WeHi4Dhw4oFmzZslisej1118vtLz+/v5XfQbPX9WsWVOffPKJ+vTpo/r16+d68vXJkyc1d+5cVa9e/brbGTRokD788EM9/vjjSkhIUKVKlfSf//wnT3fSPfHEEzp9+rTatGmjKlWq6ODBg5oyZYruvPNO27U3d955p0qVKqU333xTqamp8vT0VJs2ba76CIK8uHjxotq2bauePXsqMTFR06ZN0/333283v0888YSGDh2q7t27q127dtq2bZuWL1+e6yifI9kGDx6smTNnql+/fkpISFB4eLi+/PJLrV+/Xu+9994NL4QHcBUGgNN16dLFeHl5mfPnz19znX79+hl3d3dz8uRJ29iuXbtMr169TMWKFU3p0qVNxYoVzSOPPGJ27dqV6/0tW7Y09erVu2GWvn37mjJlyjj83tWrVxtJZv78+bmW/frrr6Z3796mUqVKxt3d3YSEhJjevXub7du33zDPFQcPHjQPPvig8fHxMRUqVDCjRo0yy5YtM5LM6tWr7fJXrVrV9vrLL7807du3NxUrVjQeHh4mLCzMDBkyxBw9etRu+x999JGpVq2aKVWqlN02q1ataqKioq6aqWrVqqZv37621zExMUaSWbt2rRk8eLApV66c8fX1NX369DGnTp2ye292drZ5/vnnTYUKFYyPj4/p0KGD2bt3b65tXi9by5YtTcuWLe3WTU5ONv379zcVKlQwHh4epn79+iYmJsZunf379xtJ5u233871mSSZcePGXfXzArciizFcdQcAACBxjREAAIANxQgAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKxK/AMec3JydOTIEZUtW7bAvwYAAAAUDmOMzp49q8qVK+f6ouTCVOKL0ZEjRxQaGurqGAAAIB8OHTqkKlWqOG1/Jb4YXXkk/v79+1W+fHkXp7m1ZWVl6bvvvlP79u3l7u7u6ji3NOaiaGE+ig7moug4ffq0IiIinP7VNiW+GF05fVa2bFn5+fm5OM2tLSsrSz4+PvLz8+MvHBdjLooW5qPoYC6KjitfzOzsy2C4+BoAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYEUxAgAAsKIYAQAAWFGMAAAArChGAAAAVqVdHQAAABSspKQknTx50tUxbkpaWppL9ksxAgCgBElKSlKt2nV0ISPd1VFuipeXl0v2SzECAKAEOXnypC5kpCuw87NyDwx1dZx8s5zar6NL3nf6filGAACUQO6BofIMqeHqGPlmcjJdsl8uvgYAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYEUxAgAAsKIYAQAAWFGMAAAArChGAAAAVhQjAAAAK4oRAACAFcUIAADAimIEAABgRTECAACwohgBAABYUYwAAACsKEYAAABWFCMAAAArihEAAIAVxQgAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYFVkitEbb7whi8Wip59+2jZ24cIFDR8+XIGBgfL19VX37t2VnJzsupAAAKBEKxLFaNOmTZo5c6YaNGhgN/7MM89o8eLFmj9/vtauXasjR46oW7duLkoJAABKOpcXo3PnzqlPnz766KOPVK5cOdt4amqqZs+ercmTJ6tNmzZq1KiRYmJi9NNPP+nnn392YWIAAFBSlXZ1gOHDhysqKkqRkZGaOHGibTwhIUFZWVmKjIy0jdWuXVthYWGKj49Xs2bNrrq9zMxMZWZm2l6npaVJkrKyspSVlVVInwJ5ceX3zzy4HnNRtDAfRUdJmIucnBx5e3vLq7RFHqWMq+PkW46LGopLi9G8efO0ZcsWbdq0KdeyY8eOycPDQwEBAXbjwcHBOnbs2DW3OWnSJE2YMCHX+OrVq+Xj43PTmXHzVqxY4eoIsGIuihbmo+go7nMxd+5c65+yXZrjZqSnh+rRGOfv12XF6NChQxo1apRWrFghLy+vAtvu2LFjNXr0aNvrtLQ0hYaGqnXr1goMDCyw/cBxWVlZWrFihdq1ayd3d3dXx7mlMRdFC/NRdJSEudi2bZtatGih4EffkEdwNVfHybeco4dcsl+XFaOEhAQdP35cd999t20sOztbP/zwgz788EMtX75cFy9eVEpKit1Ro+TkZIWEhFxzu56envL09Mw17u7uXmz/IS9pmIuig7koWpiPoqM4z4Wbm5syMjJ04ZKRyba4Ok6+mUuu2a/LilHbtm21fft2u7H+/furdu3aev755xUaGip3d3etXLlS3bt3lyQlJiYqKSlJzZs3d0VkAABQwrmsGJUtW1Z33HGH3ViZMmUUGBhoGx84cKBGjx6t8uXLy8/PTyNHjlTz5s2veeE1AADAzXD5XWnX8+6778rNzU3du3dXZmamOnTooGnTprk6FgAAKKGKVDFas2aN3WsvLy9NnTpVU6dOdU0gAABwS3H5Ax4BAACKCooRAACAFcUIAADAimIEAABgRTECAACwohgBAABYUYwAAACsKEYAAABWFCMAAAArihEAAIAVxQgAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYEUxAgAAsKIYAQAAWFGMAAAArChGAAAAVhQjAAAAK4oRAACAFcUIAADAimIEAABgRTECAACwohgBAABYUYwAAACsSufnTVlZWTp27JjS09MVFBSk8uXLF3QuAAAAp8vzEaOzZ89q+vTpatmypfz8/BQeHq46deooKChIVatW1aBBg7Rp06bCzAoAAFCo8lSMJk+erPDwcMXExCgyMlILFy7U1q1b9fvvvys+Pl7jxo3TpUuX1L59ez3wwAPas2dPYecGAAAocHk6lbZp0yb98MMPqlev3lWXN23aVAMGDNCMGTMUExOjH3/8UTVr1izQoAAAAIUtT8Vo7ty5edqYp6enhg4delOBAAAAXOWm70pLS0vTwoULtWvXroLIAwAA4DIOF6OePXvqww8/lCRlZGSocePG6tmzpxo0aKCvvvqqwAMCAAA4i8PF6IcfftDf/vY3SdLXX38tY4xSUlL0wQcfaOLEiQUeEAAAwFkcLkapqam25xYtW7ZM3bt3l4+Pj6KiorgbDQAAFGsOF6PQ0FDFx8fr/PnzWrZsmdq3by9JOnPmjLy8vAo8IAAAgLM4/OTrp59+Wn369JGvr6/CwsLUqlUrSZdPsdWvX7+g8wEAADiNw8XoySefVNOmTXXo0CG1a9dObm6XDzpVq1aNa4wAAECxlq/vSmvcuLEaNGig/fv3q3r16ipdurSioqIKOhsAAIBTOXyNUXp6ugYOHCgfHx/Vq1dPSUlJkqSRI0fqjTfeKPCAAAAAzuJwMRo7dqy2bdumNWvW2F1sHRkZqc8//7xAwwEAADiTw6fSFi5cqM8//1zNmjWTxWKxjderV09//PFHgYYDAABwJoePGJ04cUIVK1bMNX7+/Hm7ogQAAFDcOFyMGjdurG+++cb2+koZ+ve//63mzZsXXDIAAAAnc/hU2uuvv66OHTtq586dunTpkt5//33t3LlTP/30k9auXVsYGQEAAJzC4SNG999/v7Zu3apLly6pfv36+u6771SxYkXFx8erUaNGhZERAADAKfL1HKPq1avro48+KugsAAAALuVwMUpLS7vquMVikaenpzw8PG46FAAAgCs4XIwCAgKue/dZlSpV1K9fP40bN872dSEAAADFgcPFaM6cOXrhhRfUr18/NW3aVJK0ceNGffLJJ3rxxRd14sQJ/etf/5Knp6f++c9/FnhgAACAwuJwMfrkk0/0zjvvqGfPnraxLl26qH79+po5c6ZWrlypsLAwvfbaaxQjAABQrDh8ruunn37SXXfdlWv8rrvuUnx8vKTLd65d+Q41AACA4sLhYhQaGqrZs2fnGp89e7ZCQ0MlSadOnVK5cuVuPh0AAIATOXwq7V//+pd69OihpUuXqkmTJpKkzZs3a/fu3fryyy8lSZs2bVKvXr0KNikAAEAhc7gYPfjgg0pMTNTMmTOVmJgoSerYsaMWLlyo8PBwSdKwYcMKNCQAAIAz5OsBj+Hh4Zo0aVJBZwEAAHCpfBUjSUpPT1dSUpIuXrxoN96gQYObDgUAAOAKDhejEydOqH///lq6dOlVl2dnZ990KAAAAFdw+K60p59+WikpKdqwYYO8vb21bNkyffLJJ6pZs6YWLVpUGBkBAACcwuEjRqtWrdJ///tfNW7cWG5ubqpataratWsnPz8/TZo0SVFRUYWREwAAoNA5fMTo/PnzqlixoiSpXLlyOnHihCSpfv362rJlS8GmAwAAcCKHi1GtWrVst+k3bNhQM2fO1OHDhzVjxgxVqlSpwAMCAAA4i8PFaNSoUTp69Kgkady4cVq6dKnCwsL0wQcf6PXXX3doW9OnT1eDBg3k5+cnPz8/NW/e3O6i7gsXLmj48OEKDAyUr6+vunfvruTkZEcjAwAA5InD1xhFR0fb/tyoUSMdPHhQu3fvVlhYmCpUqODQtqpUqaI33nhDNWvWlDFGn3zyiR566CH98ssvqlevnp555hl98803mj9/vvz9/TVixAh169ZN69evdzQ2AADADeX7OUZX+Pj46O67787Xe7t06WL3+rXXXtP06dP1888/q0qVKpo9e7bi4uLUpk0bSVJMTIzq1Kmjn3/+Wc2aNbvZ6AAAAHYcLkbGGH355ZdavXq1jh8/rpycHLvlCxYsyFeQ7OxszZ8/X+fPn1fz5s2VkJCgrKwsRUZG2tapXbu2wsLCFB8ff81ilJmZqczMTNvrtLQ0SVJWVpaysrLylQ0F48rvn3lwPeaiaGE+io6SMBc5OTny9vaWV2mLPEoZV8fJt5ybPnSTPw7v9umnn9bMmTPVunVrBQcHy2Kx3FSA7du3q3nz5rpw4YJ8fX319ddfq27dutq6das8PDwUEBBgt35wcLCOHTt2ze1NmjRJEyZMyDW+evVq+fj43FRWFIwVK1a4OgKsmIuihfkoOor7XMydO9f6p+L70OX09FA9GuP8/TpcjP7zn/9owYIF6tSpU4EEqFWrlrZu3arU1FR9+eWX6tu3r9auXZvv7Y0dO1ajR4+2vU5LS1NoaKhat26twMDAgoiMfMrKytKKFSvUrl07ubu7uzrOLY25KFqYj6KjJMzFtm3b1KJFCwU/+oY8gqu5Ok6+5Rw95JL9OlyM/P39Va1awf2iPTw8VKNGDUmXL+betGmT3n//ffXq1UsXL15USkqK3VGj5ORkhYSEXHN7np6e8vT0zDXu7u5ebP8hL2mYi6KDuShamI+iozjPhZubmzIyMnThkpHJvrmzOq5kLrlmvw7frj9+/HhNmDBBGRkZhZFHOTk5yszMVKNGjeTu7q6VK1faliUmJiopKUnNmzcvlH0DAIBbm8NHjHr27Km5c+eqYsWKCg8Pz9WoHXn69dixY9WxY0eFhYXp7NmziouL05o1a7R8+XL5+/tr4MCBGj16tMqXLy8/Pz+NHDlSzZs35440AABQKBwuRn379lVCQoKio6Nv+uLr48eP6/HHH9fRo0fl7++vBg0aaPny5WrXrp0k6d1335Wbm5u6d++uzMxMdejQQdOmTcv3/gAAAK7H4WL0zTffaPny5br//vtveuezZ8++7nIvLy9NnTpVU6dOvel9AQAA3IjD1xiFhobKz8+vMLIAAAC4lMPF6J133tFzzz2nAwcOFEIcAAAA18nXd6Wlp6erevXq8vHxyXXx9enTpwssHAAAgDM5XIzee++9QogBAADgevm6Kw0AAKAkylMxSktLs11wfeVLWa+FC7MBAEBxladiVK5cOR09elQVK1ZUQEDAVZ9dZIyRxWJRdnbx/cI6AABwa8tTMVq1apXKly8v6fK31AMAAJREeSpGLVu2vOqfAQAAShKHn2MEAABQUlGMAAAArChGAAAAVnkqRosWLVJWVlZhZwEAAHCpPBWjhx9+WCkpKZKkUqVK6fjx44WZCQAAwCXyVIyCgoL0888/S/r/zysCAAAoafJ0u/7QoUP10EMPyWKxyGKxKCQk5Jrr8oBHAABQXOWpGI0fP16PPPKI9u7dqwcffFAxMTEKCAgo5GgAAADOlecvka1du7Zq166tcePGqUePHvLx8SnMXAAAAE6X52J0xbhx4yRJJ06cUGJioiSpVq1aCgoKKthkAAAATubwc4zS09M1YMAAVa5cWS1atFCLFi1UuXJlDRw4UOnp6YWREQAAwCkcLkbPPPOM1q5dq0WLFiklJUUpKSn673//q7Vr1+rZZ58tjIwAAABO4fCptK+++kpffvmlWrVqZRvr1KmTvL291bNnT02fPr0g8wEAADhNvk6lBQcH5xqvWLEip9IAAECx5nAxat68ucaNG6cLFy7YxjIyMjRhwgQ1b968QMMBAAA4k8On0t5//3116NBBVapUUcOGDSVJ27Ztk5eXl5YvX17gAQEAAJzF4WJ0xx13aM+ePYqNjdXu3bslSb1791afPn3k7e1d4AEBAACcxeFiJEk+Pj4aNGhQQWcBAABwKYevMQIAACipKEYAAABWFCMAAAArh4pRdna2fvjhB6WkpBRSHAAAANdxqBiVKlVK7du315kzZworDwAAgMs4fCrtjjvu0L59+wojCwAAgEs5XIwmTpyof/zjH1qyZImOHj2qtLQ0ux8AAIDiyuHnGHXq1EmS9OCDD8pisdjGjTGyWCzKzs4uuHQAAABO5HAxWr16dWHkAAAAcDmHi1HLli0LIwcAAIDL5es5Rj/++KOio6N177336vDhw5Kk//znP1q3bl2BhgMAAHAmh4vRV199pQ4dOsjb21tbtmxRZmamJCk1NVWvv/56gQcEAABwlnzdlTZjxgx99NFHcnd3t43fd9992rJlS4GGAwAAcCaHi1FiYqJatGiRa9zf358nYgMAgGLN4WIUEhKivXv35hpft26dqlWrViChAAAAXMHhYjRo0CCNGjVKGzZskMVi0ZEjRxQbG6t//OMfGjZsWGFkBAAAcAqHb9f/v//7P+Xk5Kht27ZKT09XixYt5OnpqX/84x8aOXJkYWQEAABwCoeLkcVi0QsvvKAxY8Zo7969OnfunOrWrStfX9/CyAcAAOA0DhejKzw8PFS2bFmVLVuWUgQAAEoEh68xunTpkl566SX5+/srPDxc4eHh8vf314svvqisrKzCyAgAAOAUDh8xGjlypBYsWKC33npLzZs3lyTFx8dr/PjxOnXqlKZPn17gIQEAAJzB4WIUFxenefPmqWPHjraxBg0aKDQ0VL1796YYAQCAYsvhU2menp4KDw/PNR4RESEPD4+CyAQAAOASDhejESNG6NVXX7V9R5okZWZm6rXXXtOIESMKNBwAAIAz5elUWrdu3exef//996pSpYoaNmwoSdq2bZsuXryotm3bFnxCAAAAJ8lTMfL397d73b17d7vXoaGhBZcIAADARfJUjGJiYgo7BwAAgMs5fI0RAABASeXw7fqnTp3Syy+/rNWrV+v48ePKycmxW3769OkCCwcAAOBMDhejxx57THv37tXAgQMVHBwsi8VSGLkAAACczuFi9OOPP2rdunW2O9IAAABKCoevMapdu7YyMjIKIwsAAIBLOVyMpk2bphdeeEFr167VqVOnlJaWZvcDAABQXDl8Ki0gIEBpaWlq06aN3bgxRhaLRdnZ2QUWDgAAwJkcLkZ9+vSRu7u74uLiuPgaAACUKA4Xox07duiXX35RrVq1CiMPAACAyzh8jVHjxo116NChwsgCAADgUg4fMRo5cqRGjRqlMWPGqH79+nJ3d7db3qBBgwILBwAA4EwOF6NevXpJkgYMGGAbs1gsXHwNAACKPYeL0f79+wsjBwAAgMs5fI1R1apVr/vjiEmTJqlJkyYqW7asKlasqK5duyoxMdFunQsXLmj48OEKDAyUr6+vunfvruTkZEdjAwAA3JDDR4w+/fTT6y5//PHH87yttWvXavjw4WrSpIkuXbqkf/7zn2rfvr127typMmXKSJKeeeYZffPNN5o/f778/f01YsQIdevWTevXr3c0OgAAwHU5XIxGjRpl9zorK0vp6eny8PCQj4+PQ8Vo2bJldq/nzJmjihUrKiEhQS1atFBqaqpmz56tuLg42wMlY2JiVKdOHf38889q1qyZo/EBAACuyeFidObMmVxje/bs0bBhwzRmzJibCpOamipJKl++vCQpISFBWVlZioyMtK1Tu3ZthYWFKT4+/qrFKDMzU5mZmbbXV76mJCsrS1lZWTeVDzfnyu+feXA95qJoYT6KjpIwFzk5OfL29pZXaYs8ShlXx8m3HIcbSsGwGGMK5Le2efNmRUdHa/fu3fl6f05Ojh588EGlpKRo3bp1kqS4uDj179/fruhIUtOmTdW6dWu9+eabubYzfvx4TZgwIdd4XFycfHx88pUNAAA4V3p6uh599FGlpqbKz8/PafstsD5WunRpHTlyJN/vHz58uHbs2GErRfk1duxYjR492vY6LS1NoaGhat26tQIDA29q27g5WVlZWrFihdq1a5fr+VdwLuaiaGE+io6SMBfbtm1TixYtFPzoG/IIrubqOPmWc9Q1D5N2uBgtWrTI7rUxRkePHtWHH36o++67L18hRowYoSVLluiHH35QlSpVbOMhISG6ePGiUlJSFBAQYBtPTk5WSEjIVbfl6ekpT0/PXOPu7u7F9h/ykoa5KDqYi6KF+Sg6ivNcuLm5KSMjQxcuGZns4vt9puaSa/brcDHq2rWr3WuLxaKgoCC1adNG77zzjkPbMsZo5MiR+vrrr7VmzRpFRETYLW/UqJHc3d21cuVKde/eXZKUmJiopKQkNW/e3NHoAAAA1+VwMcrJySmwnQ8fPlxxcXH673//q7Jly+rYsWOSJH9/f3l7e8vf318DBw7U6NGjVb58efn5+WnkyJFq3rw5d6QBAIAC56Jrvi+bPn26JKlVq1Z24zExMerXr58k6d1335Wbm5u6d++uzMxMdejQQdOmTXNyUgAAcCtwuBhlZ2drzpw5WrlypY4fP57rCNKqVavyvK283BDn5eWlqVOnaurUqY5GBQAAcEi+HvA4Z84cRUVF6Y477pDFUnwv7AIAAPhfDhejefPm6YsvvlCnTp0KIw8AAIDLOPwlsh4eHqpRo0ZhZAEAAHAph4vRs88+q/fffz9P1wcBAAAUJw6fSlu3bp1Wr16tpUuXql69erkegLVgwYICCwcAAOBMDhejgIAAPfzww4WRBQAAwKUcLkYxMTGFkQMAAMDlHL7GCAAAoKTKUzF64IEH9PPPP99wvbNnz+rNN9/kYYwAAKBYytOptB49eqh79+7y9/dXly5d1LhxY1WuXFleXl46c+aMdu7cqXXr1unbb79VVFSU3n777cLODQAAUODyVIwGDhyo6OhozZ8/X59//rlmzZql1NRUSZLFYlHdunXVoUMHbdq0SXXq1CnUwAAAAIUlzxdfe3p6Kjo6WtHR0ZKk1NRUZWRkKDAwMNct+wAAAMWRw3elXeHv7y9/f/+CzAIAAOBS3JUGAABgRTECAACwohgBAABYUYwAAACsHC5G1apV06lTp3KNp6SkqFq1agUSCgAAwBUcLkYHDhxQdnZ2rvHMzEwdPny4QEIBAAC4Qp5v11+0aJHtz8uXL7e7VT87O1srV65UeHh4gYYDAABwpjwXo65du0q6/KTrvn372i1zd3dXeHi43nnnnQINBwAA4Ex5LkY5OTmSpIiICG3atEkVKlQotFAAAACu4PCTr/fv318YOQAAAFwuX18JsnLlSq1cuVLHjx+3HUm64uOPPy6QYAAAAM7mcDGaMGGCXnnlFTVu3FiVKlWSxWIpjFwAAABO53AxmjFjhubMmaPHHnusMPIAAAC4jMPPMbp48aLuvffewsgCAADgUg4XoyeeeEJxcXGFkQUAAMClHD6VduHCBc2aNUvff/+9GjRoIHd3d7vlkydPLrBwAAAAzuRwMfr111915513SpJ27Nhht4wLsQEAQHHmcDFavXp1YeQAAABwOYevMQIAACipHD5i1Lp16+ueMlu1atVNBQIAAHAVh4vRleuLrsjKytLWrVu1Y8eOXF8uCwAAUJw4XIzefffdq46PHz9e586du+lAAAAArlJg1xhFR0fzPWkAAKBYK7BiFB8fLy8vr4LaHAAAgNM5fCqtW7dudq+NMTp69Kg2b96sl156qcCCAQAAOJvDxcjf39/utZubm2rVqqVXXnlF7du3L7BgAAAAzuZwMYqJiSmMHAAAAC7ncDG6IiEhQbt27ZIk1atXT3fddVeBhQIAAHAFh4vR8ePH9cgjj2jNmjUKCAiQJKWkpKh169aaN2+egoKCCjojAACAUzh8V9rIkSN19uxZ/fbbbzp9+rROnz6tHTt2KC0tTU899VRhZAQAAHAKh48YLVu2TN9//73q1KljG6tbt66mTp3KxdcAAKBYc/iIUU5Ojtzd3XONu7u7Kycnp0BCAQAAuILDxahNmzYaNWqUjhw5Yhs7fPiwnnnmGbVt27ZAwwEAADiTw8Xoww8/VFpamsLDw1W9enVVr15dERERSktL05QpUwojIwAAgFM4fI1RaGiotmzZou+//167d++WJNWpU0eRkZEFHg4AAMCZ8vUcI4vFonbt2qldu3YFnQcAAMBl8nwqbdWqVapbt67S0tJyLUtNTVW9evX0448/Fmg4AAAAZ8pzMXrvvfc0aNAg+fn55Vrm7++vIUOGaPLkyQUaDgAAwJnyXIy2bdumBx544JrL27dvr4SEhAIJBQAA4Ap5LkbJyclXfX7RFaVLl9aJEycKJBQAAIAr5LkY3XbbbdqxY8c1l//666+qVKlSgYQCAABwhTwXo06dOumll17ShQsXci3LyMjQuHHj1Llz5wINBwAA4Ex5vl3/xRdf1IIFC3T77bdrxIgRqlWrliRp9+7dmjp1qrKzs/XCCy8UWlAAAIDCludiFBwcrJ9++knDhg3T2LFjZYyRdPmZRh06dNDUqVMVHBxcaEEBAAAKm0MPeKxataq+/fZbnTlzRnv37pUxRjVr1lS5cuUKKx8AAIDT5OvJ1+XKlVOTJk0KOgsAAIBLOfwlsgAAACUVxQgAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuXFqMffvhBXbp0UeXKlWWxWLRw4UK75cYYvfzyy6pUqZK8vb0VGRmpPXv2uCYsAAAo8VxajM6fP6+GDRtq6tSpV13+1ltv6YMPPtCMGTO0YcMGlSlTRh06dNCFCxecnBQAANwK8vWVIAWlY8eO6tix41WXGWP03nvv6cUXX9RDDz0kSfr0008VHByshQsX6pFHHnFmVAAAcAsostcY7d+/X8eOHVNkZKRtzN/fX/fcc4/i4+NdmAwAAJRULj1idD3Hjh2TJAUHB9uNBwcH25ZdTWZmpjIzM22v09LSJElZWVnKysoqhKTIqyu/f+bB9ZiLooX5KDpKwlzk5OTI29tbXqUt8ihlXB0n33Jc1FCKbDHKr0mTJmnChAm5xlevXi0fHx8XJMJfrVixwtURYMVcFC3MR9FR3Odi7ty51j9luzTHzUhPD9WjMc7fb5EtRiEhIZKk5ORkVapUyTaenJysO++885rvGzt2rEaPHm17nZaWptDQULVu3VqBgYGFlhc3lpWVpRUrVqhdu3Zyd3d3dZxbGnNRtDAfRUdJmItt27apRYsWCn70DXkEV3N1nHzLOXrIJfstssUoIiJCISEhWrlypa0IpaWlacOGDRo2bNg13+fp6SlPT89c4+7u7sX2H/KShrkoOpiLooX5KDqK81y4ubkpIyNDFy4ZmWyLq+Pkm7nkmv26tBidO3dOe/futb3ev3+/tm7dqvLlyyssLExPP/20Jk6cqJo1ayoiIkIvvfSSKleurK5du7ouNAAAKLFcWow2b96s1q1b215fOQXWt29fzZkzR88995zOnz+vwYMHKyUlRffff7+WLVsmLy8vV0UGAAAlmEuLUatWrWTMta+Yt1gseuWVV/TKK684MRUAALhVFdlrjAAgL5KSknTy5ElXx7gpOTk5ro4AwIpiBKDYSkpKUq3adXQhI93VUW6Kt7e35s6dqz///FMRERGujgPc0ihGAIqtkydP6kJGugI7Pyv3wFBXx8m3UmlHJEmnTp2iGAEuRjECUOy5B4bKM6SGq2Pkm6V08b2lGihpiux3pQEAADgbxQgAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYEUxAgAAsOJLZAEA+Itt27bJza14HjvYtWuXqyMUaxQjAACs/vzzT0lSixYtlJGR4eI0cAWKEQAAVqdOnZIklX9gpLL9Krs4Tf5k7Nus1B8/c3WMYotiBADAX7iXv02lK1R3dYx8yTp1yNURirXieQIVAACgEFCMAAAArChGAAAAVhQjAAAAKy6+htMV5+eDXFGhQgWFhYW5OsZNK+5zwfNaipakpCSdPHnS1TFuSmJionx9fV0dAy5EMYLTlKTng3h5+yhx965iW45K0lygaEhKSlKt2nV0ISPd1VFuire3t+bOnevqGHAhihGcpiQ8H0S6fCvsqSXv6OTJk8W2GJWUueB5LUXHyZMndSEjXYGdn5V7YKir4+Sb+XOrqyPAxShGcLri/HyQkqa4zwXPayl63AND5RlSw9Ux8u1S2hFXR4CLFd+LCwAAAAoYxQgAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYHXLfFfa9u3b5efn5+oYN6VChQrF9ktLS6Jdu3a5OkK+JSYmytfX19Ux8BeJiYlycyue/71anP//APyvW6YYdezYURcuXHB1jJvi5e2jxN27KEculn3ujGSxKDo62tVR8s3b21tz5851dQxYZZ9PkVRVgwYNUkZGhqvjALe0W6YYlYscIhMY4eoY+ZZ16pBOLXlHJ0+epBi5WE7mOckYBXZ+Vu6Boa6Oky/mz62ujoD/kZN5XpJU/oGRyvar7OI0+ZOxb7NSf/zM1TGAm3bLFCP3cpVlCanh6hgoQdwDQ+VZTP+ZupR2xNURcBXu5W9T6QrVXR0jX7JOHXJ1BKBAFM+T2QAAAIWAYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFYUIwAAACuKEQAAgBXFCAAAwIpiBAAAYEUxAgAAsKIYAQAAWFGMAAAArChGAAAAVhQjAAAAK4oRAACAFcUIAADAimIEAABgRTECAACwohgBAABYlXZ1ADhm165dro6Qb4mJifL19XV1DAAAroliVExknzsjWSyKjo52dZR88/b21ty5c10dAwCAa6IYFRM5meckYxTY+Vm5B4a6Ok6+mD+3ujoCAADXRTEqZtwDQ+UZUsPVMfLlUtoRV0cAAOC6isXF11OnTlV4eLi8vLx0zz33aOPGja6OBAAASqAiX4w+//xzjR49WuPGjdOWLVvUsGFDdejQQcePH3d1NAAAUMIU+WI0efJkDRo0SP3791fdunU1Y8YM+fj46OOPP3Z1NAAAUMIU6WJ08eJFJSQkKDIy0jbm5uamyMhIxcfHuzAZAAAoiYr0xdcnT55Udna2goOD7caDg4O1e/fuq74nMzNTmZmZttepqamSJMuZJJnCi1ro3M4elZeXlyyn9svkZN74DUWQ29ljSk9Pl+X0QeVcvODqOPnGXBQdJWEupJIxH8xF0VFS5sJyJkmSZIyT/+1tirDDhw8bSeann36yGx8zZoxp2rTpVd8zbtw4I4kffvjhhx9++CkBP3/88YczKodNkT5iVKFCBZUqVUrJycl248nJyQoJCbnqe8aOHavRo0fbXqekpKhq1apKSkqSv79/oebF9aWlpSk0NFSHDh2Sn5+fq+Pc0piLooX5KDqYi6IjNTVVYWFhKl++vFP3W6SLkYeHhxo1aqSVK1eqa9eukqScnBytXLlSI0aMuOp7PD095enpmWvc39+ff8iLCD8/P+aiiGAuihbmo+hgLooONzfnXg5dpIuRJI0ePVp9+/ZV48aN1bRpU7333ns6f/68+vfv7+poAACghCnyxahXr146ceKEXn75ZR07dkx33nmnli1bluuCbAAAgJtV5IuRJI0YMeKap85uxNPTU+PGjbvq6TU4F3NRdDAXRQvzUXQwF0WHq+bCYoyz74MDAAAomor0Ax4BAACciWIEAABgRTECAACwohgBAABYlYhiNHXqVIWHh8vLy0v33HOPNm7caFs2evRolS9fXqGhoYqNjbV73/z589WlSxdnxy0RJk2apCZNmqhs2bKqWLGiunbtqsTERLt1Lly4oOHDhyswMFC+vr7q3r273VPMT58+rS5dusjX11d33XWXfvnlF7v3Dx8+XO+8845TPk9J8sYbb8hisejpp5+2jTEXznP48GFFR0crMDBQ3t7eql+/vjZv3mxbbozRyy+/rEqVKsnb21uRkZHas2ePbXlmZqYee+wx+fn56fbbb9f3339vt/23335bI0eOdNrnKa6ys7P10ksvKSIiQt7e3qpevbpeffVVu+/dYi4Kzw8//KAuXbqocuXKslgsWrhwod3yG/3upct/L/Xp00d+fn4KCAjQwIEDde7cOdvyAwcOqEWLFipTpoxatGihAwcO2L2/c+fO+uqrrxwP79QvICkE8+bNMx4eHubjjz82v/32mxk0aJAJCAgwycnJZtGiRSY4ONhs2rTJxMXFGS8vL3PixAljjDEpKSmmZs2a5uDBgy7+BMVThw4dTExMjNmxY4fZunWr6dSpkwkLCzPnzp2zrTN06FATGhpqVq5caTZv3myaNWtm7r33Xtvy0aNHm5YtW5rExETz9NNPm0aNGtmWxcfHm0aNGplLly459XMVdxs3bjTh4eGmQYMGZtSoUbZx5sI5Tp8+bapWrWr69etnNmzYYPbt22eWL19u9u7da1vnjTfeMP7+/mbhwoVm27Zt5sEHHzQREREmIyPDGGPMBx98YOrUqWN27Nhh3n77bRMUFGRycnKMMcbs27fP1KxZ06Smprrk8xUnr732mgkMDDRLliwx+/fvN/Pnzze+vr7m/ffft63DXBSeb7/91rzwwgtmwYIFRpL5+uuv7Zbf6HdvjDEPPPCAadiwofn555/Njz/+aGrUqGF69+5tW96tWzfzyCOPmN9//9307NnTdO/e3bZs3rx5pkuXLvnKXuyLUdOmTc3w4cNtr7Ozs03lypXNpEmTzJtvvml69eplW1axYkWzceNGY4wxgwcPNpMnT3Z63pLq+PHjRpJZu3atMeZy8XR3dzfz58+3rbNr1y4jycTHxxtjjOnYsaOZPn26McaYnTt3Gh8fH2OMMRcvXjQNGzY0mzZtcvKnKN7Onj1ratasaVasWGFatmxpK0bMhfM8//zz5v7777/m8pycHBMSEmLefvtt21hKSorx9PQ0c+fONcYYM2zYMPP8888bY4xJT083kszx48eNMZf/g2TBggWF+AlKjqioKDNgwAC7sW7dupk+ffoYY5gLZ/prMcrL737nzp1Gkt3fPUuXLjUWi8UcPnzYGGNMnTp1zNKlS40xl4tY3bp1jTHGnDlzxtSoUcMkJSXlK2+xPpV28eJFJSQkKDIy0jbm5uamyMhIxcfHq2HDhtq8ebPOnDmjhIQEZWRkqEaNGlq3bp22bNmip556yoXpS5bU1FRJsn3ZX0JCgrKysuzmpnbt2goLC1N8fLwkqWHDhlq1apUuXbqk5cuXq0GDBpKkt956S61atVLjxo2d/CmKt+HDhysqKsrudy4xF860aNEiNW7cWD169FDFihV111136aOPPrIt379/v44dO2Y3F/7+/rrnnnvs5mLdunXKyMjQ8uXLValSJVWoUEGxsbHy8vLSww8/7PTPVRzde++9WrlypX7//XdJ0rZt27Ru3Tp17NhREnPhSnn53cfHxysgIMDu757IyEi5ublpw4YNki7Pz/fff6+cnBx99913tr+3xowZo+HDhys0NDR/AfNVp4qIw4cPG0nmp59+shsfM2aMadq0qTHGmHHjxpnq1aubO+64wyxYsMBkZmaaO+64w2zevNlMmTLF3H777ebee+81O3bscMVHKBGys7NNVFSUue+++2xjsbGxxsPDI9e6TZo0Mc8995wx5vJ/IfTu3duEhYWZFi1amN9++838/vvvpmbNmubkyZNmyJAhJiIiwvTo0cOkpKQ47fMUR3PnzjV33HGH7TD0/x4xYi6cx9PT03h6epqxY8eaLVu2mJkzZxovLy8zZ84cY4wx69evN5LMkSNH7N7Xo0cP07NnT2PM5aN0Tz75pAkPDzeNGzc2P/74ozl16pSpVq2aSUpKMi+88IKpXr26ad++vfnzzz+d/hmLi+zsbPP8888bi8ViSpcubSwWi3n99ddty5kL59Ffjhjl5Xf/2muvmdtvvz3XtoKCgsy0adOMMcb8+eefJioqyoSGhpqoqCjz559/mrVr15rGjRubU6dOmR49epiIiAgzZMgQk5mZmee8xeIrQW7G+PHjNX78eNvrCRMmKDIyUu7u7po4caK2b9+uJUuW6PHHH1dCQoLrghZjw4cP144dO7Ru3TqH3ufv76+4uDi7sTZt2ujtt99WbGys9u3bp8TERA0aNEivvPIKF/9ew6FDhzRq1CitWLFCXl5e+doGc1EwcnJy1LhxY73++uuSpLvuuks7duzQjBkz1Ldv3zxtw93dXVOnTrUb69+/v5566in98ssvWrhwobZt26a33npLTz31VP4uLr0FfPHFF4qNjVVcXJzq1aunrVu36umnn1blypWZixLitttu05IlS2yvMzMz1aFDB33yySeaOHGiypYtq8TERD3wwAOaOXNmni+UL9an0ipUqKBSpUrZ3V0jScnJyQoJCcm1/u7du/XZZ5/p1Vdf1Zo1a9SiRQsFBQWpZ8+e2rJli86ePeus6CXGiBEjtGTJEq1evVpVqlSxjYeEhOjixYtKSUmxW/9acyNJMTExCggI0EMPPaQ1a9aoa9eucnd3V48ePbRmzZpC/BTFW0JCgo4fP667775bpUuXVunSpbV27Vp98MEHKl26tIKDg5kLJ6lUqZLq1q1rN1anTh0lJSVJku33nde/syRp9erV+u233zRixAitWbNGnTp1UpkyZdSzZ0/m4jrGjBmj//u//9Mjjzyi+vXr67HHHtMzzzyjSZMmSWIuXCkvv/uQkBAdP37cbvmlS5d0+vTpa87P66+/rvbt26tRo0Zas2aNunfvLnd3d3Xr1s2h+SnWxcjDw0ONGjXSypUrbWM5OTlauXKlmjdvbreuMUZDhgzR5MmT5evrq+zsbGVlZUmS7X+zs7OdF76YM8ZoxIgR+vrrr7Vq1SpFRETYLW/UqJHc3d3t5iYxMVFJSUm55kaSTpw4oVdeeUVTpkyRpFzzw9xcW9u2bbV9+3Zt3brV9tO4cWP16dPH9mfmwjnuu+++XI+t+P3331W1alVJUkREhEJCQuzmIi0tTRs2bLjqXFx5zMLMmTNVqlQp5sIB6enpcnOz/1dcqVKllJOTI4m5cKW8/O6bN2+ulJQUuzM5q1atUk5Oju65555c29y1a5fi4uL06quvSrrJv7ccPllYxMybN894enqaOXPmmJ07d5rBgwebgIAAc+zYMbv1Zs2aZXcr34YNG4yfn5+Jj483L7/8su1qduTNsGHDjL+/v1mzZo05evSo7Sc9Pd22ztChQ01YWJhZtWqV2bx5s2nevLlp3rz5Vbf36KOPmilTpthev/nmm6ZRo0Zm586dpmPHjubJJ58s9M9UkvzvNUbGMBfOsnHjRlO6dGnz2muvmT179pjY2Fjj4+NjPvvsM9s6b7zxhgkICDD//e9/za+//moeeuihXLcpX/HPf/7TPPvss7bXn3/+uQkLCzPbtm0zAwcONJ06dXLK5yqO+vbta2677Tbb7foLFiwwFSpUsF1XZwxzUZjOnj1rfvnlF/PLL78YSWby5Mnml19+sT0iJy+/+wceeMDcddddZsOGDWbdunWmZs2adrfrX5GTk2Puv/9+s3jxYtvYsGHDTFRUlNm5c6e56667zFtvvZXn7MW+GBljzJQpU0xYWJjx8PAwTZs2NT///LPd8mPHjpmqVavabvG7YsKECaZ8+fKmdu3aZsOGDc6MXOxJuupPTEyMbZ2MjAzz5JNPmnLlyhkfHx/z8MMPm6NHj+ba1rJly0zTpk1Ndna2bez8+fOmR48epmzZsqZt27YmOTnZGR+rxPhrMWIunGfx4sXmjjvuMJ6enqZ27dpm1qxZdstzcnLMSy+9ZIKDg42np6dp27atSUxMzLWd7du3mxo1atg9Gyw7O9sMGzbM+Pn5mSZNmpg9e/YU+ucprtLS0syoUaNMWFiY8fLyMtWqVTMvvPCC3UW4zEXhWb169VX/HdG3b19jTN5+96dOnTK9e/c2vr6+xs/Pz/Tv39+cPXs2175mzJhhd+DDGGOSk5NN27ZtTdmyZU2PHj3M+fPn85zdYsz/PAYUAADgFlasrzECAAAoSBQjAAAAK4oRAACAFcUIAADAimIEAABgRTECAACwohgBAABYUYwAQJe/cNpischisei99967qW21atXKtq2tW7cWSD4AzkExAnBD8fHxKlWqlKKionItW7NmjSwWS64vqZWk8PBwu5JxpSxYLBb5+/vrvvvu06pVq2zL+/Xrp65du9q9tlgsGjp0aK5tDx8+XBaLRf369bMbP3TokAYMGKDKlSvLw8NDVatW1ahRo3Tq1Kkbfs569erp6NGjGjx4sG1s9OjRKl++vEJDQxUbG2u3/vz589WlS5dc21mwYIE2btx4w/0BKHooRgBuaPbs2Ro5cqR++OEHHTly5Ka2FRMTo6NHj2r9+vWqUKGCOnfurH379l1z/dDQUM2bN08ZGRm2sQsXLiguLk5hYWF26+7bt0+NGzfWnj17NHfuXO3du1czZsywfbH06dOnr5utdOnSCgkJkY+PjyRp8eLFiouL03fffae33npLTzzxhE6ePClJSk1N1QsvvKCpU6fm2k758uUVFBSU598JgKKDYgTgus6dO6fPP/9cw4YNU1RUlObMmXNT2wsICFBISIjuuOMOTZ8+XRkZGVqxYsU117/77rsVGhqqBQsW2MYWLFigsLAw3XXXXXbrDh8+XB4eHvruu+/UsmVLhYWFqWPHjvr+++91+PBhvfDCCw5l3bVrl1q1aqXGjRurd+/e8vPz0/79+yVJzz33nIYNG5arnAEo3ihGAK7riy++UO3atVWrVi1FR0fr448/VkF9xaK3t7ck6eLFi9ddb8CAAYqJibG9/vjjj9W/f3+7dU6fPq3ly5frySeftG33ipCQEPXp00eff/65Q9kbNmyozZs368yZM0pISFBGRoZq1KihdevWacuWLXrqqafyvC0AxQPFCMB1zZ49W9HR0ZKkBx54QKmpqVq7du1Nbzc9PV0vvviiSpUqpZYtW1533ejoaK1bt04HDx7UwYMHtX79elumK/bs2SNjjOrUqXPVbdSpU0dnzpzRiRMn8pyxQ4cOio6OVpMmTdSvXz998sknKlOmjIYNG6YZM2Zo+vTpqlWrlu677z799ttved4ugKKrtKsDACi6EhMTtXHjRn399deSLl+D06tXL82ePVutWrXK1zZ79+6tUqVKKSMjQ0FBQZo9e7YaNGhw3fcEBQXZTuMZYxQVFaUKFSpcdd2COpp1xfjx4zV+/Hjb6wkTJigyMlLu7u6aOHGitm/friVLlujxxx9XQkJCge4bgPNRjABc0+zZs3Xp0iVVrlzZNmaMkaenpz788EP5+/vLz89P0uWLkQMCAuzen5KSIn9/f7uxd999V5GRkfL393foAuUBAwZoxIgRknTVC55r1Kghi8WiXbt26eGHH861fNeuXSpXrtxNXRS9e/duffbZZ/rll1/08ccfq0WLFgoKClLPnj01YMAAnT17VmXLls339gG4HqfSAFzVpUuX9Omnn+qdd97R1q1bbT/btm1T5cqVNXfuXElSzZo15ebmlutoyb59+5Samqrbb7/dbjwkJEQ1atRwuKA88MADunjxorKystShQ4dcywMDA9WuXTtNmzbN7g42STp27JhiY2PVq1cvWSwWh/Z7hTFGQ4YM0eTJk+Xr66vs7GxlZWVJku1/s7Oz87VtAEUHR4wAXNWSJUt05swZDRw4MNdRn+7du2v27NkaOnSoypYtqyeeeELPPvusSpcurfr16+vQoUN6/vnn1axZM917770FkqdUqVLatWuX7c9X8+GHH+ree+9Vhw4dNHHiREVEROi3337TmDFjdNttt+m1117L9/7//e9/KygoyPbcovvuu0/jx4/Xzz//rKVLl6pu3bq5jpgBKH44YgTgqmbPnm075fVX3bt31+bNm/Xrr79Kkt5//3317dtXzz//vOrVq6d+/fqpQYMGWrx4cb6P0FyNn5+f7dTd1dSsWVObN29WtWrV1LNnT1WvXl2DBw9W69atFR8fr/Lly+drv8nJyXrttdf0wQcf2MaaNm2qZ599VlFRUfriiy/s7poDUHxZTEFfqQgAxdD48eO1cOHCAvsKjwMHDigiIkK//PKL7rzzzgLZJoDCxxEjALDavn27fH19NW3atJvaTseOHVWvXr0CSgXAmThiBAC6/IDIK18ZEhQUdNVTiHl1+PBh2wXgYWFh8vDwKJCMAAofxQgAAMCKU2kAAABWFCMAAAArihEAAIAVxQgAAMCKYgQAAGBFMQIAALCiGAEAAFhRjAAAAKwoRgAAAFb/D60PVy7Nlq5qAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the normal images have `nan` values because\n",
    "# recall is not defined for them so we ignore them\n",
    "print(f\"MEAN\\n{aupimo_result.aupimos[labels == 1].mean().item()=}\")\n",
    "print(f\"OTHER STATISTICS\\n{stats.describe(aupimo_result.aupimos[labels == 1])}\")\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(aupimo_result.aupimos[labels == 1].numpy(), bins=np.linspace(0, 1, 11), edgecolor=\"black\")\n",
    "ax.set_ylabel(\"Count (number of images)\")\n",
    "ax.set_xlim(0, 1)\n",
    "ax.set_xlabel(\"AUPIMO [%]\")\n",
    "ax.xaxis.set_major_formatter(PercentFormatter(1))\n",
    "ax.grid()\n",
    "ax.set_title(\"AUPIMO distribution\")\n",
    "fig  # noqa: B018, RUF100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Until here we just reproduded the notebook with the basic usage of AUPIMO."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Selecting Representative Samples for Qualitative Analysis\n",
    "\n",
    "Instead of cherry picking or inspecting the 92 samples from above, we'll try to choose them smartly.\n",
    "\n",
    "Our goal here is to select a handful of samples in a meaningful way.\n",
    "\n",
    "> Notice that a random selection from the distribution above would probably miss the worst cases.\n",
    "\n",
    "We will summarize this distribution with a boxplot, then select the samples corresponding to the statistics in it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAADvCAYAAAD2Og4yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0n0lEQVR4nO3dd1hUx/4/8PfSO4giRRFQQbBgAXtNJKIiauJFRRFRExtGjcZrTG4isZdorhpLihpvxBIjJmo0iAV7ReyAqGBFrBQFC/D5/uGP88sGVFAEPb5fz7NP2Jk5M3N2AnnnnNldjYgIiIiIiFRMp6wnQERERPSqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BARvYaio6Oh0WgQHR1d1lMhUgUGHnrrLFiwABqNBo0bNy60Pjk5GRqNBt98802h9d988w00Gg2Sk5OVsjZt2kCj0SgPa2trNGzYEEuWLEFeXp7SLiQkBGZmZlr95R/r6upa6HhRUVFKv7/99luB+tOnTyMoKAiVKlWCoaEhHBwc0Lt3b5w+ffp5L4Xi5s2bGDFiBNzd3WFsbIyKFSuiUaNGGDt2LO7du1fkft4UISEhWuulp6cHR0dH9OzZE2fOnCnr6b20TZs2ISwsrKynQfRa0SvrCRCVtvDwcDg7O+PQoUM4d+4cqlevXiL9Vq5cGVOnTgXwJED873//w4ABA3D27FlMmzbtmccaGRnh3LlzOHToEBo1alRgvkZGRnjw4EGB4yIiIhAYGAhra2sMGDAALi4uSE5OxuLFi/Hbb79h1apVeP/995859p07d+Dt7Y2MjAz0798f7u7uuH37Nk6cOIGFCxdiyJAhBUKaGhgaGuKnn34CAOTk5OD8+fNYtGgR/vrrL5w5cwYODg5lPMMXt2nTJsyfP5+hh+hvGHjorZKUlIR9+/YhIiICgwYNQnh4OMaPH18ifVtaWiIoKEh5PmjQINSoUQPfffcdJk6cCH19/aceW61aNeTk5GDlypVagefBgwdYt24d/Pz8sHbtWq1jzp8/jz59+qBq1arYtWsXbGxslLoRI0agZcuW6NOnD06cOIGqVas+dezFixfj0qVL2Lt3L5o1a6ZVl5GRAQMDgyK/Bi/r/v37MDU1LZWx9PT0tNYLAJo0aYJOnTrhzz//xEcffVQq8yCi0sFbWvRWCQ8PR7ly5eDn54d//etfCA8Pf2VjmZiYoEmTJrh//z5u3rz53PaBgYFYvXq11i2wDRs2ICsrC927dy/QfubMmcjKysIPP/ygFXYAoEKFCvj+++9x//59zJgx45njnj9/Hrq6umjSpEmBOgsLCxgZGWmVHTx4EB07dkS5cuVgamoKT09PzJkzR6vN9u3b0bJlS5iamsLKygpdunRBXFycVpuwsDBoNBqcOXMGvXr1Qrly5dCiRQulfvny5fDy8oKxsTGsra3Rs2dPXL58WauPxMREdOvWDXZ2djAyMkLlypXRs2dPpKenP/Ocn8bOzg7AkzD0dxcuXEBAQACsra2Vdf3zzz+V+ri4OBgbGyM4OFjruD179kBXVxdjx45VypydndGpUyds2bIF9erVg5GREWrWrImIiIgizXHNmjXK61KhQgUEBQXh6tWrSn1ISAjmz58PAFq37Yjedgw89FYJDw/HBx98AAMDAwQGBiIxMRGHDx9+ZeNduHABurq6sLKyem7bXr16ISUlRWuT6ooVK9C2bVtUrFixQPsNGzbA2dkZLVu2LLS/Vq1awdnZWes/zIVxcnJCbm4ufvnll+fOMSoqCq1atcKZM2cwYsQIzJo1C++88w42btyotNm6dSt8fX1x48YNhIWFYdSoUdi3bx+aN2+ute8pX0BAALKysjBlyhTlqsrkyZMRHBwMV1dXzJ49GyNHjsS2bdvQqlUrpKWlAQAePXoEX19fHDhwAB9//DHmz5+PgQMH4sKFC0qb57l16xZu3bqF1NRU7N+/H5988gnKly+PTp06KW1SU1PRrFkzREZGYujQoZg8eTIePHiAzp07Y926dQAADw8PTJw4Eb/88gvWr18P4MnVqpCQELi7u2PChAla4yYmJqJHjx7o0KEDpk6dCj09PQQEBCAqKuqZ8/3555/RvXt36OrqYurUqfjoo48QERGBFi1aKOc8aNAgvPfeewCAX375RXkQvfWE6C1x5MgRASBRUVEiIpKXlyeVK1eWESNGaLVLSkoSADJz5sxC+5k5c6YAkKSkJKWsdevW4u7uLjdv3pSbN29KXFycDB8+XACIv7+/0q5v375iamqq1V/r1q2lVq1aIiLi7e0tAwYMEBGRu3fvioGBgSxbtkx27NghAGTNmjUiIpKWliYApEuXLs88586dOwsAycjIeGqb69evi42NjQAQd3d3GTx4sKxYsULS0tK02uXk5IiLi4s4OTnJ3bt3tery8vKUn+vVqycVK1aU27dvK2XHjx8XHR0dCQ4OVsrGjx8vACQwMFCrr+TkZNHV1ZXJkydrlZ88eVL09PSU8tjYWK3XpDj69u0rAAo8KlWqJDExMVptR44cKQBk9+7dSllmZqa4uLiIs7Oz5ObmiohIbm6utGjRQmxtbeXWrVsSGhoqenp6cvjwYa3+nJycBICsXbtWKUtPTxd7e3upX7++Upa/5jt27BARkUePHknFihWldu3akp2drbTbuHGjAJCvvvpKKQsNDRX+eSfSxis89NYIDw+Hra0t3nnnHQBPLvf36NEDq1atQm5u7kv3Hx8fDxsbG9jY2MDDwwPz5s2Dn58flixZUuQ+evXqhYiICDx69Ai//fYbdHV1C910nJmZCQAwNzd/Zn/59RkZGU9tY2tri+PHj2Pw4MG4e/cuFi1ahF69eqFixYqYOHEiRAQAEBsbi6SkJIwcObLAFav8WyYpKSk4duwYQkJCYG1trdR7enrivffew6ZNmwqMP3jwYK3nERERyMvLQ/fu3ZUrMLdu3YKdnR1cXV2xY8cOAE/2TAFAZGQksrKynvk6FMbIyAhRUVGIiopCZGQkvv/+e5iZmaFjx444e/as0m7Tpk1o1KiR1u02MzMzDBw4EMnJycq7unR0dPDzzz/j3r176NChAxYsWIBx48bB29u7wNgODg5a62phYYHg4GDExsbi+vXrhc73yJEjuHHjBoYOHap1m9HPzw/u7u7PvZJH9LZj4KG3Qm5uLlatWoV33nkHSUlJOHfuHM6dO4fGjRsjNTUV27ZtK3af/9wX4ezsjKioKGzduhV79uzB9evXsXHjRlSoUKHIfebvP9m8eTPCw8PRqVOnQkNNfll+8HmaogYje3t7LFy4ECkpKUhISMDcuXNhY2ODr776CosXLwbwZK8PANSuXfup/Vy8eBEAUKNGjQJ1Hh4euHXrFu7fv69V7uLiovU8MTERIgJXV1clQOY/4uLicOPGDeW4UaNG4aeffkKFChXg6+uL+fPnF3n/jq6uLnx8fODj44N27dph4MCB2Lp1K9LT0zFu3Ditc3ra+fz9nIEnm8/DwsJw+PBh1KpVC19++WWhY1evXr3Avz9ubm4AUOhtv7+PU9hc3N3dteZBRAXxXVr0Vti+fTtSUlKwatUqrFq1qkB9eHg42rVrBwDK/z1nZ2cX2lf+1YR/buY1NTWFj4/PS83T3t4ebdq0waxZs7B3794C78zKZ2lpCXt7e5w4ceKZ/Z04cQKVKlWChYVFkcbXaDRwc3ODm5sb/Pz84OrqivDwcHz44YfFPpeiMjY21nqel5cHjUaDzZs3Q1dXt0D7v79FftasWQgJCcEff/yBLVu2YPjw4Zg6dSoOHDiAypUrF3sulStXRo0aNbBr167in8j/s2XLFgDAtWvXcPv2bWUjNBGVLV7hobdCeHg4KlasiDVr1hR4BAYGYt26dUrAsbGxgYmJCRISEgrtKyEhASYmJsW6clMcvXr1wu7du2FhYYGOHTs+tV2nTp2QlJSEPXv2FFq/e/duJCcna23ALY6qVauiXLlySElJAfDk6gUAnDp16qnHODk5AUChr118fDwqVKjw3LedV6tWDSICFxcX5QrM3x//fDdZnTp18J///Ae7du3C7t27cfXqVSxatKhY5/p3OTk5Wh+26OTk9NTzya/Pt2jRIkRFRWHy5Ml49OgRBg0aVOgY586dU24V5su/jebs7FzoMc96bRMSErTmwXdlERWibLcQEb16WVlZYm5uLv379y+0fu/evQJAVq1apZR17dpVLCws5OLFi1ptL168KObm5tK1a1et8r9vPH6W521aFnmyIXn8+PGyYsUKpeyfm5ZFRM6ePSvGxsZSs2ZNuXXrllaft2/flpo1a4qJiYmcO3fumXM6cOCA3Lt3r0D5wYMHBYB07txZRJ5syi3qpmVbW1utNidPnnzqpuWbN29q9XXu3DnR1dWVXr16afWbP07+uaanp8vjx4+16jMyMkRHR0c+/fTTZ55zYesgIpKQkCB6enrSuHFjpSx/0/K+ffuUsnv37knVqlW1Ni1fuHBBzMzMpFu3biIismjRIgEgy5Yt0xrjWZuW69Wrp5Q9bdOyp6enPHjwQGm3adOmApuWx44dKwAKrBPR24y3tEj11q9fj8zMTHTu3LnQ+iZNmsDGxgbh4eHo0aMHAGDKlClo0qQJGjRogIEDB8LZ2RnJycn44YcfoNFoMGXKlFc2X0tLyyJ9Qq6rqyuWLVuG3r17o06dOgU+afnWrVtYuXKlcmXmaX755ReEh4fj/fffh5eXFwwMDBAXF4clS5bAyMgIn3/+OYAnm3IXLlwIf39/1KtXD/369YO9vT3i4+Nx+vRpREZGAnjy+UAdOnRA06ZNMWDAAGRnZ2PevHlFPq9q1aph0qRJGDduHJKTk9G1a1eYm5sjKSkJ69atw8CBA/Hpp59i+/btGDZsGAICAuDm5oacnBz88ssv0NXVRbdu3Z47Tk5ODpYvXw7gyW205ORkLFq0CHl5eVofRvnZZ59h5cqV6NChA4YPHw5ra2ssW7YMSUlJWLt2LXR0dCAi6N+/P4yNjbFw4UIAT94evnbtWowYMQI+Pj5an9zs5uaGAQMG4PDhw7C1tcWSJUuQmpqKpUuXPnW++vr6mD59Ovr164fWrVsjMDAQqampmDNnDpydnfHJJ58obb28vAAAw4cPh6+vL3R1ddGzZ8/nviZEqlbWiYvoVfP39xcjIyO5f//+U9uEhISIvr6+1pWSuLg46dGjh1SsWFH09PSkYsWK0rNnT4mLiytwfEle4SlMYVd48p04cUICAwPF3t5e9PX1xc7OTgIDA+XkyZPPnU/+8WPGjJEGDRqItbW16Onpib29vQQEBMjRo0cLtN+zZ4+89957Ym5uLqampuLp6Snz5s3TarN161Zp3ry5GBsbi4WFhfj7+8uZM2e02jztCk++tWvXSosWLcTU1FRMTU3F3d1dQkNDJSEhQUSeXFHp37+/VKtWTYyMjMTa2lreeecd2bp163PPubC3pVtYWEjbtm0LPf78+fPyr3/9S6ysrMTIyEgaNWokGzduVOrnzJlT4KqNiMilS5fEwsJCOnbsqJQ5OTmJn5+fREZGiqenpxgaGoq7u3uBtf3nFZ58q1evlvr164uhoaFYW1tL79695cqVK1ptcnJy5OOPPxYbGxvRaDR8izqRiGhE/nEjmYiIXhlnZ2fUrl1b68MaiejV46ZlIiIiUj0GHiIiIlI9Bh4iIiJSPe7hISIiItXjFR4iIiJSPQYeIiIiUr0if/BgXl4erl27BnNzc35sOREREb1SIoLMzEw4ODhAR+flr88UOfBcu3YNjo6OLz0gERERUVFdvnz5hb4M+J+KHHjMzc2VgYv6zctERERELyIjIwOOjo5K/nhZRQ48+bexLCwsGHiIiIioVJTUNhpuWiYiIiLVY+AhIiIi1WPgISIiItVj4CEiIiLVY+AhIiIi1XvrAk9ycjI0Gg2OHTtW1lNRxMfHo0mTJjAyMkK9evUKbdOmTRuMHDmyVOdFRESkFqUeeEJCQqDRaDBt2jSt8t9///2t/QTn8ePHw9TUFAkJCdi2bVuhbSIiIjBx4sRSnlnZ2bZtG5o1awZzc3PY2dlh7NixyMnJUeqjo6PRpUsX2Nvbw9TUFPXq1UN4eHgZzpiIiF5nZXKFx8jICNOnT8fdu3fLYvhX4tGjRy987Pnz59GiRQs4OTmhfPnyhbaxtrYusQ9fet0dP34cHTt2RPv27REbG4vVq1dj/fr1+Oyzz5Q2+/btg6enJ9auXYsTJ06gX79+CA4OxsaNG8tw5kRE9Loqk8Dj4+MDOzs7TJ069altwsLCCtze+e9//wtnZ2fleUhICLp27YopU6bA1tYWVlZWmDBhAnJycjBmzBhYW1ujcuXKWLp0aYH+4+Pj0axZMxgZGaF27drYuXOnVv2pU6fQoUMHmJmZwdbWFn369MGtW7eU+jZt2mDYsGEYOXIkKlSoAF9f30LPIy8vDxMmTEDlypVhaGiIevXq4a+//lLqNRoNYmJiMGHCBGg0GoSFhRXazz9vaTk7O2PSpEkIDg6GmZkZnJycsH79ety8eRNdunSBmZkZPD09ceTIEeWY27dvIzAwEJUqVYKJiQnq1KmDlStXao2TmZmJ3r17w9TUFPb29vj2228LjP3w4UN8+umnqFSpEkxNTdG4cWNER0cr9RcvXoS/vz/KlSsHU1NT1KpVC5s2bSr0vAqzevVqeHp64quvvkL16tXRunVrzJgxA/Pnz0dmZiYA4PPPP8fEiRPRrFkzVKtWDSNGjED79u0RERFR5HGIiOjtUSaBR1dXF1OmTMG8efNw5cqVl+pr+/btuHbtGnbt2oXZs2dj/Pjx6NSpE8qVK4eDBw9i8ODBGDRoUIFxxowZg9GjRyM2NhZNmzaFv78/bt++DQBIS0vDu+++i/r16+PIkSP466+/kJqaiu7du2v1sWzZMhgYGGDv3r1YtGhRofObM2cOZs2ahW+++QYnTpyAr68vOnfujMTERABASkoKatWqhdGjRyMlJQWffvppkc/922+/RfPmzREbGws/Pz/06dMHwcHBCAoKwtGjR1GtWjUEBwdDRAAADx48gJeXF/7880+cOnUKAwcORJ8+fXDo0CGlz1GjRmHv3r1Yv349oqKisHv3bhw9elRr3GHDhmH//v1YtWoVTpw4gYCAALRv3145p9DQUDx8+BC7du3CyZMnMX36dJiZmSnHOzs7PzXYAU8ClZGRkVaZsbExHjx4gJiYmKcel56eDmtr6yK/fkRE9BaRIkpPTxcAkp6eXtRDCtW3b1/p0qWLiIg0adJE+vfvLyIi69atk79PZ/z48VK3bl2tY7/99ltxcnLS6svJyUlyc3OVsho1akjLli2V5zk5OWJqaiorV64UEZGkpCQBINOmTVPaPH78WCpXrizTp08XEZGJEydKu3bttMa+fPmyAJCEhAQREWndurXUr1//uefr4OAgkydP1ipr2LChDB06VHlet25dGT9+/DP7ad26tYwYMUJ57uTkJEFBQcrzlJQUASBffvmlUrZ//34BICkpKU/t18/PT0aPHi0iIhkZGaKvry9r1qxR6tPS0sTExEQZ++LFi6KrqytXr17V6qdt27Yybtw4ERGpU6eOhIWFPXXMd999V+bNm/fU+sjISNHR0ZEVK1ZITk6OXLlyRVq2bCkAZMWKFYUes3r1ajEwMJBTp049tV8iInpzlFTuyFfk79J6FaZPn4533323WFc1/qlWrVpaXxtva2uL2rVrK891dXVRvnx53LhxQ+u4pk2bKj/r6enB29sbcXFxAJ7sIdmxY4fWVYl858+fh5ubGwDAy8vrmXPLyMjAtWvX0Lx5c63y5s2b4/jx40U8w6fz9PRUfra1tQUA1KlTp0DZjRs3YGdnh9zcXEyZMgW//vorrl69ikePHuHhw4cwMTEBAFy4cAGPHz9Go0aNlD4sLS1Ro0YN5fnJkyeRm5urvAb5Hj58qOw/Gj58OIYMGYItW7bAx8cH3bp105rr0zZm52vXrh1mzpyJwYMHo0+fPjA0NMSXX36J3bt3a611vh07dqBfv3748ccfUatWrWe/aERE9FYq07elt2rVCr6+vhg3blyBOh0dHeVWTL7Hjx8XaKevr6/1XKPRFFqWl5dX5Hndu3cP/v7+OHbsmNYjMTERrVq1UtqZmpoWuc9X4e/nmf8Ot8LK8s995syZmDNnDsaOHYsdO3bg2LFj8PX1LdaG63v37kFXVxcxMTFar01cXBzmzJkDAPjwww9x4cIF9OnTBydPnoS3tzfmzZtXrHMbNWoU0tLScOnSJdy6dQtdunQBAFStWlWr3c6dO+Hv749vv/0WwcHBxRqDiIjeHmX+OTzTpk3Dhg0bsH//fq1yGxsbXL9+XSv0lORn5xw4cED5OScnBzExMfDw8AAANGjQAKdPn4azszOqV6+u9ShOyLGwsICDgwP27t2rVb53717UrFmzZE6kGPbu3YsuXbogKCgIdevWRdWqVXH27FmlvmrVqtDX18fhw4eVsvT0dK029evXR25uLm7cuFHgtbGzs1PaOTo6YvDgwYiIiMDo0aPx448/Fnu+Go0GDg4OMDY2xsqVK+Ho6IgGDRoo9dHR0fDz88P06dMxcODAYvdPRERvjzIPPHXq1EHv3r0xd+5crfI2bdrg5s2bmDFjBs6fP4/58+dj8+bNJTbu/PnzsW7dOsTHxyM0NBR3795F//79ATzZdHvnzh0EBgbi8OHDOH/+PCIjI9GvXz/k5uYWa5wxY8Zg+vTpWL16NRISEvDZZ5/h2LFjGDFiRImdS1G5uroiKioK+/btQ1xcHAYNGoTU1FSl3tzcHH379sWYMWOwY8cOnD59GgMGDICOjo5ytcjNzQ29e/dGcHAwIiIikJSUhEOHDmHq1Kn4888/AQAjR45EZGQkkpKScPToUezYsUMJkwDQtm1bfPfdd8+c68yZM3Hy5EmcPn0aEydOxLRp0zB37lzo6uoCeHIby8/PD8OHD0e3bt1w/fp1XL9+HXfu3Cnpl42IiFSgzAMPAEyYMKHALScPDw8sWLAA8+fPR926dXHo0KGX2uvzT9OmTcO0adNQt25d7NmzB+vXr0eFChUAQLkqk5ubi3bt2qFOnToYOXIkrKysCt1D8izDhw/HqFGjMHr0aNSpUwd//fUX1q9fD1dX1xI7l6L6z3/+gwYNGsDX1xdt2rSBnZ0dunbtqtVm9uzZaNq0KTp16gQfHx80b94cHh4eWu+aWrp0KYKDgzF69GjUqFEDXbt2xeHDh1GlShUAQG5uLkJDQ+Hh4YH27dvDzc0NCxYsUI4/f/681lv8C7N582a0bNkS3t7e+PPPP/HHH39ozXXZsmXIysrC1KlTYW9vrzw++OCDl3+hiIhIdTTyz40yT5GRkQFLS0ukp6fDwsLiVc+LXhP3799HpUqVMGvWLAwYMKCsp0NERG+Jks4dZfouLXr9xMbGIj4+Ho0aNUJ6ejomTJgAAMqmYSIiojcRAw8V8M033yAhIQEGBgbw8vLC7t27ldt9REREbyIGHtJSv379Z36aMRER0Zvotdi0TERERPQqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BAREZHqMfAQERGR6jHwEBERkeox8BAREZHq6ZX1BIiInicxMRGZmZllPQ0iKgHm5uZwdXUt9XEZeIjotZaYmAg3N7eyngZRqbEz02CQlwG+j3mE6/ekrKfzSpw9e7bUQw8DDxG91vKv7CxfvhweHh5lPBuiV8847Sw8dg1Cj69+RraVusJ+XFwcgoKCyuSKLQMPEb0RPDw80KBBg7KeBtGrd00H2AV4uLsDDvXKejaqwU3LREREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6JR54srKycPToUWRlZZV010RERKQCcXFxpZ4TSjzwxMfHw8vLC/Hx8SXdNREREalAUFBQqecE3tIiIiIi1WPgISIiItVj4CEiIiLVY+AhIiIi1WPgISIiItVj4CEiIiLVY+AhIiIi1WPgISIiItVj4CEiIiLVY+AhIiIi1WPgISIiItVj4CEiIiLVY+AhIiIi1WPgISIiItXTK+kOs7OzAQBxcXEl3TURvYXy/5bk/20hojfX33+PS/t3usQDT3JyMgAgKCiopLsmordYcnIymjdvXtbTIKKXkJ8R8n8uzd/pEg88zs7OAIDly5fDw8OjpLsnordMXFwcgoKClL8tRPTm+vvvcWn/Tpd44DE2NgYAeHh4oEGDBiXdPRG9pfL/thDRm+vvv8el/TvNTctERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkeiUeeNzd3RETEwN3d/eS7pqIiIhUYPny5aWeE/RKukMTExM0aNCgpLslIiIilfDw8ICJiUmpjslbWkRERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHol/uWhREQlKSsrCwBw9OjRMp4JUekwTjsLDwBx8fHIvp5X1tMpUXFxcWU2NgMPEb3W4uPjAQAfffRRGc+EqHTYmWkwyMsA38/qhev3pKyn80qYm5uX+pgMPET0WuvatSsAwN3dHSYmJmU7GaJS1LmsJ/CKmJubw9XVtdTH1YhIkeJjRkYGLC0tkZ6eDgsLi1c9LyIiInqLlXTu4KZlIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPQYeIiIiUj0GHiIiIlI9Bh4iIiJSPb2iNhQRAEBGRsYrmwwRERER8P/zRn7+eFlFDjyZmZkAAEdHxxIZmIiIiOh5MjMzYWlp+dL9aKSI0SkvLw/Xrl2Dubk5NBrNU9tlZGTA0dERly9fhoWFxUtPkEoW1+f1xzV6vXF9Xm9cn9dfUddIRJCZmQkHBwfo6Lz8DpwiX+HR0dFB5cqVi9yxhYUF/2V7jXF9Xn9co9cb1+f1xvV5/RVljUriyk4+blomIiIi1WPgISIiItUr8cBjaGiI8ePHw9DQsKS7phLA9Xn9cY1eb1yf1xvX5/VXVmtU5E3LRERERG8q3tIiIiIi1WPgISIiItVj4CEiIiLVY+AhIiIi1XuhwDN//nw4OzvDyMgIjRs3xqFDh5S6UaNGwdraGo6OjggPD9c6bs2aNfD393+5GZNi6tSpaNiwIczNzVGxYkV07doVCQkJWm0ePHiA0NBQlC9fHmZmZujWrRtSU1OV+jt37sDf3x9mZmaoX78+YmNjtY4PDQ3FrFmzSuV81G7atGnQaDQYOXKkUsb1KXtXr15FUFAQypcvD2NjY9SpUwdHjhxR6kUEX331Fezt7WFsbAwfHx8kJiYq9Q8fPkSfPn1gYWEBNzc3bN26Vav/mTNn4uOPPy6181GT3NxcfPnll3BxcYGxsTGqVauGiRMnan23Eten9OzatQv+/v5wcHCARqPB77//rlX/vLUAnvxN6927NywsLGBlZYUBAwbg3r17Sn1ycjJatWoFU1NTtGrVCsnJyVrHd+rUCWvXrn2xE5BiWrVqlRgYGMiSJUvk9OnT8tFHH4mVlZWkpqbK+vXrxdbWVg4fPiwrVqwQIyMjuXnzpoiIpKWliaurq1y8eLG4Q9JT+Pr6ytKlS+XUqVNy7Ngx6dixo1SpUkXu3buntBk8eLA4OjrKtm3b5MiRI9KkSRNp1qyZUj9q1Chp3bq1JCQkyMiRI8XLy0up279/v3h5eUlOTk6pnpcaHTp0SJydncXT01NGjBihlHN9ytadO3fEyclJQkJC5ODBg3LhwgWJjIyUc+fOKW2mTZsmlpaW8vvvv8vx48elc+fO4uLiItnZ2SIiMnfuXPHw8JBTp07JzJkzxcbGRvLy8kRE5MKFC+Lq6irp6ellcn5vusmTJ0v58uVl48aNkpSUJGvWrBEzMzOZM2eO0obrU3o2bdokX3zxhURERAgAWbdunVb989ZCRKR9+/ZSt25dOXDggOzevVuqV68ugYGBSv0HH3wgPXv2lLNnz0r37t2lW7duSt2qVavE39//hedf7MDTqFEjCQ0NVZ7n5uaKg4ODTJ06VaZPny49evRQ6ipWrCiHDh0SEZGBAwfK7NmzX3ii9Hw3btwQALJz504ReRIy9fX1Zc2aNUqbuLg4ASD79+8XEZEOHTrIwoULRUTkzJkzYmJiIiIijx49krp168rhw4dL+SzUJzMzU1xdXSUqKkpat26tBB6uT9kbO3astGjR4qn1eXl5YmdnJzNnzlTK0tLSxNDQUFauXCkiIkOGDJGxY8eKiEhWVpYAkBs3bojIk/8piYiIeIVnoG5+fn7Sv39/rbIPPvhAevfuLSJcn7L0z8BTlLU4c+aMAND6u7V582bRaDRy9epVERHx8PCQzZs3i8iTgFWzZk0REbl7965Ur15dLl269MJzLtYtrUePHiEmJgY+Pj5KmY6ODnx8fLB//37UrVsXR44cwd27dxETE4Ps7GxUr14de/bswdGjRzF8+PAXuwxFRZKeng4AsLa2BgDExMTg8ePHWuvl7u6OKlWqYP/+/QCAunXrYvv27cjJyUFkZCQ8PT0BADNmzECbNm3g7e1dymehPqGhofDz89NaB4Dr8zpYv349vL29ERAQgIoVK6J+/fr48ccflfqkpCRcv35da40sLS3RuHFjrTXas2cPsrOzERkZCXt7e1SoUAHh4eEwMjLC+++/X+rnpRbNmjXDtm3bcPbsWQDA8ePHsWfPHnTo0AEA1+d1UpS12L9/P6ysrLT+bvn4+EBHRwcHDx4E8GS9tm7diry8PGzZskX5mzdmzBiEhobC0dHxxSdZnHR09epVASD79u3TKh8zZow0atRIRETGjx8v1apVk9q1a0tERIQ8fPhQateuLUeOHJF58+aJm5ubNGvWTE6dOvXCKY0Kys3NFT8/P2nevLlSFh4eLgYGBgXaNmzYUP7973+LyJMEHhgYKFWqVJFWrVrJ6dOn5ezZs+Lq6iq3bt2SQYMGiYuLiwQEBEhaWlqpnY9arFy5UmrXrq1c0v37FR6uT9kzNDQUQ0NDGTdunBw9elS+//57MTIykp9//llERPbu3SsA5Nq1a1rHBQQESPfu3UXkydW2oUOHirOzs3h7e8vu3bvl9u3bUrVqVbl06ZJ88cUXUq1aNWnXrp1cuXKl1M/xTZabmytjx44VjUYjenp6otFoZMqUKUo916fs4B9XeIqyFpMnTxY3N7cCfdnY2MiCBQtEROTKlSvi5+cnjo6O4ufnJ1euXJGdO3eKt7e33L59WwICAsTFxUUGDRokDx8+LNaci/xt6UUVFhaGsLAw5fnXX38NHx8f6OvrY9KkSTh58iQ2btyI4OBgxMTElPTwb63Q0FCcOnUKe/bsKdZxlpaWWLFihVbZu+++i5kzZyI8PBwXLlxAQkICPvroI0yYMIEbZIvh8uXLGDFiBKKiomBkZPRCfXB9Xq28vDx4e3tjypQpAID69evj1KlTWLRoEfr27VukPvT19TF//nytsn79+mH48OGIjY3F77//juPHj2PGjBkYPnz4i2+4fAv9+uuvCA8Px4oVK1CrVi0cO3YMI0eOhIODA9dHpSpVqoSNGzcqzx8+fAhfX18sW7YMkyZNgrm5ORISEtC+fXt8//33xdpwXqxbWhUqVICurq7Wu0gAIDU1FXZ2dgXax8fHY/ny5Zg4cSKio6PRqlUr2NjYoHv37jh69CgyMzOLMzw9xbBhw7Bx40bs2LEDlStXVsrt7Ozw6NEjpKWlabV/2noBwNKlS2FlZYUuXbogOjoaXbt2hb6+PgICAhAdHf0Kz0J9YmJicOPGDTRo0AB6enrQ09PDzp07MXfuXOjp6cHW1pbrU8bs7e1Rs2ZNrTIPDw9cunQJAJR1KOrfPADYsWMHTp8+jWHDhiE6OhodO3aEqakpunfvzjUqpjFjxuCzzz5Dz549UadOHfTp0weffPIJpk6dCoDr8zopylrY2dnhxo0bWvU5OTm4c+fOU9drypQpaNeuHby8vBAdHY1u3bpBX18fH3zwQbHXq1iBx8DAAF5eXti2bZtSlpeXh23btqFp06ZabUUEgwYNwuzZs2FmZobc3Fw8fvwYAJR/5ubmFmuypE1EMGzYMKxbtw7bt2+Hi4uLVr2Xlxf09fW11ishIQGXLl0qsF4AcPPmTUyYMAHz5s0DgAJrxvUqnrZt2+LkyZM4duyY8vD29kbv3r2Vn7k+Zat58+YFPsrh7NmzcHJyAgC4uLjAzs5Oa40yMjJw8ODBQtco/2MGvv/+e+jq6nKNXlJWVhZ0dLT/M6Wrq4u8vDwAXJ/XSVHWomnTpkhLS9O6u7N9+3bk5eWhcePGBfqMi4vDihUrMHHiRAAl8DevWDfA5MnbwgwNDeXnn3+WM2fOyMCBA8XKykquX7+u1e6HH37QejvZwYMHxcLCQvbv3y9fffWVsvOaXtyQIUPE0tJSoqOjJSUlRXlkZWUpbQYPHixVqlSR7du3y5EjR6Rp06bStGnTQvvr1auXzJs3T3k+ffp08fLykjNnzkiHDh1k6NChr/yc1O7ve3hEuD5l7dChQ6KnpyeTJ0+WxMRECQ8PFxMTE1m+fLnSZtq0aWJlZSV//PGHnDhxQrp06VLgrbb5Pv/8cxk9erTyfPXq1VKlShU5fvy4DBgwQDp27Fgq56UWffv2lUqVKilvS4+IiJAKFSooe9xEuD6lKTMzU2JjYyU2NlYAyOzZsyU2Nlb5uJmirEX79u2lfv36cvDgQdmzZ4+4urpqvS09X15enrRo0UI2bNiglA0ZMkT8/PzkzJkzUr9+fZkxY0ax5l/swCMiMm/ePKlSpYoYGBhIo0aN5MCBA1r1169fFycnJ+VtZvm+/vprsba2Fnd3dzl48OCLDE1/A6DQx9KlS5U22dnZMnToUClXrpyYmJjI+++/LykpKQX6+uuvv6RRo0aSm5urlN2/f18CAgLE3Nxc2rZtK6mpqaVxWqr2z8DD9Sl7GzZskNq1a4uhoaG4u7vLDz/8oFWfl5cnX375pdja2oqhoaG0bdtWEhISCvRz8uRJqV69utbnYOXm5sqQIUPEwsJCGjZsKImJia/8fNQkIyNDRowYIVWqVBEjIyOpWrWqfPHFF1qbVbk+pWfHjh2F/jenb9++IlK0tbh9+7YEBgaKmZmZWFhYSL9+/SQzM7PAWIsWLdK6aCIikpqaKm3bthVzc3MJCAiQ+/fvF2v+GpG/fWQlERERkQrxu7SIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIjojREWFgaNRgONRoP//ve/L9VXmzZtlL6OHTtWIvMjotcXAw+RSuzfvx+6urrw8/MrUBcdHQ2NRlPgm9kBwNnZWSs85IcAjUYDS0tLNG/eHNu3b1fqQ0JC0LVrV63nGo0GgwcPLtB3aGgoNBoNQkJCtMovX76M/v37w8HBAQYGBnBycsKIESNw+/bt555nrVq1kJKSgoEDByplo0aNgrW1NRwdHREeHq7Vfs2aNfD39y/QT0REBA4dOvTc8YhIHRh4iFRi8eLF+Pjjj7Fr1y5cu3btpfpaunQpUlJSsHfvXlSoUAGdOnXChQsXntre0dERq1atQnZ2tlL24MEDrFixAlWqVNFqe+HCBXh7eyMxMRErV67EuXPnsGjRImzbtg1NmzbFnTt3njk3PT092NnZwcTEBACwYcMGrFixAlu2bMGMGTPw4Ycf4tatWwCA9PR0fPHFF5g/f36BfqytrWFjY1Pk14SI3mwMPEQqcO/ePaxevRpDhgyBn58ffv7555fqz8rKCnZ2dqhduzYWLlyI7OxsREVFPbV9gwYN4OjoiIiICKUsIiICVapUQf369bXahoaGwsDAAFu2bEHr1q1RpUoVdOjQAVu3bsXVq1fxxRdfFGuucXFxaNOmDby9vREYGAgLCwskJSUBAP79739jyJAhBUIXEb19GHiIVODXX3+Fu7s7atSogaCgICxZsgQl9b3AxsbGAIBHjx49s13//v2xdOlS5fmSJUvQr18/rTZ37txBZGQkhg4dqvSbz87ODr1798bq1auLNfe6deviyJEjuHv3LmJiYpCdnY3q1atjz549OHr0KIYPH17kvohIvRh4iFRg8eLFCAoKAgC0b98e6enp2Llz50v3m5WVhf/85z/Q1dVF69atn9k2KCgIe/bswcWLF3Hx4kXs3btXmVO+xMREiAg8PDwK7cPDwwN3797FzZs3izxHX19fBAUFoWHDhggJCcGyZctgamqKIUOGYNGiRVi4cCFq1KiB5s2b4/Tp00Xul4jURa+sJ0BELychIQGHDh3CunXrADzZ49KjRw8sXrwYbdq0eaE+AwMDoauri+zsbNjY2GDx4sXw9PR85jE2NjbK7TQRgZ+fHypUqFBo25K6+pQvLCwMYWFhyvOvv/4aPj4+0NfXx6RJk3Dy5Els3LgRwcHBiImJKdGxiejNwMBD9IZbvHgxcnJy4ODgoJSJCAwNDfHdd9/B0tISFhYWAJ5s4rWystI6Pi0tDZaWllpl3377LXx8fGBpaVmsjb39+/fHsGHDAKDQjcLVq1eHRqNBXFwc3n///QL1cXFxKFeu3EttJo6Pj8fy5csRGxuLJUuWoFWrVrCxsUH37t3Rv39/ZGZmwtzc/IX7J6I3E29pEb3BcnJy8L///Q+zZs3CsWPHlMfx48fh4OCAlStXAgBcXV2ho6NT4OrGhQsXkJ6eDjc3N61yOzs7VK9evdjBo3379nj06BEeP34MX1/fAvXly5fHe++9hwULFmi9owsArl+/jvDwcPTo0QMajaZY4+YTEQwaNAizZ8+GmZkZcnNz8fjxYwBQ/pmbm/tCfRPRm41XeIjeYBs3bsTdu3cxYMCAAldpunXrhsWLF2Pw4MEwNzfHhx9+iNGjR0NPTw916tTB5cuXMXbsWDRp0gTNmjUrkfno6uoiLi5O+bkw3333HZo1awZfX19MmjQJLi4uOH36NMaMGYNKlSph8uTJLzz+Tz/9BBsbG+Vzd5o3b46wsDAcOHAAmzdvRs2aNQtc4SKitwOv8BC9wRYvXqzcevqnbt264ciRIzhx4gQAYM6cOejbty/Gjh2LWrVqISQkBJ6entiwYcMLX1EpjIWFhXILrTCurq44cuQIqlatiu7du6NatWoYOHAg3nnnHezfvx/W1tYvNG5qaiomT56MuXPnKmWNGjXC6NGj4efnh19//VXrXWRE9HbRSEnvHiQiekXCwsLw+++/l9hXQSQnJ8PFxQWxsbGoV69eifRJRK8nXuEhojfKyZMnYWZmhgULFrxUPx06dECtWrVKaFZE9LrjFR4iemPcuXNH+eoJGxubQm/lFdXVq1eVjdNVqlSBgYFBicyRiF5PDDxERESkerylRURERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESq93/Jtq5JqKbqVQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 700x200 with 1 Axes>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(7, 2))\n",
    "boxplot_data = ax.boxplot(\n",
    "    aupimo_result.aupimos[labels == 1].numpy(),\n",
    "    vert=False,\n",
    "    widths=0.4,\n",
    ")\n",
    "_ = ax.set_yticks([])\n",
    "ax.set_xlim(0 - (eps := 2e-2), 1 + eps)\n",
    "ax.xaxis.set_major_formatter(PercentFormatter(1))\n",
    "ax.set_xlabel(\"AUPIMO [%]\")\n",
    "ax.set_title(\"AUPIMO Scores Boxplot\")\n",
    "num_images = (labels == 1).sum().item()\n",
    "ax.annotate(\n",
    "    text=f\"Number of images: {num_images}\",\n",
    "    xy=(0.03, 0.95),\n",
    "    xycoords=\"axes fraction\",\n",
    "    xytext=(0, 0),\n",
    "    textcoords=\"offset points\",\n",
    "    annotation_clip=False,\n",
    "    verticalalignment=\"top\",\n",
    ")\n",
    "fig  # noqa: B018, RUF100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the values in the boxplot (e.g., whiskers, quartiles, etc.), we're going to use `matplotlib`'s internal function `mpl.cbook.boxplot_stats()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['mean', 'iqr', 'cilo', 'cihi', 'whishi', 'whislo', 'fliers', 'q1', 'med', 'q3'])\n"
     ]
    }
   ],
   "source": [
    "boxplot_data = mpl.cbook.boxplot_stats(aupimo_result.aupimos[labels == 1].numpy())[0]\n",
    "print(boxplot_data.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll select 5 of those and find images in the dataset that match them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  statistic  value  image_index\n",
      "0    whislo   0.00           65\n",
      "1        q1   0.53           58\n",
      "2       med   0.89           63\n",
      "3        q3   1.00           22\n",
      "4    whishi   1.00            0\n"
     ]
    }
   ],
   "source": [
    "image_selection = []\n",
    "\n",
    "for key in [\"whislo\", \"q1\", \"med\", \"q3\", \"whishi\"]:\n",
    "    value = boxplot_data[key]\n",
    "    # find the image that is closest to the value of the statistic\n",
    "    #     `[labels == 1]` is not used here so that the image's\n",
    "    #         indexes are the same as the ones in the dataset\n",
    "    #     we use `sort()` -- instead of `argmin()` -- so that\n",
    "    #         the `nan`s are not considered (they are at the end)\n",
    "    closest_image_index = (aupimo_result.aupimos - value).abs().argsort()[0]\n",
    "    image_selection.append({\"statistic\": key, \"value\": value, \"image_index\": closest_image_index.item()})\n",
    "\n",
    "image_selection = pd.DataFrame(image_selection)\n",
    "print(image_selection)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that they are sorted from the worst to the best AUPIMO score."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing the Representative Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize what the heatmaps of these samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# will be used to normalize the anomaly maps to fit a colormap\n",
    "global_vmin, global_vmax = torch.quantile(anomaly_maps, torch.tensor([0.02, 0.98]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 5, figsize=(16, 7), layout=\"constrained\")\n",
    "\n",
    "for ax_column, (_, row) in zip(axes.T, image_selection.iterrows(), strict=False):\n",
    "    ax_above, ax_below = ax_column\n",
    "    image = cv2.resize(read_image(image_paths[row.image_index]), (256, 256))\n",
    "    anomaly_map = anomaly_maps[row.image_index].numpy()\n",
    "    mask = masks[row.image_index].squeeze().numpy()\n",
    "    ax_above.imshow(image)\n",
    "    ax_above.contour(mask, levels=[0.5], colors=\"magenta\", linewidths=1)\n",
    "    ax_below.imshow(image)\n",
    "    ax_below.imshow(anomaly_map, cmap=\"jet\", vmin=global_vmin, vmax=global_vmax, alpha=0.30)\n",
    "    ax_below.contour(mask, levels=[0.5], colors=\"magenta\", linewidths=1)\n",
    "    ax_above.set_title(f\"{row.statistic}: {row.value:.0%} AUPIMO  image {row.image_index}\")\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "axes[0, 0].set_ylabel(\"Image + GT Mask\")\n",
    "axes[1, 0].set_ylabel(\"Image + GT Mask + Anomaly Map\")\n",
    "fig.text(\n",
    "    0.03,\n",
    "    -0.01,\n",
    "    \"Magenta: contours of the ground truth (GT) mask. \"\n",
    "    \"Anomaly maps colored in JET colormap with global (across all images) min-max normalization.\",\n",
    "    ha=\"left\",\n",
    "    va=\"top\",\n",
    "    fontsize=\"small\",\n",
    "    color=\"dimgray\",\n",
    ")\n",
    "\n",
    "fig.suptitle(\"Anomalous samples from AUPIMO boxplot's statistics\")\n",
    "fig  # noqa: B018, RUF100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The heatmaps give the impression that all samples are properly detected, right?\n",
    "\n",
    "Notice that the lowest AUPIMO (left) is 0, but the heatmap is (contradictorily) showing a good detection.\n",
    "\n",
    "Why is that?\n",
    "\n",
    "These heatmaps are colored with a gradient from the minimum to the maximum value in all the heatmaps from the test set.\n",
    "\n",
    "This is not taking into account the contraints (FPR restriction) in AUPIMO.\n",
    "\n",
    "Let's compare with the heatmaps from some normal images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 5, figsize=(16, 7), layout=\"constrained\")\n",
    "\n",
    "# random selection of normal images\n",
    "rng = np.random.default_rng(42)\n",
    "normal_images_selection = rng.choice(np.where(labels == 0)[0], size=5, replace=False)\n",
    "\n",
    "for ax_column, index in zip(axes.T, normal_images_selection, strict=False):\n",
    "    ax_above, ax_below = ax_column\n",
    "    image = cv2.resize(read_image(image_paths[index]), (256, 256))\n",
    "    anomaly_map = anomaly_maps[index].numpy()\n",
    "    ax_above.imshow(image)\n",
    "    ax_below.imshow(image)\n",
    "    ax_below.imshow(anomaly_map, cmap=\"jet\", vmin=global_vmin, vmax=global_vmax, alpha=0.30)\n",
    "    ax_above.set_title(f\"image {index}\")\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "axes[0, 0].set_ylabel(\"Image\")\n",
    "axes[1, 0].set_ylabel(\"Image + Anomaly Map\")\n",
    "fig.text(\n",
    "    0.03,\n",
    "    -0.01,\n",
    "    \"Anomaly maps colored in JET colormap with global (across all images) min-max normalization.\",\n",
    "    ha=\"left\",\n",
    "    va=\"top\",\n",
    "    fontsize=\"small\",\n",
    "    color=\"dimgray\",\n",
    ")\n",
    "\n",
    "fig.suptitle(\"Normal samples (test set)\")\n",
    "fig  # noqa: B018, RUF100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the normal images also have high anomaly scores (\"hot\" colors) although there is no anomaly.\n",
    "\n",
    "As a matter of fact, the heatmaps can barely differentiate between some normal and anomalous images.\n",
    "\n",
    "See the two heatmaps below for instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(7, 4), layout=\"constrained\")\n",
    "\n",
    "for ax, index in zip(axes.flatten(), [87, 65], strict=False):\n",
    "    image = cv2.resize(read_image(image_paths[index]), (256, 256))\n",
    "    anomaly_map = anomaly_maps[index].numpy()\n",
    "    mask = masks[index].squeeze().numpy()\n",
    "    ax.imshow(image)\n",
    "    ax.contour(mask, levels=[0.5], colors=\"magenta\", linewidths=1)\n",
    "    ax.imshow(anomaly_map, cmap=\"jet\", vmin=global_vmin, vmax=global_vmax, alpha=0.30)\n",
    "    ax.set_title(f\"image {index}\")\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "axes[0].set_title(f\"{axes[0].get_title()} (normal)\")\n",
    "axes[1].set_title(f\"{axes[1].get_title()} (anomalous)\")\n",
    "\n",
    "fig.text(\n",
    "    0.03,\n",
    "    -0.01,\n",
    "    \"Magenta: contours of the ground truth (GT) mask.\\n\"\n",
    "    \"Anomaly maps colored in JET colormap with global (across all images) min-max normalization.\",\n",
    "    ha=\"left\",\n",
    "    va=\"top\",\n",
    "    fontsize=\"small\",\n",
    "    color=\"dimgray\",\n",
    ")\n",
    "\n",
    "fig.suptitle(\"Normal vs. Anomalous Samples\")\n",
    "fig  # noqa: B018, RUF100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One would expect image 65 (anomalous) to a 'hotter' heatmap than image 87 (normal), but it is the opposite.\n",
    "\n",
    "This shows that the model is not doing a great job."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing the AUPIMO on the Heatmaps\n",
    "\n",
    "We will create another visualization to link the heatmaps to AUPIMO.\n",
    "\n",
    "Recall that AUPIMO computes this integral (simplified):\n",
    "\n",
    "$$\n",
    "    \\int_{\\log(L)}^{\\log(U)} \n",
    "    \\operatorname{TPR}^{i}\\left( \\operatorname{FRP^{-1}}( z ) \\right)\n",
    "    \\, \n",
    "    \\mathrm{d}\\log(z)   \n",
    "$$\n",
    "\n",
    "The integration bounds -- $L$[ower] and $U$[pper] -- are FPR values.\n",
    "\n",
    "> More details about their meaning in the next notebook.\n",
    "\n",
    "We will leverage these two bounds to create a heatmap that shows them in a gradient like this:\n",
    "\n",
    "![Visualization of AUPIMO on the heatmaps](./pimo_viz.svg)\n",
    "\n",
    "If the anomaly score is\n",
    "1. too low (below the lowest threshold of AUPIMO) $\\rightarrow$ not shown; \n",
    "2. between the bounds $\\rightarrow$ shown in a JET gradient;\n",
    "3. too high (above the highest threshold of AUPIMO) $\\rightarrow$ shown in a single color.\n",
    "\n",
    "> Technical detail: lower/upper bound of FPR correspond to the upper/lower bound of threshold.\n",
    "\n",
    "> **Why low values are not shown?**\n",
    ">\n",
    "> Because the values below the lower (threshold) bound would  _never_  be seen as \"anomalous\" by the metric.\n",
    ">\n",
    "> Analogously, high values are shown in red because they are  _always_  seen as \"anomalous\" by the metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FPR bounds\n",
      "Lower bound: 0.00001\n",
      "Upper bound: 0.00010\n",
      "Thresholds corresponding to the FPR bounds\n",
      "Lower threshold: 0.504\n",
      "Upper threshold: 0.553\n"
     ]
    }
   ],
   "source": [
    "# the fpr bounds are fixed in advance in the metric object\n",
    "print(f\"\"\"FPR bounds\n",
    "Lower bound: {aupimo.fpr_bounds[0]:.5f}\n",
    "Upper bound: {aupimo.fpr_bounds[1]:.5f}\"\"\")\n",
    "\n",
    "# their corresponding thresholds depend on the model's behavior\n",
    "# so they only show in the result object\n",
    "print(f\"\"\"Thresholds corresponding to the FPR bounds\n",
    "Lower threshold: {aupimo_result.thresh_lower_bound:.3g}\n",
    "Upper threshold: {aupimo_result.thresh_upper_bound:.3g}\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we re-sample other normal images\n",
    "#    the FPR bounds are so strict that the heatmaps in the normal images\n",
    "#    become almost invisible with this colormap\n",
    "max_anom_score_per_image = anomaly_maps.max(dim=2).values.max(dim=1).values  # noqa: PD011\n",
    "normal_images_with_highest_max_score = sorted(\n",
    "    zip(max_anom_score_per_image[labels == 0], torch.where(labels == 0)[0], strict=False),\n",
    "    reverse=True,\n",
    "    key=lambda x: x[0],\n",
    ")\n",
    "normal_images_with_highest_max_score = [idx.item() for _, idx in normal_images_with_highest_max_score[:5]]\n",
    "\n",
    "fig, axes = plt.subplots(2, 5, figsize=(16, 7), layout=\"constrained\")\n",
    "\n",
    "for ax, (_, row) in zip(axes[0], image_selection.iterrows(), strict=False):\n",
    "    image = cv2.resize(read_image(image_paths[row.image_index]), (256, 256))\n",
    "    anomaly_map = anomaly_maps[row.image_index].numpy()\n",
    "    mask = masks[row.image_index].squeeze().numpy()\n",
    "    ax.imshow(image)\n",
    "    #\n",
    "    # where the magic happens!\n",
    "    #\n",
    "    ax.imshow(\n",
    "        # anything below the lower threshold is set to `nan` so it's not shown\n",
    "        #    because such values would never be detected as anomalies with AUPIMO's contraints\n",
    "        np.where(anomaly_map < aupimo_result.thresh_lower_bound, np.nan, anomaly_map),\n",
    "        cmap=\"jet\",\n",
    "        alpha=0.50,\n",
    "        # notice that vmin/vmax changed here to use the thresholds from the result object\n",
    "        vmin=aupimo_result.thresh_lower_bound,\n",
    "        vmax=aupimo_result.thresh_upper_bound,\n",
    "    )\n",
    "    ax.contour(anomaly_map, levels=[aupimo_result.thresh_lower_bound], colors=[\"blue\"], linewidths=1)\n",
    "    ax.contour(mask, levels=[0.5], colors=\"magenta\", linewidths=1)\n",
    "    ax.set_title(f\"{row.statistic}: {row.value:.0%}AUPIMO  image {row.image_index}\")\n",
    "\n",
    "for ax, index in zip(axes[1], normal_images_with_highest_max_score, strict=False):\n",
    "    image = cv2.resize(read_image(image_paths[index]), (256, 256))\n",
    "    anomaly_map = anomaly_maps[index].numpy()\n",
    "    mask = masks[index].squeeze().numpy()\n",
    "    ax.imshow(image)\n",
    "    ax.imshow(\n",
    "        np.where(anomaly_map < aupimo_result.thresh_lower_bound, np.nan, anomaly_map),\n",
    "        cmap=\"jet\",\n",
    "        alpha=0.30,\n",
    "        vmin=aupimo_result.thresh_lower_bound,\n",
    "        vmax=aupimo_result.thresh_upper_bound,\n",
    "    )\n",
    "    ax.contour(anomaly_map, levels=[aupimo_result.thresh_lower_bound], colors=[\"blue\"], linewidths=1)\n",
    "    ax.set_title(f\"image {index}\")\n",
    "\n",
    "for ax in axes.flatten():\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "axes[0, 0].set_ylabel(\"Anomalous\")\n",
    "axes[1, 0].set_ylabel(\"Normal\")\n",
    "fig.text(\n",
    "    0.03,\n",
    "    -0.01,\n",
    "    \"Magenta: contours of the ground truth (GT) mask. \"\n",
    "    \"Anomaly maps colored in JET colormap between the thresholds in AUPIMO's integral. \"\n",
    "    \"Lower values are transparent, higher values are red.\",\n",
    "    ha=\"left\",\n",
    "    va=\"top\",\n",
    "    fontsize=\"small\",\n",
    "    color=\"dimgray\",\n",
    ")\n",
    "\n",
    "fig.suptitle(\"Visualization linked to AUPIMO's bounds\")\n",
    "fig  # noqa: B018, RUF100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the AUPIMO scores make sense with what you see in the heatmaps.\n",
    "\n",
    "The samples on the left and right are special cases: \n",
    "- left (0% AUPIMO): nothing is seen because the model completely misses the anomaly\\*;\n",
    "- right (100% AUPIMO): is practically red only because the detected the anomaly very well.  \n",
    "\n",
    "\\* Because the scores in image 65 are as low as those in normal images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cite Us\n",
    "\n",
    "AUPIMO was developed during [Google Summer of Code 2023 (GSoC 2023)](https://summerofcode.withgoogle.com/archive/2023/projects/SPMopugd) with the `anomalib` team from Intel's OpenVINO Toolkit.\n",
    "\n",
    "arXiv: [arxiv.org/abs/2401.01984](https://arxiv.org/abs/2401.01984) (accepted to BMVC 2024)\n",
    "\n",
    "Official repository: [github.com/jpcbertoldo/aupimo](https://github.com/jpcbertoldo/aupimo) (numpy-only API and numba-accelerated versions available)\n",
    "\n",
    "```bibtex\n",
    "@misc{bertoldo2024aupimo,\n",
    "      author={Joao P. C. Bertoldo and Dick Ameln and Ashwin Vaidya and Samet Akçay},\n",
    "      title={{AUPIMO: Redefining Visual Anomaly Detection Benchmarks with High Speed and Low Tolerance}}, \n",
    "      year={2024},\n",
    "      url={https://arxiv.org/abs/2401.01984}, \n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Utils\n",
    "\n",
    "Here we provide some utility functions to reproduce the techniques shown in this notebook.\n",
    "\n",
    "They are `numpy` compatible and cover edge cases not discussed here (check the examples)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Representative samples from the boxplot's statistics\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import ndarray\n",
    "from torch import Tensor\n",
    "\n",
    "\n",
    "def _validate_tensor_or_ndarray(x: Tensor | ndarray) -> ndarray:\n",
    "    if not isinstance(x, Tensor | ndarray):\n",
    "        msg = f\"Expected argument to be a tensor or ndarray, but got {type(x)}.\"\n",
    "        raise TypeError(msg)\n",
    "\n",
    "    if isinstance(x, Tensor):\n",
    "        x = x.cpu().numpy()\n",
    "\n",
    "    return x\n",
    "\n",
    "\n",
    "def _validate_values(values: ndarray) -> None:\n",
    "    if values.ndim != 1:\n",
    "        msg = f\"Expected argument `values` to be a 1D, but got {values.ndim}D.\"\n",
    "        raise ValueError(msg)\n",
    "\n",
    "\n",
    "def _validate_labels(labels: ndarray) -> ndarray:\n",
    "    if labels.ndim != 1:\n",
    "        msg = f\"Expected argument `labels` to be a 1D, but got {labels.ndim}D.\"\n",
    "        raise ValueError(msg)\n",
    "\n",
    "    # if torch.is_floating_point(labels):\n",
    "    if np.issubdtype(labels.dtype, np.floating):\n",
    "        msg = f\"Expected argument `labels` to be of int or binary types, but got float: {labels.dtype}.\"\n",
    "        raise TypeError(msg)\n",
    "\n",
    "    # check if it is binary and convert to int\n",
    "    if np.issubdtype(labels.dtype, np.bool_):\n",
    "        labels = labels.astype(int)\n",
    "\n",
    "    unique_values = np.unique(labels)\n",
    "    nor_0_nor_1 = (unique_values != 0) & (unique_values != 1)\n",
    "    if nor_0_nor_1.any():\n",
    "        msg = f\"Expected argument `labels` to have 0s and 1s as ground truth labels, but got values {unique_values}.\"\n",
    "        raise ValueError(msg)\n",
    "\n",
    "    return labels\n",
    "\n",
    "\n",
    "def boxplot_stats(\n",
    "    values: Tensor | ndarray,\n",
    "    labels: Tensor | ndarray,\n",
    "    only_label: int | None = 1,\n",
    "    flier_policy: str | None = None,\n",
    "    repeated_policy: str | None = \"avoid\",\n",
    ") -> list[dict[str, str | int | float | None]]:\n",
    "    \"\"\"Compute boxplot statistics of `values` and find the samples that are closest to them.\n",
    "\n",
    "    This function uses `matplotlib.cbook.boxplot_stats`, which is the same function used by `matplotlib.pyplot.boxplot`.\n",
    "\n",
    "    Args:\n",
    "        values (Tensor | ndarray): Values to compute boxplot statistics from.\n",
    "        labels (Tensor | ndarray): Labels of the samples (0=normal, 1=anomalous). Must have the same shape as `values`.\n",
    "        only_label (int | None): If 0 or 1, only use samples of that class. If None, use both. Defaults to 1.\n",
    "        flier_policy (str | None): What happens with the fliers ('outliers')?\n",
    "                                        - None: Do not include fliers.\n",
    "                                        - 'high': Include only high fliers.\n",
    "                                        - 'low': Include only low fliers.\n",
    "                                        - 'both': Include both high and low fliers.\n",
    "                                    Defaults to None.\n",
    "        repeated_policy (str | None): What happens if a sample has already selected [for another statistic]?\n",
    "                                        - None: Don't care, repeat the sample.\n",
    "                                        - 'avoid': Avoid selecting the same one, go to the next closest.\n",
    "                                        Defaults to 'avoid'.\n",
    "\n",
    "    Returns:\n",
    "        list[dict[str, str | int | float | None]]: List of boxplot statistics.\n",
    "            Keys:\n",
    "                - 'statistic' (str): Name of the statistic.\n",
    "                - 'value' (float): Value of the statistic (same units as `values`).\n",
    "                - 'nearest' (float): Value of the sample in `values` that is closest to the statistic.\n",
    "                            Some statistics (e.g. 'mean') are not guaranteed to be a value in `values`.\n",
    "                            This value is the actual one when they that is the case.\n",
    "                - 'index': Index in `values` that has the `nearest` value to the statistic.\n",
    "    \"\"\"\n",
    "    # operate on numpy arrays only for simplicity\n",
    "    values = _validate_tensor_or_ndarray(values)  # (N,)\n",
    "    labels = _validate_tensor_or_ndarray(labels)  # (N,)\n",
    "\n",
    "    # validate the arguments\n",
    "    _validate_values(values)\n",
    "    labels = _validate_labels(labels)\n",
    "    if values.shape != labels.shape:\n",
    "        msg = (\n",
    "            \"Expected arguments `values` and `labels` to have the same shape, \"\n",
    "            f\"but got {values.shape=} and {labels.shape=}.\"\n",
    "        )\n",
    "        raise ValueError(msg)\n",
    "    assert only_label in {None, 0, 1}, f\"Invalid argument `only_label`: {only_label}\"\n",
    "    assert flier_policy in {None, \"high\", \"low\", \"both\"}, f\"Invalid argument `flier_policy`: {flier_policy}\"\n",
    "    assert repeated_policy in {None, \"avoid\"}, f\"Invalid argument `repeated_policy`: {repeated_policy}\"\n",
    "\n",
    "    if only_label is not None and only_label not in labels:\n",
    "        msg = f\"Argument {only_label=} but `labels` does not contain this class.\"\n",
    "        raise ValueError(msg)\n",
    "\n",
    "    # only consider samples of the given label\n",
    "    # `values` and `labels` now have shape (n,) instead of (N,), where n <= N\n",
    "    label_filter_mask = (labels == only_label) if only_label is not None else np.ones_like(labels, dtype=bool)\n",
    "    values = values[label_filter_mask]  # (n,)\n",
    "    labels = labels[label_filter_mask]  # (n,)\n",
    "    indexes = np.nonzero(label_filter_mask)[0]  # (n,) values are indices in  {0, 1, ..., N-1}\n",
    "\n",
    "    indexes_selected = set()  # values in {0, 1, ..., N-1}\n",
    "\n",
    "    def append(records_: dict, statistic_: str, value_: float) -> None:\n",
    "        indices_sorted_by_distance = np.abs(values - value_).argsort()  # (n,)\n",
    "        candidate = indices_sorted_by_distance[0]  # idx that refers to {0, 1, ..., n-1}\n",
    "\n",
    "        nearest = values[candidate]\n",
    "        index = indexes[candidate]  # index has value in {0, 1, ..., N-1}\n",
    "        label = labels[candidate]\n",
    "\n",
    "        if index in indexes_selected and repeated_policy == \"avoid\":\n",
    "            for candidate in indices_sorted_by_distance:\n",
    "                index_of_candidate = indexes[candidate]\n",
    "                if index_of_candidate in indexes_selected:\n",
    "                    continue\n",
    "                # if the code reaches here, it means that `index_of_candidate` is not repeated\n",
    "                # if this is never reached, the first choice will be kept\n",
    "                nearest = values[candidate]\n",
    "                label = labels[candidate]\n",
    "                index = index_of_candidate\n",
    "                break\n",
    "\n",
    "        indexes_selected.add(index)\n",
    "\n",
    "        records_.append(\n",
    "            {\n",
    "                \"statistic\": statistic_,\n",
    "                \"value\": float(value_),\n",
    "                \"nearest\": float(nearest),\n",
    "                \"index\": int(index),\n",
    "                \"label\": int(label),\n",
    "            },\n",
    "        )\n",
    "\n",
    "    # function used in `matplotlib.boxplot`\n",
    "    boxplot_stats = mpl.cbook.boxplot_stats(values)[0]  # [0] is for the only boxplot\n",
    "\n",
    "    records = []\n",
    "    for stat, val in boxplot_stats.items():\n",
    "        if stat in {\"iqr\", \"cilo\", \"cihi\"}:\n",
    "            continue\n",
    "\n",
    "        if stat != \"fliers\":\n",
    "            append(records, stat, val)\n",
    "            continue\n",
    "\n",
    "        if flier_policy is None:\n",
    "            continue\n",
    "\n",
    "        for val_ in val:\n",
    "            stat_ = \"flierhi\" if val_ > boxplot_stats[\"med\"] else \"flierlo\"\n",
    "            if flier_policy == \"high\" and stat_ == \"flierlo\":\n",
    "                continue\n",
    "            if flier_policy == \"low\" and stat_ == \"flierhi\":\n",
    "                continue\n",
    "            # else means that they match or `fliers == \"both\"`\n",
    "            append(records, stat_, val_)\n",
    "\n",
    "    return sorted(records, key=lambda r: r[\"value\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.00     0.00     65      1\n",
      "1        q1   0.53     0.53     58      1\n",
      "2      mean   0.74     0.75      7      1\n",
      "3       med   0.89     0.89     63      1\n",
      "4        q3   1.00     1.00     22      1\n",
      "5    whishi   1.00     1.00      0      1\n"
     ]
    }
   ],
   "source": [
    "# basic usage\n",
    "boxplot_statistics = boxplot_stats(aupimo_result.aupimos, labels)\n",
    "boxplot_statistics = pd.DataFrame.from_records(boxplot_statistics)\n",
    "print(boxplot_statistics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repeated Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.00     0.00     67      1\n",
      "1        q1   0.59     0.59     58      1\n",
      "2      mean   0.78     0.79     43      1\n",
      "3       med   0.98     0.99      9      1\n",
      "4    whishi   1.00     1.00      0      1\n",
      "5        q3   1.00     1.00     36      1\n"
     ]
    }
   ],
   "source": [
    "# repeated values\n",
    "#   if the distribution is very skewed to one side,\n",
    "#   some statistics may have the same value\n",
    "#   e.g. the Q3 and the high whisker\n",
    "#\n",
    "#   let's simulate this situation\n",
    "\n",
    "# increase all values by 10% and clip to [0, 1]\n",
    "mock = torch.clip(aupimo_result.aupimos.clone() * 1.10, 0, 1)\n",
    "\n",
    "# 'avoid' is the default policy\n",
    "# notice how Q3 and the high whisker have the same value, but different indexes\n",
    "print(pd.DataFrame.from_records(boxplot_stats(mock, labels, repeated_policy=\"avoid\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.00     0.00     67      1\n",
      "1        q1   0.59     0.59     58      1\n",
      "2      mean   0.78     0.79     43      1\n",
      "3       med   0.98     0.99      9      1\n",
      "4    whishi   1.00     1.00      0      1\n",
      "5        q3   1.00     1.00      0      1\n"
     ]
    }
   ],
   "source": [
    "# this behavior can be changed to allow repeated values\n",
    "print(pd.DataFrame.from_records(boxplot_stats(mock, labels, repeated_policy=None)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAADvCAYAAAD2Og4yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3ZUlEQVR4nO3de1yP9/8/8Me7d+eziIpUKIVyqDmfNk1I01iIEDanDGNmPjtozofhY+aw7YP5TA4z2YfGyGkOyyk5J6Eck3NFSPX8/eHX9fVeoUjl8rjfbu/bvF+v1/W6Xtf7lTx2Xa/remtEREBERESkYnqlPQAiIiKiV42Bh4iIiFSPgYeIiIhUj4GHiIiIVI+Bh4iIiFSPgYeIiIhUj4GHiIiIVI+Bh4iIiFSPgYeIiIhUj4GHiKgM2rFjBzQaDXbs2FHaQyFSBQYeeuPMnz8fGo0GjRo1KrA+OTkZGo0G3377bYH13377LTQaDZKTk5Wy1q1bQ6PRKC8bGxu89dZbWLx4MXJzc5V2oaGhMDc31+kvb1tXV9cC9xcdHa30+9tvv+WrP3HiBEJCQlC5cmUYGRnBwcEBPXv2xIkTJ573USiuX7+O4cOHw93dHSYmJqhYsSIaNmyIMWPG4O7du4Xu53URGhqqM1/6+vpwdHRE9+7dcfLkydIe3kvbsGEDwsPDS3sYRGWKfmkPgKikRUREwNnZGfv378eZM2dQo0aNYum3SpUqmDJlCoDHAeK///0v+vfvj9OnT2Pq1KnP3NbY2BhnzpzB/v370bBhw3zjNTY2xoMHD/JtFxkZieDgYNjY2KB///5wcXFBcnIyFi1ahN9++w0rV67E+++//8x937p1Cz4+PkhPT0e/fv3g7u6Omzdv4ujRo1iwYAEGDx6cL6SpgZGREf7zn/8AALKzs3H27FksXLgQf/75J06ePAkHB4dSHuGL27BhA+bNm8fQQ/QEBh56oyQlJeHvv/9GZGQkBg4ciIiICIwbN65Y+rayskJISIjyfuDAgahZsya+//57TJgwAQYGBk/dtnr16sjOzsaKFSt0As+DBw+wdu1a+Pv7Y82aNTrbnD17Fr169UK1atWwc+dO2NraKnXDhw9HixYt0KtXLxw9ehTVqlV76r4XLVqECxcuYM+ePWjatKlOXXp6OgwNDQv9Gbyse/fuwczMrET2pa+vrzNfANC4cWN07NgRf/zxBz766KMSGQcRlQxe0qI3SkREBMqVKwd/f3988MEHiIiIeGX7MjU1RePGjXHv3j1cv379ue2Dg4OxatUqnUtg69evR2ZmJrp27Zqv/YwZM5CZmYkff/xRJ+wAQIUKFfDDDz/g3r17mD59+jP3e/bsWWi1WjRu3DhfnaWlJYyNjXXK9u3bhw4dOqBcuXIwMzODl5cX5syZo9Nm27ZtaNGiBczMzGBtbY1OnTohPj5ep014eDg0Gg1OnjyJHj16oFy5cmjevLlSv2zZMnh7e8PExAQ2Njbo3r07Ll68qNNHYmIiunTpAjs7OxgbG6NKlSro3r070tLSnnnMT2NnZwfgcRh60rlz5xAUFAQbGxtlXv/44w+lPj4+HiYmJujdu7fOdrt374ZWq8WYMWOUMmdnZ3Ts2BGbN29GvXr1YGxsjFq1aiEyMrJQY1y9erXyuVSoUAEhISG4fPmyUh8aGop58+YBgM5lO6I3HQMPvVEiIiLQuXNnGBoaIjg4GImJiThw4MAr29+5c+eg1WphbW393LY9evRASkqKziLV5cuXo02bNqhYsWK+9uvXr4ezszNatGhRYH8tW7aEs7Ozzj/MBXFyckJOTg5++eWX544xOjoaLVu2xMmTJzF8+HDMnDkTb7/9NqKiopQ2W7ZsgZ+fH65du4bw8HCMHDkSf//9N5o1a6az7ilPUFAQMjMzMXnyZOWsyqRJk9C7d2+4urpi1qxZGDFiBLZu3YqWLVvizp07AICsrCz4+flh7969+PjjjzFv3jwMGDAA586dU9o8z40bN3Djxg2kpqYiJiYGn3zyCcqXL4+OHTsqbVJTU9G0aVNs2rQJQ4YMwaRJk/DgwQO89957WLt2LQDAw8MDEyZMwC+//IJ169YBeHy2KjQ0FO7u7hg/frzOfhMTE9GtWze0b98eU6ZMgb6+PoKCghAdHf3M8f7888/o2rUrtFotpkyZgo8++giRkZFo3ry5cswDBw7Eu+++CwD45ZdflBfRG0+I3hAHDx4UABIdHS0iIrm5uVKlShUZPny4TrukpCQBIDNmzCiwnxkzZggASUpKUspatWol7u7ucv36dbl+/brEx8fLsGHDBIAEBAQo7fr06SNmZmY6/bVq1Upq164tIiI+Pj7Sv39/ERG5ffu2GBoaytKlS2X79u0CQFavXi0iInfu3BEA0qlTp2ce83vvvScAJD09/altrl69Kra2tgJA3N3dZdCgQbJ8+XK5c+eOTrvs7GxxcXERJycnuX37tk5dbm6u8ud69epJxYoV5ebNm0rZkSNHRE9PT3r37q2UjRs3TgBIcHCwTl/Jycmi1Wpl0qRJOuXHjh0TfX19pTwuLk7nMymKPn36CIB8r8qVK0tsbKxO2xEjRggA2bVrl1KWkZEhLi4u4uzsLDk5OSIikpOTI82bN5dKlSrJjRs3JCwsTPT19eXAgQM6/Tk5OQkAWbNmjVKWlpYm9vb2Ur9+faUsb863b98uIiJZWVlSsWJFqVOnjty/f19pFxUVJQDk66+/VsrCwsKEv96JdPEMD70xIiIiUKlSJbz99tsAHp/u79atG1auXImcnJyX7v/UqVOwtbWFra0tPDw8MHfuXPj7+2Px4sWF7qNHjx6IjIxEVlYWfvvtN2i12gIXHWdkZAAALCwsntlfXn16evpT21SqVAlHjhzBoEGDcPv2bSxcuBA9evRAxYoVMWHCBIgIACAuLg5JSUkYMWJEvjNWeZdMUlJScPjwYYSGhsLGxkap9/LywrvvvosNGzbk2/+gQYN03kdGRiI3Nxddu3ZVzsDcuHEDdnZ2cHV1xfbt2wE8XjMFAJs2bUJmZuYzP4eCGBsbIzo6GtHR0di0aRN++OEHmJubo0OHDjh9+rTSbsOGDWjYsKHO5TZzc3MMGDAAycnJyl1denp6+Pnnn3H37l20b98e8+fPx9ixY+Hj45Nv3w4ODjrzamlpid69eyMuLg5Xr14tcLwHDx7EtWvXMGTIEJ3LjP7+/nB3d3/umTyiNx0DD70RcnJysHLlSrz99ttISkrCmTNncObMGTRq1AipqanYunVrkfv857oIZ2dnREdHY8uWLdi9ezeuXr2KqKgoVKhQodB95q0/2bhxIyIiItCxY8cCQ01eWV7weZrCBiN7e3ssWLAAKSkpSEhIwHfffQdbW1t8/fXXWLRoEYDHa30AoE6dOk/t5/z58wCAmjVr5qvz8PDAjRs3cO/ePZ1yFxcXnfeJiYkQEbi6uioBMu8VHx+Pa9euKduNHDkS//nPf1ChQgX4+flh3rx5hV6/o9Vq4evrC19fX7Rt2xYDBgzAli1bkJaWhrFjx+oc09OO58ljBh4vPg8PD8eBAwdQu3ZtfPXVVwXuu0aNGvl+ftzc3ACgwMt+T+6noLG4u7vrjIOI8uNdWvRG2LZtG1JSUrBy5UqsXLkyX31ERATatm0LAMr/Pd+/f7/AvvLOJvxzMa+ZmRl8fX1fapz29vZo3bo1Zs6ciT179uS7MyuPlZUV7O3tcfTo0Wf2d/ToUVSuXBmWlpaF2r9Go4Gbmxvc3Nzg7+8PV1dXRERE4MMPPyzysRSWiYmJzvvc3FxoNBps3LgRWq02X/snb5GfOXMmQkND8b///Q+bN2/GsGHDMGXKFOzduxdVqlQp8liqVKmCmjVrYufOnUU/kP9v8+bNAIArV67g5s2bykJoIipdPMNDb4SIiAhUrFgRq1evzvcKDg7G2rVrlYBja2sLU1NTJCQkFNhXQkICTE1Ni3Tmpih69OiBXbt2wdLSEh06dHhqu44dOyIpKQm7d+8usH7Xrl1ITk7WWYBbFNWqVUO5cuWQkpIC4PHZCwA4fvz4U7dxcnICgAI/u1OnTqFChQrPve28evXqEBG4uLgoZ2CefP3zbjJPT098+eWX2LlzJ3bt2oXLly9j4cKFRTrWJ2VnZ+s8bNHJyempx5NXn2fhwoWIjo7GpEmTkJWVhYEDBxa4jzNnziiXCvPkXUZzdnYucJtnfbYJCQk64+BdWUQFKN0lRESvXmZmplhYWEi/fv0KrN+zZ48AkJUrVyplgYGBYmlpKefPn9dpe/78ebGwsJDAwECd8icXHj/L8xYtizxekDxu3DhZvny5UvbPRcsiIqdPnxYTExOpVauW3LhxQ6fPmzdvSq1atcTU1FTOnDnzzDHt3btX7t69m6983759AkDee+89EXm8KLewi5YrVaqk0+bYsWNPXbR8/fp1nb7OnDkjWq1WevToodNv3n7yjjUtLU0ePXqkU5+eni56enry6aefPvOYC5oHEZGEhATR19eXRo0aKWV5i5b//vtvpezu3btSrVo1nUXL586dE3Nzc+nSpYuIiCxcuFAAyNKlS3X28axFy/Xq1VPKnrZo2cvLSx48eKC027BhQ75Fy2PGjBEA+eaJ6E3GS1qkeuvWrUNGRgbee++9AusbN24MW1tbREREoFu3bgCAyZMno3HjxmjQoAEGDBgAZ2dnJCcn48cff4RGo8HkyZNf2XitrKwK9YRcV1dXLF26FD179oSnp2e+Jy3fuHEDK1asUM7MPM0vv/yCiIgIvP/++/D29oahoSHi4+OxePFiGBsb41//+heAx4tyFyxYgICAANSrVw99+/aFvb09Tp06hRMnTmDTpk0AHj8fqH379mjSpAn69++P+/fvY+7cuYU+rurVq2PixIkYO3YskpOTERgYCAsLCyQlJWHt2rUYMGAAPv30U2zbtg1Dhw5FUFAQ3NzckJ2djV9++QVarRZdunR57n6ys7OxbNkyAI8voyUnJ2PhwoXIzc3VeRjl559/jhUrVqB9+/YYNmwYbGxssHTpUiQlJWHNmjXQ09ODiKBfv34wMTHBggULADy+PXzNmjUYPnw4fH19dZ7c7Obmhv79++PAgQOoVKkSFi9ejNTUVCxZsuSp4zUwMMC0adPQt29ftGrVCsHBwUhNTcWcOXPg7OyMTz75RGnr7e0NABg2bBj8/Pyg1WrRvXv3534mRKpW2omL6FULCAgQY2NjuXfv3lPbhIaGioGBgc6Zkvj4eOnWrZtUrFhR9PX1pWLFitK9e3eJj4/Pt31xnuEpSEFnePIcPXpUgoODxd7eXgwMDMTOzk6Cg4Pl2LFjzx1P3vajR4+WBg0aiI2Njejr64u9vb0EBQXJoUOH8rXfvXu3vPvuu2JhYSFmZmbi5eUlc+fO1WmzZcsWadasmZiYmIilpaUEBATIyZMnddo87QxPnjVr1kjz5s3FzMxMzMzMxN3dXcLCwiQhIUFEHp9R6devn1SvXl2MjY3FxsZG3n77bdmyZctzj7mg29ItLS2lTZs2BW5/9uxZ+eCDD8Ta2lqMjY2lYcOGEhUVpdTPmTMn31kbEZELFy6IpaWldOjQQSlzcnISf39/2bRpk3h5eYmRkZG4u7vnm9t/nuHJs2rVKqlfv74YGRmJjY2N9OzZUy5duqTTJjs7Wz7++GOxtbUVjUbDW9SJREQj8o8LyURE9Mo4OzujTp06Og9rJKJXj4uWiYiISPUYeIiIiEj1GHiIiIhI9biGh4iIiFSPZ3iIiIhI9Rh4iIiISPUK/eDB3NxcXLlyBRYWFnxsOREREb1SIoKMjAw4ODhAT+/lz88UOvBcuXIFjo6OL71DIiIiosK6ePHiC30Z8D8VOvBYWFgoOy7sNy8TERERvYj09HQ4Ojoq+eNlFTrw5F3GsrS0ZOAhIiKiElFcy2i4aJmIiIhUj4GHiIiIVI+Bh4iIiFSPgYeIiIhUj4GHiIiIVO+NCzzJycnQaDQ4fPhwaQ9FcerUKTRu3BjGxsaoV69egW1at26NESNGlOi4iIiI1KLEA09oaCg0Gg2mTp2qU/7777+/sU9wHjduHMzMzJCQkICtW7cW2CYyMhITJkwo4ZGVnq1bt6Jp06awsLCAnZ0dxowZg+zsbKV+x44d6NSpE+zt7WFmZoZ69eohIiKiFEdMRERlWamc4TE2Nsa0adNw+/bt0tj9K5GVlfXC2549exbNmzeHk5MTypcvX2AbGxubYnv4Ull35MgRdOjQAe3atUNcXBxWrVqFdevW4fPPP1fa/P333/Dy8sKaNWtw9OhR9O3bF71790ZUVFQpjpyIiMqqUgk8vr6+sLOzw5QpU57aJjw8PN/lnX//+99wdnZW3oeGhiIwMBCTJ09GpUqVYG1tjfHjxyM7OxujR4+GjY0NqlSpgiVLluTr/9SpU2jatCmMjY1Rp04d/PXXXzr1x48fR/v27WFubo5KlSqhV69euHHjhlLfunVrDB06FCNGjECFChXg5+dX4HHk5uZi/PjxqFKlCoyMjFCvXj38+eefSr1Go0FsbCzGjx8PjUaD8PDwAvv55yUtZ2dnTJw4Eb1794a5uTmcnJywbt06XL9+HZ06dYK5uTm8vLxw8OBBZZubN28iODgYlStXhqmpKTw9PbFixQqd/WRkZKBnz54wMzODvb09Zs+enW/fDx8+xKefforKlSvDzMwMjRo1wo4dO5T68+fPIyAgAOXKlYOZmRlq166NDRs2FHhcBVm1ahW8vLzw9ddfo0aNGmjVqhWmT5+OefPmISMjAwDwr3/9CxMmTEDTpk1RvXp1DB8+HO3atUNkZGSh90NERG+OUgk8Wq0WkydPxty5c3Hp0qWX6mvbtm24cuUKdu7ciVmzZmHcuHHo2LEjypUrh3379mHQoEEYOHBgvv2MHj0ao0aNQlxcHJo0aYKAgADcvHkTAHDnzh288847qF+/Pg4ePIg///wTqamp6Nq1q04fS5cuhaGhIfbs2YOFCxcWOL45c+Zg5syZ+Pbbb3H06FH4+fnhvffeQ2JiIgAgJSUFtWvXxqhRo5CSkoJPP/200Mc+e/ZsNGvWDHFxcfD390evXr3Qu3dvhISE4NChQ6hevTp69+4NEQEAPHjwAN7e3vjjjz9w/PhxDBgwAL169cL+/fuVPkeOHIk9e/Zg3bp1iI6Oxq5du3Do0CGd/Q4dOhQxMTFYuXIljh49iqCgILRr1045prCwMDx8+BA7d+7EsWPHMG3aNJibmyvbOzs7PzXYAY8DlbGxsU6ZiYkJHjx4gNjY2Kdul5aWBhsbm0J/fkRE9AaRQkpLSxMAkpaWVthNCtSnTx/p1KmTiIg0btxY+vXrJyIia9eulSeHM27cOKlbt67OtrNnzxYnJyedvpycnCQnJ0cpq1mzprRo0UJ5n52dLWZmZrJixQoREUlKShIAMnXqVKXNo0ePpEqVKjJt2jQREZkwYYK0bdtWZ98XL14UAJKQkCAiIq1atZL69es/93gdHBxk0qRJOmVvvfWWDBkyRHlft25dGTdu3DP7adWqlQwfPlx57+TkJCEhIcr7lJQUASBfffWVUhYTEyMAJCUl5an9+vv7y6hRo0REJD09XQwMDGT16tVK/Z07d8TU1FTZ9/nz50Wr1crly5d1+mnTpo2MHTtWREQ8PT0lPDz8qft85513ZO7cuU+t37Rpk+jp6cny5cslOztbLl26JC1atBAAsnz58gK3WbVqlRgaGsrx48ef2i8REb0+iit35Cn0d2m9CtOmTcM777xTpLMa/1S7dm2dr42vVKkS6tSpo7zXarUoX748rl27prNdkyZNlD/r6+vDx8cH8fHxAB6vIdm+fbvOWYk8Z8+ehZubGwDA29v7mWNLT0/HlStX0KxZM53yZs2a4ciRI4U8wqfz8vJS/lypUiUAgKenZ76ya9euwc7ODjk5OZg8eTJ+/fVXXL58GVlZWXj48CFMTU0BAOfOncOjR4/QsGFDpQ8rKyvUrFlTeX/s2DHk5OQon0Gehw8fKuuPhg0bhsGDB2Pz5s3w9fVFly5ddMb6tIXZedq2bYsZM2Zg0KBB6NWrF4yMjPDVV19h165dOnOdZ/v27ejbty9++ukn1K5d+9kfGhERvZFK9bb0li1bws/PD2PHjs1Xp6enp1yKyfPo0aN87QwMDHTeazSaAstyc3MLPa67d+8iICAAhw8f1nklJiaiZcuWSjszM7NC9/kqPHmceXe4FVSWd+wzZszAnDlzMGbMGGzfvh2HDx+Gn59fkRZc3717F1qtFrGxsTqfTXx8PObMmQMA+PDDD3Hu3Dn06tULx44dg4+PD+bOnVukYxs5ciTu3LmDCxcu4MaNG+jUqRMAoFq1ajrt/vrrLwQEBGD27Nno3bt3kfZBRERvjlJ/Ds/UqVOxfv16xMTE6JTb2tri6tWrOqGnOJ+ds3fvXuXP2dnZiI2NhYeHBwCgQYMGOHHiBJydnVGjRg2dV1FCjqWlJRwcHLBnzx6d8j179qBWrVrFcyBFsGfPHnTq1AkhISGoW7cuqlWrhtOnTyv11apVg4GBAQ4cOKCUpaWl6bSpX78+cnJycO3atXyfjZ2dndLO0dERgwYNQmRkJEaNGoWffvqpyOPVaDRwcHCAiYkJVqxYAUdHRzRo0ECp37FjB/z9/TFt2jQMGDCgyP0TEdGbo9QDj6enJ3r27InvvvtOp7x169a4fv06pk+fjrNnz2LevHnYuHFjse133rx5WLt2LU6dOoWwsDDcvn0b/fr1A/B40e2tW7cQHByMAwcO4OzZs9i0aRP69u2LnJycIu1n9OjRmDZtGlatWoWEhAR8/vnnOHz4MIYPH15sx1JYrq6uiI6Oxt9//434+HgMHDgQqampSr2FhQX69OmD0aNHY/v27Thx4gT69+8PPT095WyRm5sbevbsid69eyMyMhJJSUnYv38/pkyZgj/++AMAMGLECGzatAlJSUk4dOgQtm/froRJAGjTpg2+//77Z451xowZOHbsGE6cOIEJEyZg6tSp+O6776DVagE8vozl7++PYcOGoUuXLrh69SquXr2KW7duFffHRkREKlDqgQcAxo8fn++Sk4eHB+bPn4958+ahbt262L9//0ut9fmnqVOnYurUqahbty52796NdevWoUKFCgCgnJXJyclB27Zt4enpiREjRsDa2rrANSTPMmzYMIwcORKjRo2Cp6cn/vzzT6xbtw6urq7FdiyF9eWXX6JBgwbw8/ND69atYWdnh8DAQJ02s2bNQpMmTdCxY0f4+vqiWbNm8PDw0LlrasmSJejduzdGjRqFmjVrIjAwEAcOHEDVqlUBADk5OQgLC4OHhwfatWsHNzc3zJ8/X9n+7NmzOrf4F2Tjxo1o0aIFfHx88Mcff+B///ufzliXLl2KzMxMTJkyBfb29sqrc+fOL/9BERGR6mjknwtlniI9PR1WVlZIS0uDpaXlqx4XlRH37t1D5cqVMXPmTPTv37+0h0NERG+I4s4dpXqXFpU9cXFxOHXqFBo2bIi0tDSMHz8eAJRFw0RERK8jBh7K59tvv0VCQgIMDQ3h7e2NXbt2KZf7iIiIXkcMPKSjfv36z3yaMRER0euoTCxaJiIiInqVGHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj19Et7AERE9HSJiYnIyMgo7WEQFRsLCwu4urqW+H4ZeIiIyqjExES4ubmV9jBIZezMNRjobYgfYrNw9a6UyhhOnz5d4qGHgYeIqIzKO7OzbNkyeHh4lPJoSC1M7pyGx86B6Pb1z7hvXbKBOj4+HiEhIaVy1pKBh4iojPPw8ECDBg1KexikFlf0gJ2Ah7s74FCvtEdTYrhomYiIiFSPgYeIiIhUj4GHiIiIVI+Bh4iIiFSPgYeIiIhUj4GHiIiIVI+Bh4iIiFSPgYeIiIhUj4GHiIiIVI+Bh4iIiFSPgYeIiIhUj4GHXlhmZiYOHTqEzMzM0h4KERG9RuLj40v83w4GHnphp06dgre3N06dOlXaQyEiotdISEhIif/bwcBDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqqdfnJ1lZWVh/vz5OHv2LKpXr44hQ4bA0NAQAJCTk4Ndu3YhJSUF9vb2aNGiBbRabXHu/oWU1XERERFR8Sm2wPPZZ59h9uzZyM7OVspGjx6NTz75BI0bN8aoUaOQnJys1Dk7O2PmzJno3LlzcQ2hyCIjI8vkuIiIiKh4Fcslrc8++wwzZsxA+fLl8dNPPyElJQU//fQTypcvjxkzZqBLly7w9PRETEwMMjIyEBMTA09PT3zwwQeIjIwsjiEUWWRkJD744IMyNy4iIiIqfhoRkcI0TE9Ph5WVFdLS0mBpaamUZ2VlwczMDOXLl8elS5egr/9/J40ePnwIU1NTiAgyMzNhbGys1OXm5iIwMBDHjx9HYmJiiV5GysnJQY0aNeDp6Ynff/8denr/l/tKc1yvm0OHDsHb2xuxsbFo0KBBaQ+HSHX4d4xeiSuHgR9bAQP+Ahzqleiu836mATz35/ppueNFvfQlrfnz5yM7OxsTJ07UCTsAEBMTg9zcXADAwoULMWLECKVOT08PY8eORdOmTbFr1y60bt36ZYdSaLt27UJycjJWrFihE3ZKe1yvm/v37wMA4uPjS3kkROqU93cr7+8a0evuyZ/lkv65funAc/bsWQBAx44d89WlpKTka/ekOnXq5GtXEvL2l7f/fyqtcb1u8tY+hYSElO5AiFQuOTkZzZo1K+1hEL20J9fMlvTP9UsHnurVqwMAoqKi8OGHH+rU2dvb52v3pOPHj+drVxLy9nf8+HE0bty4zIzrdePs7AwAWLZsGTw8PEp3MEQqFB8fj5CQEOXvGtHr7smf5ZL+uX7pwDNkyBCMHj0aX375JUJDQ3UuazVp0gR6enoQEQwaNEhnu9zcXEyZMgUuLi5o0aLFyw6jSFq0aAFnZ2dMnjy5wDU8pTWu142JiQkAwMPDg+sLiF6hvL9rRK+7J3+WS/rn+qXv0jI0NMQnn3yC1NRUVKlSBT/++COuXLmCH3/8EU5OTsjNzYWIoGvXrjp3QwUGBiIqKgrffvttiS8M1mq1mDlzJqKiohAYGFhmxkVERESvRrE8h2f69OkAgNmzZ2PgwIH/17m+PkaPHq08h6dp06ZKnYuLC3777bdSe95N586d8dtvv5W5cREREVHxK7YHD06fPh0TJ0586pOWO3XqVOaeaNy5c+cyOS4iIiIqXsX61RKGhoY6t54/SavVlslbvMvquIiIiKj48MtDiYiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9Rh46IW5u7sjNjYW7u7upT0UIiJ6jSxbtqzE/+0o1q+WoDeLqakpGjRoUNrDICKi14yHhwdMTU1LdJ88w0NERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHLw8lIiqjMjMzAQCHDh0q5ZGQmpjcOQ0PAPGnTuH+1dwS3Xd8fHyJ7u9JDDxERGXUqVOnAAAfffRRKY+E1MTOXIOB3ob4YWYPXL0rpTIGCwuLEt8nAw8RURkVGBgIAHB3d4epqWnpDoZU571S2q+FhQVcXV1LfL8aESlUvEtPT4eVlRXS0tJgaWn5qsdFREREb7Dizh1ctExERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKqnX9iGIgIASE9Pf2WDISIiIgL+L2/k5Y+XVejAk5GRAQBwdHQslh0TERERPU9GRgasrKxeuh+NFDI65ebm4sqVK7CwsIBGo3lqu/T0dDg6OuLixYuwtLR86QFS8eL8lH2co7KN81O2cX7KvsLOkYggIyMDDg4O0NN7+RU4hT7Do6enhypVqhS6Y0tLS/6wlWGcn7KPc1S2cX7KNs5P2VeYOSqOMzt5uGiZiIiIVI+Bh4iIiFSv2AOPkZERxo0bByMjo+LumooB56fs4xyVbZyfso3zU/aV1hwVetEyERER0euKl7SIiIhI9Rh4iIiISPUYeIiIiEj1GHiIiIhI9V4o8MybNw/Ozs4wNjZGo0aNsH//fqVu5MiRsLGxgaOjIyIiInS2W716NQICAl5uxKSYMmUK3nrrLVhYWKBixYoIDAxEQkKCTpsHDx4gLCwM5cuXh7m5Obp06YLU1FSl/tatWwgICIC5uTnq16+PuLg4ne3DwsIwc+bMEjketZs6dSo0Gg1GjBihlHF+St/ly5cREhKC8uXLw8TEBJ6enjh48KBSLyL4+uuvYW9vDxMTE/j6+iIxMVGpf/jwIXr16gVLS0u4ublhy5YtOv3PmDEDH3/8cYkdj5rk5OTgq6++gouLC0xMTFC9enVMmDBB57uVOD8lZ+fOnQgICICDgwM0Gg1+//13nfrnzQXw+Hdaz549YWlpCWtra/Tv3x93795V6pOTk9GyZUuYmZmhZcuWSE5O1tm+Y8eOWLNmzYsdgBTRypUrxdDQUBYvXiwnTpyQjz76SKytrSU1NVXWrVsnlSpVkgMHDsjy5cvF2NhYrl+/LiIid+7cEVdXVzl//nxRd0lP4efnJ0uWLJHjx4/L4cOHpUOHDlK1alW5e/eu0mbQoEHi6OgoW7dulYMHD0rjxo2ladOmSv3IkSOlVatWkpCQICNGjBBvb2+lLiYmRry9vSU7O7tEj0uN9u/fL87OzuLl5SXDhw9Xyjk/pevWrVvi5OQkoaGhsm/fPjl37pxs2rRJzpw5o7SZOnWqWFlZye+//y5HjhyR9957T1xcXOT+/fsiIvLdd9+Jh4eHHD9+XGbMmCG2traSm5srIiLnzp0TV1dXSUtLK5Xje91NmjRJypcvL1FRUZKUlCSrV68Wc3NzmTNnjtKG81NyNmzYIF988YVERkYKAFm7dq1O/fPmQkSkXbt2UrduXdm7d6/s2rVLatSoIcHBwUp9586dpXv37nL69Gnp2rWrdOnSRalbuXKlBAQEvPD4ixx4GjZsKGFhYcr7nJwccXBwkClTpsi0adOkW7duSl3FihVl//79IiIyYMAAmTVr1gsPlJ7v2rVrAkD++usvEXkcMg0MDGT16tVKm/j4eAEgMTExIiLSvn17WbBggYiInDx5UkxNTUVEJCsrS+rWrSsHDhwo4aNQn4yMDHF1dZXo6Ghp1aqVEng4P6VvzJgx0rx586fW5+bmip2dncyYMUMpu3PnjhgZGcmKFStERGTw4MEyZswYERHJzMwUAHLt2jURefw/JZGRka/wCNTN399f+vXrp1PWuXNn6dmzp4hwfkrTPwNPYebi5MmTAkDn99bGjRtFo9HI5cuXRUTEw8NDNm7cKCKPA1atWrVEROT27dtSo0YNuXDhwguPuUiXtLKyshAbGwtfX1+lTE9PD76+voiJiUHdunVx8OBB3L59G7Gxsbh//z5q1KiB3bt349ChQxg2bNiLnYaiQklLSwMA2NjYAABiY2Px6NEjnflyd3dH1apVERMTAwCoW7cutm3bhuzsbGzatAleXl4AgOnTp6N169bw8fEp4aNQn7CwMPj7++vMA8D5KQvWrVsHHx8fBAUFoWLFiqhfvz5++uknpT4pKQlXr17VmSMrKys0atRIZ452796N+/fvY9OmTbC3t0eFChUQEREBY2NjvP/++yV+XGrRtGlTbN26FadPnwYAHDlyBLt370b79u0BcH7KksLMRUxMDKytrXV+b/n6+kJPTw/79u0D8Hi+tmzZgtzcXGzevFn5nTd69GiEhYXB0dHxxQdZlHR0+fJlASB///23Tvno0aOlYcOGIiIybtw4qV69utSpU0ciIyPl4cOHUqdOHTl48KDMnTtX3NzcpGnTpnL8+PEXTmmUX05Ojvj7+0uzZs2UsoiICDE0NMzX9q233pLPPvtMRB4n8ODgYKlataq0bNlSTpw4IadPnxZXV1e5ceOGDBw4UFxcXCQoKEju3LlTYsejFitWrJA6deoop3SfPMPD+Sl9RkZGYmRkJGPHjpVDhw7JDz/8IMbGxvLzzz+LiMiePXsEgFy5ckVnu6CgIOnatauIPD7bNmTIEHF2dhYfHx/ZtWuX3Lx5U6pVqyYXLlyQL774QqpXry5t27aVS5culfgxvs5ycnJkzJgxotFoRF9fXzQajUyePFmp5/yUHvzjDE9h5mLSpEni5uaWry9bW1uZP3++iIhcunRJ/P39xdHRUfz9/eXSpUvy119/iY+Pj9y8eVOCgoLExcVFBg4cKA8fPizSmAv9bemFFR4ejvDwcOX9N998A19fXxgYGGDixIk4duwYoqKi0Lt3b8TGxhb37t9YYWFhOH78OHbv3l2k7aysrLB8+XKdsnfeeQczZsxAREQEzp07h4SEBHz00UcYP348F8gWwcWLFzF8+HBER0fD2Nj4hfrg/Lxaubm58PHxweTJkwEA9evXx/Hjx7Fw4UL06dOnUH0YGBhg3rx5OmV9+/bFsGHDEBcXh99//x1HjhzB9OnTMWzYsBdfcPkG+vXXXxEREYHly5ejdu3aOHz4MEaMGAEHBwfOj0pVrlwZUVFRyvuHDx/Cz88PS5cuxcSJE2FhYYGEhAS0a9cOP/zwQ5EWnBfpklaFChWg1Wp17iIBgNTUVNjZ2eVrf+rUKSxbtgwTJkzAjh070LJlS9ja2qJr1644dOgQMjIyirJ7eoqhQ4ciKioK27dvR5UqVZRyOzs7ZGVl4c6dOzrtnzZfALBkyRJYW1ujU6dO2LFjBwIDA2FgYICgoCDs2LHjFR6F+sTGxuLatWto0KAB9PX1oa+vj7/++gvfffcd9PX1UalSJc5PKbO3t0etWrV0yjw8PHDhwgUAUOahsL/zAGD79u04ceIEhg4dih07dqBDhw4wMzND165dOUdFNHr0aHz++efo3r07PD090atXL3zyySeYMmUKAM5PWVKYubCzs8O1a9d06rOzs3Hr1q2nztfkyZPRtm1beHt7Y8eOHejSpQsMDAzQuXPnIs9XkQKPoaEhvL29sXXrVqUsNzcXW7duRZMmTXTaiggGDhyIWbNmwdzcHDk5OXj06BEAKP/Nyckp0mBJl4hg6NChWLt2LbZt2wYXFxedem9vbxgYGOjMV0JCAi5cuJBvvgDg+vXrGD9+PObOnQsA+eaM81U0bdq0wbFjx3D48GHl5ePjg549eyp/5vyUrmbNmuV7lMPp06fh5OQEAHBxcYGdnZ3OHKWnp2Pfvn0FzlHeYwZ++OEHaLVaztFLyszMhJ6e7j9TWq0Wubm5ADg/ZUlh5qJJkya4c+eOztWdbdu2ITc3F40aNcrXZ3x8PJYvX44JEyYAKIbfeUW6ACaPbwszMjKSn3/+WU6ePCkDBgwQa2truXr1qk67H3/8Ued2sn379omlpaXExMTI119/ray8phc3ePBgsbKykh07dkhKSoryyszMVNoMGjRIqlatKtu2bZODBw9KkyZNpEmTJgX216NHD5k7d67yftq0aeLt7S0nT56U9u3by5AhQ175Mandk2t4RDg/pW3//v2ir68vkyZNksTERImIiBBTU1NZtmyZ0mbq1KlibW0t//vf/+To0aPSqVOnfLfa5vnXv/4lo0aNUt6vWrVKqlatKkeOHJH+/ftLhw4dSuS41KJPnz5SuXJl5bb0yMhIqVChgrLGTYTzU5IyMjIkLi5O4uLiBIDMmjVL4uLilMfNFGYu2rVrJ/Xr15d9+/bJ7t27xdXVVee29Dy5ubnSvHlzWb9+vVI2ePBg8ff3l5MnT0r9+vVl+vTpRRp/kQOPiMjcuXOlatWqYmhoKA0bNpS9e/fq1F+9elWcnJyU28zyfPPNN2JjYyPu7u6yb9++F9k1PQFAga8lS5Yobe7fvy9DhgyRcuXKiampqbz//vuSkpKSr68///xTGjZsKDk5OUrZvXv3JCgoSCwsLKRNmzaSmppaEoelav8MPJyf0rd+/XqpU6eOGBkZibu7u/z444869bm5ufLVV19JpUqVxMjISNq0aSMJCQn5+jl27JjUqFFD5zlYOTk5MnjwYLG0tJS33npLEhMTX/nxqEl6eroMHz5cqlatKsbGxlKtWjX54osvdBarcn5Kzvbt2wv8N6dPnz4iUri5uHnzpgQHB4u5ublYWlpK3759JSMjI9++Fi5cqHPSREQkNTVV2rRpIxYWFhIUFCT37t0r0vg1Ik88spKIiIhIhfhdWkRERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8REREpHoMPERERKR6DDxERESkegw8RPTaCA8Ph0ajgUajwb///e+X6qt169ZKX4cPHy6W8RFR2cXAQ6QSMTEx0Gq18Pf3z1e3Y8cOaDSafN/MDgDOzs464SEvBGg0GlhZWaFZs2bYtm2bUh8aGorAwECd9xqNBoMGDcrXd1hYGDQaDUJDQ3XKL168iH79+sHBwQGGhoZwcnLC8OHDcfPmzeceZ+3atZGSkoIBAwYoZSNHjoSNjQ0cHR0RERGh03716tUICAjI109kZCT279//3P0RkTow8BCpxKJFi/Dxxx9j586duHLlykv1tWTJEqSkpGDPnj2oUKECOnbsiHPnzj21vaOjI1auXIn79+8rZQ8ePMDy5ctRtWpVnbbnzp2Dj48PEhMTsWLFCpw5cwYLFy7E1q1b0aRJE9y6deuZY9PX14ednR1MTU0BAOvXr8fy5cuxefNmTJ8+HR9++CFu3LgBAEhLS8MXX3yBefPm5evHxsYGtra2hf5MiOj1xsBDpAJ3797FqlWrMHjwYPj7++Pnn39+qf6sra1hZ2eHOnXqYMGCBbh//z6io6Of2r5BgwZwdHREZGSkUhYZGYmqVauifv36Om3DwsJgaGiIzZs3o1WrVqhatSrat2+PLVu24PLly/jiiy+KNNb4+Hi0bt0aPj4+CA4OhqWlJZKSkgAAn332GQYPHpwvdBHRm4eBh0gFfv31V7i7u6NmzZoICQnB4sWLUVzfC2xiYgIAyMrKema7fv36YcmSJcr7xYsXo2/fvjptbt26hU2bNmHIkCFKv3ns7OzQs2dPrFq1qkhjr1u3Lg4ePIjbt28jNjYW9+/fR40aNbB7924cOnQIw4YNK3RfRKReDDxEKrBo0SKEhIQAANq1a4e0tDT89ddfL91vZmYmvvzyS2i1WrRq1eqZbUNCQrB7926cP38e58+fx549e5Qx5UlMTISIwMPDo8A+PDw8cPv2bVy/fr3QY/Tz80NISAjeeusthIaGYunSpTAzM8PgwYOxcOFCLFiwADVr1kSzZs1w4sSJQvdLROqiX9oDIKKXk5CQgP3792Pt2rUAHq9x6datGxYtWoTWrVu/UJ/BwcHQarW4f/8+bG1tsWjRInh5eT1zG1tbW+VymojA398fFSpUKLBtcZ19yhMeHo7w8HDl/TfffANfX18YGBhg4sSJOHbsGKKiotC7d2/ExsYW676J6PXAwEP0mlu0aBGys7Ph4OCglIkIjIyM8P3338PKygqWlpYAHi/itba21tn+zp07sLKy0imbPXs2fH19YWVlVaSFvf369cPQoUMBoMCFwjVq1IBGo0F8fDzef//9fPXx8fEoV67cSy0mPnXqFJYtW4a4uDgsXrwYLVu2hK2tLbp27Yp+/fohIyMDFhYWL9w/Eb2eeEmL6DWWnZ2N//73v5g5cyYOHz6svI4cOQIHBwesWLECAODq6go9Pb18ZzfOnTuHtLQ0uLm56ZTb2dmhRo0aRQ4e7dq1Q1ZWFh49egQ/P7989eXLl8e7776L+fPn69zRBQBXr15FREQEunXrBo1GU6T95hERDBw4ELNmzYK5uTlycnLw6NEjAFD+m5OT80J9E9HrjWd4iF5jUVFRuH37Nvr375/vLE2XLl2waNEiDBo0CBYWFvjwww8xatQo6Ovrw9PTExcvXsSYMWPQuHFjNG3atFjGo9VqER8fr/y5IN9//z2aNm0KPz8/TJw4ES4uLjhx4gRGjx6NypUrY9KkSS+8///85z+wtbVVnrvTrFkzhIeHY+/evdi4cSNq1aqV7wwXEb0ZeIaH6DW2aNEi5dLTP3Xp0gUHDx7E0aNHAQBz5sxBnz59MGbMGNSuXRuhoaHw8vLC+vXrX/iMSkEsLS2VS2gFcXV1xcGDB1GtWjV07doV1atXx4ABA/D2228jJiYGNjY2L7Tf1NRUTJo0Cd99951S1rBhQ4waNQr+/v749ddfde4iI6I3i0aKe/UgEdErEh4ejt9//73YvgoiOTkZLi4uiIuLQ7169YqlTyIqm3iGh4heK8eOHYO5uTnmz5//Uv20b98etWvXLqZREVFZxzM8RPTauHXrlvLVE7a2tgVeyiusy5cvKwunq1atCkNDw2IZIxGVTQw8REREpHq8pEVERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqsfAQ0RERKrHwENERESqx8BDREREqvf/AOgpDv2lH617AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 700x200 with 1 Axes>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# fliers\n",
    "#    if the distribution is very skewed to one side,\n",
    "#    it is possible that some extreme values are considered\n",
    "#    are considered as outliers, showing as fliers in the boxplot\n",
    "#\n",
    "#    there are two types of fliers: high and low\n",
    "#    they are defined as:\n",
    "#        - high: values > high whisker = Q3 + 1.5 * IQR\n",
    "#        - low: values < low whisker = Q1 - 1.5 * IQR\n",
    "#     where IQR = Q3 - Q1\n",
    "\n",
    "# let's artificially simulate this situation\n",
    "#     we will create a distortion in the values so that\n",
    "#     high values (close to 1) become even higher\n",
    "#     and low values (close to 0) become even lower\n",
    "\n",
    "\n",
    "def distortion(vals: Tensor) -> Tensor:\n",
    "    \"\"\"Artificial distortion to simulate a skewed distribution.\n",
    "\n",
    "    To visualize it:\n",
    "    ```\n",
    "    fig, ax = plt.subplots()\n",
    "    t = np.linspace(0, 1, 100)\n",
    "    ax.plot(t, np.clip(distortion(t), 0, 1), label=\"distortion\")\n",
    "    ax.plot(t, t, label=\"identity\", linestyle=\"--\")\n",
    "    fig\n",
    "    ```\n",
    "    \"\"\"\n",
    "    return vals + 0.12 * (vals * (1 - vals) * 4)\n",
    "\n",
    "\n",
    "mock = torch.clip(distortion(aupimo_result.aupimos.clone()), 0, 1)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 2))\n",
    "ax.boxplot(\n",
    "    mock[labels == 1].numpy(),\n",
    "    vert=False,\n",
    "    widths=0.4,\n",
    ")\n",
    "_ = ax.set_yticks([])\n",
    "ax.set_xlim(0 - (eps := 2e-2), 1 + eps)\n",
    "ax.xaxis.set_major_formatter(PercentFormatter(1))\n",
    "ax.set_xlabel(\"AUPIMO [%]\")\n",
    "ax.set_title(\"AUPIMO Scores Boxplot\")\n",
    "num_images = (labels == 1).sum().item()\n",
    "ax.annotate(\n",
    "    text=f\"Number of images: {num_images}\",\n",
    "    xy=(0.03, 0.95),\n",
    "    xycoords=\"axes fraction\",\n",
    "    xytext=(0, 0),\n",
    "    textcoords=\"offset points\",\n",
    "    annotation_clip=False,\n",
    "    verticalalignment=\"top\",\n",
    ")\n",
    "fig  # noqa: B018, RUF100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.24     0.24     44      1\n",
      "1        q1   0.65     0.65     58      1\n",
      "2      mean   0.79     0.78     29      1\n",
      "3       med   0.94     0.93     63      1\n",
      "4        q3   1.00     1.00     22      1\n",
      "5    whishi   1.00     1.00      0      1\n"
     ]
    }
   ],
   "source": [
    "# `None` is the default policy, so the fliers are not returned\n",
    "print(pd.DataFrame.from_records(boxplot_stats(mock, labels, flier_policy=None)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "with option 'low'\n",
      "  statistic  value  nearest  index  label\n",
      "0   flierlo   0.00     0.00     65      1\n",
      "1   flierlo   0.00     0.00     67      1\n",
      "2   flierlo   0.01     0.01     71      1\n",
      "3   flierlo   0.09     0.09     64      1\n",
      "4    whislo   0.24     0.24     44      1\n",
      "5        q1   0.65     0.65     58      1\n",
      "6      mean   0.79     0.78     29      1\n",
      "7       med   0.94     0.93     63      1\n",
      "8        q3   1.00     1.00     22      1\n",
      "9    whishi   1.00     1.00      0      1\n",
      "with option 'both'\n",
      "  statistic  value  nearest  index  label\n",
      "0   flierlo   0.00     0.00     65      1\n",
      "1   flierlo   0.00     0.00     67      1\n",
      "2   flierlo   0.01     0.01     71      1\n",
      "3   flierlo   0.09     0.09     64      1\n",
      "4    whislo   0.24     0.24     44      1\n",
      "5        q1   0.65     0.65     58      1\n",
      "6      mean   0.79     0.78     29      1\n",
      "7       med   0.94     0.93     63      1\n",
      "8        q3   1.00     1.00     22      1\n",
      "9    whishi   1.00     1.00      0      1\n"
     ]
    }
   ],
   "source": [
    "# one can choose to include only high or low fliers, or both\n",
    "# since there are only low fliers...\n",
    "\n",
    "# 'low' and 'both' will return the same result\n",
    "print(\"with option 'low'\")\n",
    "print(pd.DataFrame.from_records(boxplot_stats(mock, labels, flier_policy=\"low\")))\n",
    "\n",
    "print(\"with option 'both'\")\n",
    "print(pd.DataFrame.from_records(boxplot_stats(mock, labels, flier_policy=\"both\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "with option 'high'\n",
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.24     0.24     44      1\n",
      "1        q1   0.65     0.65     58      1\n",
      "2      mean   0.79     0.78     29      1\n",
      "3       med   0.94     0.93     63      1\n",
      "4        q3   1.00     1.00     22      1\n",
      "5    whishi   1.00     1.00      0      1\n"
     ]
    }
   ],
   "source": [
    "# and 'high' will return no fliers (same as `flier_policy=None` in this case)\n",
    "print(\"with option 'high'\")\n",
    "print(pd.DataFrame.from_records(boxplot_stats(mock, labels, flier_policy=\"high\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other applications and `only_label` argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "stats for the maximum anomaly score in the anomaly maps\n",
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.46     0.46     65      1\n",
      "1        q1   0.63     0.63     48      1\n",
      "2       med   0.70     0.71     10      1\n",
      "3      mean   0.73     0.73    118      1\n",
      "4        q3   0.81     0.81    115      1\n",
      "5    whishi   1.00     1.00     22      1\n"
     ]
    }
   ],
   "source": [
    "# other applications\n",
    "#    since the function is agnostic to the meaning of the values\n",
    "#    we can also use it to find representative samples\n",
    "#    with another metric or signal\n",
    "#\n",
    "#    in the last plot cell we used the maximum anomaly score per image\n",
    "#    to select normal images, so let's reuse that criterion here\n",
    "\n",
    "# recompute it for didactic purposes\n",
    "max_anom_score_per_image = anomaly_maps.max(dim=2).values.max(dim=1).values  # noqa: PD011\n",
    "print(\"stats for the maximum anomaly score in the anomaly maps\")\n",
    "print(pd.DataFrame.from_records(boxplot_stats(max_anom_score_per_image, labels)))\n",
    "# notice that the indices are not the same as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.42     0.42     90      0\n",
      "1        q1   0.43     0.43     80      0\n",
      "2       med   0.45     0.45    105      0\n",
      "3      mean   0.46     0.46     89      0\n",
      "4        q3   0.48     0.48     75      0\n",
      "5    whishi   0.52     0.52     95      0\n"
     ]
    }
   ],
   "source": [
    "# we can also use the `only_label` argument to select only the\n",
    "# samples from the normal class\n",
    "print(pd.DataFrame.from_records(boxplot_stats(max_anom_score_per_image, labels, only_label=0)))\n",
    "# notice the labels are all 0 now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  statistic  value  nearest  index  label\n",
      "0    whislo   0.42     0.42     90      0\n",
      "1        q1   0.52     0.52     95      0\n",
      "2       med   0.65     0.65     17      1\n",
      "3      mean   0.66     0.66     45      1\n",
      "4        q3   0.77     0.77    108      1\n",
      "5    whishi   1.00     1.00     22      1\n"
     ]
    }
   ],
   "source": [
    "# or we can consider data from both classes (`None` option)\n",
    "print(pd.DataFrame.from_records(boxplot_stats(max_anom_score_per_image, labels, only_label=None)))\n",
    "# notice that the labels are mixed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cite Us\n",
    "\n",
    "AUPIMO was developed during [Google Summer of Code 2023 (GSoC 2023)](https://summerofcode.withgoogle.com/archive/2023/projects/SPMopugd) with the `anomalib` team from Intel's OpenVINO Toolkit.\n",
    "\n",
    "arXiv: [arxiv.org/abs/2401.01984](https://arxiv.org/abs/2401.01984) (accepted to BMVC 2024)\n",
    "\n",
    "Official repository: [github.com/jpcbertoldo/aupimo](https://github.com/jpcbertoldo/aupimo) (numpy-only API and numba-accelerated versions available)\n",
    "\n",
    "```bibtex\n",
    "@misc{bertoldo2024aupimo,\n",
    "      author={Joao P. C. Bertoldo and Dick Ameln and Ashwin Vaidya and Samet Akçay},\n",
    "      title={{AUPIMO: Redefining Visual Anomaly Detection Benchmarks with High Speed and Low Tolerance}}, \n",
    "      year={2024},\n",
    "      url={https://arxiv.org/abs/2401.01984}, \n",
    "}\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "anomalib-dev",
   "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.14"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}