examples/cd_mmd_cifar10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum Mean Discrepancy drift detector on CIFAR-10\n",
    "\n",
    "### Method\n",
    "\n",
    "The [Maximum Mean Discrepancy (MMD)](http://jmlr.csail.mit.edu/papers/v13/gretton12a.html) detector is a kernel-based method for multivariate 2 sample testing. The MMD is a distance-based measure between 2 distributions *p* and *q* based on the mean embeddings $\\mu_{p}$ and $\\mu_{q}$ in a reproducing kernel Hilbert space $F$:\n",
    "\n",
    "\\begin{align}\n",
    "MMD(F, p, q) & = || \\mu_{p} - \\mu_{q} ||^2_{F} \\\\\n",
    "\\end{align}\n",
    "\n",
    "We can compute unbiased estimates of $MMD^2$ from the samples of the 2 distributions after applying the kernel trick. We use by default a [radial basis function kernel](https://en.wikipedia.org/wiki/Radial_basis_function_kernel), but users are free to pass their own kernel of preference to the detector. We obtain a $p$-value via a [permutation test](https://en.wikipedia.org/wiki/Resampling_(statistics)) on the values of $MMD^2$. This method is also described in [Failing Loudly: An Empirical Study of Methods for Detecting Dataset Shift](https://arxiv.org/abs/1810.11953).\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 os\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from tensorflow.keras.layers import Conv2D, Dense, Flatten, InputLayer, Reshape\n",
    "\n",
    "from alibi_detect.cd import MMDDrift\n",
    "from alibi_detect.cd.preprocess import UAE, HiddenOutput\n",
    "from alibi_detect.models.resnet import scale_by_instance\n",
    "from alibi_detect.utils.fetching import fetch_tf_model\n",
    "from alibi_detect.utils.kernels import gaussian_kernel\n",
    "from alibi_detect.utils.prediction import predict_batch\n",
    "from alibi_detect.utils.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 MMD 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": [
    "X_c = []\n",
    "n_corr = len(corruption)\n",
    "for i in range(n_corr):\n",
    "    X_c.append(X_corr[i * n_test:(i + 1) * n_test])"
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAYt0lEQVR4nO2dWYxl11WG/33nc6ea56qu6rntdntI0g7uhITECQKhKEqEEAIEeQUpDyjwQBBR8ghCvIBACN4IUkAWg5AiSGJQEhPbwnZijNvtHmvqmqtu3brzuHmoEmpZ+992msa1bf5PihSd1fvec869/z3l9e+1lrHWQggRHrHjPgEhhBuJU4hAkTiFCBSJU4hAkTiFCBSJU4hAkTjfQxhjvmyM+csH/W/fwWtZY8yZB/Fa4p1j5HMeH8aYLwD4EoDTAA4A/D2A37HW7h/neb0VY4wFcNZae/O4z+X/E3pyHhPGmC8B+H0Avw1gAMBPAJgH8G1jTMrx7xPv7hmK40biPAaMMUUAXwPwRWvtP1trO9baRQC/AGABwK8YY75qjHnGGPN1Y8wBgC8cHfv6Pa/zq8aYJWPMrjHm94wxi8aYTx3F/uffGmMWjv40/TVjzLIxZscY87v3vM6TxpjnjTH7xph1Y8yfuH4gxLuLxHk8XAGQAfB39x601lYBfBPAp48OfRbAMwAGAfz1vf/WGPMwgD8F8MsApnD49J15m/f9KIDzAJ4G8BVjzENHx3sAfhPAKICnjuK/cR/XJR4gEufxMApgx1rbdcTWj+IA8Ly19h+stX1rbeMt/+7nAfyTtfY5a20bwFcAvF0C4WvW2oa19lUArwJ4DACstS9ba1+w1naPnuB/DuDj93dp4kGh/445HnYAjBpjEg6BTh3FAWDF8xrT98attXVjzO7bvO/GPf+/DiAPAMaYcwD+CMCHAGRx+L14+e0uQvzfoifn8fA8gBaAz9970BiTB/CzAJ49OuR7Eq4DmL1nbQRg5D7P588AXMNhRrYI4MsAzH2+lnhASJzHgLW2jMOE0B8bY37GGJM0xiwA+FsAqwD+6h28zDMAPmOMuXKUvPkq7l9QBRxaOVVjzAUAv36fryMeIBLnMWGt/QMcPqH+EIfCeBGHf6Y+ba1tvYP1rwP4IoBv4PApWgWwhcMn8o/LbwH4JQAVAH8B4G/u4zXEA0abEN4nHP1JvI/DP03vHPf5iP89enK+hzHGfMYYkzXG5HD4BH4NwOLxnpV4UEic720+C2Dt6H9nAfyi1Z9C7xv0Z60QgaInpxCB4t2E8MlPf5A+VpOFIl23VdpzHt/b48UWrUqTxoYm+Xslhrm1Z5LktyfOf5M6lQ6NLb9ylcaSRb4Vde7sFI1FCbf70e8k6ZpelzsmQ2MZGps6ye9VPOH+KvR7rk1MhySS/JoP9vh93N7YorFO331tTz15ga6xLX6O3/rW92hsZoHvdoySaRpbW9lwHo9HebqmkCvQ2He+8azzovXkFCJQJE4hAkXiFCJQJE4hAkXiFCJQJE4hAsVrpSTyPC0fjQ3RWL7l3nu9VyrRNcMTPNU8eZpbEfvNPo3RIg1iGwBAvVmlsV6f2wMDxQEaGxvn15awbjvioNyja/pxfo750SyNdXr8NVsNd6zXadM16ZyvCCbOz6PF72MiFTmPjwxwO61eLfPYQZ3Gttd4+WuU4jZR3LqvLVccpGva5P760JNTiECROIUIFIlTiECROIUIFIlTiECROIUIFL+V4klfJ9M81Zwvuq2D3B5fMzE7TGNRIUdj5Ta3FRIJUtkR45fda7y1Pew9r+f5Kct5bKdOl1dNxKzbOmjWDuiaZpvH+t1RGmuWeeXP3oa7Yiie4tUxYyfc5w4AiRS3WVo1bs9kIvdnnUl7qnSa3LZp1rlt067zWuaJEf59zBTd1Scdz7NufWmNxhh6cgoRKBKnEIEicQoRKBKnEIEicQoRKN5s7cDYOI1V9t19ggAgk3dvvi4M8R4rg1M8O1b19DBPxngWL0M2UXf6fLN8t8kzmilP5tJ0eeavtMEzyhny89iqVugaGL6JOhvnWeNCjt//fsd9Ih3Ds66s7xAA9Ls8SxqL89dMJt33OB7jGdkoza95cm6axmbn5mlsaoZ/91skE726uErX1Bu86IOhJ6cQgSJxChEoEqcQgSJxChEoEqcQgSJxChEoXislneDpa+OJjU+609cHrR3+ekl+Kq0y91JSMd42P9l3//b4hje123xTtq9bUXmHp8qjHN+438y4bZHBEd6PJl/g1kHFcpul3uU2US/rvo+mzTftN8q8P08qxX/3TZLf/yyx4dIxvsm+OM77Jl14nI9xgOc7bCN+jjEyziMbcavtA1ce5efB3ufHXiGEeFeQOIUIFIlTiECROIUIFIlTiECROIUIFK+VUinzNvfGU9mxsrzkPJ7zTAuu7/K+OL0Otw5SnqqU2r67L04sy3sZeaspPBUaKU+Pm5ET3BbJDbrHOGQLvIIEMf6b2utwC6DjKe8x1n1t1S1efVTe5uMMHr58nsZGJvkoD5DTTyf5d2CwyK2q3DDvg9Xo8c+64zHOhvLuz3Nojn+/K1VemcTQk1OIQJE4hQgUiVOIQJE4hQgUiVOIQJE4hQgUv5VS41UHnRhPNS/+6DXn8Zl53myp4KncGMzxqgPraf5VLtfcAY9d0vdUYeQ953jysRM0NnpmhMbipMLBGP67ubnELa6VN3iTqeECtzAuPnLJefyl1922GADs73B7IFfgk75jcW73tFruqprsIJ8Onklz2ymX80xnt3yd6fFzHB0ccx5/7fVX6Jo3r16nMYaenEIEisQpRKBInEIEisQpRKBInEIEisQpRKB4rZR6k095bve55dCy7lhumlsKUZ/v6O+1uV8SM7xJUz7jTqNv7/FmXM0Gf6/TjyzQ2MITMzTWsrxpGHNMKmvcLrn+g/+isWrZY2+c90yHhvu6i+N8Zkja89OejvHKnw7/qFGYcTfy2mrx6phCntssuYjbcIk+P0d0uVXYI3Nlbl9foWs2b23x9yLoySlEoEicQgSKxClEoEicQgSKxClEoHiztRFpjQ8A1R0+WmFyZtZ5fOH0KbpmKOKbspdv3aGxtdt8Y/bwmDuLlySZSQBoT/IN27MXJmksluS/c7GmZ6xF19275/bLfAN7bY9s6Adw/lF+jy98+CEaW192ZxqLnpTshcvnaCxW5JnhaJBn7ZNZ9/s12+5+UACwuce/pwY8IxuP8Z5QPc8k7UrF7WJsb/GeSv0+30jP0JNTiECROIUIFIlTiECROIUIFIlTiECROIUIFL+VMsw3FKdKfGN2DO40dD7DW+NHRZ4OP/UQb+2/sbzBY5vu1PZknveVefxRbjfMkYndAGDJFG0A6MZ4z6Ibr990Ht9e3qZrJk66e9gAwIUPX6Sxwgi/x42Ge+p1scB3qacnhmkslvRsfAcvmti86b7uuXMTdE2jyws0EjGPheHbnN/nNsvO9przeGmX24tRjN97hp6cQgSKxClEoEicQgSKxClEoEicQgSKxClEoHitlEyCVxYkPanmbsfdUr/f431ZjKdCIPK01D99kdssL3/vRefxa3fv0jWXPsqtiFaSp+WTZX5tI5affwXuKckXz52la0bPclshmePWR63Oq1nG5t3nkRrg597gDhGGI17VcetH3P5aXXb32vnoBfe4CADox9w2EAD4ikFsjI9j6PS4VdjvuMeU9Hvu7z0A9A2PMfTkFCJQJE4hAkXiFCJQJE4hAkXiFCJQJE4hAsVrpUzE+U76xTpPX/d67qqDTouPJeh1eao5luZp+dlzCzS2vuhu/rWxw22P9LR7HAAA7HYPaGy8zM+/0ONNw4Yidzr/zCeepmuGp3k1SLnBLYCq4SMNWj13ZUdqzWMP1Ph9rEZ8KnrSM0LjzBNuaywzyiukdnf5eI16xzOuI8Vj6TivnMmQZTHD7cBqtUJjDD05hQgUiVOIQJE4hQgUiVOIQJE4hQgUiVOIQPFaKdUST//WqrzCgWXKyyVuRVjPjv7xOc+MkohXTTzy1GPO45eap+maeJyXWjR2uE0xkeLVINkeT7Gj5J5EvXHb3fgLAOJxPkW76GkkFe/xe9XquG2RVInPlUkl+HvtrHF744xnEnUL7vvYrHDrLuGpnjqo8fklLcs/68lBfm19cq8SKS6n6QnelI2hJ6cQgSJxChEoEqcQgSJxChEoEqcQgSJxChEoXivFZLk9MDnLm0y1Wu60d6/Dd/q3mzxlX9rgc0PGF+ZobGjEXb2R2+OX3Vpxz8EAgJkUn/XSifF5HW3DU/bT0+7X7JB0PQB0VtxNsABgu8M7WvXjvAqjkHNXx+QiXlGTSPFZIzHPHBLfKPudXbdd1V7kNpYd5hZR1nOO8cjzbEpye6ZFuoYtnD9F15w8we0vhp6cQgSKxClEoEicQgSKxClEoEicQgSKfxzDYI7GUjs8CxYV3dmzVIK/XSLOY6U13r5/fIpviu/F3RvOuwc8M9wp8d43Wz3eAymZ4ZntomeSdoYkBbMFnhlu1nnWu+Xp7eQrLmA9bqoJ/npxz4ZzePpPpUaGaGxuwJ1h7/f5vb/55iqNDU2M01grybPX1QZ/vziRTZTm3+G25a/H0JNTiECROIUIFIlTiECROIUIFIlTiECROIUIFK+VUqtxW6Hb5pu5uyRD3e1zC6DX4xu2E1k+IqF+wPscZQbcm7kTRd7D5spPfZzGXnzlFRr795d+SGOXPFOqJ4bc51LZdfcWAoCBQb4ZfXZiisYaNf6au/vuUQ1Nj6WAOP/MNne5/ZUtcBtu/ox7HINp8u/OyT4vEljc40UCieI0jdWa/LoXb9xyHr9z/RpdM7XwERpj6MkpRKBInEIEisQpRKBInEIEisQpRKBInEIEitdKaTd4X5xc1m1TAEAHbpuln+Ep76jIXy+b463s2RRtAOiTKoy7Zd6i/2yW2yxPXvoAjb38ylUaq7f4OUakR0/GM3U5FuPjHdbWNmksneZVJPMLC87jts/fK+mp6pjzjOtY95zjzTfc9/HcxSfomtPDF2ls70Xef2rPU4HUAb+23QN3P6OBoVG65tRpPgKEoSenEIEicQoRKBKnEIEicQoRKBKnEIEicQoRKF4rJQ5edZDNc+ujOOKOtfqeKckpT4v+1XUay426G0IBwMGae10mxS2FF67yyoKPPHaZxj73+c/R2OrSIo31SHVPpsAtHXgGZRfy/CPt9Xkl0dqqu4okleIVQf0uf71ExO/xxCy3xsq7bgtmZ4M38bpZ5hPTpyYXaGx1Y5HGbJ5Xzpw4f8J5fPHqHbpmY3WHxhh6cgoRKBKnEIEicQoRKBKnEIEicQoRKBKnEIHitVKyEU+jd3s8nz807N6dH2vx1HuzzWdybN31zMLgbg+6HXfzr2iKz8/YS/J5Ij94lTfx+rlP/jSN2Sav7lm+ddN5PB1xq6rV5s2npid5ZUTaM8tjv+Ju/pVJ8Rkwpsc/z80Stw56nsnWUc49V6ZR43ZJp8WrS777wxs0tljnzeHyg9wKGhhx62L2/CxdMzrBJ8Ez9OQUIlAkTiECReIUIlAkTiECReIUIlC82dpogE9X7llfjxt3pmttiW8Mbud49ref4LHNZZ7JnV1wZ8jaDZ4ZHp7hmdyrz/+IxnLf+z6NPfEIH8fQbLizpClPj6bRSb4pvl1397cBgHabFx6MDo84j/eNr18RH7nQa3t+99v8Nbvk/Xp9nkWP0nyT+soWH8cQG+GZ7b2dEo119/edxz/wMT5yYXJU2Voh3jdInEIEisQpRKBInEIEisQpRKBInEIEit9KyWdprNLkqe07b7o3c9c8m6FzWd5XpsNdG9QavO1/POneRH17cZmuOdjjm6FnLp2hsW8++xyNVVp80/aTly45j7eafFN5Nuu+LgBIJflHWiYWAMDtpchj6cSSvDAiHXlGb8T5ObaJZdLq8PvR8ozkmDvFxyBUE3xzeznGKyqGJsh3Nc2LBDabfAQIQ09OIQJF4hQiUCROIQJF4hQiUCROIQJF4hQiULxWSjrBU8Pr2ys0tnTtTefxS5f5BOJ4gvsllR5Py+cH3JOhAaDZcPfaGRnmIxyWV5ZobOrcPI2d/ODDNHZzkVfOnFpwt/Y/Pc/fq+mZGt3tcQtgfHKGxtZW3dddOuDWUgr8c+l6Rj+UPHZVOuv+ztk+t0tsl9t6qQyvgKl5JpzPnnR/LgAw/7Dbnrlb4hZdtcn7PjH05BQiUCROIQJF4hQiUCROIQJF4hQiUCROIQLFa6WU93k1RbXMKxzyWfduf+NJh6fT3AIYHuJVGOs7fNRBjTS0WjjN0+QDY0M0duvGLRq7MM+rH2IJXt3Ttu4Ue73J7ZIiub8AUOny5mXtDo9li4PO4zv7vEFWo8SbYBUL3OLKJvkzIWbctshQjlfAVHruJmkAkKvxUQ2DniqSgQne6G27te08Xu1yiwiWNyFj6MkpRKBInEIEisQpRKBInEIEisQpRKBInEIEitdKqXumCWfTPJ1/5VOfcB6/8NApumZll9sUqwe8YqVxg1spjbrbjqh0uKUzlnfPDAGA3T5vUPbG69do7GMXH6Ox0bx7Hk1ll1dMFD1VNabL56GU657KCOP+KsR44QlyOT6zJZvh1odvSnWazD3pG24D1dOe72mdX8CpKV6ls5vg71cqu78HyYhbM90Gr5xh6MkpRKBInEIEisQpRKBInEIEisQpRKB4s7XDkzwrOHX2HI09TnrtDI3yzdDFYZ79TfEkKRJ53iNmd9Odle33+Qbl5aV1GhvM8vNPjk3S2FaDv99cLuc8Hu/yQoBek2dku57p1T14xjiQEQkpw3+/G12e9Z4a99wPvpce1Zr7Xu177mHT8u9AY5+f43aD93aynknUpu3uj5TOeUZXpHlPJbrmx14hhHhXkDiFCBSJU4hAkTiFCBSJU4hAkTiFCBSvldKo882/q9W7NNbubDqPz588SdfMTozS2Pnp8zQWj/FLiFJ7zuOtFt+E3Krwaz4o83T4o+e4tZTx9PzZ33JvcB9LcNtjdZt7S3c9G+Zt0m3bAMCpSbd1UMjyDewm7ilIaPNN9okY76dTrbotk65nsvVEnvf7uVq7QWOv37lDYyfnPZv6U+7Ps0OmgwPAyhIf1cDQk1OIQJE4hQgUiVOIQJE4hQgUiVOIQJE4hQgUr5Wyu8FT9l3PNOGr19xp45Ob3H658tRlGhsd5Lv950dnaSwec6f6VzwjBuYe4mn5rVU+fuDmzf+gscEhXqFRtO7qkwpvjYTlZV5N8eYSnzg+PsKvbTTrtjfGBnlPpaFBd/8jAFhZ59ZB0WPPDA67x0LUanykxfaB2zIDgL0aH9VQ9kzthvFUupDv/sbtm3RN1OdVRgw9OYUIFIlTiECROIUIFIlTiECROIUIFIlTiEDxj2No8MqCYoantm8suif/Lt9xV6sAQPWAT3K+fOVhGhse4pOoJ0fdE6xzEW/UtVxapLH+LK/qqGb4+R/UuL3RzbirTyp9Typ/jFdMJBJzNFaqcluhywpMiNUDAAclPt18ZII3yGpUyzRWKrtjsQSvZLm7yy2/V27yypPRx/l4EF9js9XrbisrT+woAEhZNfgS4n2DxClEoEicQgSKxClEoEicQgSKxClEoHitlCjLm0yhy22WWM9tA2xu8OZTz/7jczRWHOCNpM5eOkNj2YS7amK2MEbXpD2jnN/s82oQM0VDSLW4HWFb7vvYyXgaWo3y6pLxLj+R2h6fAF0h55G3vHKj3uYNrRIRtxVyaT4BukSsmzurt+maa4u8GgSeCpjxGV7R9J/ffZHGPv6hDzmPX/7Jp+ia7//rt2iMoSenEIEicQoRKBKnEIEicQoRKBKnEIEicQoRKF4rJZnj2vVMHEdyyF2xMj/IG12tvrFBY899+1UayxZ5qjybc1tBuYhf1/gAr1RIZnmzq6Udns4/qHNbpBm5m0WVyu7KHgCotHmsucUrPrJ1bo11+sPO4/sZbi2l0rw6pt3m60pV3pDrLqlY2UtyO6pX4Nc1OcK/H9t3lmgs4Tn/E2fcDefiCW4VDuZ5JRRDT04hAkXiFCJQJE4hAkXiFCJQJE4hAsWbrbX9Oo3t7/KeOet33dnEhz68QNe0azwbt7/LN1//27+8RGPdmDsT2j7HU83THR4bKfJs7fnJizRWqvAM6lbd3f8mDj7uIhvj/ZtaKfc4AwC4/sOrNLa+5R5RMTV7mq7Zu32LxtpNPk/CgPdHisbd53/iYT7dfOiEu1cUANSavG9SLMGfTSNTvLjARu7vyH6Fa2L/wDNfg6AnpxCBInEKESgSpxCBInEKESgSpxCBInEKESheK2V/k09yvvbydRpr1lrO43EyegAARua4BdBuuF8PAO7e4K34X4B7w3wyStI1B2N8U3Zxj5/j9DjfMD9YGKWxVNL9+5g1vAfPWJa/3tgCt1nmB/hG9e++4Lak7tR4QcJOjU8qH/EUOcycmKex2Vl3D6S5aT5mYmeXf0+r4H2OAG7fFQp8zEerTyyTHr/34zO85xZDT04hAkXiFCJQJE4hAkXiFCJQJE4hAkXiFCJQjPVMLhZCHB96cgoRKBKnEIEicQoRKBKnEIEicQoRKBKnEIHy311ehR8F48kKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAcoklEQVR4nO2deXDd1X3Fz9WzdmuzVluSJVmS5VW2seUFY4MJNhgIe1kKAdKWFNI13ZdJ006TNE0nacpkBpi2JA2GAAHMYsALAe/yvi+SrX3f99WSdfuHHh3XvecVOk24Sc5nJhNxj77v/X7vvaP7fL/3fr/GWgshhH+EfdYXIIRwI3MK4SkypxCeInMK4SkypxCeInMK4Sky5y8IxpiHjTHbP+vrYBhjBowxsz7r6/hlwijPKYSfaOYUwlNkzv8FY8w1xpjjxph+Y8xPjDGvGGO+boxJMsZsMca0G2O6gz9nXRFXY4y56Yr//ltjzKbgz1HGmE3GmE5jTI8x5rAxJj2oPW6MqQo+X7Ux5uErxvde8Xj/YoypN8b0GWOOGmPWXPVcrxpjfhR8nLPGmGWf4F5rjDF/Yow5ZYzpDd5r1BX6E8aYCmNMlzHmbWPMjCs0a4wpCP58qzHmXPC5G40xf3LF791ujDkRvO/9xpjiT/+u/Gogc4bAGBMBYDOAHwKYBuDHAO4OymEAfgAgB8BMAMMAvv8JH/oxAAkAsgEkA3gSwLAxJhbA0wA2WmvjAFwL4AR5jMMAFgev6yUAP7nSSADuAPAygEQAb3+Ka7sfwC0A8gAUA3gcAIwxNwL4h6A+HUBt8PFd/DuA3w7ewwIAHwYfYwmA5wH8dvC+nwPwtjEm8hNe268UMmdoVgKYAuBpa+2YtfYNAIcAwFrbaa193Vo7ZK3tB/ANANd/wscdw+SHs8Bae9lae9Ra2xfUJgAsMMZEW2ubrbVnXQ9grd0UvIZxa+13AEQCKLriV/Zaa9+z1l4G8AKARZ/w2p621jZZa7sAvIPJPwAA8DCA5621x6y1owD+EsAqY0wuub95xph4a223tfZYcPxLAJ6z1h4M3vd/ABjF5OssrkLmDM0MAI32v6+a1QOAMSbGGPOcMabWGNMHYDeARGNM4BM87gsAtgF42RjTZIz5tjEm3Fo7COABTM6kzcaYd40xc1wPEPz6eT749bMHkzNxyhW/0nLFz0MAoowxUz7BtV0dNzX48wxMzpYAAGvtAIBOAJmOx7gXwK0Aao0xu4wxq4LjOQD+OPiVtid43dnBxxZXIXOGphlApjHGXDGWHfz/P8bkTLXCWhsPYG1w/OPfHQQQc0Vcxsc/BGfhv7PWzsPkV9fbATwa1LZZa9dj8qtjGYB/vfqigv++/DNMfsVMstYmAui94rl/FjRh0lwfX0MsJmf/xqt/0Vp72Fp7J4A0AG8CeDUo1QP4hrU28Yr/xVhrf/wzvO5fWGTO0JQCuAzgd40xU4wxdwJYHtTiMPnvzB5jzDQAX7sq9gSAB40x4cHFmPs+Fowx64wxC4OzbB8mvwZOGGPSjTF3Bj/4owAGMPk192riAIwDaAcwxRjzNwDi/5/umfFjAF80xiwO/hvxmwAOWmtrrvwlY0xEMCebYK0dw+T9fXwP/wrgSWPMCjNJrDHmNmNM3M/42n8hkTlDYK29BOAeAL8JoAfAIwC2YNI43wMQDaADwAEAW68K/yqAfADdAP4Ok4s2H5MB4DVMfnDPA9iFya+6YQD+CJOzVBcm/w37lOPStgWf7wImv2qOIPh1+2eFtfYDTN7T65j8RpEP4EHy618AUBP8uv8kJv+9CmvtEQBPYHJxqhtABYILTuJ/ok0InxJjzEEAz1prf/BZX4v45UYz5/+CMeZ6Y0xG8GvtY5hML1w9Swrx/84nWb37VacIkwsasQCqANxnrW3+bC/p/4YxZiaAc0SeZ62t+3lejwiNvtYK4Sn6WiuEp4T8WltyfwGdVlf1xtK4gxEpzvGq87tpTHfrONXmZKZRbU04z/m/eLf7BNPgliznOABMRPNvrOH7O6gWhQGqzX2RZzl6/z7ZOZ4yVEVjYgOFVJu/5n+kHf+L8hTnfgYAQGofef1L+OvRGOCbjnpDvFZJ5jzV4rrd79mclfzzsTfyOqrt27SFauvm87kpejiHat1DNc7x7GNspyVQW1xCtb3vlTrz05o5hfAUmVMIT5E5hfAUmVMIT5E5hfAUmVMITwmZSgnM4SmH/T0FVIssu+gcj47iKYDwJXzfdldJPtW2RfP0wPwPpjvH7a/xlM5QVRPVzpa1US1xTjjVul8Zplr83UXO8bIPearqiZQEqu0Z4SmM9LhRqgUuZzvHwwdbaUx4XS/VoqcOUu3UYdcR0Ek2FI05x7un3EZj1kfyNFZDcwTVarv4hqipSe73BQByxoac4ztWP0ZjhhcdpxpDM6cQniJzCuEpMqcQniJzCuEpMqcQniJzCuEpIVMpc84lUa1nURTV5g67T2G81smX1+8L8JRIIMSJhHd2naFaZ777+QrOT6MxTUf4NcYP85MiJmIp1Vbez9Msda9WOseT1/GYNz/cT7X5qfykSFwEPxlx+dK/O8dbjy6hMW2BaKqtSeKptpz4dqrVHXenPlYv3ENjss/y1FLldQeodlMT/wynj/O4/q+4T0nNeY3PdR0v8ZQO/so9rJlTCE+ROYXwFJlTCE+ROYXwFJlTCE8JuVpbw/ebI29XA9XOFrlLD2WltDjHAeDEtA1U2/lyLdVy0ngPnJXV7o3Nh6fyioPT61nlSKCLLyijriORam+9xFeAr212dVsAYo4cpjHRG66hWlLSZarFxZdRrbnavam/Kcu9yRsA1mbxje9Tjxyi2rGUXKrVF7nrT90Uv57G9F13lGobfvQI1WbE8cMFyTdyLanZXfcptYivorde5iv9DM2cQniKzCmEp8icQniKzCmEp8icQniKzCmEp4RMpSyJiaFadC5fRs+/K905vv/YXBoTN87ruax/m9eIiV3O6/NMme1e8p51yn19AHBsKk8dDJbzl2tGfAXVWkd5HZ5tfe7HvCVnDY0ZaZ5PtX2dp6iWOMFrCF2Tea9zPK2uh8ZUNvL76lx+E9UqTm+j2gZSy6jIumsLAUBdxmKqdeWcpNq8gUiqtQd4qvBQmTs1FpnMN/Rn2I1UY2jmFMJTZE4hPEXmFMJTZE4hPEXmFMJTZE4hPCVkKuXAxQtUS5gbR7WwI+70RmoTXw6vreZtEGrm8lMYqWd51+g24z4FE4dOGtOTxZfQ4y7z1ME03gUB69JXUu3ErQuc4/05vPvz5VKeHliVUEO1yhp3ewoAuJjofm+iK/gplyOD/LUvWsavY2Whu3s1AKSMuFNSly6tozHLKnlKZM9SnqKr7uiiWvKE+7QQAGy8Idc5PtrvHgeApBCfHYZmTiE8ReYUwlNkTiE8ReYUwlNkTiE8ReYUwlNCplIGyzOottgcodqmD9xL1J8r5svaCJtNpQdOdVOtvJB3xL5j0N0uYN/2ZhrTlseXvL8Qz4uJDS3Ko9oI+MmO2eHu1E3Ccd5WYT9+SLVNvSEKcoVIIeWYa93Ptew5GtNYzk/3/FYtPznTPpWf3kiGuwVIxBd4K4m8ioNUS3mBX2N6Mm8nMRHDc2MbK922ebbuOzRmaw9PO30d/+wc18wphKfInEJ4iswphKfInEJ4iswphKfInEJ4SshUyvAa3lNkn82iWkysu/hXfxRPzawpTqVa3Vne2bpxGe8YPLjNnVaYv8TdFwQA9hzif68OPctPztw1fDN/zETey+NSu/u+qyJ42mNgMz8xMYUs8wPA6Cx+0qX5192njCLL+PsyZzH/fFwYr6ZaoJb3lQn7g373eOk/0ZiOTH5qafQe3r267G1epC4+kZ/gebXKfd+nwFNtGXv5KS6GZk4hPEXmFMJTZE4hPEXmFMJTZE4hPCXkam10Gl/NatvOO/8uTnRvKM4oLKAxGXfxkvqv1b5Ftb4yvkE5K365czxygq+EDuXzlcT8ntuoVtXMN+4H2vnqcOFS99/H4x9W0piW87ytRcZdfJN97jR+b6eHjHM8Gu7XEABuj3HHAIAJ4wcSAt/lq+9dR9yfkdrEIhqzsImv1s488BOq5S/ro1p3C7dGe5b7kMPc13lrkB3JvKs4QzOnEJ4icwrhKTKnEJ4icwrhKTKnEJ4icwrhKSFTKetbeOpjT94eqrVZd7n9lUt4h+ozL/LS+Nd+yV3fBgD2/eEzVBvKd29U35m0jMbcM8jrHM3p5S9XzZwqqs2atpBq3TvdNYRmVo/QmEWf46mZ5MIHqFbIOyugfdx9HanXhfOYCL4B30YtpVrH5t08rtm90f7RXH54oG0mv46Wseup1orNVJvYWsi1whed46dG+WcA4F3FGZo5hfAUmVMIT5E5hfAUmVMIT5E5hfAUmVMITwmZSukb4yc+0tq5Nhhwpw7qDvOaM13zd1BtYT1fDr/n5vup9sLml53jg7EJNGZeFk8f7Q7wVgeDF/lJl7sT+SmSAwsHneMRw7w20qIN/ARM/yHejbx6A5VQkOZuWxDo4Cc3Jk7+OtXClr5EtYbX5lCtKcmdZlm7Op/GVHXtpFpyBX+uzBJ+0iVxI3/9X2xzp+LKa3k6EFk8NcbQzCmEp8icQniKzCmEp8icQniKzCmEp8icQnhKyFRKeB5POdRV8U7D/ftqnOPF372HxoxG8VL2aOWnMFKW8jYI6+rcp1L6p/Ou3P1N/OTM+M38Gm/dxE9vnO3/kGrrO90pgrT7+Ws1tdt9ggQAWtav49rZN6h2otedSslPOkNjco4eopqdyz87B6L3Uu3mTHea4qE8fl+vjm2h2ukwXnQreoC3FBkdLKfa1Cz3yZm8PP6+dDVSiaKZUwhPkTmF8BSZUwhPkTmF8BSZUwhPkTmF8JSQqZTzVbwwVddWXuBr1mp3kaz4d3jRpJHrx6h2cvxWql2bzq+jZY57GX120uM0JiE2iWsVB6hWlOEuagYAt9V+QLU3It1vwZY3eLphfdEKqi1o4ymApPy1VFty/LBz/GBDCo0xv3mCaqcv8WpiT5XwnjNtEe6UyaZWfsqlf5CfPJmydSrVJop5obGda6dRbdEhd5fqw9U8fRQ2ZwnVaMynjhBC/FyQOYXwFJlTCE+ROYXwFJlTCE+ROYXwFGMtL7r11a88TsX6KF74qa/CXQgrNZBMY0Y3tlKtqIova58ruY5q01LcbdYLTvJW5HVH+bL8MnOJakN5+6nWXsWLhnVmXnSOt47F0pi4w7zQ2HVF/PTGywm8TfwKBJzjBb3HaUx/22qqxc7gabOyOPc9A0BTw/vO8cxpD9OY6hyeTotM4aenRlP5SaKk3bVUOzfe6xyvOcCLeN254jzVvvXNI843RjOnEJ4icwrhKTKnEJ4icwrhKTKnEJ4ScuN72SK+0mXOP021uOnuOjxnL/MN7HN56R68k8BbP1zbXU21pMBG53j/Gb6hv8XyDfhvjVdSreB8HtXaZ/DN13bK55zjxat+RGNOBO6jWkM7fyE3FvPuyhHb3avvByLcK94AkLHqFNXCxkuptrpnEdW2pvQ7x9Ni+EGA3qO8FtDUL/MV9pY2vsp7qTiHat2lzc7xkcwaGnPsGN8Uz9DMKYSnyJxCeIrMKYSnyJxCeIrMKYSnyJxCeErIVErA8rYF50/EUG1gdY9zfM0F3v0ZuRlUGjnJ69GkDPHUQfWpXc7xeze4l+sB4K6wh6jW/ArvdvyNE5updmMCry80K7zeOV6w/09pTGFONNUypvH0RvkJnhpL6HP3C5iaW0Zjquvd7zMAtLTyv/vZs3lvgt9Yf4dz/EgPT0WsHud1jj76w59Sbel1q6i2fT5/r+ta33GOD+yuojHLn1pMNYZmTiE8ReYUwlNkTiE8ReYUwlNkTiE8ReYUwlNCplJwpp1KSxfxLs/959ylh5LS+N+Cnj0TVEtI5KcOGgZXUm16nvv0wPvHW2jMfVH8pELM5+OpFl/H0yX9FwepdqbI3cHarOI1labxFXuU17prAQHAWBG/t+xs9wmelgneGfrG63lNqPf3Tqda5xme7jm1y52i64pyt0AAgJJ4/voWfIu3Qbjwaoh204383rKb3LWYIubNpjHL0ufz5yJo5hTCU2ROITxF5hTCU2ROITxF5hTCU2ROITwlZColt54vQ3+UcpBqD38u1TneFuBPF0jiBbKijnRTrXW6++QJAGSecF9H2Px5NGZzAy/Dv2D+I1T7yhO8PcWZLfzURGy8u0x/bDN/vIZE/lpdWNBFtc833k21k5HHnOMRbWk0pmk7Py1UXMOLmu1bOpdqa+e654vbevhJomMHTlMtc/x2qpXG8JROWic/jXMp7ivO8bKw52jM4DZ+kghfdA9r5hTCU2ROITxF5hTCU2ROITxF5hTCU2ROITwlZColqaSCakWV/PRDW4u7J8qUaUdpTHPEBf54R3j37Zkl6VSzQ+6TBYH2NhqTGMt7gww0bKda0YoiquX18jTAm2Xuv48FlvcTaW9zFwUDgPULeLfptlH+nl0z5r6OmnsW0JhVB4up9kzbD6iWHeCpoMun3Smkt5qjaExEAy82V3HuGapV1/CCc63x/D2bXeK+/t6uSBoTm/Dp50HNnEJ4iswphKfInEJ4iswphKfInEJ4SsjV2rH3+YpVd567IzMARMS6a7OEjfAO1cPtfJW0Pp2vqp2q5O0Ylma6uyHHTeFtBGJz1lCt/E2+2X/zMF+Jnrv0Tqrl5Lm7Zfevu5fGJO9xb+gHgME9vJPzuSG+0Xv+UveKeGTpWRrz6iW+8T03RKvyzEGutQ+4V9jTp43QmKoBXueo4UfjVCto5au85/6Cr7xu3efupv5Ha9fSmFmx/D1jaOYUwlNkTiE8ReYUwlNkTiE8ReYUwlNkTiE8JWQqJfbhx6iW3fw+1QYb3BvmY1N57Z7EGbz1w82dmVTrrGmgWlrkkHP80Gm+hF5VytMlax5aSrV3nn+Xanmlh6nWe9M1zvGh13gKICmlgD9eZCnVCpv45vHuaHedo+Xt/O93RaM7DQQAqWM3ca2QX0fPmDvl0Dh8hsakTOFdo/ufyKfa4nLejqHnOD/0UbLKXQ9oWyQ/hJET7U6/AACrcqSZUwhPkTmF8BSZUwhPkTmF8BSZUwhPkTmF8JSQqZT+Ab5UfuEAr0eT0DvNOR5x4w4aExHDT4M0nHZ3qAaAWVm8/UBygns5f20Er0n0QjtPESU28dYEkXflUu21fby0/33pNzjHG5fupzFJlbx2T28qb12Rmsy7ke85meAcP3yRvx6XEtZRrT/9PaqVVvL0xnjRPud46mWeTvuokb9WaaXXUm3rOP9c3R+XTbWRAreWnM07Ww9V8E7lDM2cQniKzCmEp8icQniKzCmEp8icQniKzCmEp4RMpUTV8WJRo6t5UaXRGnchr96O6TQmPZa3M0i/gxca27KbnzDJvdDiHC+4gacU5o/yVgENIToolyxaTjW7mHdyLid/Hme9y5fls1KbqNYazU+sRKSVU21ut/vETfitvIjXSD0/uRFoM/w66uOpVhlY4RzvTeeFxm7PfYhq7x7bSbUlAf6Za/trfm+Dpe6007u1B2jMXTvHqIbPu4c1cwrhKTKnEJ4icwrhKTKnEJ4icwrhKTKnEJ4SMpXSOYWfIskq5yc0Ht54s3P8cpy7MBIAnK10F+MCgPiwWqqdy15Ptab4Ouf4vCm8M/Sjs5ZR7T3wImRnKrj22Arepbqr57xzfPSCu8MzAGQ8wfvKLPv2Oaq1z3SnAACgd/nbzvH4i9fTmKaq41Qzt91DtbBjvEP43HZ3obGTBfzkyc4RnloqWHiMao1x/L0u2fw7VLsYt9k5fkvjNhozeFsO1RiaOYXwFJlTCE+ROYXwFJlTCE+ROYXwFGMtr6fze999korFMXyD9fL41c7x7iy+mXtu5Nep9tIg3/jeXMRry1z8jrsEfmzgAo2xdfwaG4bdtZEA4O6n5lOt6ShfaZw3bY5zvHCCr173x/Prj+nkHbEPNfCaOZHXz3KOR/NFaIw3f0C1pUl8JfRI9CGq5Ve6N5xPJPCaVQcOF1KtJX+Qap3pE1RDH3+++mvcq8M39o3SmOhefqDiq19+xnlKQDOnEJ4icwrhKTKnEJ4icwrhKTKnEJ4icwrhKSE3vg8ddy/zA0DbQp4eePOyWytOmElj6qL5UvPtC3i65ExgJ9XyV7o7KO+a4PV+xufzzdzXdPBUROwBd0dmAEhYxjsvLxhwb3A/ljBAY7LOurthA0DTm/x9qX6Q1wOaNeZOz8SMLaAx8wK8xlRdBq+Zk3OOfw52NLjTG6aMpz3Wroil2vkTvNN3xT5e56jwAZ6eWdHnPnhQerqExiS38E3x+LJ7WDOnEJ4icwrhKTKnEJ4icwrhKTKnEJ4icwrhKSFTKbXzTlItZhdfln8lzp2CeTQ9gsY8mMc7Mnc18WXt/o5VVBsvcR+pWP8RT1OUJSVTLW4PP6nw9PTXqPZI9Q1UK+93pxXiunlH6X8r5fWFBpbyrtFLtvA0ywPp9znHj3+Jn9KpPcXbU0TXvks1O5M/5k3Tw53jZ3kZKZzf4a4VBQA/DWug2rJeXtOqKmoh1Tp63O0kEme8SWOGZ/D3jKGZUwhPkTmF8BSZUwhPkTmF8BSZUwhPkTmF8JSQqZS4i51Um5nITx2sOf0t5/im7byFQ2dbEtVuSc2g2gA/CIDYrjjneGYhX9ZOjMqi2sWv8mX5gh3uNgIAUBrgJyMyK444xwOLeffnuZG8M/cbza9TLT+cdwF/fqO7XcD6Zl4gq6OXfz4q5vECX0su8lMk4Ufdp2NsifuEEQCcLeBdqA9E8bTZCPjpmHXv88J3u9d2OceHW3ibj44w3gmeoZlTCE+ROYXwFJlTCE+ROYXwFJlTCE+ROYXwlJCpFMTz3g/9I9zXPaSwVsdBnn75aHcZ1RZ8ny/nV73Hl6jXrO9xjgfaVtKYVWv4NUZ8yLXky+6+LACwrY93NY5a5S42NlrP0y8jqbzo1hML3T1PAKBsnJ+qaW9wpxwOj+/iMbH8uTZE87SNXcRTY+9uPuEcHw7naYqD9TzFtSaDp2AG1vJUyg+feZlqjxff4RxPf+BrNOb915+nGkMzpxCeInMK4SkypxCeInMK4SkypxCeInMK4SkhUylJl3hr+e35/ESC/cB9oiJ3Ci+elRCiANKJl56lWuo9vBV8fdeNzvGkxbyXem3V7VRLKt5JtV3D7hQAAKRmPEq1wtPu+z4+nfch6RzeQTVbydvVR1zmxb9uaHIXSisr5ic+enrHqTZ0nBfxasn9D6otmd3qHN87wXvR5Azwnj7n489T7cYsPjftfHyEakOZU53jsefcJ2oAYN2cBKoxNHMK4SkypxCeInMK4SkypxCeInMK4SnGWl4r5feffYCKZyveonGNe92bnmcUPUxjLvFuDJio/z7VxkpuptpfPXCTc7w2ga9AXn+Jb4buvobXCYqo4zWQ2stfolrtmHsV0hby2jfhY+6N1wAQXXeRalV7tlKt74j7rS56lLfQOP49fs/YwA8r5I7xuMpGdwfra9N4ewTE8Q7VhyP4qmv4KM8QFM8J8X4WlDvHe8i1A0DaRfchDAD486/tcN6AZk4hPEXmFMJTZE4hPEXmFMJTZE4hPEXmFMJTQm58rzjDl38bN/H6Qilx7hoxEQ/yGjaR1W/w63iHp0sC48ep9o/mp87x1Su/R2NeXuhujwAAGZfc9X4AYO0p3l05f5RvfO9LO+gcn9LBN44nLee1jHJn8TYIE6aDamdmv+AcP3KQX/uhbJ7/ujUpkmo92XwTeyHcaZHUL4aoWXUhkWoZI4epFpHDUx+JSTzFmDzm7r49mHcLjRkJ/4hqDM2cQniKzCmEp8icQniKzCmEp8icQniKzCmEp4Q8lSKE+OzQzCmEp8icQniKzCmEp8icQniKzCmEp8icQnjKfwLFqpPFw6cZkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAWc0lEQVR4nO2dW6hl2VWG/7Fu+3oudaqrq7uSvpBIgz4oAUXQB0NUFCR5UZB4exF9EVTwEtCAPkQD4lvyqESMaMxLJGICGkFQBJ8MQgJC0pq+pTtdfS61z76t2/ShTrAo5j9S3Va6h83/PVXtedbec821/r3OGf8cY1hKCUKIeBRv9QSEEHkkTiGCInEKERSJU4igSJxCBEXiFCIoEqcQQZE4A2NmT5rZpZmV38bP+G8z+xEy9l4ze+Hb9dnCR+IMxP1CSSk9l1JappSGt3Je4q1B4hTfFsyseqvn8P8difMhcPXE+y0z+w8zW5vZn5rZTTP7vJmtzOwLZnbt6mc/YGZfMrNzM/snM/vOq9c/CeBJAH979avsb5vZ02aWvnmjm9ktM/usmZ2a2VfM7JfumcPvm9mnzezPrz7zS2b2vQ94Ct9nZl82szMz+4SZTcl5JjP7jnv+/2dm9pGrf7/XzF4wsw+Z2csAPvFG1lL8LxLnw+MnAfwogGcAvB/A5wH8DoAbuLvOv2pmzwD4KwC/fvX653BXjE1K6ecBPAfg/Ve/yv5R5jM+BeAFALcA/BSAPzSz990z/oGrnzkG8FkAH3/Auf8sgB8D8O6r+X/4QU/6Ph4DcALgKQC//AbfQ1whcT48PpZSeiWl9CKAfwbwbymlf08p7QB8BsB7APw0gL9LKf1DSqkD8McAZgB+4Fu9uZk9AeAHAXwopbRLKX0RwJ8A+IV7fuxfUkqfu/ob9ZMAvucB5/7xlNLzKaVTAH8A4IMPeNz9jAB+L6W0Tylt3+B7iCskzofHK/f8e5v5/xJ3n3hf++aLKaURwPMA3vEA738LwGlKaXXPa1+779iX7/n3BsD0Af/2e/6+97z1AMfkePXqy0g8BCTON5eXcPdXPgCAmRmAJwC8ePWSl7/3EoATMzu457Un7zn2/8IT973nS+TnNgDm9/z/sfvGlX/4EJE431w+DeAnzOyHzawG8BsA9gD+9Wr8FQDvyh2YUnr+6uc+amZTM/tuAL8I4C8ewrx+xczeaWYnAH4XwF+Tn/sigJ8xs9LMfhzADz2EzxYEifNNJKX0nwB+DsDHANzG3cDR+1NK7dWPfBTAh68iub+ZeYsPAngad59sn8Hdv+++8BCm9pcA/h7AswC+CuAj5Od+7WrO57gbRPqbh/DZgmCqhCBETPTkFCIo2sXxNsfMngTwZTL8XSml597M+YgHR7/WChEU98n5/e97D1Xuru/pcdv9Pvv64HwRlE1Nx6r5hI6lgidsDGSO/a7Nvg4Awy4/97uDHR1aHM3p2HLJx+o6fwmapqHHTBd8bHGY3Xl39z2nfI3TmL82fTfSY7yxbu/s1R/5cU2dv55VafSY1Z1LPrba8M9q+H01OtPvh/xaVTW/LlXJpfaFT/1j9uT0N6cQQZE4hQiKxClEUCROIYIicQoRFIlTiKC4Vsr0kFsAxcBjzcU2r/neuJXSzHhY22puAexabulgzM8jOWF5q/n3VTXlc1weLejYYjHjn0e+H82ZI5xyX6NjV+133ApiVYpGYhvcPYjPsYCzxs6pVZZfj8KxX9oNt7/WF9xK6af83KqK2yJW5sdS4vdp277+/QR6cgoRFIlTiKBInEIEReIUIigSpxBBkTiFCIpvpSy5PVAP3MKoiB3Rg9svpWNT9DyKjtIbLPIxe+cIJCdkX1T8u4xldQBAu+cWBpuM14Fh6PhlS856eHNM5LB6wi0FllEDfItvfceCqUpipST+joVTYLAwPv+m5lZhM3HuR7LGuxWvBrpfO9lOBD05hQiKxClEUCROIYIicQoRFIlTiKC40dpmzuvRdC2PPllDorUjj/BaxTcNjzsnMuw0fU4gxw08ojk4G+mNbKQHgP2Kr8dY88hrXZEaQhN+XrOGb6RfLpZ8Hs6m+JZE34vCiZKSaDhAA+UAgBL8Wi8WeYeA1RYCgMG5B6oZX4/JzLm/nRpZ+4tV9vXN7oIfs339fZ305BQiKBKnEEGROIUIisQpRFAkTiGCInEKERTXSjl0at90Tqi/I5vHV1tez2Ucndj74ITsvePIBuXesUv2TjuGoeLnbN6YYwPUpEVCueRrP075Zu61s63fWSkYGfUSAcyxMBqnvcbU2VR+dO04f4zTSqJwLL9iwVs1uH6PU/Ln4Di/Yf6xx0/oMaw2koeenEIEReIUIigSpxBBkTiFCIrEKURQJE4hguJaKevLNR0bnRpCHWnVkJxskM7pkjx6Y71Ta6fP1+5xmmGjmfGQ/WzJQ/bLE579MD86oGM1sRUKr2O30yKhu+RWUHLWirVyLh1PYe5YbYeHfD1mM6d2zzS/xjPP1iN1hwBg49xz5vSFqJznVklqIA1OrajeaV/C0JNTiKBInEIEReIUIigSpxBBkTiFCIrEKURQXCvl9tk5HRuJTQFwy2Q2d9o7OOX2k5NP0Tv2gJGCVhPHLlkc8LHDR7klcnCDWweFk7HSbfPz36129JjdHae9Q8etDy8bhLVW6IeWHlOW/PaZL7ldMiN2CQBMSdGtuVO4rHPaO0wrPv9uz+3AzSUvyLU6vZN9/RsvvUyPWV/yjCyGnpxCBEXiFCIoEqcQQZE4hQiKxClEUCROIYLiWim9013ZSh6+ruu8HTGd81D+uOeflTpuHQxOhsnAOlHz+liYHfLBw0e5PVDN+HpsVjy7Z32WD9lfvspD7+0ltwAWXjfyA8feoFkffD0WS349qyn/3q/mfK1KMpaMW0RNwe2vxvhY52SRbMl1AYDbL55lX3/pq6/QY3Zr9UoR4m2DxClEUCROIYIicQoRFIlTiKC40dqD4yM6VjubuRfzfFRwPuN1YPYrHp28c5qPjgFAWfAo75LUsUklrytTLvmSTJyN453T4mHc8EjjuMvPpeBTxLUTHpG9fpO3BFg4m9FLUk9necA3+x8686gO+TpWpPM5AIyWj6Du9jzauW/5+nr1lvqO3zu7HY/kbjb5pITdjm+yH0envwNBT04hgiJxChEUiVOIoEicQgRF4hQiKBKnEEFxrRS35o/TLmBB6r14m7Krgk/Fa7lw6ZTbn5Du24sDXsOmnDu74p3y/Zd7bgUVLf8ONNKZe3mN18x55Ba3S05u5DtDA0DJEgHAmzzPnGSFouHXLHkWxtapgUTqAU1JMgUA9E737cEZ23VODaEdr+G03eSv9dDz9/Paa9BjXvcRQog3BYlTiKBInEIEReIUIigSpxBBkTiFCIprpVROG4TCGStJp+F6wsPh48jtjXrHs1lsxy2MdZsPhxcDn8dJza2UCjwcPin52MEBtyM2c5L1Meffm/Pr3JIq5/zcnO4JKEmWUe8kU3Rb3kUbXjdyJ4tkMctbSHPHWkpwrBSnVcNA2nUAwDBwu2cc82OjU3PrjaAnpxBBkTiFCIrEKURQJE4hgiJxChEUiVOIoPhWipNZMIx8Bz59PyeLoah58SmvIBcaPra+k+/MXU34d1JZ83OeOMkUywm3N5rD63RsIFkwtWO/jCVf+13HWz90Hbc+RtJh27VEvK/2Cb/W5TE/t9k8b5mMldPdfHRaeTh2SQE+Zomv8dDlC3ltt7wIWVKBLyHePkicQgRF4hQiKBKnEEGROIUIisQpRFBcK2XtdGQenWJRs13eFhmcQl3TGc8GWc4foWOPHPFeHml9Lfu67bmlUDjdjgunH0o58PccWv6e3Sa/jsWWZ+lMp9yKmAz8kvaOLZJ6YlU4dlqqPXvDyVpy7p1E7pFu5AW3+pGvfTc6hdcqfj8ulnz9j46I3fM4vz8ap0AZQ09OIYIicQoRFIlTiKBInEIEReIUIigSpxBBca2U6YIX1nLqJtHSX4NjU8ApkNUYH5tWTt+TirR03zn2gNN7xbUOjIfRvayJasiPjef8/Vrj9kDp9JxpGm5XNXU+q8a1PZyvdq9HyciTN7AjreC7KV/7uuLnPFs494dzz3kF7AqyxjefuMXn4dhffA5CiJBInEIEReIUIigSpxBBkTiFCIobrZ04m3+9btM1aceQnI3jPfjG5r0T5a2c9gnDOr8huljxz0ptvj4MAJQ1/y6bTng0rqr5Mhvplj06G8e7vbOOznVJpPYNAHTkuNTw86q9sZpH+gtnHStSu2fldcN26llVTqS/BX9Pr1UDC+SOTsS+9ewNgp6cQgRF4hQiKBKnEEGROIUIisQpRFAkTiGC4lop7Z6H3uGUsu/JpmfHHQAKZxNyw+uvVE44vyEblK3g9kt3yTeVn52f0TFb87D8jWtHdGxONkTXDV+Pcskv2+hs3N/u+I7zHbG5kjnWjDvmtLWY8NYbDbPhtit6zM65T+/seR2szZZbanfOTunY5UW+zUc94ddsOdfGdyHeNkicQgRF4hQiKBKnEEGROIUIisQpRFD8ztZONkVZ87BxTcLG7HUAaCY8i6Es+HcIs20AfnJ1ya2Zgzlv74CBe0Gr9SUdK0rnvEn2Rlnyz7KCj1VOzZ+DQ35uJanD0/d8fVvHttm33GbpHevDiDW2nPK57zqeDXLutBQpnNbczZTfj0fX8ut//ZF8+w8AuPEYbynC0JNTiKBInEIEReIUIigSpxBBkTiFCIrEKURQXCtlNuMFvqqGH1oTC6ZzimetnZD9duQZDlOnY3BPiiodVDwrZVbyc3784JiOnex4NsvY84yVssjPvyTZGQAwOu0dOsfCaJ1CaYlkGRVOywK3ZpVjf3kHbtf57JO9cw/sW97ZelLyaw0nS4cVXgOAgVT4OnuNZ84M3O2h6MkpRFAkTiGCInEKERSJU4igSJxCBEXiFCIorpWycDJFvBB7QcLQ7Y6HvC87ntUxOjbLYs7nOFbEpmicrsV0BKgGfs6L+ZKOXd65oGN7UlhrVvEiWPWE2wNLJ5NoTHwdO2L3jE7WT9tza2zn2Gad857rTf4eOVtzq+p0yzNPLnpexKuv+H0wczKoDo/z17qp+DFL5/5g6MkpRFAkTiGCInEKERSJU4igSJxCBMWN1iYnItv2fCdvt81H3LZrHlXzWjWUrJUwgNLZoNws8t8965ZvePZaDPROtPnmIa8fc3xywt+TrGPtbOgvnRpCydkUb85m9IrUhNrunI7jLd+M7kzD7RrdkbpErdOxe++MbfZOB/bGqVvltPk4vHE9+/qNRx+lx0waZwM+QU9OIYIicQoRFIlTiKBInEIEReIUIigSpxBBca2Ufcc3L2+cjch7shF58ML8lVOrZnBC704YvSefZ85G+oOG1xBa1nxs/xq3lp66+RgdOzo4zA84dgNLLACAwalXNDj2l5EaQrVTb6mu+Tp6NXhGp61FSSykZsZtj4XT+Tw565icTtT1zEkgmOTHzgaevGE7Pg+GnpxCBEXiFCIoEqcQQZE4hQiKxClEUCROIYLiWinjjlspw5ZndpSkVs1iyevidI690YHbJQXpyAzwTJey4d9JbeVYGEueqZBGHpZf9zyb5ajKWymN0327c9oqDE6rg97pRL1j9tfgfJbTY4DVJAKAzutNQO6douTXpXFsuFntpDvNuE00PeT3KshtcLl36mCBrz1DT04hgiJxChEUiVOIoEicQgRF4hQiKBKnEEFxrZS+48WdCuPh8IqEr70aR0XNvycKpw1C6WRNzOf5cPjMaWdQOXaDV1jr6ICX2585cxxB1tH4epROG4Gu55aDkwyCgRT/8iwulskCAJOJU6DMqeZWjHkLJg18PZJjtY2ls1aOzdKXTodwkq3l2Udv5DmoJ6cQQZE4hQiKxClEUCROIYIicQoRFIlTiKC4Vso7b/HeD4sZD5UvDvLb9veJZ7m8enFGxy7WPKuj3fFQf086KHt9PMzxGwanW3Nf8jlOjpwO4V3eSulIx2sAMOc71SviNTqF0soib1eZkx1jTi8dOFkYidlHACpyapWTlVKZY+k4/WGSk520de7Vntg9cHvYKCtFiLcNEqcQQZE4hQiKxClEUCROIYLiRmufeebddOz4gEcgJ7N81OqyXdFjFuc8KviNs3M6trrDo2rdlkRre37MmHgEcjLl52w1X8q24xuiZyn//egEIDF6gw6lE11lI6PTzsCNyDpzvLzk9adWbT7qvRp4Esa64OH3XeXcH4VTy6h22lCwyLbXUsRdxzx6cgoRFIlTiKBInEIEReIUIigSpxBBkTiFCIprpbxy+xt0bNfyLs/Hx3nLYTrn4embJ9fpWNPwEPXtyQUdW1/mLQwnKg9rHbuh4e0YJnNeQ6h1NqOfr4m95GxSb51u3sn5vp1NuBU0JR2lvXYXnXNevZMkMLb83PotuWaDU9Mn8XmcJd4ioW35uU1Hfj0Hsv77PU9+2G2dm46gJ6cQQZE4hQiKxClEUCROIYIicQoRFIlTiKC4VspXnv0vOrZY8iySk5NF9vXHb92gxxwfH9Cxg4XTZdjJfrDiTvb1VcHD2p1TO+b2Jc+quTjjFkw98u9A2+Y/r984NWw6bivMncyZx64/Qsfq4+Ps64U5dXGcNgj7gc+fFgoCgCJ/S+52PJPlYp/vyg0Al4nbG31ybn+vVcMmf//02w09xmmYTtGTU4igSJxCBEXiFCIoEqcQQZE4hQiKxClEUFwr5dWXT+nYacGLO738Yj5ufHbOM0ieetctOnb9kSM6tpjnbRsAqEjRrWbCC4adN9wu2TbcHthuuK2waXlYnhaFqrk1U/Q8O6YruMXVGr9mrJCX5wB0ib+f10V75dgi55u8LXLuHbPjFsbW8TDmNc+sspYXZTPSwXo2ca6Zu5LsGCFESCROIYIicQoRFIlTiKBInEIEReIUIiiulbInu+8BYL/lmQDDmM8E8IocFU7n32nDp3ntmNss02k+06VxeoZUTjGx8xk/55Z0qAaA5FgpRZf/vKrnxdAmPbcA6o5/3xbcHcCGWBWVk5XSEUsBANrOKXa1dzJMLvOZROdbp1CXYxHNl7xQ14lTlG3m3AfXb70j+3pyesecnb5Gxxh6cgoRFIlTiKBInEIEReIUIigSpxBBkTiFCIprpRRLb5jretzmd+BvtzyE/vLXX6VjyyNe4Gs2d3q2zPPHHU/5+02mh3RsOuXZLBcrnnGzNW4rtAOxnVqeaTFsub3R7Pg1KxIv/rUt8mNFwa9zO/JMnNWO22bnHR/bkGynccLn4d2l05Jng1RO8a/DBbdZjq/lba5uz62U1HLLj6EnpxBBkTiFCIrEKURQJE4hgiJxChEUN1rb7nlUbbd32gWQzr+VU3Pm9ml+wzMA2LMv0LHRicY9/a7851074l20FxMeyZ06NWcOJ7ydxB2ymRsAVlW+ZtE28WjtvufX5fSC10B64Yx3Krc+fyuUJa9J1DvdpksnWaFe8qjxzZvXsq+/Y8YTAUanLUThtFVopk4LjYY/ty53+Xt/u+LR39Ud7lQw9OQUIigSpxBBkTiFCIrEKURQJE4hgiJxChEU10q5uM03eu8dK6We5d92ICX/AWC94aHm/gVeq6Z3Oih3Y37snU9zm+JwwTcoLyd8U/zRnB93TGoZAUC7zJ93e8zD8vs1t1JeO+K1ap576et07NVTVh+JWxFzJ4Hg2nXexfza9RP+nqSLuTl1n/ZOLaOx4PdcUTg9Ixz6Pn9fTRo+j9mCXzOGnpxCBEXiFCIoEqcQQZE4hQiKxClEUCROIYJiybE3hBBvHXpyChEUiVOIoEicQgRF4hQiKBKnEEGROIUIyv8AWo5/p8SBtLIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAZrUlEQVR4nO2daYzlWVnGn3P3pW7tt/aqru6uruqeXphpZmMcNoGoQEIgfpAIoiESEkw0wUiMHzQIicbEBMIHjTHgFxQ0RIMhihkMIDMMMMwA3TO9d3XtXdu9VXX37fih7iTt5DyHHqHpIz6/pJPu89a5de7//3/uuf2+531fY62FECI8Ivd7AUIINxKnEIEicQoRKBKnEIEicQoRKBKnEIEicd4DjDGLxpi33uXPvt4Yc/kuf/ZNxpiVn2x14v8KEud9xlr7TWvtwk/jtYwxnzPGfOKn8Vri/iNx3keMMbH7vQYRLhLnveMRY8yLxpiCMeazxpjUy19LjTEfM8ZsAPjsK7+qGmPOG2OeN8YcGGP+0RjzhVfuhsaYjxpjNo0x68aY3+qOfQjArwP4A2NMyRjz5e74ojHm940xPzTG7HVfL3XHa73TGPOCMaZojHnaGHPuDtvHjDGr3bVcNsa8pTv+qDHme8aYfWPMbWPMX97TK/n/FWut/vyU/wBYBHABwDSAQQDfAvAJAG8C0ALw5wCSANLdsZXuvASAWwB+F0AcwHsANAB8omt/ef7Hu/a3A6gAGOjaP/fyz75iLd8BMNFdy0sAPty1PQRgE8BjAKIAPtD9+SSABQDLACa6PzsL4Hj3788AeH/37z0AHr/f1/zn8Y92znvHZ6y1y9baXQCfBPDe7ngHwB9ba+vW2uor5jwOIAbg09baprX2SzgU1p00AXy8a/8KgBIOheTj09bate5avgzgwe74hwD8tbX2WWtt21r7dwDq3XW0cSjSB4wxcWvtorX2+h1rmDPGDFtrS9bab9/9ZRF3i8R571i+4++3cLhzAcCWtbZG5kwAWLXdLcnxOgCwY61t3fHvCg53Lx8b5OePAPho9ytt0RhTxOFuP2GtvQbg9wD8CYBNY8w/GGNefg8fBDAP4JIx5rvGmHf+mN8v/hdInPeO6Tv+PgNgrft3XxrQOoBJY4whr/PjeLUpRssAPmmt7b/jT8Za+/cAYK39vLX2SRyK2OLw6zistVette8FMNId+ydjTPZV/m7xY5A47x0fMcZMGWMGAfwRgC/cxZxncPh18neMMTFjzLsAPPoqfudtAMdexc//DYAPG2MeM4dkjTHvMMbkjDELxphfNMYkAdQAVHH4lRzGmPcZY/LW2g6AYve1Oq/i94q7QOK8d3wewFcB3ABwHYcOIS/W2gYOnUAfxOFD/z4A/4rD/wfeDX+Lw/8jFo0x/3wXv+97AH4bwGcAFABcA/CbXXMSwJ8B2Mbh1+IRAH/Ytf0ygIvGmBKATwH4Ncf/n8VPiPmf/70RoWGMeRbAX1lrP3u/1yJ+tmjnDAxjzBuNMWPdr7UfAHAOwL/d73WJnz06oRIeCwC+CCCLw6/Ev2qtXb+/SxL3A32tFSJQ9LVWiEDxfq39i099hG6rndwgnXersOEcX9vdonPaBw1qGxjjIbTm4CS1Ie4e7kT5Z1L7gL/c9e97DsL0Gmo6feIItfXH3JfYNnvpnHaLvDEAQ3m+jvGj/J5FY+5HodNu0jmxeIra9nepCRsba9RW60Sd42999BSdE6m3qe1LX/0ytR2ZPU5tuXiG2paXX3ku5BCbHqJzBrL82n/8fX/qvGnaOYUIFIlTiECROIUIFIlTiECROIUIFIlTiEDxhlJaPUk+MT9AbQN19xnoVRJiAYD8aJraRo7zcMlajYcO2GdPJ8ZDEbXaHrWVOjwUNNXLk0EmR3LUlrLuMFFhj9+aRpSvIzM8TG31dovaUHWHdNpNlnoKJLK+x4ffz2adh82SCXcIabSvn86plHaorbzP7+fttU1qqyb4+pPWrYtYb57OqVd9z6kb7ZxCBIrEKUSgSJxCBIrEKUSgSJxCBIrEKUSgeEMpbY/7OpXkYZae3j7n+NBugs4Zn+In+uM5HoowjTKfF3P/Phvh6ziolqhtIMKzHwZ6uOu93vKEDqz7fdfKPK2j1OChlE6Lh1Iae7wU0e5G0TkeSfBrNTzDs4ViCZ4n3CjzkE4m7d4v0kke/rI1dyYLAFQq/Nr3Vfh7Gxji4bt4r/t5rHrktHNrm9oY2jmFCBSJU4hAkTiFCBSJU4hAkTiFCBSvt7YvP0Ft+0XuMUz3uL14uQFeF6d3nHtrd+v80HAqwuvYZBPuddQ6vHNAu8a9v73kUDYAJFp8jbsbvBh6J1Jxjjc9h7njhns7+6Pcs92f5etH0+25rBj+vqKeBIJOi7/nCHeuIhF3ryPqmZRJ8no/M9MPUNvRaV6XaHJyhNr2iSd6a/EWn1PlNoZ2TiECReIUIlAkTiECReIUIlAkTiECReIUIlC8oZR0jLuvUx7b6NiUc3yrvkjn2Linrs8eD31kIj3c1iFvz/JQRKXhDm0AAD9CDexsr1JbLssPo1dS7lBFfojP6cvxQ9mblrdP2GrxMFEn405kMA1+rep7/FolEzwE04nzK5khYbhMhIeIEiP8AP7pB8/xdXie4Uaar9FE3c9VLs3ldPSJ89TG0M4pRKBInEIEisQpRKBInEIEisQpRKBInEIEijeUsr9X4BM7vJ7OraVrzvEBT7fg6g4vm99u8q7ASU9WSplkznQyvmwKHm5IGP5Zlk7ydQzO8GyQnn53Nk42567DBAA2wtfRavIQRrvEwwNR637NyiavfbOytU5tZx45Q235MZ6BZEjpoZTn2UkleUgkPciv/UGb11SqgT/fwz3u9U9M85DOfmmf2hjaOYUIFIlTiECROIUIFIlTiECROIUIFIlTiEDxhlKKZR5WqER49sPiC990jh89MkfnDHiKT41muYs6WuefL+ska6La4q73SIOHG0ayo9Q2+RqeKTIwx0NBsah7/Z66Wli/xcNON1/ihaRmPNks58886hx/+uISnbO9zVtX9OR4uCQS5WGKWt0d3oh52l1EwENjuSwPcWUsD1dF27ydxFj/uHP8xYvP0TkXXrxAbe//Bfe4dk4hAkXiFCJQJE4hAkXiFCJQJE4hAkXiFCJQ/Fkpnr4hlY6nS7J198mYn8jTOakOL+DU8XSvThhPNkjKncmwuHubzilU3R2eAeD0Gd5348hDvK9MhVwPADARdwGt8hpfx6Wnv0dt23sr1Da08Hpqa8J9jXtGePhlPuIrvObp9J30dAifdMeQ1ts8O2YqOkZtuTQvAJfs8EwX0+LvrdN0h26uXeFF3urX+TPM0M4pRKBInEIEisQpRKBInEIEisQpRKB4vbWZHu7N2t3mh6+PTbrrx8wdX6BzRtP8UPnidXdNIgBYv8Ft+fyAczwL7uG1Y9zzN3mSrzES559zsRo/mB0hh/CvP3eJzinv8q7iT5zjHtmFx3hrgsUl9+8bSXKv68OPzFNbuZdfj0i/JxEg4/bWVhv8PZeTPGnCgK+f5BwAANoRLo29g5pzfHWTrzHTefX7oHZOIQJF4hQiUCROIQJF4hQiUCROIQJF4hQiULyhlMQgP4yeKfBWDXEknON9Ke5CT/Vyl/fxU/zA+cYSbwmwfNvt2p731KN54tzj1JYfm6W2docX/WlGeF2iyxdfco5vLC3SOZNHeS2mhcd4B+XcEA8hDVTdh/P7c+57CQDVUf58ROO8TlOSPB8AsH5twzl+av4onVP2tNCIR/ghexvhocIa64oOYGvLnVywtsPDesciPOzE0M4pRKBInEIEisQpRKBInEIEisQpRKBInEIEijeU0hNLUlvW42puNd2n9q2nxD0iPBSR9pTUnzt9mtq+8w13W4hvrV6lc849+SBfR9xd7wcAevf4Go9Yd3YMAMQx7R6ff5LOyZ7gtZiiWX7PDioH1DZ8xJ3ZEevj72uvye/naJp/7l95wR0uAYCbJDT29pM8NNOM8LYQHevp9B3h96XR3qU223Rfx07b/dwDQMO4W4P40M4pRKBInEIEisQpRKBInEIEisQpRKBInEIEijeUMhflhZNWKrzFQLvtzsJokq7FANBu8eyBSJK70SfnZ6ltetGdJbDo6ciMCX5JVlq8JcCDe7yIV+8Ov44DUXfX7l9687vonL6Jfmrbre5Q247h7QKqbfc16V/j1yNa5mGKWppfq5Thrzn/0FnneDPDM1nKVX4/D5o8cyaX4OtIR/mzmo26WzUkDX9Od0q84zhDO6cQgSJxChEoEqcQgSJxChEoEqcQgSJxChEo3lBKpeBxUZf2qc0SV3OxwOegzUMp+elxPi/NsyZOv+4h5/iDtWN0TizKMwua27zb9BFPhkYmzrskY8cd+ti+1qRTolHebXrM8k7OWfACazXSrTm96ukqHuUFsrbXeajtUU+vlALpbRL1ZDTFYjzMUih7eqxYHi6Z6edF4EDudTTB13Fy1J195EM7pxCBInEKESgSpxCBInEKESgSpxCBInEKESjeUIovE2ByaoLaqnW3+73d5AWyGjXu1i5ubFJbfnaG2tJDw87xsV3+vlpXblBbJs4zHJrgYadGmvdKmZh2hz6aDX6tWpeXqG2LX0Z0EjxzJpdxZ85kPf1tYkleTCwS4aGIXsPDZpmlPfK7eEgnNsTvZy7B12g90ZJOnF+rMomMzS+8hs6ZnznCfxlBO6cQgSJxChEoEqcQgSJxChEoEqcQgeL11sb7+cHm2Daf2t/rLnOfinEPWDzKbftrvPbN6Dj3GrdJSRdb5N7T5gY/nL/pOSgdz3i8kz3cY5hKuA9RZ8Y9h9TL3NtZL/Gy/7bFPcClovsAfsnTkiMa515SxHmdo8TwKLVND445xztlnpCws7lMbWOe56PIHznsV3niAevcnvPUuqpY7m1maOcUIlAkTiECReIUIlAkTiECReIUIlAkTiECxRtKKZW5+7fV4K7mBqkh1OjwEEaizWsBmQwP6VT3eegj0ec+VG76++icM0+8gdpeunqR2r713Wep7ez8A9Q22udux3CwybtQ9w3wjsyTR3nooFbmr7lTdNfaqVX4PUOMh3Ru7/L2A5k+HsM4suBux4ASD2Md7fAaTYu7nhpCvbwWU6nGn+/Fq5ed41evPEXnvH72PdTG0M4pRKBInEIEisQpRKBInEIEisQpRKBInEIEir+GUJVnOPRmeNZBFe5S/I0Ud4enej1tBLJ5amu1ucs7Rlo83Njj7vXjWR6KODV7htqee/YH1FZpcFd/OjPkHE/F+edmxPCO0utrG9SWTPIQxuQMaVFhPZ/fhrdImG7y+7K+sk5t1y684ByfO/s4nTOUmKO23Qv8d+338ZYRBfBrtbLvfn4GB07SOSeOn6I2hnZOIQJF4hQiUCROIQJF4hQiUCROIQJF4hQiUPwFvsBDAKbHXb4fAMaG3EWayh0emkl6akXtrKxRW2rYHYoAgOLainM87WlL8Nwib8fw2jl3p2wAePd7f4PaVm5cpbZ2zZ31kepxt5IA4P1IzXlCAG1PVtDm2m3neCLhzpoBgE6Th8ZiPfyGjs7yAl97m+5MqN11nuVS3HO3cACA8fEFaru5xdtaRHp4sa4jC0fdr/civ88bK+4CagAAEoHRzilEoEicQgSKxClEoEicQgSKxClEoEicQgSKN5SSTXM3epM1IgEwOOh2lY/VeaZCqcEzBDZXeXijjydGoNl0u6+z4+5QDwCspHnRKuMp8HX+1GupbXKU/76lS+7XTGZ5EbJ6w9Mpe8wTgjE8NFZtuLt2pzydoU2T9y+5TQqGAUA7yx+7dM7dc6a6z0MRzRp/dr6+dInaftTkBc+yPOkK/UPuYnRTC9N0Tn6UZzsxtHMKESgSpxCBInEKESgSpxCBInEKESheb22yj3sMo5YfsI5E3B6+9VvX6ZxilntJG55V3l7iB6JnZkec480q9zIOTXLP6vPPfJ/aspd4Z+uF41PUNjLjtiV6uLtwKMfrLTXLu9RmO7yzdZp55mP84hcOeCJDu+65aTW+J7RIfaS25R7qdIIfsl8ucK+xJQkaALC9zWsxVYrbzvFH3vA2Oic/LG+tED83SJxCBIrEKUSgSJxCBIrEKUSgSJxCBIo3lJLu4R2lCzV+4vzGZXfn31qBl8aPZdxhDwCokE7ZAFCtlqhtNu6u9XJpkdd6mdjlh6GHzp6gti899TVqe1vrPLU9fNLd4qHZ4KEl30dqPMu7XtfK/FqhScIsHU9biCQP6SQzfP1pzxtoWPe9rnvaO9QtDxGNH+NtEKIxHoIpRvj6p0bd96yd5EkCt2r82X+MyEw7pxCBInEKESgSpxCBInEKESgSpxCBInEKESj+UEqMZ1qseUrZL136oXP8NY+co3MiMV6TqNTmYZtY3yC11aruTIaRQR62ubLsDgMBwJPzs9Q2/tqz/DUX3W0hAOBYzf2ag708JAJPDSGAhwdSOd66olRwt2NAm7c6SHhCGK02b9VQ2PG05ehx1+exbf6ebYuH2uIpHhI52ONtPqaPznLbA8ed4zcLV+icgqfOEUM7pxCBInEKESgSpxCBInEKESgSpxCBInEKESjeUEqxyN3opT1eHr8/4w6LRDo8syCV5K7y/AAvJra4zd3yxYa7kNfc8Vk6py/PQzOXr/Iwy0NHePZDLsa7gDcMCTl4siLgKXYF65lHimcBQKbfHWbZ3uOFrqpbvJhYb45fx0yCryMC9zMykOEZHweGZxL1tfkax5I86yo3yrtvr9TdobFCi/8uWB6WZGjnFCJQJE4hAkXiFCJQJE4hAkXiFCJQJE4hAsUbSimXC9SWS/Lwxtve+g7n+PQ0702xWefdqxP7PDxQvMozPsoV9/p3m5N0znQPX+NOh/fduHDxArW98/Rj1DZM3PmNPX7tE20eioCn03ejwzN/YNzZLJEWf8Gsp/1zJsVDB9UD3qsmSYpkdcCzOsp9PISR3eTZMSfH56ntcoyHZ27vuYt1xdPujBoAaFc9IS6Cdk4hAkXiFCJQJE4hAkXiFCJQJE4hAsXrrR0c44eXZ0+c5rYBd5n+Dnitl8Ekr+sznS9TW6SHexO3WVmcjrszMQDcuLVMbfnMMLVF89w7uVwtUtsMaQkQrXu8e21+HX0ft23PgfNEzO19T1g+p9riiQzJBOmUDSDquWelkjvZotjg3tN6gT/GVU9ixFLrGrUVhvm9jjXc7zuT5Z3gG0nuNWZo5xQiUCROIQJF4hQiUCROIQJF4hQiUCROIQLFG0qpVPgB5cXSLWqrN92u8vwIP3CeBj9InwefF0nzt5Cd3neO1+rcrX3rgB84393jtXvOz7+Z2lIZvsbimjveMxzlh9QrntL+qwV+CNym+P08ls87x3MZfpjbxDyPT5R/7scMrwdUOnCHUlotfs9GoxPU9uIBD2M9t/EDahs5wpMVcgl3+KtZ5ffl5i0etgEpV6SdU4hAkTiFCBSJU4hAkTiFCBSJU4hAkTiFCBRvKGVng9fM2WjxrIkXLv3IOT534gSd8/DZ89SWBM+MGADPnIkQd/5Giocbpk7xbITbKzzMcvna16htYGCG2mrWHXY62OeZG0uLPIx1eZXbRgZ5SGo46e46PtjDr28uxzNPtrd51+jeNM/g6SftMMqlHJ2zdcBbg+zWS9S2vc/n5T2tK0rk2V+78TydM9DxFHciaOcUIlAkTiECReIUIlAkTiECReIUIlAkTiECxRtKKVV5AaehlLuIFwBcWXR3Q752k7dOONjnBZwef+IstSUMz5pIwR0WmYrwEMBWYpXa2lM8BFBJ8SJk2+WL1FYjbQt2cjxLpzw9QG17PfyWNkjxLABoseJfHV5MrLLHX29ogodtqp4QRqHoziKJxPl9Xi3ygm1PL1+htr4HT1Jb2vCsoKUr7iJwA54MnqzlGUEM7ZxCBIrEKUSgSJxCBIrEKUSgSJxCBIrEKUSgeEMpGV9xpxYvZpRqu1/2RxvX6Zziv/AwSx9vQYG5Mwt8Heh1jo+ChyKSaFGbTfAQAMZ51kHE0ybDkGJj+ykeWhoZ5gWthlrj1FbddRc8A4AD0nm5p85DRJUmfwZiFV7EKxN3F8gCgAIJ3dxcvUTnXPLYyhmezTI3eZzavv91nmX07ofdxdzmTz9E5zz/wlPUxtDOKUSgSJxCBIrEKUSgSJxCBIrEKUSgSJxCBIo3lBLN8oyEkqfleGsg4xxf6D9G5yy/dIPa/vM/vk1tmV6eKRKZcWdaxODuCwIAGUxR2xh4TCdheCioluKfgfWUOwSzs79O5+w2lqitvcl7g+QrPITULI84x4sJ970EgESKv16jwUNLhQOeRbJK2svf7uOPamGfPwNTCXc4DQC2b16ltqxn/TNz7jCRBX8GhnuGqI2hnVOIQJE4hQgUiVOIQJE4hQgUiVOIQPF6a22H14gp7Hi6+K66a6yce4y3Y2iV+VK2d3j7hK/8+zeo7S2/8jrneGOKt5Lw+dRS4J7LURyltgb4yfdduL2rsV7uUe4Br99kE/yg+rXnebuArW13h+3RKV6/qXhjkdoaFX5w31je6iA+5a77NL5wms7JzvAowH6NP8ORGF/H4PhbqK2eZAkEPGliZ5/rhaGdU4hAkTiFCBSJU4hAkTiFCBSJU4hAkTiFCBRvKGXvNq+Zc/G5H1JbpVxxjpsUD0UMTvPwRqPqfj0AWL7K2yc8g+86x2Pv8bQ6GOLvuc8TaMmA1+5htYwAYBju5IIaeJ2dfk8NpP5efktPvJHX03nq2f9yjj9f5nWftmvXqG1scI7apma4bXrKXR9pfIjXTdov8fDRLrgN5NoDQM5zUL0Kd3gmBl5za2SS2xjaOYUIFIlTiECROIUIFIlTiECROIUIFIlTiEAx1vJaKUKI+4d2TiECReIUIlAkTiECReIUIlAkTiECReIUIlD+G+iSy7Tyfs/LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAAD3CAYAAADmIkO7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAOdklEQVR4nO3dW2xcx33H8d/sLpfLu0jdKIqkZEuyZEm246R2JKd2nRpw3b7kpUAaGAgC1Cha9KUoiha9oMhD2regDdCHAkHgFnCDGm6CXtAgLVo0TlwkrlsnTmzHEWRKMiVKJEWueOdyV5w+iC5cdWbp84dl/oF+Py+SZnj+O7vcnw55zs5MiDEKgD+lnR4AgDTCCThFOAGnCCfgFOEEnCKcgFOE05EQwqMhhJ98AHU+F0J46YMYE3YO4XQkxvidGOPxD/MxQwifDyE892E+Jt4fwgk4RTh3QAjhYgjhd0MIb4YQ6iGEZ0MItRDC4yGEy1tfcySEMB9C+OjWv0dCCLMhhMe3/j0QQvhKCOFqCOFKCOELIYRy5vG+FEKYDCEshhD+K4Tw6Fb7U5J+T9KnQwjLIYTXitbGnUM4d87Tkn5O0hFJ90j6g/d2xhjflvQ7kp4LIXRLelbSX8YYv7X1JX8hqSXpqKQHJT0p6ZnMY70i6SOShiR9VdILIYRajPGbkv5Y0vMxxt4Y4wOG2rhDCOfO+bMY42SMcV7SH0n6zO1fEGP8sqTzkl6WdEDS70tSCGG/pF+Q9BsxxpUY44ykP5H0S6kHijE+F2OcizG2YoxflNQpKfm7bdHauHMqOz2A/8cm3/P3S5JGMl/3ZUl/L+lXYoyNrbZDkjokXQ0hvPt1pdtq/o8Qwm9J+uWtx4iS+iXtyTxeodq4cwjnzhl7z9/HJU3d/gUhhF5JfyrpK5I+H0L42taZdlJSQ9KeGGOr3YNs/X7525KekPRGjHEzhFCX9G7ybp+W9L5r487ix9qd8+shhNEQwpBu/bj6fOJrviTpP2OMz0j6R0l/LkkxxquS/lnSF0MI/SGE0tYFpJ9J1OjTrd8fZyVVQgh/qFtnzndNSzocQigZauMOIpw756u6FYIJSW9L+sJ7O0MIn5L0lKRf22r6TUkfDSE8vfXvz0qqSnpTUl3S3+jW76W3+ydJ35R0Trd+fF7X//4R9YWtP+dCCK8WrI07KDDZ+sMXQrgo6ZkY47/s9FjgF2dOwCnCCTjFj7WAU5w5Aafa3ud86hc/aTqt1heXk+0z09OWcto7lrtfLm1UuwrXazVtPy1MfPe1bN/YqYOmmvv37Uq2ry7bxjh2dMh0XK27lmwP/+c26Puzvpr/KO70VPH3waNnT5nG8bXnv5HtOziW+9xHe6uLa8n2rsG9pnovvvBvIdXOmRNwinACThFOwCnCCThFOAGnCCfgFOEEnCKcgFOEE3CKcAJOEU7AKcIJOEU4Aafazkqp7U7PmNjOUEc12R7Duqne3vH92b6ZenqGQDut9eLHSFIp3Mz2dfd1mmpWK73J9vm1a6Z6reaA6bjJty4n2/v32t4DvQP54zY38q9jzkBvetbMdpbqq9m+xpBtxs2BQ4eT7fW5FVO9HM6cgFOEE3CKcAJOEU7AKcIJONX+am3NdgVysLcv2R67N0z11hqb2b7uSvE1hNaaS6ZxdFQ7sn2L04umms3F5PIxqnbatsPsrKWv/m5nYDj93Erl9Pi2Ezcb2b5qV/51zOnrTr+ntnP6ofuzfcMHbYvYl/u7k+0raxdM9XI4cwJOEU7AKcIJOEU4AacIJ+AU4QScIpyAU4QTcIpwAk4RTsApwgk4RTgBpwgn4FTbWSmtjfzMgnYWl9M7W5dtS7aosZ4/sNwqvh7NTePzCiE/Q6OjZptFMnS4P9ne3W+bhbGwZntu2kw/t7mJSVO5u0+OZft6B3oK1xvcZds1+tDJu7N9G5stU83BgfT3pvdY/rEsOHMCThFOwCnCCThFOAGnCCfgFOEEnGp7K2V20bZo1YU330i2n3zglKneUHd+obHV+fxy+9lj6ulbPdsZHstvC3Hv40dMNTt60t+Ca+eum+pd+fGM6bhHPvHT6XGsvmOqN9A/lO3r6F4oXK/SXXwhN0ka3JW/JRVuFl9oTJL6M8/tH77+kqleDmdOwCnCCThFOAGnCCfgFOEEnCKcgFOEE3CKcAJOEU7AKcIJOEU4AacIJ+BU2w++q2LL7v7R8WT7fWceMtVbvbGS7XvrldcK1zt6cLdpHOMP5deI6RiwfYh6dir9AfeJVyZM9c4+9bDpuLmVqWT7mceOmeoNDufXCeqp5Hcqz5lZuGoaR60jP45Spf3bP+fq1Fyy/Z2JK6Z6OZw5AacIJ+AU4QScIpyAU4QTcIpwAk4RTsApwgk4RTgBpwgn4BThBJwinIBThBNwqu3H8rur+W0Q2lnPbDbdUbHN3Bg5dDDbd+knFwvXWyrZdn+uDfZm+/oWbDMc1mfTWxOcvv9eU73uoeK7RkvS0Hh6h+1qn20bhO5qNdv3wxeL75Z98sm7TOOQ8ruibwTb+3v5xoV0vca6qV4OZ07AKcIJOEU4AacIJ+AU4QScIpyAU22v/5dXmqaijdW1ZPvSjSVTva6e/OX84x8rfsuhft12e6C1ll9obNeS7TZRo5K+5XD05IOmeh29tv9vl5o3ku2dbxdfjEuSVluz2b5jIyOF68VqMI2jnnkvStJAr+171hnSt2cqHR/suY4zJ+AU4QScIpyAU4QTcIpwAk4RTsApwgk4RTgBpwgn4BThBJwinIBThBNwinACTrWdldLZa5u90TuQXmRqqZ5ezGo75cwsAEkaPjxauN74um2Br+pafobDRnXVVHNkby3Z3pi2bbO+etk2e2Ogty/Z3tdlWzCssZlfDK3SLD7bqXnV9vqWlVltTpI6bTNu1pvp4+5/+H5TvRzOnIBThBNwinACThFOwCnCCTjVfg2hLtsaK31D6St/5VLZVG9hdi7bt3v/3sL19lRsVyBn6vlxLK/ariaODg8l2/fsSl/F3U5HybYtxNx8Pdm+vtIy1VvZyF9hv/vI0cL1NlZsr+9QayPbt34zP8Z2lubTa2E1W7Yx5nDmBJwinIBThBNwinACThFOwCnCCThFOAGnCCfgFOEEnCKcgFOEE3CKcAJOEU7AqbZTGErRth7N4J70TIsO607CpfQsF0laWiy+W3ZntK0dc/bhx7J93375e6aa1fJAsr2v27h2T2PddNzo6MFke6Vs+541N/OzWRaX5wvX27V33DSO0Zn8rJQ3Z5dNNa8vpcc/0N9vqpfDmRNwinACThFOwCnCCThFOAGnCCfgFOEEnCKcgFOEE3CKcAJOEU7AKcIJOEU4Aafazkrp7bHNjNg7Npxs3zTOmHjrRz/O9g0M7ytcr9Vve17f+eFr2b5HHj5jqvn97/17sr1UrprqjRwo/npI0lJmr5eOim1WyvzcdLavq7+7cL3pK5dM43j94sVs39SmbW+T8XvS7+/yTdv+NjmcOQGnCCfgFOEEnCKcgFOEE3Cq/dXaPbtNRRfr6XV96nPXTPVa1fxaRmtLK4XrrVRtO2wvbzazfS+/nr+S287xE/ck22s123o0G03bFfHB/vRaRqsbtp2tQ+zM920Wf/2jbLtQLzbbjL+jy1Rz5Ua65skzR0z1cjhzAk4RTsApwgk4RTgBpwgn4BThBJwinIBThBNwinACThFOwCnCCThFOAGnCCfgVNtZKbF101T0By99N9n+wMfvM9VbWm6ze3Wp7VNImpmdNY3jxEdOZfvO/+iCqea4DiTb9w3a1jlqNWzHzU5fTbaXNvM7Q7dTbvPffmO5+No9zY2GaRx7DqR3WZekcrTNdDn18QeT7ZfqtvdADmdOwCnCCThFOAGnCCfgFOEEnCKcgFNt70Ms1eumog+cTi9addeYbcGwSk9+0ap3popflj9y8qRpHBfP5S+V3zW831Sz0pFevGyzuWaq12rabg/U+tILfM3NTpnq9dTabLmwWXwRslbZ9nrsWswvAFcdGDPVvLSY3hqi0WYBOAvOnIBThBNwinACThFOwCnCCThFOAGnCCfgFOEEnCKcgFOEE3CKcAJOEU7AqbYffB8ZtX2Y++Thfcn20fFhU73jjfQ6O5L0+jtvF643NTltGkdPm52Qu/bsMdUsZXZsXqsvmOrNLhs/ID6QHn9PZ81Ub2Rf/ns9NTVZuN7qiu1D5Wuz+bWHLsfi45CkntquZHsMH+y5jjMn4BThBJwinIBThBNwinACThFOwCnCCThFOAGnCCfgFOEEnCKcgFOEE3CKcAJOtZ2VMjGRXnZ+O99/9cVk+2M/+4ip3oOZ7R0k6fT4kcL1NsvFZ7JI0tVSfnuK63XbDIdSZ3on6nht0VTvjfMXTcf9/JmzyfZ9u9MzjLYzcTm/dUVPZ7VwvdJmp2kc0/X86zgbyqaalf2ZrSvOnzfVy+HMCThFOAGnCCfgFOEEnCKcgFOEE3Cq7a2UUkzvuryd6QvXk+1//ezfmeo1nn4823fixInC9e4buds0jnI5f3ugfmPJVHOhmd7lOVSK326QpIHqIdNxq7qZbF9atC00Nrg7v4v5Un2ucL2Ja5dN43j1Uv648dF7TTU1nb49c+yQbafsHM6cgFOEE3CKcAJOEU7AKcIJOEU4AacIJ+AU4QScIpyAU4QTcIpwAk4RTsApwgk41X5WSp9tVkp1oDvZPjebnq2ynb/9q3/N9n3mV4sv/HT2dHoxq+38VF9+a/lzV8+Zai6sbSTbb8zPmur11JdNxy3fSG8v39Fr+/+7p9yR7ZtZyC+UljO5YZsdM3zoQLZv8fI1U80nP51+/5SUfg2tOHMCThFOwCnCCThFOAGnCCfgVNurtTOTV0xFJy+kr4Ld8zHb+jaTr09l+77x9W8Xrlfr7zKN4+DgwWzfseHjppo3VueT7as9g6Z6rU7bVc3/eOnVZPv6vO21urnWzPYNH8y/jjmfeOIJ0ziuL+evDHf3pu8qbCfW0ue01YWGqV4OZ07AKcIJOEU4AacIJ+AU4QScIpyAU4QTcIpwAk4RTsApwgk4RTgBpwgn4BThBJwKMcadHgOABM6cgFOEE3CKcAJOEU7AKcIJOEU4Aaf+G6YGD1JinIJ1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "i = 4\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",
    "We are trying to detect data drift on high-dimensional (*32x32x3*) data using a multivariate MMD permutation test. 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: **UAE** (Untrained AutoEncoder), **BBSDs** (black-box shift detection using the classifier's softmax outputs) and **PCA**. \n",
    "\n",
    "#### Untrained AutoEncoder\n",
    "\n",
    "First we try UAE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:alibi_detect.utils.saving:Directory my_path does not exist and is now created.\n"
     ]
    }
   ],
   "source": [
    "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",
    "uae = UAE(encoder_net=encoder_net)\n",
    "preprocess_kwargs = {'model': uae, 'batch_size': 128}\n",
    "\n",
    "# initialise drift detector\n",
    "cd = MMDDrift(\n",
    "    p_val=.05,          # p-value for permutation test \n",
    "    X_ref=X_ref,        # reference data to test against\n",
    "    preprocess_kwargs=preprocess_kwargs,\n",
    "    kernel=gaussian_kernel,  # use the default Gaussian kernel in MMD\n",
    "    kernel_kwargs={'sigma': np.array([1.])},\n",
    "    chunk_size=1000,\n",
    "    n_permutations=5    # nb of permutations in the test, set to 5 for runtime\n",
    ")                       # purposes; should be much higher for a real test\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 optional `chunk_size` variable will be used to compute the maximum mean discrepancy distance between the 2 samples in chunks using [dask](https://dask.org/) to avoid potential out-of-memory errors. In terms of speed, the optimal `chunk_size` is application and hardware dependent, so it is often worth to test a few different values, including *None*. *None* means that the computation is done in-memory in NumPy.\n",
    "\n",
    "Let's check whether the detector thinks drift occurred within the original test set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? No!\n"
     ]
    }
   ],
   "source": [
    "preds_h0 = cd.predict(X_h0, return_p_val=True)\n",
    "labels = ['No!', 'Yes!']\n",
    "print('Drift? {}'.format(labels[preds_h0['data']['is_drift']]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, no drift occurred. The p-value of the permutation test is above the $0.05$ threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8\n"
     ]
    }
   ],
   "source": [
    "print(preds_h0['data']['p_val'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now check the predictions on the perturbed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corruption type: gaussian_noise\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: motion_blur\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: brightness\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: pixelate\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for x, c in zip(X_c, corruption):\n",
    "    preds = cd.predict(x, return_p_val=True)\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('')"
   ]
  },
  {
   "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",
    "for i in range(n_corr):\n",
    "    X_c[i] = scale_by_instance(X_c[i])\n",
    "X_ref = scale_by_instance(X_ref)\n",
    "X_h0 = scale_by_instance(X_h0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialisation of the drift 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": [
    "# use output softmax layer\n",
    "preprocess_kwargs = {'model': HiddenOutput(model=clf, layer=-1), 'batch_size': 128}\n",
    "\n",
    "cd = MMDDrift(\n",
    "    p_val=.05,\n",
    "    X_ref=X_ref,\n",
    "    preprocess_kwargs=preprocess_kwargs,\n",
    "    kernel_kwargs={'sigma': np.array([1.])},\n",
    "    chunk_size=1000,\n",
    "    n_permutations=5\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no drift on the original held out test set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? No!\n",
      "0.4\n"
     ]
    }
   ],
   "source": [
    "preds_h0 = cd.predict(X_h0)\n",
    "print('Drift? {}'.format(labels[preds_h0['data']['is_drift']]))\n",
    "print(preds_h0['data']['p_val'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare this with the perturbed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Corruption type: gaussian_noise\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: motion_blur\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: brightness\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: pixelate\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for x, c in zip(X_c, corruption):\n",
    "    preds = cd.predict(x)\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('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kernel bandwidth\n",
    "\n",
    "So far we have defined a specific bandwidth `sigma` for the Gaussian kernel. We can however also sum over a number of different kernel bandwidths or infer `sigma` from *X_ref* and *X* using the following heuristic: compute the pairwise distances between each of the instances in *X_ref* and *X*, and set `sigma` to the median distance.\n",
    "\n",
    "Let's first try a range of bandwidths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd = MMDDrift(\n",
    "    p_val=.05,\n",
    "    X_ref=X_ref,\n",
    "    preprocess_kwargs=preprocess_kwargs,\n",
    "    kernel_kwargs={'sigma': np.array([.5, 1., 5.])},\n",
    "    chunk_size=1000,\n",
    "    n_permutations=5\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original test set sample\n",
      "Drift? No!\n",
      "0.2\n",
      "\n",
      "Corruption type: gaussian_noise\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: motion_blur\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: brightness\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n",
      "Corruption type: pixelate\n",
      "Drift? Yes!\n",
      "Feature-wise p-values:\n",
      "0.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "preds_h0 = cd.predict(X_h0)\n",
    "print('Original test set sample')\n",
    "print('Drift? {}'.format(labels[preds_h0['data']['is_drift']]))\n",
    "print(preds_h0['data']['p_val'])\n",
    "print('')\n",
    "\n",
    "for x, c in zip(X_c, corruption):\n",
    "    preds = cd.predict(x)\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('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A bandwidth can also be inferred from *X_ref* and *X* using the heuristic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd = MMDDrift(\n",
    "    p_val=.05,\n",
    "    X_ref=X_ref,\n",
    "    preprocess_kwargs=preprocess_kwargs,\n",
    "    chunk_size=1000,\n",
    "    n_permutations=5\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? No!\n",
      "0.2\n"
     ]
    }
   ],
   "source": [
    "preds_h0 = cd.predict(X_h0)\n",
    "print('Drift? {}'.format(labels[preds_h0['data']['is_drift']]))\n",
    "print(preds_h0['data']['p_val'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inferred bandwidth: 1.4132\n"
     ]
    }
   ],
   "source": [
    "print('Inferred bandwidth: {:.4f}'.format(cd.permutation_test.keywords['sigma'].item()))"
   ]
  },
  {
   "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": 23,
   "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 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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd.X_ref = cd.preprocess_fn(X_ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift? Yes!\n",
      "0.0\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'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}