doc/source/examples/cd_ks_cifar10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kolmogorov-Smirnov data drift detector on CIFAR-10\n",
    "\n",
    "### Method\n",
    "\n",
    "The drift detector applies feature-wise two-sample [Kolmogorov-Smirnov](https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test) (K-S) tests. For multivariate data, the obtained p-values for each feature are aggregated either via the [Bonferroni](https://mathworld.wolfram.com/BonferroniCorrection.html) or the [False Discovery Rate](http://www.math.tau.ac.il/~ybenja/MyPapers/benjamini_hochberg1995.pdf) (FDR) correction. The Bonferroni correction is more conservative and controls for the probability of at least one false positive. The FDR correction on the other hand allows for an expected fraction of false positives to occur.\n",
    "\n",
    "For high-dimensional data, we typically want to reduce the dimensionality before computing the feature-wise univariate K-S tests and aggregating those via the chosen correction method. Following suggestions in [Failing Loudly: An Empirical Study of Methods for Detecting Dataset Shift](https://arxiv.org/abs/1810.11953), we incorporate Untrained AutoEncoders (UAE) and black-box shift detection using the classifier's softmax outputs ([BBSDs](https://arxiv.org/abs/1802.03916)) as out-of-the box preprocessing methods and note that [PCA](https://en.wikipedia.org/wiki/Principal_component_analysis) can also be easily implemented using `scikit-learn`. Preprocessing methods which do not rely on the classifier will usually pick up drift in the input data, while BBSDs focuses on label shift. The [adversarial detector](https://arxiv.org/abs/2002.09364) which is part of the library can also be transformed into a drift detector picking up drift that reduces the performance of the classification model. We can therefore combine different preprocessing techniques to figure out if there is drift which hurts the model performance, and whether this drift can be classified as input drift or label shift.\n",
    "\n",
    "### Backend\n",
    "\n",
    "The method works with both the **PyTorch** and **TensorFlow** frameworks for the optional preprocessing step. Alibi Detect does however not install PyTorch for you. \n",
    "Check the [PyTorch docs](https://pytorch.org/) how to do this.\n",
    "\n",
    "\n",
    "### Dataset\n",
    "\n",
    "[CIFAR10](https://www.cs.toronto.edu/~kriz/cifar.html) consists of 60,000 32 by 32 RGB images equally distributed over 10 classes. We evaluate the drift detector on the CIFAR-10-C dataset ([Hendrycks & Dietterich, 2019](https://arxiv.org/abs/1903.12261)). The instances in\n",
    "CIFAR-10-C have been corrupted and perturbed by various types of noise, blur, brightness etc. at different levels of severity, leading to a gradual decline in the classification model performance. We also check for drift against the original test set with class imbalances. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import tensorflow as tf\n",
    "\n",
    "from alibi_detect.cd import KSDrift\n",
    "from alibi_detect.models.tensorflow import scale_by_instance\n",
    "from alibi_detect.utils.fetching import fetch_tf_model, fetch_detector\n",
    "from alibi_detect.saving import save_detector, load_detector\n",
    "from alibi_detect.datasets import fetch_cifar10c, corruption_types_cifar10c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data\n",
    "\n",
    "Original CIFAR-10 data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "(X_train, y_train), (X_test, y_test) = tf.keras.datasets.cifar10.load_data()\n",
    "X_train = X_train.astype('float32') / 255\n",
    "X_test = X_test.astype('float32') / 255\n",
    "y_train = y_train.astype('int64').reshape(-1,)\n",
    "y_test = y_test.astype('int64').reshape(-1,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For CIFAR-10-C, we can select from the following corruption types at 5 severity levels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['brightness', 'contrast', 'defocus_blur', 'elastic_transform', 'fog', 'frost', 'gaussian_blur', 'gaussian_noise', 'glass_blur', 'impulse_noise', 'jpeg_compression', 'motion_blur', 'pixelate', 'saturate', 'shot_noise', 'snow', 'spatter', 'speckle_noise', 'zoom_blur']\n"
     ]
    }
   ],
   "source": [
    "corruptions = corruption_types_cifar10c()\n",
    "print(corruptions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's pick a subset of the corruptions at corruption level 5. Each corruption type consists of perturbations on all of the original test set images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "corruption = ['gaussian_noise', 'motion_blur', 'brightness', 'pixelate']\n",
    "X_corr, y_corr = fetch_cifar10c(corruption=corruption, severity=5, return_X_y=True)\n",
    "X_corr = X_corr.astype('float32') / 255"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split the original test set in a reference dataset and a dataset which should not be rejected under the *H<sub>0</sub>* of the K-S test. We also split the corrupted data by corruption type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5000, 32, 32, 3) (5000, 32, 32, 3)\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "n_test = X_test.shape[0]\n",
    "idx = np.random.choice(n_test, size=n_test // 2, replace=False)\n",
    "idx_h0 = np.delete(np.arange(n_test), idx, axis=0)\n",
    "X_ref,y_ref = X_test[idx], y_test[idx]\n",
    "X_h0, y_h0 = X_test[idx_h0], y_test[idx_h0]\n",
    "print(X_ref.shape, X_h0.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Class Ref H0\n",
      "0     472 528\n",
      "1     510 490\n",
      "2     498 502\n",
      "3     492 508\n",
      "4     501 499\n",
      "5     495 505\n",
      "6     493 507\n",
      "7     501 499\n",
      "8     516 484\n",
      "9     522 478\n"
     ]
    }
   ],
   "source": [
    "# check that the classes are more or less balanced\n",
    "classes, counts_ref = np.unique(y_ref, return_counts=True)\n",
    "counts_h0 = np.unique(y_h0, return_counts=True)[1]\n",
    "print('Class Ref H0')\n",
    "for cl, cref, ch0 in zip(classes, counts_ref, counts_h0):\n",
    "    assert cref + ch0 == n_test // 10\n",
    "    print('{}     {} {}'.format(cl, cref, ch0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_corr = len(corruption)\n",
    "X_c = [X_corr[i * n_test:(i + 1) * n_test] for i in range(n_corr)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualise the same instance for each corruption type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVUElEQVR4nO2dWYxl11WG/33One+tW7fmnmx3t7vjMbYcxZkYFQMyEiFSggIyUx54AKQ8oMADQUTJIwjxAgIhEA9gJEAWREqEopDwgIzMkBA7sYkdW3bP1UNVdVV11Z3P3TxUCbVg/6vdRbtrdef/pEjRWb3PPefc899d3v9ea4UYI4QQ/sj2+wKEEGkkTiGcInEK4RSJUwinSJxCOEXiFMIpEucdRAjhMyGEP7vV//ZtnCuGEE7cinOJt0+Qz7l/hBA+CeDTAO4HsAng7wH8ZoxxfR8v6/8QQogATsYY39jva/leQjPnPhFC+DSA3wHwGwCmAXwAwH0A/jGEUEn8+9LtvUKx30ic+0AIoQ3g8wA+FWP8coxxFGM8BeATAI4C+LkQwudCCM+FEJ4NIWwC+OTusWevO88vhBBOhxBWQwi/HUI4FUL4kd3Y//zbEMLR3T9NfzGEcCaEsBJC+K3rzvO+EMILIYT1EMJyCOEPUz8Q4vYice4PHwJQA/B31x+MMW4B+AcAP7p76KMAngPQAfBX1//bEMLDAP4IwM8COIid2ffwDT73+wE8AOApAJ8NITy0e7wA8GsA5gF8cDf+qzd/W+JWInHuD/MAVmKM40RseTcOAC/EGL8QY5zEGHv/69/9FIAvxhifjzEOAXwWwI0WED4fY+zFGF8C8BKAxwEgxviNGOO/xhjHuzP4nwD4ob3dmrhV6L9j9ocVAPMhhFJCoAd34wBw1jjHoevjMcZuCGH1Bp978br/3wXQAoAQwrsA/D6A9wJoYOe9+MaNbkK8s2jm3B9eADAA8LHrD4YQWgB+HMDXdg9ZM+EygCPXja0DmNvj9fwxgFexsyLbBvAZAGGP5xK3CIlzH4gxbmBnQegPQghPhxDKIYSjAP4WwDkAf/k2TvMcgI+EED60u3jzOexdUFPYsXK2QggPAviVPZ5H3EIkzn0ixvi72Jmhfg87wvg37PyZ+lSMcfA2xr8C4FMA/ho7s+gWgMvYmZFvll8H8AyAawD+FMDf7OEc4hajTQh3Cbt/Eq9j50/Tt/b5csQtQDPnHUwI4SMhhEYIoYmdGfjbAE7t71WJW4XEeWfzUQAXdv93EsDPRP0pdNegP2uFcIpmTiGcYm5CWFlZodPqeJza3LJDCHenRebmvqw/doyYOYz8TEdjVMYG3ejDwoSHSCwaLlEw5ph34i/DvbwH1nUsLS0lT6iZUwinSJxCOEXiFMIpEqcQTpE4hXCKxCmEU0wrJc/z23UddwRurBSDMClozDQVsvS9TaxEl2i8H9GwPjJ+JQHMZrGu/s62UhiaOYVwisQphFMkTiGcInEK4RSJUwinSJxCOMW0Uqzl3+/FPNDbec/mcr11HZFnfJiuCLVF+O/3YMQzk0rlMv+wgl9jHvbyjI17doKsFCHuIiROIZwicQrhFIlTCKdInEI4xVyttVYM74RN4Iw7fqXZePSFtcI+4QPHk/SK52jMN9K//uabNLZ0YJHGJsMhjS3MziSP16p89XdyB3yfe9GLZk4hnCJxCuEUiVMIp0icQjhF4hTCKRKnEE55Rza+38k2i8Ve7+vWWzf8OvJyhcYKo65Pbyvdc3d9Y5uOubSyRmP1qSaNzU1N0VgW0vOF1XKBtXD4f2HZiLf+05Jo5hTCKRKnEE6ROIVwisQphFMkTiGcInEK4RTTSslIiX7AznC4nRjuwA36D6Sx7JJsj1ZKYSy+T0g2SJ7z383hcERjV1Y3aWxzu09jvUE6+2S7m7ZYACCrNmhsu8czT1oN/sWMSYgbRKbr8Y5wu6xCzZxCOEXiFMIpEqcQTpE4hXCKxCmEUyROIZxiWinb3R4PTvhyeIl0xI7GmLzEuyRbsWCU72c2SzbZ229SZuUjGMvrWwNuYbCMlXqJfzV9ow3CsmGlXL7KY6yD9Yh5GwC617b4ZxkZK+fOL9PYwyePJ4/ff/QIHZNHo5u32brCeA8st4TErE4S5rtDxwghXCJxCuEUiVMIp0icQjhF4hTCKRKnEE4xrZT1Hs9IaDV4AaeslO5rUUy4BWC6G8YqdG7EMuKlhGyPv0l7LGp2cfk8jc3OziaP12s8D2PQ79JYo8rHHViYp7FIHvJ2l9tAzQr/rGGf23B5xgtybQ3S79zY7NvDX2O7uJp1zj2M2mMzcoZmTiGcInEK4RSJUwinSJxCOEXiFMIp5mptqT1HY4Wx4jnKyEb1wDcoW7FiwmOZtYJKYnEvxYVg1ysyyi1hPOSr3oFt2jZWtjtGq4PRyLi3nHeHbrTSLRKs1dqQV40YfyDVOr+OQB7kmLRpAIBodWPY43dmFaBiV2+f7ubfOc2cQjhF4hTCKRKnEE6ROIVwisQphFMkTiGcYlopf/4Xz9JYMOoBlcnG99ZUjY45cexeGnvysYdprGT8vLCaRWbHbmt93dgNPTasjxmyuR0AKtX0M2Eb0QGgUuEWxtwMr7cUwWMlsom9YtQyQpl/n/0xfx7rm1d5bGMjefzaxjodM7JqXRmFfebmOjR28kS6lhEAlCvpZ2K5JcwistDMKYRTJE4hnCJxCuEUiVMIp0icQjhF4hTCKaaV0jMyEoY9HiuT5fdr6VVyAEDDWLIvHnqQxvqRd1DOiJVSrdTpGGs5vLAsGMNmmZ5doDHaLdvI+hmSbtgAkBt1fWBkdrAzTozsjFOn36Sx85cv09ja6iqN9XppW6QYcGtmaHTRHgx4vaUj9yzR2L338PYPTWKlWJksljXG0MwphFMkTiGcInEK4RSJUwinSJxCOEXiFMIpppXyiY99nMYGRiZAs562KoKx1Fyny9NAMAo4bW4a3ZrHo+TxcolnU5TqPBaNDtu9EV/OjxN+bxmxTFhmDwCUjOsol40WA9nNW0Ejwz7qT9LPFwCa7RaNzXQ6NFYM0+es5dz+Wl/lHt2586do7MSxEzSWZ4a1R55JbthpascgxF2ExCmEUyROIZwicQrhFIlTCKdInEI4xbRSJiMj+8HQNVvob1V4j496jRet6vW5XdId8T4qp948lTxeMbJS7j12H429dfYCjX3py1+jsVHGbZEa6UTdMJ5H07B7ptttGutMp/uhAMATTzyWPL4wP0PH3H/kMI1lgds9uZEdM+yn+8qUDGujt8gLqB062OGxwwdprCj4e9Xtpu0eZiECZkIQRTOnEE6ROIVwisQphFMkTiGcInEK4RSJUwinmFbKF774FRqbjHhGQoZ0hkar0qBjpgwL4OhJXmxpYY5nP8wdTPdfmZ1fpGNqTW5TrH/nNI29/J2zNNYzUhJYgknJyOCZMq7xxL3cCvrg+95DY3PNtM3SzPkrEo2aVcMhL8g1LtJ2CQB0SU+UUcHft3qDP49Oh9t3ly5eorGVlTX+ec20ZbJ0gL9XjQa3xubb6WevmVMIp0icQjhF4hTCKRKnEE6ROIVwirla+/VvvkxjtTIv+z8cpDeqlyv8t+D9H3iSxk6f5yuhq8s0hEcfeSR5vGJsHO8OeC2gsrEZ/Yn3pDeOA0C/x1cnK+X0V3Dy+DE65pGHHqCxQ/MdGms3+MbsST9932cvXqFjLl/lHaqXV/i47a1tGltfX08eH474M2SdpgHeORwAijFfER+N+Gpzo5NeXX0U6fcNAKaNpIPjB9LtOjRzCuEUiVMIp0icQjhF4hTCKRKnEE6ROIVwimmlXDnHN3rPzvDaMoePpDcAP/zYSTqmXOW7qF958d9pbKnGl8pbIV0H5vIK91+a7Wkam2vzz/rJp3+QxjKjgMz0dPrz5ufm6Ji1Nd4Z+q3Tr9PYxjqvxbS5cS15/Nom7wy9vs0tkbVN3iJhbCRNlMvpekuVKq/DlOXG823z96pjtIWYWeTWR7WRTuCo1Hlix5bRCZ6hmVMIp0icQjhF4hTCKRKnEE6ROIVwisQphFNMK+X8d/+LxjaNzsU/8WO/nDz+9NNP0TFf/Sder2iRZAEAwGLDaPFQSi+j14xW2UvTvJbRlBGrGXVsxkY9IJY1MS74NV587TyNnbnM6+IMR0Yto1r6OU5N8VYHizVuHYxIh+obUa6kLZPcsEus2NQUf3fapHbPzjm5BbO1nbaXLl1aoWP6fW5J4b2PJw9r5hTCKRKnEE6ROIVwisQphFMkTiGcInEK4RTTSul3edbBux9/lMY+/NSHk8fnOjzT4vveb2R1ZEZrgjIvutVupe2BvMJtj5LR9Toa1zEhLSgAYOMqzyJpl9LXP6H9wYHjD/Bnv3jkXTS2dpVnpUyRDI1Rwe85RP7bXs749U8m3Cbq99PZG1vbW3RMnPAu1FtdPu7sMs9O6ve49THqpq/R6obdaPL3lKGZUwinSJxCOEXiFMIpEqcQTpE4hXCKxCmEU0wr5fiD6d3yAPDTP/9LNNYt0pkFr73BMyYmgRdwqhkZMCOjvfLaOlnanvBl8qLo0VgwntYEvJfHtc108SwAyC+lszcuXL5MxwwGPONj0uc9PppGBs+br59LHn/rzBk6JpT4dzY7z22z4YA/q42NdGGw1RWe8RENCyPLuG0TjFizzi21DsngqRm9dHpb/L1iaOYUwikSpxBOkTiFcIrEKYRTJE4hnGKu1n78mWdobObAERp76eX0yt/QqCszNDZDF8Ym8DgxassgvZIbjJo+hVG7JxrjMvNnzuigPE5/3soqX9kej/nKn7EAiU67Q2PDYXoFdW2VJz8g59/LygpvPzAY8esfk7YFxZAnFuRGZ+tGjXdgr1p1icb83oZ99h7zVeN6kydbMDRzCuEUiVMIp0icQjhF4hTCKRKnEE6ROIVwimmlfPPFr9PYt779Io0FpDcN5znfKF0yagHlJWsZmp8zJ0v9pQr/TaoZnbJZ12UAqFT59WdGXaI8ps/ZrvDO4VnVSATI+XJ+v+Cb4sfE7amQLs4AMOryDezdbV6vaDjm4wLrem14VUOjzlFBWicAwPY1fh0Nw55ZmE4//5LRkoN0mTDRzCmEUyROIZwicQrhFIlTCKdInEI4ReIUwimmlfL8P3+Vxrqb6zRWKaeX3+sN3knYupQ88lg0fl+yMrNSeN2hGuk0Ddg1YipGl+dSg9fTqVWm0+fLDNvJ+EkNNX5vIRjZMYN01seAZIkAwGjEM0UmRvdwGNdRYhk8RnsHVPmzmm5aMf5etepGNks5fW/lwLOuQsFtG4ZmTiGcInEK4RSJUwinSJxCOEXiFMIpEqcQTjGtlKWFNo0t967QWFGsJ4+3Z2f5hRjtGDZXrtLYtU1egGpUpJf6J0ZWRDQKjZkY1kelvsg/r5x+xmOj90NmeCkNIwOmWed2TzEiGSsTbnugyq8jWHaVkfFRJ3bVLOlSDgBHWtyiO3JwnsaMJBIM+ryFRhbT9lIp5/fcafPvhX7OTY8QQtwWJE4hnCJxCuEUiVMIp0icQjhF4hTCKaaVEke8ONJ0k+/av9ZPLzWPii065oEHH+HXcZBbMFdWVmns8mq6G/IW63gNoNu1ul7zAlmTMc/eaJbSmScA8OBj9yePXzC6YV8xMoJ6Q24t9fq8RwnrK1Mt8++5aRQ86zS5dbDQ6dDYgUMHksdPHF6iYxarPGNlyyg0trbG7cDcKALXaKaLr7Wm+D3PzfGCbQzNnEI4ReIUwikSpxBOkTiFcIrEKYRTzNXa1QvpDtUAUIz46mSP1IHpnj1Dx8warRrma3zTc3nAV1frpM1zL+ebuWPkK7JW52KrLk63l141BoAfeDK9Sv3IQ++mY86cOU1jq+s8SWBA6gQBoBvcS0btnnrG73neqLfUafLvsyDP+OIKf3deW1mmsWB0tm4v8tpO9TbfTN+YSl//7Dw/X2uar9gzNHMK4RSJUwinSJxCOEXiFMIpEqcQTpE4hXCKaaUcMDacnzvDbZbxgNgRgdsUb333NRrbqPDaN9avy/YkXR5/e8zL5k+Mze1grQIA5IHXj7Hq0fznv3wlefyHm7x79aNGl+feNLcAJmNuBYVx+r77Q26ZbRgtBljSAQCcfvUSja300hvV+2X+fOuL/D2dOdChsWqbv1e50Y6hMZ2u+1RtcIso5KbUkmjmFMIpEqcQTpE4hXCKxCmEUyROIZwicQrhFHN9956T99DYplGbZfscW0bny+F9w8JYG/MWCRWjbcGQZJgU0cguiXtrxxCi1VGaj3vjW/+RPH72Grd7FjJeqyZGbvcUhgWzRTJ4LpLWAwDwhpERdM5oedFt8O9s6p6DyeNLx+6jY2od3jYEmfGK5/x5tFrcymqQjJWszDNxYrj5eVAzpxBOkTiFcIrEKYRTJE4hnCJxCuEUiVMIp5hWSnuG7/ZfWOLdmpeJlWI4CmYD5YFRWGtkjGOWSYE9dq82iEbGinXjo166RcL2Cm8VkFU7NJYPuPVxwXiOLyJtfbxR4s9qu8WLsjWP8PYDC4cO0djcQrrtQrXJM0iGxrOPhjVWLfHiZbkVy9OxvGR0IydjLDRzCuEUiVMIp0icQjhF4hTCKRKnEE6ROIVwimml1I0eJVWjF0aZdAUuRnxZ20jqwNjoQwLLFmHDrA8zsjosJkbqSTRiW5P09b86NLqKV3hWyqt9XjzrlTHver1Gil3N3nOMjjl4lFsiHaM4XNUoXpZN0s9qZFgieYkX48qNTJFShY8LGf/OiiJtSQXje86UlSLE3YPEKYRTJE4hnCJxCuEUiVMIp0icQjjFtFJGRtGt7R7v/zHVqSWP97d50aeCWAoAUBjL0IXlfJBgMOp72bkznGjYM9Hok7GdpZ/x88MNOuZ01yiG1uDPqrTEC7YdOLyQPH5sYZ6OmZvmbdYzwy7ZNrJI+sQ2KxlZIjXD1qsZ/UtKlfR7CgC1Os+CqdbS48plnqWzFzRzCuEUiVMIp0icQjhF4hTCKRKnEE65wWotX13NK3zFbWYhvUI2avGNxmNjU7wRwshY5Y1ktZZ0HgAABGO11trYbG1uR4mv4pVKZKO30Vl5MM03lR+f5rWdZmZ524JWO/0qtBp8lbRa469P3+iiPTRqGUWy4pmXjVfVevZGrGxsfLdqCJXJtbDaQsANakwRNHMK4RSJUwinSJxCOEXiFMIpEqcQTpE4hXCKaaXkZb4M3ZnlG5tbZPN1MeTLyZaVMi4Mu8SwPjLS1TgYv0mZVQcm40vlWcnYcF7m910nS/ZTU3zD9lJrmsZaVV5fqGnUHqpU0xbG0NjLvUVqRQFAz0iasBIZasR2qhjJA5YlYrVBCEanb6tD+HCY7jpeqfBu5JWy2jEIcdcgcQrhFIlTCKdInEI4ReIUwikSpxBOCdaSsRBi/9DMKYRTJE4hnCJxCuEUiVMIp0icQjhF4hTCKf8NBnlQTNnC8YwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcU0lEQVR4nO2deXCd13nenxf7QgAX+w5i5U4QpLhI1L5YUmzZpGJVlqs4SsZ17dRJmyaeTDudNPVM3Ewz8UzjyR9VHduVLdtSVMWO5MgUbWuluO8rSJAgQWIn9vViPf0DV1OWOc8XORNXh/bzm/EYOg/Od7/v3u/Bd3ne876vOecghAiPpA/7BIQQfmROIQJF5hQiUGROIQJF5hQiUGROIQJF5rxFMLOnzWz3h30eDDObMLP6D/s8fpkwxTmFCBM9OYUIFJnzH8HMNpnZMTMbN7OXzOxFM/tTM8s3sx+Z2XUzG078XHXDvCtm9tAN//1fzOz5xM8ZZva8mQ2a2YiZHTKz0oT2W2bWnni9y2b29A3je2443l+a2TUzGzOzI2Z2902v9Tdm9u3Ecc6Y2eYPcK1XzOxLZnbSzEYT15pxg/45M7toZkNm9oqZVdygOTNrTPz8UTM7m3jtLjP70g2/95iZHU9c914za/75P5VfDWTOCMwsDcAPAPwvAAUAvg/g8YScBOBbAJYDqAEwDeCvPuChnwGQB6AaQCGALwCYNrNsAF8D8GvOuRwA2wEcJ8c4BKAlcV7fA/DSjUYC8AkALwCIAXjl5zi3JwE8CqAOQDOA3wIAM3sAwJ8l9HIAHYnj+/gGgM8nrmEdgDcSx9gI4JsAPp+47mcBvGJm6R/w3H6lkDmjuR1ACoCvOefmnHN/C+AgADjnBp1zLzvnppxz4wC+AuDeD3jcOSzdnI3OuQXn3BHn3FhCWwSwzswynXM9zrkzvgM4555PnMO8c+6rANIBrLzhV/Y4515zzi0A+A6ADR/w3L7mnOt2zg0BeBVLfwAA4GkA33TOHXXOzQD4jwDuMLNacn1rzCzXOTfsnDuaGP/XAJ51zh1IXPdzAGaw9D6Lm5A5o6kA0OX+31WzawBgZllm9qyZdZjZGIB3AMTMLPkDHPc7AF4H8IKZdZvZn5tZqnNuEsCnsPQk7TGzvzezVb4DJL5+nkt8/RzB0pO46IZf6b3h5ykAGWaW8gHO7eZ5yxI/V2DpaQkAcM5NABgEUOk5xicBfBRAh5m9bWZ3JMaXA/jDxFfakcR5VyeOLW5C5oymB0ClmdkNY9WJ//9DLD2ptjnncgHckxh//3cnAWTdMK/s/R8ST+EvO+fWYOmr62MAfjOhve6c+wiWvjq2Avj6zSeV+PflH2HpK2a+cy4GYPSG1/5F0I0lc71/DtlYevp33fyLzrlDzrkdAEoA/BDA3ySkawC+4pyL3fC/LOfc93+B533LInNGsw/AAoDfNbMUM9sBYGtCy8HSvzNHzKwAwJ/cNPc4gKfMLDWxGPPE+4KZ3W9m6xNP2TEsfQ1cNLNSM9uRuPFnAExg6WvuzeQAmAdwHUCKmf1nALn/PJdM+T6A3zazlsS/Ef8rgAPOuSs3/pKZpSVisnnOuTksXd/71/B1AF8ws222RLaZfczMcn7B535LInNG4JybBfDrAD4LYATAbwD4EZaM898BZAIYALAfwK6bpv8xgAYAwwC+jKVFm/cpA/C/sXTjngPwNpa+6iYB+AMsPaWGsPRv2N/xnNrride7gKWvmnEkvm7/onDO/RRL1/Qylr5RNAB4ivz6ZwBcSXzd/wKW/r0K59xhAJ/D0uLUMICLSCw4iX+INiH8nJjZAQD/wzn3rQ/7XMQvN3py/iOY2b1mVpb4WvsMlsILNz8lhfhn54Os3v2qsxJLCxrZANoBPOGc6/lwT+mfhpnVADhL5DXOuav/P89HRKOvtUIEir7WChEokV9rHS7Sx+rQ/w3b/QMK25Z5x9v3z9A59ev3UG1hjr/Wuxv535eqEf+usMYuHvMentxLtfbGB6h2W8kU1ebODVPtRHqfd3xzWymdM1Q2RrW8vPNUS67dSbXzl/3jk1UdfgHAyr5iqrk078amJW1xC9W6h3/oHa88sY7O6XyqkWqr9rdTDSt4htvxKb6ssKKsxDvevriJzhlOG6La3Sjwxqf15BQiUGROIQJF5hQiUGROIQJF5hQiUGROIQIlMpRi/XyJ2lInqNbW5A+ZNDWN8heb40kVZ1PWUi1zmG/WyX17wC/s5KcxNLudarcNkOMBcLM8dHDp0J183l3+9zFe50uTTGhpPGU0o3Qnf60FKqGXRIIyL/XTOefhDykAQFLBHNd4ZAx1aY96x5dV8FDEqsusIAPQVcY/z8nXOqmWf6f/PABgfu4t73hVJp2CyUsFXGzwD+vJKUSgyJxCBIrMKUSgyJxCBIrMKUSgyJxCBEpkKCU+cZ1qMyN8Ob8+luodv35+ns65topnKmyK82X04Rm+ZJ/5Sf9S/zWconOmRldQbXd1EdUqh3goaLqMv821/Su9471reDZFqZul2u5OHlpqThqkWlbtcu94Tpx/LkkR731xjGcg5VMFOLfgD2+Mrc+mc7bN3EG1gfkaqm3Y3Eq1cX6roqsn5h0frZ+kc4oa+Pkz9OQUIlBkTiECReYUIlBkTiECReYUIlCiq+8NHeDiwDYqtY33esezm/nhKhyvi3PpYh7V5lbzGjfT5q8703KEr9INJPtXmgEg3sI3o2e285XLgwv8b+CKBv+Kcl4XX6G+kOWvOwQAK3r57uviAb76fvpe/6rsez+N0zkFKbzI/N01PBFgzEjBIgDXpv6Vd7x+BT+P2CK/5uSf+OtZAcDutXzje1kRX2G/a8S/Ej0d44kA0xO+rhpLFJRnq4aQELcSMqcQgSJzChEoMqcQgSJzChEoMqcQgRK58X2kgIdLYm18Xl2jf6n8EnhNov40rjVk8iXvPcbPccOb/nlD9/ON40loolo1+GbuzulXqPZwjJ9ja5d/+b2umodSeBMEoPhsLdVmank9oJp3/BvmU2Z4K4lXBnlNqHX1/D3OXIxRrbbM3+S6fJCHuGbG06jW+WtUwhOdVVzM4Tf4scv+e2R9XkT46Cq/v1HuH9aTU4hAkTmFCBSZU4hAkTmFCBSZU4hAkTmFCJTorJThTioOX+XL0JljXd7xjLujIje89g0OtFCpfxuvtXPxsL9LdWPqQ3SO68zg51EzTaWF9RGxpTEe/JhJ8bd4iCXfzeekv0W1zKtbqXa+6xLVplLHveN5i/7PEgD2zG2gWk4H7/3waB0P6ZTc4X+vXm7llYfuzeSvlR3xmWVd53WfZsZ4+4SFRv9n9ualQjqnpIiHe7bmQVkpQtxKyJxCBIrMKUSgyJxCBIrMKUSgyJxCBEpkVsr8KF/yTi/oplpHhr9VQ+PwejrnYESN/jsiOiGX9IxQLbl5p3d8qJu3kmgt5cWiPl4Yo9r1EZJaAKB4kLdB2FeU5R2vGyOtpgH05/A367sD+6lWPclDDh05/lDW2lnegiL1IA/NvFJ5kmrnHWnlDOD3Z/2dvhsmY3ROWhUvnnXhYkQRryae6VJewkOFe+EvGlZUeY7OqZ+PakLhbxmhJ6cQgSJzChEoMqcQgSJzChEoMqcQgSJzChEokaGUVxb4zvzs2gqqrcAh73jHEH+5XPCwDVBKla4SnnVQecHfR6XQ1tI5TUnnqdYZ93ehBoCi5NNUa8vihcHm2m7zC5t5COD6Xh7u2QgeLplt4O/jylLS5+MVHhIZKOYhjN/Z9BmqnXpjN9V2r/SHUrbdxu+3NvgzagAgY76OalPnaqk2msb76dQU+7Odko7voHOK7lFnayF+aZA5hQgUmVOIQJE5hQgUmVOIQImsITTqfkzF+BSvcz9DFqYGTvHOysnr+UbvyjG+UjfMy8Agc8LfLbsqi9dzwYS/HQAATOWepVrKzBqqnRngZfoXk/xv1qp8viLblsqTDmoO+Lt5A0BrjK8AN6ZPescvFfFuze2TfKV8+zzvOD4xT1aoAeS2+9+PopGrdE7mEzyqcHyMv4+bcvnncjGidchBl+4/3nVez6oqmSdGLCtcphpCQtxKyJxCBIrMKUSgyJxCBIrMKUSgyJxCBErkxvc/+6K/tgkAXNvxJaq1nLndOz7dyUMR9Z9dTrW8uSeoNpwVp1pTkr9uy+HCeTpn8xjv5HwhdxXV6oaep1o8fRPVCpx/U3z2GZ4IMFvP2yCghnebzorxVhPjy/wbxFPHefjr8Tz+mWX18/f40BgPwRyd9Ic3FgZ4p+yKb/D6PK818oSEY2v47Z+/eJxqW0r990GshIcl06b9dYei0JNTiECROYUIFJlTiECROYUIFJlTiECROYUIlMhQyrXNx6iW38ozATqGfuAdz1l2F50ze4D3XCh70p8xAQAF+3jWxMSGIe/4ppm/o3NOzfw21Sr9DY0BAFbOO1GvBg85zC6QVg2pvF7RVGcf1aZLeZ2glusRXcz7/OdxoYFn6Zz6C/4+Ht7mrwUEAPPmTcIAACT1+ltlDF7nIa5LjofoZuM8FFTUycMbA6X/jmqtPf77cWsXDz3OtPBMorRKf9aVnpxCBIrMKUSgyJxCBIrMKUSgyJxCBIrMKUSgRBb4+tbl/VScP8CLZA086N+133KCd5SuKjtFtaJ5Hra5PMLDLFmN/m7NGX28jcB8MX+tdTM8M+JaBj9m+hTvDn1xRZd3vK+HtxG4K6IFxUBEcGz1ed5Z/HjtCe94yyBvXdE2yDM+9nfxcEllJT//g3P+rJoNqZV0zuG+Eao9uJ+HS67fV021fMe7XtfF/dddnc+LvGFzVNRSBb6EuKWQOYUIFJlTiECROYUIFJlTiECROYUIlMhQyj50U3FTP0/RmCjx97uY2Ms7PF+d4svQ2zb7s0sA4Pzb56h29JR/Ob906xyds7aKh4i+W8g7KO999kf8mBNZVGvP9Wfj5C7wHh/LiqmEh4p42Kb0dD3Vpn/PHwq6PY/3qTk7wft/jPmjWACALTkvUG1mfJt3fBTX+HkM30O1dZM8G8Q9wouQZY/ysF9Bnj/MMnSVh9oKaqqoBkChFCFuJWROIQJF5hQiUGROIQJF5hQiUGROIQIlssDX3q+8SLUjl3nhp+tl/n4Xg5kb6ZzR+89TrbOVF8jaOM77ZOQtf8M7XvTwI3TO/nbebjzphzxcMnCMt3R/tvUBqm3N+Z53/PDEvXRO8cJJqqVv/zjVkjbwImp/nHSndzweERKZz/RnsgDAumrez+W1kV+nWvMRf2jsnYf4/bZl+VGqXe3h99XMIX/YBgBSK3j46+KVMe/4Qxt4BkznZf8cAGiuy/WO68kpRKDInEIEiswpRKDInEIEiswpRKBErtY+t9/fVgEA0t/YSrWsnVe84y0ZfDNxfKiBap3NvKR+T41/pQsAbivzr0Amo4nOqZtop1rO9tVUOzXN69FUrOH1aDZ3+o+Z80W+kvjM5L+hWvpqvpl70zTfmH2i11+75+i5PDrn0UXeBuGvs3hCQvUEr/v0Qrm/o/fsX/Ou3DbKb+OyOX+XdQB4ftUBqo2N8df7xOFL3vFDX+X3QOYwb6/RXNfiHdeTU4hAkTmFCBSZU4hAkTmFCBSZU4hAkTmFCJTIGkL1sYeouGIzbxdQVtriHV+zZYHOGTzIa98MFvAl+4JkfwgAAB4s2+EdPx5rpXMeSObL4fFUHhIZKeDXVnJfjGrlR/zzMnJ5ksBgRGn/vr/7FtVy+/hm9P6OI97xl0p5CKDn6j6qlU7wthYT2ceplnvNX2dqbuNH+Wu1HKba5F4eCnpwrpdqy7f+T6rtmtjvHW9Zt4nOGc3cTbWnVz+mGkJC3ErInEIEiswpRKDInEIEiswpRKDInEIESmRWyuXSYapNdP2Mav/hj/7KO/4Ha3mNlb+s8O/0B4CSM7w1QW0fv4T5XH+Nmzuz+XVhWTOVKqr50vv65bzmz2Q673rd3+LPFFnbxs9x78U9VNvVy3s1rItoMTCXVeMd7y7jNZV2tPDwUVIG71594S946436Hf5jXpjqo3O2n+chopk7/F3WAQBx/zUDwFjWc1R7vOwZ7/ihiuP8pc7zex8k2UlPTiECReYUIlBkTiECReYUIlBkTiECReYUIlAiQym4Ukml5AcfpdqT5U95x6cmCumch8ELa40+wsv+ny3nZfOzlvk7Oc+m8/O4q41nUwzfxguNTbdOUG00g3fLbhj1Z00c2BajczrQRbVNGf6iZgBQODBCtblMf5fqP+ni3bwP5ZVQbX0/D0m5Hfz8rzb6b8n6Il6crO0Iz+BJT+Xhr8FZnp20+CrPTprM+qp3vKyOh5Z662NUA+7zjurJKUSgyJxCBIrMKUSgyJxCBIrMKUSgyJxCBEpkKKWmaoRq/z7jE/yg6/3hgT0vnqJzunfwzsX5XTwzYluhP1wCALn9/tDHRCPvJ3J0nodLOna/R7W3SnjPlnuOrKXaaz1/7x3f/BwvFpXfwM+/YfkLVHt1mBdza0jzF1g7/YK/4BYA9IJ/nmOxAaqhlmfcNFzwH/PHEZk9U8m8L8u2l/l7dSmL3zujG65RrWrxIb+Q5O/zAgApu+JUA4lK6skpRKDInEIEiswpRKDInEIEiswpRKBErtY+/clPUy3nCb46+dp3v+Md75/hnaEvf5u3VYhHdC5eX8lrywzN+jdtT73BV9VG4nyD9YbqMapdebOWarfNvEu18sJB7/ir/d4K/QCAxQleJ+j07ItU2zLPW168e9m/cd8iPrNL2WVUyznB3+O5Yb5B/PIa/8b94yMX6Zz7W/lK6PjtfAN+bLiNarUFpVRLa/afy/Suo3TOilWbqcbQk1OIQJE5hQgUmVOIQJE5hQgUmVOIQJE5hQiUyM7Wj3z6P1Ext5uHFUbG7/YL67vpnMlVGVRb9dZpqh1b4GGW6ZKPecdL8g/SORlT+VRrXs3rBOW6LVQbaOCb6Yt+stw7fvY3eIgo9uoFqmXwfeqYWnGean3JD3jH6+d4SGTZHt7pu3cL3xRfFedtId5+3X/Mu1en0zlnz3Its5iHREYyeQuQynv5vJpRf82i7R/lHdgPDfP79Et3f1GdrYW4lZA5hQgUmVOIQJE5hQgUmVOIQJE5hQiUyFCKrbmfi908iwRZ/iX7tT38cJN1fOn6ypS/VQAA5A/MUS27qc47Pj4UUR8mk2eDlBWuoFrfbIxqyUW8/cBUnn/e7dm81cHE2HGqFczyTs4Tefza2sk5Nl78OJ3z3oWXqNZSyrtXt1ZRCSW9/tpD+e28rUJ387+lWt7qY1RrTL9MtYV1POtq9B1/PaPyNTyjqWktP97vPfxlhVKEuJWQOYUIFJlTiECROYUIFJlTiECROYUIlMgCXw8X8xL4u4t5J2q8M+UdHiu9h06ZXLGfH6+Hd5ue6+OhlM5WstTPoxSY7+AZEwUdV6l2usR/zUsiD8Gsz/SHpF7ccoDOcb08g6e+mGszl3gWRk26P9OidfFndE5pPKIz9ENHuNZeS7Wrk7d5xzM3tdA5GxavUK1ygBc1a6njWUbxQX7vH4/5WzwsVvN7cV3bw1QDkfTkFCJQZE4hAkXmFCJQZE4hAkXmFCJQZE4hAiUylHJimmdv/Iv5rVR7KeMn3vG6tbww1XDeeqqt2s//hjT9KS8y1f0V/+WdrOTdjvv6efbAqUyeXZI2xs+j/M5Cqi07keUdX+ji8Z6Vpbw3SN84r/A1CB5WuJLq/8xq53nBs9psnl5yYZSHdJrnbqfaqkf8mTMpeTwLam6eFyFL6q2k2snObKp1N/JjNo34s50qL6XROVPN/P5g6MkpRKDInEIEiswpRKDInEIEiswpRKBE1hCqyNlCxZ6qK/yo58nKqztBpzTd4V8BA4DKogqqDQ/z06i/4t/ofQZ8w/OF63xzO/wLqwCAGl4CCdOdfDP955/8fe94VeEmOmd+aohqbXO8rs/VI6RNBoCZ6/6N3qWF/DNLfqKeaveO8/YUyT18JTpW6l8t/27By3RO0p7X+PEqf5NqmUNnqLa9cR3Vrq31t3/YcJbfp907eTfyT215QDWEhLiVkDmFCBSZU4hAkTmFCBSZU4hAkTmFCJTIje8pK3jNGUzwGEaze9M7PvYv+eHmvsfDA+c28Y7MlUd5+4EfFPk7Da8a5OeBlBoq1dTyMMvYsUyqzYCHHEbP+ZMLmj/DN45fS4lR7c97P0W1XWt4S4MNE/5N8W+e5CGRnkpeR+r6q7wm1IFJ3ll8/LR/M/2jTbw+TzzG74HhZfwW35L7GNXemuehj4d7j3rHM3byDtsFC7w9BUNPTiECReYUIlBkTiECReYUIlBkTiECReYUIlAiQylbV99HteoOvlS+94J/yb7oezykMAAetqk8ykvjd4F3a8aAf2m7axWvs4MeHi4Z4E2SMYWoGjE8++Hr+/x1lVJbH6VzntnCw1i7sqeptiEiLLIn95ve8fRB3kpi3HiNqZ++10S1ijwertq83N8dui/V36YBAOzTF6l2/nUeghnO5OG72jQeb7OGmHe89SS/Fzdv4ZlQDD05hQgUmVOIQJE5hQgUmVOIQJE5hQgUmVOIQIkMpWwpuk61Hx/soZrBv4w+AJ4VkVXOz6OrJ6JrdCHvrrxxMM87fqwtYlm7nIdEpsob+bxFf6YCACCFh5A2nV3mHU+vOUnnjM5splrNAH8/3r30BtXe2eRvu7Avq4POKZvirRrGl52m2mz55/i8bf54VUkhL5KWdZxXV5sq2EC1wo/4s5YAoDqbh9sOdfhDQU9X30HnHCvpo9pqMq4npxCBInMKESgypxCBInMKESgypxCBInMKESiRoZSTaQ1Uq83moZT2bH/oo3tynB+vhzci4WWkgPggD4scy/O/XmYj717trvEshtk0fo6LJ3nPmcxsfgXs1QZ7ebjkxBn+Pl6a54XGXo3IjhlL9oed7sn3h3oAID/O36vGnR+nWnE971I9e9Hfo2RXDn+ttQv8vVoz4+/YDQCZu3jRrf51PISUu9rfLftYB+84nrIvImvpcf+wnpxCBIrMKUSgyJxCBIrMKUSgyJxCBIrMKUSgRIZSkmZ4Vspz6bx1eG65f6l/4SovqDTUzAtTxS/yMEVGWhrVsrP9xaJWHZmkcw6kL1Btsf9dquXzCAaGI9rV9037QxiH7DCd0z7Pi27NFfIwUX3OR6iWlOXPmujO55kWj89tpNqrc7y/Td19/M3KhP/8PzLIM0hOtPG+MiXrvki18bPPU23+Ph5GXN7pv39iFSf48WZ55gxDT04hAkXmFCJQZE4hAkXmFCJQZE4hAiVytXayjddtuW9jNdV6D/jL49cs8hXZ5F6+wbpurX/VFQD2veffhAwAGXWXvePvpfBaQJvz+cpwyry/dQIAZPAFYAyn30m1huwr3vEDd/EuycnxHKo9kVxLtWMb+ebxum/4V4Dj2/nncux2vtH7qeV8c37/t89RrafUv6I8Ob2ezsn5LG+d0PzfeN2kgZ38mPFufsyCppXe8avtPLrRvNhJNWzb6h3Wk1OIQJE5hQgUmVOIQJE5hQgUmVOIQJE5hQiUyFBKVglfDr/2Hm/zvDLFvxH54gr+ct2FPITRNMY3WH/ssQNUOzu63Ds+urWYzsm5vI1qUwtfpVruDA9vpO3op1ryaf/fx6J1/L169MgnqZZZ9DOq5Y/yztaxb/pbZeQc43WHftrGOzmnzPFQRHylPxQBAPFR/4b5ihTeDTvlb1uoNvC7+6hmb1VR7cEneZjl9DJ/faHiBX68sTi/ZoaenEIEiswpRKDInEIEiswpRKDInEIEiswpRKCYczyEIYT48NCTU4hAkTmFCBSZU4hAkTmFCBSZU4hAkTmFCJT/A1VUgS933+TGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUKElEQVR4nO2dW6hl6VWFx9xr388+97p1dar6aksUbAPxRR8MXlCQzouCxNuL6IuggpeABvQhGhDfkkclYsRog3SImnhpVFQCPhkiiaTptk26q+mu6rqcOtd9/X2oEyyKf4yu2p6qninH91S15vnXXmedNc6qmuOfc0YpBcaYfLTe7QswxtSxOI1JisVpTFIsTmOSYnEakxSL05ikWJzGJMXiTExEXIyIvYho7uNn/HdE/ACJfSAiXr9fn200Fmci7hRKKeXrpZRRKWX+bl6XeXewOM19ISLa7/Y1fLNjcZ4Ax2+8X4uIL0XEfkT8YUScjYjPR8RuRLwYEZvHX/vBiPhyRNyIiH+KiPceH/8UgIsA/vL4n7K/HhGPR0T5xoMeEecj4rMRcS0iXo6In7vtGn47Ip6PiD8+/swvR8T77/Jb+K6I+EpEXI+IT0ZEn3yfJSKevu3vfxQRHz3+8wci4vWI+HBEvAngk8vcS/O/WJwnx48C+EEAzwB4DsDnAfwGgNO4dZ9/MSKeAfBpAL98fPxzuCXGbinlpwF8HcBzx/+U/b3KZ/wZgNcBnAfwYwB+NyK+77b4B4+/ZgPAZwF84i6v/ScB/BCAp46v/yN3ue5OzgHYAvAYgJ9f8hzmGIvz5Ph4KeWtUsolAP8C4N9KKf9eSjkC8AKA9wH4cQB/XUr5+1LKFMDvAxgA+O53OnlEXADwPQA+XEo5KqV8EcAfAPiZ277sX0spnzv+P+qnADx7l9f+iVLKa6WUawB+B8CH7nLdnSwA/FYpZVxKOVzyHOYYi/PkeOu2Px9W/j7CrTfe175xsJSyAPAagEfv4vznAVwrpezeduxrd6x987Y/HwDo3+X//V6745zn72JNjSvHv4zMCWBxPljewK1/8gEAIiIAXABw6fiQqt97A8BWRKzeduzibWv/L1y445xvkK87ADC87e/n7oi7/vAEsTgfLM8D+JGI+P6I6AD4FQBjAF84jr8F4MnawlLKa8df97GI6EfEdwD4WQB/cgLX9QsR8Z6I2ALwmwD+nHzdFwH8REQ0EfHDAL73BD7bECzOB0gp5asAfgrAxwG8jVuJo+dKKZPjL/kYgI8cZ3J/tXKKDwF4HLfebC/g1v/vXjyBS/tTAH8H4L8AvALgo+Trfun4mm/gVhLpMyfw2YYQ7oRgTE785jQmKd7F8ZATERcBfIWEv62U8vUHeT3m7vE/a41JinxzXr58mSr3YRX1LXfjAbLEbVz2zi+zbvm7IT4tTvhKHvCjyJ6RZTVx9uzZ6gn9f05jkmJxGpMUi9OYpFicxiTF4jQmKRanMUmRVkrT8L5StlJOCnIf5e1d7t4XZUeQUCl8TUijhV9jLGWliE/6JrdSGH5zGpMUi9OYpFicxiTF4jQmKRanMUmxOI1JirRSVGr4YbVS7sf3tZw5c/LXoc64IJaJvh/KLuGrQpzzQRtZJ4mtFGP+n2BxGpMUi9OYpFicxiTF4jQmKTJb++A3gZ8sWTLKOt9Zj6prX8wXNDZVsRmfwTuZ1mOz2Yxfx4Kfb2U4oLFumxdU9Hvd6vGWeBbvx095mWf/pPXiN6cxSbE4jUmKxWlMUixOY5JicRqTFIvTmKQsvfH9YWXZdLhu+aOiZMO5WDEX28OnC77yYMJtkcPDcfX4bDqla5qGX8dwwK2UpuGPXZD3hf65LPmcinNmMBH95jQmKRanMUmxOI1JisVpTFIsTmOSYnEakxRppchKgBO2WbJUFihkDx5hYcwW914pcjThFsbe4RGPHdQtEQAYT1WFSf362y1hl/TrFSQAMBHVMT1xI+n0h3ufJHFXUboqQUWW35zGJMXiNCYpFqcxSbE4jUmKxWlMUixOY5IirZTxeEJjymYJkn7Xa/jviWUrEtRUZsZCNdYSdslENLsai2qQI1L1sU+qRADg6s09Grt+c59/lrBnQL439du73+HRm6tDGjuztU5jp7c2qseHgx5dowdl33tF0DuxTC3LMp/kN6cxSbE4jUmKxWlMUixOY5JicRqTFIvTmKRoK0UkgDuiSVNDrI8ict6i+AGijxRCJanJxy0Kr5hYzMU8ETE3ZCIaYY1FrNWq38duh88TGfX793w+AJiIWSlzEpuL+9EWv9qbDq9YmQdfyMy7nrDa5BtGOCknXXiiLJ1lqq785jQmKRanMUmxOI1JisVpTFIsTmOSoscx9FdobC4zXWwTNc+ShoipzegtEWMR1f9I90bisRCxtkjj9br1H8Ggy7Odwx7P1pIB1QCAUGMQSJZX1Q60xITqtkjlthseYyGeJwda6ufSEs+OOqfK8tLAyaZr/eY0JikWpzFJsTiNSYrFaUxSLE5jkmJxGpMUaaV8+vm/oLHVFd4jZpv0iHnk9DZd88jpLRrbWh3RmISkr2VWW6X5Q1kR3HNodzo01iEbxKPFbYq+8DfmIrYQv4tn5KaMxWb/wyM+FmL/4IDGDg54D6TJuH7OEMUKynZaW+V24ObGGl8nnm/WC6uoHkJLbLL3m9OYpFicxiTF4jQmKRanMUmxOI1JisVpTFKklXLj+k0amx7xcQEtMsl50HB7YH3AKy1WRRWG7C9E+s4omwKiv02I62/A0/mNnAhQ/wbE5AdZpTMVYyEOx4c0dnO/bn1cu7FD11y5eo2fb3eXxsZiMnchz05XVMAoW+/0qU0a64jKmRXxPDbiOThJ/OY0JikWpzFJsTiNSYrFaUxSLE5jkmJxGpOUUA2tvvAfr9JgI5oZdUmKWjW66oiOSlF4ZcRiwWPMMlGNropIk6uqg4WoBmkJ64al5VW6Xk36noufpxoncUCssb0Dbr/si1iZ8yoSzEWDtVl9nSrqEIVEGI24JXL2LK+SWhtxe4Y1KFPXqGLb6xvVsN+cxiTF4jQmKRanMUmxOI1JisVpTFIsTmOSIqtSmsK12xHp/IbZCqTiAAD2j1TKnlfH7O/zZlF7pMnUWAwU2dnj13Hprcs0dn2HV2H0RFXNiFRUrK/xpmZbmxs0dvoUtwdOnRJN1Dbr1Rtb67wJVpGlMyIkbBY2WVx9VhHNvxrhs/R7PRqbyWFA9ZiaASO9FILfnMYkxeI0JikWpzFJsTiNSYrFaUxSLE5jkiKtlM/81d/QWK/Dl66u1K2DDZGW39jisdU1Pu+is1KfywIAvVb9OvZucGvmzavcLvnqy5do7I03+bpej89K2SRzYM6JuTLlIrcOtkf8Pq40vAnZqE1+nqqZmKhymUwmfN10yteRWBETboqoCJoJ22Z+TVkwvCpouDKoHh/IpmD8PXhqbbV63G9OY5JicRqTFIvTmKRYnMYkxeI0JikyW/ul/3yJxkZDnplaX69v5j539hRdM2/xzeizBc/ujYY8k9sf1LNqayK798g5/lmt4Bm8px9/nMbWxWTus6fr9+TRc2fomm0xkbkr5lNMD/gYhJfeulI9fn3nBl1zfYdnvXdFQcKYTK8GgBnJALPRGsA7TA7v8s3tvb4oSBA/s83tepHAxgZ3DtosGw7gyXOnq8f95jQmKRanMUmxOI1JisVpTFIsTmOSYnEakxTdQ2jBNy8Pe7xd/eqwftrWgvfnufz6qzS2ozYNb3N7ZpPYFANhvzz92Dkae+aJR2msESMe1KbnIBvLDw/5vXrl1Zdp7Pq16zR2c5fbG7t7+9Xj+8J+OZrw6eZz0S+qEVOq+/269TEgthgADNvcLhn0xbrV+oZzABiMuJUS7XoBwXjG7cCjKS8SYPjNaUxSLE5jkmJxGpMUi9OYpFicxiTF4jQmKdJK6S54qvzsBrdSnnryQvV4h2fQsXuDWwDDLu99s73Kr2NrtW6ZjMQ4gxWRXm9E9cNcjARgU6MB4OYusTAm3ErZEVUdY2Hp9Le47bRy+nz1uLSIxMRuMa0DjajQ6HbrsY64992uqEoR6yD6Ek2m/Ge2S35mV6/yZ3ix4DYL3v9s9bDfnMYkxeI0JikWpzFJsTiNSYrFaUxSLE5jkiKtlMcuXqSx73z2fTT23m//1urxpsXz6wvR2r8nUuV90cCpwywY0Wp/ISyAWeHXOCcTmQGg0xfjGAb1plDDDd7g68x7xBiEqZgaLRqbtYhloiZKL+Y8NpejGrhNcXhYn0Z+dFQ/DgA7uzfE+eq2BwAcicof9TyC2GYhnm9VmcTwm9OYpFicxiTF4jQmKRanMUmxOI1JisVpTFKklXLmwjM0NtisVzEAwNVdkmIXNkVLTF1uzUQ5yx5PeRcyY2UhKkhmwhKZznnDs7G0B3hsPKlf44ForDUTzaJmosmUmkR9cFT/vCMxoTpE6clCNPhidgkAHJHYbMqvoy0sjI4oheqTCphbMTFjhVTVtKRdwm0nht+cxiTF4jQmKRanMUmxOI1JisVpTFL0OIbBFo29eolPNT46ulY9PlObqFUyS2zYXiIJhiKytSrLuCg8Ezqf80zoVGQ8J5N6llRtDg9x/SFS4oWMfgCAMbnGyZRP+l6Im6965izm/Jzse2uLid2truhlJPocYcbfTYsWv8fzYNO3+WfpTC5Zc88rjDEPBIvTmKRYnMYkxeI0JikWpzFJsTiNSYq0Uv75H1+ksaYlNqo39U3DDZkIDACN2PjetHkPnqYRrfhZa3/Zvp9veG6LMQJqbEEE/966pJ9Rp8/XzBeiP4+wKWbC0lkQL2sxFRvpxWb02Yxfh7KyIurXMRNWymzGr2M+EdZSn1sf7eDPQdPtV4932vxdJx4Pit+cxiTF4jQmKRanMUmxOI1JisVpTFIsTmOSIhO8//C3L9BYFK7rdrueam4LK6UV/FJaYnyCmpLc6dUtk25P9Ifp8UnZ3X59UjYAdLs81unyc7aJFdQSfXEi7kPlDKk+UWukJSJ+7TfCFumSnj89MdJiOOTP1VDYcD1i+QFAt+H3sY26ddMSFU0xd1WKMQ8NFqcxSbE4jUmKxWlMUixOY5JicRqTFGmlrHAHALMpb0DVNPXKgqYR1RSHvLLg8IC375cNqGiTKdEwLHjqPYhFBACtNrdSWh0eizar4BGVM8o+ElU6aroym6zQdPiabp9fx3DAbYq1FX4fB8P6vdoWazbXeGxlhd8PUZyEdoc3Lxv06jerP1QVUtzuYfjNaUxSLE5jkmJxGpMUi9OYpFicxiTF4jQmKdJKEdlwtESw26untockTQ4Ajfg9sbfD57LsiNj+wX71+Ew0rQKZgwEA7S5Pr/eGvHKmOxrwz+vWY7tH3KqaTA9pbAa+ritsomGvnupfW+E/5601/vM8tbFGY6dFbHOt7t+NemIOSeF22rxwi64RpTMrA/69jTZG9TUjvqZHNKHwm9OYpFicxiTF4jQmKRanMUmxOI1JiszW7ly5QmNNW2Qnu/UM37y/R9cMyRoAaInN7W3w7GqnVc+uRptnXdkIBwDYJFk6AHj0/Hkae+KJb6Gxxx97msYYuzd5hlpN0Va/ibtk53tfjBjoiz5HbXLvAT1O4mBcn4q+v1+fAA4AIH2HAGC4scpj65s0NtpYp7HV9Xq2eWWVPx9MEwq/OY1JisVpTFIsTmOSYnEakxSL05ikWJzGJEVaKQG+Qfxgl/f1OSy79TUtYb+ITcghev5MxfiBCRkXUMT5GtFSv+yLkRE3+Ob2lSt8Q/QFsln6qa3T/LOGfBP1pMvv41jYLPvj+mb6nZv14gEAuHYkYmMe2xdWymJARmiIzfIrwhIZbIjN6KIgoSPucWdQX9d0eJ+gWGK0td+cxiTF4jQmKRanMUmxOI1JisVpTFIsTmOSIvO7w1U+j2E85r1qpkdkSrLozzMV9oayUgpE9UMhVSliGoOa5KxGP+zvc2vpxmVe3XP1pVeqx88MeVVEV9yPXXAr6EZw2+lt8rO50uLn2+lya2xGegEBQG97g8ZGZ+qxlU1ul/RWxFTxvrBExITzbp/bLNQyEVZhUaO+CX5zGpMUi9OYpFicxiTF4jQmKRanMUmxOI1Jih7HMOTp8L0urzoYEytlQawNgNseAOQgauGk0KD6jSSKUrCYibb/snKGxw6Jc7MnJn1DpOXfFHbVpS6/jiuD+qMwWeeWQnebN88andqiseH6Bo31R/UmWaqhnPLGQjQhU+dsmnuPtUTlSUucj6655xXGmAeCxWlMUixOY5JicRqTFIvTmKRYnMYkRVopnS6fhNwSMzRAssZlwX2PhbJLlJeiQtJmqTMT3syhiF0XjcZeWwhbZFavZrkq5rnMOvzev81/ZDgSk6h7Z+vWx8ZZ3mhsdZNXzvSGS1R1AAhiOYSyPTr8MW7aPBYt8QzL57FOS1g6LVelGPPwYHEakxSL05ikWJzGJMXiNCYpFqcxSdEDHBpeotHwbDj6w3raezbj5xMui4wpeyZYNYg430ykww9UYzAxZv1INcnq1C2YN/p8zco6rxbqbvJKkY1TvEnW1lY9Nhrx87VFE6xFlz9aRVRoRKfuBXX6vBlXn8wueadYVzT46gi7h9kz0ppZAr85jUmKxWlMUixOY5JicRqTFIvTmKTIbC3LugLA2iZvc98f1DU/m/IM5FxtihcxMrz6OFg/rDYhq4xbW2QZVZHAUGQaR6RP09oaz5Jur/EN5+tr9R48ADAUk5x7/Xp2sumIURgiQ636RclihaZ+/1UPnhZZAwAttWFe9hAS5xR9iRhqbAj9nHteYYx5IFicxiTF4jQmKRanMUmxOI1JisVpTFKklbKxxdPyYs8w5tP6SAAxNBoq817EOIMQtkiL2CKNsEtUCr3T4an3vrRSeGyVbB5fE6MwVvpianSXb9huic3chdgKE2E3FGUpiPsoLSky0kCNR1h2w7my6OZz3hNqPqvHIoTd4x5Cxjw8WJzGJMXiNCYpFqcxSbE4jUmKxWlMUqIoD8MY867hN6cxSbE4jUmKxWlMUixOY5JicRqTFIvTmKT8D6J0+XDsJU83AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVA0lEQVR4nO2da4xd11XH/+vc5zzuzJ0Zjx+TsWPHdlzc1pAoSUtE21BQRUnDS/1AIailhSoCJJCKqBBCoCqV4AsSVT+AUJXypTS8BDIqD4H40JI0rQtR2zzrpHb8fo8973vvuZsPc4NMtP/L9tSP5fT/k0YanzX73H3OPf+7r9faay1LKUEIEY/iVk9ACJFH4hQiKBKnEEGROIUIisQpRFAkTiGCInHeAMzskJn9+FX+7bvM7KWr/NuHzOzo9zY7cbsgcd5iUkpfTintuR7nMrPPm9nj1+Nc4tYjcd5CzKx6q+cg4iJx3jjuN7PnzeyCmT1hZs3Xv5aa2SfN7CSAJ974VdXM7jWz/zGzeTP7GzN78o2roZl9wsxOm9kJM/vlwbGPA/hFAL9jZgtmtn9w/JCZ/baZfdPMLg7O17zsXB8ws2fNbM7MnjKzfZfZPmlmxwZzecnMfmxw/AEzO2Bml8zslJn9yQ29k9+vpJT0c51/ABwC8G0AWwFMAvgvAI8DeAhAD8AfA2gAGBocOzoYVwdwGMBvAqgB+DkAHQCPD+yvj//UwP6TAJYATAzsn3/9b98wl68BmBnM5QUAjw1s9wA4DeAdACoAPjz4+waAPQCOAJgZ/O12ADsHvz8N4JcGv48CeOetvudvxh+tnDeOz6aUjqSUzgP4NIAPDY73AfxBSmk1pbT8hjHvBFAF8JmUUjel9PdYE9bldAF8amD/EoAFrAnJ4zMppeODuewH8EOD4x8H8OcppWdSSmVK6S8BrA7mUWJNpHvNrJZSOpRSeuWyOewysw0ppYWU0lev+q6Iq0bivHEcuez3w1hbuQDgTEpphYyZAXAsDZakzHkA4FxKqXfZv5ewtnp5nCR/fyeATwy+0s6Z2RzWVvuZlNJBAL8F4A8BnDazL5rZ69fwMQB3A3jRzL5uZh+4wuuLdSBx3ji2Xvb7NgDHB797aUAnANxhZkbOcyWuNcXoCIBPp5Tal/0Mp5T+CgBSSl9IKf0I1kScsPZ1HCml76SUPgRg4+DY35rZyDW+trgCEueN49fNbNbMJgH8HoAnr2LM01j7OvkbZlY1s58G8MA1vOYpAHddw9//BYDHzOwdtsaImT1sZi0z22Nm7zWzBoAVAMtY+0oOM3vUzKZTSn0Ac4Nz9a/hdcVVIHHeOL4A4N8AvArgFaw5hFxSSh2sOYE+hrWH/lEA/4S1/wdeDZ/D2v8R58zsH67i9Q4A+FUAnwVwAcBBAB8ZmBsA/gjAWax9Ld4I4HcHtp8A8JyZLQD4UwA/n/n/s/gesf//3xsRDTN7BsCfpZSeuNVzETcXrZzBMLP3mNnmwdfaDwPYB+BfbvW8xM1HO1TisQfAXwMYwdpX4g+mlE7c2imJW4G+1goRFH2tFSIoV/paq2U1Il7QwnnH3GEVdjo+ykrns92ZR8V96nrZo31Y9jgAFCCTv33IXpxWTiGCInEKERSJU4igSJxCBEXiFCIoEqcQQdEOodsR7yO11+XDjIcjiDcfpTPCi7O5e1ucmE6lyF9c8X0Y1dPKKURQJE4hgiJxChEUiVOIoEicQgRF4hQiKAql3I7015mW4trIo9Dnn9+lc7rkRW28jBWtF/+H7oQQQZE4hQiKxClEUCROIYIicQoRFHlrb0e8j1SycfxKA7ukHlDP2fne8byuXlkfZ1xBXq9ijofavebblzfnVQnxJkDiFCIoEqcQQZE4hQiKxClEUCROIYKiUMptibMZ3YlT9J1xKySE0eElidBzNrc3nFCK2zyBTv/mriNeasHNmolWTiGCInEKERSJU4igSJxCBEXiFCIoEqcQQblCKMVxKHu+5puIU+IGxfX+7HFr93DK9czDaxrtTGOFtagGsJJvGr1mIyGTjjPG1pcA49oKEp7xOkkUN+BhvO7PzrrmIIQIicQpRFAkTiGCInEKERSJU4igSJxCBMUNpXhZDK6qma/f6wZw7U2XB/PwsjDyqRZF6b2Yg9eu2fH1r+vlvCiWM42LHW5b8DJMyDk7XmvrDp/k3CKfyMRonZ+zTp4sJ5Wlts5QSuE8/n7Di7zVbSVRuu28s4e1cgoRFIlTiKBInEIEReIUIigSpxBBkTiFCIobSlly3OiFF1VIec17WQyF4yr3sg68UErRX2fI5DpTFNc3a8Kcz9S6c8lN5x6zhJuaM8bLtimqTT7QYZUc9x5U736sH++Zy1M6YyokXOKhlVOIoEicQgRF4hQiKBKnEEGROIUIiuutvehslPZK6jNPbsV5Nc8r6GyTRt3bFU83qq/Ti1t1GwlQCrIBfw0vG+Dax4w4LaU9T25JTtlzPr473m302jE4zwHzd646Du9EogMA/KSJdS5NNXa+9Z2OopVTiKBInEIEReIUIigSpxBBkTiFCIrEKURQ3FDKfx98jdoKp5BNvZJ3No+1+WborRsnqa3NfNcA3FAE7Vvg+NedcIm3sbnvbcDnr+bPhY7gY2peVMG5VSxHwAseeUEgJwoHp8MDLbXTcV6s6oRZvISKmvP0152N6lUSGkvO++JuiidPiFZOIYIicQoRFIlTiKBInEIEReIUIigSpxBB8WsILa5QW2eJ2+rV/GkvXeSv1SJjAKC9YYwPrKyjPo/XVsGhdFs/cPp97s9nzvf1zdAf5yVvsMiYV4bJy1hZdiaysswqBQHLy8vZ4+UyD8B0lnjrh9XOErVt376J2jZNtKiN3cfk3H0vzMLQyilEUCROIYIicQoRFIlTiKBInEIEReIUIiiWnLDCq5e4b9ic2EGVeI0Lx5vMGhoDwGiD24ZqTvGskrjfEw/blE5WStcJmHQdN3qve+2hFC8m4jTRdu+xF2ahBb6c93nJec+6XiNnL4uE2ZyO3SWPzGB+icfvtmwZp7a6kwlVr5KO6XwIKl4mEXnXtHIKERSJU4igSJxCBEXiFCIoEqcQQZE4hQiKm5XiVWJy3cYkclB1qj55Xa87Xtdrp8lKn7Vr9spWeaEDJw1j/388xccZ98sPN/OdYFpNHj9qDfNCaRPjPINnss0zLcaG83OsOSGFUa9/iRPS8d5rGiZywmnJmeNYi4dLKl7YyXn2e+T5cTvpONfMirJp5RQiKBKnEEGROIUIisQpRFAkTiGCInEKERQ3lHLg+WPUlro8LlIv8ikEY7VhOqblhAA2zXJby7uCft7HzjIwAL9zyUVeKwpfPnCE2uacvjIkkoKm09hkaoyHUt6+/U5qe+jBe6mtQcIiXraQVyfNuWT3JrNzeuESr308jaYBcB5hdJ1QCsvIchKa3DBLU6EUIW4vJE4hgiJxChEUiVOIoEicQgTF9db+59PfpraRGnEzAuj1LmWPN4b4Z8GPvvt+aivmueuvdDyoU+38pmevjUDP80A649794D5qW1rkRW6a9fxbsHfXDjpm48QItY3xt8X1vLLu0Il33cCqc68WHFfowsIitc1dmMse7/T4Paw3+GNca3DPds9pl91x3LVDk/kEgi3TU3RM1blXbGu+Vk4hgiJxChEUiVOIoEicQgRF4hQiKBKnEEFxQyknvnuY2mamJ6ht166N2eP33L+bjqkP893QL3zza9S2o8GLy0y+7a3Z40MbuMvbC6VUuVce73vwB6jNa13BauZ49W08vE3g80s8HHFxbj5//AKPVZ1b5CGR004b817HSZogfRAaTb7zvajyi54Y4zdyot2mtulNvN5SQQoreQkV3gZ8+jrXPkQIcTOQOIUIisQpRFAkTiGCInEKERSJU4iguKGUl596ntoWt45S20c++Fj2+O6dM3TMayd5DZ6tJAsAADbWeShliLRJbjitlbsl/7xyGmIj1fk5e05P6UQK6iTnc7PjuOWPnZ2jthUnjaRo5DNdmrOTdMw2Z46z6wgdADwU5HbzdpYYLyRV8cY5RX9Yt+9VpyaRuwySR1grpxBBkTiFCIrEKURQJE4hgiJxChEUiVOIoLihlDPHeNbB+x9+G7XdtWM2e7wo+WvNbthKbds3c1vDOSdzlS87Rau6zvm8mvpVp3pWSUI6ANAjVcOYux7w5zg81Ka2itPSgHWi9ubRX8e9vxKsjYOXLVQ4tq4TSimda1uZ59k4vcX8A9Qr+Q0ZbTmtuUfz7Ua0cgoRFIlTiKBInEIEReIUIigSpxBBkTiFCIobSnnX+3+Q2h796K9Q24Xl/HFzXO9eBKN0QhieG51+8jgu9J7X0di5WzXHU54S/wxkYZElb47OzbJ8U3EA/pvNLvvQiVPOCXlsZmoDz2bpdPgkL8zlC4OdOX2WjklOCKPKmsAAKBzb2PAQtW0gGTxDQ/whWL5ERLH2atmjWjmFCIrEKURQJE4hgiJxChEUiVOIoLje2o/+2i9QGy/sD5w7mT/uONXgOBmx6nhQS8dWIZ7cwtkM7ZbN9+rROB5lbxzzDp+fn6NjypJ7/irOPZ4ab1Pbaif/jp45zZMfknPRx47x7IKVHp8/21Tedzy8Faez9dgIb/U9ZHxtqqzwa1st88WCzAlHDLecXh4ErZxCBEXiFCIoEqcQQZE4hQiKxClEUCROIYLihlKe+cYBaqsUfGi/l980XHGK2FSqfNNwUeFu6JSccxJXf73JP5OGhvhrsa7LANBo8vkn59oSiUmNF206plbnrTA6Ne7OX3TiTl0SQqqPDPMxCzygNj9/idpWSz6uyrpeOz0XllnMDEDp1AJanOfzGGvy53t4In//G+P82XEeK4pWTiGCInEKERSJU4igSJxCBEXiFCIoEqcQQXFDKU/u/3dqO3Fsjtra9bz7fXKCd6jul3wqtS63laSdAQBUhvOhlKEWTxNpOT7v0VEeEmkM8ZBDak5RW7M6nj0+ZDxsU3cyYMy5tsIpuNRZzmd9LJMsEcCvBVQ63cPNmUeDdAFPBb/o/jC/V5va3Dbe5s/VVItns4wO5a9tuMpbW9f6Xh5XHq2cQgRF4hQiKBKnEEGROIUIisQpRFAkTiGC4oZS7rk7XyYeANr1M9Q2vzCXPb79LbxEf63kLu+z371Aba8c4gWo5k7mXf2ri9yt3fX6IHg0+PzbExv569Xz93jZ6f1QdbJqpod4G4GplpNhspLPWLGShz1Si8+jGOEhnbFhfm3jI/lw1Ww73wIBAPa2eYhu57YN1OYMQ7czT22NSj681ODRF0y129xI0MopRFAkTiGCInEKERSJU4igSJxCBEXiFCIobiilabw40r4ZniXw6kq+F0YaXqBj7tv9Vmqr7eIhmOOnzlHbIdINef44L4L12mv8mufmeIGsU6d49sbdO/OZJwBw73t2Zo8/d5678l86N0dt5xd4aOlch/coaVg+9OH1Gpke5eGj2Uke0rlzQ5vatm3fnD2+785NdMx2kn0EAIuLvNDY3BwPBzaaTubS2ET2+GiLZy1NTbWpjaGVU4igSJxCBEXiFCIoEqcQQZE4hQiK663tHD1Kbd1Lh6jt8GS+Ds/CEe4t3Avu+dvY458h7SXundzWzHtXj0/wzdynT3OP7MQE9/KecloCLHTyXmMA+Nn35r3Uv7/37XTM4cOHqe3MBZ4ksLLi9A8nG9xrTu2eceeatzjdptvD3JNbku7QJ89+h445+CJ/Ts3Z7D82yxMSmhM86aNJzjmxgdeKchqwU7RyChEUiVOIoEicQgRF4hQiKBKnEEGROIUIihtKmb2LtyY4+hzfIP7qKtmobtz1fvif/5Xa5vvc9V70nU7O9Xw9oHLLLB1TP8Od3ne2+Pz33ME3Sv/j1/km9he+8qXs8ff1eNhj5jwPzaDq9GroOfWRWGuFxO/HapVf82mnXcfhU8ep7Wwnv1F9xWmPMLRjhtomtrWprTbhtNAY4tJojufDLJUG10t/HeugVk4hgiJxChEUiVOIoEicQgRF4hQiKBKnEEFxQynTu7mL+vTiCWp75GjeZf/yy9z1vnLuG9R2vrOL2urGy/R3Uj5jpTx+kY7ZcvZFavMw8LL/jzzyw9R28O8+lz1+/KvP0jEzKzx8BHM+b6s88weNfAjmZMFbVxzs8rpJR0l2CQAsbR6ltta+Pdnjm/a+hY5pTPMaU6ng15wqTp2gUT7HobF8H4ceeBjLwF+LoZVTiKBInEIEReIUIigSpxBBkTiFCIrEKURQ3FBKa2Ka2jZu2kZtF44eyR7PO8nX6DuJFqttnhnR7XKXd0k6WJdYX7jEI4FfwMu2n9pmjuRL+C/2+b2H8SwdrPLwxpKTsfJsNX+PD47zbJDFWd5mYmTvVmqb3rmb2qY25zOGGmP8fV4Bz7YpE7eNOPej4tkq+TWN5yytD62cQgRF4hQiKBKnEEGROIUIisQpRFAkTiGC4oZShpu8X0S9ybv4durPZ483x3ivlPRTfB7PmNPjw8kGaaV874q7E3flYz8Pe3j0H3mE2l4ynsFzX/dA9viLXV5AbbzFM3FeXOTdmp8jxbMA4PzufAhj8oH76Zgte7dTW3s7DwXVRvhzVennQxhdJ05RONk21RovulWp8zCRl9xTlvnwDAuxAEChrBQh3jxInEIEReIUIigSpxBBkTiFCIrEKURQ/LbzZZfaLi5zl/1IO59BkGo8JLLQn6C2lnFX+en0LWory3yBrwulUyDrZ3hoxmM1neTGCs/gOf5wPpvlK99aoGMOl9wtf373Jmqr3vcgtW3esSV7fMdYm46ZcmxFPV8ECwDmnTXhEikMVqvyMaNNHhJpDPN+KNU6D7M0h3i4qqg4/WiuI1o5hQiKxClEUCROIYIicQoRFIlTiKC43trVkm++Rp3XqmlPE0/dKN8sP+bsbG71+WfIbJ/XFyrKvAd4rM8vuwbuGTZnN/SqY5uqci91o5pvM7C4k3sgV6vco3zXyB3UNjHNa/6MTuQ9wKND/H1pDHGv5UqPe8RXnVYN/Vr+/hc151E1Z1O5Y6s6G98LxzsMWrOI34+0jgpDWjmFCIrEKURQJE4hgiJxChEUiVOIoEicQgTFDaUUNV7KfnKSl8dPw/n6ManDPwtKJ5QyRWq2AEDfcV8XRT50U3HCJRUnJFI4XZK9O1l1NvyPk7L/rRbfOD5Tb1PbaJW3ahip8lBWvZm/ts4Iv78Lw/y15p3kglXj5xyu5sMbzYpTJ6jOr8ucTepW8Pc6JS/0ce1rmjpbC/EmQuIUIigSpxBBkTiFCIrEKURQJE4hgmK+y1gIcavQyilEUCROIYIicQoRFIlTiKBInEIEReIUIij/C15d8mj0C8E2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAANNUlEQVR4nO3dW4xdV33H8d//nDnnzP3mmfgyHhNjagtCaBQeaCtwI0QFqkjhJZWq9DGtVMETAqSWPhQJKJUKgqpSW0VVi7iIvrRAWySgagMRVdMmqoKCIdQmCZOLL2M7Y4/nfs7iwafVyKy1J/vfGc+/6fcjjeTZa/Z/1jlzfmdZe5+1lqWUBCCexn53AEAe4QSCIpxAUIQTCIpwAkERTiAowhmImb3NzJ7ehTqPmNlDu9En7B/CGUhK6dGU0qnb+TvN7Fkze8ft/J14ZQgnEBTh3Af90ep3zeyMmV01s78ys0Ezu8/Mnu//zAkzu2Jm9/a/P2Jml8zsvv73v2Bm/2pmL5vZk/99PPO7TpjZP5vZZTNbNLMvmtlkv+3zko5J+nszWzazD9epjT2WUuLrNn9JelbSU5LmJU1L+q6kj0m6T9Lz237utySdkTQs6RuS/rh/fE7SZUm/qptvsL/S/3623/6IpIf6/35dv70jaVbSdyR95pa+vGPb95W1+bp9X4yc++dPU0oLKaUrkj4u6Tdu/YGU0sOSzkp6TNJhSR/pN/2mpK+nlL6eUuqllL4l6XHdDNStNc6mlL6VUlpPKV2S9GlJv1zRr1dcG3trYL878P/YwrZ/PyfpSOHnHpb0NUm/nVJa7x97jaQHzOz+bT/XkvQvt55sZgclfVbS2ySN6eZoeLWiX6+4NvYWI+f+md/272OSXrz1B8xsVNJnJP2lpD8ws+l+04Kkz6eUJrd9jaSUPpn5PZ+QlCTdnVIa182R0ba13zotqU5t7CHCuX/eZ2ZH+4H7iKS/yfzMZyU9nlJ6SNI/Svrz/vEvSLrfzN5pZs1tF5OOZmqMSVqWtGRmc5I+dEv7BUmv3fZ9ndrYQ4Rz/3xJ0jcl/VjSOd28IPQ/zOw9kt4l6Xf6hz4g6V4zezCltCDpPZJ+T9Il3RztPqT83/Ojku6VtKSbAf/bW9r/UNLv96/MfrBmbewh61+hw21kZs/q5tXUf9rvviAu3g2BoAgnEBT/rQWCYuQEgqr8EMLi4qJrWH21jsZmtvMP7bdez3VaauQfW5LzMffK5zUavj567MVLcbdfBzMzM9mCjJxAUIQTCIpwAkERTiAowgkERTiBoAgnEBThBIIinEBQhBMIinACQRFOICjCCQRVOSvF++n7V+uslNv5uLzPvZl3xkcze7TrLdeoeN93PI/2M4sEvvoxcgJBEU4gKMIJBEU4gaAIJxDUnlytjb7Wjveq6148rt2+Atxr5K+67mRlrZs9/pPnz7vqDY50im3TIyO1602O1T9HktIeXOW9Xa9vRk4gKMIJBEU4gaAIJxAU4QSCIpxAUIQTCIpwAkERTiAowgkERTiBoAgnEBThBIKqnpXiLHo756Qkxy8zZw8bFbMRvHMfur3dnTXx3PlF13kvX1/PHl+6tuyqN761WWybGBmtX9D5oqo+LfasK0ZOICjCCQRFOIGgCCcQFOEEgiKcQFCVt1KWVtZcRQdb7UKLb23/Vqu8aFXTc1W757yEXtG2vLrqqrnZyy+sNdBqueotLq24zlvdzPfjRuH4Ti4vXCi2ra3lb9tU+aU33+XqhwrPr/S/uCVSOK+xy2uJMXICQRFOICjCCQRFOIGgCCcQFOEEgiKcQFCEEwiKcAJBEU4gKMIJBEU4gaAqP/j+7z98xlW03cyXnZgcctU7eexIsW3IszN08r0n9RrlD0pvbG24ajYKXem0yjtDVzk5d9B1Xnsov3P02pbv09xdKz/Haxv1Jwmce+GSqx+pt1Vsmzkw6ao5NTpc+GWuckWMnEBQhBMIinACQRFOICjCCQRFOIGgCCcQFOEEgiKcQFCEEwiKcAJBEU4gKMIJBFU5K+Xbj37XVbS0LP3s9Lir3p2//t5yY6/+bJBO2zc7Zqti9f6hsTFXTStszd1zbhUwOu6bzdIrvE+vXF9y1fv+uXPFtosXyls1lNx42bfD9sry9WLb6dNvcdV885tenz3u22ykjJETCIpwAkERTiAowgkERTiBoAgnEJSligWyvv2f51xLFg218ndoRjqDnnLqtMq7E1++uFC7Xqvju5Wy1i1fLH/qRz921WwO5HcBHx329XF6asJ13tH5o9njXc8CapKur5QX8Wo47jn01nw3KhYWniu23fXGU66a4+P5v81Q27cb+cz4ePa+GSMnEBThBIIinEBQhBMIinACQRFOICjCCQRFOIGgCCcQFOEEgiKcQFCEEwiKcAJBVS7w9ZWvfdNV1Lr5GQkHxqZc9U6+8USxbXZ6tHa90alZVz+e/I8ni21/8dd/56qZCquhjbSbrnrvPP1W13kn5o9lj4+N+mYSTbbKC40t36i/WNdaqr9VvSTd9Ybya+fCSxddNc+fzx8/On/EVW9mPL/wHSMnEBThBIIinEBQhBMIinACQVVerX3ssSdcRacm82usnHqDb4uBx5/4t2Lbqbn6V8gGOvl1e3YyMTFSbHvwgXe7ag4V1lU6cfxOV7122nKd9/SZH2aPX7h42VXvufMvFdtWbtyoXS85NzsYrdgmY2hw2FVzeCJf0wZ9r6tT84eyxxk5gaAIJxAU4QSCIpxAUIQTCIpwAkERTiAowgkERTiBoAgnEBThBIIinEBQhBMIqnJn67n5e1zbGn/qT/4oe/zw3B2ecnrm7NPFtvnpA7XrHX7NcVc/xqYmi21dK+++XWVgMD8r5crL11z1nvqBd4ftfD8aTd9Mi0bF+36vV/9ltbW14eqHVQw/vd6mq+bK9fysmrX1FVe99z/4Xna2Bv4vIZxAUIQTCIpwAkERTiAowgkERTiBoAgnEBThBIIinEBQhBMIinACQRFOIKjKvVLe/q5fcxU9eOxk9vj6um8WwPE3nS62FTaGrrS46ppso0vr68W2tXXfLJKtXn7/mIWXCtsn72Td99g2u/n+P/v8C65603eUdw+/fq3+c/XCwoKrH6by3jGdjm/38NnRiezxwSHfDJ4SRk4gKMIJBEU4gaAIJxAU4QSCqrxaO3XwoKvoV776SPZ4z/te0ChfBTO1apdrtXxX6ZoD5fNabd+VupTyNVeWr7rqWcO3llG3mz/v+tKqq975F8trGa2v19/Zulfo307anfJrrrnl2y173fLPyeCgb+f2EkZOICjCCQRFOIGgCCcQFOEEgiKcQFCEEwiKcAJBEU4gKMIJBEU4gaAIJxAU4QSCqpyV8vCffcpV1Hr5zA8N59de2UkzVcz4aFY+hHw/Rn0zSEZHRoptg2MzrprDI4fy9Zr1Z9tI0sCwb2bEZmF9pNUb9WeQSNJmxa7RDau/zlGj7Xs+JsY7xbah6SFXzXYj//puyTfbqYSREwiKcAJBEU4gKMIJBEU4gaAIJxBU5X2ImYn6tykk6erVK9njh+YOu+oNbJYvvZ+/eKF2vfPPXHL1Y3OzfHug1fHdJjp85z3Z45dX11z15Lyc31D+Fkyn5buFMTJWvl118MB07Xp3zx9x9WN2fLDYtrl13VWz08q/DqanDrjqlTByAkERTiAowgkERTiBoAgnEBThBIIinEBQhBMIinACQRFOICjCCQRFOIGgKj/Zvn5t0VV0dTW/7swLZ8+46h2fKH+geDzV33n5mm24+mHKr7MjSd2N8ofiq5z++fxkgPvv+UVXvSuXfX+zXqH/3d6Wq163W36ubqyv1K730tJ/ufqx2hgrto0fusNV88Dh/HpRh+bnXPVKGDmBoAgnEBThBIIinEBQhBMIinACQRFOICjCCQRFOIGgCCcQFOEEgiKcQFCEEwiqclbK3Gt9S+AvPp6ffbK8cdVV7+xSedn8ptV/f9ns+WalKFXsyGw9V8lH/+HL2eMPnDnnqvd6822hca2Zf2xnnFsWfG9tudjWPTxVu97xt9zr6sfAcHk38qHxcluVidn8rJRmu7z1gwcjJxAU4QSCIpxAUIQTCIpwAkERTiCoyuvus4d8O1E3mz/IHk+9ilsRFVbVLbZZ1e2NAt9Njx3kN4be0Y2V/GJXL16+6Kr3I8fiWZL0RDv/HG/OH3LVm7z7dcW26YP1X1cDgx1XP9Qsjz8Dbd+u3Y1GvmZzwHcbq/h7drUagF1DOIGgCCcQFOEEgiKcQFCEEwiKcAJBEU4gKMIJBEU4gaAIJxAU4QSCIpxAUJUfo9/q+bZSnzgwlD2+urzmqrdZMeWj65jp0vA9LFVNPelZ01XxamGOzOdWLrjq9e4Yd513+OTx7PG75o+56o3MzBbbrqn+38w6vhkkYxMTxbahkWFXzVa7nT3ebPpeAyWMnEBQhBMIinACQRFOICjCCQRVebV2wLm6/NHjB7LHN1a2XPU2Nsur/nRT/cV7msn3ntSo2PrBWr4rdYPt/HkT0/W3LJCknzs07zpvcjJ/5bI96nsRrHbLl8S7jrV2vOv9WMUV1N1eQ2i3MXICQRFOICjCCQRFOIGgCCcQFOEEgiKcQFCEEwiKcAJBEU4gKMIJBEU4gaAIJxCUJcfO0AD2HiMnEBThBIIinEBQhBMIinACQRFOIKifAvCg83TBXtKoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "i = 1\n",
    "\n",
    "n_test = X_test.shape[0]\n",
    "plt.title('Original')\n",
    "plt.axis('off')\n",
    "plt.imshow(X_test[i])\n",
    "plt.show()\n",
    "for _ in range(len(corruption)):\n",
    "    plt.title(corruption[_])\n",
    "    plt.axis('off')\n",
    "    plt.imshow(X_corr[n_test * _+ i])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also verify that the performance of a classification model on CIFAR-10 drops significantly on this perturbed dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test set accuracy:\n",
      "Original 0.9278\n",
      "gaussian_noise 0.2208\n",
      "motion_blur 0.6339\n",
      "brightness 0.8913\n",
      "pixelate 0.3666\n"
     ]
    }
   ],
   "source": [
    "dataset = 'cifar10'\n",
    "model = 'resnet32'\n",
    "clf = fetch_tf_model(dataset, model)\n",
    "acc = clf.evaluate(scale_by_instance(X_test), y_test, batch_size=128, verbose=0)[1]\n",
    "print('Test set accuracy:')\n",
    "print('Original {:.4f}'.format(acc))\n",
    "clf_accuracy = {'original': acc}\n",
    "for _ in range(len(corruption)):\n",
    "    acc = clf.evaluate(scale_by_instance(X_c[_]), y_test, batch_size=128, verbose=0)[1]\n",
    "    clf_accuracy[corruption[_]] = acc\n",
    "    print('{} {:.4f}'.format(corruption[_], acc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the drop in performance, it is important that we detect the harmful data drift!\n",
    "\n",
    "### Detect drift\n",
    "\n",
    "First we try a drift detector using the **TensorFlow** framework for the preprocessing step. We are trying to detect data drift on high-dimensional (*32x32x3*) data using feature-wise univariate tests. It therefore makes sense to apply dimensionality reduction first. Some dimensionality reduction methods also used in [Failing Loudly: An Empirical Study of Methods for Detecting Dataset Shift](https://arxiv.org/pdf/1810.11953.pdf) are readily available: a randomly initialized encoder (**UAE** or Untrained AutoEncoder in the paper), **BBSDs** (black-box shift detection using the classifier's softmax outputs) and **PCA**.\n",
    "\n",
    "#### Random encoder\n",
    "\n",
    "First we try the randomly initialized encoder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:Compiled the loaded model, but the compiled metrics have yet to be built. `model.compile_metrics` will be empty until you train or evaluate the model.\n",
      "WARNING:tensorflow:No training configuration found in the save file, so the model was *not* compiled. Compile it manually.\n"
     ]
    }
   ],
   "source": [
    "from functools import partial\n",
    "from tensorflow.keras.layers import Conv2D, Dense, Flatten, InputLayer, Reshape\n",
    "from alibi_detect.cd.tensorflow import preprocess_drift\n",
    "\n",
    "tf.random.set_seed(0)\n",
    "\n",
    "# define encoder\n",
    "encoding_dim = 32\n",
    "encoder_net = tf.keras.Sequential(\n",
    "  [\n",
    "      InputLayer(input_shape=(32, 32, 3)),\n",
    "      Conv2D(64, 4, strides=2, padding='same', activation=tf.nn.relu),\n",
    "      Conv2D(128, 4, strides=2, padding='same', activation=tf.nn.relu),\n",
    "      Conv2D(512, 4, strides=2, padding='same', activation=tf.nn.relu),\n",
    "      Flatten(),\n",
    "      Dense(encoding_dim,)\n",
    "  ]\n",
    ")\n",
    "\n",
    "# define preprocessing function\n",
    "preprocess_fn = partial(preprocess_drift, model=encoder_net, batch_size=512)\n",
    "\n",
    "# initialise drift detector\n",
    "p_val = .05\n",
    "cd = KSDrift(X_ref, p_val=p_val, preprocess_fn=preprocess_fn)\n",
    "\n",
    "# we can also save/load an initialised detector\n",
    "filepath = 'my_path'  # change to directory where detector is saved\n",
    "save_detector(cd, filepath)\n",
    "cd = load_detector(filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value used by the detector for the multivariate data with *encoding_dim* features is equal to *p_val / encoding_dim* because of the [Bonferroni correction](https://mathworld.wolfram.com/BonferroniCorrection.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert cd.p_val / cd.n_features == p_val / encoding_dim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check whether the detector thinks drift occurred on the different test sets and time the prediction calls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from timeit import default_timer as timer\n",
    "\n",
    "labels = ['No!', 'Yes!']\n",
    "\n",
    "def make_predictions(cd, x_h0, x_corr, corruption):\n",
    "    t = timer()\n",
    "    preds = cd.predict(x_h0)\n",
    "    dt = timer() - t\n",
    "    print('No corruption')\n",
    "    print('Drift? {}'.format(labels[preds['data']['is_drift']]))\n",
    "    print('Feature-wise p-values:')\n",
    "    print(preds['data']['p_val'])\n",
    "    print(f'Time (s) {dt:.3f}')\n",
    "    \n",
    "    if isinstance(x_corr, list):\n",
    "        for x, c in zip(x_corr, corruption):\n",
    "            t = timer()\n",
    "            preds = cd.predict(x)\n",
    "            dt = timer() - t\n",
    "            print('')\n",
    "            print(f'Corruption type: {c}')\n",
    "            print('Drift? {}'.format(labels[preds['data']['is_drift']]))\n",
    "            print('Feature-wise p-values:')\n",
    "            print(preds['data']['p_val'])\n",
    "            print(f'Time (s) {dt:.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No corruption\n",
      "Drift? No!\n",
      "Feature-wise p-values:\n",
      "[0.9386024  0.13979132 0.6384489  0.05413922 0.37460664 0.25598603\n",
      " 0.87304014 0.47553554 0.11587767 0.67217577 0.47553554 0.7388285\n",
      " 0.08215971 0.14635575 0.3114053  0.3114053  0.60482025 0.36134896\n",
      " 0.8023182  0.21715216 0.24582714 0.46030036 0.11587767 0.44532147\n",
      " 0.25598603 0.58811766 0.5550683  0.95480835 0.8598946  0.23597081\n",
      " 0.8975547  0.68899393]\n",
      "Time (s) 1.172\n",
      "\n",
      "Corruption type: gaussian_noise\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[4.85834153e-03 7.20506581e-03 5.44517934e-02 9.87569049e-09\n",
      " 3.35018486e-01 8.05620551e-02 6.66609779e-03 2.68237293e-01\n",
      " 1.52247362e-02 1.01558706e-02 1.78680534e-03 1.04267694e-01\n",
      " 4.93385670e-08 1.35106135e-10 1.04696119e-04 1.35730659e-06\n",
      " 2.87180692e-01 3.79266362e-06 3.45018925e-04 1.96636513e-01\n",
      " 1.86571106e-03 5.92635339e-03 4.70917694e-10 5.92635339e-03\n",
      " 5.07743537e-01 5.31427140e-05 3.80059540e-01 1.13354892e-01\n",
      " 2.75738519e-02 7.75579622e-07 3.23252240e-03 2.02312917e-02]\n",
      "Time (s) 2.516\n",
      "\n",
      "Corruption type: motion_blur\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[3.39037769e-07 1.16525307e-01 5.04726835e-04 3.81079665e-03\n",
      " 6.31192625e-01 5.28989534e-04 3.61990853e-04 1.57829020e-02\n",
      " 1.94784126e-03 1.26909809e-02 4.46249526e-09 6.99155149e-04\n",
      " 3.79746925e-04 5.88651128e-21 1.35596551e-07 2.00218983e-05\n",
      " 7.15865940e-02 7.28750820e-05 1.04267694e-01 1.10198918e-04\n",
      " 2.22608112e-04 1.52403876e-01 6.41064299e-03 3.15323919e-02\n",
      " 3.04985344e-02 8.97102946e-05 6.54255822e-02 2.03331537e-03\n",
      " 1.15137536e-03 8.04463718e-10 9.62164486e-04 3.45018925e-04]\n",
      "Time (s) 2.622\n",
      "\n",
      "Corruption type: brightness\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[0.0000000e+00 0.0000000e+00 4.0479114e-29 0.0000000e+00 0.0000000e+00\n",
      " 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00\n",
      " 0.0000000e+00 0.0000000e+00 0.0000000e+00 2.3823024e-21 2.9582986e-38\n",
      " 0.0000000e+00 0.0000000e+00 0.0000000e+00 7.1651735e-34 0.0000000e+00\n",
      " 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00\n",
      " 3.8567345e-05 0.0000000e+00 0.0000000e+00 0.0000000e+00 0.0000000e+00\n",
      " 0.0000000e+00 0.0000000e+00]\n",
      "Time (s) 2.572\n",
      "\n",
      "Corruption type: pixelate\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[3.35018486e-01 2.06576437e-01 1.96636513e-01 8.40378553e-03\n",
      " 5.92454255e-01 9.30991620e-02 1.64748102e-01 5.16919196e-01\n",
      " 2.40608733e-02 4.54285979e-01 1.91894709e-04 1.33487985e-01\n",
      " 1.63820793e-03 1.78680534e-03 1.13354892e-01 9.04688612e-02\n",
      " 8.27570856e-01 6.15745559e-02 6.54255822e-02 9.06871445e-03\n",
      " 2.38713458e-01 6.89552963e-01 1.07227206e-01 8.29487666e-02\n",
      " 3.42268527e-01 1.37110472e-01 3.64637136e-01 3.00327957e-01\n",
      " 3.72297794e-01 9.06871445e-03 4.98639137e-01 9.78103094e-03]\n",
      "Time (s) 2.554\n"
     ]
    }
   ],
   "source": [
    "make_predictions(cd, X_h0, X_c, corruption)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, drift was only detected on the corrupted datasets. The feature-wise p-values for each univariate K-S test per (encoded) feature before multivariate correction show that most of them are well above the $0.05$ threshold for *H0* and below for the corrupted datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### BBSDs\n",
    "\n",
    "For **BBSDs**, we use the classifier's softmax outputs for black-box shift detection. This method is based on [Detecting and Correcting for Label Shift with Black Box Predictors](https://arxiv.org/abs/1802.03916). The ResNet classifier is trained on data standardised by instance so we need to rescale the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = scale_by_instance(X_train)\n",
    "X_test = scale_by_instance(X_test)\n",
    "X_ref = scale_by_instance(X_ref)\n",
    "X_h0 = scale_by_instance(X_h0)\n",
    "X_c = [scale_by_instance(X_c[i]) for i in range(n_corr)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we initialize the detector. Here we use the output of the softmax layer to detect the drift, but other hidden layers can be extracted as well by setting *'layer'* to the index of the desired hidden layer in the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from alibi_detect.cd.tensorflow import HiddenOutput\n",
    "\n",
    "# define preprocessing function, we use the\n",
    "preprocess_fn = partial(preprocess_drift, model=HiddenOutput(clf, layer=-1), batch_size=128)\n",
    "\n",
    "cd = KSDrift(X_ref, p_val=p_val, preprocess_fn=preprocess_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we can see that the p-value used by the detector for the multivariate data with 10 features (number of CIFAR-10 classes) is equal to *p_val / 10* because of the [Bonferroni correction](https://mathworld.wolfram.com/BonferroniCorrection.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert cd.p_val / cd.n_features == p_val / 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no drift on the original held out test set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No corruption\n",
      "Drift? No!\n",
      "Feature-wise p-values:\n",
      "[0.11587767 0.5226477  0.19109942 0.19949944 0.49101472 0.722359\n",
      " 0.12151605 0.41617486 0.8320209  0.75510186]\n",
      "Time (s) 7.608\n",
      "\n",
      "Corruption type: gaussian_noise\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "Time (s) 16.233\n",
      "\n",
      "Corruption type: motion_blur\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "Time (s) 15.688\n",
      "\n",
      "Corruption type: brightness\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[0.0000000e+00 3.8790170e-15 2.2549014e-33 4.6733894e-07 2.1857751e-15\n",
      " 1.2091652e-05 2.3977423e-30 1.0099583e-09 4.3286997e-12 3.8117909e-17]\n",
      "Time (s) 14.733\n",
      "\n",
      "Corruption type: pixelate\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "Time (s) 13.781\n"
     ]
    }
   ],
   "source": [
    "make_predictions(cd, X_h0, X_c, corruption)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Label drift\n",
    "\n",
    "We can also check what happens when we introduce class imbalances between the reference data *X_ref* and the tested data *X_imb*. The reference data will use $75$% of the instances of the first 5 classes and only $25$% of the last 5. The data used for drift testing then uses respectively $25$% and $75$% of the test instances for the first and last 5 classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5000, 32, 32, 3) (5000, 32, 32, 3) (5000,)\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(0)\n",
    "# get index for each class in the test set\n",
    "num_classes = len(np.unique(y_test))\n",
    "idx_by_class = [np.where(y_test == c)[0] for c in range(num_classes)]\n",
    "# sample imbalanced data for different classes for X_ref and X_imb\n",
    "perc_ref = .75\n",
    "perc_ref_by_class = [perc_ref if c < 5 else 1 - perc_ref for c in range(num_classes)]\n",
    "n_by_class = n_test // num_classes\n",
    "X_ref = []\n",
    "X_imb, y_imb = [], []\n",
    "for _ in range(num_classes):\n",
    "    idx_class_ref = np.random.choice(n_by_class, size=int(perc_ref_by_class[_] * n_by_class), replace=False)\n",
    "    idx_ref = idx_by_class[_][idx_class_ref]\n",
    "    idx_class_imb = np.delete(np.arange(n_by_class), idx_class_ref, axis=0)\n",
    "    idx_imb = idx_by_class[_][idx_class_imb]\n",
    "    assert not np.array_equal(idx_ref, idx_imb)\n",
    "    X_ref.append(X_test[idx_ref])\n",
    "    X_imb.append(X_test[idx_imb])\n",
    "    y_imb.append(y_test[idx_imb])\n",
    "X_ref = np.concatenate(X_ref)\n",
    "X_imb = np.concatenate(X_imb)\n",
    "y_imb = np.concatenate(y_imb)\n",
    "print(X_ref.shape, X_imb.shape, y_imb.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Update reference dataset for the detector and make predictions. Note that we store the preprocessed reference data since the `preprocess_at_init` kwarg is by default True:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd.x_ref = cd.preprocess_fn(X_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? Yes!\n",
      "[5.2598997e-20 1.1312397e-20 5.8646589e-29 9.2977640e-18 6.4071548e-23\n",
      " 1.4155961e-15 5.0236095e-19 1.6651963e-20 4.2726706e-21 3.5123729e-21]\n"
     ]
    }
   ],
   "source": [
    "preds_imb = cd.predict(X_imb)\n",
    "print('Drift? {}'.format(labels[preds_imb['data']['is_drift']]))\n",
    "print(preds_imb['data']['p_val'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update reference data\n",
    "\n",
    "So far we have kept the reference data the same throughout the experiments. It is possible however that we want to test a new batch against the last *N* instances or against a batch of instances of fixed size where we give each instance we have seen up until now the same chance of being in the reference batch ([reservoir sampling](https://en.wikipedia.org/wiki/Reservoir_sampling)). The `update_x_ref` argument allows you to change the reference data update rule. It is a Dict which takes as key the update rule (*'last'* for last *N* instances or *'reservoir_sampling'*) and as value the batch size *N* of the reference data. You can also save the detector after the prediction calls to save the updated reference data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 7500\n",
    "cd = KSDrift(X_ref, p_val=.05, preprocess_fn=preprocess_fn, update_x_ref={'reservoir_sampling': N})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reference data is now updated with each `predict` call. Say we start with our imbalanced reference set and make a prediction on the remaining test set data *X_imb*, then the drift detector will figure out data drift has occurred."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? Yes!\n"
     ]
    }
   ],
   "source": [
    "preds_imb = cd.predict(X_imb)\n",
    "print('Drift? {}'.format(labels[preds_imb['data']['is_drift']]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now see that the reference data consists of *N* instances, obtained through reservoir sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert cd.x_ref.shape[0] == N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then draw a random sample from the training set and compare it with the updated reference data. This still highlights that there is data drift but will update the reference data again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "perc_train = .5\n",
    "n_train = X_train.shape[0]\n",
    "idx_train = np.random.choice(n_train, size=int(perc_train * n_train), replace=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? Yes!\n"
     ]
    }
   ],
   "source": [
    "preds_train = cd.predict(X_train[idx_train])\n",
    "print('Drift? {}'.format(labels[preds_train['data']['is_drift']]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we draw a new sample from the training set, it highlights that it is not drifting anymore against the reservoir in *X_ref*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? No!\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(1)\n",
    "perc_train = .1\n",
    "idx_train = np.random.choice(n_train, size=int(perc_train * n_train), replace=False)\n",
    "preds_train = cd.predict(X_train[idx_train])\n",
    "print('Drift? {}'.format(labels[preds_train['data']['is_drift']]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate correction mechanism\n",
    "\n",
    "Instead of the Bonferroni correction for multivariate data, we can also use the less conservative [False Discovery Rate](http://www.math.tau.ac.il/~ybenja/MyPapers/benjamini_hochberg1995.pdf) (FDR) correction. See [here](https://riffyn.com/riffyn-blog/2017/10/29/false-discovery-rate) or [here](https://matthew-brett.github.io/teaching/fdr.html) for nice explanations. While the Bonferroni correction controls the probability of at least one false positive, the FDR correction controls for an expected amount of false positives. The `p_val` argument at initialisation time can be interpreted as the acceptable q-value when the FDR correction is applied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? Yes!\n"
     ]
    }
   ],
   "source": [
    "cd = KSDrift(X_ref, p_val=.05, preprocess_fn=preprocess_fn, correction='fdr')\n",
    "\n",
    "preds_imb = cd.predict(X_imb)\n",
    "print('Drift? {}'.format(labels[preds_imb['data']['is_drift']]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adversarial autoencoder as a malicious drift detector\n",
    "\n",
    "We can leverage the adversarial scores obtained from an [adversarial autoencoder](https://arxiv.org/abs/2002.09364)  trained on normal data and transform it into a data drift detector. The score function of the adversarial autoencoder becomes the preprocessing function for the drift detector. The K-S test is then a simple univariate test on the adversarial scores. Importantly, an adversarial drift detector flags **malicious data drift**. We can fetch the pretrained adversarial detector from a [Google Cloud Bucket](https://console.cloud.google.com/storage/browser/seldon-models/alibi-detect/ad/cifar10/resnet32) or train one from scratch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_pretrained = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from tensorflow.keras.regularizers import l1\n",
    "from tensorflow.keras.layers import Conv2DTranspose\n",
    "from alibi_detect.ad import AdversarialAE\n",
    "\n",
    "# change filepath to (absolute) directory where model is downloaded\n",
    "filepath = os.path.join(os.getcwd(), 'my_path')\n",
    "detector_type = 'adversarial'\n",
    "detector_name = 'base'\n",
    "filepath = os.path.join(filepath, detector_name)\n",
    "if load_pretrained:\n",
    "    ad = fetch_detector(filepath, detector_type, dataset, detector_name, model=model)\n",
    "else:  # train detector from scratch\n",
    "    # define encoder and decoder networks\n",
    "    encoder_net = tf.keras.Sequential(\n",
    "            [\n",
    "                InputLayer(input_shape=(32, 32, 3)),\n",
    "                Conv2D(32, 4, strides=2, padding='same', \n",
    "                       activation=tf.nn.relu, kernel_regularizer=l1(1e-5)),\n",
    "                Conv2D(64, 4, strides=2, padding='same', \n",
    "                       activation=tf.nn.relu, kernel_regularizer=l1(1e-5)),\n",
    "                Conv2D(256, 4, strides=2, padding='same', \n",
    "                       activation=tf.nn.relu, kernel_regularizer=l1(1e-5)),\n",
    "                Flatten(),\n",
    "                Dense(40)\n",
    "            ]\n",
    "        )\n",
    "    \n",
    "    decoder_net = tf.keras.Sequential(\n",
    "        [\n",
    "                InputLayer(input_shape=(40,)),\n",
    "                Dense(4 * 4 * 128, activation=tf.nn.relu),\n",
    "                Reshape(target_shape=(4, 4, 128)),\n",
    "                Conv2DTranspose(256, 4, strides=2, padding='same', \n",
    "                                activation=tf.nn.relu, kernel_regularizer=l1(1e-5)),\n",
    "                Conv2DTranspose(64, 4, strides=2, padding='same', \n",
    "                                activation=tf.nn.relu, kernel_regularizer=l1(1e-5)),\n",
    "                Conv2DTranspose(3, 4, strides=2, padding='same', \n",
    "                                activation=None, kernel_regularizer=l1(1e-5))\n",
    "            ]\n",
    "        )\n",
    "    \n",
    "    # initialise and train detector\n",
    "    ad = AdversarialAE(encoder_net=encoder_net, decoder_net=decoder_net, model=clf)\n",
    "    ad.fit(X_train, epochs=50, batch_size=128, verbose=True)\n",
    "    \n",
    "    # save the trained adversarial detector\n",
    "    save_detector(ad, filepath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialise the drift detector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "idx = np.random.choice(n_test, size=n_test // 2, replace=False)\n",
    "X_ref = scale_by_instance(X_test[idx])\n",
    "\n",
    "# adversarial score fn = preprocess step\n",
    "preprocess_fn = partial(ad.score, batch_size=128)\n",
    "\n",
    "cd = KSDrift(X_ref, p_val=.05, preprocess_fn=preprocess_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make drift predictions on the original test set and corrupted data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H0: Accuracy 0.9286 -- Drift? No!\n",
      "imbalance: Accuracy 0.9282 -- Drift? No!\n",
      "gaussian_noise: Accuracy 0.2208 -- Drift? Yes!\n",
      "motion_blur: Accuracy 0.6339 -- Drift? Yes!\n",
      "brightness: Accuracy 0.8913 -- Drift? Yes!\n",
      "pixelate: Accuracy 0.3666 -- Drift? Yes!\n"
     ]
    }
   ],
   "source": [
    "clf_accuracy['h0'] = clf.evaluate(X_h0, y_h0, batch_size=128, verbose=0)[1]\n",
    "preds_h0 = cd.predict(X_h0)\n",
    "print('H0: Accuracy {:.4f} -- Drift? {}'.format(\n",
    "    clf_accuracy['h0'], labels[preds_h0['data']['is_drift']]))\n",
    "clf_accuracy['imb'] = clf.evaluate(X_imb, y_imb, batch_size=128, verbose=0)[1]\n",
    "preds_imb = cd.predict(X_imb)\n",
    "print('imbalance: Accuracy {:.4f} -- Drift? {}'.format(\n",
    "    clf_accuracy['imb'], labels[preds_imb['data']['is_drift']]))\n",
    "for x, c in zip(X_c, corruption):\n",
    "    preds = cd.predict(x)\n",
    "    print('{}: Accuracy {:.4f} -- Drift? {}'.format(\n",
    "        c, clf_accuracy[c],labels[preds['data']['is_drift']]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While *X_imb* clearly exhibits input data drift due to the introduced class imbalances, it is not flagged by the adversarial drift detector since the performance of the classifier is not affected and the drift is not malicious. We can visualise this by plotting the adversarial scores together with the harmfulness of the data corruption as reflected by the drop in classifier accuracy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "adv_scores = {}\n",
    "score = ad.score(X_ref, batch_size=128)\n",
    "adv_scores['original'] = {'mean': score.mean(), 'std': score.std()}\n",
    "score = ad.score(X_h0, batch_size=128)\n",
    "adv_scores['h0'] = {'mean': score.mean(), 'std': score.std()}\n",
    "score = ad.score(X_imb, batch_size=128)\n",
    "adv_scores['imb'] = {'mean': score.mean(), 'std': score.std()}\n",
    "\n",
    "for x, c in zip(X_c, corruption):\n",
    "    score_x = ad.score(x, batch_size=128)\n",
    "    adv_scores[c] = {'mean': score_x.mean(), 'std': score_x.std()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAE8CAYAAACYW1XSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABCH0lEQVR4nO3debxc8/3H8ddbdhFr1JYgtSYhsUS0tqQIsVSsRS2xVPjZFa1qcaWlqBZVbcUWayxBGktbtSQoRRCCBEFICGKXRJDk8/vj853k5LrL3OTeOXfufJ6PxzzunJlzZj73zMz5nO/3fBeZGSGEEEK5WCrvAEIIIYSGiMQVQgihrETiCiGEUFYicYUQQigrkbhCCCGUlUhcIYQQykokrkTScEm/yzuO+kh6WVL/Itc1Ses2bUTlR9IYST/LO468SZoiace848hLQ35LoXlp8YkrHaQ+ldQu71gag5n1NLMxS/o6knpKekDSJ5I+k/SspF0bIcQWQdJhKfHvn3cs5UDSspIulfSOpJmS3kjLnfOODWo+MW2s31K193k5/f8zJc2TNCezfGZjvlcla9GJS9LawLaAAXvk8P6S1Cj7WFLrxnidjHuA/wCrAt8DTgS+aMw3aIKYS2kw8AlwaCnftBz3maS2wENAT2AgsCzwQ+BjoO9ivF6rastls09SMlzGzJYBHgOOLyyb2fl5x9dimFmLvQFnA/8F/gTcW+25TYHngC+B24Bbgd+l5yYCu2fWbQ3MADZLyz8AngA+A14A+mfWHQOcl973K2Bd4DDgzfRebwEHpXXXAR7Gf+AfATcDy2deawrwS+BF4OsUxxRgx/R8X+DJFMd04C9A28z2Bqxbw37pnJ5bvo59NwgYjyezN4CB6fHVgdH4QX0ycFRmmypgJHBT2u5nwHLANSm+d4HfAa3S+usCY4HP0/9/Wx3x3AG8n9Z9FOiZeW44cAVwX9rHTwHrZJ4fAExK2/4lvefP6nivtYD5wD7AXGDVzHOtgDPTPvkSeBbomp7riZ8MfAJ8AJyZie93mdfoD0yr53M+I/MerwB7VYvxKPx7Wnh+M+B04M5q6/0ZuKyW/3MK8Ku0/afAdUD79NxLwI8z67ZJn9GmNbzOz9L/u0wd+7Q7/tv4DHgZ2KPa5/c34H5gFrBjLftkke9zdr8W9mn6bD5K2xd+Z0OAb4FvgJnAPZn/v/BbagdcCryXbpcC7aq99qnAh/h3+fAijj9jSN+zuvYnsHb634ak954OnJZZd6nM9+Fj4HZgxTyOqc3llnsATfrP+YH1WGDz9MVdJT3eFngbOCV9gfZNzxd+BGcDN2deZzdgYrq/Rvry7Jq+UAPS8srp+THAO/hBrDV+4P4C2CA9vxrpoIsfuAekH83K+AH50sz7TsGTR1egQ+axwo9tczyJtk5f/onAyZnta0tcAl4H7gX2LOyXzPN98YP8gPQ/rgFsmJ57FPgr0B7YBE/o26fnqtJ+3DNt1wG4G7gS6IiX7J4Gjk7rjwB+ndZtD2xTx2d5BNCJhQeY8ZnnhrPw7L41fgJwa3quM35w3zd91qfgyaiuxHUW8HS6PwE4NfPc6emxDdJ+7A2slGKbjh/c2qflLTPx1Ze4qn/O++EnCUsB++MH9NUyz70LbJFiWBdPtqul9ZZP67XGD7Sb1/J/TsEPqF2BFfGTrcJv4BdkTiTwE5kJtbzOrcD1dezPNvhv8Uz8t7d9+kw2yOyfz4GtM9+FmvZJfYlrLn6S2g7ol/bFBjV9BjX8loYC/8O/oyvjJ6a/rfbaQ9P/siswG1ihnuPPGBYmrlr3JwsT1wj8d7Ix/rsqxHZSiq1L+t+uBEbkfXzN85Z7AE32j8E2+EG0c1qeBJyS7m+Hn9kos/4TmR/BuumHtXRavhk4O93/JXBjtff6NzA43R8DDM081xE/y9yn8AOsI+Y9geczy1OAI6qts+DHVsP2JwN3Z5ZrTFzpuS546eMNvHTxKLBeeu5K4JIatukKzAM6ZR77PTA83a8CHs08twp+ttwh89iBwCPp/g3AMKBLAz/b5dP/tlxaHg5cnXl+V2BSun8o8L/Mc8LPnutKXK+TTgDwEskLmedeBQbVsM2B2c+u2nPDqT9xHVFbPGmd8YX3Td+3k2pZ75+kUjCwO/BKHa85BTim2n57I91fHf8NLJuWRwK/qOV1/gNcUMf7bIuXlpfKPDYCqMrsnxtqiK36d7+YxNUx8/ztwFk1fQbVf0v472DXzHM7A1Myr/0V0Drz/IfAD+r5zMawMHHVuj9ZmLg2zGx7EXBNuj8R2CHz3Gr4sa11Xe/fkm8t+RrXYOABM/soLd+SHgP/Er1r6VuQvF24Y2aT8S/LjyUtjV8fuyU9vRawX2rQ8Jmkz/AkuVrmtaZmXmsWfsZ8DDBd0n2SNgSQtIqkWyW9K+kLvIqt+sXsqdRC0vqS7pX0ftr+/Bq2r5GZTTOz481snfQ/zcITCXiCeqOGzVYHPjGzLzOPvY2XyGqKdy38DHV6Zl9diZ/Vgp+FCng6XdQ+opb/s5WkC9IF/y/wAw7V/tf3M/dnA8tkYs5+Hkbd+3RroBteigD/3DeWtElarm3f1PZ4sRaJSdKhksZn9ttGLPx/63qv64GD0/2DgRsb8L5v4/sLM3sPL4HtI2l5YBf8BK4mH7Po97+61YGpZja/2nvV9r2p67G6fJp+b9n3WL3IbVcncwyoYduPzWxuZjn7HatXkfuzxs8C/x3dnfkuTMRPIFcp9v1bmhaZuCR1AH4C9EsH9ffxKqLeknrjVTprSFJmszWrvcwI/Cx6EH7WOjk9PhUvcS2fuXU0swsy22YTImb2bzMbgP+4JwFXpafOT+tubGbL4geabEzfea1q/pZeb720/Zk1bF8vM5uKXyPaKPM/rlPDqu8BK0rqlHlsTbzaqqZ4p+Ilrs6ZfbWsmfVM7/u+mR1lZqsDRwN/raX5/k/xz2FHvOp17fR4Mf/rdPxA7xv4Z9619tUZnF53fPrePJV5vPA/1bRvpgLfr+U1ZwFLZ5ZXrWGdBftN0lr4d+R4YCUzWx6v0iv8v7XFADAK6CVpI7zEVVuyKcjuizXxz7igkAT3A540s+znnPUgsLOkjrU8/x7QtVpDpbq+N7U9Npu69+MK1WLI/j91/Y4KMa5Vy7aNpb79WdtnMRXYpdoxp30dn0eL1yITF17lNg/ogV+H2QS/OPwYXnX0JF6tcKKkNpL25rutn24FdgL+j4WlLfBS0Y8l7ZxKAu0l9ZfUpaZAUqlqUPpBfY1fHC6ceXZKy59LWgO/ftIQnfDrZzNTKe7/itlI0gqSzpW0rqSlUpPlI/B6dPDGFIdL2iE9v4akDVOCewL4ffq/ewFHpn3yHWY2HXgA+GNqLr2UpHUk9Utx7JfZb5/iB5f5NbxUJ3zffYwfuBrSOus+oKekvVPrtBOpOXEgqT1+wjOEhd+bTYATgJ+m7a8GfitpvdRqtJeklfDrhatJOllSO0mdJG2ZXno8sKukFSWtilfp1qUjvi9mpLgOZ+FJBSmG0yRtnmJYNyU7zGwOXg11C36d7p163us4SV0krYhfb7wt89wovNHHSSwsjdfkRvzgeqekDdPnvJKkM+VdLJ7Ck84v0u+tP/BjFpZqizUe/xxaSRqIX8eq7lxJbSVtiyfuO9LjH1D7iQX4iepvJK2cfg9nU8v3egmMou79eZakpSX1BA5n4Wfxd+C8wmecYhzUyLGVlZaauAYD15nZO+ms/n0zex+/pnMQfnDcG2/t9wlelXdX9gXSQfdJYCsyP+Z08B6El25m4D/Y06l9Xy4F/Bw/e/oE/7EVEsy5+Bf5c/wAe1dNL1CH0/DSyJf4Gfptda++wDd4qeVBPPG9hCeGwwDM7Gn8h3NJim0sC89GD0zbvoc3vDjHzB6s470OxS/IF1qujWRhtdIWwFOSZuItFU8yszdreI0b8KqTd9Pr/K+GdWqUqor3Ay7AE996eJVNTfbEr2XcUO17cy3e0GEgfvH/djwhf4En+Q6p+nQAfkB+H79O9qP0ujfirU+npO3q/JzM7BXgj/j37wP8Yv1/M8/fgbdcvQX/7EfhjSsKrk/b1FdNSHqNB/BWr2/grT4L7/MVcCdedVrrd9PMvsZLw5Pw611f4I1wOgNPmdk3+H7ZBW9J91fgUDObVER8WSel1/kM/x2Pqvb8+/h37D28pHlM5j2uAXqk6rbq24H/3+PwVowT8BbHjTogQRH7cyzeiOUh4GIzeyA9fhn++3hA0pf493/LGravGFr0Mk8IodxJWhNPIqua2RL1zZN0NrC+mR1c78o5SqW4m8ysxpqP5qKm/Snvb/oW0KbadbRQi7Lp2BdCqF+6jvRzvDvAkiatFfGq4EMaI7ZKF/uz8bTUqsIQKk66jvoFXmV5zhK+1lF4Nfg/zezRRgivojX3/SnpWkkfSnqplucPkvSipAmSnkiN3HITVYUhhFDhJG2HNxS7wcw2quH5rfBBGD6VtAveBy+362xRVRhCCBXOzB5N19pqe/6JzGJhFI/cRFVhCCGEhjgSH6ElN2VV4lpqqaWsQ4cOeYcRQghlZfbs2YY38S8YZmbDGvo6kn6EJ65tGiu2xVFWiatDhw7MmjWr/hVDCCEsIOkrM+uzhK/RC+/8vouZfdw4kS2eqCoMIYRQp9Q38C7gEDN7Le94yqrEFUIIofFJGoGPgt9Z0jS8O0UbADP7Oz4E1kr4eKIAc5e0BLckyqo5fMeOHS2qCkMIoWEkzTaz2gZBLjtR4gqhiXz77bdMmzaNOXPm5B1KqBDt27enS5cutGnTJu9QmlQkrhCayLRp0+jUqRNrr702i86gE0LjMzM+/vhjpk2bRrdu3fIOp0lF44wQmsicOXNYaaWVImmFkpDESiutVBEl/EhcITShSFqhlCrl+xaJK4QW7rzzzqNnz5706tWLTTbZhKeeeqr+jSrYqFGjkMSkSQ2dLiyUSlzjCqFE1j7jvkZ9vSkX7FbvOk8++ST33nsvzz33HO3ateOjjz7im2++Wez3nDt3Lq1bl+awMXHD7o36et0nTSxqvREjRrDNNtswYsQIzj333EaNoWDevHm0atWqSV67EkSJK4QWbPr06XTu3Jl27doB0LlzZ1ZffXWeeeYZttpqK3r37k3fvn358ssvmTNnDocffjgbb7wxm266KY888ggAw4cPZ4899mD77bdnhx12YNasWRxxxBH07duXTTfdlH/84x95/ouNaubMmTz++ONcc8013HrrrYAnmdNOO42NNtqIXr16cfnllwPUuA+HDx/O8ccfv+D1dt99d8aMGQPAMsssw6mnnkrv3r158sknGTp0KFtssQUbbbQRQ4YModA1afLkyey444707t2bzTbbjDfeeINDDz2UUaNGLXjdgw46qEXt94aKxBVCHaqqqpC04FZVVZV3SA2y0047MXXqVNZff32OPfZYxo4dyzfffMP+++/PZZddxgsvvMCDDz5Ihw4duOKKK5DEhAkTGDFiBIMHD15wof+5555j5MiRjB07lvPOO4/tt9+ep59+mkceeYTTTz+9xQzF9o9//IOBAwey/vrrs9JKK/Hss88ybNgwpkyZwvjx43nxxRc56KCDat2HdZk1axZbbrklL7zwAttssw3HH388zzzzDC+99BJfffUV9957L+BJ6bjjjuOFF17giSeeYLXVVuPII49k+PDhAHz++ec88cQT7LZb/SXuliqqCkOoQ1VV1YIz5sLfcrLMMsvw7LPP8thjj/HII4+w//778+tf/5rVVluNLbbYAoBll10WgMcff5wTTjgBgA033JC11lqL117z0X0GDBjAiiuuCMADDzzA6NGjufjiiwFvPfnOO+/QvXvjVu3lYcSIEZx00kkAHHDAAYwYMYK33nqLY445ZkEV6YorrsiECRNq3Id1adWqFfvss8+C5UceeYSLLrqI2bNn88knn9CzZ0/69+/Pu+++y1577QV4vyyAfv36ceyxxzJjxgzuvPNO9tlnn5JV2TZHlfufh1AhWrVqRf/+/enfvz8bb7wxV1xxRYNfo2PHhYMumBl33nknG2ywQWOGmbtPPvmEhx9+mAkTJiCJefPmIWlBcipG69atmT9//oLlbNP09u3bL7iuNWfOHI499ljGjRtH165dqaqqqrcZ+6GHHspNN93ErbfeynXXXdfA/65liarCEFqwV199lddff33B8vjx4+nevTvTp0/nmWeeAeDLL79k7ty5bLvtttx8880AvPbaa7zzzjs1Jqedd96Zyy+/fME1meeff74E/0nTGzlyJIcccghvv/02U6ZMYerUqXTr1o3evXtz5ZVXMnfuXMAT3AYbbFDjPlx77bUZP3488+fPZ+rUqTz99NM1vlchSXXu3JmZM2cycuRIADp16kSXLl0WXM/6+uuvmT17NgCHHXYYl156KQA9evRoqt1QFnJNXJJOkfSypJckjZDUPs94QmhpZs6cyeDBg+nRowe9evXilVdeYejQodx2222ccMIJ9O7dmwEDBiwoAcyfP5+NN96Y/fffn+HDhy9o1JF11lln8e2339KrVy969uzJWWedlcN/1vhGjBixoIquYJ999mH69Omsueaa9OrVi969e3PLLbfQtm3bGvfh1ltvTbdu3ejRowcnnngim222WY3vtfzyy3PUUUex0UYbsfPOOy9Sqrvxxhv585//TK9evdhqq614//33AVhllVXo3r07hx9+eNPthDKR2yC7ktYAHgd6mNlXkm4H7jez4bVtE4Pshjz0798faPg1rokTJ7aI6z6heZg9ezYbb7wxzz33HMstt1yt69X0vWtpg+zmXVXYGuggqTWwNPBezvGEEEKz8+CDD9K9e3dOOOGEOpNWpcitcYaZvSvpYuAd4CvgATN7oPp6koYAQwDatm1b2iBDCKEZ2HHHHXn77bfzDqPZyK3EJWkFYBDQDVgd6Cjp4OrrmdkwM+tjZn0quflnCCEEl2dV4Y7AW2Y2w8y+xaeF3irHeEJodOU0UWsof5Xyfcszcb0D/EDS0vIhjXcAihtMLIQy0L59ez7++OOKOZiEfBXm4yp0Wm7J8rzG9ZSkkcBzwFzgeWBYXvGE0Ni6dOnCtGnTmDFjRt6hhApRmAG5pcv1opGZnQOck2cMITSVNm3atPiZaEPIQ97N4UMIIYQGicQVQggVTtK1kj6U9FItz0vSnyVNlvSipJqHBCmRSFwhhBCGAwPreH4XYL10GwL8rQQx1SoSVwghLIFyn7MNwMweBT6pY5VBwA3m/gcsL2m10kT3XdGjN4QQlkCZzNnWWtK4zPIwM2tIK+41gKmZ5WnpsemNEVxDReIKIYSWb66Z9ck7iMYSVYUhhBDq8y7QNbPcJT2Wi0hcIYQQ6jMaODS1LvwB8LmZ5VJNCFFVGEIIFU/SCKA/0FnSNHxgiDYAZvZ34H5gV2AyMBvIdTbLSFwhhFDhzOzAep434LgShVOvqCoMIYRQViJxhRBCKCuRuEIIIZSVSFwhhBDKSiSuEEIIZSVaFRZp4obdl2j77pNKP7lzOcYM5Rt3CKE0osQVQgihrETiCiGEUFYicYUQQigrkbhCCCGUlUhcIYQQykq0KgwhlJVodRqixBVCCKGsROIKIYRQViJxhRBCKCuRuEIIIZSVSFwhhBDKSiSuEEIIZSUSVwghhLISiSuEEEJZicQVQgihrOSauCQtL2mkpEmSJkr6YZ7xhBBCJZI0UNKrkiZLOqOG59eU9Iik5yW9KGnXPOIsyLvEdRnwLzPbEOgNxFgsIYRQQpJaAVcAuwA9gAMl9ai22m+A281sU+AA4K+ljXJRuSUuScsB2wHXAJjZN2b2WV7xhBBCheoLTDazN83sG+BWYFC1dQxYNt1fDnivhPF9R54lrm7ADOC6VPy8WlLH6itJGiJpnKRxc+fOLX2UIYRQ/loXjqPpNiTz3BrA1MzytPRYVhVwsKRpwP3ACU0abT3yTFytgc2Av6Xi5yzgO3WrZjbMzPqYWZ/WrWMw+xBCWAxzC8fRdBvWwO0PBIabWRdgV+BGSbnljzwT1zRgmpk9lZZH4okshBBC6bwLdM0sd0mPZR0J3A5gZk8C7YHOJYmuBrklLjN7H5gqaYP00A7AK3nFE0IIFeoZYD1J3SS1xRtfjK62zjv4MRpJ3fHENaOkUWbkXfd2AnBz2llvAofnHE8IIVQUM5sr6Xjg30Ar4Foze1nSUGCcmY0GTgWuknQK3lDjMDOzvGLONXGZ2XigT54xhBBCpTOz+/FGF9nHzs7cfwXYutRx1SbvflwhhBBCg0TiCiGEUFbqTVyS1pf0kKSX0nIvSb9p+tBCCCGE7yqmxHUV8CvgWwAzexFvdRJCCCGUXDGJa2kze7raYzGERQghhFwUk7g+krQO3gQSSfsC05s0qhBCCKEWxTSHPw4YBmwo6V3gLeCgJo0qhBBCqEWdiSsNd3+sme2YBsBdysy+LE1oIYQQwnfVmbjMbJ6kbdL9WaUJKYQQQqhdMVWFz0saDdyBj+AOgJnd1WRRhRBCCLUoJnG1Bz4Gts88ZkAkrhBCCCVXb+Iysxj4NoQQQrNRzMgZXSTdLenDdLtTUpdSBBdCCCFUV0w/ruvwuVlWT7d70mMhhGaqqqoKSQtuVVVVeYcUQqMp5hrXymaWTVTDJZ3cRPGEEBpBVVUVY8aMAVjwN4SWopgS18eSDpbUKt0OxhtrhBBCCCVXTOI6AvgJ8D4+1NO+xEzFIYQQclJMq8K3gT1KEEsIIYRQr2JaFV4vafnM8gqSrm3SqEIIIYRaFNM4o5eZfVZYMLNPJW3adCGFEEpl4obdl2j77pMmNlIkIRSvmMS1lKQVzOxTAEkrFrldCGVj7TPuq/W599/8uN51plywW6PHFEKoWTGNM/4IPCnpt5J+BzwBXNS0YYUQQigVSQMlvSppsqQzalnnJ5JekfSypFtKHWNWMY0zbpA0Dh+r0IC9zeyVJo8shBBCk0vTV10BDACmAc9IGp09zktaD/gVsHW6XPS9fKJ1tZa4JC0tqQ1A+gf+A7QFNixRbCGEEJpeX2Cymb1pZt8AtwKDqq1zFHBF4ZKRmX24JG8o6ceSiqnxq1FdG/4LWDu9ybrAk8D3geMkXbC4bxhCCKHkWksal7kNyTy3BjA1szwtPZa1PrC+pP9K+p+kgUsYz/7A65IuktTgwlBdVYUrmNnr6f5gYISZnSCpLfAsUGM9aAghhGZnrpn1WYLtWwPrAf2BLsCjkjbOtjhvCDM7WNKywIH4MIKGj4E7wsy+rG/7ukpclrm/PV5VSCpKzl+cYEMIITQ77wJdM8td0mNZ04DRZvatmb0FvIYnssVmZl8AI/GqydWAvYDnJJ1Q37Z1Ja4XJV0s6RRgXeABgGxn5BBCCGXvGWA9Sd1SjdoB+IwgWaPw0haSOuNVh28u7htK2kPS3cAYoA3Q18x2AXoDp9a3fV1VhUcBJ+HXuXYys9np8R7AxYsbcAghhObDzOZKOh74N9AKuNbMXpY0FBhnZqPTcztJegWYB5xuZksy2Po+wCVm9mi1WGZLOrK+jWtNXGb2FfCdRhhm9gTelyuEEEILYGb3A/dXe+zszH0Dfp5ujaEKH7QdAEkdgFXMbIqZPVTfxovdHLGxpKlSnpd0b96xhBBCKIk7WLStxLz0WFFyT1x4dWQMeBZChYtZmytK69TQD1jQ6K9t0Rs3SUhFktQF2A04j8YrgoYQylBzn7U5xrNsVDMk7ZGunyFpEPBRsRvXmrgk3cOiTeIXYWaNMUfXpcAvgE51xDEEGALQtm3RCTmEEELzdQxws6S/AMI7QB9a7MZ1lbiatOWgpN2BD83sWUn9a1vPzIYBwwA6duxYayINIYRQHszsDeAHkpZJyzMbsn2t17jMbGxdtyWMG2BrYA9JU/AOaNtLuqkRXrcixPWAEEI5k7QbcCzwc0lnSzq7vm0K6r3GlUYF/j3ef6t94XEz+/5ixLqAmf0KH22YVOI6zcwOXpLXrCTN/XpACCHURtLfgaWBHwFXA/sCTxe7fTGtCq8D/gbMTW9yA9BiSkZRcgkhhJLbyswOBT41s3OBH+KjcRSlmMTVIXUIk5m9bWZVeEvARmNmY8xs98Z8zWJVVVXRr18/+vXrh5lF4mpCcZIQQkjmpL+zJa0OfIuPV1iUYprDf53mTXk9DQvyLrBMg8MMFS+qN0MIyT1p3Ns/AM/hLdivKnbjYhLXSXhd5InAb/GR4gc3OMwQQggVLxWEHkpTotyZRk1qb2afF/sa9SYuM3sm3Z0JHL44gYYQQggAZjZf0hXApmn5a+DrhrxGXR2QLzWzk2vriNxIHZBDCCFUnock7QPclQbwbZC6Slw3pr8xhUkIIYTGdDQ+zN9cSXPw0TPMzJYtZuO6pjV5VlIrYIiZHdQooYYQQqh4ZlbrMH/FqPMal5nNk7SWpLbZkXxDCCGExSVpu5oerz6xZG2KaVX4JvBfSaOBWZk3+FNREYYQQgiLOj1zvz3QF3gWb7Ver2IS1xvpthR1jOIeQgghFMPMfpxdltQVny2kKMU0hz+34WGFEEIIRZsGdC925WIG2V0ZnzOrJ4sOsltUkS6EEELIknQ5C7tZLQVsgo+gUZRixiq8GZgEdAPOBaYAz9S1QQghhPIhaaCkVyVNlnRGHevtI8kk9VnCtxyHX9N6FngS+GVDZgcp5hrXSmZ2jaST0jxcYyVF4gohhBYgdXu6AhiAV9k9I2m0mb1Sbb1O+BCATzXC244E5pjZvEIMkpY2s9nFbFxMievb9He6pN0kbQqsuHixhhBCaGb6ApPN7M3U7elWYFAN6/0WuJCFI7sviYeADpnlDsCDxW5cTOL6naTlgFOB0/BJv05pSIQhhBCarTWAqZnlaemxBSRtBnQ1s/sa6T3bm9nMwkK6v3SxGxfTqvDedPdzfCLJEEII5aW1pHGZ5WFmNqyYDdNo7n8CDmvEeGZJ2szMnkvvsTnwVbEbF9Oq8CLgd+lF/wX0Ak4xsxYzC3IIIbRwc82stgYV7wJdM8td0mMFnYCNgDGSAFYFRkvaw8yyybAhTgbukPQePk7hqsD+xW5cTOOMnczsF5L2wlsU7g08CkTiCiGE8vcMsJ6kbnjCOgD4aeHJNE9W58KypDHAaUuQtDCzZyRtCGyQHnrVzL6ta5usYq5xtUl/dwPuaMhkXyGEEJo3M5sLHA/8G5gI3G5mL0saKqlJpq+SdBzQ0cxeMrOXgGUkHVvs9sWUuEZLmoRXFf5f6pDcGK1KQgghNANmdj9wf7XHzq5l3f6N8JZHmdkVmdf8VNJRwF+L2bjOxJUuyt0D/AH4PI0WP5uam0qGEEKoxcQNix7RqEbdJ01spEiahVaSVJhEMvUla1vsxvVNazJf0hVmtmnmsVlkRokPIYQQGuhfwG2SrkzLRwP/LHbjYq5xPZSG+dDiRBdCCCFU80vgYeCYdJvAoh2S61RM4joauAP4RtIXkr6U9MXiRBpCCCGY2Xx86Kgp+Mgd2+MNQ4pSTAfkmIMrhBDCEpO0PnBgun0E3AZgZg0a3KLeEpfcwZLOSstdJfVteMghhBAq3CS8dLW7mW1jZpcD8xr6IsVUFf4V+CELO6TNxEcSDiGEEBpib2A68IikqyTtgI+c0SDF9OPa0sw2k/Q8LGhvX3SzxRBC01n7jNrHPH3/zY/rXafoZlwhNAIzGwWMktQR71Z1MvA9SX8D7jazB4p5naKmNUlt7Avt7VcG5i9O0CGEEIKZzTKzW8zsx/jYiM/jLQ2LUkzi+jNwN54VzwMeB85fnGBDCCGELDP71MyGmdkOxW5TTKvCmyU9CxTqIvc0sxbVhTuEEEL5KGZakz8Dt2bHlWoMkroCNwCr4NWQw8zsssZ8jxBCCC1PMVWFzwK/kfSGpIsl1TanS0PNBU41sx7AD4DjJPVopNcOIYTQQtWbuMzsejPbFdgCeBW4UNLrS/rGZja9MPulmX2J95peo+6tQgghVLpimsMXrAtsCKxFA4bmKIaktYFN8SFAqj83BBgC0LZttMIPIYRKV8zIGRelEtZQ4CWgT2rC2CgkLQPcCZxsZt8ZAzG1NuljZn1at25Ing0hhNASFZMJ3gB+aGYfNfabS2qDJ62bzeyuxn79EEIILU+tiUvSZunuM8CaktbMPl+4PrW40jQp1wATzexPS/JaIYQQKkddJa4/pr/tgT7AC3g/rl7AOHz8wiWxNXAIMEHS+PTYmWkK6RBCCKFGtSauwjDzku4CNjOzCWl5I6BqSd/YzB5nMQZXDCGEUNmKuca1QSFpAZjZS5K6N2FMIYQWLAYGDkuqmMT1oqSrgZvS8kHAi00XUgjlZ+KGS3Yu131SjKIW8iNpIHAZ0Aq42swuqPb8z4Gf4QNHzACOMLO3Sx5oUszIGYcDLwMnpdvLwGFNGFMIIYQSSbN/XAHsAvQADqxhFKPn8a5QvYCRwEWljXJRxYycMcfMLjGzvcxsL3wIqGgFGEIILUNfYLKZvWlm3wC34nNlLWBmj5jZ7LT4P3wqktwU1aNX0qbAgcBPgLeA6HNVIktyPSCuBYQQirAGMDWzPA3Yso71jyTnw0td/bjWx5PVgcBHwG2ACq0Ny01cEA4hVLDWksZlloeZ2bCGvoikg/HuUf0aLbLFUFeJaxLwGLC7mU0GkHRKSaIKZS1OEkJoduaaWW0ze7wLdM0sd0mPLULSjsCvgX5m9nXjh1i8uq5x7Q1MBx6RdJWkwkSSIYQQWo5ngPUkdZPUFjgAGJ1dIV0uuhLYw8w+zCHGRdSauMxslJkdgI8I/whwMvA9SX+TtFOJ4gshhNCEzGwucDzwb3zmj9vN7GVJQyXtkVb7A7AMcIek8ZJG1/JyJVFv4wwzmwXcAtwiaQVgP+CXwANNHFsIIYQSSEPt3V/tsbMz93cseVB1KKYf1wJm9mmaZmSHpgoohBBCqEuDElcIIYSQt0hcIYQQykokrhBCCGUlElcIIYSyEokrhDp89vjNfD31Jb6e+hJvX7g7nz1+c94hhVDxihqrMIRKtfw2B7H8NgflHUYIISNKXCGEEMpKJK4QQghlJRJXCCGEshKJK4QQQlmJxBVCCEsgWp6WXrQqDCGEJRAtT0svSlwhhBDKSiSuEEIIZSUSVwghhLISiSuEEEJZicQVQgihrETiCiGEUFYicYUQQigruSYuSQMlvSppsqQz8owlhBBCecgtcUlqBVwB7AL0AA6U1COveEIIoVLVV4iQ1E7Sben5pyStnUOYC+RZ4uoLTDazN83sG+BWYFCO8YQQQsUpshBxJPCpma0LXAJcWNooF5Vn4loDmJpZnpYeCyGEUDrFFCIGAden+yOBHSSphDEuQmaWzxtL+wIDzexnafkQYEszO77aekOAIQBt27bd/Ouvv270WPr37w/AmDFjGv21m1I5xl2OMZerctzX5RhzOZD0DTAh89AwMxuWnqv3WCzppbTOtLT8Rlrno1L9D1l5DrL7LtA1s9wlPbaItHOHAXTs2DGfLBtCCOVtrpn1yTuIxpJnVeEzwHqSuklqCxwAjM4xnhBCqETFFCIWrCOpNbAc8HFJoqtBbonLzOYCxwP/BiYCt5vZy3nFE0IIFaqYQsRoYHC6vy/wsOV1nYmc5+Mys/uB+/OMIYQQKpmZzZVUKES0Aq41s5clDQXGmdlo4BrgRkmTgU/w5JabmEgyhBAqXE2FCDM7O3N/DrBfqeOqTQz5FEILVFVVxdixYxk7diySqKqqyjukEBpNlLhCaIGqqqoiWYUWK0pcIYQQykokrhBCCGUlElcIIYSyEokrhBBCWYnEFUIIoaxE4gohhFBWInGFEEIoK5G4QgjNQnSaDsWKDsghhGYhOk2HYkWJq0zF2WkIoVJVfOIq1wRQVVWFmS24lUPc5bqvQwjNi3KcUqXBOnbsaLNmzco7jBBCKCuSZptZx7zjaCwVX+IKIYRQXiJxhRBCKCuRuEIIIZSVSFwhhBDKSiSuEEIIZSUSVwghhLISiSuEEEJZicQVQgihVpJWlPQfSa+nvyvUsM4mkp6U9LKkFyXt35QxReIKIYRQlzOAh8xsPeChtFzdbOBQM+sJDAQulbR8UwUUiSuEEEJdBgHXp/vXA3tWX8HMXjOz19P994APgZWbKqAYHT6EEFq+1pLGZZaHmdmwIrddxcymp/vvA6vUtbKkvkBb4I2Gh1mcSFwhhNDyzTWzPrU9KelBYNUanvp1dsHMTFKtA9xKWg24ERhsZvMXN9j6ROIKIYQKZ2Y71vacpA8krWZm01Ni+rCW9ZYF7gN+bWb/a6JQgbjGFUIIoW6jgcHp/mDgH9VXkNQWuBu4wcxGNnVAkbhCCCHU5QJggKTXgR3TMpL6SLo6rfMTYDvgMEnj022Tpgoo5uMKIYQWLubjagSS/iBpUuqodndTtvcPIYTQsuRVVfgfYCMz6wW8BvwqpzhCCCGUmVwSl5k9YGZz0+L/gC55xBFCCKH8NIfm8EcAt9X2pKQhwJC0aJK+KklU39UamFvvWs1LOcYM5Rl3OcYM5Rl3OcYM+cbdIaf3bRJN1jijrg5tZvaPtM6vgT7A3tbMW4lIGldXB77mqBxjhvKMuxxjhvKMuxxjhvKNuzlqshJXXR3aACQdBuwO7NDck1YIIYTmI5eqQkkDgV8A/cxsdh4xhBBCKE95tSr8C9AJ+E/qqPb3nOJoiGIHpGxOyjFmKM+4yzFmKM+4yzFmKN+4m52y6oAcQgghxJBPIYQQykokrhBCCGUlElcjkbRM3jGEEEIliMTVCFLSGibpoLxjKZDUVVJbSR3TcnzWIYQWIQ5mjcOAfwH7Sdo372Ak7Qb8E2+9eZ2kDcxsfrkkL0kq5rFK1Rz3haStJW2bdxy1KewzST3zjiUsubI4kDV3ZjYLn0RtBHB4nslLUhd8vpzjgbOAp4ExknqWQ/KS1LbQIV3SupI2AJ8yPN/ISi9zsP2hpG0l/Qia376QtD7+nZuRdyy1SVPO7wLcK2mzvOOpLvNZbyfpWEm7xeWH2jXrg1hzlz3zNbMvzew24CbyTV4fAE/go+5/aGYX4weVByStb2bzc4qrXpI2Bn4idxJwJzBC0tWS1so5vJJLB9vdgSuAXsBlkn6Sc1iLSCWYUcBDZjapuZ4YSeoFXArsY2bPSVpD0nKSWuUcGrDgsx6A9/X6GLgZGNxc92feYqcsJklKX7ZdJV0p6feSNjazEfiX7hBJB5QwnsJn2RZYBjiicGZuZpcBlwFnSmrfHKuako2BPfHS4vbA1sBWgIDT8gsrH2meuuOAgcAXwJfA2OZyMJO0rJm9DIwDdpO0THM6Mar2PTfgLqCrpN8A9+AnmZvkENoi0olaO2D/dHsTeAu4qzntz+akWfwAylFKWrsBv8UTVTe8dLClmd2C/0gOk7RqKRJF5gveDjgHOEHSCZlVbge+Br5uhlVNAkj77Q6gP7A8sIKZzcEP3tukmQIqQkpOs4Fp+LToRwGHm9kHwC55X6uRtC5wvqStzOxQ4HXgVknt0/O5nhxlTix/JGlvfD+uARwKvAHsgieHXKoNJS0tae20uC4+/N6z+AnaFcBeZjZd0kGS+ucRY3MWiWsxSeoA9AUOAlYAVscT2HBJfczsevxA835TJgpJWxVKdilRPQscBlwLDJV0SroG0R/YHE8IzUbhAFNYTtWtf8FLF9tJWiMlr1GU51QWRctc59ga2NnMvsGT1wXAEDN7LT13Pl4KzdMMPLafSOprZj8FPgXul9Q+75OjlLT2BC4B5pjZp3jy3zfVinQGtgUm5RTi2sDPJZ0L/Bn4Hv797gn8ysymSOoN/Jo4Tn9HDPnUANUPspKWwxPBCOBIM5so6SlgFaCXmX1Rgph2ww/01+MTcv4e2Amff6crfpb5DvAj4DAzm9DUMS2OdE2rF/AScA3wA3yutnnpsSOAPcxsYm5BNqFMCWEX4HL8pOcxSSsDp+MnHncAhwC/MbPRJYytlZnNS/c3A9qa2f/S9//nwIrAVWb2oqTbgYvMbFyp4quJpBWBW/HvzXS8SrCvmf1NUj/8ZOCCwhRLOcTXCvgjcDhwvpldKKkN8Cd8HNeV8JPhc0v5WZeLSFwNJGl7YC284cN9kjrhieM3+FnUnsD1ZvZiCWMagH/h/2dmR6X68v1SnK3xH8hSpUikxcqeBEjaFI9/OLARsBpePdgzPf4Y8FczeyufaJuOpNXN7L3CfbxkebSZPS+pD7Ah8Ah+rW8eMCMltEVOopowvs7Aw8BuZjZV0un4NbdfmdnT6TrctcCy+Fx7TzV1TMVISfUBfIb1TsB8YDfgSuBCoJuZvVKq/ZiJK/u9/yl+srY+fswozFO4Cf5Zz8sjxnIQiasIkpZKTck3B0YC/8BLN6+a2a8lXQW0BwYAR5nZPTnEOAi4CjjRzG5N10gGAxvgZ5aflTqm2lT78f4A2AL4wMxul9QVOAYvLZ6MJ683zezdvOJtSvKZES4rlCQl/RH/32ekv7OAaWZ2eo4xXgoMAn5oZu9LOgUvAf4+lbz2xUsOp5jZaznFWCix9sUbYnyU/p4KjDSzsem7diRwTKEEmVOMm7LwJGS6pKPwE97zganA7mb211LHV1bMLG613PDWeUun+/3w5rTbp+VeeAOM09LyqkD3wvE4p3h3A14EDkjLSwGd8t6PmfiEl/wKy0cDE1LMo/DGGOAnBZfiiTiXfVni/fJ9/Ppoa2BT4A/49ZfW+GSrV+axH4BWmfsX4gfV1dLyKelz+xUwHtiuGezH3YHngJOAl4Gtqz03AU8KecTWKv3dBb+u9ivg7cwx40i8dDgFGJD3vmzut9wDaK43/NpVFbByWh6MNxg4Ji23wptv/xu4MO94M3Hvkg4w++YdSw2xtc7c3x24v5Bc8VLs7zLJa/XCvm/pt/Rdehq/TpR9vF86EP84x9iWz9yvwlvnFZLXT/Bq6NwPtPg1oYfw68uH4X0Zv5eSfyd8ZJs90rolOwnIfofxk90XgHWAvYD3gZnAJun5DYDN8t6X5XDLPYDmfMMbNnQF9k7Lh+Ade7dIy63Sl7FZfdnwKsvv5x1HtZhWBkbjpa7WeJP9zzP7cs2UvC7JHixb4o2FVfTr4c2xlW734X2LlPbHNcCg7DYliG2VQiICdsVrFW7HWzmCt3KbAqyZlts0k325IvA3vOTyGLB+eny39N1bOrt+iWJrnfbflWm5HdAd2A54Lj12IX79rVkdQ5r7LfcAmuONRatIBuPVWHum5Z8Br5Cphohb0ft16ZRUl0/L56V92zstr423BGvxJS38msbTeMOH64ADU8K6B2+l2gpYNq1bqqTVCm+FNwI4Fi+19MOHDvsTcFxa74/4CC3tsr+VEu+/BQkr89gwvIn+Wml5G7zE2iPHz3kdvFbmj5nHTgSGpvuD8JLijnl/J8vplnsAze2W+UGsDxyY7h+MX4PYKy0fi/duXy7veMvtlg7Y04Dl0oH6THxop83T863zjK9E+2AFvKVg93Tw3wMvLfRNyeNxYOOcYlsTv/Y4Argu8/gg4F5glbS8QTPYj7umfTUMb+CzLvDXlAhOwKvl9mgGcXZLJyiXpOV9gRvTCcGLeNeZkp2gtIRbdGzLyLT66Y/38zhd0p5mdhNeR76XpJ+Yt/jpZ2af5xhuWTKzUXipdRzehPr3wES8M2Y7vNqkpWuH/++Y2dfAWLz14C5mNs/MtrGc+tuZ2Tt4leVTQJ/UiRfzptrz8WG4wEfKyE1qPXgscC5+7flgoAc+8sSdeGfok8xsdKlH8ch0JF9D0lrm3TiOAHpKOs/MRuLHk+XwLgQvQvMbPLk5i+bw1cinZrgK/wEMwnuzjzVvYv4zYEf8B/FBjmGWPUm74lVOPzSzzyStZGYf5x1XU8icEHUGPjWzeZLOwhsP/MnM3pK0P95J/CRgruXQXLtazN/D+wL2wZPYY/gMCAeY2XM5x9YZuAV4y8yOTh13T8GvR/8LeMDMvs05xj3xY8hHeA3D+fg4olcBL5jZaZl1o59WA1V8iUvSOmlolYIfAreb2b34YK8TgCGp5HU18PNIWkvOzO7HmwQ/lH64LTJpwSKjvN+Jj9LfB7/uMR0YmfpFXYgPqvp13kkLwMw+xK8/Po9XaV2C9xEsedJK4/rtkO73wrsKjAJ2kLR7SlIXAx8CP8arYnMjH0fyFLyj9qP4MWWmmU3Bq2H7SupRKJlF0mq4ii9xSSo0S33BzGZL+jGesE4xs1fSOg/jfS+uMLOXCx2S84u65ZCPKD4z7ziaUhom6QLgDLyV20b4mfd4vCPvinhn9rE5hVgrSavgJa+HC7+HHGLoiFcpb4w38DkIH8bsQPx60d/M7P40jNJaZvZmDjFmO9WvD+yNT09yBHCQmb0paRMzGy8fy3FOqWNsSVrnHUCe0pft7jSu2cuSDsOvN/wQv561DN5keym8X9FewMuRtBpPBSSt1fHqv1mptPJcGpdxCD6+5J15nHFnqi/bAvPNbG5NJ2Rm9oGka81sdtpuwbiFpWJmsyT9E09YT5rZ5BTL/fjoGKdJapOuw5U8aaUYLY3MsSp+QrIlXnV5YEpaOwO/lbSXtdBRYEqpYqsKMz/cXfHWXZfgLZLWxVsQzsNHbb4B+D+8f82yaiYTz4Wy8TV+fWglSYNhwfxoT+NDW62UR1Dpu78H3kfrCkk/sBpmyE6JarbSdCWlTFrVYnkWL8V8K+maFMsM4Eng70AuyUCuEOd2+ESVU1Jck4GdJB2OdycYGkmrcVRs4ko/3L54nXNrM/szPir3LXhnxQvwxhkD8SbCQ4EbmsP1h1A+0rW724Grga0kHZIevxg43sw+yiMuSRvgI7uPwK/j3iBp22zyKpSu5APpjlKJZqGWtIKk1VIsO0u6EO+KMhYffPl7kobJB6M9E3jcSjwafWEfpdLyRul61QS8ChMzuwjvk7gSPjL9yWZ2b6lbOLZUFVtVKGlZfFT3z80H4JSZ/V3SPGC0pP3N7NG03t7AfuazvYbQIGb2haT78ObkO6eEMJx0kCu11MDhPGCM+fxnSJoFXC3paDMbI6l1qj5cDm9UUmVmb5cgto54a7wvJL0EXISfUJ4hqauZ/SaVXK/BO26fZWl0/VKRDwS9u6SrWTj7wgx8jM2+kj7GR9iZADyYrQ6PhhiNo6IaZ2SqB9fBxy9bGm/i+8t0ICmsdwwwyczGpOU2eTevDeUvNePeHRhnZi/lGMcyePXaCniNw7vpd3E0PqTTJmb2SSpp/QPva/R4CePbF5+PbUXgv2Z2jaQuwG3AI2b2m7Teauajq5d6apKV8WveH+AjdXzNwgGS78CnU5mNjwRzopnl2uetJaqoxAUL+ldU4Z0WnwS+wUtUF5vZtdXWjf4VoV6ZE6KO+BxKc7KPV1u3nXmn45I1dMjE1wdPVtPxYcuuAubgExm+m9btYmbT0rXcM4HHCidwJYgzO2HlVniT8nnAGeYzAq8G/BMvxZxWx0s1ZYytzWxuun8PPgr9xYUqX0nn4PPi/VvSqmb2fh5xtnQVdY1L0kr4BfGDzGxbfMSGz/BEdqqk/8uuH0krFCMlhUH4iBPDUomh8PiCaxrpwPy1pA7p+ZJcL01x7IYnqgPw6rfz8BFM2gHnphINZjYtbTYfH6JoTCliTMl1nqTekq7Exxi8EO/Au6ekNc1sOj7M052liKkmqfp0m7Q/R+KtBw9OrUfBS4kD0/0ZecRYCSoqceGjYHQCOqflm/CpBNbFf8TNclr70LxJWg8/IfoTXrV2vnx22wXJq1pDh3vSNk0Z04qpapLUIvBYvG/ikcCheJ+o4/A+i8uThqAqMFeSrgqZEuGOeNVlPzy5voQ3luoG/DQlr/fM7MlSxFU9xvR3K/wE4Kf4PhR+MnBQ+mzvwBvjlLQFZqWpqMYZZva5pDuBfpI+NrOXJN2Nz2E1Lq5jhYaS1B1vPPCCmY1Oj30B/CVdG70enzxzXrWGDk123SNdwzocuD3FMEfSfLxaHLxj7BX42IhzUkOkPGYEbmU+NqPJZxe/Hu+r9V+8JHM5nljb4nN/5dYVJdMK+TzgcPOZn9fFO5Rvgg+UvCzwu0JVcGg6lVbiAr/A2wa4RNJQ/Af8z0haYXGY2UT8OscmktaXd+L9D359Zqik1TMlrXuBc8zssSaOaSaeBL4FTkwlrvvxVoNrpSrw1sA6KcmVXKpa21veARp8wNm7U9Xk7fisz93wktcT+Ezjb+URa8ZyeF+t7dPy28Bb+Mg7BwCjImmVRsU1zgCQ1AkfHWNdYLyZPZFzSKFMZKq1NsUPZM+nkvwl+KC55wKTzfsgdTazj+R9fo7Hv2uPNmFs7fC5zKbJB8ndHtgBr3K7Bh9+6DTgWny4pJPN7J9NFU89sfYBvsIbirQDlgEeBI42s3+ldS7Bm5g/gQ8GMD/v687pWuYf8Wb4IyT1wwcv2N7MPssztkpSkYkrhCUhH3HiXHz0i2WBa83sP5Iuwmc1/qWZvVZtmw5m9lUTxiT8ZKwn3kBgO2B//HrRTvisxX9O6ywNfJnHtaJCrCn5r46P4fgift1oAHAUPlfVZOBS/JrhCmb2izxirYl8PNOb8Wbv84GbCtXEoTQqsaowhAaRtGyq6kPSRvjYgzvh3Sn6AvtL2iUdXN/GE8MimjJppdc3PDntAvwC+JeZzTSz+/CJFdcCTgfeMLMH8kpahVgl7YSX+kbjgw7/FE9gf8KT15n4eI7jgQ0ldcy20MyTmd2Dz/+1LvCMpTm/mkt8laCiGmeE0FCpWvl0/FooeBPnn+MlmxPxGZ2PAs5JpaqTc4hRqRXge6lv0Sx8XM3tzOzRdGBth08C2b7U8VUnHy1/EHCrmT0m6VNgMN5C7zq8JNMKLzVego9aMyuveGuS9ukc4FpJb5jZXXnHVEmiqjCEeqTrRe3x60WjzOxTSccBbczs0tT0fTe8RdnEEsdWqHbbHG909BXwKnA2XqNyA94NpCvwSuoLVXKZ7gDCB8z9Bi+1vJHi/xF+HXAs3jBjqfT82OrVrs2JpAH4/5DLqPSVKqoKQ6iFFo76PQPvs7MHsFtqiTcZ+IOk04FzgKtLnbRgQbXbzviArvvhI0vsBVxGGnUCvxZneSStVGIlJa3t8M65l+ODz+5SaGxhZo/gszOMNZ9M8yvgmuactADM7D+RtEovSlwh1CBTkukHrGNm10raB08KD5rZ8DR6wo/wqeIfyCHGpfAO9XcCfzKfTHFzvBPsqcA9eJVmGyvx6OkpvqWBf+FJ9GV8XNDn8Knst8WvEQ01s79U2y6GWgt1imtcIdQgJa1d8ZZtQ9Jjd0qaCRwqqQ3e+u3+tG7JDraZ92qDj7n5FPBVqo57Vj5R5eHAaDN7oYbtSsJ8Hq9L8FLfLOBIM3siddx9B9gKOFPSymZ2Tma7SFqhTlFVGEINJK2AdyI+wHyajx9J+iU+LNh1+PWulTJVXaVMCCYfLPouPLEOBHZmYWvG2XgzbVXfrlQxZt7zbnzE+c1ZtOPuNOANvMHIf0odVyhvUeIKIclUDy6TGmA8A1wu6VW8ldt84M9mtq+kl6zEI39n4lseOAzvS2TANvgYfx0ktU7LZ1kaxTxvZvagpMPwa4JvpI67n+FTvPzRzN6K6sHQEJG4QmCRpLAHsIek3+BNsb8AHjazp9NoD6dLalvqpAULSlpbApsBz5rZLSn2WXgDkaWA4cAtZvbf5pQMzOxuSXOB6yUdgE+nco6ZfZ6ebxZxhvIQjTNCSFLT5guB4woddOVjD86XtAs+wOpQMxtV4rgKSXUrvJpyMj681C/waeu/lXQwXiW3hZVoVPfFIWlvYChwlJk92ZySaygfkbhCxcskp/Pw+Z9uw1sLDgRex5uaH4nPxjs6j4NtKmn9Dvi5mU2Q9Ft8OpKRwBMpea1haULI5kzSimb2Sd5xhPIVjTNCgO7p79N4S7f78MFdX8QH0v0Mv2aUS9JKlsOT6YC0PBT4BB9xYhuAckhaAJG0wpKKEleoWJkOxq8BD5nZ0fKZgNuZ2RuSNsEnG93TzCbnFWeBfGTy84HzzOyW1BBjKHCzmb2cb3QhlE4krlCxUiOLb9LoDk/jozYck57bGR9N/efmA9U2C6lv2W+By81seM7hhJCLSFyhIknaGOiNl7Smp1EeJgD/NrNjU2OMWdaE82ctrtTy8QJgR+ADiyniQ4WJ5vChomSuUfXGh2+aL+lhM3tf0n7AOElfmtkv8420dula25NmNiPvWELIQySuUBEyCWsNYJqZ3STpE+Cg9Pxd+FT3l+LzVzVrkbRCJYuqwlAx0qC4ZwKP483eL8WHIfopPm3J1sCB5nNERf+iEJqpSFyhIkj6AXA1PoHhqcCWwGPAr4AVgPWAb83sidyCDCEUJfpxhRZP0prA+vhsxevjSescoBtwMTDfzMZG0gqhPETiCi1SmmkXSb2B0/CqwBn4aBhHmNm9wKfAsvikhiGEMhGNM0KLlMb2+zGetNoAG+LTfPQF3kvzaXXH54iKzrshlJG4xhVaJEmr4DMD/8zMJkk6FlgZT2IDgfeAm8zs9hzDDCEshihxhZbqG7wqvFANeBXwd7xq8DLg3jTnVrQeDKHMxDWu0CKZ2ad4iWt7SRuZ2bfA7UBHvAn8V3nGF0JYfJG4Qkt2G9AWn3n3PHzswSq8FLYBxASGIZSjuMYVWjRJy+JTlfQG7geWxqsNB5jZB3nGFkJYPJG4QsWQ9CPg98DRZvZC3vGEEBZPJK5QMSStBrQ1s7fzjiWEsPgicYUQQigr0TgjhBBCWYnEFUIIoaxE4gohhFBWInGFEEIoK5G4QgghlJVIXCGEEMrK/wOnggqe0v0PLAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu = [v['mean'] for _, v in adv_scores.items()]\n",
    "stdev = [v['std'] for _, v in adv_scores.items()]\n",
    "xlabels = list(adv_scores.keys())\n",
    "acc = [clf_accuracy[label] for label in xlabels]\n",
    "xticks = np.arange(len(mu))\n",
    "\n",
    "width = .35\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "\n",
    "p1 = ax.bar(xticks, mu, width, yerr=stdev, capsize=2)\n",
    "color = 'tab:red'\n",
    "p2 = ax2.bar(xticks + width, acc, width, color=color)\n",
    "\n",
    "ax.set_title('Adversarial Scores and Accuracy by Corruption Type')\n",
    "ax.set_xticks(xticks + width / 2)\n",
    "ax.set_xticklabels(xlabels, rotation=45)\n",
    "ax.legend((p1[0], p2[0]), ('Score', 'Accuracy'), loc='upper right', ncol=2)\n",
    "ax.set_ylabel('Adversarial Score')\n",
    "\n",
    "color = 'tab:red'\n",
    "ax2.set_ylabel('Accuracy')\n",
    "ax2.set_ylim((-.26,1.2))\n",
    "ax.set_ylim((-2,9))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can therefore **use the scores of the detector itself to quantify the harmfulness of the drift**! We can generalise this to all the corruptions at each severity level in CIFAR-10-C:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accuracy(y_true: np.ndarray, y_pred: np.ndarray) -> float:\n",
    "    return (y_true == y_pred).astype(int).sum() / y_true.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Severity: 1 of 5\n",
      "Loading corrupted dataset...\n",
      "Preprocess data...\n",
      "Make predictions on corrupted dataset...\n",
      "Compute adversarial scores on corrupted dataset...\n",
      "Get labels for malicious corruptions...\n",
      "\n",
      "Severity: 2 of 5\n",
      "Loading corrupted dataset...\n",
      "Preprocess data...\n",
      "Make predictions on corrupted dataset...\n",
      "Compute adversarial scores on corrupted dataset...\n",
      "Get labels for malicious corruptions...\n",
      "\n",
      "Severity: 3 of 5\n",
      "Loading corrupted dataset...\n",
      "Preprocess data...\n",
      "Make predictions on corrupted dataset...\n",
      "Compute adversarial scores on corrupted dataset...\n",
      "Get labels for malicious corruptions...\n",
      "\n",
      "Severity: 4 of 5\n",
      "Loading corrupted dataset...\n",
      "Preprocess data...\n",
      "Make predictions on corrupted dataset...\n",
      "Compute adversarial scores on corrupted dataset...\n",
      "Get labels for malicious corruptions...\n",
      "\n",
      "Severity: 5 of 5\n",
      "Loading corrupted dataset...\n",
      "Preprocess data...\n",
      "Make predictions on corrupted dataset...\n",
      "Compute adversarial scores on corrupted dataset...\n",
      "Get labels for malicious corruptions...\n"
     ]
    }
   ],
   "source": [
    "from alibi_detect.utils.tensorflow import predict_batch\n",
    "\n",
    "severities = [1, 2, 3, 4, 5]\n",
    "\n",
    "score_drift = {\n",
    "    1: {'all': [], 'harm': [], 'noharm': [], 'acc': 0},\n",
    "    2: {'all': [], 'harm': [], 'noharm': [], 'acc': 0},\n",
    "    3: {'all': [], 'harm': [], 'noharm': [], 'acc': 0},\n",
    "    4: {'all': [], 'harm': [], 'noharm': [], 'acc': 0},\n",
    "    5: {'all': [], 'harm': [], 'noharm': [], 'acc': 0},\n",
    "}\n",
    "\n",
    "y_pred = predict_batch(X_test, clf, batch_size=256).argmax(axis=1)\n",
    "score_x = ad.score(X_test, batch_size=256)\n",
    "\n",
    "for s in severities:\n",
    "    print('\\nSeverity: {} of {}'.format(s, len(severities)))\n",
    "    \n",
    "    print('Loading corrupted dataset...')\n",
    "    X_corr, y_corr = fetch_cifar10c(corruption=corruptions, severity=s, return_X_y=True)\n",
    "    X_corr = X_corr.astype('float32')\n",
    "    \n",
    "    print('Preprocess data...')\n",
    "    X_corr = scale_by_instance(X_corr)\n",
    "    \n",
    "    print('Make predictions on corrupted dataset...')\n",
    "    y_pred_corr = predict_batch(X_corr, clf, batch_size=256).argmax(axis=1)\n",
    "    \n",
    "    print('Compute adversarial scores on corrupted dataset...')\n",
    "    score_corr = ad.score(X_corr, batch_size=256)\n",
    "    \n",
    "    print('Get labels for malicious corruptions...')\n",
    "    labels_corr = np.zeros(score_corr.shape[0])\n",
    "    repeat = y_corr.shape[0] // y_test.shape[0]\n",
    "    y_pred_repeat = np.tile(y_pred, (repeat,))\n",
    "    # malicious/harmful corruption: original prediction correct but \n",
    "    # prediction on corrupted data incorrect\n",
    "    idx_orig_right = np.where(y_pred_repeat == y_corr)[0]\n",
    "    idx_corr_wrong = np.where(y_pred_corr != y_corr)[0]\n",
    "    idx_harmful = np.intersect1d(idx_orig_right, idx_corr_wrong)\n",
    "    labels_corr[idx_harmful] = 1\n",
    "    labels = np.concatenate([np.zeros(X_test.shape[0]), labels_corr]).astype(int)\n",
    "    # harmless corruption: original prediction correct and prediction\n",
    "    # on corrupted data correct\n",
    "    idx_corr_right = np.where(y_pred_corr == y_corr)[0]\n",
    "    idx_harmless = np.intersect1d(idx_orig_right, idx_corr_right)\n",
    "    \n",
    "    score_drift[s]['all'] = score_corr\n",
    "    score_drift[s]['harm'] = score_corr[idx_harmful]\n",
    "    score_drift[s]['noharm'] = score_corr[idx_harmless]\n",
    "    score_drift[s]['acc'] = accuracy(y_corr, y_pred_corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compute mean scores and standard deviations per severity level and plot the results. The plot shows the mean adversarial scores (lhs) and ResNet-32 accuracies (rhs) for increasing data corruption severity levels. Level 0 corresponds to the original test set. Harmful scores  are scores from instances which have been flipped from the correct to an incorrect prediction because of the corruption. Not harmful means that the prediction was unchanged after the corruption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_noharm, std_noharm = [], []\n",
    "mu_harm, std_harm = [], []\n",
    "acc = [clf_accuracy['original']]\n",
    "for k, v in score_drift.items():\n",
    "    mu_noharm.append(v['noharm'].mean())\n",
    "    std_noharm.append(v['noharm'].std())\n",
    "    mu_harm.append(v['harm'].mean())\n",
    "    std_harm.append(v['harm'].std())\n",
    "    acc.append(v['acc'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAEWCAYAAAA6maO/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABG20lEQVR4nO3dd3gU5fbA8e/ZTSX03pOA6A428MfFrogNUSkqShFBBa4F7FjuVcxFVFRs2BFpiiCiIipXLxbUq3gBFUWYtdBDJ3SQkuT9/TETXELKhmQz2c35PE+e7PQzsztz5n3nnRkxxqCUUkpVdD6vA1BKKaXCoQlLKaVUVNCEpZRSKipowlJKKRUVNGEppZSKCpqwlFJKRYWYTlgiMkFERngdR3FEZLGIdAhzXCMiR0U2ougjInNEZIDXcXhNRFaIyHlex+GVkuxL0UhE+ojIf7yOwytRmbDcg9NWEUn0OpayYIw51hgzp7TzEZFjReQ/IrJFRLaJyPci0rkMQowJItLfTfhXeR1LNBCR6iLyjIisEpFdIrLU7a7rdWxQ8AlpWe1LBSyrQuxbxpjJxpgLQuKqVCewUZewRCQNOBMwQBcPli8iUibbTUTiymI+IT4AZgMNgfrALcCOslxABGIuT/2ALcA15bnQaNxmIpIAfAYcC3QCqgOnAllA+yOYnz9fd7Rtk4jvW8WJwm1W9owxUfUHDAO+AZ4CPsw3rC3wA7ATeAuYCoxwh9nAJSHjxgGbgJPc7lOAb4FtwE9Ah5Bx5wAPu8v9EzgK6A8sc5e1HOjjjtsS+Bxnx94MTAZqhsxrBXAP8DOwz41jBXCeO7w9MNeNYx3wPJAQMr0Bjipgu9R1h9UsYtt1BRbi7GhLgU5u/8bATJyD+R/AwJBpMoDpwBvudAOAGsBrbnxrgBGA3x3/KOBLYLu7/m8VEc/bwHp33K+AY0OGTQBeAD5yt/H/gJYhw88Hgu60z7vLHFDEslKBXOByIBtoGDLMD/zD3SY7ge+BZu6wY3EOVFuADcA/QuIbETKPDkBmMd/zvSHLWAJ0zxfjQJzfad7wk4ChwDv5xhsNPFvIeq4A7nOn3wqMB5LcYb8Al4aMG+9+R20LmM8Ad32rFrFNLZx9YxuwGOiS7/t7CZgF7AbOK2SbHPJ7Dt2uedvU/W42u9Pn7WeDgAPAfmAX8EHI+uftS4nAM8Ba9+8ZIDHfvO8ENuL8lq8tZD3D2bcuwdm3tuEcR05w+98DTM837rPAaPdzUftSf5xjztM4x5MRbr//usO/cuPa7W6Dq0ryHUfjn+cBlDhg54B6E/B/7g+2gds/AVgJ3O5+SVe4w/N+/MOAySHzuRiw3c9N3B9EZ5xS5/ludz13+BxgFc7BK879ke0AjnGHN8I92OIcsM93d5Z67o/qmZDlrnB/2M2A5AJ2sv/DSZ5xQBrOAey2kOkLS1gC/A58CHTL2y4hw9vjHNzPd9exCRAI+eG/CCQBbXASeUd3WIa7Hbu50yUD7wGvACk4Z5vzgL+7408B/umOmwScUcR3eR1Qjb8OLAtDhk3gr7P5OJzEP9UdVhfnoH6F+13fjpOEikpYDwDz3M+LgDtDhg11+x3jbscTgTpubOtwDmpJbvfJIfEVl7Dyf889cE4OfDgHl91Ao5Bha4C/uTEchZNkG7nj1XTHi8M5wP5fIeu5Aueg1QyojXPAy9sH7ibkBALnBGZRIfOZCkwsYnvG4+yL/8DZ9zq638kxIdtnO3B6yG+hoG1SXMLKxjk5TQTOdrfFMQV9BwXsS8OB73B+o/VwEslD+eY93F2XzsAeoNYR7Ftt3e/kZJyTn35uHInud7gHqOaO68f5TZ3idhe1L/V3Yxzifu/JhCSsQrZf2N9xNP55HkCJgoUzcA6edd3uIHC7+/ksnLMoCRn/25Af/1HuDlXF7Z4MDHM/3wO8nm9ZnwD93M9zgOEhw1JwzqQuz9vxioi5G/BjSPcK4Lp84xzcyQqY/jbgvZDuAhOWO6wpTmljKU5p4iuglTvsFeDpAqZpBuTk7VBuv0eBCe7nDOCrkGENcM6Ok0P69QK+cD9PAsYATUv43dZ0162G2z0BGBsyvDMQdD9fA3wXMkxwzpaLSli/4yZ+nBLITyHDfgW6FjBNr9DvLt+wCRSfsK4rLB53nIV5y3V/b7cWMt6/cUu9OGfyS4qY5wrghnzbban7uTHOPlDd7Z4O3F3IfGYDI4tYzpk4pWNfSL8pQEbI9plUQGz5f/vhJKyUkOHTgAcK+g7y70s4+0HnkGEXAitC5v0nEBcyfCNuIinhvvUSbiLM95s62/38X+Aa9/P5Id9HcftSf2BVvvn2p+iEFfZ3HI1/0XYNqx/wH2PMZrf7TbcfOF/UGuN+S66VeR+MMX/glFYuFZEqONe/3nQHpwI93Iup20RkG05ybBQyr9Uh89qNc4Z8A7BORD4SkQCAiDQQkakiskZEduBUpeW/SL2aQojI0SLyoYisd6d/pIDpC2SMyTTGDDbGtHTXaTdOAgEnMS0tYLLGwBZjzM6QfitxSmAFxZuKc0a6LmRbvYJzdgjOGZ4A89wWW9cVsp5+ERnpXsjfgXOgId+6rg/5vAeoGhJz6PdhKHqbng6k45QawPnejxeRNm53YdumsP7hOiQmEblGRBaGbLfj+Gt9i1rWROBq9/PVwOslWO5KnO2FMWYtTonrchGpCVyEc+JWkCwO/f3n1xhYbYzJzbeswn43RfUrylZ3fwtdRuMwp21MyDGggGmzjDHZId2hv7FDFLNvpQJ35jt+NAtZ1ps4iQigN4ced4ral6CE26uE33HUiZqEJSLJwJXA2e7BfD1OVdCJInIiTjG7iYhIyGTN881mCs4PpyvOWeofbv/VOCWsmiF/KcaYkSHThiZCjDGfGGPOx9mpg8Cr7qBH3HGPN8ZUxznAhMZ02LzyecmdXyt3+n8UMH2xjDGrca4BHReyji0LGHUtUFtEqoX0a45TPVVQvKtxzgrrhmyr6saYY93lrjfGDDTGNAb+DrxYSCum3jjfw3k4Vaxpbv9w1nUdzgHBmcD5zpsVPjr93PkudH83/wvpn7dOBW2b1UCLQua5G6gS0t2wgHEObjcRScX5jQwG6hhjauJU3eWtb2ExAMwAThCR43BKWMUdgEK3RXOc7zhPXvLrAcw1xoR+z6E+BS4UkZRChq8FmuVrgFTU76awfnsoejvWyhdD6PoUtR/lxZhayLRHrJB96+F8x48qxpgp7vC3gQ4i0hTozl8Jq8h9KW9xRxBiuN9x1ImahIVTtZYDtMa5ztIG56Lv1zhVRHNxqg9uEZF4EbmMw1szTQUuAG7krx8NOKWgS0XkQvfMP0lE8n5gh3FLUV3dHWkfzgXPvDPNam73dhFpgnN9pCSq4Vwf2+WW2m4MZyIRqSUi/xKRo0TE5zY9vg6nDh+cC7vXisi57vAmIhJwd75vgUfd9T4BuN7dJocxxqwD/gM86TZ79olISxE5242jR8h224qzw+UWMKtqONsuC+eA9Ug46+n6CDhWRC5zW07dQsEJAxFJwjnRGcRfv5s2ONcFervTjwUeEpFWbivQE0SkDs41i0YicpuIJIpINRE52Z31QqCziNQWkYY4VbdFScHZFpvcuK7lrwMebgx3icj/uTEc5SY5jDF7cap23sS5DreqmGXdLCJNRaQ2zvXEt0KGzcBpzHErf5UQCvI6zgH1HREJuN9zHRH5hzjNuf+Hk2zudve3DsCl/FWKDddCnO/BLyKdcK5T5fcvEUkQkTNxEvbbbv8NFH5CAc4J6v0iUs/dH4ZRyO+6KGHsW68CN4jIye53lyIiF+edBBpjNuFcVhgPLDfG2G7/IvelMBW0DWYQ3nccdaIpYfUDxhtjVrln8euNMetx6pX74BwUL8Op492CU2X3bugM3B/IXOA0QnZi96DdFac0swlnRx1K4dvHB9yBc7a2BWcny0ss/8L5sWzHObC+W9AMinAXTuljJ86O8FbRox+0H6eU8ilOwvsFJyH0BzDGzAOuxWlxtB2nVV3e2Wcvd9q1OBeBHzTGfFrEsq7BudCe1xJtOn9VH/0N+J+I7MJpeXirMWZZAfOYhFNFs8adz3cFjFMgt0q4BzASJ+G1wqkGKUg3nGsVk/L9bsbhXMjuhHNRfxrOwWMHTnJPdqtJz8c5EK/HuQ52jjvf13Fak65wpyvyezLGLAGexPn9bQCOD43ZGPM2TkvUN3G++xk4jSbyTHSnKa46EHce/8FpxboUp3VZ3nL+BN7BqSIt9LdpjNmHU/oN4lzP2oHTIKAu8D9jzH6c7XIRTiu0F3Gu0wTDiC/Ure58tuHsxzPyDV+P8xtbi1OyvCFkGa8Brd3qtPzTgbPeC3BaJS7CaUF8JA8SKG7fWoDTwvN5N9Y/8oaFeBNne76Zr39R+1I4MoCJ7ja40o0nrO84Gsmhl3yUUhWRiDTHSR4NjTGluv9HRIYBRxtjri52ZA+5pbY3jDEF1nSowkXLd1xSeiOaUhWce53oDpxm/aVNVrVxqnz7lkVsquKJ5e/Y0ypBEakpItNFJCgitoic6mU8SlU07nXSHThVkw+Wcl4Dcaq7/22M+aoMwlMVTKx/x55WCYrIROBrY8xYcR4FU8UYs82zgJRSSlVYniUsEamB00KohdELaUoppYrh5TWsdJwWeePFuY/qe5wWZaE3CSIig3CaJAP8X5UqVVAV0969ewFISkryOJLSi6V1gdhan1haFyif9dmzZ48xxkRTq/ACeVnCaofTlPl0Y8z/RORZYIcx5oHCpklJSTG7d+8ubLDyWIcOHQCYM2eOp3GUhVhaF4it9YmldYHyWR8R2WOMKewm8KjhZcbNxHn2Wt5TB6bj3L+klFJKHcazhOXevLlaRI5xe52Lc/OcUkopdRiv78MaAkx2Wwguw3kSg1JKKXWYqHrShV7Dqtii4drCgQMHyMzMPHihuzDr1zsPim/YsMBHFEadWFqfWFoXKNv1SUpKomnTpsTHxx/SP1auYXldwlKqXGVmZlKtWjXS0tIQKfzB8D6fU1t+zDHHFDpONIml9YmldYGyWx9jDFlZWWRmZpKenl4WoVU4Ud/MUamS2Lt3L3Xq1CkyWSkVjUSEOnXqFFt7EM00YalKR5OVilWx/tvWhHUEjDFsfPIp9i1f7nUoSilVaeg1rCOQvW4d26ZPZ9vbb9Ps5ZdIbtPG65DUkcqoUWDvI76akLG92FFEhDvuuIMnn3wSgFGjRrFr1y4yMjIKnWbGjBkcffTRtG7d+vBFZmRQtWpV7rrrroP90tLSWLBgAXXr1i35OhSjV69eLF68mGuvvZbbb7+9wHHmzJnDqFGj+PDDD8t8+Vd8d4XzIew3qBVtUb9FxY5TtWpVdu3adbB7woQJLFiwgOeff75sggjx9ddfc8MNNxAfH8/cuXNJTk4ucLwOHTowatQo2rVrV+YxVFRawjoC8Y0bkzZ1Cr7q1VnZrz87P/vM65BUFElMTOTdd99l8+bNYU8zY8YMliyJ3G2KOTk5YY23fv165s+fz88//1xoslKHy87ODnvcyZMnc99997Fw4cJCk1VlpQnrCCWkppI2dQqJxxxD5pBb2DJ5stchqSgRFxfHoEGDePrppw8btmLFCjp27MgJJ5zAueeey6pVq/j222+ZOXMmQ4cOpU2bNixdurREy+vWrRuXXXYZl1xyCWPGjDnYv2rVqtx5552ceOKJzJ07l6pVqzJ06FCOPfZYzjvvPObNm0eHDh1o0aIFM2fOBOCCCy5gzZo1tGnThq+//poOHTqwYMECADZv3kxaWtqRb5go9cEHH3DyySfTtm1bzjvvPDZs2AA4Jd++ffty+umn07dvXzIyMujXrx9nnnkmqampvPvuu9x9991ceumlDBgwgAMHDjB27FimTZvGAw88QJ8+fZgzZw6XXHLJwWUNHjyYCRMmeLSm3tOEVQpxtWuTOmE8Vc8+mw0PjWDjk09icnO9DktFgZtvvpnJkyezffuhVYhDhgyhX79+/Pzzz/Tp04dbbrmF0047jS5duvDEE0+wcOFCWrZsedj8nn76adq0aXPwb+3atQeHjRs3jnfffZfp06czevRosrKyANi9ezcnn3wyP/30E2eccQa7d++mY8eOLF68mGrVqnH//fcze/Zs3nvvPYYNGwbAzJkzadmyJQsXLuTMM8+M4BaqWP78889Dtm/e9gA444wz+O677/jxxx/p2bMnjz/++MFhS5Ys4dNPP2XKlCkALF26lM8//5yZM2dy9dVXc8455/DBBx+QlJTERx99xIABAw5+15P1JPgweg2rlHxVqtD0udGsHzGCrFfHcmDdeho98jC+hASvQ1MVWPXq1bnmmmsYPXr0IdU+c+fO5d133wWgb9++3H333WHN7/bbbz/sGlae0aNHM3XqVADWrVvH77//Tp06dfD7/Vx++eUHx0tISKBTp04AHH/88SQmJhIfH8/xxx/PihUrjnRVY0JycjILFy482J13DQuce/uuuuoq1q1bx/79+w+5B6pLly6HfL8XXXTRwW2ak5NDp06d+O233zj66KMr/TYOh5awyoDExdHwwQepd8cd7PjwQ1YPHETOjlK9yVxVArfddhuvvfYakXx6y5w5c/j000+ZOnUq77//Pm3btj3kdRZ+v//guPHx8QebRft8PhITEw9+LuwaTFxcHLlurUIs3/9TlCFDhjB48GAWLVrEK6+8csh2SEk59OESods0dHuLSIHbOHT7QuXdxnk0YZUREaHuoIE0fvwx9vzwAyv7XM2Bdeu8DktVYLVr1+bKK6/ktddeO9jvtNNOO1gamjx58sFqt2rVqrFz584SL2P79u3UqlWL5ORkli1bxnfflVHTOldaWhrff/89ANOnTy/TeUeL7du306RJEwAmTpxYpvNOTU1lyZIl7Nu3j23btvFZJW/gpVWCZaxGly7E1atH5pBbWNGzF83GjCHpmKO9DksVppBm6L/++isQ+cf/3HnnnYc0jX7uuee49tpreeKJJ6hXrx7jx48HoGfPngwcOJDRo0czffr0Aq9jFaRTp068/PLLdO7cmfT0dE455ZQyjf+uu+7iyiuvZMyYMVx88cVlOu/CTD/FSYwV5dFMGRkZ9OjRg1q1atGxY0eWl+H9mc2aNePKK6/kuOOOIz09nbZt25bZvKORPvw2Qvb++iurBw4id88emj7/HCllfKCoiKLh4be2bWNZVrHjlVfCKi+xtD6xtC5Q9utT0G88Vh5+q1WCEZJ0zDGkvTWV+EYNWTVwENs/KPsbKJVSqjLRhBVB8Y0akTp5MlXatmXt0KFsfvVVoqlEq5RSFYkmrAjzV69Os7GvUr1zZzY9+RQbHnoIE+ZTBVRk6EmDilWx/tvWRhflwJeQQONRTxDfuBFZY1/jwIaNNBn1BD597Eq5S0pKIisrS18xomJO3vuwkpKSSjytHbA6Ac8CfmCsFbRH5hueCowD6gFbgKutoJ3pDusH3O+OOsIK2mXbVDKEJqxyIj4f9e+6i7iGjdjw8MOs6n8tTV9+ibhatbwOrVJp2rQpmZmZbNq0qcjx8t4CmxsjTy6JpfWJpXWBsl2fvDcOl4QdsPzAC8D5QCYw3w5YM62gHfrwylHAJCtoT7QDVkfgUaCvHbBqAw8C7QADfO9Ou7XUK1MATVjlrPbVfYhrUJ+1dw1lZc9eNHt1DAnNm3sdVqURHx8f1ttYb7zxRqBit3gsiVhan1haF6gQ69Me+MMK2ssA7IA1FegKhCas1sAd7ucvgBnu5wuB2VbQ3uJOOxvoBEyJRKB6DcsD1c8/n+bjx5OzbRsrevbiz0XFv95AKaVKIU5EFoT8DQoZ1gRYHdKd6fYL9RNwmfu5O1DNDlh1wpy2zGjC8kiVk9qSOmUKvipVWHlNP3Z+8YXXISmlYle2MaZdyN+Y4ic5xF3A2XbA+hE4G1gDlHvrMU1YHkpske68oqRFCzJvHszWqW95HZJSqvJZAzQL6W7q9jvICtprraB9mRW02wL/dPttC2fasuTpNSwRWQHsxMnU2caYyvPqTFdc3bqkTppI5h13sD4jgwPr11Hv1lu1BZtSqrzMB1rZASsdJ9n0BHqHjmAHrLrAFito5wL34bQYBPgEeMQOWHmtxy5wh0dERShhnWOMaVMZk1UeX0oKzV54gZo9riDr5VdYd+99mP37vQ5LKVUJWEE7GxiMk3xsYJoVtBfbAWu4HbC6uKN1AH61A9ZvQAPgYXfaLcBDOElvPjA8rwFGJGgrwQpC4uJoOHw4cY0asXn0c2Rv2kST0c/ir1rV69CUUjHOCtqzgFn5+g0L+TwdKPBx/FbQHsdfJa6I8rqEZYD/iMj3+VqtHCQig/JathT2Tp5YISLUu+kmGj3yCLvnzWPl1X05sGGj12EppVSF4HXCOsMYcxJwEXCziJyVfwRjzJi8li1xcZWjQFjzsu40e+klDqxaxYqePdn3++9eh6SUUp7zNGEZY9a4/zcC7+HcwKaAqmeeQeobr2OyD7Ciz9XsnjfP65CUUspTniUsEUkRkWp5n3Fal/ziVTwVUVLr1qRNmUpc3bqsvn4AO2bNKn4ipZSKUV6WsBoA/xWRn4B5wEfGmI89jKdCSmjahLQ3J5N0wgmsueNOssZPiPknMiulVEE8uyhkjFkGnOjV8qOJv2ZNmo97jbV338PGxx7jwLq1NLjnHsTv9zo0pZQqN143ulBh8iUm0uTpp6jdrx9bJ73OmtvvIHfvXq/DUkqpcqMJK4qIz0eD++6l/r33sHP2bFZddz3ZWyPyFH+llKpwNGFFoTr9+9Pk6afY+8svrOzdh/2ZmV6HpJRSEacJK0pV79SJ5uNeI3vLFucVJb8s9jokpZSKKE1YUaxKu3akvTkZX0ICK6+5hl1ffeV1SEopFTGasKJcYsuWpE6dQkJaKqtvvIlt0wt83JdSSkU9TVgxIL5+fVInvU7Kqaey7v4H2PTc83qvllIq5mjCihH+qik0e+lFanTvzuYXXmDd/fdjDhzwOiyllCozleNpspWExMfT6JGHiW/UiM0vvkj2xk00feZpfCkpXoemVNiOn3h8kcOXrV8W1niL+i0qs5hKI9bWx0tawooxIkK9W4bQ8KHh7P72W1b2vYbsTZu8DksppUpNS1gxqlaPHsTXr0/mbbezomcvmr06hsQWLbwOS0WInsWrykBLWDGs6tlnkzppErl797KiV2/2fP+91yEppdQR0xJWjEs+/jjSpk5h9YCBrLr2Oho/8QTVL7zA67CUUhWIHbA6Ac8CfmCsFbRH5hveHJgI1HTHudcK2rPsgJUG2MCv7qjfWUH7hkjFqSWsSiChWTNSp04hqXVr1tx2G1smTfI6JKVUBWEHLD/wAs6b31sDveyA1TrfaPcD06yg3RboCbwYMmypFbTbuH8RS1agJaxKI65WLZpPGM/aoUPZ8MijHFi3nvpD70J8lfOcRa/5KHVQe+APK2gvA7AD1lSgK7AkZBwDVHc/1wDWlmuErsp5tKqkfElJNHnmGWr16cOW8eNZc+ed5O7b53VYSilvNQFWh3Rnuv1CZQBX2wErE5gFDAkZlm4HrB/tgPWlHbDOjGSgmrAqGfH7aXD/P6k/dCg7//0xq68fQM727V6HpZSKrDgRWRDyN6iE0/cCJlhBuynQGXjdDlg+YB3Q3K0qvAN40w5Y1YuYT6lowqqERIQ6119H41Gj2PPTT6zo3YcDa9Z4HZZSKnKyjTHtQv7GhAxbAzQL6W7q9gt1PTANwArac4EkoK4VtPdZQTvL7f89sBQ4OlIroQmrEqtxycU0f/VVsjduZPkVPdjxyX+8DkkpVf7mA63sgJVuB6wEnEYVM/ONswo4F8AOWBZOwtpkB6x6bqMN7IDVAmgFLItUoJqwKrmUU04m7a2pxDduzJpbb2XNXUO1ilCpSsQK2tnAYOATnCbq06ygvdgOWMPtgNXFHe1OYKAdsH4CpgD9raBtgLOAn+2AtRCYDtxgBe0tkYpVWwkqElu0IG3qFDaPGcPml15mz7x5NBrxEFXPOuuQ8bRlnVKxyQras3AaU4T2GxbyeQlwegHTvQO8E/EAXVrCUoDz4Nx6N99M2tSp+GtUZ/Wgv7PugWHk7NrtdWhKKQVUgIQlIn4R+VFEPvQ6FgXJxx1L2vTp1BlwPdumT2d5167snjfP67CUUsr7hAXcilNvqioIX2Ii9e+6i9TJb4Dfz6p+/dnw6EjiD+hLIZVS3vE0YYlIU+BiYKyXcaiCVTnpJFrMeI9avXqyZeJEHh+XQ8u1mrSUysjIQEQO/mVkZHgdUqXgdaOLZ4C7gWqFjeDe4DYIICEhoXyiUgf5qlSh4bBhVD33XDbfNoARk3KYcaow/QwfOX7xOjylPJGRkcGcOXMADv5XkedZCUtELgE2GmOKfOeFMWZM3s1ucXFe59fKq+rpp3PXAD9fHydc/q3hkYk5NNuopS2lVPnxskrwdKCLiKwApgIdReQND+NRxdiTJLx4iZ/HL/dRayc8Nj6HrnNzkVxNXEqpyPMsYRlj7jPGNDXGpOHcWf25MeZqr+JR4VtwtI87B/qZf7TQZ04uD72eQ6MsTVpKqciqCK0Eo45ecIWdVYSnu/l4touPxlvg8XE5XEVV9KqWUipSKsRFIWPMHGCOx2GETS+4ukT45lhhSXPh77NyGbqsNmdThQnbDZtraOpSSpUtLWGpUttaTRh5pY+HyOJYEhg1NodzfsoFo9WEyqG1EqosVIgSlooBIrzPbuazl0cbNuHGWbm0/0145SIf26pqaauy01oJVRa0hKXK1FpyGN7bz/jzfBy/wvDkqzmcaud6HZZSKgZowlJlzojw77/5uPs6P+trwe0zcrl1Rg5V92gVYUloNZpSh9KEpSJmbR3hgWv8TDnbx8m/Gp4cm8NJv0eutBVrB/iMjAzOPvtszj77bIwxUb8+SpWWJqxKLtIH+Vyf8N5pPu7r72dHFbh3ei5r//lPcnbtKtPlgB7glYp12uiikiuvi+ErGwj39ffT47+5dH9vBrvnzqXxI4+QcsopEVumUiq2aAlLlZvsOGFKBz9pb07Gl5DIqv7Xsn7Ew+T++afXoSmlooAmLFXuktu0If29d6nVty9b33iD5d26s+fHH70OSylVwWmVoPKELzmZhv/8B9XOPZd1//gHK/tcTZ3rr6fukMH49DUySpUrO2B1Ap4F/MBYK2iPzDe8OTARqOmOc68VtGe5w+4DrgdygFusoP1JpOLUEpbyVMopJ5M+831qXNadrFdfZcUVPdhr6wuolSovdsDyAy8AFwGtgV52wGqdb7T7gWlW0G6L87DyF91pW7vdxwKdgBfd+UWEJizlOX/VqjQeMYKmL79E9tYtLL/yKja//DImO9vr0JSqDNoDf1hBe5kVtPfjvO6pa75xDFDd/VwDWOt+7gpMtYL2PitoLwf+cOcXEZqwVIVRrUMHWsycSfXzz2fTM8+yoldv9i1b5nVYSsWCOBFZEPI3KGRYE2B1SHem2y9UBnC1HbAygVnAkBJMW2Y0YakKJa5WLZo89SRNnn6KA6tXs7z7ZWyZOBGTq493UqoUsvPe3O7+jSnh9L2ACVbQbgp0Bl63A1a55w9NWKpCqn7RRbT4YCYpp57KhkdHsqpff/ZnZnodllKxaA3QLKS7qdsv1PXANAAraM8FkoC6YU5bZjRhqQorrl49mr70Io0eeYS9ts3yLl3ZOm0aRl9bolRZmg+0sgNWuh2wEnAaUczMN84q4FwAO2BZOAlrkzteTztgJdoBKx1oBcyLVKCasFSFJiLUvKw7LWa+T9IJJ7B+2IOs/vvfObBho9ehKRUTrKCdDQwGPgFsnNaAi+2ANdwOWF3c0e4EBtoB6ydgCtDfCtrGCtqLcUpeS4CPgZutoJ0TqVj1PiwVFeIbN6b5uNfY+uYUNo4axbIuXWh4//1Uv+RiRPR9W0qVhntP1ax8/YaFfF4CnF7ItA8DD0c0QJeWsFTUEJ+P2lf3ocWM90hMT2ft0KGsufU2srds8To0pVQ50ISlok5CWhqpk9+g3p13sOuLL1h2aRd2fvaZ12EppSJME5aKSuL3U3fgQNKmTyeufn0ybx7M2nvuJVmbvysVszRhqaiWdMzRpL81lbo33cj2Dz/k/s1ZBPbt8zospVQB7IB1aWnu3/IsYYlIkojME5GfRGSxiPzLq1hUdJOEBOrdcgtpU6ewT4Rbt24j8/bb2Z8ZsdtBlFJH5irgdztgPW4HrEBJJ/ayhLUP6GiMORFoA3QSEX2bnzpiyccfz6N16/BhSgq7vpjDss6d2fjkUxF5u7FSquSsoH010BZYCkywA9ZcO2ANsgNWtXCm9yxhGUfekSTe/dM7QlWpHBDho2pVafnxv6l+USeyXn2VpRd2cm44zonY7SFKqTBZQXsHMB3nIbuNgO7AD3bAGlLkhHh8DUtE/CKyENgIzDbG/M/LeFTsiG/YkMaPPUba29NISE1l/bAHWX7Z5eyeO9fr0JSqtOyA1cUOWO8Bc3AKKe2toH0RcCLOzclF8vTGYWNMDtBGRGoC74nIccaYX0LHcZ8qPAggQV/sp0oo+fjjSZ38Bjs/+YSNT4xi1bXXUfWcc2jUyrCujt5wrFQ5uxx42graX4X2tIL2HjtgXV/cxBWilaAxZhvwBc4LwPIPG5P3hOG4OH0whyo5EaF6p060mPUR9e68gz3z5vHk2Bz6zc4h5U+thVaqHGUQ8qxBO2Al2wErDcAK2sXeTOllK8F6bskKEUkGzgeCXsWjYp8vMZG6AwfS8pOP+eIE4aLvDc+9nMNF83Px52jiUqocvA2E3iyZ4/YLi5clrEbAFyLyM87TgmcbYz70MB5VScTVrcurF/kZep2fZQ2Faz/NZdTYHE76PRf0SfBKRVKc+1ZjANzPYV/rCTthiUiyiBxTwuAKZYz52RjT1hhzgjHmOGPM8LKat1LhWF1fGNHTx8grnN3g3um53D81l2YbNWkpFSGbQp4Ajx2wugKbw504rItCInIpMAonE6aLSBtguDGmS5ETKlXRifBDK+GnFsIFPxh6/DeXJ8blMINavMx2r6NTKtbcAEy2A9bzgACrgWvCnTjcVgwZQHucpogYYxaKSHqJwlSqAsvxC//+m/DVccIV/82ly4KqXEAK78/NZdbfhANx2qJQqdKygvZS4BQ7YFV1u0t0V3+4CeuAMWZ7vvcOab2Jijm7k4WJ5/sZt2ANt1KTPnOqcP6PMPkcH3MDAvruLaVKxQ5YFwPHAkl2wALACtphXRIK9xrWYhHpDfhFpJWIPAd8eyTBKhUNVpLNHWxmeC8fexLh9hm5DH8jh5Zr9TxNqSNlB6yXcZ4nOASnSrAHkBru9OEmrCE4GXEf8CawHbitJIEqFY1+SfNxz7V+XrnIR8Mt8OjEHAbPzKH2Dk1cSh2B06ygfQ2w1Qra/wJOBY4Od+JiqwRFxA98ZIw5B/jnEYepVJQyPuGzNsK3ltBtbi4XzzOc/GsOM08WcnvswVelitchKlUqdsDqBDwL+IGxVtAemW/408A5bmcVoL4VtGu6w3KARe6wVVbQLqox3l73/x47YDUGsnBucQpLsQnLGJMjIrkiUsMYo82mVKX1Z6IwpYOfT9sY+szJpcc3hqWdLqLe7bdTo2sXxFchHhyjVInYAcsPvIDz8IZMYL4dsGZaQXtJ3jhW0L49ZPwhOE9cz/OnFbTbhLm4D+yAVRN4AvgBpy3Eq+HGGu4etgtYJCKvicjovL9wF6JULNlUU3imm5/7+/qJa9iQdffdx4oeV7Jn/nyvQ1PqSLQH/rCC9jL3Rt6pQNcixu8FTCnpQtwXN35mBe1tVtB+B+faVcAK2sPCnUe4Cetd4AHgK+D7kD+lKq3fmgppU6fQ+InHyc7KYmXfa8gccgv7V63yOjSl8osTkQUhf4NChjXBuR8qT6bb7zB2wEoF0oHPQ3on2QFrgR2wvrMDVrfCArCCdi5OSS6ve58VtEtUaxdWs3ZjzEQRSeCvi2O/GmMOlGRBSsUi8fmocemlVDvvPLLGjyfr1bHsmjOHWn37UvfGG/BXC+u9dEpFWrYxpl0ZzKcnMN0K2qEvl0u1gvYaO2C1AD63A9Yi936rgnxmB6zLgXetoF3ilkthlbBEpAPwO052fBH4TUTOKunClIpVvuRk6t10Ey0//pjql1zClvHjWXrBhWydMgWTne11eEoVZQ3QLKS7qduvID3JVx1oBe017v9lOA+XaHv4ZAf9Hedht/vsgLXDDlg77YC1I9xAw71x+EngAmPMrwAicrQb9P+FuyClKoP4BvVp/Ogj1Lq6DxsfHcn6fw1n65tvUv/ue6h65hleh6dUQeYDreyAlY6TqHoCvfOPZAesAFALmBvSrxawxwra++yAVRc4HXi8sAVZQbtUVQ7hJqz4vGQFYIz5TUTiS7NgpWJZ8rHH0vz1SeycPZuNT4xi9cCBpJx1Jg3uuYfEli29Dk+pg6ygnW0HrMHAJzjN2sdZQXuxHbCGAwusoD3THbUnMDVfVZ4FvGIHrFycGruRoa0L87MDVoE1c/lf6FiYcBPWAhEZC7zhdvcBFoQ5rVKVkohQ/YILqNqhA1tff4PNL73Esi5dqXXVVdQdMpi4WrW8DlEpAKygPQuYla/fsHzdGQVM9y1wfAkWNTTkcxJOC8XvgY7hTBxuwroRuBm4xe3+GudallKqGL6EBOpcfx01undj03PPsXXqVLZ/+CF1b7yR2n16Iwlhvw5IqahmBe1LQ7vtgNUMeCbc6cNt1h4HPGuMucwYcxkwGqfoqJQKU1zt2jR68EFavD+D5BNOYONjj7H00kvZ+emnGH1xpKqcMnGqFcMSbgnrM+A8nBuIAZKB/wCnlSg0pRSJrVrRfOyr7PrqKzY89jiZg4dQpX17Gtx3L0lW2PuuUlHHDljP8debPnxAG5wnXoQl3ISVZIw5+N4SY8wuEdEHqClVClXPOouUU09l67RpbH7ueZZfdjk1LutOvVtvJb5+fa/DUyoSQts+ZANTrKD9TbgTh1sluFtETsrrEJF2wJ/hLkQpVTCJj6d2nz60/ORjavfrx/aZH7C000VsfvllcvfuLX4GSkWX6cAbVtCeaAXtycB3dsAKu/ATbsK6DXhbRL4Wka9xnjU1uMShKqUK5K9Rgwb33kPLD2aSctqpbHrmWZZ27ky7P/8Evb6lYsdnOJeU8iQDn4Y7cZEJS0T+JiINjTHzgQDwFnAA+BhYXvJYlVJFSUhLo9nzz9N8wgT8NWpy/fYdDNm6jezNm70OTamykGQF7YOXl9zPZVbCegXY734+FfgHzuOZtgJjShanUipcKaecTPr0t5lSvRpH7d/Psu7d2f3dd16HpVRp7bYD1sHLS3bA+j9KcHmpuEYXfmPMFvfzVcAYY8w7wDsisrCkkYYSkWbAJKABTquRMcaYZ0szT6Viifj9fFWlCn/Ex/NQteqsuvY66t54I3Vvvgnx610lKirdBrxtB6y1gAANcXJLWIorYflFJC+pncuhj5QPt4VhYbKBO40xrYFTgJtFpHUp56lUzFkbH0/69Lep0bUrm198kVX9r+XAhg1eh6VUiVlBO+/y0o3ADYBlBe2wX1VVXMKaAnwpIu/jFNu+BhCRo4BSvX3YGLPOGPOD+3knYFPIO1iUqux8VarQeOSjNHr0Uf785ReWd+vOrq/CevyaUhWGHbBuBlKsoP2LFbR/AaraAeumcKcvMmEZYx4G7gQmAGeYv27H9wFDjizkw4lIGs4j6f9XVvNUKhbV7N6N9HemE1evHqsH/Z2No0ZhDuir6VTUGGgF7W15HVbQ3goMDHfiYpu1G2O+M8a8Z4zZHdLvt7zSUWmJSFXgHeA2Y8xh70URkUF5b8nM1vcKKUViixakTXuLmlddRdbY11jZ9xrqbtem7yoq+O2AJXkddsDyA2E/TDPc+7Aiwn1FyTvAZGPMuwWNY4wZY4xpZ4xpFxdX2stmSsUGX1ISjf6VQZOnnmTf77/z+Lgc2v2W63VYShXnY+AtO2Cdawesc3EuO/073Ik9S1giIsBrgG2MecqrOJSKZtU7dyb93XfYUBPufieXfrNziMvW0paqsO7Babx3g/u3iENvJC6SlyWs04G+QEcRWej+dfYwHqWiUkJqKg/09fNRO+HiBYYRr+fQYIsmLVXxWEE7F6etwgqcd2F1xGlwFxbP6tiMMf/FaYevlCql7Dhh4vl+lqTmcuNHuTw2PodXOvuYa3la668UAHbAOhro5f5txnlqElbQPqck89Ffs1IxZP7RPu6+zs/qenD7jFwGfpxD/AEtbSnPBXFKU5dYQfsMK2g/B+SUdCbaikGpGLO5hpDRx0/Pr3Lp+p3h6Mwc7iCOlWgrW1UwO2B1Ap7FeTHvWCtoj8w3/GkgrzRUBahvBe2a7rB+wP3usBFW0J5YwCIuA3oCX9gB62OcB6iXuIZNS1hKxaAcvzD5HD+PXOmj1i54g4ZcTIrXYakKyG1a/gJwEdAa6GUHrEOeOmQF7dutoN3GCtptgOeAd91pawMPAifjXJN60A5YtfIvwwraM6yg3RPnKRdf4Dyiqb4dsF6yA9YF4caqCUupGLawpY+h1/tZwn7+RR1u+jCHxP1aRagO0R74wwray6ygvR+n9NO1iPF74TRHB7gQmG0F7S3uTcCzgU6FTWgF7d1W0H7TCtqXAk2BH3FaDoZFqwSVinFbqwk3sZHrqcGARTU4am0OT3fzs7q+tnmqROJEJPRtv2OMMXlv3GgCrA4ZlolTYjqMHbBSgXT+eq5sQdOG9Yg9N8GNoQRv/tCEpVQlkAOMYTvretXilpm5PDoxh/Hn+/jsRAHRxFWgjBpFD1+xO7zx0puXTTylk22MaVcG8+kJTLeCdokbTJQFrRJUqhJZnObj7uv92E2Fv/87l1vfzyV5n1YRVnJrgGYh3U3dfgXpyV/VgSWdttS0hKVUJbM9RXikp4+ucw1XfZVLy3U5PN3dz/KGWtKqpOYDreyAlY6TbHoCvfOPZAesAFALmBvS+xPgkZCGFhcA90UqUC1hKVUJGRFmnOYjo4+fuBwYMSmHTgty+euFDKqysIJ2NjAYJ/nYwDQraC+2A9ZwO2B1CRm1JzDVCtomZNotwEM4SW8+MNztFxFawlKqEvu1mXD3dX5u/iiX62bnkjlkCI0ffhh/jWKuy6iYYgXtWcCsfP2G5evOKGTaccC4iAUXotIkrLR7PyrT+a1fllXm810x8uIym5dS4dpVRXjsCh8Xzzf0m/Mly7tfRpOnniS5TZuSzSi2GimoCkirBJVSIMJH7X2kvTkZRFhxdV+yXnsNk6uvLFEVR6UpYSlVIZVDqWTDexvY8+seAH7p/wv1utajQfcGBY6bfMIJpL/3Luvuf4CNT4xi97x5NB45krhahz28wBMlWRcVezRhqehSwQ7w0aBB9wYlit9fvTpNnn2GrVOmsPHRkSzv2o0mT46iyt/+FsEow1PSdanoYu23FmmasCoDvbZQIrF2UDwSIkLt3r2p0rYta267nZX9+lNvyGDqDBqE+P1ehxcz9LdWMnoNSylVqCTLIu2dd6jeuTObnh3NqgEDyN60yeuwVCWlCUspVSR/1RQaP/E4jR4ewZ8/LmRZt+7s/vZbr8NSlZAmLFUm8uri9/y6h1/6/8KG9zZ4HZIqQyJCzcsvJ/3tafhr1WTV9QPY+MwzmGx9x5YqP3oNS5UJrYuvHBJbtSL97bdZP2IEWS+/wp75C2jy5CjiGzb0OjRVCWgJSylVIr7kZBo//DCNH3+MvbbN8m7d2TlnjtdhqUpAE5ZS6ojU6NKF9HemE9ewIZk33MiGH6tjPHnphKosNGEppY5YYno6aW9NpVbvXmz5tSorPq/L/l3a7F1FhiYspVSp+BITaThsGE1O38L+HXEs/6QeO1YneR2WikGeJiwRGSciG0XkFy/jUEqVXvVme0m/cBMJ1bJZ801t1n9fnVytIlRlyOsS1gSgk8cxKKXKSELVHNLO3UztY3ax9feqrPi0Lvt3ahWhKhueJixjzFdAxF72pZQqf+KHBm130PTMLLJ3O1WE21cmex2WigFel7CKJSKDRGSBiCzI1psUlYoa1ZrsI/3CjSTWPMDaubXoE1efeMTrsFQUq/A3DhtjxgBjAFJSUvT93UpFkfiUXFI7ZrFpUTXOsGuQ5ksie+9W4pL0PVuq5Cp8wlJKRTfxQf0Td/Kvn3YwML4RKz+vQ/NzsohP1qRVUdgBqxPwLOAHxlpBe2QB41wJZAAG+MkK2r3d/jnAIne0VVbQ7hKpOCt8laBSKjYszt3D8/vXkL3Hz8rP6nJgtx5+KgI7YPmBF4CLgNZALztgtc43TivgPuB0K2gfC9wWMvhPK2i3cf8ilqzA+2btU4C5wDEikiki13sZj1Iqsv4we2nWIYucfT5W6k3GFUV74A8raC+zgvZ+YCrQNd84A4EXrKC9FcAK2hvLOUbA4ypBY0wvL5evlCp/VeoeoPk5Waz6og4rP69D6jlZJFTTG7YiLE5EFoR0j3HbBwA0AVaHDMsETs43/dEAdsD6BqfaMMMK2h+7w5LsgLUAyAZGWkF7RlkHn0fL5Eqpcpdc+wCpHTdjsoWVn9dl3w4taUVYtjGmXcjfmOInOUQc0AroAPQCXrUDVk13WKoVtNsBvYFn7IDVsqyCzk8TllLKE0m1smneMQuTi5O0tmsbMI+sAZqFdDd1+4XKBGZaQfuAFbSXA7/hJDCsoL3G/b8MmAO0jVSgmrCUUp5JqplN6rlZCLDy8zrs3apJywPzgVZ2wEq3A1YC0BOYmW+cGTilK+yAVReninCZHbBq2QErMaT/6cCSSAWqCUsp5anE6tmknrsZ8cHKL+qSvk5vtyxPVtDOBgYDnwA2MM0K2ovtgDXcDlh5rf4+AbLsgLUE+AIYagXtLMACFtgB6ye3/0graEcsYenpjFLKcwnVckg9dzOrPq/DsCk5PHyVnz+a6FMxyosVtGcBs/L1Gxby2QB3uH+h43wLHF8eMYKWsJRSFURC1RxSz81iZzLcPzWHY1ZrSUsdShOWUqrCiE/J4cGr/WyrCv98K4fWK/VpGOovmrCUUhXK1mrCg338bKoB903L5fjlmrSUQxOWUqrC2V5V+FdvP+trwT1v59JmqSYtpQlLKVVB7UhxklZmPRj6Ti7/97smrcpOE5ZSqsLaVUUY3svPivpw57u5nBzUpFWZacJSSlVoe5KEEb38/NEIbpuRy+mLNWlVVpqwlFIV3p+JwsM9/QSbwZAPcjn7Z01alZEmLKVUVNiXIDx6pZ9fUoUbP8ql40JNWpWNJiylVNTYHy88doWPn1oIN/w7lwu+16RVmWjCUkpFlQPxwhOX+5jfShjwn1w6z9OkVVlowlJKRZ3sOOGp7j6+O0bo/1kuXedq0qoMNGEppaJSjl94ppuPbyyhz5xcLv+vJq1Yp09rV0pFrVyfMLqLj2x/Lld9nUtcjuGts3wg+qT3WKQJSykV1YxPePESJ2ld/q0hLjeXyR208igWacJSSkU9I8KYi3xk+3Lp+p0hLjsX088gWtKKKXoaopSKCUaE1y708dHfhIsXGNYPH47J1etascTThCUinUTkVxH5Q0Tu9TIWpVQMEGHiuT7eP0XYNmUq64YN06QVQzyrEhQRP/ACcD6QCcwXkZnGmCVexaSUigEiTO7g4/qTBrH5xZfgQDaNHnkY8fu9jqzCsgNWJ+BZwA+MtYL2yALGuRLIAAzwkxW0e7v9+wH3u6ONsIL2xEjF6WUJqz3whzFmmTFmPzAV6OphPEqpWCFCvVtuoe4tQ9j+/vusHXo35sABr6OqkOyAlVd4uAhoDfSyA1brfOO0Au4DTreC9rHAbW7/2sCDwMk4x/QH7YBVK1KxetnoogmwOqQ7E2elDyEig4BBAAkJCUe8sBUjLz7iaQvS4bsnAJhTxvONiIztRQ+f08Edb06Roy0qk2BKKZbWBWJrfSrgutS76SZ8CQlsHPUkJjubJqOeQMI9jlTA9YmQ9sAfVtBeBmAHrLzCQ2ht10DgBStobwWwgvZGt/+FwGwraG9xp50NdAKmRCLQCt9K0BgzBhgDkJKSYjwORykVZeoMGIDEx7Ph0ZFk3nqAJs8+g68UJ78xKJzCw9EAdsD6BqfaMMMK2h8XMm2TSAXqZZXgGqBZSHdTt59SSpWp2v360WDYA+z64gsybx5M7t69XodU3uJEZEHI36CSTg+0AjoAvYBX7YBVs4xjDCsIr8wHWolIOk6i6gn09jAepVQMq927NxIfz/phD7L6xhtp9uKL+JKTvQ6rvGQbY9oVMiycwkMm8D8raB8AltsB6zecBLYGJ4mFTjunLAIuiGclLGNMNjAY+ASwgWnGmMVexaOUin21evSg0aOPsOd/81g96O/k7t7tdUgVwXyglR2w0u2AlYBTeJiZb5wZuInJDlh1caoIl+Ecvy+wA1Ytt7HFBW6/iPD0PixjzCxjzNHGmJbGmIe9jEUpVTnU7NaNxo8/zp4ffmDVgIHk7NzpdUiesoL2YYUHK2gvtgPWcDtgdXFH+wTIsgPWEuALYKgVtLPcxhYP4SS9+cDwvAYYkVDhG10opVRZq3HJxUh8PGvuvJNV1w+g+atj8Neo4XVYnrGC9ixgVr5+w0I+G+AO9y//tOOAcZGOEfTRTEqpSqr6hRfQdPSz7LNtVl57Ldlbt3odkiqGJiylVKVVrWNHmr7wPPv/WMqqfv3JzsryOiRVBE1YSqlKrepZZ9Hs5ZfYv2oVK/v148DGjcVPpDyhCUspVemlnHYazca8woG161h1TT8ObNjgdUiqAJqwlFIKSGnfnuZjXyV70yZWXt2XA2v0OQYVjSYspZRyVTnpJJqPe42cbdtY2fca9q9eXfxEqtxowlJKqRDJJ55I8wnjyd2920laK1Z4HZJyacJSSql8ko89luYTJ2D27WNl32vYt3Sp1yEpNGEppVSBkgIBUidNxBjDymv6sfe337wOqdLThKWUUoVIbNWK1EmTEL+fVdf0Y+8SfSG6lzRhKaVUERJbpJP6+iQkOZmV/a/lz0VR8ErGGKUJSymlipGQmkrq66/jr16dVddex54ff/Q6pEpJE5ZSSoUhoWkTUl+fhL9ObVZfP4Cj9u/3OqRKRxOWUkqFKb5RI1InvU5cw4YM3rqVVvs0aZUnTVhKKVUC8Q3qkzppIkvjE9jm10NoedKtrZRSJTTi+ee5YtHPTP/mG0SEjIwMr0OqFPQFjkopVUIZGRmapDygJSyllFJRQROWUiriMjIy+PLLL/nyyy+1Ck0dMa0SVKqCyjvIA4gIDz74YNQe6LUKrWKzA1Yn4FnAD4y1gvbIfMP7A08Aee9ced4K2mPdYTlA3t3Uq6yg3SVScWrCUjEjlg7woAd5VT7sgOUHXgDOBzKB+XbAmmkF7fzPoXrLCtqDC5jFn1bQbhPhMAGPqgRFpIeILBaRXBFp50UMpRFL1Ruxti7GmIN/0bwuSpWj9sAfVtBeZgXt/cBUoKvHMRVIjDHlv1ARC8gFXgHuMsYsCGe6lJQUs3v37ojGppRSsUZE9vNXtR3AGGPMGAA7YF0BdLKC9gC3uy9wcmhpyq0SfBTYBPwG3G4F7dXusGxgIZANjLSC9oxIrYcnJSxjjG2M+dWLZSulVCWUbYxpF/I3poTTfwCkWUH7BGA2MDFkWKoVtNsBvYFn7IDVsoxiPkyFv4YlIoOAQW6nEZE/vYznCMThnHnEAl2XiiuW1ieW1gUqxvokFzFsDdAspLspfzWuAMAK2lkhnWOBx0OGrXH/L7MD1hygLRCRN15GLGGJyKdAwwIG/dMY836483HPBEp6NlBhiMgCY0zUXacriK5LxRVL6xNL6wJRsT7zgVZ2wErHSVQ9cUpLB9kBq5EVtNe5nV0A2+1fC9hjBe19dsCqC5xOSDIraxFLWMaY8yI1b6WUUmXDCtrZdsAaDHyC06x9nBW0F9sBaziwwAraM4Fb7IDVBaekuAXonzc58IodsHJxLjGNLKB1YZnxpNHFwYWLzKEEjS6iURScXYVN16XiiqX1iaV1gdhbHy951ay9u4hkAqcCH4nIJ17EUU6itjqzALouFVcsrU8srQvE3vp4xtMSllJKKRUufZagUkqpqKAJSymlVFTQhBUhItJJRH4VkT9E5F6v4ykNERknIhtF5BevYyktEWkmIl+IyBL38WC3eh3TkRKRJBGZJyI/uevyL69jKi0R8YvIjyLyodexlJaIrBCRRSKyUERitmFZedJrWBEgIn6cx5ccfJgk0MsYE7HmnpEkImcBu4BJxpjjvI6nNESkEdDIGPODiFQDvge6ReN3IyICpBhjdolIPPBf4FZjzHceh3bEROQOoB1Q3RhzidfxlIaIrADaGWM2ex1LrNASVmS0B/4wxiwzxlToh0mGwxjzFc69F1HPGLPOGPOD+3knzg2QTbyN6sgYxy63M979i9ozUBFpClyM8yQFpQ6jCSsymgCrQ7ozidKDYiwTkTScx8j8z+NQjphbhbYQ2AjMNsZE7boAzwB34zwYOxYY4D8i8r37iDlVSpqwVKUkIlWBd4DbjDE7vI7nSBljcowxbXCe/9ZeRKKyylZELgE2GmO+9zqWMnSGMeYk4CLgZrdqXZWCJqzIKPZhkso77vWed4DJxph3vY6nLBhjtgFfAJ08DuVInQ50ca/7TAU6isgb3oZUOsaYNe7/jcB7OJcKVClowoqM+UArEUkXkQSch0nO9DgmxcGGCq8BtjHmKa/jKQ0RqSciNd3PyTiNfIKeBnWEjDH3GWOaGmPScPaXz40xV3sc1hETkRS3UQ8ikgJcAER9K1uvacKKAGNMNpD3MEkbmGaMWextVEdORKYAc4FjRCRTRK73OqZSOB3oi3MGv9D96+x1UEeoEfCFiPyMc5I02xgT9c3BY0QD4L8i8hMwD/jIGPOxxzFFPW3WrpRSKipoCUsppVRU0ISllFIqKmjCUkopFRU0YSmllIoKmrCUUkpFBU1YqsIQkYYiMlVElrqPs5klIkeXcwwdROS0kO4bROSaMpivT0RGi8gv7hO854tIemnnW8wyh4vIee7n20SkSiSXp1SkabN2VSG4N/R+C0w0xrzs9jsR56ndX4cxfZx7/1uB3SWIIwPYZYwZVdJpi5lvL+By4EpjTK77oNfdxpitZbmckOX5jTE5Id0r0CeHqyinJSxVUZwDHMhLVgDGmJ+MMV+L44mQ0slVcLA09LWIzASWFNCdFvoOLxG5y01IiMgcEXnWvXH4FxFp7z4M9wbgdrf/mSKSISJ3udO0EZHvRORnEXlPRGqFzOsx991Uv4nImQWsXyNgnTEm1123zLxkJSIXiMhcEflBRN4WkarivE/t7ZDYO+S9I6qg8d3+K9w4fgB6iMgEEblCRG4BGuPcZPyFiFwnIs+EzHugiDxdqm9PqXKgCUtVFMfhvJuqIJcBbYATgfOAJ8R5rxXASTjvgDq6kO6iVHEfHHsTMM4YswJ4GXjaGNOmgJLdJOAeY8wJwCLgwZBhccaY9sBt+frnmQZc6ibCJ0WkLYCI1AXuB85zH5S6ALgD+BQ42X2sD8BVwNQixs+TZYw5yRgzNa+HMWY0sBY4xxhzTkgs8e4o1wLjwtheSnkqzusAlArDGcAUt4prg4h8CfwN2AHMM8YsDxk3f3dRpoDzvi8RqZ73XL6CiEgNoKYx5ku310Tg7ZBR8h6i+z2Qln96Y0ymiBwDdHT/PhORHkAy0Br4xqkVJQGYa4zJFpGPcRLLdJz3RN0NnF3Q+CGLequ4lXZf+Pg5cImI2EC8MWZRcdMp5TVNWKqiWAxccQTT7S6iO5tDaxGS8o2b/wJuaS7o7nP/51DIfmWM2Qf8G/i3iGwAugH/wXkGYK8CJpmK80zKLcACY8xO91pfYePD4dujMGOBf+A8LHd8mNMo5SmtElQVxedAooS86E5ETnCvB30NXCXOywrrAWfhPFC0OBuA+iJSR0QSgfyvXM+7FnYGsN0Ysx3YCVTLPyN32NaQ61N9gS/zj1cYETlJRBq7n33ACcBK4DvgdBE5yh2WEtIy8kucKs6BOMmLYsYvyiHr5b7osRnQG7ekqVRFpwlLVQjGaa7aHThPnGbti4FHgfU47xL6GfgJJ7HdbYxZH8Y8DwDDcZLbbA5/9cZeEfkR57pV3hPoPwC65zW6yDd+P5zrZz/jXFMbXoJVrA984DYC+Rmn9Pe8MWYT0B+Y4s53LhBw488BPsR5AeCHbr9Cxy/GGOBjEfkipN804JtItVRUqqxps3ZVKYnIHOAuY8wCr2Pxitvq8GljzGdex6JUOLSEpVQlIyI1ReQ34E9NViqaaAlLKaVUVNASllJKqaigCUsppVRU0ISllFIqKmjCUkopFRU0YSmllIoK/w91/Pcc4raHnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_labels = ['0', '1', '2', '3', '4', '5']\n",
    "\n",
    "N = 6\n",
    "ind = np.arange(N)\n",
    "width = .35\n",
    "\n",
    "fig_bar_cd, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "\n",
    "p0 = ax.bar(ind[0], score_x.mean(), yerr=score_x.std(), capsize=2)\n",
    "p1 = ax.bar(ind[1:], mu_noharm, width, yerr=std_noharm, capsize=2)\n",
    "p2 = ax.bar(ind[1:] + width, mu_harm, width, yerr=std_harm, capsize=2)\n",
    "\n",
    "ax.set_title('Adversarial Scores and Accuracy by Corruption Severity')\n",
    "ax.set_xticks(ind + width / 2)\n",
    "ax.set_xticklabels(plot_labels)\n",
    "ax.set_ylim((-1,6))\n",
    "ax.legend((p1[0], p2[0]), ('Not Harmful', 'Harmful'), loc='upper right', ncol=2)\n",
    "ax.set_ylabel('Score')\n",
    "ax.set_xlabel('Corruption Severity')\n",
    "\n",
    "color = 'tab:red'\n",
    "ax2.set_ylabel('Accuracy', color=color)\n",
    "ax2.plot(acc, color=color)\n",
    "ax2.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "ffba93b5284319fb7a107c8eacae647f441487dcc7e0323a4c0d3feb66ea8c5e"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}