doc/source/examples/cd_online_wine.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Online Drift Detection on the Wine Quality Dataset\n",
    "\n",
    "In the context of deployed models, data (model queries) usually arrive sequentially and we wish to detect it as soon as possible after its occurence. One approach is to perform a test for drift every $W$ time-steps, using the $W$ samples that have arrived since the last test. Such a strategy could be implemented using any of the offline detectors implemented in `alibi-detect`, but being both sensitive to slight drift and responsive to severe drift is difficult. If the window size $W$ is too small then slight drift will be undetectable. If it is too large then the delay between test-points hampers responsiveness to severe drift.\n",
    "\n",
    "An alternative strategy is to perform a test each time data arrives. However the usual offline methods are not applicable because the process for computing p-values is too expensive and doesn't account for correlated test outcomes when using overlapping windows of test data. \n",
    "\n",
    "Online detectors instead work by computing the test-statistic once using the first $W$ data points and then updating the test-statistic sequentially at low cost. When no drift has occured the test-statistic fluctuates around its expected value and once drift occurs the test-statistic starts to drift upwards. When it exceeds some preconfigured threshold value, drift is detected.\n",
    "\n",
    "Unlike offline detectors which require the specification of a threshold p-value (a false positive rate), the online detectors in `alibi-detect` require the specification of an expected run-time (ERT) (an inverted FPR). This is the number of time-steps that we insist our detectors, on average, should run for in the absense of drift before making a false detection. Usually we would like the ERT to be large, however this results in insensitive detectors which are slow to respond when drift does occur. There is a tradeoff between the expected run time and the expected detection delay. \n",
    "\n",
    "To target the desired ERT, thresholds are configured during an initial configuration phase via simulation. This configuration process is only suitable when the amount reference data (most likely the training data of the model of interest) is relatively large (ideally around an order of magnitude larger than the desired ERT). Configuration can be expensive (less so with a GPU) but allows the detector to operate at low-cost during deployment. \n",
    "\n",
    "This notebook demonstrates online drift detection using two different two-sample distance metrics for the test-statistic, the maximum mean discrepency (MMD) and least-squared density difference (LSDD), both of which can be updated sequentially at low cost. \n",
    "\n",
    "### Backend\n",
    "\n",
    "The online detectors are implemented in both the *PyTorch* and *TensorFlow* frameworks with support for CPU and GPU. Various preprocessing steps are also supported out-of-the box in Alibi Detect for both frameworks and an example will be given in this notebook. Alibi Detect does however not install PyTorch for you. Check the [PyTorch docs](https://pytorch.org/) how to do this. \n",
    "\n",
    "### Dataset\n",
    "\n",
    "The [Wine Quality Data Set](https://archive.ics.uci.edu/ml/datasets/wine+quality) consists of 4898 and 1599 samples of white and red wine respectively. Each sample has an associated quality (as determined by experts) and 11 numeric features indicating its acidity, density, pH etc. We consider the regression problem of tring to predict the quality of white wine samples given these features. We will then consider whether the model remains suitable for predicting the quality of red wine samples or whether the associated change in the underlying distribution should be considered as drift."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Online detection with MMD and Pytorch\n",
    "\n",
    "The Maximum Mean Discepency (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",
    "$$\n",
    "MMD(F, p, q) = || \\mu_{p} - \\mu_{q} ||^2_{F}\n",
    "$$\n",
    "\n",
    "Given reference samples $\\{X_i\\}_{i=1}^{N}$ and test samples $\\{Y_i\\}_{i=t}^{t+W}$ we may compute an unbiased estimate $\\widehat{MMD}^2(F, \\{X_i\\}_{i=1}^N, \\{Y_i\\}_{i=t}^{t+W})$ of the squared MMD between the two underlying distributions. Depending on the size of the reference and test windows, $N$ and $W$ respectively, this can be relatively expensive. However, once computed it is possible to update the statistic to estimate to the squared MMD between the distributions underlying $\\{X_i\\}_{i=1}^{N}$ and $\\{Y_i\\}_{i=t+1}^{t+1+W}$ at a very low cost, making it suitable for online drift detection.\n",
    "\n",
    "By default we use 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import torch\n",
    "import tensorflow as tf\n",
    "import pandas as pd\n",
    "import scipy\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "np.random.seed(0)\n",
    "torch.manual_seed(0)\n",
    "tf.random.set_seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data\n",
    "\n",
    "First we load in the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fixed acidity</th>\n",
       "      <th>volatile acidity</th>\n",
       "      <th>citric acid</th>\n",
       "      <th>residual sugar</th>\n",
       "      <th>chlorides</th>\n",
       "      <th>free sulfur dioxide</th>\n",
       "      <th>total sulfur dioxide</th>\n",
       "      <th>density</th>\n",
       "      <th>pH</th>\n",
       "      <th>sulphates</th>\n",
       "      <th>alcohol</th>\n",
       "      <th>quality</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "      <td>4898.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>6.854788</td>\n",
       "      <td>0.278241</td>\n",
       "      <td>0.334192</td>\n",
       "      <td>6.391415</td>\n",
       "      <td>0.045772</td>\n",
       "      <td>35.308085</td>\n",
       "      <td>138.360657</td>\n",
       "      <td>0.994027</td>\n",
       "      <td>3.188267</td>\n",
       "      <td>0.489847</td>\n",
       "      <td>10.514267</td>\n",
       "      <td>5.877909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.843868</td>\n",
       "      <td>0.100795</td>\n",
       "      <td>0.121020</td>\n",
       "      <td>5.072058</td>\n",
       "      <td>0.021848</td>\n",
       "      <td>17.007137</td>\n",
       "      <td>42.498065</td>\n",
       "      <td>0.002991</td>\n",
       "      <td>0.151001</td>\n",
       "      <td>0.114126</td>\n",
       "      <td>1.230621</td>\n",
       "      <td>0.885639</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>3.800000</td>\n",
       "      <td>0.080000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.600000</td>\n",
       "      <td>0.009000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>9.000000</td>\n",
       "      <td>0.987110</td>\n",
       "      <td>2.720000</td>\n",
       "      <td>0.220000</td>\n",
       "      <td>8.000000</td>\n",
       "      <td>3.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>6.300000</td>\n",
       "      <td>0.210000</td>\n",
       "      <td>0.270000</td>\n",
       "      <td>1.700000</td>\n",
       "      <td>0.036000</td>\n",
       "      <td>23.000000</td>\n",
       "      <td>108.000000</td>\n",
       "      <td>0.991723</td>\n",
       "      <td>3.090000</td>\n",
       "      <td>0.410000</td>\n",
       "      <td>9.500000</td>\n",
       "      <td>5.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>6.800000</td>\n",
       "      <td>0.260000</td>\n",
       "      <td>0.320000</td>\n",
       "      <td>5.200000</td>\n",
       "      <td>0.043000</td>\n",
       "      <td>34.000000</td>\n",
       "      <td>134.000000</td>\n",
       "      <td>0.993740</td>\n",
       "      <td>3.180000</td>\n",
       "      <td>0.470000</td>\n",
       "      <td>10.400000</td>\n",
       "      <td>6.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>7.300000</td>\n",
       "      <td>0.320000</td>\n",
       "      <td>0.390000</td>\n",
       "      <td>9.900000</td>\n",
       "      <td>0.050000</td>\n",
       "      <td>46.000000</td>\n",
       "      <td>167.000000</td>\n",
       "      <td>0.996100</td>\n",
       "      <td>3.280000</td>\n",
       "      <td>0.550000</td>\n",
       "      <td>11.400000</td>\n",
       "      <td>6.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>14.200000</td>\n",
       "      <td>1.100000</td>\n",
       "      <td>1.660000</td>\n",
       "      <td>65.800000</td>\n",
       "      <td>0.346000</td>\n",
       "      <td>289.000000</td>\n",
       "      <td>440.000000</td>\n",
       "      <td>1.038980</td>\n",
       "      <td>3.820000</td>\n",
       "      <td>1.080000</td>\n",
       "      <td>14.200000</td>\n",
       "      <td>9.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       fixed acidity  volatile acidity  citric acid  residual sugar  \\\n",
       "count    4898.000000       4898.000000  4898.000000     4898.000000   \n",
       "mean        6.854788          0.278241     0.334192        6.391415   \n",
       "std         0.843868          0.100795     0.121020        5.072058   \n",
       "min         3.800000          0.080000     0.000000        0.600000   \n",
       "25%         6.300000          0.210000     0.270000        1.700000   \n",
       "50%         6.800000          0.260000     0.320000        5.200000   \n",
       "75%         7.300000          0.320000     0.390000        9.900000   \n",
       "max        14.200000          1.100000     1.660000       65.800000   \n",
       "\n",
       "         chlorides  free sulfur dioxide  total sulfur dioxide      density  \\\n",
       "count  4898.000000          4898.000000           4898.000000  4898.000000   \n",
       "mean      0.045772            35.308085            138.360657     0.994027   \n",
       "std       0.021848            17.007137             42.498065     0.002991   \n",
       "min       0.009000             2.000000              9.000000     0.987110   \n",
       "25%       0.036000            23.000000            108.000000     0.991723   \n",
       "50%       0.043000            34.000000            134.000000     0.993740   \n",
       "75%       0.050000            46.000000            167.000000     0.996100   \n",
       "max       0.346000           289.000000            440.000000     1.038980   \n",
       "\n",
       "                pH    sulphates      alcohol      quality  \n",
       "count  4898.000000  4898.000000  4898.000000  4898.000000  \n",
       "mean      3.188267     0.489847    10.514267     5.877909  \n",
       "std       0.151001     0.114126     1.230621     0.885639  \n",
       "min       2.720000     0.220000     8.000000     3.000000  \n",
       "25%       3.090000     0.410000     9.500000     5.000000  \n",
       "50%       3.180000     0.470000    10.400000     6.000000  \n",
       "75%       3.280000     0.550000    11.400000     6.000000  \n",
       "max       3.820000     1.080000    14.200000     9.000000  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "red = pd.read_csv(\n",
    "    \"https://storage.googleapis.com/seldon-datasets/wine_quality/winequality-red.csv\", sep=';'\n",
    ")\n",
    "white = pd.read_csv(\n",
    "    \"https://storage.googleapis.com/seldon-datasets/wine_quality/winequality-white.csv\", sep=';'\n",
    ")\n",
    "white.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the data for both red and white wine samples take the same format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fixed acidity</th>\n",
       "      <th>volatile acidity</th>\n",
       "      <th>citric acid</th>\n",
       "      <th>residual sugar</th>\n",
       "      <th>chlorides</th>\n",
       "      <th>free sulfur dioxide</th>\n",
       "      <th>total sulfur dioxide</th>\n",
       "      <th>density</th>\n",
       "      <th>pH</th>\n",
       "      <th>sulphates</th>\n",
       "      <th>alcohol</th>\n",
       "      <th>quality</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "      <td>1599.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>8.319637</td>\n",
       "      <td>0.527821</td>\n",
       "      <td>0.270976</td>\n",
       "      <td>2.538806</td>\n",
       "      <td>0.087467</td>\n",
       "      <td>15.874922</td>\n",
       "      <td>46.467792</td>\n",
       "      <td>0.996747</td>\n",
       "      <td>3.311113</td>\n",
       "      <td>0.658149</td>\n",
       "      <td>10.422983</td>\n",
       "      <td>5.636023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>1.741096</td>\n",
       "      <td>0.179060</td>\n",
       "      <td>0.194801</td>\n",
       "      <td>1.409928</td>\n",
       "      <td>0.047065</td>\n",
       "      <td>10.460157</td>\n",
       "      <td>32.895324</td>\n",
       "      <td>0.001887</td>\n",
       "      <td>0.154386</td>\n",
       "      <td>0.169507</td>\n",
       "      <td>1.065668</td>\n",
       "      <td>0.807569</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>4.600000</td>\n",
       "      <td>0.120000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.900000</td>\n",
       "      <td>0.012000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>0.990070</td>\n",
       "      <td>2.740000</td>\n",
       "      <td>0.330000</td>\n",
       "      <td>8.400000</td>\n",
       "      <td>3.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>7.100000</td>\n",
       "      <td>0.390000</td>\n",
       "      <td>0.090000</td>\n",
       "      <td>1.900000</td>\n",
       "      <td>0.070000</td>\n",
       "      <td>7.000000</td>\n",
       "      <td>22.000000</td>\n",
       "      <td>0.995600</td>\n",
       "      <td>3.210000</td>\n",
       "      <td>0.550000</td>\n",
       "      <td>9.500000</td>\n",
       "      <td>5.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>7.900000</td>\n",
       "      <td>0.520000</td>\n",
       "      <td>0.260000</td>\n",
       "      <td>2.200000</td>\n",
       "      <td>0.079000</td>\n",
       "      <td>14.000000</td>\n",
       "      <td>38.000000</td>\n",
       "      <td>0.996750</td>\n",
       "      <td>3.310000</td>\n",
       "      <td>0.620000</td>\n",
       "      <td>10.200000</td>\n",
       "      <td>6.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>9.200000</td>\n",
       "      <td>0.640000</td>\n",
       "      <td>0.420000</td>\n",
       "      <td>2.600000</td>\n",
       "      <td>0.090000</td>\n",
       "      <td>21.000000</td>\n",
       "      <td>62.000000</td>\n",
       "      <td>0.997835</td>\n",
       "      <td>3.400000</td>\n",
       "      <td>0.730000</td>\n",
       "      <td>11.100000</td>\n",
       "      <td>6.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>15.900000</td>\n",
       "      <td>1.580000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>15.500000</td>\n",
       "      <td>0.611000</td>\n",
       "      <td>72.000000</td>\n",
       "      <td>289.000000</td>\n",
       "      <td>1.003690</td>\n",
       "      <td>4.010000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>14.900000</td>\n",
       "      <td>8.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       fixed acidity  volatile acidity  citric acid  residual sugar  \\\n",
       "count    1599.000000       1599.000000  1599.000000     1599.000000   \n",
       "mean        8.319637          0.527821     0.270976        2.538806   \n",
       "std         1.741096          0.179060     0.194801        1.409928   \n",
       "min         4.600000          0.120000     0.000000        0.900000   \n",
       "25%         7.100000          0.390000     0.090000        1.900000   \n",
       "50%         7.900000          0.520000     0.260000        2.200000   \n",
       "75%         9.200000          0.640000     0.420000        2.600000   \n",
       "max        15.900000          1.580000     1.000000       15.500000   \n",
       "\n",
       "         chlorides  free sulfur dioxide  total sulfur dioxide      density  \\\n",
       "count  1599.000000          1599.000000           1599.000000  1599.000000   \n",
       "mean      0.087467            15.874922             46.467792     0.996747   \n",
       "std       0.047065            10.460157             32.895324     0.001887   \n",
       "min       0.012000             1.000000              6.000000     0.990070   \n",
       "25%       0.070000             7.000000             22.000000     0.995600   \n",
       "50%       0.079000            14.000000             38.000000     0.996750   \n",
       "75%       0.090000            21.000000             62.000000     0.997835   \n",
       "max       0.611000            72.000000            289.000000     1.003690   \n",
       "\n",
       "                pH    sulphates      alcohol      quality  \n",
       "count  1599.000000  1599.000000  1599.000000  1599.000000  \n",
       "mean      3.311113     0.658149    10.422983     5.636023  \n",
       "std       0.154386     0.169507     1.065668     0.807569  \n",
       "min       2.740000     0.330000     8.400000     3.000000  \n",
       "25%       3.210000     0.550000     9.500000     5.000000  \n",
       "50%       3.310000     0.620000    10.200000     6.000000  \n",
       "75%       3.400000     0.730000    11.100000     6.000000  \n",
       "max       4.010000     2.000000    14.900000     8.000000  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "red.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shuffle and normalise the data such that each feature takes a value in \\[0,1\\], as does the quality we seek to predict. We assue that our model was trained on white wine samples, which therefore forms the reference distribution, and that red wine samples can be considered to be drawn from a drifted distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "white, red = np.asarray(white, np.float32), np.asarray(red, np.float32)\n",
    "n_white, n_red = white.shape[0], red.shape[0]\n",
    "\n",
    "col_maxes = white.max(axis=0)\n",
    "white, red = white / col_maxes, red / col_maxes\n",
    "white, red = white[np.random.permutation(n_white)], red[np.random.permutation(n_red)]\n",
    "X = white[:, :-1]\n",
    "X_corr = red[:, :-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although it may not be necessary on this relatively low-dimensional data for which individual features are semantically meaningful, we demonstrate how [principle component analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis) can be performed as a preprocessing stage to project raw data onto a lower dimensional representation which more concisely captures the factors of variation in the data. As not to bias the detector it is necessary to fit the projection using a split of the data which isn't then passed as reference data. We additionally split off some white wine samples to act as undrifted data during deployment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = X[:(n_white//2)]\n",
    "X_ref = X[(n_white//2):(3*n_white//4)]\n",
    "X_h0 = X[(3*n_white//4):]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define a PCA object to be used as a preprocessing function to project the 11-D data onto a 2-D representation. We learn the first 2 principal components on the training split of the reference data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PCA(n_components=2)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca = PCA(2)\n",
    "pca.fit(X_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully the learned preprocessing step has learned a projection such that in the lower dimensional space the two samples are distinguishable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAACb8UlEQVR4nOz9eZRk2VnYi/72GeLEPOVcmVlz1tTVc/WgFt2NhAYQuJF5sJCRbbhgY56NWebdf8zlLuPFXbyFH/j52euxDMbWk8EywpaxaEBWS6InSd1d3VU91diVlTXlnBnzfMb9/tgRGVnVNXZldVdR51crVkSccUdU5Pft/Y1CSklISEhIyN2L9lEPICQkJCTkoyVUBCEhISF3OaEiCAkJCbnLCRVBSEhIyF1OqAhCQkJC7nKMj3oAH4TBwUG5devWj3oYISEhIXcUhw8fLkgphy7dfkcqgq1bt3Lo0KGPehghISEhdxRCiPOX2x6ahkJCQkLuckJFEBISEnKXEyqCkJCQkLucO9JHEBIScufiui5zc3N0Op2Peih/Y4lGo0xMTGCa5nUdHyqCkJCQD5W5uTlSqRRbt25FCPFRD+dvHFJKisUic3NzbNu27brOCRVBSEjIraFUgulpqFQgm4WpKcjn6XQ6oRK4hQghGBgYYHV19brPCX0EISEhG0+pBAcPgm3DwIB6PnhQbYdQCdxibvT7DRVBSEjIxjM9DYmEegjRfz09/VGPLOQyhIogJCRk4zl+HL79bfhP/wn+/M9hdhbicWUmugNIJpOX3f77v//7/NEf/REAX/7yl1lYWNjQ+z777LP89m//9oZe83oIfQQhISEby8wMvPKKEvwjI9BowDe/CU8/Ddu3f9Sjuyl+6Zd+ae31l7/8Zfbv38+mTZs27PrPPPMMzzzzzIZd73oJVwQhISEby0svwY4dEImA60IqBZalfARTUzd8uVK7xMG5gzx3+jkOzh2k1C7d1PB+53d+h3/37/4dAL/6q7/KJz/5SQCef/55vvjFL64d9+u//uvcf//9PP744ywvLwPwL//lv+R3f/d3+drXvsahQ4f44he/yAMPPEC73ebw4cM8/fTTPPzww3z2s59lcXHxovv6vs+2bduQUlKpVNB1nZdffhmAp556iunpab785S/zy7/8ywD83M/9HL/yK7/CE088wfbt2/naf/kvUK1Co8Hv/Kt/xSOPPMJ9993Hb/zGb9zU9wGhIggJCdloFhdh0yalDExTrQiyWUgmIZ+/oUv1lIDt2QzEBrA9+6aVwZNPPsl3v/tdAA4dOkSj0cB1Xb773e/y1FNPAdBsNnn88cd55513eOqpp/jDP/zDi67xkz/5kxw4cICvfOUrvP322xiGwT/9p/+Ur33taxw+fJif//mf59d//dcvOkfXdXbv3s3x48f53ve+x0MPPcR3v/tdbNtmdnaWqcsoycXFRb734ov85Z/8Cf/8X/wLMAy+9Z3vMH3iBK+/8gpvv/02hw8fXlMoH5TQNBQSErKxjI2pmWsup4Q/QLmsnMU3yHRxmoSZIBFR5/aep4vTPDbx2Aca3sMPP8zhw4ep1WpYlsVDDz3EoUOH+O53v7u2UohEIvzYj/3Y2vHf/va3r3rN9957j6NHj/LpT38aULP/sbGx9x335JNP8vLLL3P27Fl+7dd+jT/8wz/k6aef5pFHHrnsdT//+c+jOQ779u1juRsO+q0XXuBbL7zAgw8/DJpGo9Fgenp6TYl9EEJFEBISsrE8/TR85SvqdSajlEK5DJ/73A1fqtKpMBAbuGhb3IxTbBc/8PBM02Tbtm18+ctf5oknnuC+++7jhRde4PTp0+zdu3ftmF4Ipq7reJ531WtKKbnnnnt49dVXr3rcU089xb//9/+ehYUFfvM3f5Pf+Z3f4cUXX+TJJ5+87PGWZYHvg2EgpVy716/97/87/+hnf1Z9vxtAaBoKCQnZWHbsgC9+Ua0AZmfV8xe/qLbfINlolpbbumhby22RjWZvaohPPvkkv/u7v8tTTz3Fk08+ye///u/z4IMP3lD8fSqVol6vA7B7925WV1fXFIHruhw7dux95zz66KO88soraJpGNBrlgQce4A/+4A+uPpvXdQiCtbef/dSn+NIf/RGNdhuA+fl5VlZWrnvclyNcEYSEhGw8O3Z8IMF/KVMDUxycOwiolUDLbdF0m+wf2X9T133yySf5rd/6LT72sY+RSCSIRqNXnJVfiZ/7uZ/jl37pl4jFYrz66qt87Wtf41d+5VeoVqt4nsc/+2f/jHvuueeicyzLYnJykscff3xtHH/yJ3/Cvffee+UbRaPQbK69/cwnPsGJ48f52Kc+BUKQTCb5L//lvzA8PHxD41+P6C037iQOHDggw8Y0ISF3JidOnFgzwVwPpXaJ6eI0lU6FbDTL1MAU+diNOZ3veDwPOh1lJtJ1pRyMq8/jL/c9CyEOSykPXHpsuCIICQm5rcnH8h/YMfw3BsPoO95vAaGPICQkJOQuJ1QEISEhIXc5oSIICQkJucsJFUFISEjIXU6oCEJCQkLuckJFEBISEnId9ArO3SgfVWnpGyEMHw0JCblrkVIipUTTbt2c+KMqLX0jbMinF0J8SQixIoQ4eoX9Qgjx74QQp4UQ7wohHlq372eFENPdx89uxHhCQkL+BtFre/nccxe1u/ygnDt3jt27d/P3//7fZ//+/czOzvI7v/M7ly3r/Fu/9Vvs2rWLH/iBH+C9995737VuurT01762dq0rjeHDYKPU4JeBH77K/h8BprqPXwT+PYAQIg/8BvAY8CjwG0KI3AaNKSQk5E7nGr2PPyjT09P843/8jzl27Bjvvfce09PTvP766xeVdT58+DBf/epXefvtt/nGN77BG2+88b7rfODS0t/7Hn/5l3/JP//n/xyAb33rW5cdw4fFhpiGpJQvCyG2XuWQHwf+SKp6Fq8JIbJCiDHgB4FvSylLAEKIb6MUyp9sxLhCQkLucNb3Pob+8/Q0PPbBs423bNmyVu/nW9/6Ft/61rd48MEHAdbKOtfrdf723/7bxONxgCuadz5QaWlNU6Wluw1vrjSGmyktfSN8WM7icWB23fu57rYrbX8fQohfFEIcEkIcWu3W5Q4JCfkbTqWiWl6uZwN6HyfW9UaQUvJrv/ZrvP3227z99tucPn2aX/iFX7juaz311FN897vf5fXXX+dzn/sclUrl2qWl1917I8Zws9wxUUNSyv8gpTwgpTwwNDT0UQ8nJCTkwyCbhdbFZahptdT2DeKzn/0sX/rSl2g0GkC/rPNTTz3F17/+ddrtNvV6nb/4i7+47PkfqLT0dY7hw+LDihqaBybXvZ/obptHmYfWb3/xQxpTSEjI7c7UlPIJgFoJtFqqJPP+mytDvZ7PfOYznDhxgo997GMAa2WdH3roIX76p3+a+++/n+Hh4Suaej5QaenrHMPNlJa+ETasDHXXR/CXUsr3/Q8JIX4U+GXgcyjH8L+TUj7adRYfBnpRRG8CD/d8BlciLEMdEnLncqNlqCmVlE+gUlErgampG+59fDfyoZehFkL8CWpmPyiEmENFApkAUsrfB76BUgKngRbwv3X3lYQQ/xfQc8f/5rWUQEhIyF1GPn9TjuGQa7NRUUN/5xr7JfBPrrDvS8CXNmIcISEh1yCcXYdchjvGWRwSEnKT3KKY/A/CndgZ8U7iRr/fUBGEhNwtrI/JF6L/enr6Qx1GNBqlWCyGyuAWIaWkWCwSjUav+5yw1lBIyN1CpaJWAuuJx6FY/FCHMTExwdzcHGE+0K0jGo0yMTFx3ceHiiAk5G6hF5O/Lplqo2PyrwfTNNm2bduHes+QqxOahkJC7hamplQMfrMJUvZfX6YeTsjdRagIQkLuFnphmJalzEGWpd6HUUN3PaFpKCTkbiKMyQ+5DOGKICQkJOQuJ1QEISEhIXc5oSIICQkJucsJFUFISEjIXU6oCEJCQkLuckJFEBISEnKXEyqCkJCQkLucMI8gJCTk1hGWvb4jCFcEISEht4bbqOx1yNUJVwQhITdCOMO9ftaXvYb+8/R0mN18mxGuCEJCrpdwhntjVCqqzPV64nG1PeS2IlQEISHXy23S2OWOoVf2ej0fQdnrkGsTKoKQkOslnOHeGGHZ6zuGUBGEhFwv4Qz3xgjLXt8xbIizWAjxw8C/BXTgP0opf/uS/f8G+ET3bRwYllJmu/t84Eh33wUp5TMbMaaQkA1nakr5BECtBFotNcPdv/+jHdftTFj2+o7gphWBEEIHfg/4NDAHvCGEeFZKebx3jJTyV9cd/0+BB9ddoi2lfOBmxxEScsvpCbXpaTXDzWaVEvibNMMNo6LuSjZiRfAocFpKeQZACPFV4MeB41c4/u8Av7EB9w0JuXlKJTh0CE6cUO/37YOHH76y8PubPMPtRUUlEioqqtVS70Nzzt94NsJHMA7Mrns/1932PoQQW4BtwPPrNkeFEIeEEK8JIT5/pZsIIX6xe9yh1dXVDRh2yF1PqQTf+Q68/TbEYurx1lvw7W/fnSGhYVTUXcuH7Sz+AvA1KaW/btsWKeUB4GeA/48QYsflTpRS/gcp5QEp5YGhoaEPY6whf9OZnlYCP5dTNv94XL0ul+9O4RdGRd21bIQimAcm172f6G67HF8A/mT9BinlfPf5DPAiF/sPQkJuHZUKuC5Eo/1t0ajadjcKvzAq6q5lIxTBG8CUEGKbECKCEvbPXnqQEGIPkANeXbctJ4Swuq8HgY9zZd9CSMjGks2CaUKn09/W6ahtd6PwC+P+71puWhFIKT3gl4HngBPAf5NSHhNC/KYQYn0o6BeAr0op5bpte4FDQoh3gBeA314fbRQSckvpRcSUy2rm22qp17nc3Sn8wrj/uxZxsVy+Mzhw4IA8dOjQRz2MkL8J3GjUUEjIHYwQ4nDXJ3sRYfXRkLubfB4+8xn1CAm5SwlLTISEhITc5YSKICQkJOQuJ1QEISEhIXc5oSIICQkJucsJFUFISEjIXU6oCEJCQkLucsLw0ZCQu52w9PRdT7giCAm5m+mVnrZtVXrattX7u7H66l1MuCIICbmbWV96GvrP09M313chXGXcUYQrgpCQu5lbUXo6XGXccYSKICTkbuZWlJ4OG9zccYSmoZDrZ2YGXnoJFhdhbAyefhp2XLaPUMidwtSUmq2DWgm0Wqr09P79N36tnjnoO9+B8XHYvBkymf61i8WNG3fIhhKuCEKuj5kZ+MpXlJCYnFTPX/mK2n670TNNPPfc3WOS+KCfeaNKT683B42PQ70OR49Ctar2hw1ubmtCRRByfbz0kqrTn8uBpvVfv/TSRz2yi7kb7dM3+5l7yuCzn/3g/QfWm4M2b1aNbYSACxfCBjd3AKEiCLk+Fhf7y/wemYzafjtxp9mnN2L1cjt85vVO50xGmZZSKZif768y4O5bqd0hhIog5PoYG+sv83tUq2r7R8XlhOid1IB9o1Yvt8NnvtTpnMnA9u3wqU9drATuppXaHUSoCEKuj6efVm0cy2UIgv7rp5/+aMZzJSGqaXdOA/aNmslns7C0pGzyr76qnpeWPtzPfK1+x7fDqiXkioSKIOT62LEDvvhF9cc7O6uev/jFjy5q6EqCZb0Qut0bsG/UTH5wEN58E2o1SKfV85tvqu0fFtdyOt8Oq5aQKxKGj4ZcPzt2XF7w98IGZ2eVEEqllMPwVmaTVipqJbCeeBzabSWApqeVQMpmlb36dsxq7ZlTetm88MFWL4UCPPSQ+rw9ZbBtm9r+YSrqnjK4HBv1WUNuCaEiCLk5eiYa31eOQU1ToYOxmNr3QaNQrsXVBMvVBNLtxEbF8FcqMDp6sb9Gyg83bv9aJSU2Ml8hZMPZENOQEOKHhRDvCSFOCyH++WX2/5wQYlUI8Xb38Q/W7ftZIcR09/GzGzGekA+RnommWFR/4Kap7NMvvgjnz8Phw7fmvteySd8sH0YuwkbF8N+K7OAb4Xqc3hv1WUNuCTe9IhBC6MDvAZ8G5oA3hBDPSimPX3Lon0opf/mSc/PAbwAHAAkc7p5bvtlxhXxI9Ew09TroOpw5o1YHKyvguvDWWyp6ZL2JYiMKkvUEy60wAfUEWyKhPlurpd7fCsG1EauXj3q2fb2F6+6UldpdyEasCB4FTkspz0gpHeCrwI9f57mfBb4tpSx1hf+3gR/egDGFfFisn42+9prKNH73XRVZlEwqofDss/3Z4Y2ETF5tVn4rq1teT4TL7ZS9/FHPtkNH8B3PRiiCcWB23fu57rZL+b8JId4VQnxNCDF5g+cihPhFIcQhIcSh1dXVDRh2yPVQapc4OHeQ504/x8G5g5Talwi8qSllCioUVDhpq6VWBO22yjOYmlLve0L0esMIr6YwbnX28LUE2+2YvbwR2cEflI/aNBVy03xY4aN/AWyVUt6HmvX/5xu9gJTyP0gpD0gpDwwNDW34AEPeT08J2J7NQGwA27MvVga9Wfnp00rwDw2B56nnbFYJeMNQYYw9IXq9s8erKYwbjEm/pjK7lGsJtkvvXy4rn8j/8X/Al750e9ZfupXcan9NyC1nIxTBPDC57v1Ed9saUsqilNLuvv2PwMPXe27IR8d0cZqEmSARSSCEIBFJkDATTBenL54Vp1LKB7BzJ+zaBVu2wNatal+7rRRBT4he7+zxagrjBkwR11Rml+Nagm39/efm4H/9L7XqSaVu72J81+KjLlwX8pGxEeGjbwBTQohtKCH+BeBn1h8ghBiTUvaK0jwDnOi+fg74fwohct33nwF+bQPGFLIBVDoVBmIXx+rHzTjFdhHm182K02lwHCXwNU05ipeXldlk+3blRO4J0et1bF4tPLRWg0OH1OojlVLVUA3j/cqkVGLhe19nolzCGBjC3jZJIqfqJU0Xp3ls4gqOy2s5oteP7fBhNQbLUmPIdX/KL720cTH8H0a3r5t1kIeO4Duam1YEUkpPCPHLKKGuA1+SUh4TQvwmcEhK+SzwK0KIZwAPKAE/1z23JIT4v1DKBOA3pZRh8ZHbhGw0S8ttkYj0hXHLbZGNZi9O6JqcVGUNolEl3D/xCSW4JiaUmWi94LqWkD18GP7sz1T0URDApz8NDzzQVxjj43DsGBw5ogRvPA5nz8LeverYHjMz8Oyz6EvHiI6M42IQr1RpPbifeDatlNnVuJpgW6/MVleV8O901IoIVJ2d2XWur5sR5B9WBNOtalkZckewIQllUspvAN+4ZNu/WPf617jCTF9K+SXgSxsxjpCNZWpgioNzSuDFzTgtt0XTbbJ/ZD9kp/uz4l61yVOnlPAeGoInnriyoLqSkD18GP7Nv1ECb/duVdn0a19T5qVHH1X3OHRIOaanpvpmol7oau9+pZKKVDIMjE0TOK0m0Qst7M3jWGdnqe3fppTZB6U7/uqRwzTcCp2lZdg5xbAJKbi4GN/NCvIPS0BfKVM7bCZzVxBmFodckXwsz2MTjzFdnKbYLpKNZtk/sp98LP9+E49hKN/Ajc5U18+W/+RPlKDr1ciZmKBpQOm9wxz/kV1kW9Pc8/YhktlBdc/eca2WymruMT2tbPaDg4w6JjOustdHimVsPJrusFJmH4CZ0gwvnXuJmfIMTafJD/3E4zz4/HFaToeZ4jQ79CFSTQc+97n+WG5EkF+6epidVf6W9dwKAR2WgLirCRVByFXpKYP379iAhK71s2XDgPfeU7b2Tgc2baJpwJzWIL1QYiA2QMttcb52ns3xGKlLncXrqVSUkuh0SMWS7MjtYKmxRGdhDrHzCR6beEwpsxtkpjTDV979CrlYDkMYSCR/5Z8g+ukH2Xl0gcziIssjUVJf/Pm+f+BGZtqXWz3MziqT2/ryEVcS0Ddjgvqok9JCPlJCRRBy/VyuZ/EHNU/MzKhQy5UVJbSi0b6tvVSCIKCU1og1bYLxsbWoJX/3LoonZkhZSXVOp6ME34MP9q+dzaqs5jNnAEhFE6S0ARjNwA98Hj6AEgB46dxL5GI5ctEcs9VZBmID1J063zUWyP+dH0dKSbFdZOd6J/GVZtpCKMG7Xmj3VjJnzypzVyqlvufpaeWQv5qAvlET1OWUxuUUO7x/nGE00N84wjLUIdfHRvYs7l2rXFbO38VFeOcdJWQ6nTWbvzh7hlitSfVzn1g71X/wASojaSUwKxX1vGULPPxw//pTU8pnsH27qn20vKwijJ555qaE2GJjkYyloo7iZhzHd0hGkhTaBWCdI3093VDUammRI0tHOHjqBU4dfZn68uz7E9KOH1ffjeMov4vjqLFns9cOzbyR3IorJcTBxUlpcPslzoXcEsIVQcj1sb5nMfSfv/EN5ci9kRlj71pSqpm7pinB12qpiKN334XVVSJGlAv/20/AffvWTq0nTeKf/EFoZa58z/VmK9NUEUUbMJMdS45RtavkojlGk6OcLp+m7tQZiA7QdJp9R/p68nnK9+3ixKvPkmr6pPODdLwix/Qie/VtZHpCG1Ri3uioqtwKSsmdP6++p3374JFHrvwZbsQEdb1+izCS6K4hVAQhF3MlO/PioloJrEfX1Qzx/vv75ojvfEcJ+SC4smLoXcs01QxY15WwKxTg3nvhV38VDANL2qxs8kg4zYujlnY+dm3zzi2Ia39669N85d2vAJCxMgzGBjlTPsMDow9gGVbfkX4Jp2QB78DDOJEEDpB64VW8ZIQL1VnujXb7QPd8HkGgoqQ8D06eVN/Lpk392fiuXep7uvT/50acvderNMJIoruGUBGE9LmanVnT4JvfVDP4bBb27FGROqOjfeHjeXDunLrOgQNXtlP3+h/ncsqpatsqLDQa7bea1HUyjz3GYzEuH7V0qz7/JUqw1L1/pVMhG83yo7t+lHeW3mG2NstYcoyfuuen2JG/euLYpYl5fiZF3LYpB/X+Qa2W+i6Gh5Wgfe89JXRHRtT35HnKdPT1ryt/yNRUXzk89tiNOXuvV2mEkUR3DaEiCOnTMwV4nkraqtdVNM/KiiosVy4rBdFuw1//NUQi8LPrWkjMzioh4Th9O3XvulNTKk/g+HGVhHXhgjJ39LKSt22D++5TPoJ33lF5BG+8QT6b5bGpKZi4QeF/oxE0l1GCtb/6n0x7S8Q1i3R+kPKkSymh8/m9n78hZXRpYp69bRL99UOkU2k14+8J7aefVrkY27ap7OleBFU2qxL2CgWlLHVd/f/s39/3Azz22PVHcV2qNKan+5/92DE1jh07wkiiu4hQEYT0qVSU4D92TNmpMxkl9L/+dSVkNm+GEyfUbD6bVY/1poN6XSmHVKq/LR5Xq4Rz55S9O5tV1+mtHqpVdc+dO9XqIptVs+FCQcXPX7qquB4B380q7uUS4LrX7pZ2qT3c8yideItsPELw8INo7Q6DR8/A/u1XL09xGXaJQU4cepZI08fID1KeHKS1fyuPtHLvF9q5nBpLEKjH/v1KwcZi6rvI5fo+hNlZuOeevqnmes1h630oR47AK68owb9pk/r/+MpX+v2o75S2nyE3RagI7mQ2qgZNLyz09dfVjG/Pnv51hFCz0xdfVPtcV5krRkaUcHj5ZSVsp6aUQK9UlFmox9KSak5Tr/frEoEyAQWBmvEeONBvYrO8rISb511+VXGtEMl1WcW9XALOnFERRFdzcl5qD5+dpRHXSaHTEoIgHkMDcrMFLiTN6/9uSyVy757inswUF2IFarUC+XfL7P2hZ8iMX8ak1BPSvc9qGP3Vga73zTLRqBLaH9RU07vPsWNqJXZpEECvVlJYQ+iuIAwfvRMpleBb34Lf+z0lPA3jg4f2rQ8Lve8+df5rrykB3m6r2fnCgpo5FotKQJ85A9/9rhKepqnMPf/tvylhtWWLGo+Uyin85ptqlRCJKOF/9KgSPkKoe3qeGne9roSO5ynT0Hp6lUWvJ0SyF4ufzapjYjH16DlYr8SlVVHrdWJGjFa0P1cKYlGqK+eZr81ff0nr7pgz+THuHb2Xx3Z9gj1bHyY3V7j6eesrevZWB089pb7Ddls9DOPmyz0vLqqV33oyGbU95K4hVAR3Gj1b9pkzalbesxd73lVr8l+R9WGh6bRa+sfjSoDbtopcEULZ8dttJSxtW81GQQnXRx5R5/VWA72Y95UVeOghZeKJRPrCv9FQrz1PCWDT7CuwREIdI0R/jL1Z76Xlp6tV9T185zt9Jbguq3iNaFQpgqvNnC8tPW0YjPhRasMZ2m4bKSXl0jxH7XmGE8PXX9L6MiWzZ90S33zjq/zrv/x1nv3Pv87Cf///XV6J95TB3/t7SsHm82q11GsF2jPd3Iyppue4X8/6WkkhdwWhaehOoWcGeu01JWirVWXT7QnMS+3FV6JnBjpzRgm8EyeUQFlfLbPXg/jcObUiAKV0Wi0lUEEJt3Zbxb5HoyraxbbV/p4p4bnn1LXicbWq6K0yQJk7BgbUTHfTJhWB1GgoYR2LqePWO1L371fO5l756Z7JKplUSWm9FZGuK0XQzSomGlXCeH0p7MtxacmM7dtJ5nLsySe44BSoV5apF5eYfOwTDKaUkOw5f6/qM7gk8ma2OsdfH32WjIhy/zmNiuHxV+3X+JFYjIkr+THWj63d7kcNbYSt/umn1YoQ1EqgWlVBAb1aSSF3BaEiuBNYH9GiaeqxvKwUwuCgEozvvafMRNPT8OUvK5v4T/xEP+O2VIL/+l/VH73vq2uNjKgolTNnVMXQe+5RgrVS6beaHBxUq4FWS83qx8aUoHRdpSw2b1avT57sVwTtKa2TJ9Ux1ao6x3HU62QSdu+mFtPpfP1/4M80EAODJAf3koymYGiIeqXA3MnvU42BsWsvO9plcuWyUia2rRRYLzP54x/vr4hsu59VXCio70nXry+r+FJ7eKlEZnqaeysmbN/L81OzpMe2XnTKWn+GK3FJ5M27Z14hF0RIxFL4EZN4zMJzahyuHmdi8lNX9mNcr63+Rv1GO3Yox/BLL6nJwNiYUgIb1Ush5I4gVAR3Autt471wy82b+7P4mRklFE+cUILRdZXd/tVXVXjn/v1KCfzZnymhGImoWfbyspox1+tqdn7ihEpYqtWU0KlWVU+BI0eUYOnNxD1PzR4jEQCaTotqbRHn+9Ok/rRAstLGyuXUeE+dUoL/scdUNJFtww/8ALXRAWbfeoF0Nklm5jzBYoHm7BzBMz8JFcmF+RMkOj7xdIbGSo2F+f+BmZgg2euHvLqqVg6uqx5HjypFFgQbl1V8ifBNzHHl/gxXuUb5vl3MH36J9ulF3q2eYvjeJ0ieXsWPqu8vaSZZbi7ffLLWlfJArpSE1mPHjlDw3+UIKeVHPYYb5sCBA/LQoUMf9TA+PHomFiGUcO41gSmVlAIoFFTLxF5vgPl5JaSTSSUYo1FlwllZ6c/4pVSCR9eVwAwCtb8nRMbHlbBYXVVOydOn+6GemYyy/ZdKtJ02F7bkiJcbDLz8Bp2RATwNcss1zGpVxcRHIkp4/+APKjNQPM5pexnzzFmSZxfw8lkEwGoBvd3B3jpBY+cWrEiU2NvHMUoV7LiFZpqMDEwqf8bKivo+2m3Ys4eGdGm+c5iOpeE+/ACDP/ITZPc/fPXv9Qbptb1MmImLMp2vVs300nP++7H/TtWu8plChoQ08WMWNadGzIjx45Of6tcS+iD06gKtTwBbXFRK8eGHL84FCFtJ3pUIIQ5LKQ9cuj10Ft8JrI9o6TWBCQI1I960CX7mZ5SwHRxUwjoaVYI+nVa2+UZDCcx0N4HJstSqwnXVvmZTzdYzGaUomk1lLjp6VK0Wmk2lLLZtgyefhMcfX+sbUI3rGIkk6ZOnccaGENkMVsumk4gopVEq9TuVtdtqvPU65YE4uhXD2bYZb3wMd3wM5/57CByHoFrFyA9gzi2idWyCRAyr1cFvNfoO4p7PwrJon52m/vK30SoV5L37CZoNlv/w31A5enhD/xt6Jbktw6LYLmIZ1jVLWl/a9/ljmz+G4zu8mayjt9q0qgXq7RoPZ/ZdNQKop1CuGq10uV7OhULfFHitYnQhdy2haeij4kZsuVdrAjM9rWaBQ0NKaNu2+sO3bWWW6RV1MwwlAJaX1esgUIJUSqVEpqfVs2Eo00+rpcbTixgqFtW+aFQVhVtagtFRCnsHaXfKpColKnu3knEEyY6DEzPVamN5WY09FlPXSachCMhF4mjnLhDs3I7W6aCXK1AoEi2Uibg+pTNniVyYR0gwimW01VX0dAq0buXRXrKYYdB+81VIJXC3TCLTSUySuEDhf/3Z2qqg1C5dVCpiamDqA5WquGJ/hitwaXmJyfQkf2v33+Ll8y/zTizKVBF+NP4Am4a2XfE3sH5V0evLcHDu4PuV0OVKQhQK/QY+PcJ6QSGXECqCj4JLbbmnT/e7c+3Y0U/x73G1JjA9JXHfffDVryr7vuf1TUn79ikzSq+sca+Wja4rBWFZ6r71ugq5tO2+41XX4exZmqZBS7PxFk9jzZ8hMjxKMpfDrpTJvnqB1kNT+LEoA++cRu84aC2HmKFDq3utUkkpnPl5pZj+yT9hOG7SfO555NlzGFJD1usY80tYvgaNDqyuYp6aIbJaRsgAN2phbdqsavX3HOQTExCL0Tz3LsbwGM7mTWtfmZHN4Z4/p77u6xGkt6hB/OX6Pg/EBvjC/i9ct0JZv6qAq0QrXa4kRC+Kaj0fUb2gjVLGIRtPqAg+CtY7f+fmVNauZfUThNan+Pe4UtRIPq+cgSdPKpPNK6/0o3o+9SlVtuGNN5QS8DwlSINAHbtjBzz/PDgOrqHhAZ5lYJba6HYHc3aWjudQ1NqQjJBathG+TWtpFt03CVpNojGdoVfeQTiCWKWJLyS4DnorUONbnxBlGEoZPP886QMH0D7zY3h/+lX8wCVSbRFJpokMZSEIGDs+B8sFpOuBaRLVTMwz55RpaXxc9USu19UKY2IzzbhFpNUh+u5J9NUSstPCHRuGgwc5E6+SSF1FkN7CBvFX7ft8nVy6quhd69JopVIMzkzoeKfeItOGTZP7yDzzjFoZNpvXVS/oVgnr617VhHwkhIrgo2B9OYPDh5V9PpVS9vpLU/xLJRU7f+KE2r5vn3L8rRdQhYLa9tRTqmaPrqtw0gsXlPDfuVPdp2c22LpVmZKSSXjwQZw3XqPj2Zhtl0izjdax8QIfWm2a2Sippkdk2UVISaAJopUOzVQEEbPQNZ1NZ0sUBxOUR9LkFisI3cB0HTU2TevbpzVNrTqkpHnuNEWtjZWKYooEEUfDSmXUuLJZos89B24AQod0VimUXs2gTZvgp35q7eOn9mzD/Z3fJDE9jwZoxTIycEjkhuHdd0mcO0R28zbc8THsbZP4ucxFgrR65DBzzfPUbI9kI0U2mqXWXKT+winEY4/flDC8at/nLuuFryY0JBIp5Zogvtyq4tJopTVBm0oQf/zjrLotzrlNHtuUI5+7er2g3v1nq7PM1maZGphiNDHKUnOJV2ZfYTI9yWRm8qa+h+te1YR8JGyIIhBC/DDwbwEd+I9Syt++ZP//A/gHgAesAj8vpTzf3ecDR7qHXpBSPrMRY7qt6dlye1U+43EVnTMyovZnMio0tFRSWbPnzvWX8m+9pbZ/+tP9P+YLF/oZuysr/QigVksJ39deUyGmY2PKxn/hgpqZex4MD1PPxIkugSYDhJT40RhCBMi2TWKxiZeMgyYQXoBVbeNrIIVAs6IgfaSmkWk4NHI59Vk8CY6nxiZEv6aQrivnrtNmqdzAkBp6fpBWJkFjuMGwNUgsQI2t5/TslaewLPW91Lulm9e1T8xOTWGMTeG+dwFqdWQ6gTWxn5gXwDe/iXX/FuxalYjnk/n6N5HxGLXxIfSPP0JpvMTM6YPoQ8NkzAzFdpHXZl9lz8BuRjoGS93s4ZuZuV7Nr7B+pmxoBocXD4OEh8ceXstc3jWwi1PFU8CVVxXXFLRXiERaf/+G08DQDM6Uz+AHPuer5zE0g4bTWBvLB/0erndVE/LRcNOKQAihA78HfBqYA94QQjwrpTy+7rC3gANSypYQ4v8O/L+An+7ua0spH7jZcdwR9OzQs7NKoC8tKWFWraoVQSKhhLnr9nvVlkpqleD7/VLQxaLa9pnPqP1zc8rsYhhKobz4Yj8EdGZG5QF4nrpX7z6mqbZ/9rP4jovMpPANgyAeR0Yj+I5N5NwsUhMEnoeuaeApc48WgOGqBiouPu2YTrzpEF9YQaCjWzEl9F1XKSPHUQqp+9w5cwpj3w6iro+zZQKz2ca1fOqNArG2UH6KXE4dH4ko5RWJ9COeWi31/bmu+hznz5NcKsAnf0R9F/G4ukZ3FjxqajTPn8J0XPyBYRpbRgmaVXY9d4hTgJ4fJOEJgoig0qmSiqaoV1cZHNl9y2eu6wX40ZWj5KI5kDBXn2P/sBL0hVbhmquKqwnaq5l71t+/4TTIWlk6fofDC4eZzEwS1aNU7epNfw/Xs6oJ+ejYiBXBo8BpKeUZACHEV4EfB9YUgZTyhXXHvwb83Q24751FrzRyva4cum++qYT78LAS8LFYP3onk1HZnadPK2EnpTo/Gu2XYj54UNX16VXlPHKkX8Ihm1WCcH5eCf1aTY2hV764WFTXqVbh+HGsjoOoNxB+QJCKY3sOtaDBkCEwPB3DdpC6BlLQyzoxXR8CG+JRPEPiIdG8AMOMYNqO+myAC0oJAdLUEYaOMTtPKpHAGRvCSyVo3buHxGtv0nYdqHaUySoaVUrPstQNSyW1/aGHlFLRdWg0aB/8Pu3Zc4jVVTr33UMmliXeaKjjazVotYi/+Q5moUAnl8T3fDLlCsPGI8RiDtZzf4390z9J+/XXWSjXOVo7xbBI4LZdqjvHKc6+ChJqTu2ygvRmberrBXjdrpOJZkBC1Vb1f3rC/FrRSlcStEKIq9rm198/ZaXoeB2iZpRCu8DUwBQdt0PKSl00lmsxU5rhpXMvsdhYZCw5xtNbn94QX0nIrWMj8gjGgXWFapjrbrsSvwD8r3Xvo0KIQ0KI14QQn7/SSUKIX+wed2h1dfWmBvyh0yuNbNtKOPW6eAWBUgwf/7gSzOfPq209R3GvINuFC0ow9uzk8bh6fP3rynRUKCjBGI8rE9HwsGoFOTysrqnryjTTbqv727ba3m5Dp0Mkk0MSEOg6QbOJU14m0u5goqP5PoFpgoAAn04iigBM20NzPMymjRGNUh8bgFgM0+tm/gqhlACAJkCALyVNr4WtScpjGbRAYhQrSF2n9PBe9GxXgBqGipz6sR/rRx0ND8Mv/IL6HnbuhDNnsP/Tf8T/678mMj2DWaqROPQOhdIc9tKC+rzVqnoWAjORIOVrDKzWyXqCmC8gGmXgxHlWmiscHAdbl0w4UZa9Gt8eqrEccdGFztHVoyzUFzCEcVGhuZ5Zxfbs6y9Cdwk9AQ59Qdzx+sL3emfNUwNTNN0mTaeJlHKth7JAXJTHkIgkyDR9Fr79dXjuOcaPz+OsLgEqtLXttZVyiA5Q6VRoe20m05PXPZaZ0gxfefcrNN0mk+lJmm6Tr7z7Fcrt8g3nYIR8eNx0ZrEQ4ieBH5ZS/oPu+78HPCal/OXLHPt3gV8GnpZS2t1t41LKeSHEduB54IeklDNXu+dtk1l8vSGH/+N/wJ/+aV9g96pv9so9jI0pZdBsqoStnj235yN47jnlPxBCKY5MRgn9XnmFXk2dJ55QpqGZGTXrbzRUlJCU/YqerquUDajVQiQCo6PY0sOtVnG8DoHjYGoGQtOIlup4lgG6hhAaMdvD8X2E76+tDrx0nCCZwItGyBoJOHcOV9fxAhcCkIaGF9HRHR83Gce1dOZ2j2E98DDG/vsRZ88SVMts2vEAVGs0Tr6L7XWwrATpWA5tZYWKBZVtI7RXF1iwS+x/4RiJhoORTGKgQ6uFHYvQHsxgJTMMp0b6K6qxMbW6qFb7mdTbt9MazLMcdfmjfQ4nt2cYjA+i6zrThWkm0hOMpcawfZum22R7bju5aI79w/tpOk0sQ61UbM++aBbe23fp7P1KK4f1Nno3cC/yEZi6uZa5DFxz5XG5e7wx/wYDsQFEtzihXq4Sf/MIZcPlsakfpFpe4tT5N3EPPERkaJSlxhLTpWmy0SyVToWp/BSjydHryqIG+NKbX6LpNpWJq0u5UyZhJvj5h37+6n9PIbecK2UWb4RpaB5Y39V8orvt0gF8Cvh11ikBACnlfPf5jBDiReBB4KqK4LbgekMOZ2bgz/9cCZ+eEigWVVhno6GEuu8rhTI5eXFmaT6vQkCXlpRQi0RUC0foF3ObnOzXpn/55X6ryL174XvfU6/jcTXG1VWlBKTsVy31fbhwAWvLFqz8EK3zZxG2i9Q97GwaJ2kRadlgB2BaXbOMwNUtzHabAIlv6EggMbsM+QFcw8DHRwrQBTgE+IHEDAQyEaW0bYRBM8WZlEYkY5ANXGqPP8Cs5qGdPcLgzjFGj5wnaJSY2Qz2viGyZpJzE3FOX3iJ+05VSDU9mhHwvSZxGUFLJQhiqufBmw+M8vFdP0Tq290+De22Mov18idcl44mqa5ewNm1lfxqBW1HjoX6ApPZSdKxNGkrzUprhbgRZ+/gXpJm8n3mGuD6wjqvETrZs/+3O20eHH0QicSTHkkjuWY6uZ7Qy8uZjy41GVlnZ2lFNOKpIRCCTH6MXTzE7IUV5pMmQ4khntj8xJqSutF+0YuNRSbTkzScBsvNZZpOk5gRY1kuX/W8kI+WjVAEbwBTQohtKAXwBeBn1h8ghHgQ+APUymFl3fYc0JJS2kKIQeDjKEfy7c+lrQ3Xd9JaH6Hx0ksqoUcIJYh6TVN6UUKplBLOzeblK2Tm8/B3/65SMr6vzErf/rZahezdq0wlY2PK1PPOO+peDz2klIauq+v3Moi1dZZAXe9n6Pbi+2MxhBegSQi8AKtax49GcSM+eJKIbYOmIQwDrdMiECANE73doaULolHVY9eNR9BaHaQRwSVA73SItkECerGCSBm085LCQIzHzCFqfoAZT9GuL2Al4yzbZdL5BJaRZ3bQp1WfI5sa4rmFQ4wNpLEiDpIqcRs84eMZEj2ZR7dd7HSCxc//EPMLGnu2bVPlKDxPfS++r1ZPlkU1ZdHcu4eSW2M0iLMttw0AUzeJ6TFs3+axcfX/6PjOFc01V7PL92bmNbt2UUSPJz3OV85zqniKxydUeOrV7P89JfBBQi8vtc17xVWaKZP9mUmqnSoXqrM07BrZZsAj449cU7Fci7HkGAv1BYotZf5JRpIU20UEglK7FJqCblNuWhFIKT0hxC8Dz6HCR78kpTwmhPhN4JCU8lngd4Ak8N+7S9RemOhe4A+EEAHKX/Hbl0Qb3b5c2toQLp+6v7iobNpnz6rSyNWqmqUKoVpCCqH233fflStA9pLGvvpVJdgyGaU8ymWVY9AL0dymhBknT/bNQYbRL9LWqx4aj/dDOqWERAK300G2Ghi9sE9dIH1B4DoIz8NQvl/cIEC6DroEX9eQQmK4AYmWhxeLEvElrYkRYoUq0rHRV0sYEnwBnZiBbjtMTC9z7hOb0E2L19/7a8amNpP3BC2vhTUxSe71d9GOnsPyNVKJFq18hNOPpSkvLuCm8pzbmmPsbJF008ORAX7gEV1apRMRnNi2heym7SxkffbUTJVg12zSiUaw3SZ+PkHzEx9nYVOKYT2Ns7hMamgc27OJGBEadoNNqU28tfQWA/EBYkbsInPNYmOR6eI0k+lJ0tE0s9VZXN/F9V1M3cTUTZKRJFE9etHs/WOTHwOUE/joylGiehRNaFcMy1xv5jlZOLkWQbT2U7tOx+2leQzjuTz7Y8MAHF05SsyIkQssWongusNDr+Ygf3rr0/z2936buBknqSep23Vs3+bpLU+HOQO3MRuSRyCl/AbwjUu2/Yt1rz91hfNeAe7diDF86FyursvlUvfHxpRA3rFD2fH371fJXpOTKrt3cFDNznt9A3pc6n/oNXJJJJRymZtTgrxU6veyHRhQimZkRNUD6tUlcrtu21isr4R65iEpcXsF6Ax97fZCghQakaaNCCToAj9qodsOSKlm977S30EihjQtRKcD0ThePkPZt8mcq6hSRVGTdkRHCvBTEYTtEj0xjfOTT1JKCsaOnyR37jX2xCXNTJyR4+fRaw2ausCQbYZkHOoBgzLK6+k2wzmTXZuyxOfr+I6P5vogfPx8hupTj7JQmuYHCjF4513YvZtOqUi1toIwLM49uodCuk3zwnmC8T042yZpp2PsyG/hfPU8EknMjPHje34cS7c4unIUz/dIRpJcqF24yG6+1FhiujjNUGKIiK5KSi81lrh/9P6LZu+D8UGmi9M8Mv4Is9VZYkYMJKSt9GVn95eakizd4vDiYQ5sOkDGUm0lbyT0Mt+Gx+aBChCdgFKZY9U5YpEoCRdEu4P24H4SpnFNYX144TB//M4f03AbjCXGuGf4Hkrt0poC2ZHfwROTT3CmdIblxjKD8UE+vuXjTKQmwpyB25gws/iDcrm6LpdL3e91gMrl1Iw9n1cC+8ABZbbpOZlBXW92VlUMXVxUheWmppQZ6S/+QvkFhFD+hfvuU9VBZ2dVNu6+fSq6yLKUmadS6begBHW9nsO43VZmIimV76Ib1SN8X73WdRAC3QswA4lrGEhdQ5gmOA5SE2i+VLoE8JMpRLWG8APKozlWO0WSrTq+rmEJcDRBBLCRtKSLO5ggQZRG2iL3317hWMJkddMIQ2cLbP/2NJ1cGmd8lKJbJdL0SbR9Mt8/zvwntvG94BSzfo35+7Yh9TliSy4xPU1nfBNBOsHWY/OYlRrjBxdhZCfN0UFOxaq00lHqhslgoYD28JM0tg7x3dVzPGrcw/JYCjdwGU2OsiO/A13oa0lcD44+uBbueHjxMFMDU4wlVYeyYrvIeHqctJVmMj3JbG2Wcwvn6LgdPrH9E2tCeyo/xStzr9B0mtQ6NSzDouN12JFXK8BLZ/eXJoftGtjFoYVDnCqc4sCmAzcWenk5X5YQtBolMmYKP5PG3r1DZVtLeVVhPVOa4Q8O/QExM7bmB3jx3Iv84NYfvEiB7Bvax47cjvc50cOcgduXUBF8UK5WCG49l+sA9Q//4cVmoN4fq+8rIX7+vBL65bKKOOoJcNtWDWlWVlTBtfvvV9t9XymAXE69X15eaxqDlOr66bSa9eu6Mgv5/lpzd0+TaLqO5thIBMLzVJRJ0I0L8jyk1JC+11UCAYEAN25hCgPRqOMS0MhH8RtlBjsmyWIT2g7SB8N2ILAQOhiBpBn1qY0lyT3/KpG2j1Mr0W7YrNodJoWk7Xcgk4RSnWHHpD2coxITRBJpfqqyhVlrgaXyCucODJGP7WFvfAuR96ZpRQSGZjJVMQjOznA8YVBuOdhBB2HqQEBZs2kJj0knhpcdY27vOIPDm+isLDC61GT89AybJvcx3Zl5X30iP/CpLZxjR6OIXq0Ts89h7djOeWeJaqdKzIixKbWJ1dYqR1eOsn94Pxkrg6mba6GTAQGBDNS+6OVn95cmh2WiGR4ee5ijq0dvyHELXN6XNTpKrNJg6b6dN5Tg9dK5l9A1fS0KKW2lATi2coxcrB8lFOYM3HmEiuBmuN72gdfqANX7Yz17th8WWigoRVCvq/e2rRRAsah8C7quhPvjj6tjcjmlRFot9TBNdd7SUr/UdK/uT08RdFtfikYdPxJB+r4y/QBSSgQqISwAkIHyC0jQu/rBjVtUp7bhN2vIICDadpHtZje5TKIJDd8MMF1JENhIU8cQkKjAn31SZ//R08ylwI4apKXB1rk6bUNjoO4hi02GmhEqpke80qS9dTdjIzsoGxaj2QkeTyQwrBjnG/NEvvd9NpVddCtGa7XDfNTDSccxF5bRNufpdGycwGGsKZibGuHotjibkpswdZOh1BifzT8CFw7CeKK/ujv4EqnHnyBYJ2szTZ/Ka89xfHgHkVQW3fGw3nyXYFuK2MQuYmaMXCxHw20gEFyoXGB7bvtFYZc9IWloBlLKywrJyyWHmbrJ4xOP37iN/Qq+rM3lFLNuU729grC+1BcwU55hOD6M4ztr4bPJSJLZ2uxFCuR66iuF3F6EiuCjYr0P4ORJalObqb71MnJ+nszJ80R9iWVa/WijXs0dz1MF6IIAHnlEmY8iEeVD+NM/VZUmheiXly4U+lnLW7YoBdBo9LONTRMxOopYWlxTAqBMPtLQcA0dw3HxNcBT0QAAUoDVcpC1Fs7mrTRqq6QLi7TyScxGBzuiE/UlnrI6YQQQsX1sS+PCpiSZCytYDcH2xRq6YZBKaOTbAjwwNINSp86kTBBUinh2m9l6AucdF384w87cHtofu5/snz/HzpMruAsliprEFB61ap1YVZDKjeGsLJIczhNYaSoLZ/BllLkHtlNul8lFc4ylxpQAWz9rnp2lefD7JE6+RvnlF1j8qR8md48Ku5bT09iWgUiomP+CaCGDJoOLPtZWi7bXRtM0ntzyJJV2hfn6PHuH9l4kBK9HSH6QGfUVHbhX8GWlRzfz2MTUFcdxuZDXptMkHU3T8ToARPQIpXaJpJlkauDihjofJOJoowjLXd84oSL4KFhvCioUaL99mNbX/oh0qUqQz+PHLOzFJTRXYvpBv8xEKqWcv5qmnMz33qvMRgMDKst4ZkYdB8qR7HmQTuO4Ds2JYTwdooUKUQJMw1C+As/DCgLsZAoqNZVBrGtIw0AXAikE4IEv19LQA8DTQQYezeIi0U6btOfh6wItmsBaqSKDACl9Yo5yLKMJFT+qaeTKDe7rQHkgwUTNx7d0qkaV4bpLrONR25IltlLBaRfQ7TbFLQPIHTsYlib5d8+T2ZZHj80RJJPIWBRdSqTj4KRTlPQmBgaRZpn2+ACOIRgsdiCa5dVPTHE66zIYGWR7bju6pisBduwNJSz/9E/pHHyFekzD3DpGvFFm7NkXeM+uURrLMdmR7NrxcepOk5bbZCA2QC63A3dliaPNFYYSQ+zI7yBjZchH8+wd2ntZYXglIblegBmaQcfv0Pba15xRXzVP4Sq+rKsJ68sVsXts/DGeP/c8O3I7aLttFhuL+IHPPzrwj24bQRuWu/5ghIrgVnKlzOPpaaUEzpyBIKASNIkVKugBSM9H79gQtXDtFmajoRSAlP2s4mhUmYksSxWW++pXVfIY9H0GhqEiggwdO5AYiyvEKhX0Vht8iWsamJmsul67jVWtqXOEoB2NoKE6kxmeg0QipCodK7q30ALwA0msYVMeAq3SJt5oYTgtPKFh6ALhB2tDQgikJnADH9EJSGHi1j3OTyQZrLnkSzbCDfAjOqJap5WN4+MRxKJg6Jgdm03RQbxKQKm6RD7/OPqbNYLZCzj5jCo/LWGkbWCbOpVYQH04xfLuYbRcHn/7NnZv2orWDf0cSgz1Z4r1OnzjG7CwQDNpYQgYOjVPZM9W2rkEo4ff492nJ3hy24NksBjJjwLKfFavLLPn3s9RzZQot8scWz6GqZvkYjk+vePT1/9TudwM/DoyeUEJbV/6nK2cpW7XSVkpBuIDF1cevZYva904povTfOfMdxhPjbM5s3nNlzE1MEXH7yClZLGxyCObHuHprU+vOb3Xn/9Rzcavt9z1Rz3O241QEdwqrpZ5XKn0TTYLCzQSEWLpNIFtg64T5DIYto2dTRJvd2f2PUev6/bNRLUa/PVfw//8n2o10GsAo+sqoxgI7A4MZIlUKmi2qwrHGZpaVdi2WjlEIqo5TSKuzEatNq4Jwg8wfejoEEFN6CVgm2D6IFB+gPzMItG2j28AHR9XSOIuGEHXLiSBIKBj6RieJOZCMx0h3/apjGWYjTTYvWBTzUdJiRiu7mN6AXrHxdQNvHoLcfgtjiWPkdTiiJNHiW/aTLbZomoJrA7UExa0mxhNF9+zmdueITMyiomBE/jEIwmG4kM8MfnE+//gFxfVd2DbNI0Ax27jtytUZjvUtz7MznqSU9mteKNbEUfPoAFBLIpTq5BxdYKd22ksX+Bs+Sw1p0Y6ksb2bQ4tHLqor8DVCtXdTL3+C9ULLNQXaDgNzpTPUGwViRpRPr754+rc6/RlrVdG46lx6k6doytH2ZzZTNWustpcJR/L8/m9n7+s0LwdZuPXU+76dhjn7UaoCG4VV8s8zmZVtdCREWg2seIp3HgUIxbFTyUJMkkl5A0Tyg01w08mlQO4UukL72efVeWYCwUl2HtO4K4SQAg0JJGyCu30hnKIlkoQw3VxOx384irOYI6IkIhWS4WJOiqBzAPaJhg++CiZHugof4GmEekExF0b0S1dJH2IeAGGJhASgiBA0zUQGgEBngaaLmgbQKeFJyBWrDLkaURiSUqDKZLLTdqGZLBuI5yAjumyEJOklur4qQjupigrUZtDR/+c+88pm/xQy2fVttGjMbKZOMgAv1HnhChS1Sw+ldzDnpU46T1TlOCirF9DMwhOPM9KzGWXV8ArV2jFDIKRDLGOx/mlWY6kM4ynx6maOuzfTm62gLeyTCuhs+WHnuFQ5wzFdpGdAzuJGlGK7SLvLL2D53s8vfXpiwQNXL5cxFJjiZXGCu+svKNCS7M7eGLzE8TNS5rRX4a6XafYKnJ09ShJM8lIYoRiu8h3znyHz0197qIZO1x5NrxeGW3ObuboylFabouXzr/E9tx2TN1kODl8RaF5OzSfuZ5y17fDOG83QkVwq7ha5vEjj6hicm+9BSsrDHk2hXSU+EoVo+PQ1sFOWGRdU+UH9EpAzM+r172evS+/rFYejqO29zKFu7hS4puayvD1A2i0IIDAtcFzCXzQXJfYwjIegBeovABNKQEdaBtgaKD56r2nARICAmX2D1RKeCAgEqh9hoC2paN5AbpuYNgufiAxpHIyxzxoGwHlGHScNpvqOjXTJrfg4osImaUqdiSCEY1Q0ztEWgGWDQOVJu3CDMFQlHZskGY+iVMqcT7vMiATZNoeK0GVlBZDT2TJNT22j9+DkUhzvH6eicMvcGRrfE0ITxen+avpv+ILaY3BIMuxPTXyby1DxCQhA2TEIN12efXRDD9gZXh408NMF6e5kDTJRveyZ2CKXCzP8df+kqyVVYliQLVTJR/Ls9hYXKv4CUrQAO8TQjW7xkvnXmKluUI+midjZThTPsNCfYGf2PcT1/yppa00p4qniIgIcTOOG7jEjBhpPc1L515iR34HpXaJQwuHODR/iPPV8+zK7+KB0QcuymxeP5vOWBn2D+/nhTMv0LAbKlciM0nGytB0mpcVmrdD85nrcbTfDuO83QgVwa2gVFJC+8gRlew1Oals8eszj0dHlUkikSBWshlMj1ETJp1Wk2i9RWz/AaxCtwnNyooS9omEUgK9lo2uq1YOsZjyG3Sdv6DCPqWhITQdGQT4AvRqAzeisol1X2UPa4DnB/iWhgYIT5WD0H3wdUi4yh+ABMcAw4NmFCxf6Zy1XOSu/hFA4EscAgzLxDZ1Ym0HA3VNG3XOoC+JOTDQcEh7kJMa9YyGk9SJClVuohlxWMxZDFUdBksOrbjBStbCQLD/dJVTQzoTvsEL+yyKSZ3hhs49czYRLcAazWI4HbYvtCDnQzLLG8e+w5nETioLZxicq9AuLLAvm2Vu2xCbDy6SqbVoJ2MMl5pk8Dh+YAszn3yQ2JZNBDK4onPVqjYZXFwi0/JxUnGmY0Ui2RSO14/CulyhuqpdZbY6y7tL7zJTniFrZbFMC1MzSUQS1J068zVVv7E3i79QvUDdrq8J5qmBKSYzk5iaSdSI0nJbxMwY2WSWpJlksbFIqV3iOzPf4VzlnPIhRFJcqF7ACRweG3+MhJlgujj9vtl0xsqQiWbYkd9xUYkLN3B5a+6t960obofmM9cTlXU7jPN24+5WBNdbRvpGr3nwoKqf32tCc+SIyiNotVS8/+//vioFMTCgpOnoKLFOh9gFCfc+psxAyaSqI1Tq1rbXdXWuYahze+UieiUjel3HhFArAwChoQkNIf01803QnZlrUpl7gq4NX7cDhOyafQL1w+joKuwTASJQSsCQIB2QXQ3QKxLVXSisvU+2AzzdRjrdFUPv+9HUdTQPUjY4AjpRgRc1QdcwWg628NF1nbYpaCQjmB2bfNTENCLkGy6gEcQ0hooOzz80TDVhk2k47PSynN2awKxW2FMsslUbItMu4lcP09q/hTdapwhmffbMVKmbkjed8wyurGLNLVLJbcZoWXiWxanNGtaee3nr8S042QRx32W+Ns9zp597v2OxVOKRuYCTnTIiNUDMcdk2u8K72+rs2PHI2s/i0kJ1nvQ4unKUIAiYq81RbpUxNZOm0ySiR4gaUdKRNHWnvmbT9qXPQn0BTWjUnTpRI0qpXWLXwC4G4gNIJJPxSRzPIXFhiXvfO8dILaB08F/j7oiS2zxKoVUgHU3jeA51u85sbZZ7hu6h2C7yyPgjF82ml5pLzNXm1grpTWZUkeHDi4dJR9Lvs6/fLolk1wpdvV3GeTtx9yqC6y0jfaOs9w3E4yqbeHVVRQiNjCjBfeGCEvTNpsoQ1jQlxF97TSmP0VH13jSVzT8W6/sGfL+fJdxSDU1wHHVdXVcrkIWFbgN0wHPAk2u9A6xuTTlJz9avBHNvdUAAZvfYuAOuAY4OMR9MqQS65amVw1pOwbrnXoipIQFPHbvmXwDMQJ3XLYKNKcFoS1zXxrNM7MCllDbopKK0fIeUr7Nn3iFuS4oZaCRUbetBR5LUBG/uTjJv6ggSbJnL4EQDhPQZO1WlmREsDCUZLDRovnqBxb0uo9OnmZMmdS0gHklQ9T20coVTVsDio+NUlzskV0o8WK2y/7UZvnVPjDNahUc2PYImNNzA7dfWaUP9T79C/MJ5DG+RudwyJcun0V6h8G6FghUgpWRrdiu6pl9UUvp85Ty+73Oucg4pJUOJIWpOjZpdY9/QPlIRVel0LDm2ZtM+WzlL3IwTM2K03TbFdpFt2W0UWgWe2f0M//bgv+V0+TT3VuM8cnAFJ51gcv8nKS6dYdv/OsvyMz/EQiqB4zlE9Ah1u444fYbEc8cZLlbJ75nnY4/ez3sxj7OVs8zV5nhs/DFWWivU7BpHlo/gBR5IVfaiZ/aar83zr1/518TNOJrQ1kyHve5kt5sDNkx4ez93ryK43jLSN8p630Amox5Swve/r8pLnD3br43f66I1PKz6DA8NKeF/7pwyLY2NKfPR4GC/4qhtKyXQ7bxFNquUQ7mstq+sgGEgbBvRnYdL+oK3O8FHAFYATtfG39vfe+2jfhy6BxZclENw6Y+mFxjUO6Z3v95x+rpj1isB1r0Wnko4a+kejXiShU1Jdp9rEDtXIdmRuAZEXZfRQgCaQUZKtKSOVqkSDBrsG9xHpbCI06ww1hbMjSVItDokzp+mYrc4nfK492SaZszgSKxOx7eJGTEabgPXNtD1JKJaZdt8Ey2Rp5TQSZXKbH6vQvLeKaJGFDdwOVM+w/bcds6cPoQxU2b+7NswMsLOpol/5i3O5NuQTvGAtoXXvA6vzb5Gxa7wd/b/HUD5CZYaS7x47kUabkN1PfNtNfP2WiT0BKeKpzheOM6OzA5+5r6fodAqMBAbWGtl2XAbLNWXWGmuIKVEIhlNjvKFe77A8dXjTH7/25wTNsOZUU6W3qMaVAgiLvFX32DkJz7BTGkG27cZWqyz55WTBKPj5HbdB80m2T/7Kx774hchAyOJERKRBCOdEWZrs6w2V5mvz/NjUz+2FlI6W53lxfMvogud8eQ47668i+M7/K3df4uB2ACniqfIxXK3nZD9KBPebkfuXkVwvWWkL2VmRtUNWlxUgvrppy8uH3GlqqS96y8tqRl8r1lKo6HKRi8tqRl9r4Vluw3Hj6vrlUpKAfQUgW2rc7NZFT1UqSiTkWmqfZdwqeDtCW2JCgu93LE6fYG+/vzL/WDWKxp5yfv1x1yqBNbfMwLIept4TGAgGK9KnFQM4fsUPVWmIqYZZJsu5aRO3fBJNQP+/X8uUsjFOLv1GMVdE7iehtVeYTHhEDQKOOUOZ0ciGIMjjLcCOqUFkgOCdtqg7jQI8AmMGLbmk1mtkcgMoMWSpKRJcmKcVLLNaCXC9FCLmBGj6TY5OHuQ7dOrLBsZNuXTRG2bM+1FarrHprJHK2mRGd3CroFhCq0C762+x//5/P9JINUsWQjBQGyAhdoCpXaJTtAhY2XQhc5KawU0eGj0Ie4ZuYeV5gq6ptNyW6SsFMV2kfnaPG23je3bHJw/SKlV4vN7Ps+uwV2MpEaI+G9S2hTvFriL4PguBa2JdeocJ1fHGIoPsdhcZNvRCyTGtrF5872krCREu/8hL71E5anxvuM4qnwFUkq+P/t9VporvHT+JQqtArPVWfKxPNvz21ltrzIYH8T2bN5afIsf3/PjwN0djXOncPcqgustI72emZl+JdHJSTWb/8pX+j2G4cpVSfftU69nZ1WJiHJZCe10WpmOpFTPtZoS+NGoOi8IlIBvtfqvQSmTpaW1165lIQNXxXAaRtcBAK5lYNpef+bNxUL7coKZdcdei26bgjUzUbDu9frreKgf25XuKVEmqnRHkqnZzGeiJIs1lmMSjIC8rxFpdtB8yNYCEtEImmURRJJk3CjZZcFZrcTL1hIpXyc2V2eyJHEMjWjToexd4EQmTWXIZPNSCycRo2Z4mJ2AVcslPzDA2GoTbzCH6fp47RpH8wItnoVKmYFOHl5/nfaFIzQMB83J81Z+iaN2h61zDeyIIDB1UgWbYnWV+f37adgNjqwcIR/LU+/USUVTPH/meRKRBAEB5yvnSZgJYmYMJ3AQmmD3wG6GEkP86O4fpdqpkjAT2L7NYmORQqvA9y58DyklXuCxLbsNS7ewIzZvLLzBSFLN3DcPD6CVV3BMh1KrzGJ9gUwbtNFxpUS8Np/Z8Rl+OJEhsWP3xQ2LMhmYnSUbveeyDtV0NM2fnfwzNDRMzeRc5RyL9UV2DuykbtdJRpKYuslyQ3Uku9ujce4UNqJ5/Z3J1JQStM1mv4lLs3lxq8hLeemlfmnnVkutHhYX4Utf6jt1e8k7lqX2W5ZqKgMq8evll1WZ6UhEKYpeM/kgUIqgWu1fu+cHWF1Vf6zJpFIAvbG6rmoSLwS0m+A4qneAF6x5Zw3Pu6JAX2/KuRnW+wouVQK97T1l4F7hfr2w1bYJEdcjXaiTL3cYKnZoWALDl0ghcAAhJUYAnmmg+RB4Liumg78wR7xh06mV2VywSbRcylZAzJbsWZUs6G0KCY3lrEkTl8EmxBM5nt+u8V+Gllgw2uTqDi3hsbxliHIkYH55mrnGAiPHzjGzeJxyTJCQOtZSAW2lQNF0OT6qsexUCFZWWLV83tueZd5sc7p8mpgeQwYSIQRtp03drVNoFVhtrlLtVFlqLFFsFem4HUzNpGbXWGmtcHz1OEII4macSqey1oQ+G81i+zYdr0OlU8H2bTShUelUmK3NUrfrVB69D69QQFRq1DoVBh2TZNNh4eGdHNh0gC/e+0Um0hMEo6OcOfc27y6/y3RxhrrdUL+/sTGmBqZouk2aThMpJU2nSdNtslRfImEmiEfiaxVIdU3nVOEUiUgCx3do2A0G44NAGI1zp3D3rgiut4z0ehYX1Uqg0VA9hKNR5QCem7vY0bw+k3NmRiV++b4qD63rasZuGMq802z2wz57jeVdVx2vaep4UNsajX5kUCymzumep4SrRAoN4XtrDlrNf//H6Ani9auEa60OLsf6S1/tfEE3JwFoWxDxIOG//3gtgHgHWpakZHpMeAE5H1qmx2pcx2pLDMATUI9rRAKXthdQs1vMD8TIVFpsFRaVwOHNEZ198y7blzVmhjUW8xHumbVxyxXiiRwJIw96wKxfxzIsYmaM5nCK0vQ8+uAgK81lWtJFa7XYNLSbc26BVdkkLuK0TY0j6TZbTlXoWG2SxEgl8yzkLV7YrjFrn8eYmaflthhJjmAZFn7gU7WrOL5DO2gjpUQX6v/WDVyK7SJpmSZuxBmIDdB0mtTtOkvNJd5Zfofz5fPM1edYbizTdtsIISh1SuzM7USiBPW5yjm2ZLYwNyQ5/cQoB2ZaJJZL1PNJ3vvUfpxNeWJuk7gZ51zlHO2pBFtOVUiLPO24zYXzb7NVZkn83D+8rEM1EUnwzdPfhADyiTxbslsYiA/w5sKbzJRn+NjEx/o+gs1/a0153M3ROHcKd68igOsvI91jbEyFgr72mpo55XKqP8C2bcrEND2tHLs9H0Iqpcw3Y2Nqe6ullIfvq+3NphLkjcbFncSE6DeXh/450O8z7LpgGKp1JF2hGoDsLgVU9f33C/1LX19t29XwuNix3Lvn1ZC6Chn1gLamIogkKkqpZ7KSEpJNjy1LHsKHjAsxF6oRn0JUmbFdHRy/zUDbhmaAG4eRJZ94x6OoO1gmxFo+NUsQ9zzSNZisCpIdjQsZh6DTYFBL8ta4Bo7D09PgB8uUx3UiAxYPvHOahxsayzvGeP3pe7ErVYqyhRd42J5NuVMmHjiM1sqU4m08rY0VTRC1LKKmjmV4a/H8TaeJqZu03BZtv41A0HJbRPQIuqaDBE3TyJpZml4TUzOZSE2wd3AvnvR4/szzvDH/BulomlKzRMttUWgViOgRCu0CK40VEpEEWSvLydWTbM5sxgs8Ynv28/aUy1JziXKnjOeXqJ5XUUeO76xFMRU/uYfMa++QW+qgbRrn1NMHeLBr5lzvUD28cJjfP/T7lFtlEpEEDafBu0vvct/ofWzJbuFk8SSvL7xOPpbn8cnHiZtxLMO666Nx7hTubkVwo+Tz8NxzSlrlu8XKXnsNHnhAmXmOHFE+gJ4P4Z131LbPfEYJ9lRKrRCqVfWIRtW1ek7gSKRfaVTT1DmOo+5VrfazhjWtrzQ8F2nq4PprwtRHCeWe6ae3OrjU8Xuj+N1r+oCtq3pDPafztX5IPuBEIHBVZnFH614v6CsVDZV5rAUw2ICEozKZIy5kA2hZcGpQMNSUTFSAIEBakK9DuuVRtyDbhoEWNA0oJCWahCfOw1zGZilrImJRLMdn3i6gz3nkJ0bY/V4Zs+1wbuUEu0omxVwSc2iIWKnOZ//qPU6PRXDGYiwmoxRaqyw1lrjnfIcLGYNj41FMzSQVaTGu53i0lUPfPETcjLMlu4VvTn+TSrsCAtp2G8uw1lYIUVN5Z6WUtN02ET3CT+z5CfYN71tL/DpXPcdYaozF+iKBppLaqp0qdbtOIpKg7baJm3HqTp0gCBhNjvLUlqf4q1N/xZGVIxRbRc6UzzAUHyJqqPu9s/wOm9ObefHci+wd2kvrp36Yo50iFyoXiNlv4cxNvq820h+/88fEzTh7h/ZyYvUEXuCRiWZ4d/FddF3nmd3P8CM7f2QtJn9nfieFVoE35t/40Iu6hQXlbpy710fwQXjjDXjwQSWYGw3lWNu9G44dU7P906f7PoSeWWdoSCWGHTmiwkJ7DWciEeUo1jTlR9D1ftVQ0wRNw41EaKYTVC2wAx/X99UKomda8jzwAwLTVDZ4vW+nD9YNe30E0M3QE/oAEb8vvHsz+isR0C1Z4ffNUK5QJh4D9Vi7tgQ9gKStchEMXyW/6VKtDDYXJcsZncWsoBlRgj/hgKtBKQpRX6067luCXSuwvdhdgRjw1oSOLSQymWLA0RlpSIZnVqDVBiQfm3HIz5cx55fxZs9hCYOK5TPy3gL3v3SS8bemsUsr6C2bXMPlvXyAhq7s/16Ls51FqFSwPZtkJLlm/ml6alVgGAYxM8ZAdABDM9CERkSLMJ4eZzwzzlhqjIbbWKsgOpgYpON12JbdRt2pI4Nus5/AxxAGKSuFrulEjSjD8WHKTpmEmeBM+QyjyVGm8lNrq46TxZOcLp8mHomzM7+ThfoCUSNKtVOl6TWZr83jSx8Nba3sRKmt/F7TxWmaTnMtOmh7fjuGZlBulym2izw0+hAfn/z4Wl6BH/g8+96z2J7NQGzgfde7lfSS7z6Ke9/JhCuCG2F+XpmBNm1SfoFe16+5OWXmAaUcesTj6pjvf1+tFFxXCf96XQl921bHWJbyE7Tba+e7tSpORKO9fRIRT+A12iAkjq4joxZGo4XmeKqWUODjGQLp992wl2r4m1UGvXN7ZiDjBq7VU1IE6tE2IOn2lVVvxeIBEQm2UMpFAtJQqwIzANNTmc7tqk/Egdkc1NvqPKlDvtNNWPMh5qj3MV+tXnYVwT3nMeKZCL9MEHiUsxZvDzo4QK7iY7RdSlHJRF0QGA4FfRW/DamWz1Jesv1Ug6HTPqeGBFLX2V3WOG8aFHAItICRpk7Wb/LQEZsgbXN2yCBqRtkW3cb2/HYqnQotp8X56nkCGZAyU2SsDAkzQb4No4sNJs4dY3RHnLdTTb6jq+zhil1hU2oTs7VZbN9GCEHMjKFpmioDEctgaAYRXdUaemvpLR4cfZAL1Qu8Mf8Gpmau7ZurzTGRmsD2bQxhsFhfZLY6S6VTIaJH1HguKcJW6VQYTY6qmkPRNHuH9qoQ1MYiyUiSz+z4zFpeAUChVcAP/KsWdbtVs/awoNwHY0MUgRDih4F/i5I3/1FK+duX7LeAPwIeBorAT0spz3X3/RrwCyh58CtSyuc2Yky3hPFxNZsfHFQZwb06QCMjytewuNj3HYAyBZ04oYT9+laSvcSxbFYpBVAz/N45QuBKF3dsCNP1CZaWqGVimKYgWm/jew6BqRHxdQJN0BYeMSRCA+FfWeDfzIpgvUN5fQLa5fDWndMzJem9TDagZXZzB7y+Mujt9lDCXKy7jiZAaqryqeyVqAjUDk1KYj7EOqou0pk8DDQg4YHeVkqnE1Pmom3FAE1ro/sqY3rFdZgoeMQ9wWrEw4iAKXVMx6OsefjVDg09QHOgNpJGZix812akavP9XTH2lCX3nGvz1oTGABEeXJSc3ymwMwm8RoX8QpnlHWmCwQh1p87O/E7mq/PEIjEEguXGMoZusMmNsn+2SSQ9jm1pnFo8SvaMQ2dTmvr8WWrLL2IlY5CHUhw1cxcatU6NZCRJy2khdcnuwd203BarzVWePfks//Pk/0QgSFpJ/MCnYTfINH3K09/kEX0ILXeUs1kbUzd5qBonYzs0Ekf5fkfjvFan6TWp2lU0obFveB8vnX8JgKSVJB6Jk4/n2Tewj+/Pfh8v8EiYCUZTo5yvnmdLZstFv4n1YaQ3Wgb6RpRGWFDug3HTikAIoQO/B3wamAPeEEI8K6U8vu6wXwDKUsqdQogvAP8K+GkhxD7gC8A9wCbgO0KIXVLKy8S6fESsr0e0c6dqYAJKsMdiSqD/4A8qs1Emo1pFgnq9sKCOGR5W/oCBAfXcaKjr9lYGUuIMDVC5dwfa0jKm49OJDSDv2YPmBch338F0fWKNNi4BpgeOaeBJG6Pjk0EJzV6s/s1yNUVyLYew333U4mpWb9nKvGN0zTsNE1KOMg05ujL7iO4+g76De82v0S1TIYTa3rD6CXETFYnlKrNQ1FXX1ANIuVCJdR3PLiSr0DRhvBJQSgo6lsHJEYi1PDoRFZ11dEjQ0eGReR9HCDqGYMeyixZI5gcNUg2XqNRYSpgEDZtY0+ZE3uLhcw6fewc6KUFpapzh7fex1FhmngWsqMYjSxpuwybVOkNmxGc56/LZ+3+UlJWi3C6zUF/Af+1VjjTPAMtoQmMoMcRuPcvWl4/wnVyZdkJgthtMzTi8Pi5I57Nr/0mmMHE8h3K7TOFMgdfnXletJCUU2gUEgvOV8yAh2wyYWBSUTcHiWAarMscnS4P4fhM5mqKR0nHrNfyX/5r4gfuIZHK8vfg2g/FBEpEET295muMrx5mrzqFpGg+PPcxwYpgXz71I1IiqUFMk87V5dg/svuh3sT6M9EZm7evrLBWaBY4sH+GV2Vd4astTeIF3Wxa+uxPZCLnxKHBaSnkGQAjxVeDHgfWK4MeBf9l9/TXg/yuEEN3tX5VS2sBZIcTp7vVe3YBx3TyX1iOKxeDJJ+H8eWXvz+fhB34Atm9XZp/VVWXuOXtWnb+4qFYRS0vKfGSa6hq9ZLJuITk7laAuO3gxC2NqJ3ahyGwWBnybzEoVp+c78CVWx0bTdKTjogeSwBCqUTyqFERPGfSKv30QrpTw1Xu+2nUlyuErBFRMGHTUBduGEvppVykTB1XwTuueZKByDFyNtaJ4Pf9Dz3cQoBRGvq3qIOmBqojqC3V83VLbbBM6ploF9MxZ6ZYyKy0nfBqG4NR4Er3eoKkLthUDAk3j2DBUopJtFY8n5jQ6QlJKKtPU1HyL2aEIPhrtZITxlqDjBJzOBMQDnbwhyNVc9I5P1IiyKTlONgL7piuc2qwRrTikLhzhR5KDbNpqwc5JTpdO8+bSm4zMnqRoOSQ9g47Xodwuoy17JKsdOqMpkmaSglcgFtG4tx6luHWCXDxHtV1lvjbPSmsFU5gkrSQrrRUcz8ENXIQUSCRO4KCjs7uiU9Ed3IhFzIzTEE065xfRhcaJtE3KS+HiEjMDNi1U8Ma2U2gVmC5OK59FfJDBxCBPbn0SgcDSLc5WzrJ3aC/VTpVyu0zLbfHo+KO8Nvcac7W5tXN00a+1dCOz9l73tTPlM8SMGCPJEebr8/zBoT/gmd3PMJocvS0L391pbIQiGAdm172fAy41xq0dI6X0hBBVYKC7/bVLzh3fgDFtDJerR/TAA8oM9NhjSknYtnLaHjumhPzevSpfYHZWCe9CQZmQymWlOFZXlS9g506VjWzbLG3NETtygkjLpv5Dj2DOL2FZAecrS+w526AtfFLFMr70VZkG18cCnIhOEDGg46B7SlSb9KN7NoreLL/n1L0aOuBLSDQhqishqqFm6W2zbwKKAJ7s1y7qJZTp9PshxNatC12h/AwxB2Ld92hKCQQa1KJQjKvS2YEBMU8j0AIi3eJ4QkIrAiNNiAQejx6rcGHARMvFKWZdRpsBjvSIOjBZgtVowGoMJurQbnssJwSpmkMkIljIauQ6AbWkRVrE6BiSFSmJBg7MnGR+GFzf5b5CnKYhycwVqBs+c7GAAdMi+53vsppJcmLlBG/Ov8k+0SDWFpSlp0I7hY5espmNQ6XtU3caatVoRcjWO8zX5qnYFXR0BuIDOL5DwkzQ8TqUWiWSkSS6pqOjU3bKCAS6pjPkmCxGJaPJEVZbqwzFh0gJqHTKtFyBQKPttfDNJIMdnbeaiwR+QLFdxNRVaewd+R3KzOQ0lBPbrjMQHWAwNojMSRbqCzRd5VgejA9yoXqBt5fe5v7R+5kuTjM1MHVDs/ZKp0KhWSBmxNb6PXTcDrqmU2wXGUuNvW9FERaUu3HuGGexEOIXgV8E2Lx584dz02vVI+rt7ymBWEyZe956S4WPtlqqOFwyqZzKhYJyHvca0BcKYJpoh9+m5XaIlGpUNiWJ50ZIRTNQu0AxrZOqgNQ0jEDZSAIkOhBxfHD8983+r2W+uVF6wtpdFzJ6OXo2/V7oquF3ext0lxMRv1+qGvqmpl6ZaoEyF5le/x6OUMLdMdR22TWJtKMQ91RfBF3A6QzMZSHvanTcgIgTrPVR8DSISjX2SHdMO4tweMzFd3zmcwZLVsBPHlOO5mRHfd5NHhwbhJE25NoSNwIFy8fwQUZjWJqB1ukwlzeJ6FE2FwPEyirxTZvRfZ1YqclyxMaJCBxdJyIMzsgi6arFmy/9V17iNZpuk/ODJo/OSYpBG9eUmK6Hr0EhDi4eru9hoBNzApYjKnLI8z3qbh1DU0XranoNgcCXPn7gI5FE9AhJM0nbaSszSjzCGBncwKfm1MhEM1xoLSFlQEQfACExtQibzUHOUSGqb+ZM9Qwtt6WS2xpLGJrBvcP3slhfXKt/1PE6xMwYHa9D02muhc5Opiep2TVysRwCsRbBs2tgF4cXDlNql3ADF1MzycfyfGrHp973m8pGsxxZPsJIcmRtW6VTYTg+TN2ur21bv6IIC8rdOBuhCOaByXXvJ7rbLnfMnBDCADIop/H1nAuAlPI/AP8B4MCBAzdbFeH6WFePqNqpcqE6S6u6SiKRZ1O7RL63f2lJzfJbLbUaCIK1fsR4XemVTCrz0dSUKibXrWvUlj7J46epj+ZxhvKIdgfeO4jp+AwsLCLNCOlCE8MXeLEoge1ieRcXlrsVMcDdKtLAusS0dXkDV/IhuL0Q1m5Hs2j3uaP3w0171+gFEvVWHFqgupz1zEEe0IyomXzEh3jXhITshqJ2m+dUo9CJKlPTYixAjyrhX43CcBPSNtQstd/yVV6CocPWCnxrj04lLnj6VMC5nGCgJWlFlYlJD2D/CsxmIWODEUgayQhLGYOJqo8pLGYmU3QsSdP3OJnVGKz5jNgGhYjghdECW89XWU3peI6PpVtkA4t345LZc23qI3V8fBYtn++NCbaUJINtKEfhL3bB7pIygbVMiLo+puvz7rhOBEmxVcSXPslIknqnTkCAQBDIgEpQIapH0TUdTdOQSNVfYCzB5PmAkmximgZBvUEx6hPVoyQ9jab0mDDzdKpFXh8XRFaOcap0iq2ZrUxkJtCExtGVo+wZ3EPaStN0mwzEBpgpz9DxOwQyUK0/ZcBkepLZ2iwxI0bUVGGqvZn7mfIZ1gqjrz1d/k96amCKV2ZfUb4AK0vHU6uBmBkjZaXWjgv9ADfHRiiCN4ApIcQ2lBD/AvAzlxzzLPCzKNv/TwLPSymlEOJZ4L8KIf7fKGfxFPD6BoxpY+gWkKvaNY7WZ0i4GjnPpLB5mNm5gzw+sYvc9w4pM1A0qsJAV1bUquCFF9TMv7caaLdhz55+f+F6HWwbu7iIFo0RaXUo7pwgM1cgc/Issm3jpuJYno/RaYPvo8kIWuDjxCNEnGCtG9nVcLk4uex6EajZeDUO+aYyvxhBf9/lCFAmGyvoC/teAxyBEvS9cQj6Cqa3QlgLGe0+fNHNPtYg2epHRPlA3FUnSh1sA4pJnbmUz4Wcxp/cK/jFQ5L75wN0HwKpnMx6oPwL9QiczqnniUpAK6YzXoMzGYmtCwpJyZYymK5SHGNViAbgRiT7liSptstL+2KMNQM8IQgCH9N2sQPB9x8appZq0HSa+F6Nnz7rEGto1GMGhu/iOTZHoy3m/DYBSnDagc1KHFYuaU9cicHOEmvK4cgweOk4vqtKVKQjaZpOE1e6yvwjdHx8NKnR9tr43ZgLXdMZjg9T8du8NNpiZwkmHJNStMHp3TEyVoZH2wl2ugZGaoCDw1CzXNrVWXRNp+E2aHttkmaSbDTLdHGaB8ceZGpgiuniNB2vQ82ukbJSpCIphhPDZKIZjq8eJ2Nl6HidNaG9Prx1Z37n2me9UvvLfCzPM7uf4dn3nmW5ucxgfJADmw4wXZxmIDaAlDL0A2wAN60Iujb/XwaeQ8mcL0kpjwkhfhM4JKV8FvhPwB93ncEllLKge9x/QzmWPeCffGQRQ1fqVvbYY8x+7+uk6y7GwBCt/ZNYuQwJp8kpWeCxXE7lFpw+rQRzPK6EvOOo0hKFgnrdLQnB3Jw6futWmJ4mcn6WZi6Fr4GYmyc5V8KstzHcADcmaA3G8TRIzq3iGy6eruNFDAzPResloV2FGw0ZXT/bN6WaRS+nVabu1BVyctbP5QyplIaUfcdxx4S4rX5svdm/xftNRB5KGQSo8M62rpLIjK4JyBHKDGTrKpJIop6PDsGLW31ahmBLXdnDl5Ie8UEYrysFMNjsd1srx1WWczkGLUOybdVjPg3DLZWNvJqAVAeGfYjYkAQKSajFJEYgmahr/OgZkz+/32KwBdl6g1rM5GBeUpLLOGUHiUTEBc/ts/jRkx5jLY3VmM98VsOTDucGTXTNQwiBjhLgl1KOwxuXKAfNa5I1sowkR1TJC9/G9E3QwPM9NDRiZgxTU1VqNU0jkAHLzWV8fGKJGHMjwxTNGC2nhSY0qprG20MRBmIDdLwOw6lt7EuO887SOyQiCRbri7y7+C57h/eyM7eT1dbqWpROzybfoxfl03SaJCNJKnYFiWRHXpWtaLmqJLvruxxdObqWPDeRmqDttS/7+9qR38HPPvCzF4WQPjr+KIVWIfQDbBAb4iOQUn4D+MYl2/7Futcd4KeucO5vAb+1EeP4wFyjW9n8vnEGYvfhrKv/s2aTlKieBA8+2G8uUyqpaKGhIbUSqFZVrkGppMJOu8lj7aE8S9tHEbrA0CMMLRVJlprg+QSWha8J9JUCi4NRRlM6ibYPPhi+xA18jG7TYP0quuDSekDrt12O3vaewznuwFwaxhpXXlWsZQ6j/AKOrmbqZqAcxDODsKmqTDSG2+09sO5ePZ9GBLWC8TS1koi66lptAzrdxDKne7DV/UCFOMzkoBaBfFtSt+ATp30y7YDhulJIrqYc0FEfKtGuDyKA5QT4MmCkBc/dm+DvHWySaUPJgn2+ylQuW92wV0sdb/genUSM+5bAjQW8sdng3UkLJ5NkrjZHXMaxNIuG3cDG5njKZPlewa5SQNYRtOI60zmdVtLC8DvYnv0+s8iVFINAEBDQtJuMJEaodqr4vvIHWLqFLnQsw8LzPXLRHD4+NVv5DgYTgzScxloyWtpKo6Gx0FhABhIpJSvNFWzPJm2lcQOXLdkt1O06Q/Ehih1lf4+ZMe4Zued9Qnd9rL+hGXT8DslIknKnzNTA1Nrqpek2GU+Pc3jxMLlobm3FcHjxMA+MPXCFX9jl7f495RJy89wxzuJbyjW6lV01yiGLUhyZjMoXmJpa60DWSkap5Szc+AhadoB8NoncsoXW8XeRp99lYfsQ1o4tZI68h501sBptAhGgBxLNtokGAW7UJFdo0I4IGrEolmaQqnbQpVQtB7oZW1dyEAd0wytRJpLrrQ3US+iydYh6kOtcWRHYmvINJN11G7tC29Nha9fW3TAh6YPTDQ+N8H6F1PMb2Lrqa+yYaiZ/dBgmat1yE0Axpmb6p/N9k0+uAxcyPpNlqMQEpYSkXoG0qYR4y1S2/pijrjNVhoUkHBwPOLY5wh8Df/vNJoNtOJdWK5ndBXA0WEoLpIDAlwxXHDqGIGZl2R4bIDV/jhecVWJJNRNv2A1UvjJ4eKxEBUubAiwsTCPA9V1i6ORiOQqtAoEfdL9zQdJM0nSb7/uODQwMzUBKVdL6dPE0lmERNaNohoYlLKp+Fc/3EJoqWT3bmEUGSsk4noMvfSxhUWwX0YXOhdoFYkYMIQRSShpug82pzdSdOrZnE8iARCTBYHyQveZeUlaKLdktHNh0AOgL/9nqLLO1WaYGphhNjK6Zaj65/ZMA74vgObxw+OLiV92HuOkiKCEflFARwDWjg64am2yX+2Wmq1WVLDY8TDOXZEarkSubRN0IncDjbCrAr18gvXUTicNFOnPnmBtKkvrYbjadWkR0bDRP4usCoesILyBS75DWAmbGYxixBCv5FL7vsvvoEtG2i2doaFJVbluvDHqRRAIlTBHK5i2C66gS2j3fMSDjQHxdtM/lzjUDJaybhsr+1bv1grRuOrLsGv3j3RDStN3vi3wpAuVr9wRYGpzNqmt3DGUz17qTZ8tVjuSZAUg74ERNvrZb8NmTLoWoxJCSybZBRA+oWwFby8rcZHqwkIJGBKIO7C7Cy1tVKeiTW2P83pDJtlWPgWKbsarPbE7wwKpO0he4IiDwfXQJpYEYy6LDGXeZocwA985VaCci5G2bBd3jZE6Fs8ruPwAbG3yIRqLd+kRt0mYa13Cp2BUMDBWfr6midAAODkb3z7QXDaSh4eGRMTM0O00s08IRDhEtgi99JhITqtGNFHiBh0RSs2t4gYffrUtebpcxhKHKYnttnECFoAohiBpRVporbMlsoek2CQjQNX1tVp6P5S/KDm44DQzN4Ez5DHEzTsZS5SbWh3OuJ5ABD489zFxtjmqnSspK8fDYw6ofcshHQqgI4Jrdyq7Y7LqNyiSemlK+gGZTNan/xCc4L4pYp1tEvECZefZOsVI5T2augIxFWByO4QkdZzDLe9uHcCMGg6fmQCgTSKzjIwKfQACaQEZMpISBC6vorg+eS6Cphi2B7CZxrZvxe5oSnpFu3Z1mXKMpAjKdvqO2N8O/VMgL+uYUy1Ez/l4eQS/MU6w7dy36R4N6QiPbBs0L8Lqx/joq8ifRAVPvVy5l3TV6KxAdpSSKETibEiylJcN1VbH00DjsW1UrlLYOZwbgz/fBuxNRmmmLltviM0KSdFX8fwcPMxBsrqiV03wcMh3YXoZWFKYHlaP5E7MGhSGop3wqcY1DEzrmlhzxeocfOOOzv+AyVnFZjaNKO6R0FmIeSzmDwdggm4IEmxZmeHmzzXLUxGwHfGxe4/vjPuVLbPwRI0LSSK4VZ7MDGws1di/waLgNQJmHBjqCrQXBYEdQjsF7OZ9GIgICDGkQ0SIQBUu3iOpRWp4qb70zv5OZ8gzxSJy6XVclKJAIBE7ggKYUn67ptL02US1K020iECw0Fnhi8xOcKp1aE/Cf2PoJmm6TXCy3Vk20ZtfWsoMbTkNF9PgdZquzZIYzVy3rkI1msT2b/cN9527Pp3AjhFVGN45QEcCV20vu7/9QLxub/O7BvklpbAzuvVdlE6+s0KRK8v57Wfnkk6S/9TLGUoGBmXPEi1V0uYAYzZO0EpwYS7Oid0jqHdoGxIXAkBLX1DGiEYTQEYHHcD3Ao0202cHwJIk2BEI1anG6sfq9uj4SNYs2A2UWQodKQidXCdYiddZXEu3VBOohUWGcHkqAxvy+4O9lLveUSU8p+N1zErbEMwS6q8xCgab2W546UHQ1ju9frHh6Zi43ANuC5RQUkoK4JxmvKid0OQbzaSW8693sNisS55GlgGrV5XTOZDbl8pNH1RjbJgxXJVJAxYKsrTKei3G14jiXFwgEY03JVBneSguEEGzLbSMfz7PinsSOVPizBy02FXWmigGb64LVgRize8cYHNrEXG2e4VWflYFBGpEVYppO29IIpMvOknL2rrf5BzKg5bYYTgyzM7+T1c4q87V5kKwdo6GRbLk8NC8QiRSFRIukq3FgPuDwuEsjESEXzxHIgIgWwfZsPj75cTzfYyAxwN6hvSTnkpwsnqTQLGDqav3V9ttIJB23A6ZSBl7g4QhlxpJCYmgG5yrn+KFtP4QudGzfxjKUoorq0bUV8cG5g3xs8mMAa7kEvTBRuHo450Zk/95ovaKQqxMqAvhg3crg8ial0VEwTYJ77mHVs1XN+HumGP29P6Yjm3SGsjSFJG77nPncD5AQbaIXjlMay9LMJIiZYAcufqeNMDSseApjtUiyITE7Dq4eEHGU1DcEBLqGGXTND732lHRn6rqGbgdIDbJ1l0ivJIOvEqxc+iaa9UXlejNzze8XgXMMCHpZuuuOXQv31FTYaOBD4AW4BnQiGhoS4UtMDfwI6N20it49oB9JRKCygitRpbgy7QBbU6sbCWyuwPQQ5FwdU9Mo5GMMm0nSdQdruc7H33OoWEDEoK1JIo5PxIdCTEUJZTuCQCiH8ngNtlUEuyoG0ZbDRMPBtBLEpM49ERNyBuVCknrWYlU0Wc04LN07zGgdhheq2NLjwuoxjLZDcw4W9w6x1YpT7VSxPZuGCSNtrauc/e7/i4HjOSTjSQzNYKG+wFh6jCNLRwhkgClMAtVdiJ0laJiSZDzGqJGh3CljBw5TJUF9y04qnQp1t47dsfGFz3cvfJfB+CARI8I9w8qZe7xwnGQkScfv0HE7aN1/QghMzaTpNYlq0TUndK1TIxPJMFOe4e/e/3fRhb62ErY9m7OVs2tRPjEjxnRxmkfGH2EyPcnRlaN0/A6pSOqancmuuMK+AQEeVhndWEJF0ONGu5XBVU1K62c9+fMLFB67l3a9hN7u0DYkrYiJOTvP6o8ewNszhhnNUX+7SOrQSTqWRiubwncdRgsVNNPAFwGGAxFpoAsPT0gQOn7UwnA9DKlBp7NmpvEEaE6gSjaYJlHp4+sSX5P4ERXpKiRIV/0IHA3sqCDZkgSiqyR8ZfOne82ltArHRPbzEyRKeGuBMieJQBL1oZrQwdTxIhFkq4HXNV2luk2LfQMaqAqkAEJTCuDcgM5sVjBZ9Ii6oJmwklHnDtdhogrnhmDIi+D7Jk8ea+DpGsWISbopGW8FnJlIEG3YBK7PeApcE5ayOuNNjXgnwPQDdCn5e4cDXMNhOS0YWA546r+XeW/U5NQOSSAFj85C48CDnKPCcmMFBNSzBqbrURcOiWqHZtzk5I4EduCQNLNoXfVmVYpUoqpjnI6+Jmx1dB4ae4jV1irnq+eJGlHcwFXKQkJEjyAQ5DoexRjECchH82hCpxTAYCtgrlNVGcSejy98InqEmBHD8R2iZpRis8ibC2/SclvKuYzAMixM3SSmx6g5yl8gZNdUhKoomtWzSKHifi3d4uFND5OP5ZmtzjJfnyduxNeifBzfoVAtsG9oH2krzfbcdqZL0yQiievqTHa92b9XMv+EVUY3llAR3AxXMSmtn/U487NoW7awddeDAJytnGWmMM3YaoV7hu5Z+6EXxrIsPrqZydkaqXIL3dWxsyn0epu4p9FJCITrEw0CTM8nMAWGJ8GMqJBTQ9XvIVA5AAI1mw4EiECAoRMxDBrCwzM94o6kGlN9gl1TmWY0S2Xe0vULTA8JhuqSTEc5e2tRFdYZVUVTsSMgNUHcldiWjo2PZ8DqUAzDh3jLxZSCxayG8H2alnLyJpx+eGpHg9kBnaWcgWFYzAzrTK2W8TVVWvpkXiPbgZXukmcmL7CqAUM1D0fXWBxLMOTr5FtL1CMabq1DrOOzasGr23V+aMYnUpEsZgVxHXYvSywPWhFBKyoYbmsk7IDFlCTbFGxZcWh4K9jRGOm5VSJbEgwnh3ACF69W5XTM5siWOFFjUjluPYs9M1VWa0UiyTQTWo5Vv8qJEUk6kkQTGi23hS51klaSUruEF3g8OPqgso2bSTpuB0MYeNLDEAbtWISk61Lv1GlEGviBR9STlCzVzcz2bTzpETWiRPQIAQGO7/Be4T0Wa4v4qOY1CPqrAaHhChdLt9A0ZQwMgoCklWRTehMZK0PDafC5XZ8jbaXXBPlCfYFThVPU7BoSyXh6nFw0x1hyDMtQUUhDiSGe2PzEhpplrmb+CauMbiyhIrgZrmFSWpv13HtMKQhLOcPuG7mXbeRYztepBh478jsot8ucHPk+Ozppzg/mSSyXmSj7ROptVqXDWNnFdH2klMpJLCFwPRAaumZBp4GrqxA8oYH0u92LNYEmUUlogURzApwkNFIJmh2btq7hRXziDY+oI9G7SWSeAV4iRqDZdKKCZkKnYYHVdhlsqp7DCU9guBIfSS2m4yQsKrrLu7tz6JrBlnMVzuYDFjIRlvIGnzvhMZP2yDV8IhoI06Qa9cg3JU1LoPkBnbjGYGDhaSBcFSU0VgkQumAhA62YwXLWYGdTI+6CMDXswCUajRONJElValSiGic3RRgra+TrHi9t09leEWQ6UEhZvBWHPUsunQhE2z7NqIbuS+K+QCdGIRJlcqFJJGNRu3ACI7YVMThIpTyP1mgztzXGRGoCT3pUO1VWoz7BthTxc3X22CbG4CZmxmN47hwyUP2B01YaL/DIxXK4vsvm9GaklMqUM3IPr82+RtttqwYzRoQLozqPXDARRpK6XWNIJMlam3h+qImhSVVSQggykQwdv0PLaZGKpKh2qjScBkkzieM5RI0onqmK2bXdNkIKsvEsAoGneei6TsyI4QYuhjDYO7yXj09+nEqnAihhfLZylgu1C2SsDKZucrp0mqSZ5Ae3/eAtNcNczfwTVhndWEJFcLNcj0np6afhK19RrzMZqFZJNR1SX/x5dnYbhZfaJf717hdINGymztVJmylIuNQsg4qokjBiZFarqtqmlDjZBPF6G9+M4JkauqbCQyUSqXeduUK9d+MWmu1059MBdVMgWy1qUQ3T0GkZBsmmKnQm6ZV10FkYjGLUHQoJwblNMfbN2SRc5ayNdU06HRN8XcPpRgOd2JJgIa2Tr3VYyRrUfIdyXJB0NQzHY3PT560Jg4jjs3/FJe8pJ7DheuyogF+o4SVcihmL8YKNJ9TqwY5I4q7G8U0aw1oSd1OKxJLz/2/vTGPjyq4D/d169erVvhf3VSQlkdrV7Ja723a3nXa3HQdOOxN4ErQTO3ESZDwzCJDJIAYMzI8MEnjWzI8ZzIw9WZxBgEniJE5jbEOON9kORp2WuttuLW6RlChR3JfaWHu9uvPjFimKokSqKZIieT+AYL3HV++dOiTvuefcc88haWeRtoMGM0xbwMViCTwBC49ZxmkUcDgk01GTfzjspOoykB4Pv/h6gduRGr6ageV0YAoo1yqESzBtWTTlBabDTcZRI9vZTNtMgZv5UTyJCNmBftqnbvHkDyboKljMBA0u9nl4zZrD0RvG03QYw2HQRhPeQhtD80N4XV6K1SISSdQdXS4FMZYeQwhB0S5iOkzy5KnJGqVqCWekhQm/yZO5MJWFOfyNbVxPOKiWJ7CqRWpV9dv0urw4bDW7NwzVNtNrekkVU3gtL12+LsbSY6rhvXRhOSyCriA1ajT7m3EZLoJWkFK1RHe0m+c6n6NQLTCzOMPZ4bOMZ8ap1Wp0hbuWjUnQCuI1vcilHtpbxIPCP49inUFzB20IHhEPTGXr6YGPfhT+5m/g4kW16/jnfk6drxP1RHn+9M9xXX6LwK0LuCqSxbCHiWoJq/UI1Ys3WGwIg9OFKOSR1SqVxjhuDMzFAotBD57FIkZFLRrbhorblwSI/CI1Q/U5TpmSvM/JhFvSM1/jRwnJwflafbAWlF0GLqdFLmBRFGVyYRcTcRem6abkqVERNVxFG7PsIGUYuMuqIurtjhC3QpKY4UfkalQsPzgXcZYFx8aKLMQNrrR7OH6rSPdCjdc7TMK2pGu6QsEJtmFgSpv2hRpiIcd42GDR52TcL/EWbfxVBzdiDrUj92A3taSTCA4C71zlUArsWIqCE5ItUdw+J88PTVJzukk3NdEhc+QrEqftxEp0UotNs1gqE5+tUrBcmFVJQFi4DBt/JEGFIiFXgIpbUDlxiKHMDbBcTPU14hub4mcupFmw4Jq3RmBR8tOvV5g9Iplr8nMzfRMDY7kr2Onm07zY+yJfH/o6l6YvYTjUesFcfo6pxSlcThcJT4KQFcIhHPgtP4YwOBA9gOWweL2wgLOjBb/bS8ku4UPNigWCYqXIYmVRPc/lZ7GySMQfoSXQwkR2ggZXAyErhDPipC/ax3x+nnQpzenm0xxvPI7P8vHO7DsMJ4cZiA/wYu+LFKoF3ph4g9PNp4l5Yrw9/TbpUhq/y0/Cm8BlqCyl+eI8QSu4pf9T64V/dJXRR4c2BI+AdVPZFhZUMbqXXrqzljAzo86vyEwadB/ANS9IdTRQRWKVbOJXp+lIlnAJF/mgi6mDLUTmslQyaSKhJsxL16BYwmlaSFGg6lR9CmwhsMqqII9VgppZo0SZbMyP1xegLTVL3hS4hKAqbPLRAHalhgRmmuMYtQr5fIWLR+L4TB8WTsozZbBMikGboO0nGbbI59MEshWSnQnylRQNCyXmA2UcJai4BUbVIEWNeUuVzRhqddM+scjx8Sr5kJdvDHpIegRnhnJ4qw5c1RpGDdoWbDI+g6THwU9aDKp+D3/1bJSnihHaaz5KrirexSLXWkOUS3kaijmiC2VGz3Rhzc0z0egh1tCFJ5slmM0SiDbhD8TJWjDT6KdzPsPVPiftBZOWmylMKfnR8Wa6bZOOmptMKMCN/gDzrgoll0FbycWc288TI3myXhMzEiKfncAVCONxtvArGSf/JnqDgBUg4omQr6qWkT2xHiSSqcUpuiPduAwXqVKKyfQkIXeIVDFFsVrE5/IRcqswz7GGYwgEByIH+O7od/GYHmzbpmpX8Tq99Cf6mVqcIuaOMZefYzY/S9gTJlaLEfFG6Iv1MZmdxJaqPaXL6aIr3EWjv3G5GuhYZoywO0xHqAOP6eH57uep1qrMLM5wuvk0zYFmABK+BFFPFCklpmEu7y3oi/TRHmpf+x/iEaHDP9uHNgSPgHVT2dYpYbFE5PYcJ3vfy3i0HeP1C4RmFpBlN2alhCMUwV+r0ZW3yBzpxvr+OQJXhhHZRUSxjE9KyhjYToGjBkLaFFwOAkWJMNSGtKrDJpS3WSDDgXnJTMikJ+vCrNbwGCZ50yZWEkyEXKR8bsxaiAOtXfjfvEzA6aWzYOKswGzQIuszKZcLGDWBIxiiGAthjc5RKORJBSAT99G4KGlNSgyc+OZthhMlZmJ+ZsImi07J6OEm2m7M0Vl105Mp47bL5C0V4nIA/qJNV0qQjLiZFS5ivhjJg710NJ6g/1vvUDlUwL48SagoKYS9XGkOEHrnJmmfQdZXI1osEnYFSXa4iSwWmXMu8L12C19bjKZuN8/etKlVQrxz3CDjhlbbjzMlKba04Dh5lD6jxPjNc3SacRx+VaQtmipht3aSLmU4FD/Mk82DeJ1uhn/0PZ585klmcjPk7BwO4aAz0rlcx6cr3EW5WsZn+Wj0NzKXm6Naqy43W3EZruUdyLlKDgcOmgJN/MFP/QGvvvMq8/l5pnJTtIRaiPvivHjgRZLFJG6nm4XCAg6Hg3xZbSgLWAGe6XhGeSCGQX+sn5gvxuXZyyRcCZr9zSSLSZKFJLP5WT7Q/QFe7HkRgLPDZ+8Kx7QH25nITDCcHKY/0Q+okE17qJ2+WN+W/U/Bo0kz1WyM/WMI7ldd9BGwbirbeg1ulm+UIhRrIhRthjlg7PuUyk5SssjioT7cxQquS1dwVXJ4iw6qpoGzpmrPIEB4vDiqqlOLWSjhcLqQPkHRNHHlCrhsibVYpuowkT4vfqcTV75KwXITnS8QtWsYbjfHRotULYPpT/wMs5kJAp19+KWLTLJIZGQCQ5jM+Ry0FEwolBkPFCjduo7b8nKtHXLOMrF0hUhZEM/VMKWDibgHv+kiIE2mo15i4WaMSIzOXIKmyUVipRRFh4OK5cApK2pjnIRA3ibTEOZ2m5cPLoQ5FT/C4oUbtP7oBnkTrOYO0o4KlUKeci5N2JZUwmG6yzVKBqRb42TKs8QXCtxscWEYBuVqmXRbnAuHGpiPdC+HOl7ofoGT4VNMf/tV0naOsC/O6eAh7Eyac01FwpZFJuzDun0LXyrFs44uwtffwg74mLfKdEe6afQ3ciBygJGFEUzDJFlMslBYUJMDoTaULe2ircgKxxqOMZ2bJl9RlUC9Hi9up5uIJ8JAYoAPHvggJ5tPMjQ/xPnb57EMi4Oxg4TcIdLFNNfmr+ExPQwkBpBI0sU02VKWp9ue5pn2ZxjPjJMtZ2nwNdAcaGY6N43P8hHzqUqjyWJyuSQE3BuOCblDHGs8hhSSq3NXKVQKDCQGllNLtxod/tke9ochWKe66GZZN5VtnRIWd24UJp2c4lZlHsf0ZRLlNIH2NoJuDyW/m6w7T7ipiUSmhOw5SDo1hTNXwJkrqNLPpgvT66KSy1D2WGTiXgzDRXB+kRqqbIXEJlqtMdMVwpfJ45A23oZ2zEoaMhlq1QoO28ZrhonN5ghNZRHRONWJ2xRklVzUT6QsMTI2M16buM+L9MJYg4NwqIlTPxrDn6oyEzRIhkxKUSeNk4vMCIOQP4YZS5Cwi7gGTnLI28jIEUH+y3+NXaviFA6qbpOyLZEOgVGrkfMavNMTIuKOcua2JHEsTtuTTxN640+ZnBzGbo+y6KwxI2cJFisYITczRzpoKvipUuZWeY5IEUqGE+ehw7w30cRoapSSXaLB18BieRFDGMuhjkhrD5Gf+dTypEE01MgMNvA+q8oPbv2ATH83H/rKG8iSm1ysjL9UxDO/gO9omNLsJL6GFqZz01hOi1K1RMKbIOFLUK1VlzNvqrUqpmEymhrlcPwwhWqBheICSOgL9RH3xmn2Ny97CEuD4VKoZKn4nNPhpDPcueHdtGeHz9IWbHtgjZ+1wjH5Sp6eSA/vbX/v8rlr89eIeCJ6dr5H2B+GYIOhmXfLWv885dkpTuQjcPmsqmuQTKoyFPcpYQGQbIsz/I2zmIEI/nACmb/EfClHpHWQjlBU9UfujcE//AN0xxCLi1QCAUqWEwMDs1Am0xRBEqFglzCzOWrZDLZdU20uDUHN6cBhGHjmU5SCfoI1E9fNGYqNUXwH+sA0se0iRY8T/40xcJqUz1/kHW+BWshF2QxgLZQoNcYoh1yMNEYoHu3Hn5viyKyB0WBTszIEK0VC6SqZU8eYOCbwpJO0HTxBKNGG3dbCW7OXuC4XsA88QeE3f42J8f9Iw0QKp+ECjwBZo+AxudEXI9Z8gI+kGjh84jC15oPcSo/hDkuco3n80w5KUTctE1nCOcl4e5iQbRLqP018IUXp6g9oCDZz68Mn6G2PkyvnyJayZMvZ5TLG94Q6VmSCtRQWGLv9GiMz7+ByuDDsGre7ohxKOqguZkgGy1Tf8wQNrhLmyE2s5m7m8nMALFYWebbzWYKuIKliiiZfEx2hDubyc7QEWvjk8U8ymhplKDnE8cRx3C43EXeErkgX/bF+qvLuImybDZVspMbPWs+IeCJYhqV38e5h9och2GhoZmQEzp1T9YKam1Xa54rMnvux+p8nVhCcHJcEoxbE6gO/EGogLxTuW8LimpzDHjyN//Y8ti+PbIjjKBRIz9zCcfkS5elJCmE/tEUIWAY+G9yhmCrXWauRdaUQPh/umSTlQ51Up6fwD+UwiqV6n2SBM+DHNC1clTJFp4/5JwaIXbiCr6EFt8ME0427bMPAU2Qvv8n4xDU8NUnZcuJ0OJgzcoi4j2JLhNvHO7ltFugL+umaNMibgrDPT8YFox2tVHKLNEaCZFua+MCkm/iHP07aUeHqrYt0OeOMHIwzl58jGTFwf+6f4/sff4UvlUVWbcrlAtN+wdvHEvR7uxhIe7F7e7g0cwmP04P/QC8FA+w3LtKcdTPtdnHlkJtyOMjhxGE8vjBl08Jq/CnME08z0NC8rOfJxUluLNxYboRysvkkgy2Daw6oUU+Ug7GDfPUnX8VtuEmUTbzxJkYO+Ih5EswX5mjx+2jIe/hkxyl+6PMymh4l5Arx4d4P0x5UC6o9kR5mcjM0B5rpT/TflVV2pOEIpXo5kiVy5Rx+571F2DYTKtno4uvqZ5wdPovXvLt6nt7Fu7fYH4ZgI6GZkRGV6x+JqMbz6bQ6fuWVhzIGgAo7Rd13eyBNTaqV5QM8kFQxRSzRRL6hGU4fpXDsELE/+UvM8/9I1uem1tKAw+ujYgimC7M0vGcQ/9ANJXsySb4xhKsiKXW3Y5kear2HcWaqOG/dxoMDLCeUKlCtYUajmNE4AWcMgjEYHVNB+bY21UHNtrkehtAEGKYLn6xSrdoEbZNJd4mQbbNg1QhXVR56THpwhCMIdwaXQ5U96G07QQdhov3PE2opwvAw6Z/8I6FYiOr7nuZQVzuHUINeqaHE8O944dw54mmbks9iLuLiTLyNU4c/iK8XLqdH8FgePKaHcmsT/vkkxcMDjDb4cJseQqkZPCcHCfoTJNPTXK/NEys5GL3wLWLHzhBr61M7fIXBKyde2fBMei4/x6mmUxjCIN4yjjdzE2fVJl1M0RPppc/bohot9J7i4OkzvNz/8nIW2VIrRcNh8HL/y2s+c7uyY96tR7Ey9JkuphnLjDGbm10uR63DQ7uf/WEINlBdlHPnlBGIRNTx0vdz5zZkCO5iIx7IGt7H6rWGclc72ZY4pRMHCPtiGP4ApUQMKlVq2TRTURe9hw+DaYLTSdVrU8kXsAcOIcoVnJMzVAIeXA4DvHVjlM1CJqM2tsXjMD6uzkt5J621pQUyGYa7QxwptCFsm0h+kVRqCttl4q7WKNcq3HLm+Lh1hGPRYzg6A8xOX6cSayBidHKy4yAB06u8lWRe3b+3l/FIjkjNwrp+CzsUxI6oksWFaoGnn/8kFw4e5vzsVQAGEgMcb3mCUD0Ft/A35/FHGpFOCU6DSmsTnkqFmJAMdL2HTEOImzLN1OQIzvOv0/X8c4RbWggkJ5j64XcoPl2kqWPgoTNPUsUUfdE+Ls9eZrYljGdiFmt8krRZpSl2XIX9OjvV3xkPP+BuZ3bMu/EolgxVppRhJDmCQzgwDZMGf4Ou+LlH2B+GYCPVRScnlSewklBINaZ/WNbzQNbyPr70JfoPtHF79G1Mw4Xdf4iFgS78mTSFgYP4/Q3YS60yazWs7CJXPv5+entfgrNnIRbDX8ow+sP/i7tmY4aDlKhSXGzDX6jC7AIUi2rwd7nAMNTgBaqnMig5nE5lLF56CdMxzJvdnfS/PkLV2Yw1E2A2M4mjVMZ1/CT/Kn6GlsHnVS+GaIFAvgKnP6g6tA0NwfS00rsQyhvy+fAvBsnbZXzCg3VjjHwktLywHvVEebHnRZXKuJTlNfz6nSyvM2fID1/Hv5DGDgXIvv8MtaaIalvZfZQgcIx2Rt66QrVzgFqsDYBIrA2v04dzTnL0mYcPq6yMrY9lxhg61UGT16B3ukzAdqg2pU88cdff08MOuI9zdsySbF+9+lUqdoWEL0F7qJ2QFbpv03nN7mJ/GAJYvxREc7MaCJc8AVDHzc33fw+snZa6ngey2vswTZiaIjg8TPsHnmUmO4391pvEUymaW/sx8zOU7DJupwWAI7PIYix4T1ZSyBeiu/s02b//GpVcFjORIG57sZ59vyo3euuWGvTb2pRMY2MQCKjnd3Up+bxedf5XfoUjCyP8+Y//nHLgKCe+8SaOfAlvrIXDL/0Czf1Pqs80N6f0eubM3bo4depOiu7Zs+q+QEdIlSzG6cafzqxdsnhllpfTCW++Cd/8JgdPDvCPbVEWE013FuXbojw1LpUsdV1XZyaRT7/nrl+TKxAiN/UujDp3ZsQ+08eRxBHy4W5y3Tka287APpkJRz1RWoOtHG88rtKV6+i1gr3B/jEE67FGPSCSSfjpn77/ex6UlvogD2S19zE9rb7bNoFQgkAoAbGCahzQncDxVorR5BzFUBRPvgwTk+SPdNF/8RaMo2b5165BJkNwcp5gzwl1z0gEhs7DYgHcbmUA+vuhUoHRUTV4FgrKECyxwvj1RHt45fgrnBs9x9X+W3g/cJLTLYM0h9RM+55w1wpju1xyYyFFa2mc9mSFULSZkDvE0Yaj3J68xoK7tnbJ4qUsr2oVLl9WHkZjI8HxWZ6yI1xzF5n3FFQI5egLBHtY1nXa7eB6X5zU3CX8diNNvkb8lp9yNo0nvo5Rvw8PDN1s4f6Uxw1d8XPvIra6cNRWMDg4KC9cuPDob/ywWUOvvaYygVaGgHK5dReF+eM/VtcteQQ//rFqcRkOw/PPq3NSqsHl4EHo7SX7za+TvHGZRUvgCkVInHiGUKQJpqbUQBQOw/CwenZnpzI0mYwybrdv3/kc8/MQDMLHP64G2G98445XUCopI7jWAvlan7XejY3W1rsGwZUlN7yml/LsFOaFNzjYeVrJvOQh3W8fRz3UxeXLyoPxeJQ+0mk4fvy++l16riOZovDD71H1uqlYTjrNBGaxTO9HXiHS+pDrPQ9i5URgpef3iPanPG6s/r0uLWrrNYLdgxDiopRycPX5TXkEQogo8BdAFzAKfEJKmVx1zUngvwNBVIn735dS/kX9Z38KPAek65d/Wkr51mZk2hQ9PQ+3MLzRtNTVrPY+qlUVv+/uvnNNsagG53AYenoI/LN/SQDuHpDTaWVAnE41AGUyara/NAi98YYyCktrE6mUMj6WRbotwa3MbSonGkn8ZIx4UuIJhO+fJbU63DU1pe5/+vQ93tBQ/u6SG1ZDM6XB04zdmiFUM9fvALe0xpLNKv0s6SMQeKB+l0t9tDWQ+4CHzKU3kDOTzMXdvO8jv/pojQBs+f6Uxw1d8mHvstnQ0OeAb0spvyCE+Fz9+HdXXZMHfllKOSSEaAEuCiHOSilT9Z//aynlVzYpx86w0R3Dq+npUQPuuXMqHt/VpQbsUkm9H9SgvSITZZmVxmdsTM2Wq1XlVUQiagY9MaGydW7fVhlBbW13BvdajcK573D11kVcgQjurl5m29oYzybpbz1F6H6GcPWC+8yMMgJLaygrBsFU7N6SG65EE+N+k6O9Lz1YN3DH6DidKnQlhPre0/NA/a4s9eFrbMfX2I6UkvnC/KM3AvDuJwK7mMd5UVvz7tmsIfhZ4Pn66y8D32OVIZBSXlvxekIIMQMkgNQmn73zbCQt9X6s9j4WFuDCBbiqUifXykQB7jY+SzPmq1fV+Y4OFR7K59WaQKWi7nvy5J33p9NMdUQJLOTx3UohKhW8pknO72LomOQen3ElKxfcl8I3K6kPguHWTeadLz3H6VT6jcfhyJE7ns999LvtMex3OxHQaB4zHJt8f6OUcrL+egpofNDFQoinABcwsuL07wshfiyE+EMhhPWA9/6GEOKCEOLC7OzsJsV+RCwNWJalZoFLseuHiQ8vxZlff10N6r/0S/BbvwUf+tDa9+nrU4NhLgd+v5qVJpPKCPj9avYfDCrvYiksVKnU8/mTkEwyebwby3AhkKh6dRLLcJEpZTYu99IguJIV/ZpzlRyT2UnennmbTClzV975QmFh/ftHo0oHn/2sMorV6rr6XXpurpxDSrmckbRlVTJX/i5SKfU7/P73VYhuYQOfUaN5TFh3sVgI8S2gaY0ffR74spQyvOLapJQyssa1CCGaUR7Dp6SU51ecm0IZhy8CI1LK31tP6C1bLN5u3u1i41Kmyq1bKvxTq6n4+VII5ehRNXu2LDWbXrUAfvX813DcuIlHGNg+L+XWJgq1Ck6Pj6Mf+9VHIvtCYYGvXv0qC4WFe/LOLae1ZeGFBzYI2pIH1j25116DREIZB9Pc04vGmt3Lu14sllK+8ICbTgshmqWUk/VBfeY+1wWBrwGfXzIC9XsveRMlIcSfAL+znjx7ine72LjkiSzl7q8ciFaHUKLRe0JQHT+ZZMRVoeZ14y5XcFy6QulAEwfca9n7dWR4QL/mncg73/YYdjSqPLnnnrs7RAR7dtFYs/fY7BrBq8CngC/Uv//d6guEEC7gb4E/W70ovMKICOBl4NIm5dldPIrFxmgUXnwRBgfv5LP7/ffPyhkawtfaSbddZKqaJitzeL1+BjIWgcMdDyf/Opv09k3e+T5cNNbsLTZrCL4A/KUQ4jPATeATAEKIQeA3pZS/Vj/3fiAmhPh0/X1LaaJ/LoRIoHqtvwX85ibl2V2EwyoNc35eLfwGAmpASSQe/l7r7ZxeIpWCvj4Cly8T8LaoReVCQWUBrc5Q2iT7ptWgXjTW7HL0hrKdZGQEvvQllfLpcKhYv8sFv/7r9V3BW7BjdWkfQrWq0k+zWRVK6ulRi7OPmG2P2e8E+2xjmWb3siUbyjSb5Pp1NSMX4k5WjGXBW2+pAWUrOqotpbz6fGo9YWnQeuKJR/KRVrMv8s43UtRQo3mM0YZgJ7lyRZV89q5o+pHPww9/qEpAbMWOVT1obQ0bDc1pNI8h2hA8jhQKdxsHeLSLj3rQ0mg0K9jshjLNZujvV2sAhYIqqlYoqOOBgftu1tJoNJpHjTYEO8ngoKozZNvKANi2Ov7oR+/sWJXyzutHnNWj0Wg0oENDO0s0Ci+8sHZ20FLWkI7jazSaLUYbgp3mfvF6HcfXaDTbhA4NaTQazT5HGwKNRqPZ52hDoNFoNPscbQg0Go1mn6MNgUaj0exztCHQaDSafY42BBqNRrPP0YZAo9Fo9jnaEGg0Gs0+RxsCjUaj2edoQ6DRaDT7HF1r6HFjYWFrWlRqNBrNfdAewePEUu/bUkm1qCyV1PHCwk5LptFo9jDaEDxODA2ptpQ+n+pjvPR6aGinJdNoNHuYTRkCIURUCPH3Qoih+vfIfa6zhRBv1b9eXXG+WwjxmhBiWAjxF0II12bk2fWkUmu3qEyldkIajUazT9isR/A54NtSyj7g2/XjtShIKU/Wvz624vy/A/5QStkLJIHPbFKe3U04rFtUajSabWezhuBngS/XX38ZeHmjbxRCCOCDwFfezfv3JH19ukWlRqPZdjZrCBqllJP111NA432ucwshLgghzgshXq6fiwEpKWW1fnwbaL3fg4QQv1G/x4XZ2dlNiv2YstSVzLJUi0rLUsc6a0ij0Wwh66aPCiG+BTSt8aPPrzyQUkohhLzPbTqllONCiAPAd4QQbwPphxFUSvlF4IsAg4OD93vO7ke3qNRoNNvMuoZASvnC/X4mhJgWQjRLKSeFEM3AzH3uMV7/fl0I8T3gFPDXQFgI4ax7BW3A+Lv4DBqNRqPZBJsNDb0KfKr++lPA362+QAgREUJY9ddx4FngipRSAt8Ffv5B79doNBrN1rJZQ/AF4ENCiCHghfoxQohBIcT/ql/TD1wQQvwINfB/QUp5pf6z3wV+WwgxjFoz+KNNyqPRaDSah0SoifnuYnBwUF64cGGnxdBoNJpdhRDiopRycPV5vbNYo9Fo9jm70iMQQswCN+uHcWBuB8V53ND6uBetk7vR+riX/aKTTillYvXJXWkIViKEuLCWq7Nf0fq4F62Tu9H6uJf9rhMdGtJoNJp9jjYEGo1Gs8/ZC4bgizstwGOG1se9aJ3cjdbHvexrnez6NQKNRqPRbI694BFoNBqNZhNoQ6DRaDT7nF1nCDbaFa1+bVAIcVsI8V+3U8btZCP6EEKcFEL8PyHEZSHEj4UQ/3QnZN1KhBAfFkK8U+92d0+DJCGEVe+CN1zvite1A2JuKxvQyW8LIa7U/ya+LYTo3Ak5t4v19LHiun8ihJBCiH2TTrrrDAEb74oG8G+B72+LVDvHRvSRB35ZSnkE+DDwX4QQ4e0TcWsRQhjAfwM+AgwAvyiEGFh12WeAZL0b3h+iuuPtWTaokzeBQSnlcVSDqH+/vVJuHxvUB0KIAPBbwGvbK+HOshsNwYa6ogkhnkA1yvnm9oi1Y6yrDynlNSnlUP31BKpc+D27C3cxTwHDUsrrUsoy8H9QelnJSj19Bfipepe8vcq6OpFSfldKudQb9TyqFPxeZSN/I6Amj/8OKG6ncDvNbjQE63ZFE0I4gP8E/M52CrZDbLRLHABCiKcAFzCy1YJtI63A2IrjtbrdLV9T73+RRlW83atsRCcr+QzwjS2VaGdZVx9CiNNAu5Tya9sp2OPAuo1pdoJH0BXts8DXpZS398Kk7xF1iaPePOh/A5+SUtYerZSa3YoQ4pPAIPDcTsuyU9Qnj/8Z+PQOi7IjPJaG4BF0RXsaeJ8Q4rOAH3AJIRallA9aT3hseRRd4oQQQeBrwOellOe3SNSdYhxoX3G8Vre7pWtuCyGcQAiY3x7xdoSN6AQhxAuoCcVzUsrSNsm2E6ynjwBwFPheffLYBLwqhPiYlHLP17zfjaGhdbuiSSlfkVJ2SCm7UOGhP9utRmADbKRLnAv4W5QevrKNsm0XrwN9Qoju+mf9BZReVrJSTz8PfEfu7d2U6+pECHEK+J/Ax6SUa04g9hAP1IeUMi2ljEspu+rjxnmUXva8EYDdaQg20hVtP7ERfXwCeD/waSHEW/Wvkzsi7RZQj/n/C+AscBX4SynlZSHE7wkhPla/7I+AWL0b3m/z4GyzXc8GdfIfUB7zX9X/JlYbzz3DBvWxb9ElJjQajWafsxs9Ao1Go9E8QrQh0Gg0mn2ONgQajUazz9GGQKPRaPY52hBoNBrNPkcbAo1Go9nnaEOg0Wg0+5z/D+EFOoSuo5loAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "enc_h0 = pca.transform(X_h0)\n",
    "enc_h1 = pca.transform(X_corr)\n",
    "\n",
    "plt.scatter(enc_h0[:,0], enc_h0[:,1], alpha=0.2, color='green', label='white wine')\n",
    "plt.scatter(enc_h1[:,0], enc_h1[:,1], alpha=0.2, color='red', label='red wine')\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define our online drift detector. We specify an expected run-time (in the absence of drift) of 50 time-steps, and a window size of 10 time-steps. Upon initialising the detector thresholds will be computed using 2500 boostrap samples. These values of `ert`, `window_size` and `n_bootstraps` are lower than a typical use-case in order to demonstrate the average behaviour of the detector over a large number of runs in a reasonable time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No GPU detected, fall back on CPU.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  3%|▎         | 74/2500 [00:00<00:03, 736.80it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating permutations of kernel matrix..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2500/2500 [00:04<00:00, 599.06it/s]\n",
      "Computing thresholds: 100%|██████████| 10/10 [00:02<00:00,  3.59it/s]\n"
     ]
    }
   ],
   "source": [
    "from alibi_detect.cd import MMDDriftOnline\n",
    "\n",
    "ert = 50\n",
    "window_size = 10\n",
    "\n",
    "cd = MMDDriftOnline(\n",
    "    X_ref, ert, window_size, backend='pytorch', preprocess_fn=pca.transform, n_bootstraps=2500\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define a function which will simulate a single run and return the run-time. Note how the detector acts on single instances at a time, the run-time is considered as the time elapsed after the test-window has been filled, and that the detector is stateful and must be reset between detections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_run(cd, X, window_size):\n",
    "    n = X.shape[0]\n",
    "    perm = np.random.permutation(n)\n",
    "    t = 0\n",
    "    cd.reset_state()\n",
    "    while True:\n",
    "        pred = cd.predict(X[perm[t%n]])\n",
    "        if pred['data']['is_drift'] == 1:\n",
    "            return t\n",
    "        else:\n",
    "            t += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we look at the distribution of run-times when operating on the held-out data from the reference distribution of white wine samples. We report the average run-time, however note that the targeted run-time distribution, a Geometric distribution with mean `ert`, is very high variance so the empirical average may not be that close to `ert` over a relatively small number of runs. We can see that the detector accurately targets the desired Geometric distribution however by inspecting the linearity of a [Q-Q plot](https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average run-time under no-drift: 48.516\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyAklEQVR4nO3dd5hURdbH8e8BRUBcEWRdRBgMyIphETHumnNE9zWD4OKKi3FXDCgmVMxhzYqKoLQBwyq6JkDMEQmKKMoqqIiCIAgSJJz3j6qWZrpnpmeYnu6e+X2ep5++XX3v7bo09KGqbtUxd0dERCRVvXxXQERECo+Cg4iIpFFwEBGRNAoOIiKSRsFBRETSKDiIiEgaBQepk8zMzWyzKh471cz2KeO9Xc1scqZ9zexCM7uvajWuVP32MLNvc/05UrspOEjRiD+0i8xsgZn9YGaDzaxJvuuVyt3fcPf2Zbx3lbv/HcDM2sYAtUZVPsfMTjSz5fHP4mczG29mh1ThPIPN7Mqq1EFqNwUHKTaHunsToBPQGbio9A5V/cEtQu/EP4umwP3AMDNbL79VktpCwUGKkrtPB14AtoLfuolOM7MvgC9i2clmNsXM5pjZcDPbsNRpDjKzL83sRzO73szqxeM2NbNXzGx2fC9hZk1LHbu9mU0ys5/M7AEzaxiPLbNLx8wuM7Oh8eXr8Xlu/N//7rGeW6fs/3szW2hmLSr4s1gBDAIaAZtm+NwtzOxVM5trZp+Y2WGxvBfQFTgv1uHZ8j5H6hYFBylKZtYaOAgYl1J8OLAj0MHM9gKuBo4GWgLTgEdLneYIQuujE9AF6Jk8fTx2Q2ALoDVwWaljuwL7E36MNydDC6YCu8Xnpu7exN1fi/XrlrLPccAod59V3oliS+nvwAJiYEx5b03gWeBl4PfAGUDCzNq7+0AgAVwX63BoJa9BajEFByk2T5vZXOBN4DXgqpT3rnb3Oe6+iPDjPcjdx7r7EuACYGcza5uy/7Vx/6+BfxN+jHH3Ke4+wt2XxB/mm4DdS9Xjdnf/xt3nAAOSx66mIcBxZmbx9QnAQ+Xsv1P8s/g+fv4R7j6v9D5AE+Aad//V3V8Bnqum+kotVlf6ZqX2ONzdR5bx3jcp2xsCY5Mv3H2Bmc0GWgFTM+w/LR6DmW0A3ALsCqxD+E/UT+V81m/Hrg53f8/MFgJ7mNkMYDNgeDmHvOvuf6ngtBsC38Sup6RphD8HkTKp5SC1SeoSw98BJckXZrY20ByYnrJP65TtNvEYCK0RB7Z2998RunqMVZV1bFXqmmpI/LwTgCfcfXElz1vad0Dr5HhK1IaVfw5allkyUnCQ2uoR4G9m1tHM1iL84L/n7lNT9jnXzNaL4xdnAY/F8nUI/ffzzKwVcG6G859mZhuZWTOgX8qx2ZoFrAA2KVU+lDAW0g14sJLnzOQ9YCFh0HlNM9sDOJSV4y8/ZKiDiIKD1E6x6+li4ElgBmHg+NhSuz0DfAiMB/5LuB0UoD9hkHpeLH8qw0c8TBjk/RL4H1CpuQLuvpAwVvFWvItop1j+DaE7zIE3KnPOMj7nV0IwOBD4EbgT6O7un8Vd7icM4M81s6dX9/Ok9jAl+xEpLGY2CPjO3St7B5RItdGAtEgBiXdT/RXYNs9VkTpO3UoiBcLMrgAmAte7+1f5ro/UbepWEhGRNGo5iIhImqIec1h//fW9bdu2+a6GiEhR+fDDD39093LX7Crq4NC2bVvGjBmT72qIiBQVM5tW0T4561Yys4Zm9r6ZTYgrQfaP5YPN7Ku4/vx4M+sYy83Mbo2raH5kZp1yVTcRESlfLlsOS4C94po2awJvmtkL8b1z3f2JUvsfCLSLjx2Bu+KziIjUsJy1HDxYEF+uGR/l3RrVBXgwHvcu0NTMWuaqfiIiUrac3q1kZvXNbDwwExjh7u/FtwbErqOb47o3EFaJTF3p8lu0cqSISF7kNDi4+3J37whsBOxgZlsR1tX/I7A90Aw4vzLnNLNeZjbGzMbMmlVuDhQREamiGpnn4O5zgdHAAe4+I3YdLQEeAHaIu01n1WWQN2LV5ZWT5xro7p3dvXOLFuXeiSUiUuskEtC2LdSrF54Tidx8Ti7vVmqRzLtrZo2AfYHPkuMIMdvV4YTlAiAkNeke71raCZjn7jNyVT8RkWKTSECvXjBtGriH5169chMgcnm3UktgiJnVJwShYe7+XEzc3oKQPGU88I+4//OEnMBTCOvP/y2HdRMRKTr9+sHChauWLVwYyrt2rd7PyllwcPePyLCypLvvVcb+DpyWq/qIiBS7r7+uXPnq0NpKIiJFok2bypWvDgUHEZEiMWAANG68alnjxqG8uik4iIgUia5dYeBAKCkBs/A8cGD1jzdAkS+8JyJS13TtmptgUJpaDiIikkbBQURE0ig4iIhIGgUHERFJo+AgIiJpFBxERCSNgoOIiKRRcBARkTQKDiIikkbBQURE0ig4iIhIGgUHERFJo+AgIiJpFBxERCSNgoOIiKRRcBARkTQ5Cw5m1tDM3jezCWb2iZn1j+Ubm9l7ZjbFzB4zswaxfK34ekp8v22u6iYikkhA27ZQr154TiTyXaPCksuWwxJgL3f/E9AROMDMdgKuBW52982An4CT4v4nAT/F8pvjfiIi1S6RgF69YNo0cA/PvXopQKTKWXDwYEF8uWZ8OLAX8EQsHwIcHre7xNfE9/c2M8tV/USk7urXDxYuXLVs4cJQLkFOxxzMrL6ZjQdmAiOA/wFz3X1Z3OVboFXcbgV8AxDfnwc0z3DOXmY2xszGzJo1K5fVF5Fa6uuvK1deF+U0OLj7cnfvCGwE7AD8sRrOOdDdO7t75xYtWqzu6USkDmrTpnLldVGN3K3k7nOB0cDOQFMzWyO+tREwPW5PB1oDxPfXBWbXRP1EpG4ZMAAaN161rHHjUF4Uli2DpUtz+hG5vFuphZk1jduNgH2BTwlB4si4Ww/gmbg9PL4mvv+Ku3uu6icidVfXrjBwIJSUgFl4HjgwlBc0d/jvf2GbbeD223P6UWtUvEuVtQSGmFl9QhAa5u7Pmdkk4FEzuxIYB9wf978feMjMpgBzgGNzWDcRqeO6di2CYJBqwgTo0wdGjYJ27aB9+5x+XM6Cg7t/BGybofxLwvhD6fLFwFG5qo+ISFH67ju4+GJ44AFYbz245Rb4xz+gQYOcfmwuWw4iIlJVv/wCN9wA110XxhjOPjvca7veejXy8QoOIiKFZPlyePBBuOii0Go46ii45hrYZJMarYbWVhIRKRQjR8J220HPnuG+2rfegmHDajwwgIKDiEj+ffopHHII7LsvzJsHjz4Kb78Nu+yStyopOIiI5MvMmXDqqbD11vDGG2F84dNP4Zhjwj22eaQxBxGRmrZ4Mfz733DVVWFRp9694dJLYf31812z3yg4iIjUlBUrQpfRBReEhZwOOyy0FnI8Z6Eq1K0kIlIT3nwTdtopzLxr3hxeeQWeeaYgAwMoOIiI5NaUKfB//we77hpuTR0yBMaMgT33zHfNyqXgICKSC3PmwL/+BR06wEsvweWXw+efQ/fuIf1cgdOYg4hIdfr1V7jjDrjiinBbas+eITC0bJnvmlVK4YcvEZFi4A5PPhlaCmefDdtvD+PHw733Fl1gAAUHEZHV98EHsNtucOSR0LAhvPBC6Eraeut816zKFBxERKpq2rRw99EOO4TxhHvuCa2FAw7Id81Wm8YcREQqa948uPrqMJHNLKyWev75sM46+a5ZtVFwEBHJ1rJlYQzh0kth1iw44YSQW7R163zXrNqpW0lEpCKp6TlPPTUMOo8ZE5bWroWBARQcRETKN2FCWC31kENCroWnn4bRo8PS2rWYgoOISCbffRfmKGy7bRhkvvVWmDgRunTJ+4qpNUFjDiIiqX75Ba6/PjyWLYM+fcKAc9Om+a5ZjcpZy8HMWpvZaDObZGafmNlZsfwyM5tuZuPj46CUYy4wsylmNtnM9s9V3URE0ixfDoMGQbt20L9/6Eb69NMQJOpYYIDcthyWAX3cfayZrQN8aGYj4ns3u/sNqTubWQfgWGBLYENgpJlt7u7Lc1hHEZGQnrNPH/joo7By6hNP5DULWyHIWcvB3We4+9i4PR/4FGhVziFdgEfdfYm7fwVMAXbIVf1ERJg0CQ4+OAw4//xzQaTnLBQ1MiBtZm2BbYH3YtHpZvaRmQ0ys/ViWSvgm5TDviVDMDGzXmY2xszGzJo1K5fVFpHaaubMkH1tm23grbdC11GBpOcsFDkPDmbWBHgS+Ke7/wzcBWwKdARmADdW5nzuPtDdO7t75xYtWlR3dUWkNlu0KMxs3mwzuO++MGdhyhQ455ywJpL8Jqd3K5nZmoTAkHD3pwDc/YeU9+8FnosvpwOps0k2imUiIqtnxQp45BG48MKCT89ZKHJ5t5IB9wOfuvtNKeWpa9ceAUyM28OBY81sLTPbGGgHvJ+r+olIHfHGG2GQuVs3WH/9MIGtgNNzFopcthz+DJwAfGxm42PZhcBxZtYRcGAqcAqAu39iZsOASYQ7nU7TnUoiUmVffBEWw/vPf6BVq5Ces1u3osjCVghyFhzc/U0g08jO8+UcMwAYkKs6iUgdMGdOyMJ2xx3QoEHYPvtsaNw43zUrKpohLSK1Q+n0nCedFNJz/uEP+a5ZUVL7SkSKW+n0nDvsEBbLGzhQgWE1KDiISPF6//2V6TkbNYIXXwyPrbbKd82KXqWCg5nVM7Pf5aoyIiJZmTYNjj8edtwxDDwPHAjjxsH+WpKtulQYHMzsYTP7nZmtTbjtdJKZnZv7qomIlDJvHvTtG25DffppuOiiEBxOPhnW0BBqdcqm5dAhzmw+HHgB2Jhwi6qISM1YtgzuvDPMbL722rDMxeefh8HnWpS3uZBkExzWjDOdDweGu/tSwhwFEZHccofnnoOtt4bTTgtjCR9+GOYsbLRRvmtXq2UTHO4hTFZbG3jdzEqAn3NZKRERxo+HffaBQw8Ny1888wy88gp06pTvmtUJFQYHd7/V3Vu5+0EeTAP2rIG6iUhdNH06/O1vIQhMmAC33RbScx52mFZMrUHZDEhvYGb3m9kL8XUHoEfOayYidcuCBXDppbD55vDwwyH5zpQpcPrpsOaa+a5dnZNNt9Jg4CVCdjaAz4F/5qg+IlLXLF8O998fgsLll4f0nJ99VmfTcxaKbILD+u4+DFgB4O7LAC2IJyKrb8SI0H30979D27YhC9tjj8HGG+e7ZnVeNsHhFzNrTrxDycx2AubltFYiUrsl03Putx/Mnw/DhoWMbDvvnO+aSZTNrJGzCbkWNjWzt4AWwJE5rZWI1E4zZ4ZxhXvvhSZNQtfRGWfAWmvlu2ZSSoXBwd3HmtnuQHvCEtyT41wHEZHsLFoE//53SNG5aFFIz3nJJSH5jhSkbO5W6g4cD2wHdCIk6+me64qJSC2wYgUMHRqWu7jwQthrL/jkE7j11goDQyIRhiHq1QvPiUSN1FiibLqVtk/ZbgjsDYwFHsxJjUSkdnjjjbCE9pgxYdD5wQdhjz2yOjSRgF69YOHC8HratPAaoGvX3FRXVmXulVsJw8yaAo+6+wE5qVEldO7c2ceMGZPvaohIqtT0nBttBFddFX7RK5Ges23bEBBKKymBqVOrraZ1lpl96O6dy9unKvkcfiEsvicistKcOfDPf4akOyNGwJVXwuTJcMIJlc7b/PXXlSuX6ldht5KZPcvKhfbqAR2AYbmslIgUkSVLVqbn/PnnMGehf//VysLWpk3mlkObNqtRT6mUbMYcbkjZXgZMc/dvKzrIzFoTxiU2IASXge5+i5k1Ax4D2hIW9Dva3X8yMwNuAQ4CFgInuvvYSlyLiNSkZHrO88+HL7+EAw4It6ZWQxa2AQNWHXMAaNw4lEvNyGbhvddSHm9lExiiZUAfd+8A7AScFtdl6guMcvd2wKj4GuBAoF189ALuquS1iEhNee892HVXOOqo8Kv90kvwwgvVlp6za9eQ3K2kJKy1V1ISXmswuuaU2XIws/lkzttggLt7uelC3X0GMCNuzzezT4FWQBdgj7jbEOBV4PxY/qCHEfJ3zaypmbWM5xGRQjB1KlxwATz6KGywQfjF7tkT6tev9o/q2lXBIJ/KDA7uXm3plcysLbAt8B6wQcoP/veEbicIgeOblMO+jWWrBAcz60VoWdBGHZAiNWPevHDX0S23hMHliy+Gc89VFrZaLOukq2b2e8I8BwDcPav7BsysCfAk8E93/9lS1mN3dzezSt1L6+4DgYEQbmWtzLEiUklLl4bWwWWXwezZ0L17uAtJWdhqvWxmSB9mZl8AXwGvEQaRX8jm5DG96JNAwt2fisU/mFnL+H5LYGYsnw60Tjl8o1gmIjUtmZ5zm21CPoWttgqT2QYPziowaHZz8cvm5uMrCAPKn7v7xoQZ0u9WdFC8++h+4FN3vynlreGsTBbUA3gmpby7BTsB8zTeIJIH48atVnrO5OzmadNCjEnOblaAKC7ZBIel7j4bqGdm9dx9NFDuzLroz8AJwF5mNj4+DgKuAfaNrZF94muA54EvgSnAvcCplbwWEVkd06fDiSfCdtutVnrOfv1WvQUVwut+/aq3upJb2Yw5zI3jBq8DCTObSZglXS53f5NwZ1Mme2fY34HTsqiPiFSnBQvguuvghhtCVrZzzgmL5FUxC5tmN9cOZbYczOwoM2tIuMV0IfAv4EXgf8ChNVM9EcmZZHrOdu3C7ObDDgvpOa+7brXSc5Z1E6FuLiwu5XUrHQ98DdwNHED4z/0Qd781djOJSLFKTc+58cbwzjth7kI1pOccMCDMi0ul2c3Fp8zg4O5HAJsBI4EzgG/N7O6Y+EdEitEnn8BBB6Wn59xpp2r7CM1urh2yXrI75pE+kjBQ3MzdW1dwSM5pyW6RLP3ww8r0nOusEyaxnX660nPWUdks2Z3VJDgzWw/4K3AM0Ax4YvWrJyI5t2gR3HxzSM+5eHEICJdcAs2b57tmUuDKW1upCXAEcBxh6YvhhDkPr3plMwSJSM1asQIefjjcdfTNN3D44XDttbD55vmumRSJ8gakpwL7A3cCbdz9FHcfrcAgUuBefx123DEk2fn97+HVV0NWtioGBs12rpvKCw6t3b2buz/n7ktrrEYiUjVffAFHHAG77w7ffw8PPQTvvx9el6GiH37Ndq67yrtbaVFNVkREqmj2bDjrrJCec+TIcM/o559Dt27lpufM5odfs53rrqrkkBaRQrBkCdx4I2y2Gdx+e8irMGVKGGdo1Agov2WQzQ+/ZjvXXQoOIsXGHR5/PLQUzjknzFGYMAHuuSck4Ikqahlk88Ov2c51V3l3Kz1L5kxwALj7YTmpkYiU7d13oU8fePtt2HrrkJ5zv/0y7lpey6Br1/ADP21a+nGpP/zK5Vx3lddyuAG4kZDHYRFhpdR7gQWE9ZVEpIY8/e+pPLv2sbDzzsx870vePeleGDeOxKz9yuw2qqhlkM0yF5rtXIe5e7kPYEw2Zfl4bLfddi5Sq82d6xMPOc8XsZb/QiPvz8W+NvO9cWP33r3dGzd2D51G4dG4sfvQoeHQkpJV30s+SkpWnn7o0PDaLDwnj5XaLZvf8AqXzzCzT4GD3f3L+Hpj4Hl33yK3YatiWj5Daq2U9JwrfpzNg3TnIq5kOiuzsNWvHxZWLa2kBKZOXTnmULpLSP/zl+paPuNfwKtm9iUhP0MJcEo11E9ESkum5zz3XJg8Gfbck86jb2Qc26btmikwwMpuo2QA6NcvlLVpE7qMFBgkGxUGB3d/0czaAX+MRZ+5+5LcVkukbkgkVv54H/iHcdxkfWj/3Wg+oz3XtxjOXj0PYc6XBhkGjstqOaQOKHftqmAgVVPhraxm1hg4Fzjd3ScAbczskJzXTKQWyTTfINnts3TadAb5iTw7Yzuaffcxp3E7W/Mxg2YdSq9TjIMOyjxw3KuX8iZIDlU0KAE8BpwHTIyvGwPjKzquJh4akJZiMHRo5oHjNs3me38u9l9o5Itp4Ndwnv+OuRkHkMsaONaAslQF1TQgPcbdO5vZOHffNpZNcPc/5TRqZUED0lKoUruL6tVbtfunHsv5Gw9wBRfTku95hGO5gKuZRtuM5zILi6yKVJdsBqSzmSH9q5k1Ik6IM7NNgQrHHMxskJnNNLOJKWWXmdl0MxsfHwelvHeBmU0xs8lmtn8W9RIpSKVnJqcGhn15mXFsy32czJdswo68y/E8UmZgAM1GlvzIJjhcCrwItDazBDCK0M1UkcGE3NOl3ezuHePjeQAz6wAcC2wZj7nTzOpn8RkiBSfTzOQOfMLzHMjL7M/a/MKRPE6XZm8ysfGO5Z5LYwiSL+UGBzOrBySzwJ0IPAJ0dvdXKzqxu78OzMmyHl2AR919ibt/BUwBdsjyWJGCkjoz+ff8wN2cwkdsw068Sx9uoAOTeKHxkdxyq6XNPu7dW7ORpTCUeyuru68ws/PcfRjw32r6zNPNrDswBujj7j8BrYB3U/b5NpalMbNeQC+ANmpvSwFq0wZ+mLaIs7mJvlxDQxZzO6czoN4l/OjN0+Yb6MdfClE23UojzewcM2ttZs2Sjyp+3l3ApkBHYAZh7aZKcfeB7t7Z3Tu3aNGiitUQyZEVK3j4wIf4wjZnABcxgn3Zkk+4sPEt3Pxgc1asCLOXFRCk0GUzQ/qY+HxaSpkDm1T2w9z9h+S2md0LPBdfTgdap+y6USwTKQqJBAzv8xrn/tCHXfiQaS224+h6CZ6YuRtt2sBAzUyWIpPNDOmNq+vDzKylu8+IL48AkncyDQceNrObgA2BdsD71fW5Irk0/IbPWafv+Ty2/Gm+pjXdeIinFxzPPffWY5gCghSpCoNDnCF9NtDG3XvFpTTau/tzFRz3CLAHsL6ZfUu462kPM+tIaHlMJa7R5O6fmNkwYBKwDDjN3ctYOUakQMyeDf37c+Btd7GYhlzIAG7mXyymESxamTdBpBhlMwnuMeBDoLu7bxWDxdvu3rEG6lcuTYKTvFiyBG67Da68EubP554VJ3MJ/ZnJBqvspslrUqiqaxLcpu5+HbAUwN0XElZnFalb3HnjzMf5uskWcO65jFi4C39e5yP+wd1pgQE0eU2KW85mSIvUFokEHNHyXd6q9xd2ve1o5i5rwr68zH5Ln+fteVtmPEaT16TY5XKGtEjRe/rmr2jQ41j+8/3ObMKXnMR9bMs4RrJvmcdo8prUBtncrTTCzMYCOxG6k85y9x9zXjORfJo7F666igOvv4Xl1Kc/l3A95/ILTco9zCzMYxApdmUGBzPrVKooeQtqGzNr4+5jc1ctkTxZupQP/n4Pmw69jKYr5vAwPbiIK/ku84T9NBpnkNqivJZDcvZyQ6AzMIHQctiGsPTFzrmtmkjNSQx1Rp/9LOfMOo/tmcwo9uIcbmB8hvScZdE4g9QmZY45uPue7r4nocXQKS5ZsR2wLZq9LLVAMjtbJxvLhifsxX2zugBwKMPZh5FlBgaL9+o1bx4eWiRPaqNsls9o7+4fJ1+4+0Qz2yKHdRLJqUQCzjoLGs7+lgH04wQeYjbNOZU7uJeTWcaaGY8zI23RPJHaKpvg8LGZ3QcMja+7Ah/lrkoi1Sc1I1uzZrB4MfDLAs7nWvpwI/VYwXWcx9VcwM+sW+Z5Sko00Cx1SzbB4USgN3BWfP06YXVVkYKWzMiWTLzz0+zl9GQQV3Axf+CHCtNzJmksQeqicoNDzMb2Qhx7uLlmqiRSPVIzsu3HS9zAOWzNRN5iF7rwDO9TdhY2s5Dis6RE3UhSN5U7CS4ufrfCzMpub4sUoEQi5HDekom8wAG8xAE0ZiFH8jh/4c2MgSE50FxSAg89FIKDci9IXZVNt9ICwrjDCOCXZKG7n5mzWolUUXKweY3Z33MPl3AS9/Mzv+NsbuQOTuNX1sp4XPPmcMstCgQiSdkEh6fiQ6SgJRJw1skLOWXRzb+l57yNM7iCi5lD84zHKCiIZJZNcHgM2CxuT3H3xTmsj0iVJB5awcs9EozzC2nNtzzFEZzPtUyh3Sr7NW8Oc+bollSRipS3fMYawFVAT2AaYXZ0azN7AOjn7ktrpooi5Rtx0WtscXUfuvqHfEBnupLgDXZL20+3o4pkr7wB6euBZsDG7r6du3cCNgWaAjfUQN1E0iRnNZtBe/ucp+1w9h2wB+uvmElXhrIj72UMDLodVaRyygsOhwAnu/v8ZIG7/0yY83BQrismUlpy3sKCaT9yC2cykS3Zi1e4gKtoz2Qepiue4a908+Za2kKkssobc3DPkEPU3ZebWfm5RUVyoP+FS+i98DYu4krWYT4D6cVlXJYxCxtA/fowZIiCgkhVlNdymGRm3UsXmlk34LPcVUmkFHcYNowXv96CGziXt/gzW/Mxp3JXmYGhcWMFBpHVUV7L4TTgKTPrCXwYyzoDjYAjKjqxmQ0idE3NdPetYlkzwt1PbYGpwNHu/pOZGXALobtqIXCi8kUIAO+8A336wDvvsGTNbdhn6QhGsU+5h2hWs8jqK2/J7unuviNwOeGHfCpwubvv4O7ZLNk9GDigVFlfYJS7tyOkG+0byw8E2sVHL7R2k3z1FRxzDOyyC0ydyjsn38/u64wtNzA0aABDh2pWs0h1yCZN6CvAK5U9sbu/bmZtSxV3AfaI20OAV4HzY/mDcYzjXTNramYt3X0GUrfMnRv+23/rrWHQ4JJLeKzNufQ8s8lv6yRloslsItUrm0lw1WmDlB/87+G3DuNWwDcp+30by9KCg5n1IrQuaKOcjLXH0qVw993Qv3+YpdajB091upJe/Vsxe3bmQzRvQSR3yl14L5diK6HSdz25+8CYla5zixYtclAzqVHu8Mwz/NxmKzjzTEbN/hPb+lhs8AP835llBwYIORpEJDdqOjj8YGYtAeLzzFg+HWidst9GKBVp7Td2LOy1Fxx+ONO/r8fBPBfTc3bM6nA1HEVyp6aDw3CgR9zuATyTUt7dgp2AeRpvqMW+/RZ69IDOnVn84URO4w624SOe52DCKi0V04xnkdzKWXAws0eAd4D2ZvatmZ0EXAPsa2ZfAPvE1wDPA18CU4B7gVNzVS/Jj0QCtmwznyvsYha23pzFDz7GNX4eG8yfwp2cWmbe5kzq19eMZ5Fcy9mAtLsfV8Zbe2fY1wnzKqTIJfMppI4V1GcZPRnEKC7hD/zAwxzHhVxVYXrOTBo0gEGDFBhEcq2m71aSWuzUU+GuUjNUUtNzvsmfK0zPWR7dripSc/J2t5LUDokErL9+WCU1NTCkpudsxCL+jyfYlTcqHRiaNw8T29zhxx8VGERqiloOUmmZuo6SNuB7Lo/pOeexLv/iJu7gNJbSIOvzN28eAoGI5I+Cg1RKpq4jgEYs5Gxu4nyuZS2WcCtncgUX8xPNKnX+xo1D15GI5Je6lSRriUSYxJzKWMEJPMhk2nMlF/My+9GBSZzNzVkHhnrxb2FJie5CEikUajlI1vr1C33/SbvzKjfSh+0Yy/tsz/E8zJvsWubxGlAWKR5qOUiFkoPO06aF15szmafpwqvsSQtm0ZWh7MS7aYGhd+8QTJIPDSiLFA8FBylTIgFNmkC3bmHwuTk/citnMJGt2JPR9OXqjOk5k3cY3XlnHisvIqtF3UqS0T77wKhRYbsBSziTW+nHAJqw4Lf0nLP4/W/7q8tIpHZRcKjjEgk45RT45ZdM7zpH8TjXcj4bM5XnOJjzuI5P6bDqXsooLlLrKDjUYamtg9J24h1u4mx25l0msA37kDk9Z0lJjispInmhMYc6qqzAsDFf8hhH8w67UMI0enI/ncicnrNBA62MKlJbKTjUMclB5tKBYV3mcj3n8ClbcDD/5TIuZXM+5wF6soL6aedp3lwL4InUZupWqkMyzW5eg6X8g7u5lP40Yw6DOZGLuYLvaJV2fMOGcN99CggidYGCQx2QSEDPnvDrr6mlzmEM5zrOoz2fM5K9OYcbmFBGFra994aRI2uitiJSCNStVMslEtC9+6qBoRMfMpo9eYbDWU59DuY59mVExsCQnLOgwCBSt6jlUIslEnDCCStvNd2IbxhAP7rzEDNpQW/u5F5OZnn8a6DWgYgkqeVQS6TmVUg+unULgaEJ87mCi/iczTmaYVxNX9rxBXfT+7fA0Lu3AoOIrKSWQ5ErawltWJme8/KYnjPB8VzIVXzNqpMTevfWUhcisqq8BAczmwrMB5YDy9y9s5k1Ax4D2gJTgaPd/ad81K8YJBLQowcsX575/f15kRs4h634hDf4C4cxnA/YYZV9dPeRiJQln91Ke7p7R3fvHF/3BUa5eztgVHwtGZx6augyyhQYtuJjXmR/XuRAGrKYv/Iku/F6WmDo3RsWLVJgEJHMCmnMoQswJG4PAQ7PX1UK06mnhsQ4mbqRNuB7BnIy4+nI9nzAP7mZDkziP/wVsFX2VTeSiFQkX2MODrxsZg7c4+4DgQ3cfUZ8/3tggzzVrSBtuSVMmpRe3oiF9OFGzudaGvArt3AWV3JRmVnY9t5bgUFEKpav4PAXd59uZr8HRpjZZ6lvurvHwJHGzHoBvQDatGmT+5oWgH32SQ8MIT3nQwygHxsxnSf5K+dzLf9js4zn0JLaIlIZeelWcvfp8Xkm8B9gB+AHM2sJEJ9nlnHsQHfv7O6dW7RoUVNVzotEAtZaK30dpD0YzRg6M4QT+Y4N2ZXXOZInfwsMTZqEiWvKwiYiVVXjwcHM1jazdZLbwH7ARGA40CPu1gN4pqbrViiSQaFbt1VnNm/OZJ7hMEazF82ZzfEkVknPuffeIRDMn69AICKrJx/dShsA/zGz5Oc/7O4vmtkHwDAzOwmYBhydh7rlXaZ5C835kUvpzz+4m0U0oi9XcwtnsZhGv+2jQWYRqU41Hhzc/UvgTxnKZwN713R9CknpHAtrsZgzuK3c9JwlJSGngloKIlKdNEO6AJReAwmcoxnGNfT9LT3nuVzPZ2zx2zENGsCSJXmprojUAYU0z6HOSSRgjTVWroEEsDNv8za78BjHMo912ZuRHMpzqwQGs5BoR0QkVxQc8iA1KCRnOSfTc77NnylhGn9jENvxIa+U6mlr2BAeekjdSCKSW+pWqmGlB5yb8hP9GMAZ3MYy1uBSLuMGzmEha69y3BprwODBCgoiUjMUHGpI6YXy1mApvbmLS+nPevzEA/yNi7mCGWyYdqzyLIhITVO3Uo6ldyE5XXiaT9iSWzmLsXRiW8bxd+5PCwz16ysLm4jkh1oOOZRIhKCQtB1juIFz2IPXmMQWHMR/eYEDKb0wHkCHDvDJJzVXVxGRVGo55EhyWW0I6Tkf5ATGsD0dmMQ/uItt+IgXOIjSgWGNNUJrQYFBRPJJLYdqkmlmcxPm05drOJubMJyr6cvVXMB8fpd2fP36MGSIBpxFpDAoOKym0rOaIaTnPIn7uZxL2ICZDKUr/RiQlp4zSQPOIlJo1K1URYlEmIxWOjDsz4uMpyP38A8m057teZ8TGJoxMGjAWUQKlYJDJSUSIRtb6kAzrJqecy2W8FeeZHdeYwzbZzxP796wbJm6kUSkMKlbqRIyjSv8gRlcziX0ZBDzWJd/cjN3cipLaZDxHBtuCNOn10BlRURWg4JDlkqPLTTmF/pwI+dxXVbpOUG3p4pI8VC3UhbWW29lYDBW0J0hfM7mXM6lvMgBdGASfbipzMBQUqLbU0WkuKjlUIEGDWDp0rC9B6O5kT50YhzvsQPH8Bhv8ZeMxyn5jogUM7UcSkkud2EWHkuXQns++y09ZzPmcBwPszPvpAWGpk1X5mxWYBCRYqbgkGKffVZdRnt9ZnEbpzORrdid1zifa/gjn/Eox+Gl/uh694affspDpUVEckDdSmROz3kmt9KPAazNL9zDKVzGZfxIi7RjmzZVUBCR2qdOthwSCWjbdmXX0crA4BzDo3zGH7mO83md3diajzmdOzIGhg4dFBhEpHYquOBgZgeY2WQzm2Jmfav7/OutF7qOpk1btXxn3uYdduZRjmMuTdmbkRzGs6uk50ylu49EpDYrqG4lM6sP3AHsC3wLfGBmw919UnWcv1UrmDt31bJN+B/X0JejeILpbMiJPMBDnMAK6mc8x5prwq+/VkdtREQKV6G1HHYAprj7l+7+K/Ao0KW6Tv7ddyu3m/ITN9CHT9mCg3ieS+jP5nzOEE4sMzB06KDAICJ1Q0G1HIBWwDcpr78Fdkzdwcx6Ab0A2rRpU6UPOZDneYgTWI+fGERPLuYKvqdlxn014CwidVGhtRwq5O4D3b2zu3du0SJ9kDgbn7M577Ej2zKOk7mvzMCw994KDCJSNxVacJgOtE55vVEsqxYbxhTN/2MzDuZ5PuJPGffr3TtMZNNS2iJSVxVat9IHQDsz25gQFI4Fjq+uk0+fHgalU8ceQAviiYiUVlDBwd2XmdnpwEtAfWCQu1frz7aWyxYRqVhBBQcAd38eeD7f9RARqcsKbcxBREQKgIKDiIikUXAQEZE0Cg4iIpJGwUFERNKYu+e7DlVmZrOAaRXumG594Mdqrk6+1bZr0vUUvtp2TXXpekrcvdwlJoo6OFSVmY1x9875rkd1qm3XpOspfLXtmnQ9q1K3koiIpFFwEBGRNHU1OAzMdwVyoLZdk66n8NW2a9L1pKiTYw4iIlK+utpyEBGRcig4iIhImjoXHMzsADObbGZTzKxvvutTFWY21cw+NrPxZjYmljUzsxFm9kV8Xi/f9SyPmQ0ys5lmNjGlLOM1WHBr/M4+MrNO+at5ZmVcz2VmNj1+T+PN7KCU9y6I1zPZzPbPT63LZmatzWy0mU0ys0/M7KxYXpTfUTnXU8zfUUMze9/MJsRr6h/LNzaz92LdHzOzBrF8rfh6Sny/bbkf4O515kHIEfE/YBOgATAB6JDvelXhOqYC65cquw7oG7f7Atfmu54VXMNuQCdgYkXXABwEvAAYsBPwXr7rn+X1XAack2HfDvHv3lrAxvHvZP18X0OpOrYEOsXtdYDPY72L8jsq53qK+TsyoEncXhN4L/7ZDwOOjeV3A73j9qnA3XH7WOCx8s5f11oOOwBT3P1Ld/8VeBTokuc6VZcuwJC4PQQ4PH9VqZi7vw7MKVVc1jV0AR704F2gqZllTvydJ2VcT1m6AI+6+xJ3/wqYQvi7WTDcfYa7j43b84FPgVYU6XdUzvWUpRi+I3f3BfHlmvHhwF7AE7G89HeU/O6eAPY2Myvr/HUtOLQCvkl5/S3l/wUpVA68bGYfmlmvWLaBu8+I298DG+SnaqulrGso5u/t9NjNMiilq6+orid2P2xL+J9p0X9Hpa4Hivg7MrP6ZjYemAmMILRw5rr7srhLar1/u6b4/jygeVnnrmvBobb4i7t3Ag4ETjOz3VLf9NBuLOp7lGvDNQB3AZsCHYEZwI15rU0VmFkT4Engn+7+c+p7xfgdZbieov6O3H25u3cENiK0bP5YXeeua8FhOtA65fVGsayouPv0+DwT+A/hL8UPyWZ8fJ6ZvxpWWVnXUJTfm7v/EP/xrgDuZWW3RFFcj5mtSfghTbj7U7G4aL+jTNdT7N9RkrvPBUYDOxO69JIpoFPr/ds1xffXBWaXdc66Fhw+ANrF0fwGhEGZ4XmuU6WY2dpmtk5yG9gPmEi4jh5xtx7AM/mp4Wop6xqGA93jHTE7AfNSujYKVqk+9yMI3xOE6zk23j2yMdAOeL+m61ee2Bd9P/Cpu9+U8lZRfkdlXU+Rf0ctzKxp3G4E7EsYSxkNHBl3K/0dJb+7I4FXYusvs3yPuNf0g3BXxeeEvrl++a5PFeq/CeEuignAJ8lrIPQdjgK+AEYCzfJd1wqu4xFCM34poV/0pLKugXBXxh3xO/sY6Jzv+md5PQ/F+n4U/2G2TNm/X7yeycCB+a5/huv5C6HL6CNgfHwcVKzfUTnXU8zf0TbAuFj3icAlsXwTQiCbAjwOrBXLG8bXU+L7m5R3fi2fISIiaepat5KIiGRBwUFERNIoOIiISBoFBxERSaPgICIiaRQcpCCYWfOUlTG/T1kpc66ZTarhuhxuZh1SXl9uZvtU4TxtLWWV1ppmZheWev12fM5rvaQ4KDhIQXD32e7e0cNSAHcDN8ftjsCK6v68lBmkmRxOWJUzWbdL3H1kddehBqwSHNx9l3xVRIqPgoMUg/pmdm9cs/7lOBsUM9vUzF6MCxC+YWZ/jOVtzeyVuJjaKDNrE8sHm9ndZvYecF2m481sF+Aw4PrYctk0HndkPMf2ZvZ2XEP/fTNbJ37eG2Y2Nj7K/RGOs4hvt5AnYKSZPZ9y/qlmtn7c7mxmr8btHczsHTMbFz+/fSw/0cyeitfxhZldF8uvARrFa0jEsgUZ6lLfzK43sw/in9cpsbylmb0ej59oZruu5ncoxSbfs/z00KP0g5Q19oG2wDKgY3w9DOgWt0cB7eL2joTlAACeBXrE7Z7A03F7MPAccV3+co4fDByZUp/BhOUGGgBfAtvH8t8BawCNgYaxrB0wJqXuEzNc318JK2jWBzYE5iY/j5RcHUBn4NXUz4rb+wBPxu0TY53WJcyAnQa0ju8tKPW5C0rXC+gFXBS31wLGEPIX9GHl7Pv6wDr5/nuhR80+ymtaixSKr9x9fNz+EGhrYXXNXYDHbeWS9GvF550JP8AQlke4LuVcj7v78gqOL0t7YIa7fwDgcZXSuMbV7WbWEVgObF7BeXYDHnH35cB3ZvZKBftD+PEfYmbtCMtArJny3ih3nxfrMgkoYdXlpsuzH7BNsuUSP6cdYR2yQRYWq3s65c9f6ggFBykGS1K2lwONCF2icz2MS1TGL/G5qsdn8i/gB+BP8byLV+Ncy1jZ3dswpfwKYLS7H2EhH8GrKe+V/vOpzL9rA85w95fS3ghLwR8MDDazm9z9wUqcV4qcxhykKMX/tX9lZkfBb/34f4pvv01YcRegK/BGJY+fT0glWdpkoKWZbR+PWcdWLn08w8OyzycQumHK8zpwTOzvbwnsmfLeVGC7uP1/KeXrsnLp5RMrOH/S0vg///K8BPRO7mdmm1tY+bcE+MHd7wXuI6RAlTpEwUGKWVfgJDNLrlCbTPl6BvA3M/uI8GN9ViWPfxQ4Nw7+bprc2UNq2WOA2+IxIwj/u78T6BHL/sjK1klZ/kNY1XQS8CDwTsp7/YFbzGwMoRWQdB1wtZmNI/uWwUDgo+SAdBnui/UYG29vvSeefw9gQvy8Y4BbsvxMqSW0KqtInpnZYOA5d3+ion1FaopaDiIikkYtBxERSaOWg4iIpFFwEBGRNAoOIiKSRsFBRETSKDiIiEia/weti2O1G4md5AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_runs = 250\n",
    "times_h0 = [time_run(cd, X_h0, window_size) for _ in range(n_runs)]\n",
    "print(f\"Average run-time under no-drift: {np.mean(times_h0)}\")\n",
    "_ = scipy.stats.probplot(np.array(times_h0), dist=scipy.stats.geom, sparams=1/ert, plot=plt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run the detector in an identical manner but on data from the drifted distribution of red wine samples the average run-time is much lower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average run-time under drift: 6.004\n"
     ]
    }
   ],
   "source": [
    "n_runs = 250\n",
    "times_h1 = [time_run(cd, X_corr, window_size) for _ in range(n_runs)]\n",
    "print(f\"Average run-time under drift: {np.mean(times_h1)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Online detection with LSDD and TensorFlow\n",
    "\n",
    "Here we address the same problem but using the least squares density difference (LSDD) as the two-sample distance in a manner similar to [Bu et al. (2017)](https://ieeexplore.ieee.org/abstract/document/7890493). The LSDD between two distributions $p$ and $q$ on $\\mathcal{X}$ is defined as $$LSDD(p,q) = \\int_{\\mathcal{X}} (p(x)-q(x))^2 \\,dx$$ and also has an empirical estimate $\\widehat{LSDD}(\\{X_i\\}_{i=1}^N, \\{Y_i\\}_{i=t}^{t+W})$ that can be updated at low cost as the test window is updated to $\\{Y_i\\}_{i=t+1}^{t+1+W}$.\n",
    "\n",
    "We additionally show that TensorFlow can also be used as the backend and that sometimes it is not necessary to perform preprocessing, making definition of the drift detector simpler. Moreover, in the absence of a learned preprocessing stage we may use all of the reference data available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ref = np.concatenate([X_train, X_ref], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we define the LSDD-based online drift detector, again with an `ert` of 50 and `window_size` of 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Computing thresholds: 100%|██████████| 9/9 [00:09<00:00,  1.00s/it]\n"
     ]
    }
   ],
   "source": [
    "from alibi_detect.cd import LSDDDriftOnline\n",
    "\n",
    "cd = LSDDDriftOnline(\n",
    "    X_ref, ert, window_size, backend='tensorflow', n_bootstraps=2500,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run this new detector on the held out reference data and again see that in the absence of drift the distribution of run-times follows a Geometric distribution with mean `ert`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average run-time under no-drift: 51.648\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyGklEQVR4nO3debyOdf7H8deHRKKkZEq2ZIgW6VTaN9VkmlK/tonStFhTaZkWM+01lUk7UolyKoYaS8uE0h5hkHVSCAkpSiScz++P73W4neU+9znOfe77Puf9fDzux7mv731d1/25zs39Od/l+n7N3REREYlVKdUBiIhI+lFyEBGRfJQcREQkHyUHERHJR8lBRETyUXIQEZF8lBykQjIzN7MDSnjsIjNrW8hrx5vZ/IL2NbPbzey5kkVcrPhOMrOlyX4fKd+UHCRjRF+0G8xsnZmtMLPBZlYj1XHFcvcP3b1ZIa894O5XAZhZoyhB7VSS9zGzy81sS/S7+MnMppvZWSU4z2Azu68kMUj5puQgmeZP7l4DaA1kAX/Lu0NJv3Az0KfR76IW8Dww3Mz2SG1IUl4oOUhGcvdlwFvAQbC1maiHmX0JfBmVXW1mC8zsBzMbbWb75jlNOzP72sy+N7M+ZlYpOq6Jmb1rZquj17LNrFaeY48wszlm9qOZvWBm1aJjC23SMbO7zGxotPlB9HNN9Nf/iVGcB8fsv7eZrTezOkX8LnKAQcAuQJMC3vdAM5toZmvMbLaZnR2VdwY6AH+NYhgT732kYlFykIxkZvWBdsB/Y4rbA0cBLczsFOAfwIXAPsBi4NU8pzmXUPtoDZwDXJF7+ujYfYEDgfrAXXmO7QCcQfgy/j0F1GCKcEL0s5a713D396P4Osbs82dggruvineiqKZ0FbCOKDHGvFYFGAO8A+wN9ASyzayZuw8EsoGHoxj+VMxrkHJMyUEyzb/NbA3wEfA+8EDMa/9w9x/cfQPhy3uQu09z943AbcDRZtYoZv+Hov2/AR4jfBnj7gvcfZy7b4y+mPsCJ+aJ4yl3X+LuPwD35x67g4YAfzYzi7YvBV6Ks3+b6HfxXfT+57r72rz7ADWAB939N3d/FxhbSvFKOVZR2mal/Gjv7uMLeW1JzPN9gWm5G+6+zsxWA/WARQXsvzg6BjOrCzwOHA/UJPwR9WOc99p67I5w90lmth44ycyWAwcAo+Mc8pm7H1fEafcFlkRNT7kWE34PIoVSzUHKk9gphr8FGuZumNmuwJ7Asph96sc8bxAdA6E24sDB7r4boanH2F5hx5Yk1lhDove7FBjh7r8W87x5fQvUz+1PiTRg2+9B0zJLgZQcpLx6BfiLmbUys6qEL/xJ7r4oZp+bzWyPqP/iOmBYVF6T0H6/1szqATcXcP4eZrafmdUGesccm6hVQA6wf57yoYS+kI7Ai8U8Z0EmAesJnc5VzOwk4E9s639ZUUAMIkoOUj5FTU9/B0YCywkdxxfn2W0UMBWYDrxBGA4KcDehk3ptVP5aAW/xMqGT92vgK6BY9wq4+3pCX8XH0SiiNlH5EkJzmAMfFuechbzPb4RkcCbwPdAPuMzd50W7PE/owF9jZv/e0feT8sO02I9IejGzQcC37l7cEVAipUYd0iJpJBpNdR5wWIpDkQpOzUoiacLM7gVmAX3cfWGq45GKTc1KIiKSj2oOIiKST0b3Oey1117eqFGjVIchIpJRpk6d+r27x52zK6OTQ6NGjZgyZUqqwxARyShmtriofdSsJCIi+Sg5iIhIPkoOIiKSj5KDiIjko+QgIiL5KDmIiGSQ7Gxo1AgqVQo/s7OT8z4ZPZRVRKQiyc6Gzp1h/fqwvXhx2Abo0KF030s1BxGRDNG797bEkGv9+lBe2pQcREQyxDffFK98Ryg5iIhkiAYNile+I5QcREQyxP33Q/Xq25dVrx7KS5uSg4hIhujQAQYOhIYNwSz8HDiw9DujQaOVREQySocOyUkGeanmICIi+Sg5iIhIPkoOIiKSj5KDiEgm2bAB/v53+OijpL6NOqRFRDLFO+9A9+7w1VdhuNJxxyXtrVRzEBFJd999B3/+M5xxRphxb9w4uOeepL6lkoOISLrasgX69YPmzeG11+DOO2HmTGjbNulvrWYlEZF0NH06dOkCkyfDKaeEJNGsWZm9vWoOIiLpZN06uPFGyMqChQvhpZdg/PgyTQygmoOISPr497+hZ09YuhSuvhoefBBq105JKKo5iIik2jffwDnnwLnnQq1a8PHHYdKkFCUGUHIQEUmdTZvgn/+EAw8MI5AeegimTYNjjkl1ZMlLDmZWzcwmm9kMM5ttZndH5Y3NbJKZLTCzYWa2c1ReNdpeEL3eKFmxiYik3GefhX6Fm28OHc5z5sBf/wpVqqQ6MiC5NYeNwCnufijQCviDmbUBHgIedfcDgB+BK6P9rwR+jMofjfYTESlffvwRunULtYPVq2HkSBg9Gho1SnVk20lacvBgXbRZJXo4cAowIiofArSPnp8TbRO9fqqZWbLiExEpU+7w8svhnoWBA+Haa2HuXDjvvHC3c5pJap+DmVU2s+nASmAc8BWwxt03R7ssBepFz+sBSwCi19cCexZwzs5mNsXMpqxatSqZ4YuIlI4vv4TTTw8LMTRoAJ9/Do89BjVrpjqyQiU1Obj7FndvBewHHAk0L4VzDnT3LHfPqlOnzo6eTkQkeTZuDNNcHHwwTJoETz4Z+hpat051ZEUqk/sc3H2Nmb0HHA3UMrOdotrBfsCyaLdlQH1gqZntBOwOrC6L+ERESt1774W+hfnz4cIL4dFHYd99Ux1VwpI5WqmOmdWKnu8CnAbMBd4Dzo926wSMip6PjraJXn/X3T1Z8YmIJMWqVdCpUxiB9Ntv8NZbMGxYRiUGSG7NYR9giJlVJiSh4e4+1szmAK+a2X3Af4Hno/2fB14yswXAD8DFSYxNRKR05eTAoEFhOOq6dXD77dC7N1SvnurISiRpycHdZwKHFVD+NaH/IW/5r8AFyYpHRCRpZs2Crl3Dnc3HHw8DBkCLFqmOaofoDmkRkZJavx5uvRUOOwzmzQs1h/ffz/jEAJp4T0SkZN54A665BhYtgssvhz59YK+9Uh1VqVHNQUSkOJYtg/PPh7POgmrVYOJEeOGFcpUYQMlBRCQxW7bAE0+ESfLeeAPuuy8syHPiiamOLCnUrCQiUpQpU0KH89Sp4U7nfv2gSZNUR5VUqjmIiBTmp5/CHEhHHRWak159Fd5+u9wnBlDNQUQkP3cYMQKuuw6++y7c6Xz//WEhngpCyUFEJNbChdCjR7izuVWrsHTnkfluzSr31KwkIgJhqosHH4SWLeGDD6Bv3zB7agVMDKCag4gIfPRR6HCePTus4/z441C/fqqjSinVHESk4lq9Gq66Kkx58fPPYUW2116r8IkBlBxEpCJyhyFDwqpsgwfDTTeFNZz/9KdUR5Y21KwkIhXLvHlh9NHEidCmDTzzDBxySKqjSjuqOYhIxbBhA9xxR0gE06eHmVM//liJoRCqOYhI+TduXKgtfPVVWMf5kUegbt1UR5XWVHMQkfLru+/gkkvClBeVKoUkMXSoEkMClBxEpPzJyYH+/UOH88iRcOedMHMmtG2b6sgyhpqVRKR8mT493LMwaVJYx7lfP2jWLNVRZRzVHESkfFi3Dm68EbKy4Ouv4aWXYPx4JYYSUs1BRDLfqFHQsycsWQJXXx2mwahdO9VRZbSk1RzMrL6ZvWdmc8xstpldF5XfZWbLzGx69GgXc8xtZrbAzOab2RnJik1EyolvvoH27cNj993DNBgDByoxlIJk1hw2Aze6+zQzqwlMNbNx0WuPuvs/Y3c2sxbAxUBLYF9gvJn93t23JDFGEclEmzeH+Y/uvDN0Pj/0EPTqBVWqpDqyciNpycHdlwPLo+c/m9lcoF6cQ84BXnX3jcBCM1sAHAl8mqwYRSQDffZZ6HCeMQP++Ed46ilo1CjVUZU7ZdIhbWaNgMOASVHRNWY208wGmdkeUVk9YEnMYUspIJmYWWczm2JmU1atWpXMsEUknaxZE25kO+YY+P77MER1zBglhiRJenIwsxrASOB6d/8J6A80AVoRahaPFOd87j7Q3bPcPatOnTqlHa6IpBt3eOWVcM/CwIFh2c65c+G888CsxKfNzg55pVKl8DM7u9QiLheSOlrJzKoQEkO2u78G4O4rYl5/FhgbbS4DYufJ3S8qE5GKasEC6N493NmclQVvvgmtW+/wabOzoXNnWL8+bC9eHLYhzK4hyR2tZMDzwFx37xtTvk/MbucCs6Lno4GLzayqmTUGmgKTkxWfiKSxjRvh3nvhoINCH8OTT4afpZAYAHr33pYYcq1fH8olSGbN4VjgUuALM5seld0O/NnMWgEOLAK6ALj7bDMbDswhjHTqoZFKIhXQxImhw3n+fLjwQnj0Udh331J9i2++KV55RZTM0UofAQU1CL4Z55j7gfuTFZOIpLFVq8KiOy++CI0bw1tvwR/+kJS3atAgNCUVVC6Bps8QkdTKyYHnnw8dzi+/DLfdBrNmJS0xANx/P1Svvn1Z9eqhXAIlBxFJndmz4cQTwzrOLVuGSfMeeCD/N3cp69AhDHxq2DAMeGrYMGyrM3obJQcRKXvr14caQqtWYe3m558PfQ0tW27dJdlDTTt0gEWLQsVl0SIlhryUHEQkrpJ8Scc95s03QxJ48EHo2DF0PF9xRdg55vjOnUO/gPu2oaa6F6HsmLunOoYSy8rK8ilTpqQ6DJFyK+/9ABBafOI1wRR2zEsPLuO8D66HESNC/8KAAaFJqQCNGhXcYdywYfgrX3aMmU1196x4+6jmIFKBFVUrKMn9AHmPqcQWrlz/BKdddyCMHQv33Rf6FgpJDKChpulA6zmIVFCJ3CVcki/p2NdaM5Vn6EIWU/mPn84Zs/pBkyZFxqahpqlXrJqDmVUys92SFYyIlJ1EagWFfRnH+5Ju0ABq8hOPcR2TOZJ6LOMiXqVLg7cTSgygoabpoMjkYGYvm9luZrYrYaqLOWZ2c/JDE5FkSqRWUOwvaXdePGcE8+xAevIkA+jKgcxlbPWLuP+BxCfJ01DT1Euk5tAimk21PfAW0JgwLYaIZLBEagXF+pJeuBDOOosTnriAag325tzffUZPe5paDWuV6ItdQ01TK5E+hyrR7KrtgafcfZOZZe4QJxEBwl//BY0qylsr6NChiC/mTZvgkUfgnntCz3bfvtTu2ZNRO6lLM5MlUnN4hjBB3q7AB2bWEPgpmUGJSPKVStPNRx/BYYeFG9rOOCOss9CrFygxZLwik4O7P+Hu9dy9nQeLgZPLIDYR2UFFDVUtcdPN6tVhyovjj4eff4ZRo+D116F+/aKPlYyQSId0XTN73szeirZbAJ2SHpmI7JCk3GXsHmZNbd4cBg8Os6jOng1nn11aYUuaSKRZaTDwHyB3QvX/AdcnKR4RKSWlvqDNvHlwyinQqRMccABMmwZ9+kCNGjscq6SfRJLDXu4+HMgBcPfNgBbhEUlzpXaX8a+/wh13wKGHhjubBwyAjz+GQw7Z0RAljSWSHH4xsz0JK7dhZm2AtUmNSkR2WEluYMtn3Dg4+OCwZOcFF4TaQ5cu202SJ+VTIp/wDYT1nZuY2cfAi0DPpEYlIgmJ1+G8Q3cZf/cdXHIJnH56GMo0bhwMHQp165Zi9JLOihxv5u7TzOxEoBlh2c/57r4p6ZGJSFxFzY2UO/Kod+/QlNSgQUgMcUck5eSE8ay33gobNoTmpNtug2rVknotkn6KnLLbzC4rqNzdX0xKRMWgKbulIiv1aa1nzAhNRpMmwcknQ//+0KzZDkYp6SiRKbsTuVPliJjn1YBTgWmE5iURSZFS63Betw7uugseewxq1w5DVTt2DM1JUmEl0qy0Xf+CmdUCXi3qODOrT0ggdQmd2QPd/XEzqw0MAxoR7ry+0N1/NDMDHgfaAeuBy919WnEuRqQiKZVprUeNgp49YckSuPrqsDpb7dqlFqNkrpIMOfiFMPleUTYDN7p7C6AN0CO6ge5WYIK7NwUmRNsAZwJNo0dnoH8JYhOpMHaow3nJEmjfPjx23z1MgzFwoBKDbFVkzcHMxhANYyUkkxbA8KKOc/flwPLo+c9mNheoB5wDnBTtNgSYCNwSlb/ooRPkMzOrZWb7ROcRkTxK1OG8eTM88UToaM7JgYceCnMhValSJjFL5kikz+GfMc83A4vdfWlx3sTMGgGHAZOAujFf+N8Rmp0gJI4lMYctjcq2Sw5m1plQs6CBloWSCq7IGVNjTZoUOpxnzIA//hGeeir0aosUIJGJ996PeXxcgsRQAxgJXB+tCxF7bmdbrSQh7j7Q3bPcPatOnTrFOVQkYxU1gV5ca9ZA9+5w9NHw/fcwciSMGaPEIHEVWnMws58p+IvbCN/rRS4XGq0DMRLIdvfXouIVuc1FZrYPsDIqXwbETum4X1QmUqElstZzgdzh1VdDs9GqVXDtteFO55o1kx6zZL5Caw7uXtPddyvgUTPBxGDA88Bcd+8b89Jots3q2gkYFVN+mQVtgLXqbxAp4QR6CxaE9RUuuSRMo/3552GoqhKDJCjhFTnMbG/CfQ4AuHtRo6mPJSwn+oWZTY/KbgceBIab2ZXAYuDC6LU3CcNYFxCGsv4l0dhEyrNi3c+wcSM8/HDomd55Z3jySejWDSpXTmqMUv4kMlrpbOARwpTdK4GGwFygZbzj3P0jQhNUQU4tYH8HehQVj0hFk/D9DBMnQteuMH9+mCTvscdg333zHyiSgETuc7iXcJ/C/9y9MeGL/bOkRiUiWxV5P8OqVWGNhZNPht9+gzffhOHDlRhkhySSHDa5+2qgkplVcvf3gLhzcohI6Sl0rec/58Dzz4dV2V5+OUyQN2sWnHlmqkOWciCRPoc10XDUD4BsM1tJuEtaRMpIvvsZZs+GE7uGO5uPPz5MktcybkuvSLEUWnMwswvMrBrhzuX1QC/gbeAr4E9lE56IbGf9+lBDaNUK5swJNYeJE5UYpNTFqzlcAjxNWD/6FeA/7j6kTKISkfzeegt69ICFC+Hyy8P6zXvtleqopJyKd5/DucABwHjCym9LzWxAtPCPiJSVb78No4/atYOqVUNN4YUXlBgkqeJ2SLv7T+4+xN3PBA4C/gs8YWZL4h0nIonJOy1G9+7btvdvuIXPOz0ZOpzHjAl3N0+fDifq7zNJvoRugjOzPYDzgIuA2sCIZAYlUhEUNC1G/2ii+tZM5ZlvupD14lS+Pfh09n3taTjggNQFKxVOvA7pGmZ2qZm9CcwhDF+9F2jg7r3KKkCR8qqgaTFq8hOPcR2TOZL9WMrFvMIxa99WYpAyF6/msIgwOqkfoTN6U5lEJFJBbD/9hXMer/EE17IPy+lPN3pzP2uphRpxJRXi9TnUd/eO7j5WiUGkeOL1JeROuZ07/UVDFjGWsxjJ+axkb47mU67hadZSCyjmsp8ipaTQmoO7byjLQETKi3h9CbnbnTvDXzpuYo8X+nLbprvJoRK96MuT9GRLzH/LhJf9FCllJVlDWkTyiK0pdOqUvy8hr1brP6bnC4dx76Zb+XCXM2jBXF5v2IvO3XbKP01Goiu9iZSihKfsFpGC5a0pbNlS+L578AMPcQtX8xyLNzWAUaM44+yzKWr+e5GyFm8luDHEWcLT3c9OSkQiGaagUUf5OZfyEo9wI3vwI324iRfq38mcs2uURYgixRav5vDP6Od5wO+AodH2n4EVyQxKJJMUthhPrt8zn/504xTe41Pa0JUBLKh+KAP/UTbxiZREvOkz3nf394Fj3f0idx8TPS4Bji+7EEXSU24/gxdSv65e6Vfu4Q5mcghtdv4vt9cewHF8zNqGh6ovQdJeIn0Ou5rZ/u7+NYCZNQZ2TW5YIuktbz9DXn+sOp6Xa3VjtxULwjrOffvyQN26PFC2YYqUWCLJoRcw0cy+Jiz72RDoktSoRNJcYf0Me7OCZ6rfQPv1L0PNA2DoOGjbtuwDFNlBRSYHd3/bzJoCzaOiee6+MblhiaS3vP0MRg6dGciD3EqtzRvgjjvCugvVqqUmQJEdVOR9DmZWHbgZuMbdZwANzOysBI4bZGYrzWxWTNldZrbMzKZHj3Yxr91mZgvMbL6ZnVHC6xFJqoL6GQ5hBh9zLAPoxpyqrWHmTLj7biUGyWiJ3AT3AvAbcHS0vQy4L4HjBgN/KKD8UXdvFT3eBDCzFsDFQMvomH5mVjmB9xApM7n9DIsXh+1dWUcfbmIqh9OEr7hq5xdZ+NwEaNYstYGKlIJEkkMTd38Y2ATg7usJfQ9xufsHwA8JxnEO8Kq7b3T3hcAC4MgEjxUpE7H9DH9iNHNowU08wiCu4LT95nHyoEvp0LHI/xoiGSGRDunfzGwXohvizKwJsCN9DteY2WXAFOBGd/8RqAd8FrPP0qgsHzPrDHQGaKAZyaQMffMN7McSnqQn7RnFFxzEsbzCp3YsOZo5VcqZRGoOdxKm7q5vZtnABOCvJXy//kAToBWwHHikuCdw94HunuXuWXXq1ClhGCLFtHkzd9fqy1wO5HTe4RYepDXT+IRjNWuqlEtxaw5mVgnIXQWuDaE56Tp3/74kb+buW++sNrNngbHR5jKgfsyu+0VlIqk3eTJ06cLff5zO25Xa0S3nKRbRGNCsqVJ+FbWGdA7wV3df7e5vRGs7lCgxAJjZPjGb5wK5I5lGAxebWdXoJrumwOSSvo9IqVizJizE0KYNrFwJI0aweshYvGFjzZoq5V4ifQ7jzewmYBjwS26hu8ftbDazV4CTgL3MbCmheeokM2tF6L9YRHQznbvPNrPhhOVINwM93D3O3JYiyZM91Pms1zBu/74Xe7OS9w++lmvX3MOcC3ajQYNQU1BCkPLOvLCJYXJ3MFtYQLG7+/7JCSlxWVlZPmXKlFSHIeVAdnYYjbTT4gU8TQ/O4B0+J4uuDGAah2+3b/XqqjFIZjOzqe6eFW+fIjuk3b1xAY+UJwaR0pKdDddcvZGOi+9jFgdFy3Q+SRs+y5cYIAxn7d07BYGKlKEim5WiO6RvABq4e+doKo1m7j62iENFMsKYGyfyyYZuHMg8hnMB1/MYy9k37jFFTdMtkumKc4f0MdF2ondIi6S1Ef1XMaLG5by64mSqspEzeZOLGF5kYgA0fFXKvaTdIS2SrrJfyuG6Gs9zcvfmnPNLNg9wGwcxi7c5M6HjNXxVKoJEkkNp3yEtkjJjH5pNw8tP4vFfrmIOLWjFdHrzABuovt1+Fv3507AhdOsWfmr4qlQkiQxlzXuH9LHA5ckMSqS05I5CWrV4Pfftch/XbOjDT+zGFTzPYC7HC/j7qGFDDVcVSWQ9h3FmNo1SuENapCzlzqJ6wvq3eJce7L9hIYPpxM304XsKnnqlYUNYtKhs4xRJR4UmBzNrnadoefSzgZk1cPdpyQtLZMdkZ8Otl33LCznXcyH/Yi7NOYn3eJ+TCj1GfQki28SrOeROilcNyAJmEGoOhxBmVD26kONEUurll7Yw7cp+zM7pzc78xt+4lz7czG9ULfSYPfeExx9XU5JIrkKTg7ufDGBmrwGt3f2LaPsg4K4yiU6kmN68bxpN/96FS5jCfzidHjzNVxxQ6P6VK8OQIUoKInkl0iHdLDcxALj7LDM7MIkxiRTbsGd/Ys21f+eqX59iFXW4mFcYxkXEG3WtaTBECpfIUNYvzOw5MzspejwLzEx2YCKxctduNoOddgo/K1UCM+f/bCTHdT6Qq399kgF0pTnzGMbF5E0MZqH5SENSRYqWSM3hcqAbcF20/QFh0R6RpMvOhuuug9Wrt5VtiebrbeCLeIprOIs3mM6hnMdrTOaoAs+jWoJI8RS12E9l4K2o/+HRsglJJMgdipq7bnOundjEDfTlTu4mh0rcwCM8wbVsKeSfc+XKSgwixRU3Obj7FjPLMbPd3X1tWQUlFVdBNYVYx/AxA+jKwcziddpzHY+zhMInOlKNQaRkEmlWWkfodxjH9ov9XJu0qKRCKSohAOzBDzzELVzNc3xDfc5mFGM4O+55NTxVpOQSSQ6vRQ+RUte9O/SP24PldGQoj3AjtfmBPtzE3dzJL9Qo9AhNfyGy4xJJDsNg60DxBe7+axLjkQokOxsGDCj89d8zn/504xTe41PacBrjmMmhVKoE5IRRR7kLGaqWIFK64k2fsRPwAHAFsJgwLrC+mb0A9Hb3TWUTopRXvXtv+3KPVZVfuY1/cCsPsp7qdGEAz3I1TiX23BO+18xeIkkX7z6HPkBtoLG7H+7urYEmQC3gn2UQm5RzBa2mdirjmckh3Mk9jOB8mjOPgXTBqUT16qF2ICLJFy85nAVc7e4/5xa4+0+Eex7aJTswKf9iV1PbmxW8REfGcxqG05ZxdCSb7yvVBXTTmkhZi5cc3D1/pd/dtxAt/BOPmQ0ys5VmNiumrLaZjTOzL6Ofe0TlZmZPmNkCM5tZwIywUs5kZ8O6dWDk0JlnmEdzLuBf3M0dHL3rF/xlaFvcww1v7mEabSUGkbITLznMMbPL8haaWUdgXgLnHgz8IU/ZrcAEd28KTIi2Ac4EmkaPzugO7HIt9+a2fVfP5GOO5Rm6Mp1WnFhrBgcMvZvv11VTIhBJsXijlXoAr5nZFcDUqCwL2AU4t6gTu/sHZtYoT/E5sHVC/SHAROCWqPzFqKbymZnVMrN93H05Uu7cd9sv3LX+LnrxKD+yB5fyIkPpSMPdTUlBJE3Em7J7GXCUmZ0CtIyK33T3CTvwfnVjvvC/A+pGz+sBS2L2WxqV5UsOZtaZULugQYPC74yVNDV6NG8v6UlDvuFZruIWHuJHagMFd1CLSGokskzou8C7pf3G7u5mVmTfRQHHDQQGAmRlZRX7eEmN159YQvVbr+WMDf/mJw7iWD7iE47dbh/lepH0kchNcKVpRW5zkZntA6yMypcB9WP22y8qk0y3eTNTL3+S07L/TiVyuIUH6csNbKbKdruZaYlOkXSSyHoOpWk00Cl63gkYFVN+WTRqqQ2wVv0N5cDkyXDEERyefQPvcyItmc3D3JIvMUAYkaT+BpH0kbTkYGavAJ8CzcxsqZldCTwInGZmXwJto22AN4GvgQXAs0D3ZMUlZWDtWujRA2/ThuUzVvJ/jOAsxrKIxoUe0rBhGcYnIkVKWrOSu/+5kJdOLWBfJ4yOkkzmDsOGQa9e5KxYyRN+LXdwDz+zW9zD1KQkkn7Kus9ByqlRfb+i5q3dOWXTO0zhcLowlmkcXuRxZtC1q5qURNKNkoPsmI0bmdGxD6ePuJ9NVKEnT9CP7uRQuchDNbW2SPpScpCSe/996NqVQ+fNYzgX0ItH+ZZ6RR7WsGGYDkNE0ldZj1aS8uD77+Hyy+Gkk1j8v19pxxtcxPCEEsPOO6t/QSQTKDlI4txh0CA2Nm7GpiHZPMBtHJgzm7cSnKS3Rg0YNEjNSCKZQM1Kkpg5c0LP8Ycf8rkdRxcGMGfrrCpF69YN+vVLYnwiUqpUc5D41q8PS7a1agWzZ3Ptrs9zgr+fcGLYc08YOlSJQSTTqOYghXv7bejeHRYu5JWqnbj2hz58T52EDlWns0hmU81B8vv2W7joIjjzTL77sSon8R6XbByccGJQp7NI5lNykG22bIGnnoIDD2TTyFH8jXtpuGY6729dgqNoe+6pTmeR8kDJQcjOhlNqTePzndpAz56889NRHLhlFvfzN36japHHV6oU+hXcwyhXJQaRzKfkUMENe+5nVl92PePWHkF9lnAxr3AG/+ErDkjoeDN48UUlBJHyRh3SFZU7vP46J3S9lro53zKArtzOA6ylVrFOo3mRRMonJYeKaNEi6NkTxo5lBYfSnpFM5qhinWLPPeHxx5UYRMorNStVJJs2wcMPs7l5S9aNfY8beIQsphQrMXTrpr4FkYpANYeK4pNPWHNRF2otncUY2nMdj7OExBdtrlYNnntOCUGkolDNoRzLzob6u/7AQOsMxx7LT0vXcjajOI/XE04MuXc4b9igxCBSkSg5lEPdu4OZ83bHl5i6vjlXMIg+3EQL5jCGs4s83kxDU0UqOjUrlSPZ2dCpEzTZMp/xdOdU3uUzjuI0xjGTQxM+j0YgiYhqDuVAdjZUrQpXdvyVv2+5k5kcwuFMpSv9OYZPipUYNHuqiECKag5mtgj4GdgCbHb3LDOrDQwDGgGLgAvd/cdUxJcJuneH/v23bZ/CBPrTjd/zJdlcwo08wgp+l/D5NDRVRGKlsuZwsru3cvesaPtWYIK7NwUmRNtSgJYttyWGvVnBS3RkAm0xnNN4h45kJ5wYNDRVRAqSTs1K5wBDoudDgPapCyU9ZWeH1dTmzAEjh6sZyDyacwH/4m7u4GC+YDynFXmeSpW2JQU1IYlIQVKVHBx4x8ymmlnnqKyuuy+Pnn8H1C3oQDPrbGZTzGzKqlWryiLWlMlNBmbh0bEj/PILHMxMPuI4BtKF6bTiUGZwF3ezkWqFnit3SKp7mHxVSUFE4klVcjjO3VsDZwI9zOyE2Bfd3QkJJB93H+juWe6eVadOYusLZJowFHVbMshVnV94iL8yjdY05Usu5UVO4V3m07zQc6nZSERKIiUd0u6+LPq50sxeB44EVpjZPu6+3Mz2AVamIrZUa9kyNBvldRZjeIpraMg3PMtV3MJD/EjtQs9z6qkwfnwSAxWRcq3Maw5mtquZ1cx9DpwOzAJGA52i3ToBo8o6tlRq2zbUFvImhv1YwkjOYwxn8xO7cSwf0ZlnC00MNWqE5iMlBhHZEamoOdQFXjez3Pd/2d3fNrPPgeFmdiWwGLgwBbGVuexsuPTS0PQTqzKb6cmT3MMdVGYLt/AgfbmBzVTJdw7dmyAipa3Mk4O7fw3578py99XAqWUdTyoVlhiOYDLP0IXDmM4btOManmIRjfMdbwYvvaS+BBEpfek0lLXCyL2juWPH7RPDbqzlSa7hM9qwNys5n39xFmMLTAw1aigxiEjyaG6lMta2LUyYkLfUuZDhPMb17M1KnqQnf+defma3rXtoymwRKUuqOZSB3JqCWf7EsD9f8RZnMoyLWUY9jmQy1/P41sSw006aMltEyp6SQ5JlZ4fmo99+2768Cr9xO/czi4M4hk/oyRMcxSSmcfjWfbp1C4u3KSmISFlTs1KSXXVV/rITeJ/+dKMFcxnOBfTiUb6l3tbXd94ZNm4swyBFRPJQzSFJsrNDM9Kvv24r25PvGcRfeJ+T2IUNtOMNLmL4donBDAYNSkHAIiIxlBySILcpaRvnLwxiPs3oyFAe4DZaMpu3aLfdcRqBJCLpQs1KpSg7G664Yvv+hQOZwwC6cgIf8iHH0ZUBzKHldsfpJjYRSTdKDqUkb22hGhv4G/dxM334mZpcyXO8wF/wmMqa5j8SkXSlZqVS0L379onhDN5mFgfRmwd4mUtozjwGceV2iaFbNyUGEUlfqjnsgLw3tP2O5TzG9VzEcObRjJN4j/c5Kd9xakYSkXSn5FACeafVrsQWujKAB7idqmzkb9xLH27mN6pud9y++8KyZWUcrIhICahZqZiqV98+MRzGND7laJ7mGiZxFAcxi/v5W77E0K2bEoOIZA7VHBJQ0AI8NfiZe7iDa3mCVdThYl5hGBcBtt1+lSvD5s1lF6uISGlQzSGOghfgcc7lNeZyINfxOM/QhebMYxgXkzcxAAwZUlbRioiUHtUcCrHHHrBmzfZlDVjMU1zDnxjLdA7l/xjJZI4q8PjKlUNi0A1tIpKJVHMoQN7EsBObuIk+zKEFJ/MeN/AIWUwpMDHkzqK6ebMSg4hkLiWHGLnNSLGJ4Wg+YSqH04e/Mo7TaMEcHuUGtuSpdFWrFpKCZlEVkfJAzUpAvXrw7bfbl+3BD/yD2+jCQL6hPufwb0ZzToHHDx2qhCAi5UuFTw477xz+2t/G6UA2fbmB2vzAP7mRu7iLX6iR79iddoLBg5UYRKT8qXDJoeBlOoOm/I/+dONU3uUzjuI0xjGTQ7fbRzeyiUhFkHZ9Dmb2BzObb2YLzOzW0jx3YYmhKr9yJ3fxBQdzOFPpSn+O4ZN8ieHUU5UYRKRiSKuag5lVBp4GTgOWAp+b2Wh3nxP/yMQUlBhOYQL96cbv+ZJsLuFGHmEFv8u3n2ZQFZGKJN1qDkcCC9z9a3f/DXgVCukF3kF1WMlLdGQCbTGc03iHjmTnSwxmocNZiUFEKpJ0Sw71gCUx20ujsq3MrLOZTTGzKatWrSrRm/yBt5hPMy7gX9zNHRzMF4zntHz7desGOTnqcBaRiietmpUS4e4DgYEAWVlZXpJzfElTPqMNvXiU+TTP97qGpopIRZduNYdlQP2Y7f2islLRrVv4+RUH0I638iWGypWVGEREIP2Sw+dAUzNrbGY7AxcDo0vr5P36hQRRuXLYrlw5bLuHh6a8EBEJ0qpZyd03m9k1wH+AysAgd59dmu/Rr59WYRMRKUpaJQcAd38TeDPVcYiIVGTp1qwkIiJpQMlBRETyUXIQEZF8lBxERCQfJQcREcnH3Et0k3FaMLNVwOISHLoX8H0ph5Nq5e2adD3pr7xdU0W6nobuXifewRmdHErKzKa4e1aq4yhN5e2adD3pr7xdk65ne2pWEhGRfJQcREQkn4qaHAamOoAkKG/XpOtJf+XtmnQ9MSpkn4OIiMRXUWsOIiISh5KDiIjkU+GSg5n9wczmm9kCM7s11fGUhJktMrMvzGy6mU2Jymqb2Tgz+zL6uUeq44zHzAaZ2UozmxVTVuA1WPBE9JnNNLPWqYu8YIVcz11mtiz6nKabWbuY126Lrme+mZ2RmqgLZ2b1zew9M5tjZrPN7LqoPCM/ozjXk8mfUTUzm2xmM6Jrujsqb2xmk6LYh0Vr42BmVaPtBdHrjeK+gbtXmAdhjYivgP2BnYEZQItUx1WC61gE7JWn7GHg1uj5rcBDqY6ziGs4AWgNzCrqGoB2wFuAAW2ASamOP8HruQu4qYB9W0T/9qoCjaN/k5VTfQ15YtwHaB09rwn8L4o7Iz+jONeTyZ+RATWi51WASdHvfjhwcVQ+AOgWPe8ODIieXwwMi3f+ilZzOBJY4O5fu/tvwKvAOSmOqbScAwyJng8B2qculKK5+wfAD3mKC7uGc4AXPfgMqGVm+5RJoAkq5HoKcw7wqrtvdPeFwALCv8204e7L3X1a9PxnYC5Qjwz9jOJcT2Ey4TNyd18XbVaJHg6cAoyIyvN+Rrmf3QjgVDOzws5f0ZJDPWBJzPZS4v8DSVcOvGNmU82sc1RW192XR8+/A+qmJrQdUtg1ZPLndk3UzDIopqkvo64nan44jPCXacZ/RnmuBzL4MzKzymY2HVgJjCPUcNa4++Zol9i4t15T9PpaYM/Czl3RkkN5cZy7twbOBHqY2QmxL3qoN2b0GOXycA1Af6AJ0ApYDjyS0mhKwMxqACOB6939p9jXMvEzKuB6Mvozcvct7t4K2I9Qs2leWueuaMlhGVA/Znu/qCyjuPuy6OdK4HXCP4oVudX46OfK1EVYYoVdQ0Z+bu6+IvrPmwM8y7ZmiYy4HjOrQvgizXb316LijP2MCrqeTP+Mcrn7GuA94GhCk17uEtCxcW+9puj13YHVhZ2zoiWHz4GmUW/+zoROmdEpjqlYzGxXM6uZ+xw4HZhFuI5O0W6dgFGpiXCHFHYNo4HLohExbYC1MU0baStPm/u5hM8JwvVcHI0eaQw0BSaXdXzxRG3RzwNz3b1vzEsZ+RkVdj0Z/hnVMbNa0fNdgNMIfSnvAedHu+X9jHI/u/OBd6PaX8FS3eNe1g/CqIr/Edrmeqc6nhLEvz9hFMUMYHbuNRDaDicAXwLjgdqpjrWI63iFUI3fRGgXvbKwayCMyng6+sy+ALJSHX+C1/NSFO/M6D/mPjH7946uZz5wZqrjL+B6jiM0Gc0EpkePdpn6GcW5nkz+jA4B/hvFPgu4Iyrfn5DIFgD/AqpG5dWi7QXR6/vHO7+mzxARkXwqWrOSiIgkQMlBRETyUXIQEZF8lBxERCQfJQcREclHyUHSgpntGTMz5ncxM2WuMbM5ZRxLezNrEbN9j5m1LcF5GlnMLK1lzcxuz7P9SfQzpXFJZlBykLTg7qvdvZWHqQAGAI9Gz1sBOaX9fjF3kBakPWFWztzY7nD38aUdQxnYLjm4+zGpCkQyj5KDZILKZvZsNGf9O9HdoJhZEzN7O5qA8EMzax6VNzKzd6PJ1CaYWYOofLCZDTCzScDDBR1vZscAZwN9oppLk+i486NzHGFmn0Rz6E82s5rR+31oZtOiR9wv4egu4qcsrBMw3szejDn/IjPbK3qeZWYTo+dHmtmnZvbf6P2bReWXm9lr0XV8aWYPR+UPArtE15Adla0rIJbKZtbHzD6Pfl9dovJ9zOyD6PhZZnb8Dn6GkmlSfZefHnrkfRAzxz7QCNgMtIq2hwMdo+cTgKbR86MI0wEAjAE6Rc+vAP4dPR8MjCWalz/O8YOB82PiGUyYbmBn4GvgiKh8N2AnoDpQLSprCkyJiX1WAdd3HmEGzcrAvsCa3PcjZq0OIAuYGPte0fO2wMjo+eVRTLsT7oBdDNSPXluX533X5Y0L6Az8LXpeFZhCWL/gRrbdfV8ZqJnqfxd6lO0jXtVaJF0sdPfp0fOpQCMLs2seA/zLtk1JXzX6eTThCxjC9AgPx5zrX+6+pYjjC9MMWO7unwN4NEtpNMfVU2bWCtgC/L6I85wAvOLuW4BvzezdIvaH8OU/xMyaEqaBqBLz2gR3XxvFMgdoyPbTTcdzOnBIbs0lep+mhHnIBlmYrO7fMb9/qSCUHCQTbIx5vgXYhdAkusZDv0Rx/BL9LOnxBekFrAAOjc776w6cazPbmnurxZTfC7zn7udaWI9gYsxreX8/xfl/bUBPd/9PvhfCVPB/BAabWV93f7EY55UMpz4HyUjRX+0LzewC2NqOf2j08ieEGXcBOgAfFvP4nwlLSeY1H9jHzI6Ijqlp26Y+Xu5h2udLCc0w8XwAXBS19+8DnBzz2iLg8Oj5/8WU7862qZcvL+L8uTZFf/nH8x+gW+5+ZvZ7CzP/NgRWuPuzwHOEJVClAlFykEzWAbjSzHJnqM1d8rUn8Bczm0n4sr6umMe/Ctwcdf42yd3Zw9KyFwFPRseMI/x13w/oFJU1Z1vtpDCvE2Y1nQO8CHwa89rdwONmNoVQC8j1MPAPM/svidcMBgIzczukC/FcFMe0aHjrM9H5TwJmRO93EfB4gu8p5YRmZRVJMTMbDIx19xFF7StSVlRzEBGRfFRzEBGRfFRzEBGRfJQcREQkHyUHERHJR8lBRETyUXIQEZF8/h/bBiQbwofaQwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_runs = 250\n",
    "times_h0 = [time_run(cd, X_h0, window_size) for _ in range(n_runs)]\n",
    "print(f\"Average run-time under no-drift: {np.mean(times_h0)}\")\n",
    "_ = scipy.stats.probplot(np.array(times_h0), dist=scipy.stats.geom, sparams=1/ert, plot=plt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And when drift has occured the detector is very fast to respond."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average run-time under drift: 4.328\n"
     ]
    }
   ],
   "source": [
    "n_runs = 250\n",
    "times_h1 = [time_run(cd, X_corr, window_size) for _ in range(n_runs)]\n",
    "print(f\"Average run-time under drift: {np.mean(times_h1)}\")"
   ]
  }
 ],
 "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"
  },
  "metadata": {
   "interpreter": {
    "hash": "ffba93b5284319fb7a107c8eacae647f441487dcc7e0323a4c0d3feb66ea8c5e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}