intro_to_drift_detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3efcd04f",
   "metadata": {
    "id": "6a288462"
   },
   "source": [
    "# ODSC 2022: An Introduction to Drift Detection\n",
    "\n",
    "To load this notebook in Google Colab [CLICK HERE](https://colab.research.google.com/github/ascillitoe/odsc_workshop/blob/main/intro_to_drift_detection.ipynb).\n",
    "\n",
    "The notebook will run most efficiently on a GPU. To select a GPU runtime on Colab navigate to ***Runtime -> Change runtime type*** and select ***GPU*** under ***Hardware accelerator***.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this workshop, you'll learn about *drift detection* through a number of tasks:\n",
    "\n",
    "1. [Classifying newsgroups](#1.-Classifying-newsgroups) - you'll train a machine learning model to classify the topic of newsgroup posts.\n",
    "2. [Performance under drift](#2.-Performance-under-drift) - You'll examine the performance of your classifier on test data that has *drifted* from the data it was trained on.\n",
    "3. [Detecting drift](#3.-Detecting-drift) - You'll use the open-source [Alibi Detect](https://docs.seldon.io/projects/alibi-detect/en/latest/) library to detect a number of types of drift in a principled manner. \n",
    "4. [Accounting for context](#4.-Accounting-for-context) - To finish, you'll learn how to use a powerful drift detector which is able to take relevent context into account.\n",
    "\n",
    "We'll also leave you with some [homework](#Homework!) to cement your newfound knowledge. After this, you'll have a good working knowledge drift detection, and will be able to apply it to your own ML tasks!\n",
    "\n",
    "**Note:** This notebook has some sections left blank for you to fill in (marked with `???`). If you'd rather follow a completed notebook, this can be accessed [here](https://github.com/ascillitoe/odsc_workshop/blob/main/intro_to_drift_detection_master.ipynb), or loaded into Google Colab by [clicking here](https://colab.research.google.com/github/ascillitoe/odsc_workshop/blob/main/intro_to_drift_detection_master.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15aa76b1",
   "metadata": {
    "id": "bf3c9ceb"
   },
   "source": [
    "## 0. Getting ready\n",
    "\n",
    "This workshop relies on a number of files located in a [github repo](https://github.com/ascillitoe/odsc_workshop). If running locally, start by cloning the repo, and moving to the workshop directory:\n",
    "\n",
    "```\n",
    "git clone https://github.com/ascillitoe/odsc_workshop.git\n",
    "cd odsc_workshop\n",
    "```\n",
    "\n",
    "If running on Google Colab, the following code will fetch the necessary files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "507974f5",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "09fa8f99",
    "outputId": "76af0e96-cf70-402f-9964-2cabfd44e467"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "  import google.colab\n",
    "  IN_COLAB = True\n",
    "except:\n",
    "  IN_COLAB = False\n",
    "\n",
    "if IN_COLAB:\n",
    "    import zipfile\n",
    "    import shutil\n",
    "    import os\n",
    "\n",
    "    !wget -nc https://github.com/ascillitoe/odsc_workshop/archive/refs/heads/main.zip\n",
    "\n",
    "    pz = open('main.zip', 'rb')\n",
    "    packz = zipfile.ZipFile(pz)\n",
    "    packz.extractall()\n",
    "    pz.close()\n",
    "\n",
    "    srcdir = 'odsc_workshop-main/'\n",
    "    for filepath in os.listdir(srcdir):\n",
    "        if not os.path.exists(filepath):\n",
    "            shutil.copyfile(os.path.join(srcdir, filepath), filepath)\n",
    "    shutil.rmtree(srcdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa5a4d63",
   "metadata": {
    "id": "e4617fcb"
   },
   "source": [
    "### Software\n",
    "\n",
    "In this workshop we'll make use of state-of-the-art drift detectors from the open-source [Alibi Detect](https://docs.seldon.io/projects/alibi-detect/en/latest/) library. To install Alibi Detect, and the other dependencies used in this workshop, run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "466d70bd",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "7b1f5060",
    "outputId": "ee1747c7-4c36-4785-8fb8-e9f6f194003e"
   },
   "outputs": [],
   "source": [
    "!pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daf6f4c7",
   "metadata": {
    "id": "7f65feb2"
   },
   "source": [
    "### Download data and sentence transformer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df71af48",
   "metadata": {},
   "source": [
    "To download the [twenty newsgroups dataset](https://archive.ics.uci.edu/ml/datasets/Twenty+Newsgroups) and the [sentence transformer model](https://huggingface.co/sentence-transformers/paraphrase-MiniLM-L6-v2) (more on these later), run the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e8997a7",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "70af6555",
    "outputId": "418b4ac5-20aa-4c0d-ddd3-109c42988113"
   },
   "outputs": [],
   "source": [
    "# Fetch preqrequisites\n",
    "from workshop_utilities import fetch_prerequisites\n",
    "dataset, sentence_transformer = fetch_prerequisites()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b052fe58",
   "metadata": {
    "id": "43d2fd99"
   },
   "source": [
    "## 1. Classifying newsgroups\n",
    "\n",
    "### The data\n",
    "\n",
    "The dataset we'll explore in this workshop is the [20 newsgroup dataset](https://archive.ics.uci.edu/ml/datasets/Twenty+Newsgroups). [Usenet newsgroups](https://en.wikipedia.org/wiki/Usenet_newsgroup) have been around since the dawn of the internet, and can be thought of as the very first online social network! \n",
    "\n",
    "The dataset should contain about 18,000 newsgroup posts across 20 topics, including politics, science, sports and religion. Check this now! (*Note: The actual posts are stored in a list `dataset.data`, and a list of topic names is stored in `dataset.target_names`*):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb40f101",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'{len(dataset.data)} posts')\n",
    "print(f'{len(dataset.target_names)} categories:')\n",
    "classes = dataset.target_names\n",
    "classes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c926c45",
   "metadata": {},
   "source": [
    "To get an idea of what a post looks like, print out an instance and it's topic (*Note: The topic labels are stored as integers in `dataset.targets`. These can be used to index `dataset.target_names`*):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d45c9e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 0  # instance index\n",
    "post = dataset.data[i]\n",
    "topic_idx = dataset.target[i]\n",
    "topic = classes[topic_idx]\n",
    "\n",
    "print(post)\n",
    "print(f'Topic: {topic}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66fa082d",
   "metadata": {},
   "source": [
    "So, each data instance looks something like an email. We wish to train a machine learning model that can take in future instances (without labels), and predict their topic. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b594e1c",
   "metadata": {
    "id": "9d27f700"
   },
   "source": [
    "### Visualising the embeddings\n",
    "\n",
    "For classification, we'll use a simple [MLP](https://en.wikipedia.org/wiki/Multilayer_perceptron) model. But before that, we'll transform the posts into *semantically meaningful sentence embeddings* (*semantically meaningful* means that semantically similar sentences are close in vector space).\n",
    "\n",
    "To extract embeddings, we'll use the pre-trained `'paraphrase-MiniLM-L6-v2'` model from the [SentenceTransformers](https://www.sbert.net/index.html) package. To see this model in action, you can checkout the [hosted HuggingFace API](https://huggingface.co/sentence-transformers/paraphrase-MiniLM-L6-v2). The `sentence_transformer` should already have been loaded earlier, so we just need to pass it to the `EmbeddingModel` wrapper class via the `model` kwarg:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c268f0",
   "metadata": {
    "id": "1f5ea9e8"
   },
   "outputs": [],
   "source": [
    "from workshop_utilities import EmbeddingModel\n",
    "emb_model = EmbeddingModel(model=sentence_transformer)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23dd5093",
   "metadata": {},
   "source": [
    "The sentence embeddings are high dimensional (384-d to be precise!). To visualise them, we need to first reduce their dimensionality. For this, we will train a [UMAP](https://umap-learn.readthedocs.io/en/latest/) (Uniform Manifold Approximation and Projection) model. UMAP is able to take advantage of our data labels, so we'll need to pass it the data labels (from `dataset.target`) in addition to the embeddings.\n",
    "\n",
    "We won't need all of the 18000+ posts for training the UMAP model, so subsample the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "924968f0",
   "metadata": {
    "id": "8a81218b"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from workshop_utilities import set_seed\n",
    "set_seed(2022)  # This will ensure reproducibility (at least on CPU!)\n",
    "\n",
    "n_train = 2000  # A larger sample here will give a better UMAP projection, but can reduce if too slow\n",
    "\n",
    "idx_train = np.random.choice(len(dataset.data), size=n_train, replace=False)\n",
    "x_train, y_train = [dataset.data[_] for _ in idx_train], dataset.target[idx_train]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4be61ee",
   "metadata": {},
   "source": [
    "The list of training posts (`x_train`) is passed through the `emb_model` to transform them into embeddings. Print out the shape of `emb_train` here to check its dimensions (they should be `(n_train, 384)`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c360dd5",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "593a8820",
    "outputId": "7b0237e5-0b1a-442a-a867-448a4b35ddba"
   },
   "outputs": [],
   "source": [
    "emb_train = emb_model(x_train)\n",
    "emb_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4768afd7",
   "metadata": {},
   "source": [
    "The UMAP model is provided for you in `workshop_utilities`. To train it, the embeddings `emb_train` and labels `y_train` are passed to the model's `.fit()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f106e3c6",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "e87b4ca5",
    "outputId": "0e9c2c70-943f-4e2a-d6f0-ea731eb38e22"
   },
   "outputs": [],
   "source": [
    "from workshop_utilities import UMAPModel\n",
    "umap_model = UMAPModel()\n",
    "umap_model.fit(emb_train, y=y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c437444",
   "metadata": {},
   "source": [
    "Now use the UMAP model to obtain a 2-d representation `dr_train`of the 384-d embeddings `emb_train`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc9e8200",
   "metadata": {},
   "outputs": [],
   "source": [
    "dr_train = umap_model.predict(emb_train)\n",
    "dr_train.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11740a36",
   "metadata": {},
   "source": [
    "We now have a 2-d representation of the sentence embeddings, from which we can visually inspect the various newsgroup topics. The `plot_clusters` utility function is provided for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b05b60a",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 607
    },
    "id": "3606812e",
    "outputId": "cbb4c179-34b0-412c-8076-017be098b62f"
   },
   "outputs": [],
   "source": [
    "from workshop_utilities import plot_clusters\n",
    "plot_clusters(dr_train, y_train, classes, title='Training data: clustered news topics')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80e46140",
   "metadata": {},
   "source": [
    "Hopefully your data is nicely clustered by topics! The clustering is partly down to the effectiveness of the sentence transformer in extracting semantically meaningful embeddings, but it is also aided by the UMAP algorithm which has a clustering effect when used in a supervised manner(see [here](https://umap-learn.readthedocs.io/en/latest/supervised.html#using-labels-to-separate-classes-supervised-umap))."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8f0ace8",
   "metadata": {
    "id": "799ebc30"
   },
   "source": [
    "### Training a classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a44c27c",
   "metadata": {
    "id": "d394b2dc"
   },
   "source": [
    "Now to train a classifier!\n",
    "\n",
    "In the interests of time, the `Classifier` model has been defined for you in `workshop_utilities`. The classifier is initialised with the `EmbeddingModel` you initialised earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6f4e02f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from workshop_utilities import Classifier\n",
    "import torch\n",
    "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "clf = Classifier(embedding=emb_model).to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b313c4e3",
   "metadata": {},
   "source": [
    "The `Classifier` consists of a pre-trained sentence transformer model (the model you passed to `embedding`), and an MLP head. These two components are stored in the attributes `.embedding_model` and `.head`. Print out the attributes to have a look for yourself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dde1353",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(clf.embedding_model)\n",
    "print('')\n",
    "print(clf.head)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e931865",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "48f64f05",
    "outputId": "2c3e5148-6dbd-487d-8af0-18d85a673224"
   },
   "source": [
    "To demonstrate *drift*, we want to start by only training the classifer on posts from two topics; cryptology (11) and space (14). The utility function `split_data` splits a dataset with a specified number of instances per class (i.e. per topic). We can use this to generate *train* and *test* data consisting of only cryptology and space posts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b6c7adb",
   "metadata": {
    "id": "98f3a16f"
   },
   "outputs": [],
   "source": [
    "from workshop_utilities import split_data\n",
    "\n",
    "n_classes = len(classes)\n",
    "n_train_c = [0] * n_classes  # Sample 0 from all topics (will update below!) \n",
    "n_test_c = [0] * n_classes\n",
    "\n",
    "n_train_c[11], n_train_c[14] = 200, 200  # subsample topics 11 & 14 only\n",
    "n_test_c[11], n_test_c[14] = 100, 100  # subsample topics 11 & 14 w/ sample proportions as above\n",
    "\n",
    "(x_train, y_train), (x_test, y_test), _ = split_data(dataset.data, dataset.target, n_train_c, n_test_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34dab2e1",
   "metadata": {},
   "source": [
    "Check `y_train`, does it consist of only 11's and 14's?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "029ba0fb",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "deb8d8d2",
    "outputId": "90e326e7-7023-4f68-f7ce-8ed8e6a6439e"
   },
   "outputs": [],
   "source": [
    "np.unique(y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "463ab323",
   "metadata": {},
   "source": [
    "Use the `plot_clusters` function to visualise the training data. Is it as expected?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8711ecad",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 607
    },
    "id": "bff64cbd",
    "outputId": "125facde-5b8d-4fdf-f877-ff1618401adf"
   },
   "outputs": [],
   "source": [
    "emb_train = emb_model(x_train)\n",
    "plot_clusters(emb_train, y_train, classes, dr_model=umap_model, title='Training data: clustered news topics')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6c5f66a",
   "metadata": {},
   "source": [
    "If everything looks good, `x_train` and `y_train` can now be used to train the classifier! Note that training only involves the MLP head, the pre-trained embedding backbone is left untouched.\n",
    "\n",
    "In the interests of time, the classifier has actually been trained for you. If you'd rather perform the training yourself, set `TRAIN_CLF = True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db7d97c6",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "9a679e3b",
    "outputId": "feec2a86-2365-42a1-b561-7509ff4146a6"
   },
   "outputs": [],
   "source": [
    "from workshop_utilities import train_model, eval_model\n",
    "\n",
    "TRAIN_CLF = False  # Set to TRUE to train classifier, otherwise it will be loaded from disk\n",
    "filepath = 'classifier'\n",
    "if TRAIN_CLF:\n",
    "    # Train model\n",
    "    train_model(clf, x_train, y_train, epochs=5, shuffle=True)\n",
    "    clf.eval()\n",
    "    # Save model\n",
    "    torch.save(clf.state_dict(), filepath)\n",
    "    # Check train and test accuracy\n",
    "    _, _ = eval_model(clf, x_train, y_train, shuffle=True)\n",
    "    _, _ = eval_model(clf, x_test, y_test, shuffle=True)\n",
    "else:\n",
    "    # Load model\n",
    "    clf.load_state_dict(torch.load(filepath, map_location=device))\n",
    "    clf = clf\n",
    "    clf.eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4ada9dd",
   "metadata": {
    "id": "b63002d4"
   },
   "source": [
    "### Testing the classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b0174d5",
   "metadata": {},
   "source": [
    "We can now run the classifier on previously unseen instances from the *test* set. If you use `plot_clusters` on the test set, you'll see that it's distribution is qualitatively similar to the *training set*; **the instances are drawn from the same distribution, but the instances are not the same**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b17a7d4e",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 607
    },
    "id": "6ddda901",
    "outputId": "c6ca43ee-5574-4894-e220-323ec12d80d5"
   },
   "outputs": [],
   "source": [
    "emb_test = emb_model(x_test)\n",
    "plot_clusters(emb_test, y_test, classes, dr_model=umap_model, title='Test data: clustered news topics')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5d7c50e",
   "metadata": {
    "id": "0ef79da1"
   },
   "source": [
    "To convince yourself the classifier is working, you can run it on a few instances and examine its predictions (*Note: The classifier outputs logits, so we have to run .argmax(1) on the output to obtain a class prediction)*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d971576b",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 35
    },
    "id": "c34d86d4",
    "outputId": "1f41a5a5-6537-4f18-ee79-ffb41ec77325"
   },
   "outputs": [],
   "source": [
    "idx = 42\n",
    "class_pred = clf([x_test[idx]]).argmax(1)\n",
    "classes[class_pred]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0ba85b7",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "4ec4bf0c",
    "outputId": "bfd9eff6-14d9-46d6-b983-d33d91b4872e"
   },
   "outputs": [],
   "source": [
    "print(classes[y_test[idx]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef8d33be",
   "metadata": {
    "id": "6cf7b6d3"
   },
   "source": [
    "## 2. Performance under drift\n",
    "\n",
    "To artificially introduce drift, lets now introduce a new topic into our test data. Use the `split_data` function to create two datasets:\n",
    "\n",
    "1. `(x_nodrift, y_nodrift)` - A data set with the same topic prevalence as the original test data.\n",
    "2. `(x_drift, y_drift)` - A data set with the same number of instances of topics 11 and 14, but also 100 instances of topic 18 (talk.politics.misc)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a0376ee",
   "metadata": {
    "id": "6efc6219"
   },
   "outputs": [],
   "source": [
    "n_classes = len(classes)\n",
    "n_nodrift_c = n_test_c\n",
    "n_drift_c = [0] * n_classes\n",
    "\n",
    "n_drift_c[11], n_drift_c[14], n_drift_c[18] = 100, 100, 100  # \n",
    "\n",
    "(x_nodrift, y_nodrift), (x_drift, y_drift), _ = split_data(dataset.data, dataset.target, n_nodrift_c, n_drift_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e666507",
   "metadata": {},
   "source": [
    "Use `plot_clusters` to check the new *_drift* data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9df85dfe",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 607
    },
    "id": "d8d21d87",
    "outputId": "b722ef58-a761-4e19-c780-ea26ccac728b"
   },
   "outputs": [],
   "source": [
    "emb_drift = emb_model(x_drift)\n",
    "plot_clusters(emb_drift, y_drift, classes, dr_model=umap_model, title='Drifted data: clustered news topics')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "859ff7de",
   "metadata": {
    "id": "ef1c392e"
   },
   "source": [
    "It should be clear that the new `_drift` dataset distribution is quite different to the previous test data distribution, in that there is a new topic cluster. To see the effect of this *data drift* on our classifier, use the `eval_model` function on the *_nodrift* and *_drift* data sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5d292fe",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "535bc9e4",
    "outputId": "a894fe07-e757-4a47-cfbb-6133c11e1cdf"
   },
   "outputs": [],
   "source": [
    "_, _ = eval_model(clf, x_nodrift, y_nodrift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aec9ec54",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "dbebce6a",
    "outputId": "75268688-0f65-403e-8b31-4fab52051c3d"
   },
   "outputs": [],
   "source": [
    "_, _ = eval_model(clf, x_drift, y_drift)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bfa2595",
   "metadata": {},
   "source": [
    "Unsurprisingly, the classifier does not perform well when a new topic not seen during training is introduced!"
   ]
  },
  {
   "attachments": {
    "drift_detection_pipeline.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABMoAAAC+CAYAAADED+P9AAAACXBIWXMAABibAAAYmwFJdYOUAAAAGXRFWHRTb2Z0d2FyZQB3d3cuaW5rc2NhcGUub3Jnm+48GgAAIABJREFUeJzs3Xd4VGX2wPHve2cmmRSSkBBAQVqAkAIroiAEqSIiywoIqKsgLIvguhbUFRAEbCxg+dkbq1JdFSvYQJo0EUVpSQAh9KUYSkLaJDP3/f1xM5MEEpoJk8D5PM99yNy5971nAiQz5573vIpykDh98j2gXi+PsYTwJ22aXZP/NmaJv+MQQgghhBBCCCHEhWcvt5GU0r3irzxRbuMJcQFlu/KMJdtTQ/0dhxBCCCGEEEIIIfyn3BJlSkPrejG55TWeEBfSkewTNkmUCSGEEEIIIYQQl7bzTpQ1nz51POjm1iMdg9bq/V9Wh5dXYEJcCAF2u+7XonWmv+MQQgghhBBCCCGE/513okyhuwYHBCZFVwsr0Nq0oxX5Hk9geQYnREU6npOtcgryNSCJMiGEEEIIIYQQQvyBqZcK6kfWKLitZduMAtMT5CooCAsJdOpyjE2ICrVixxZW7tzm7zCEEEIIIYQQQghRSZx1oizho4kBKsfZ2/tYQ3RGbo5tw/7dTlPj8GgPTrtDVUyYQpSfyJBq+rKwCH+HIYQQQgghhBBCiErmrBNlBe7QsAA8H3ofaw37Mo7x8cafivclk0SZqPSuqdeIXglXSfWjEEIIIYQQQgghSjjnqZd9W1yT0axWHZfLnR9hoBwOu0MrUFprpZTkyUTl9u6Py/wdghBCCCGEEEIIISqpMyfKtFYA5txXFPm5OAybdtrtWmsTQykC7HYUoLVGEmWislKnfCGEEEIIIYQQQghR0mkTZYnTJ9/DjKmvF9/34fo1EawvcZikHkSlN6RNRxpFRvs7DCGEEEIIIYQQQlRiZ64oU0r3ir/yhNZKeUx3Nbvdrg0UKAVoqSETlVqOy8XCbZv9HYYQQgghhBBCCCGqgDMmypSG1vVick2tjZx8VzWnw4HNsFl9ydAoKSgTldjRnGxJlAkhhBBCCCGEEOKsnJIoS3xvykSUTgRQ0Ai0mrNuZQQoNBpDGUophZIpl6ISC7Q76NP8an+HIYS4RCS89+9rlFKjLtT1tMG05EGjF1yo6wkhhBBCCHGpOCVRZthUl2C7s110tbACtLYDuE0CQAMaExMkSSYqsaO52eTmu6C5vyMRQlwqDIw6Gm6pH1nDZShDV+S1dh457ESzpCKvIYQQQgghxKWqlKmXivpRNQpuu/LajALTE+QqKAgLCXRab/q1ViiZbCkqtxU7trBy5zZ/hyGEuATd1rJtZmhAoFmR1xj/7ScBmgrNxQkhhBBCCHHJsid8NDFA5Th7e3doraMzcnNsG/bvdpoah0d7cNodCvC9LZdEmaiMooKrcVl4hL/DEEJcIuJnPtvGMD31ATT6GlCkHNwXGGizV0gWK6F2XZfdZpMMmRBCCCGEEBXIXuAODQvA86F3h9aafceP8vHxo+HFjpPcmKj0WteLoVd4S3+HIYS4RNhM8wGNur34vvnJv4ZV1PViomv9HiqJMiGEEEIIISqUb+pl3xbXZDSrVcflchdEGOBw2B3WypZaK6UkTyYqt3d+XObvEIQQl6C6EVF5g665LtNjmoEud354cICz3Kdd/nb4f2ruhrXyi1gIIYQQQogLwJcocxg2HWS3a7TGUBBgt6OwKswkUSYqO0OKHoUQF4rW1g+cmc9iKHDa7dqtTa21SaC9lNaf58n7U81hs5XbmEIIIYQQQojTK7939EIIIcRFLnH6lAJmTLWD1bdzz7EjjP/mY2exQ4zyulaDqJr6b607yFRLIYQQQgghLiBJlAkhhBBnSaFUQu26eQ2jauR7TB2EUnaHYVjlX+VYgL1uz87yGUiUasSIEZ1M0+zh7ziEqMqUUivfeuut+f6OQwghhChvkigTQpxJMJDj7yCEqCzqR0bmt64Xk5tXUBCgwe50ODRohabcMmVpR34n05VXLmOJU3k8nrZKqUf9HYcQVZnW2glIokwIIcRFRxJlQojTGQa8CiwHZgGfA5l+jUiICyh++pQhNvRT3sca07ZoW3Lo8h1bQlFKFVaRlVuTxNb1Y+jYKK68hhNnEBwcbMbExBT4Ow4hqpI9e/Y4MjIyym2auRBCCFHZSKJMCHE6LwAngE7A9UAu8DXwIfBl4WMhLlo2ZcRqbdapH1kDW1H7sQr5gJh29DAHT0geuqIMHz48HqtCFq11HYDIyEhPz549s/0amBBVyGWXXeaePn162K+//uo889FCCCFE1SSJMiHE6ZjAOuB24M9Af+Bm4BasJNliYC7wMTI9U1zEbr+qnRnsCACtDcqzhKyYCd9+UgGjCi+t9RzgyuL79u3b55g8eXKkn0ISosp56KGHjvo7hspg6NChkU6nU5ZcFxWuoKAg8+2335bKZyEuMEmUCSHOxlFgZuEWhZUoGwT0xEqgvQzMw0qaLQDy/ROmEEIIIUTFstlshwsKCmz+jkNc/LTWN2O9xxZCXECSKBNCnKsjwNuFW12spFl/4E5gIHAMa1rmXOAbwO2fMIUQ4lT9+/c/kZSUpFwuV6jdbtfBwcHa3zEJUdn94x//MLSW/ypCCCEuDZIoE0L8EfuAlwq3+kBvrKTZwMLtCFZPs5lY0zTlXbYQ4oIYOnRofZvNtsz7WGt9uVKKr776KmTx4sVorVGF/BimEJVaWFgYo0ePNv0dR2U1ePDgjIYNG8q0OFHuJk+eHJmbmyuLZgjhJ5IoE0KUl90UJc3igQHArRQlzfYBn2JVmq1CkmZCiIrlABp4H3jzYTk5OUZOjrRUFOJsuN1SFH461atXN2vXru3xdxzi4iP3cITwL0mUCSEqQgowsXBLwKoy+ytwf+G2B/gcK2m20i8RCiEuGbfeemumw+EI0VrbHA6HttvtaK0VyIcRIUqzc+dOFi9e7O8whBBCCL+QRJkQl5i496Y2MRR9CM55MXnAxAvRdD+5cJtIUdJsIEVJs53AR8AMIPUCxCOEuEjdfffdnYEeAKZpVvfu379/v91utysAm82mDMOQBJkQZWjcuDGNGzeWRJkQQohLliTKhLjE2AwVhPY8RU7QgwnTp4xIHjzqQq6k402aPQm0w0qa9QdGFW4pWFVm7wPbLmBcQoiLQ1vgXwCGUdTaZeXKlcH+CkiIqsbtdtOkSRN/hyGEEEL4jSTKhLgEaaUyleYyBTObz5i6HdMcsmnI6E0XMAQTa8rlSmAkRUmz24AJhZs3aTYTSLuAsQkhqriQkBCzUaNGbiAAwG63a0BKyIQ4jZ07d3L8+HF/hyGEEEL4nSTKhDiDxOmTH1EYDfwdR3nRmijQoYWfGcO11lehjJWJ06d8ohzmw5vuGHPsAodUPGn2EFZFyECshQC8SbN1wCysKZoHLnB8QohKbPjw4fGAt2KsjtaayMhIT48ePXJcLlcAQEhIiPb2JAPpSybEyWrXrs3MmTP5+eef/R2KEEII4XeSKBPijIwBGn2Nv6MoPxpVsrBCaa1DleIO7bZtwFq10l88FCXN7gduwKo06w28CLwA/IBVafYBcMg/YQohKgut9RzgyuL79u7d65g6dWpEsV0GQogyPfroo/4OQQghhKg0JFEmxBlot/1GbO4Af8dRbrSKV4Y5F4gEUJABrDcMY8jGQf/a6d/gSnAB8ws3J9ANK2nWF0gCngfWYCXN5gDp/glTCCGEEEJUpN27d9tWrVrl6N+/f57D4fhDY3399dcBNWrU0K1bty4op/CEEBcZSZQJcQbJf3/oqL9jKE+JM56tqUxsGn0CpdJBD948ZPRyf8d1BnkUJc1GAD2BQVgVZ0nAFOA7rKTZZ8AJ/4QphPCX/v37n0hKSlIulyvUbrfr4OBgrbU2QKZaClGWe+65x98hiNNYsmRJwN/+9rewMx332muvZfbs2fO8VzLfuXOnbfbs2c577rknt0aNGub5jgNw4sQJ9eKLLwYPGDAgLzY21nOu569atcqxatUqx6OPPppTfP+KFSscAwcODO/Zs6crPDxcn81YU6dODW7Xrl1B+/btSyTERo8eXa1Nmzb5kigTQpRFEmVCXIK0AkOr0bHBDd+aO2DAOb+J8bMcrITYXCAC+AtWpVl34M/AW8Ciwuc/AbL9E6YQ4kJSSmEYRolNa+17TgghqpqYmBjPI4884nsfs3XrVvurr74aPHz48JzExES3d3+zZs3+0Hu53377zTZ+/PjQW265xfVHE2UZGRlq/PjxofHx8e7zSZQtWbIk4Lnnngs5OVHWunXrgldeeSXT6XSe9ViTJk0Kfeihh7JPTpSNGTMmq3bt2n/odQohLm6SKBPiEpNnC9wVWuBpsH7IyOMb/R3MH3cca1XMmVhTSf+MlTS7sfDrN4HFhc9/AZz33VYhhBBCiAupfv36nn/+85+53scLFy4MePXVV4N79OiRf/PNN7uKH3vgwAFjz549toiICLO0BJVpmuzYscN29OhRIzw83GzSpInHZrOdd2x79uyxHThwwAgJCdFNmjTxBAYGnlWVV2ZmpkpLS7PZ7XbO5bymTZt6mjZtmnvy/oMHDxp79uyxOZ1O3bBhQ0+1atXOON7tt9/uKm2/2+1m48aNdtM0adCggflHk4ZCiKpLmtsKcYnZfuf9meuHjLwY138/ipUQ6wU0AB4EfsGapvkRVuN/7/N/rLmFEEIIIUQlkJ6ebvTs2TOiTp060e3bt49s1qxZjWuuuSZy9+7dvizYpk2b7PHx8VFNmzat0aVLl+pxcXE1atSoEZ2enm58+eWXAd27d68OkJCQEKWUqqWUqpWbm1tqKe6hQ4eMa6+9NrJ+/fo1unTpUr158+ZRkZGR0WvXrnXs2LHDdsUVV0QD9OvXL8I71rp16+wAvXr1ioiKiqrZunXrqObNm0dFREREP/nkkyHesceNGxcyfvz40MzMTOU9NzExMQpg9uzZTqVUrYyMDAXgcrlU7969wy+//PLojh07Vm/ZsmVUVFRU9Jw5c5wAwcHBNTMyMtSECRNCvWPNmDHDCdCiRYuoYcOGVSv+ul5++eWg6Ojo6FatWkW1bds2qnbt2tGjRo0KLc+/KyFE1SGJsnP0r+596FMrpsxtfN87/B2iEAL2Y63e2R6oj5U0SwbuBOYBBylKmkllrRBCCCGqpN69e4dv2bLFvmLFiqMul+vQli1b0gsKCujXr1+495j77ruvWlRUlHno0KHfs7OzD2dmZh5+8803MwMDA3X37t3zP/nkk+MAq1evPpqenv57enr6706ns9TKrMcffzzk8OHDxvbt29Ozs7MP5+TkHP7444+P16hRw2zYsKFn06ZNRwDee++9DO9YV155pRugc+fO+Zs2bUrPzc09lJWVdXjcuHHZEyZMCF24cGEAwNixY3PGjBmTXa1aNe09d9WqVaX2Cn7jjTeClixZEvjDDz8czc3NPZyXl3fo+++/PxYbG+sG2Ldv3+9hYWF69OjR2d6xbrvttlIryaZPn+584IEHwu644468Q4cO/Z6Xl3fol19+OdKmTRvpYSbEJUo+IJ6jm4YO4tihw6fs3/D9SjauWI0zJKSUs4QQfrQXK2n2Elal2c3AwGJbOvApMAtYBZzVFAAhhBBCCH/6/vvvHatWrQpYuHDhsaSkpAKA2NhYz0svvXSiU6dOkZs3b7YnJia69+/fb3Tr1i2/Zs2aJkC1atX0rbfe6ksahYaGaoDw8HAdFRV12umGe/futcXExHhiYmI8AEFBQbpHjx6+1hYRERG+a5w81kMPPeTrOxYSEqLHjh2b/eGHHzoXLlwYcMMNN+QHBQXpoKAgrZTiLOIwIiMjzVatWhUAOBwO2rZt60tsRUZGaqUUQUFBZ3xNkydPDmnbtm3+q6++6lsMqkWLFu4WLVq4T3eeEOLiJYmyc9R5QJ9T9u3YsJm5L7xKWGR1/vbUOD9EJYQ4S7soSpolYPUzuw24u3Dbi7Vq5lwkaSaEEEKISmz16tUOgHnz5gV+8803Ad79eXl5CiA1NdWWmJjo7tu3r2vq1Kkha9eudfTu3dvVo0cPV6tWrc4rCdSvX7+8YcOGhcfGxkb17dvX1aNHD1dSUlLB2fQ7W79+vf3ZZ58N3rx5sz0jI8MAOHz4sNG0adNzbpbWp08f1yuvvBJcr169GrfccovrxhtvdHXt2rWgrEq4smRlZalt27bZJ0yYkHWuMQghLl4y9fIPOnLgEP++625M02Ts7GnUblDP3yEJIc5OMjARaAYkAk8ALuB+YAWwk6Lpm0IIIYQQlUp2dray2WzY7fYSySGn06lHjhyZXb9+fRNgypQpWZ988snxJk2aeF577bXgq6++OqpDhw7Vs7KyznlJ4KFDh+YtXbr0aFJSUsGHH37o7NixY2SLFi2iivdEK01aWpqtU6dOkUePHjXGjRuXPXPmzIwPP/wwo1mzZu78/PxzjqN9+/YF69atO9q/f3/X0qVLA3r16lW9UaNGNX744Ydz6kObm5urtNaEhYXJzVEhhI9UlP0BuVnZPP3Xv3H00O88+p9Xadqqpb9DEkKcn2SKEmcJFE3LvL9wS8GqMvsA2OKfEIUQQgghitSrV8/0eDyMGDEit7SVLovr27evq2/fvi6wKtD69OkT8fbbbwc99NBDOd5qMK3PLlfUsWPHgo4dOxYA/PTTT/Zu3bpVnzRpUvBbb711wjCMUseaP39+QH5+PvPnzz9utxd9BE1PTzfq1q3rmxppGMZZx9G8eXP3Sy+9dAJgx44dtptuuili1KhRocuXLz9WeIg+01hRUVFmWFiY3rJly/kvASqEuOhIRdl58rg9TB16L7tStnDX+FFc27N7uYz7268bGPPn/gxOaM3nr08rlzGFEOckGRgNXAFcB7wMRAITgFSKEmqN/RSfEKKScLvdLFq06Izb9u3b/R1qpXTvvffy4osv+juMP2z//v3ceOONpKenn9Xxbreb/v37s2bNmgqOTFzsevfu7QoNDdVjx44N9XhK5skOHjxoePft37+/xGe+bt265TudTp2XlwdArVq1TIB9+/ad8bPhyWNdffXV7ssuu8z0TveMjo42bTbbKcd5PB7ldrvxTrkEmDVrlnPv3r0lElQ1a9Y0s7OzlXd1y7ONIyYmxtOsWTO3y+XynVezZk3zf//732lfk2EY3HzzzXlz5swJOjlZVnwsIcSlRSrKztN/xj7B+mUr6HbHrdx8z9/LcdynaNW1I2NnDyIwyFlu4wohzpkJrCzcHgLaYvU0ux0raTaBokqzGVhTNYUQl5Ds7Gy6det2xuNGjRrF5MmTL0BEVcvatWtxu8unV3ZsbCxjx45l0KBB5TJecZs3b+bmm2/m888/p3nz5qc8P2bMGKKioqhRo8ZZjWe322nevDkjR45k9erVKCWfxcX5qVmzpvnOO+9kDBo0KLx58+ZRnTt3zjcMg23bttmWLl0acPTo0d9DQ0N1XFxcjdatWxckJCS4nU6n/vbbbwNDQkL07bff7gKIi4tzJyQkuO+8887wVq1aFTidTj137twMh+PUWYw33XRTRFBQEFdddVVBRESEXrVqlWPnzp22d955JxOspvo33XSTa+zYsdW+/PLLwJCQEP3CCy9k9evXz/XEE0+Etm/fvnqvXr1cu3btsq1YsSIgISHBfdL4rpCQEN2yZcuo+Ph4d0xMjMdbNVbcww8/XG3Tpk32du3aFURHR5spKSm2L7/80jl9+vQM7zG9e/d2vfjiiyGpqan2yMhIc+TIkTneSrjinn/++ayff/7Zcc0110TdfPPNeTVr1jRTU1PtNWrUMGfNmpVZDn9VQlRZulMnO4cPW43aTXOf2rLlhxLPN2vWFLu9KaYZhNY5KjX1K78EWs4kUXYePn3lLb6dPoeWnTswfOpTpR5TkJ/P0YOHqFXvCnZuTiHzyFH+1NFqdZSdkcmu5FScISE0iG+GzWFHa83hPfs4kLaT2kMHkp2RiTMk+EK+LCFE2TyUTJp1BgZhraA5AXgc+AErafYhcNA/YQohLqTw8PASU4Ty8vIICgri8ccf58knn/RjZJeenTt3kpGRceYDz4PL5SItLQ2Xy3XKc2lpacyZM4cffvihlDPLds899/D000+zcOFCuncvn1kJ4uLXpEkTz5NPPpkVHx/vSy4NGDDA1bJlyyPvvvtu0NatW22BgYFcddVV7vHjx2eHhIRogP/+97/Hly5dGrB3716bUoo+ffrkjRgxIrd27domgM1mY82aNUc/++yzwH379hlut1uVlcB9/fXXT3z99dcBaWlp9vT0dNq0aVMwbdq0zKZNm/pK2j799NPjn3/+eeCOHTts+fn5Kjg4WNeuXdtcu3btkZdffjl427ZttiZNmnief/75I8uXLw9wOBy+H6R16tQxU1NT07/++uvAgwcPGtHR0SZAy5Yt3U8++WRWYGAgAOPGjcueN29eQHJysn3r1q22K664wly9evXRNm3a+BJhU6dOzeratWv+pk2b7Lm5uSoqKkoD3HvvvTl169b1xRsdHW3+9NNPR6dNmxa0evVqx759+2zNmzd3Dxw4MK98/uaEuHB0QsJDaH2vb4dSxzBNF0qlodQaHI45av3642c94IEDQdhsHwFgGHOxPvOgQREfPxO4E9P0XmsPUF8nJNyCd/ZNcvJzyvosVaVIouwc/fDlt8ye9Bz1YpvwyNsvY7OXPp195+YUHu9zB9f27M7ulC0Eh1XjTx3b8817s5n9zLM0afknMtKPoAyDcXP+Q7XI6rzxyFhysrKY99a7hH7wCfc89wy16l9xgV+hEOIMPMCiws0JdMOqNOsDJAEvUJQ0+y9w2D9hCiH8bdOmTXzzzTc8/PDDeHsALViwgPXr1zNs2DAiIyMB+PHHH1mzZg0PPPCA79x169bxwQcfcOTIEerVq8eQIUOoX7/+Ga9pmibz589n6dKlnDhxgho1atC+fXt69eoFWEmfZcuWsXz5cg4dOkTt2rXp1KkT119/fYlx3njjDRISEoiOjubdd98lMzOTG264gVtuuQW3283MmTNZs2YNderU4R//+AfR0dG+cz/66CMcDgdXX301b731FgcOHKBFixbcfffdBAUFnTb+goICPvzwQ1atWkVBQQFt27Zl4MCBBAQElHnOlClTME2TRYsWkZOTA8Att9xC48bWe/S9e/cyY8YM0tLSiIiIoH///rRt27bEGKtWrWLevHmkp6cTFhZGq1at6N+/P+np6cyYMQOAGTNmsHjxYgDuu+8+goODmTZtGjExMbRu3brEeHv27GHmzJns3LmToKAgmjRpQv/+/bn88ssBiI6Opnv37rz99tuSKBNnrWHDhp7HH388++T9TZo08fz73/8uc9XGnj175vfs2TP/dGOHhobqs0kMJSUlFSQlJZ1SlVWc3W6nX79+p2SWY2NjPa+99lqJ6rA77rjjlGvWqVPHHDZsWG7xfQkJCe7i1WeJiYnuxMTEM5akdu/ePb979+4lXvvw4cNzTz4uJCREP/jggzkPPvjgmYYUorKrDjTyPdIarMR3O7S+k/z8J3Vc3ECVmvr1WY0WGGjidqcVjlX0uSYxsS2meWfho1/Q+lPgaOFxdwHWG4+EhP8jObnKJcqkR9k5cBcU8OK9D6NNk4Ska/l2xvt8+spbJbYVn833HZ+fl8cVsY15cdnXTJr3Ift+287sZ57lmS8+YOLcmbyw5Esa/6k5syc9hyMggIlzZxJSrRp3TRjDxLkzJUkmROWXB8zHqi6rCfwFmAO0BF4E/odVhXY3EOanGIUQfmKaJqNGjeLnn3/27Xv88ccZPXo0CxYs8O178cUXmTdvnu/xf/7zH1q3bs3KlSsJCAjgv//9L4mJiaxateqM1xw5ciR33HEHR44coXr16uzevZt77y26sfz9999z1113sWPHDoKDg/npp5+44YYbeOqpkhXykydP5umnn6Zbt24cOHCAnTt30q9fP95880369+/Pf/7zH9xuN2+++Sbt2rUrUW31zjvvMGnSJK699lp27tyJaZqMGzeOrl27kp9f9mf1rKwsOnTowD//+U9ycnJQSvHoo49y/fXXl3me1pp169ahtWbv3r2sW7eOdevWcfy4dbN86dKlxMfH8/HHHxMUFMTWrVtJSkrirbfe8o3x7rvv0qFDB7Zt20ZkZCRHjx7lscceY+/evWRlZZGamgrAli1bfON7p4x+9dVXdO7cuURMW7ZsITExkfnz5xMeHk5BQQGzZ8/ms88+K3Fc586dWbBgAQUFp805CCGEEOdH62loPQh4AthcuDcSpT7ViYlXl3la48aBOjGxhU5MvEJt3JitUlJiVEpKjEpN/afvII8nxve1Us+r1NRnVGrqGxX0Si44qSg7B6bHQ35h08tv3p1V6jHNk67luj69fI9v+ltRr4y13yyibtMmeNwedmyw/p3G/CmRT166aP49CXEpy8VKms0H7sOaltkfuAGr0uxl4DusSrNPgTLv/AohLg4tWrQgOjqaxYsX06ZNG44dO8Yvv/xCmzZtWLx4Mbfffjtaa5YsWcL9998PwO+//859993H7bffzqxZs1BKkZ2dzXXXXcff//53UlJSyuxppbXmvffeY8qUKSWSY8WTTG3btmXPnj0lKrSef/55xo0bx8iRIwkNDfXtX7NmDRs3bqRBgwYA9OvXjwceeIC7777bl/RJTk6mefPmfPHFFwwYMMB37s8//8wnn3xC3759ARg+fDjt2rXjzTff9L3Wk40fP57t27ezceNG6tWrB8Do0aNp0aIF06ZNK/GavJRSfPTRRwQEBDBkyBDuu+8+33N5eXn89a9/5S9/+QuzZs3CuxrfpEmTGDlyJP369SMqKop33nmHoUOH8vbbb/vONU0T0zSx2+1MnjyZq6++mn//+99cfXXR54rMzEw2b958SlwfffQRNWrUYPXq1b5KwpP/HgBatmxJdnY2GzdupFWrVqV+T4QQQojzptQvKiVlFoDu1Olpfv/9DbT+OxCIx/MM0B1Ax8d/W/i1xjBuxDTnYJo1gDd1bOyj2GzeXn1zVUowIg+dAAAgAElEQVTKAB0XtwClbvBdR+s5Oj5+DvAtcGOJGLR26fh4UGqTSk5uUdEvubxIRdk5CHA6+ezQjtNuT346p8TxwdWK3nAeOXiIY4cO8+mrb/m2jSt/IK5NmclcIUTVlAHMxCo5rg3chZUkuxGr8f9hiirRQvwUoxCigiml6NKli2+63rJly4iOjub+++/37du0aROHDx+ma9euAKxcuZK8vDwee+wxX0IsJCSEBx98kC1btrBnz57TXi8oKIhFixaxbds23/7iSbFq1aoREBDAjh07WLZsGXPnzsVut5OXl8dvv/1WYrw+ffr4kmQAHTt2JD8/v8QU0YSEBGrUqFHiegCNGzf2JckArr32Wjp06MDChQvLjH/OnDkMGjTIlyQDiImJoVu3bqxcubLM88qyaNEiDh48yJgxY3xJMoARI0aQm5vrq/QLDQ3l559/Zu3atb5jDMPAbj/9/eSDBw+itS4x7RSsv6/Dhw8zb968EpV2J08frVmzJmCtmimEEEJUJLVsmRubbWzRDtVZ1617cj8EhWnOBWr4HpcuCyg+DTwbOFa4/xhQvFT6GHAMravUwhhSUVaBTr7jW/OKOlSvVZN/TXvFTxEJIfzgGFbSbCYQBfQEBgI3AX8G3gC+AmYBC4DT9hARQlQtXbp04YEHHiA3N5clS5bQpUsXunbtyp133sn27dtZvHgx4eHhvkqltLQ0DMOgSZMmJcaJjY0FYMeOHdSvX59BgwaRkpLie/6FF16gQ4cOvPXWWwwdOpTY2Fji4+Pp0aMH999/vy/5dODAAfr168ePP/5Iw4YNCQ8P91U6HTxYch0Sbz8tr2rVqgFQp06dU/afOFFyUTpvf7DimjRpUub00czMTA4fPsybb77Je++9V+K5nJwcrrzyylLPO50dO3YAcN1115Vahbdr1y7Ammbav39/2rRpwxVXXMGNN97IiBEjuOqqq047vjcJdnICbNiwYSxcuJC+ffsSHh5O586dueOOO+jXr1+J47znlbZIgBBCCFHe1MaNh3V8/GGsljEOQkNrAbtOOuxn7Pa7cDgOkZd3WanjpKbeohMSBqO19Qtb67+r1NQPvM/r+Ph5eHuUKVVbJSdXuc83kii7gNr0uIH3J/8f306fQ/e7/mpNp8jIZGdyKont2vg7PCFExTtCUdKsLnAL1vTMfoV/HgO+xJqe+S0l78YIIaqgrl27kpeXx6pVq1i8eDGPPPIItWrVIjExkcWLF7NkyRI6duzoq14KDQ3FNE3y8vJwOBy+cbxN6r1TI/v06UP79u19z3sb/ffu3ZubbrqJ5cuXs2TJEmbPns27777L5s2bufzyyxk5ciQZGRns37+fWrVqAdb0ycTExBIreAIlqrCKK2vqZ3HeeIvLzc0lJKT0Ilrva73vvvu4++67T3neu9LdufB+TxctWkT16tVPeb5GDeuGecuWLdm2bRs//fQTS5Ys4YMPPuDdd9/l+++/JykpqczxvecfO3asxP6wsDAWLFhAWloa3333HfPnz2fAgAFMmDCBCRMm+I47etTqeXxyRZoQQghRETQYKFUN7+97pU69U6P1A2rjxn2Fj/bo2NhqFy7CykMSZRXEZrMTEl6yd3ftBvV4bNbbTHvsCWZPeg5ncBC5Wdn8ZcRQX6IsJDzsjKX+QoiLwj7gpcKtHtaqmf2xqs0GYiXVvsZKmn1NFVxWWQhhTR1s0KABM2fOJDU11TfFsmvXrnz77bd8//33JRrpN23aFLBWYbzxxqI2HytXrsRms/kqzfr06VPmNQMCArj++uu5/vrrGTZsGI0aNWLp0qXccccd/Prrr/Tq1cuXJPNeq7xt2rSJ7OxsX2LMNE3WrFlzymqTXkFBQSQkJPDLL7/QqFGjUo85HafTeUpTfG+V3m+//cZtt9122vMNw6BNmza0adOGhx9+mLp16/LFF1+QlJSE0+kETu0xVrt2bWrWrMnWrVtLHbNRo0YMHz6c4cOHc+utt/LZZ5+VSJRt3boVwzBo3rz5Ob9eIYQQ4pzFxbVGa+90y3RSUw+dckxwcOm/1C4x0qOsgsT8KZF3Nqw+Zf+fOiTx6sqFvP3zcqZ88xmztv7CrY8UNbV9fc0S6VkmxKVnD1bCrD3Wcs6jsfqYDQTmAbuLPX/mUg4hRKXStWtX3n//fRo3buyr/OratSvz58/nxIkTvuQZWNMEmzVrxsMPP8z27dsBa+XGF154gf79+5daGeV1/Phxnn/+ed8USq01K1asAKyEHVjJm2+++YbDh60V3levXs3TTz9d7q/5xIkTPPjgg7hcLgoKChg/fjw7duxg6NChZZ4zduxYvvvuO8aNG0dWlrXeSW5uLvPnz+fbb7897fUaN27MggUL2Lp1K2lpaeTk5NCmTRu6devGyJEj+e677zBNE7B6gj3//PO+lTGfffbZEj3W1q9fT2Zmpu97Vq9ePQICAvj000/ZsWMHaWlpeDwelFJ07tyZH374oUQs77//PitXrvStjHn48GFSU1NPSQD+8MMPtGjRgqioqLP5lgohhBDnTcfHt0Gpmb4dSs1WYJ58mFq3rjxmtBTd4M/IsJ3muEpLEmV+EhxWjajLamHYquS/GyFExdkJTAHigUSs5ZyzgfuBFVh9BLxJMyFEFdClSxc8Hk+JhFiHDh1QSlG7dm0SEhJ8++12O59//rmvT1loaChdunShXbt2vPHG6VfJNk2T5557jssuu4zq1atTvXp1RowYwVNPPcW1114LwJQpU8jMzKROnTpUr16dnj17Mn78+HJ/zT169ODIkSOEhYURERHB1KlTee655+jUqVOZ59x+++288cYbvP7664SFhREZGUlwcDBDhgwhM/P0PYCfffZZ0tLSaNasGTExMSxatAiwVqDs0KED3bt3Jzg4mJCQEOrWrcvMmUWfFT7++GNiY2MJDg6mZs2atG3blttvv92X1KtWrRovvPAC06dPp3HjxsTExPimTQ4ePJjly5dz6FDRTfnNmzfTqVMngoKCqFWrFnXq1CEwMJCXXnrJd4zb7ebzzz9n8ODB5/qtFUIIIc7WAzoubqlOSEgDfgC8DVCTyc2dWGFX1bpolZqwsNd1QsIonZh4w2nOqHRkjp8QQlReyYXbRCABa2rmnVhJs/uxkmbzsFbS/MUvEQohfAIDAzl69Khvqp7XgAED6NGjB0FBRYtLhYWFkZ6ejmEYp/T8io2NZePGjb4VMWNiYmjYsOEZrx8ZGcmBAwfYvn07+/fvJzAwkLi4OMLDw33HtGjRgq1bt7Jp0yYKCgpo1aoVgYGB3HLLLb7+ZwAbN248pUn9X//6V26++eZT+oX98ssvpxzrcDj45JNP2Lp1K/v37ycxMdG3yqPX4sWLT2k3MWLECIYMGcKGDRvIysqiVq1axMbGnrEtRdeuXfntt99wuVzk5OT4XktERAQffvghL7/8MqmpqTgcDho0aFBiQYIff/yR/fv3s3PnTsBadKD41FSAe++9l3vvvZfs7Gzy8/OJiIgA4IYbbiA2NpYZM2bw6KOPAjBp0iRGjx7Nli1byM7Opm7duqcszvDll1+Sn5/PXXfdddrXJYQQQvwBzVCqGUU9SAvQ+n0CAx9UKSkZFXZVw5iB1sOAAGBw4fXfAMpe+rqSkUSZEEJUDcWTZq2AQViJM2/SLAWrn9n7wLbShxCi8kuYObl98qDRK/0dx/lQSpU6NdJut5e6v3gCq7SxWrRocV5xNG7cuNRVJ72Cg4Np06bkIkInx1dabAEBAackxMo61is2Nta3YufJwsLCSt0fGBhI69atyxzzdAIDA0tt/F+rVq1Tkl/F1alT55TVPEsTEhJSYkECwzB46aWXuO2227jnnnt8K4OGhYWV+Rq01kyaNInx48f7Em7C764AbJy6+psQQlQtWi9EqSzfY9M0MYzjKLWbvLxf1fbtv5dy1rsotbTU8ex2F1ZbGNB6i2+/x7MOm83abxgbip+ikpN/0s2bx+F2d8EwrP4CSv36R17WhSaJMiGEqHrWFW4jgXZYCbNbgQmFmzdpNgvY4acYhTgvylTLEmZMWeQwzRHrh4zZ5e94hDiTrl27snbt2rNejMk0TT744APq1atXwZGJcxCDVenwMtbv0Wz/hiOEEOdHpaSswGrXci7nfFTmc8nJ+VhtYUru37JlE7CpzPM2bUoD0s4ljspEepQJIUTVZQIrgQeAOsB1WG/ya2K90d+OVYU2CrjcTzEKcU405BhaXe9RxvrmM6a8HfvOlEtyWfKqbOzYsTz88MP+DuOCatiwYYmptadjs9lo1KiRrHJe+WQD92JVlY1APicJIcQlS35Di0uERqNBg2maqsDjxqNN3KbJ16nrQ/PcBfJmSFR1Hqyk2UrgUeAGrEqz3sBkYBJWE8+5wAfAqctBiyrnTzMmNfWY6uJKgmpsWmkbEK5NBjrs9PnTzCmPb0jLfZuJE09enUlUQh06dPB3CFXOu4MGUZCXV2Jfz3HjqHue028vhIMHD0YAnfwdRzm5EmvqpbNwm4JVtT0EOHUZeyHEpWA2cBz450n7ve+zewDlsUKkqIQkUSYuDbpwAzxo7TZN09Sm8pgmP+3Z+TuAYWCijRw/RilEeXEB8ws3J9ANK2nWB0gCngfWUNTTrLReBReVK9/7v4j1Q0Ye93cc5c1j2h9C6eH+jqOiKIUTjdODfqV5w+Abm3300S1zBwzwnPlMISqXcau+Jjn9YJnPN137I0Z+yc9b36WsZ0glTpT98ssvV2H1y7wYhRVuK7H6gL7q33CEEH4wHtgALAY+K9zXAmtFekmSXeQkUSYuFRqlNEopA5UX7AjINFAhSmnP+kH/OvNSYkJUXXkUJc2CgOuBgcDNFCXNlmL1M/scyPRPmBXLrVyfJb43uY5W/BdtfJY85NH1/o6pfHg+UMrY7O8oypPWPIuV4AXrFkeWQn2Qa3M+IkkyUVXtP5HB9uPpRAeFnva4/3W5FleEtchBXN3LLkRo5y0mJmbbvn373vJ3HOWkCfB3ILjwcS7WVMwJwJv+CkoI4VdpwDjgNeB74AQwHXgbWIaVTL8Va2X6A8DMwj+9OmEl1EKBfcBHnLl3cAjwCtYCI6FY74fCCv8MxKpw/csfe1nibEiiTFy0ihWRYWqNqa3pl9n5LmPv8SOOTJfLwJquJsSlIpeipFkE1i/a/kB3rARaHrAIq9LsUyCr9GGqHq3UCgVjFYxTSj+QOH2qR8EClDnrhJm3dNeQiXlnHqXy2TxkzDKsN2sXjcT3pkxC4QQyQG8A42+b7npUFqUQVZ7dMFgy4B+lPjf+idcooIBJDz5GncTECxzZ+YmNjd33/fffXyyVVp2AO7FuKJ0APgQe4SK9eSSEOGuvYM3IeBXYipW8GgNEYyWt0oAFWMmy9cC1wE6gH1aS/SVgL1YyvjNnTpRlY73/vg4rMVZcBqU01RcVQxJl4uKltbUBLrdb5ea7lMfjYevhA4FbDx8IBFBKHfNrjEL4z3GsO18zgUjgz1hJsxsLv34Tq9R8LvAxUKWnJSvDWILHvB8I1xBu9S3kdrQxIBTnq8CD/o5RWJRSTq30boXx9013/WuRv+MRQlwyIrBuPAxDVowWQlhM4G9YUzBvwUp25WD1/t2N9b7ZW5sxDXgM62dId2AO8NR5XHMUVqLt5LLiZGDVeYwnzoMkysTFSykTpQCMALstMzggoAAI1TAfrZ4BMAzl9muMQlQORylKmtUA+mL1nemJlTR7CasKbS7wLVWgJ0P/jz6yJeekNbNhJGmlr8fUnYHwkw47AfrjzYNHjWTIaH+EKUqhtb6/Rv3c/yzrPFF+PgshLpSDWB9sJTkvhDhZGlaVaSOKFvdoAxzDauzvFQF4W/qsAN7Amjb5BdbPlrOdqdEEq3qseKIsgz94U1fHx3+H9RoOqJSU9n9krKpOx8dvxqog/lWlpPTz7Y+L64lh1AFJlIlLhvIYynADpoL0zRdNfyIhyl06Vu+Ft4H6wACs/gsDC7d04BPgYazy8Eqh1VsTg/Mcgdcqm3EdWrdPzdnZ1kCFFK526wZ+BW0DVR1AozM1anrK4FEPMHiUn6MXxW0eMkr6AQkhLrQthZsQQpSmgJI3iqtTsh8ZwM8UJdtnYlWc3YY1bTMM6yb056e5hhMYjTW10wC2A/UAB7AW+OkPvQKlrkDrRoXjnTMdHz8P6FU4VqBKTs7/Q/Gc7lpxceNQ6qnCa12nkpNXlvMlGgAhKHWoxF6lHkTr60ESZeKiYVW86sI+ZBrILchXeW7r59ncDT+G2JUtKDs/L8CPQQpR1ewGni3c6gO9saZn3opV/v2KvwJrPPvlsCCPq7U2zfYoneRCXacg0PpRoHKAX7VmJehVDgJXrB8y8njijKkvo/U/QZ9QqKnJg0c946/4hRBCCCFElbUVOMzpe4Z9X7jdW3jceMpOlF2HdZO6GdY0z2FYlWxbgQCkRcgFJ4kycVFR1sqWgMKmVI5hZeODcvPz56L1bgCt1I/+jVKIKmk31hTMl7D6m/2FC5goi31nyuUBhkrSSrcHlYQ7t6UGAwWgToBegVarlGGszM0IXLH9/vtdpQyzUGvuxVBjku8a9fqFil0IIYQQQlxUXsNqS/IjVmsSE7gKiMKaankn1k3lZKwKrtrA/lLGiQCeAP4JuAq/foai6rUHsVaqTzmX4HRCQmvgIaA2Wm9AqX+XeWyrVsHk5AxGqW6F8R8DvkGpd1Vycr5u0MBJcPBLQIvir1/Hx5sotV8lJz8JoFu0CMHjGVxYkRUFHEWpr4HpJ1ef6djYWOz2oWjdAghGqf8Bi1Vy8jQdHz8RuKnoYP0vHR8/EACl7lPJyfkaDOLj+xV+b64AstF6FYbxukpOPlriWnFx3YARGEYUWq8mL++sFkSQRJm4aBS27tdoFGgC7Y4sh81mA4K04p3kwaOX+jtGIS4S+syH/DHN3/13I9NQ7UElKWtFzkbad1l9AMUnoFYp7Vm5aZfrVyZONM80prJ7Vul8dUfyXaM/qNjohRBCCFGadVsOOX7ZetgR3zDSndSiToVN3RKinKUDwcUeL8SaYTEWeB1QWNVf3tkKVwBPA3Wxkl5LgBEnjdkLq49ZHazKs7uBbScdMxv4+lwC1c2adUTrBRStmtkRrXtz6iqa6MTEWuTmLkGp+JOe+gta36nr1u1GYKC9MLZiJ+q/F361CXhSJyTUxu1eAsSddNzNwJ26bt3uat++XAAdFzcIpaahdUCx40CpLlgLIvQHisfzF99XTudI3aqVJifnU6w+ykWUuhGth+lmzTqqLVt2Aej4+P7AB4BRuMhfRwIDe3IWebCzSpTpiv9MJMQfogHTNFWBx41Hm7hNk69T14fmud2Gv2MTQpxep6UT7Ud2B/5JK1t70EloOmuooYoOSUMxS2m10oO5KmXw6OTzuc6mO8Ycw/plKYQQQogLKCu3QI17c2XY4p/3OL37WjWrlf/6v64/Huy0y4dNUdk9Xsq+zws3G+A56bl/F26lPXcZ1qyMW7Cqt4ZjJYjK+n9wtIz9pTOMV7GSYhqtH8dmW4XW96N1n1OONc2XsZJSJvB3lPqiMKn2DpBEePhjOJ1Pkpc3HNO8F6WsqjKt78Uw3L7YtH4FiEMpDzAUmI9p9kWpacB1hIWNBibo2NiGhfsCUCoXGAesxKp8610Y1UTgDqxqMVDqWax+bVCtWj6HDo1EKW+S7Dny86ficMSh1JdAPQzjDaCHTkgIKIzLAPJR6mEgBaXGo/UpScNCx7x/njZR5kuQyY8uUQWYaO02TdPUpvKYJj/t2fk7gGFgujFy/B2fEMLSYuazIR48LQ2tkrSmffpurgPCC+/0uEFtAP2+0mqlK9BYuu2vj6T7OWQhhBBCnKc8l0f9Y+qiiF+3HQ6Iqx9VcNsNsTlfrUoLWptyMGDkS0vD33y023GlzjyOEJXUyYmwsp5TWL3HnsVq7j8Xa8rl4fIKRDduHA0kFj5coVJTnwHQV165nvz8nlj9zqxjGzRw4k1GwRa03ovWVwF7UOogVvLqVrVu3ePA2zo+/s94p18axn+80yl13bpBeKu+tE5F6/1Y01B3AYeAWih1KzABw+jvi0HrKSol5YVi4c8DUCkpc3VcXCxKeWObV7yZv46Pv63wSxdKLcLh+FPh41+BDsANunnz6phmfaBW4bU+USkprwLoZs32YRhby/gW7vb+eYaKMqVRkiYTVYFSCpUX7AjINFAhSmnP+kH/anjm84QQFc3beN80zeuVor1pmtcoVEDhL5dsYL238b4zn+Xrho/K8GvAQgghhCg3Y99aGfbrtsMBSS0ud736SNfjdptBr/YxriFPfxuxeuP/Ar9Yvt3Zu2PjPH/HKUQFSsRq1t8WK4F0K/BtuV8lMLBB4Y1n0NqXDFLr1x/X8fEHsBbnsoSE1CxWWRWPUt+VMmL9UvaVFBpai6IEXGKp42htjaPUFb59Sq0549ilsVbvBAhE6285NctuUFBQF6Ua+J4zjKJVjbds+Y34eDelza5Uanfh2KUnyjRgasmPiapDo8nJzzP2HEt3nHC5DE6f2RdCVKBSGu9fpa1sNqAzgeVKqUWm0qtw5q1NHjBRepQUoykq5Da19lV3/+/4MbvT4bD6MMqddyGqvJ8/+ogN8+cD4HaVtv6IEFXfvBU7nAt/3OVsWq+6+8WHumTYbVZXFIfd0E8Ma3ei35j5kS9/9Eu1Hu0augIdNvkAKi42Dqym+k8Ufv028AhwokKuZhhH8BR+DFYqwrvbWgBLRVAyx5Pl+0qpTcCrpYx45v+TppmFUdjtSOuNGMZrpRzj7SVc/HXXPuPYpdH6BBANZKLUv8qI6QB2e0Sxx0VfN25cjbJbkE3HZvua3NwTpRygcXtM5fZ45G24qDI8HpMthw8Ebjl8IBBAKXXsTOcIIcpXwvTJLygYANTxJnc0epdSzNamXqFsauXmgaO2oJS8ET6dwjcxpulRufku5XK7AZi1blV1f4YlhCjbjuNHmLbpB+uBhhPH3FSLLHqbnZZx5JRz0nftYvuqVRcqRHF6VwHtsZpyu/0cy0Xj9+M5xr9n/lgtwG7Tk+65LsN5UiIspk6E+y/XxeR9uuy3oM+WbXfe1i0211+xClEBkrB6j8UBG7GmXa6t0Ctu2rSbhIQMtA4HbtSJiTFq8+YdJCQMLNzno5KTj+q4uI0o1QKtr8Dj+U5t3brT+7xOTGyH1pcXnaCOF0u0Xel9LWrbtnQdH78Zq5rsCjyehd5m+gA6ISEJw7CSYlovQakxhV8/qmNjv1Nbt/5PgyIh4QaVnLwAAMMoupZpttIJCWuLrZy5DGgEhGGaO1Vqqq+CTcfFNQGuV9u2peumTVOx263KMaX668TEKWrz5kM4HPed5js4FdNsTWDgqZk0p92RqdGSJBNVioJQ4Avt0ZMATJuWNzlCXGAKIxq0y9t4Xxnqu42D/rWzxEGDRvkpuipDo1RhzZhyBwcEZDpstgCgGtrsrE3rJoDHbu7zc5xCiGLSc7P4Ki3FerAnGPYGQ7t0ild/2o3S1xe64sorufJmqxVLxOWXl3qMqHBhWH2D/oXViHqhf8O5OLz0wa+hWTkFxj/7tcyKrVe91PfmQ3slZn+xYnvQ7AUpwZIoExeJcOBJrP5jLqxqsklAhc+gUODRWv8fVkP8aphmSuGUy7pAHuA86ZRHUOobtI7AZtus4+OXo5QNrWMxzXoo9QbwceGxa4GBAGj9o46PB6XeUMnJ/0CpR4Cv0Lo6hrFZx8evwKraaorW9bCq1T5RqamLdHz8F1i90eKx2dJ0fPwOrOoyD1Cz2LUKX5R6Ea1f1AkJ61Ry8tV4PE9js90MRKHUtzo+fiVaZ2IYTdC6KfAz8Ibati1dx8W9h1LDgLpovVPHxx8BolHKg9a2U76BWjcBWkEpJWdKKY/M6RBVjVaYoI4kDx21wd+xCHGp2rwr5y4mTjTPfKQ4GwqlDWW4FcoG4NA6+dehY373d1xCiLINbtqGj5bvI7/A5P7aPbi5U4zvubLeXUfHxNDurrsuTIDidLKwPkx+DGwBBgMp/gyoKkvZdcQ+f9WOoMtrhHoG/zmhzEW16tUO87RNvNy1csP+wHVbDjlaNatVcCHjFKKc9cKqTK0LLMda0XLLac8obykpTxEfHwbci7X6ZShaDwbuBxoDvl7AKjX1O52YeD2m+QyG0Qatb0RrUCoXrRdgml/5xnW5phEQ0ADogbVyJ2idA6CSkxfohIRuKPUMWrcBbrQuUDgOfO0bR6kBaD0G63tzGVZ/NA9aL/Mdkpy8VsfF3YtSw1CqDlrbMc0TAGrr1p06Lq5t4WqYNwIdCnuReYAf0XqW71q5ufcTFBSAYdyJ1kEopQpX/5wGBHvHLCaLwpUvfYmyLzavC/k6dcPJGUYhqoTs/LyAMx8lhKhQkiT7wzzaJLcgX3k8ptqfcdQ+dclX4QUe96l3vIQQfpVdkM+uzKO+x7szrY4Pv6zKwGNq3B6TN+duok/7pkSElrUKfdXzxBNPBO/bt89x5MgRTpw40RnrA+HF4HIguPDrasDVwBpgJvAYkOmnuKqsZ2f/VM00NQ/9tdWJM/Ue69upSe7KDfsDP1n6W5AkykQVdRnwMtAPOI6VBJrG2fT4KmcKTFJSHtYtWoynoKAGtWrtV8uWuYHZpR6/efMyIEnXrRtE9eqX4XJlqm3bTllxXm3f7sLqr/ZIqeMkJy8F2ulWrYLJz69d5jjWFMongCd0QkJtPB4nWVmH1L59JSpKVWrq65TxO0alpv4G9NYJCQGY5mVAPikphxSU+Cyidu3KAwbrRo0eIDg4jM2b/6eshFrdUsdNSenj/dpe7URWjivAOSmvIJ+8AumnLKoubb2hEUKIKktrPDalXECgqfWBE3k5vrtiRkF+th9DE0IUk3LkIH9b8EHJnVl2Nm84jmlan4s8psnrc9fz2JA2foiwYsydO9fpdvtm0DUv3Ch6YAIAACAASURBVC5GCgjCmob5MzDdr9FUMWtTDjp+Tj0U8KfG0QXd2zQ440oVnVtd4YoMd5qLf9rjzB7sPhESZJdepqKqUFjTEf8PiAS+BEYA+/0ZFIDauDEbyCY19eyO37cvl3370v7wddetywHOahyVnHzwD13LSrrtPuNxaWkZFKukOxv2dcMn5gBjzzM2IcT/s3ff4VFV6QPHv2dKZiYhhBR6SSYEFVFRwUIRwbqIFcWyP+wKltV1FSGwu8LaAMUtdnF3sbCKoq697KoUKQqCImVBQhICBCEQ0jPJzNzz++PMhBASQpkwgbyf58kjc++5975JdjN33vue9wghRITYFAGXw1mmFHYUuatuzpT3ZyGasb7tu3Jyu84AfPxWPgW1WtBUVQeZPednrj7/WDK6tGnoFEeUefPm7VqwYEHsqlWrXH6//+VZs2b9IdoxRUh/4J9AeNGUEuAr4C7gkD7ItUQv/ntFK4A7hvcua2wsgMNu46J+6ZUzP18T9+mibPeIc4+RXmXiSJCBWcVyCJAPDAf+HdWIRMQ0tCymEEIIIeqlm6yOfltpsePT//3YaltJkbw/CxFFhb4Kvtva8EPqDaFVLPt39nL7iWcyZ9kmXtuWR9Da86+DBh7+x2Jemzh097GLFrH1f/9j84ojr61qcnKybt26tfZ4PHg8nkpge7RjipAizOeiUmATcAOwLKoRHaGWrd3mXLrml5gT0lP8A3t33u/pSlcO6eGb+fmauA/mZ3kkUSaaOSdwP2b6oBOTLBuD+fshjhJyIy6EEEIciCacEFJYUV61NC+nAMDSutFSciFE01i/q4Cx8z/ar7H+gMVjr35LtT+4175g0OJ/uTuZ98Nmzj7FtERZ+dlnfPevf0U0XhERVcB9wBtEoa/Q0eLlD36KAxh9xUn7VU0WltGlTaCXN8X/04YC56Ztpfau7eP3/j+UaDJ33nlnot/vPzpKXyPI4XBUv/jii7WnUfbHJMZ6AauA25H2P0clSZQJIYQQ+0upplywwKa1/vjHG8dc24TXEEIcgPNTj+Gszt0b3N8rpQOvf76GbTsbXNQPX3WQqTOX0P/ETjgdtprtiV260C4jg449e0Y0ZnFQfgJSgYZ/kaJR2VuK7YtXbnX16JoYGHxq1wNufn1Rf69vdc4O54ffbHDffdXJ0pfzMAoGg2NsNtuEaMfR3FiWtQ44LvSyPWZKNsAfgCcAWXziKCWJMiGEEKIR+cW7lMvhbNIKA60xLWGFEM1Gr+QOXNFj3/3qO5/XhquGHFPz+uZHP+fnvF1MzzyfXukpNdvttj3/D97zvPO45KGHIhuwOFiFjQ8RjXn98zWxltb8+oLjDirhOGyA1/eXWcviP12cI4ky0RxtA0YB3wE/RzkW0cQkUSaEEEI04rWl3ygkjSWEqEec27nH63BCLM7jpHVcTDRCEuKwK6v0q88W5XgS4l3WxWd19x3MOZITPNYZvTpULfwp37Uya4fzxIwUqdY5zI455pjq22677YBWBzwaffvtt+733nsvvp5dr9ezTRyFJFEmhIi0ROC8Wq/9mMa4K4BAvUeIxryGWdL4nmgH0tJUeVo94yovOmwNnQM68OPhupYQQggRKW99uc5T7vOrW84/scLttB90BfZF/dN9C3/Kd324YINbEmWHn91u1wkJCU3ZZuKI4PF4pE9hCyeJMiFEpHmBt4FyoBqzGkwrIBe4ApBEwIGLR5KMUfHz1XduAWZHOw4hhBCiubK05p2v1sU67DauO//YQ+rzdsHpaVWPv/Kd/uK7XPe4608rddhtjR8khBARJn95hBBN5RYgCZPk6YepNJsR1YiOXFdgfp5CCCGEEM3K19/nuTYXlNnP6dvN1yE57pCqkdwuuz73tG6+XSU+21ff57kiFaMQQhwISZQJIQ6Hb4F/AScBcUAXYArQFZMAegd4stb4ocA/gc+AF4CT65zvUuB3QGfgz8DHwN+AbnXG3QjcBHQHngY+AXqH9mUAfwltmwlcRf09qM7HLAP9KfAKcDNgr7U/ARgL/Bv4AHgw9D3WdlzoWh8B7wKTgeNr7e8W2vY+8F7oe+pba/9twK9rvT4bmISp1JsIfAg8D9S3dNopwMuYn9GjQBvgT8BZ9YwVQgghhDggMz//XyzAyAt7RmTV0BHnHFMJ8O6c9Z5InE8IIQ6UJMqEEIdLAJOIsgEdgHGY5Nm9mB5m4dWNJmMSSm5Mgq176L+1+56dB9wPLMBUrP0IDAOWYJJvYVfWGpcIZGGmgvYFlmOSYD9gpqG/jUlQ1fYY8J/QOb8DSoGHgNjQ/g7A98DdwHpgbeh6X4fiBzgGWIpJEi4GVgG9QvECtAvFfW7oXCswCcBra8VxDXBZrddnAmOA/4bOuwYYFDpPpzrjFmGSg0sx02K/xiTzTkcIIYQQ4hCsy9vlWL5uW0wvb4r/lGPbRaSn2MnHtPP36JoY+G7VVtfm7aX2xo8QTW3s2LGtpkyZEtv4yObPsiyefvppT//+/RNTU1NTzj333MQlS5Y4Gz9StCTSo0wIcTgcj6mI+h8m2RQWwCStwv23+gGZobFv1hr3JqZirFetbV2AuzAVZwB/xSSq/sSe0xRPxCTWvqq17VsgHziD3Qm61cDDwKuYxFs/YAImMfZIrWN/C4QbfE7DJP5OBEpC254PnevmUGzDgApMUq72dIRw9dpgIAVID42ru78hcZjeWeHk3lQgL3Tdx0Lb/oz5mQ/E9IsD8/M5pZFzCyGEEC2JG6hi9/u72E+vf7YmVmv49YXHRaSaLOzyQRmVT/5rafy7c9Z7fnvNqWWRPLc4cCtWrHB06NDhqGjyX1VVpd588033vffeW5mWlhacMmVK7MUXX9xmy5YtBU6n5MuEIRVlQoim8jSwAZO8+QmTJBpVZ8wL7Nmk/ipgOzCrzriZmGRb61rbfMA/ar3eETpucJ1j/8eeSTIPpppqOruTZGASbVat44eH9j9R53wW5kbaCVyOmY5ZUmv/RuAbTDUXQAGmV9u9QNta43St/XbM9M0u9exviAZeqvV6F6ZarXvotRM4DdMXrrrWuBcQQgghRG1nYqqy5UHSASguq7J98W2uu20bj3VRf68vkue+4uzulS6nQ38wP8ttWZK/PERDgYj2eysvL1evv/66+7HHHot7/vnnPStXrmywAGf+/PnOH3/80bFmzRrHE088ETtz5kx3eXl5Yw+E67V8+XLHI488Evf888/vc1rutm3bbP/4xz/czz//vGfHjh02rTVvvvmmu6qqSrndbr148eJd1113na9fv37+MWPGVBQUFNi2b98uuRFRQ/7HIIRoKh9gqpwmAhdipv0trDNmU53XaZipiOFkVPjr49D+2j3I8tkzAQSQExpTu0w/r86YVEy1Vnad7aWYJJ039LorJulVVfcbC2mLqep6uE6sGvhV6Dpgknf/xPRk24q5Eb8fiAntn4OZbvogu5OKEzH9x/allD0TfQBl7O6P1glTNby5zpht7P1zE0IIIVq644F5mPftdlGO5Yjw7tz1bl91QA0fckxFpFenjI9z6XNP61pVUFRpn7t8kzT1PzSTMPeYV0TiZKWlpap3797JL730kqegoEAtWLAg5p577olvaPyTTz4Zd8cdd7QePnx4woYNG+xTpkyJ69+/f1JJSUmjybLi4mL11ltvuW688cbWnTp1anv22WcnLVmyxJmenh5s6JglS5Y4TzrppOT33nvP/eGHH7pOOeWUpE8++SRm4sSJcS6XSyu1+7KWZTFt2rS4M888s7pz585HRcWciAyZeimEaCpfYfp+7UugzmsfpgrttAbG167cqu8NOT40pvabZ91rFDVwvMJUrO0KvS7HNL5vSPjJaSamOq2uQK3/jsYkws7GTMUMN/O/LTRmAmZ651mYpOIYzNPtofu4fmOPV3eExiTV2R7P7iSdEEIIIQwf5j3zSsz773OYdgUNPTBr8d6fl+Wx2xVXDelR2RTnv3JIj4pPF+W43/5qneecvt3k93DwnJjk7wzgD5g2HT8d7MkWL17sLC4uVuvXr99VO+m0L+vXr7f//PPPO5OTk63KykrVu3fvpL/+9a+xDz30UN2HvgC8+OKLnpkzZ7qXLFkSk56eHrjwwgurZ8yYUTx48GC/y+Xa5z3w/fff32rQoEHVs2fPLgbo169f4p133tl65MiRe1U9jh49On7t2rX2efPm7dr7TKIlk0SZEKI5WYRpYt8e029sX1IwyaY1tbYNxky13JdtQCFwDvB6re39ME36w+f7FnMjcSqm8X9dhcC60HmmNnJNMAm8j0JfAeCSOvsrMQsH/AeT5HoYU/V7sE+3yjFJx0vYc4rqpQd5PiGEEIdg2/r1/O/LL8lfvTraoey37du3k5KSYisvN59l33rrrcRAIEAwGCQQCIwC/i+6EUaMk90P0ByYB2f3YlbOvgX4PDphNV8/rNvmzMkvdgzs3bmqQ3Jck1TinH58R/8x3RIDi1bmu9ZvKnL06Nqm7sPPlqQt5oHrwQg/+E3ATC+ej3mgfQemBcgBycjICJaUlKiRI0e2vu6663znnHOOPzY2dp/Jq6FDh1YlJydbAB6PRw8fPrxq3rx5Mew9OwKAr7/+OmbRokUx3bp1C1599dW+YcOGVZ922ml+m63xysWVK1c677nnnuLw6wsvvLD6T3/6U6vrrrtuj0TZ/PnznTNnzvSsW7du59HSf01EjiTKhBDNyavA74B3MStJfodJXoUb8v+h1tgqTAJoJCb5dTfmBqKxm3YNPAuMxyTm3gQygJeBXEwiC+ANzMqcb2BuJBZjnjRfhulL5sM0+n8L04/tKcx00M6YqrAtmCmjN2KmOs4BfsH0EBuI6ScGps9ZMuYmPHz8eZgFAQ71TfsxzNPDvwHvYVbgHBuKRxp+CCHEYZS/ejVfPPlktMM4IA6HgxNOOIGff/4ZgE6dOgWLi4tt5eXlqrKycmdRUVFWlEOMlDaYBYNqfzYKYto01G3VIIB35qz3AAwf3DTVZGEjLzy+/KGXFya8+unq2EdHDyhp/Iij1omYRZwOlcIkzC7FzODoi2k9st/S09ODc+fO3fXss896brnlloTy8nI1duzY8okTJ9ab9AJISkra474zOTnZ2rFjR4PlaG+//XZxfn5+6UcffeT64osvYp5++um4mJgYfd5551UPHTq0aujQodUpKSl73ScHAgFKS0tV7esde+yxwbS0tOBJJ520R6J11apVjs6dO1vdunVrcBqnaLkkUSaEaE7KgCHAi5gnXeHHRqWYhv61bcD0QVuNaVBajent9cZ+XOcRzNPiZ9k9bfJ7zHSL8A1fOSZh9U9MkitsLSahB2ZqqQPTf+yeOmPCr1uFrpdYa//XmOmYYKZBPg78vdb+JZiVPw/VK5ifzRh2l9mPxCw2sDMC5xdCCHGA3K1bk37mmXQ87rhoh9KopKQk5s6da91xxx02rTX3339/yZw5c2J/+OEHN/Du9OnTfxvtGCNkMOYhXRLmXqQE8z798T6OabHKKwPqyyV57qQEtzWkT9cmnRJ5yVnpvmdmL2/12eIc933XnlqWkuBpqZU/q4DrDvLYyZg+wGEl7K4oO6AkWVi/fv38/fr18weDQd5880339ddfnzBy5Ehf9+7d6006bdiwoXb/YLKysuzdunXb5++yU6dO1ujRoytHjx5dGQgEmD9/vvPTTz91Pfnkk3GLFi1yvvDCC6V1j3E4HLRt29batGlTTenZ3Llznbt27dorKTdo0CB/QkKCrKgq6iWJMiFEpC3HPK3al+/3MWYTpo9XG0xD/F2YJvj+esZOwazi2A3TpLS4zv6GphkGMJVrfwTSMdMo6za9B1Nhdg6m3L0T5mZia50xb2Cq0tIwSbGtmKmTYc+FYuyEuQHfwp5JqrcxTwg7YqaT/sLeNy3n13k9lfqne9YdB2ZlzNqrY56KmWIiT8ibSGZa2jiUarCpbSRp+M/UnJz5h+NaQojISO7WjetffDHaYYi9xWLuIyYDf6b++w4BfLJog7uiyq+uHNIj4k3863LYbVx3Qc+Kp99eHv/3D1bGZd5w+l7JkRaivlXh99dYzH1qEebe9iZgxcEGsmrVKkcgEODkk08O2O120tPTg0opPB5Pg7MV/vvf/8bMnz/fOWjQIP+qVascb731luell16qe99eY/PmzbbKyso9Pit07dq1JnG2r2tdfPHFVS+88ELsNddcU7Vt2zbbG2+84XE6nfqdd95xXXXVVTWJ3aysLPvcuXNj/u///i+iK7aKo4MkyoQQzVURuxvv70sxsPIgr1HG/jUzLWDfPRw0ZsXNhliYRFx9ybjw8fmhr0hKxiT6PgUqgBMwFXS5mB5sognY7PYxNru9jdPhaNJeKpU+n0uZykdJlAkRRT/8vJ2CXRU1r0vLzcLC363eyi87d89EaiOTe5qzILtbLuxoZGyL9/68DWba5ZAehyXB0COxJKCClbzz9c+eWy89obxtm9iWWlV2sPyY/13fRQSmb27cuNE2cuTIBJvNRkJCgi4sLLRNmzattFOnTg3+XkaMGOG76667WhcVFamCggLbnXfeWXH11Vc3WI04fPjwNkuXLnU2tP+OO+6oqK+iDODRRx8tu+yyy9p06NChrdaaO+64o6JXr16BW2+9NaGwsLB01KhRlQDr1q2zf/bZZ7LAlaiXJMqEEOLo5cJUsyVjerq5MEnFy5An5U1JnXTiib6+ffs22KsjEl555ZXEar/8GoWItqLSKh54eh6xbidKQWWVyZG//MFKbDZFVXWQU49tx296RTlQsS/fhL5EIzZsKXKsyi5wnpTRtjqjy+Fprl9RVszWnz6gwynXqpc/WBk34cYzWmpV2cF6BDPV8qD7yX3xxRc1D6+HDRtWvXPnzoItW7bYfD6f6tq1q+V2u/fZ+7Zz587WzJkzd2ZlZdk7dOhgxcfH73P8kiVLCg821o4dO1pLliwpLCwsVD6fT4UTeDfffPMeid1x48ZVjBs3rqL+s4iWThJlQogj0ZOYPmZi3/Ixy4GnYqZ9FmCmqIoIGtO+fZw9NjajZoPW9oqKCtuOHTua5D02MTExaLfbZTEGIfbTF7lrWfrLpgM6ZnvF/n8OH9KnK32Oa8eK9TsIBHcXVIQTZg67jfE3nkHp93MPKAYhmqO3v/zZozVcMfjwVJOFFa7/moz+N+h356z33HhRr4rObVtJjeb+i3ivPZvNRteuXQ+osk8pRY8ePQ7b7y3U0F/ul8RBkUSZEOJIdGCfeFo2CzMtdF9TQ8UhsMfFnaK0rqlEsLRm7dq1rF27NrYprnf1NdfsaJOQIB8QhNhPy7Zt5q11PzTpNR66pT9Xjv9wr+2uGDvDB/Ugo0sbfvi+SUMQosn5A5b6dHG2O87t1Bf18x7WRJkVqGbQqV0Dny7c4Hxy5tJWf/3dkAb7W4nm5cQTTwykpqbKfYs4okiiTAghhIiAc845pzg+Pj4YCATa2G02ZXc4IvoUs2D7drVg4cLGFsoQQjRg+vlX0zk+4YCOSXTtX747vXMCw4f04IP5WVRV7/48aLfZuPPKk7CCQbQlbZXEke2/Sze6ikqrbJefnVEZ647se9z+OMbbQf8vZ2fgq+/z3N+t3lp5Rq+O1Yc7BnHgHn/8cVlZUhxxJFEmhBBCREBScnIgKTEx4PP5cDgcOByRfYv1V8vnASEORce41nSLT2yy8//26lP5ZMEGwt2p3TEO7r+uD//+zWjWL1jQZNcV4nD599z1boCrhvQ46F5Xh0TZ1Fl9e+icD5fyxOtLW81+/JJCm02eH4nDqhVmMTBxlJNEmRBCCHGAxnm9f7ArdatSCsCNUnz+6aeJdrtda61toe0RuXtPTErSF1xwgfTYEKKZax0Xw2+v6cuf31yKzx8kqZWTc7q7ec93WGeoCdEk8gvK7EvW/OLq3qVNoHePdtFZSUaDy9PKfnqvTv4lq/OdL3+wMm70FSc16cI5Lc2OHTvsH3/8cVy044i2vLy8+lbcdGBWGs8C7sb0/hVHKUmUCSGEEAfIBu2cMTFdjunRo0qDXWuN3W5XCpTGNKyNhM2bN1NaVqaQZrRCHBGuPvcYZn6+hk3bi+m24CX+/Elmzb5z7rmHU4cPx+FyRTFCIQ7Ou3PWuy1Lc/mgjOhUk9XSu1eGWpe303rp/Z/iBvfpWnVst8TDsvpmS1BQUOD46KOPWkU7jmaqNbAZGAEMAR4AXotqRKLJSKJMCCGEOAhut9vfr3//Mm1ZTl9VldvtdmulFFprFalEmd/vZ9v27RE5lxDi0GiteWfs2EbHnR2I45sqJynFe66hEpeYSHJqalOFJ0STsbTmwwUbPE6HXV82qHv0SySV3XbJ2b38Mz9e7nropYWt33h4WKHdLlMwD4VSKltrPS/acTRDtVeLLwQuBS4BXgBeBW4BRgE/H/7Qmr/UCakdtaU9AHlT8rL355i0B9M6WHYrNnRMDnUeFncf072d3+FvBWBz2/JzJ+U2yd+kQ06UTeje/QRtWddFIhghoknZbG8+vmHDqmjHIcSRaoLX2z4mNnbXpNWrj8pmWuPS0zNtWseHXvbz+XzO7779thVgsywLh8Ox+y79EBNlXbp0oWPHjod0DiGag52V5SzdFt2FivNKdu332MriYqxg/Yuzaa1Z/u67+3WeY+2maqz3JZcwYto0AGx2+37HIQ67rpjPRbJCdD0Wr8yP+WVnuf3cvt18ifHuqK1KoWt9Xo6NS3Cednwn/9I1+c6n3vi+1djrT5O+UYfgpZde+gfwj2jHcYT4CPgG+BPwG+BH4AngMSA605JrSRuXdq1SajgAmpKc7JzRzCYIkDE+o21QB58D0EovyJ2c+3RTxmLTthnAhQCMwBGOY1+UXU23Y78EoM+oPjHLpi/zA6SOT+1n07bXLKwMO+b91FZlG5j2YNoGm912AkDQFly98fGNWyMR+yEnyizL6glM8Hg8vkg9QRficNJaU1lZ6bYsawUgiTIhDlJQ6+d9FRWXZHq9P6HUl1rrhSoYXDAlL2//P6U2Y3alHrA5HG1inM4A4ATs2Tk5tT/5RuRNsKysDKfTKYkycVRYX7SDB+d9GO0w9tv0a6/ll3XrGh13+nXX4UlofAXNTscfj91ZX6sb0cx0B/4DPAv8EZC+V7W8N3e9B+CKwVFq4l8PDfQ+IcOWt7XImvnFmrgTuqf4L+rvrWr0QCEiowj4LfAOMB2YCFwO3A4sjWJcKJvqhWaEeQFp3dMW5JL7CkBAB+IUakRoaNSTegfCpm3TgQygWmn1qVZ6uwqoX2x22xCNfgNAaXUTptLvkEVs6uU111xTEhMTIz1UxBHHX12tZrzyijvacQhxpFNKzVEQp6E/WvdRAHZ7cFx6+goF3yitFwS0XvBkbu4v0Y71IKmTTjzR17dv3/Jqvz9eW5bH5XJpIKLTLV97TdpdiKPPBanHMqCzN6oxpHgi15/6rNtuI8Ub3e9HRFw5cAdwPeZD74tA1KqnmotN20rtX32f5+7aLj54Vu/OUa0YV3U+aSplV8MGnxSc+dFS9ae/L044pltiYUaXNtKvTBxO3wCnAJnAeGAR8Dzwe5rJ6pgKNSnjnow3s57J2u9Ecvsx7eO2Tdt2UA8MMu7JcB3ItcK6/K6LZ/NfNlfGumN/XeordQKEq8kGTxrs2Ojb2DM09OvsqdlXhI9Lz0w//WDibMxBJcp+7/WmBqANgIJUDezcudPhdDolUSaOKCkpKfJmKkSETMnJeRZ4dgTY0zMyjlPB4ABgoNJ6MPBbDb+1K0Wm17sVWAB8qe32hVOzstbQDJvVj2nfPs4eG5tRs0Fre0VFhW3Hjh2OQCBgQ2ucMTGhXfqQG/gnJiZis9kO6RxCNFe9UjowvMdJh/26Ozdu5N8TJgAwaz/GF24y00SvnDqV4887r8Fx+1NNdqTasWMHxcXF5OfnU1xczJw5c1R+fj5aaxYtWuRct26dLT8/n9LS0i6YhtZHg5MBO+AJfU0G7gNuwnzwbbGenLm0VTCouXFYr3KbLbqzh3Q9l7fHeOxXnntC9ZtfrHD9ZtpXbV75468KOyTHtfgEpzisfMAk4F3gZeBeTB+zOzCVqtGkgdRAq8CdwF/3NbB7ZvdztNZ/0EoPAGK8md5y4Eul1R+zp2avbOxC3nHeG1BMDhLs5M30bkLx+/rGpWemX6nRT4ReZirUORp9ExofkFjhq/iLHfs5AIMnDT52Y9XGyzdWbnwCRXgWxyBvpncDgELN1+iLwudWWj3pzfQ+BIDFRTlP5DReIt6Ag0qUBeFJZVZ7qPlk89FHHyUdbBBCRMuoUaO2RTsGIY42syFIVtZqYDWmHJ3Mbt3Std0+UMEAlBqI1lcBI1QwSKbXuw1YitYLtM22cFdi4nfTly2Lejm4Mzb2VG2WAQdMM+O1a9eydu3a2FrDIvap4dprr6V169aROp0QAqgqLWXD4sUHfFxMbCyxiYlNEFHzt2DBAmbN2p1WfOutt2r+zn388cfxtYYOD30djVqHvhZgplc9E91wouO/Sze65izf5D4uNdE/4pxjoj7tsm5FWVirNknOC87MqP7Pt1kxo6b8N/HVh361K5q91ESLtRLoD9wGTAO+AGYDdwMFUYlI8y6Kq5RWE44de+w/qqi/0Cs9M/1KC+ttFOEntsVAAnCZVvr81AmpZ218fOPyhi7jHe+9FM0r7L4vbovmVcwU1T1YWPEKlR56+ZRGdwVA1QTXHkgHKM0vVSpZxWula5dwx4b3a/TPQLta+9qGvlBKxTQU7/446KmXySkp1QMHDCizLCsmGAy2crlc8sdIHDFyc3PVDz/8IE31hDhMpuTlZQPZhJbRfjAtrYNNqdOUUgPQeiBwIUpdrLQmqbCwLNPr/RalFgILqqqrF/5l8+ao3aCfe+65Ra1atbICgUAbu82m7A6HRmulISLTLQsKCliwYEEEIhXi6Pbqbbexbf36AzomlBwREgAAIABJREFU6Dc5d7vDwbA//nG/j+t84okHdJ2jSe/evUlNTWXBggV4PB4uv/xyPWvWLKW15uyzz67YsGGDc/Pmzc6SkpK5K1eufCfa8UZID0xvofCDkErMVMzwFMwWp6S8Sk1+5bt4u10x6fYBJdGuJmtMareujsHVgeq5y3NjRk/+ss0//nDhrvhYme0kDjsL85D4Y0zPwxHAeZipmdMPdzAa/XcUpyitulfZqu6nvv5dk7Dh42+ADahSqLOyp2Qv9Y7zXo3iLSDWFrQ9gfk+GvJHwkkyxWU5k3M+9GZ6RwBvNxJiO+AepdU8oFt9Azwuz7/Ky8s/U3a11ZxefRYg8H8A8e54f4Wv4lpMJR9Kq7sCKjALIC87r6SRa+/TQSfKXE6n1b59e38wGLRXV1fj8XgOJQ4hDqudO3ZEOwQhWrRQn7KPQl9M6tWrla+y8kxgIFoPCP33PACX0xnI9HpXACZx5vd/9ZfNmwsPV6yJSUnBpMTEgM/nw+Fw4HCYt85ITLcE8PujXjwnxBGh+Jdf2LXp4FbQtDkc9Lv++ghHdHTyer306NGDXbt2kZCQwN13361XrlyptNbccMMNvjlz5th++OEHJ/DTypUrn4t2vBEyGNObzIPpKzQLGAMc0getI5XPH1R3T/uqTUFRpf2Gi3qV9/ImN4tWJfVNvayhoEePNIfP7/d/u3KL86aHP0t87sFzi2QapoiSfEzF7QhMwuwl4ArgTiD3cAVhV3a/tvRDWul/KdQDCvWprtPtJKMyo2dQBTuHXn6RPSV7KUDO1Jy3vZneR4EeKM7uNalXzOpJ9axsPwI7GtNbQfNjzpScDwFypuTM9mZ6/wf03OuY3V7OmZLzbOjf9U7vXD1pdXXGPRm7gnHB0CW0P2/K7oXC0jPTy8Pfk2WzKvImR2YRsf1OlI1PT8+yKWXmhSrVtmDHDtdbs2alaG3+ZCmlpLGKaPa6duum+/fvL0+XhDg0EX+sPGn16jLgy9AXkwYPdvg2beqtg8GBSqkBCs7V0Ae41+V0Wpnp6WvReoGGhQ6Y91hOzsbGrjEhNbVj0G6/Y2p29sRIxy+EODBV5eWU7dx5QMeEq8MGjR5N+hlnHNCxSvr/ica1Ab4CRmMqoFskS2vGPj2v9Y8/F8QMOKlT1f3X9mkWDcmh4amXtUZw4vE9bDabrXrRik0xIyd+mvTsg+cWHZea1CwSfaJFmg38F5iKqVpdDTyMmZoZPBwBZHuyZ3l93rFAbzQT6t7FW1i1+wxsr3P4NkzFraOoqige2OuNOyMjIymog2aao9rr+O3sI1GmlGpwOme07XeiTGvdvX3HjtXJSUl+rbULpZTdZrPp0AemSK32JURTyc7OprxcVvsW4iClYpa9Di8p/WFTXmzS3LkBYFno6291FwhA60HAKAWjglCzQIDWeqGy2xdM2bBhOXUWCAjabGcprX8/LjW149SNG0fX3S+EOHxWffYZ74wde1DHdjr+eI4dPDiyAYmWbitwPvB1tAOJpgpfQI17bn7C3OWbXL28Kf6/3Dek2G5vPp/x9llRVkOpXsdmOBLi3NWfLVofc9MjXyQ9fseA4nP6djvgVfiEiJAiTAL+X5jpl1OAazGJs++b/OqTsNR4NUFr/YlW+rK6u22WbWvQXlOtlV7rOJvyKW+oWqtiS9aWvfqNAWRNztrhzfRWYKau1+4lpuq83ptp4H9ILCxLhbN/OnIP8w9o6mV6errv+OOPr6yurrYD9piYGA2gtY5InxYhmtLOA3xyLYSgG3A1cA3QN7RtJ/AW8M/DGUh9CwSM7dq1k83hGACcF14gQCk1Assi0+vdDiwJLxDg8XiW+Corh6K1Xdnt12WmpaW409KuDiXkhBBCtGzrQl8t1paCMvs9T33dZv2mXY5e3hT/C+POLfK4HEfmAyUFXbp2cYy4wF39769Wx9z31zltLhuUUZl5/RmlcZ4j9HsS0XYcEAf8yJ6VYN0BP5C3H+eYD5wCjAPGA4uB54EJmJ6ITSZ7cvan3kzvN8BZdfdlPZm1wZvpXQWcoFCDveO9t2LjS1WlbtFoMyVT8TGzG6yA02gWozgX6OEd5/2Dtul3lFa/poG+YxFWszifQl2ZPi69wFJWee6U3LmHctIGE2UTunc/QVvWdeHXGsjNzXWVlpTYtdYOQNntdhMPgCTKRDOVkZFBYgtduUqIg5AMXAncgFm5R2FWvnkdUz7+BbB3f4IoeGLTpnxMTLPBLBDgUGqghrO0UgOV1kPDCwT4KipKAScAWrdCqfN9eXkfjOrT5/LmsMKmEC1V1969GT179gEdY7PbGx8khNgvWsOsL9d6npn9Q6vS8mrb0H5e38OjB5S4nfZml1BqfOrlntokpjhuuPQ0/38WrrG9Py/Ls3TNLzGPjB5QfFrPDvK+Lw7UP4F+wC3AjFrb/wxsAn6zn+epBCYB72Ea0N8LXAzcgZmi2WRsyjbe0lb9q0dZ3IqN/wAJaP5O0GS/QvLsQfsD+zq3Umq8Rs8DPCgeUVo9gvle82jiZFlloHJprCN2K9ARuFgrfbFCFWOm0x+0BhNllmX1BCZ4PB5fuFpsV1GRq7i4OEbvTotJdkw0W1prysrKSE5OlkSZEPuWhHmTHgH8CvPeUAl8gklCvQNURC26/RRaIOCd0Bdj2rePs8fFnaIsawBK/Qo4u9bwVmg9OKmw8KvfdelyYTRX1RSiRVMKu+Og15YSQhyCn9bvcD4xc0n8iqwCp8vp0L+7rk/pLRef0Gzf7/dv6uWe7DGx9mHn9NEbsvOqv1qaHXPrY18kDT61m+/eq08pz+jSRqrKxYHYgklyvQmHPGXwJ0zi7TbgKeA/mHvuu4Daq851BgbQ+OqRNbSlVyulzBMo2+6eYRsmb1iYlpk2WaEyQpu+C+/LeSJniff33t46qH+rUAMxD87zlVJfV/uq/5bzl5x9LqKVPSV7adr4tCFKqz9iepqt10o/orQ6D+gNQC+0ebRNbuh7RSm1VyWeRi9UKB/Asl3LLICEqgSrMK5wNoBWeo/pqtumbStPH5d+JjZu1VqnAy4i8Lml0TuTa665piQmJkZXVla6nTEx2mE3Txe01JCJZq66uppXXnkl2mEI0Vy1AS7FJMcuxFRb+YDPMW9e72FW/jpiTdu2rRxYACwYn56+Eq1P1ns+XbKhVL+YmJjHgd9FJ0ohhBDi8FqyZqtz+vs/tfpu9S8xAP1O7FT1x5vPLO3aPv6wNBc/3CytlNeb6rihY0r1/KU/2+Ysy3PP+2GT+1dnpFXefdXJ5d06tD4qv28Rca9heovdCfylgTHtMVVmZ2GmaM4F7gfqW4nRwrQT+Q/wIuaefBAwNnSt8DVPAj4A9qvPXu7U3FmYlXv33jcld0JDx+U8lrMxFOtByZ2c+x3mwXtt3+01zkyJnNvgeabkTq27bdn0ZX5MO5h6ZU/NzgMiuljXHomyzLS0NFvoQ4QFaQA7d+50OJ1OXV1VhcPhwOFwoLVGKyXlZKLZSk5OjnYIQjRHHuA84HpMksyFeROfg5la+T5QErXompBlWecrpRI0lCpQWqmfbUq9amn9/tTs7P3pKyGEEEIcsfILyuwffrPB/cmiHHfu1mIHQC9vin/U5SeWt5RG9zHuOMcFZ5+iS4uKquYtW+/4dHGO5/Nvcz2nHNuu+pKB3SuH9vNWxbqlh5loUDWmomwa8HegtJ4xbwA2zP22EzNl8x/A8H2cNxczo2ME8BzwKqY/8MeYFd8dmEqzhpJzognskSjTSj1pwVW1t3300UdJtV4qZLqlOAKMGjUq2iEI0Vy4MSt5jQCuAFphkmPfYirH3gAKohbdYaKUytDwrrKsNwI+339C1WZCCCHEUam0vEqtyNrhXLQyP+bbVVtdWZt3ObQGh93GgJM6VV0/9PiKASd1bhY9Rw8ny0LFtW7jvHjIabq4eGfVstWb7cvWbotZtnZbzNTXl+jBp3atGnBSp6p+J3aqbpcYa0U7XtHszAQeBB7AJM1qOxE4B9PgPzu07U7MCu6pwMZGzj0b+BKzKubtwFB2515+D7zEEdAK5Wix19TLlJSU6gEDB5ZZlhUTDARaudxuC0BrbZMMmWjuNm7cyPLly6MdhhDR5gIuwCTHLgNaY8q7F2PehN8CfoladFEwJSfnkmjHIMSRrPiXX3jtIB5ClfurOatkF3mxX/CMJ65me8Wu+mahCCH2V2VVQO0q8dkKS6pUYWmlbfuuSnveLyX2jb+U2Ndv2uXcvL3UrkO1UXa74qTubf3nnZ7qu3hgui8lwdPiE0AWqPiEZOfg/skMsar8m7Zut5as3Oz8bHGO+7PFOW6AjslxweNSk/w9uiYGuraPD3Zp1yqYlOCxUhLcVus4l1SetUwWJmk1E1P9VVsPoIjdSTKAFUAgtK+xRBmYKZqjQ+MH19oei0nOPXIwQYsDt1eiLCYmxmrfrp0/GAzaq6ur8Xg8gGmMrqQrmWjmdu7cGe0QhIgWO6Yp6AjgOqBtaPsaTH+D14H86IQmhDjSBaqryV+16qCOTQCqKJA/QKLF+v2LC1p/+M0Gz6GcI8Zp10HLUsFg4/kZl9OhT0hv6++ZluTvc1z76oG9O1VLYqdh2uayd+nc1d61c1cqK8r9hbuKdc6WQltufpF9zvJN7jnLNx3yNcaOPL30+qE9pRro6PAhphl/3X5flUBMnW1OzD36gfzuzwJOZc+ZfB5M/7BnMMk40cQc49PTs2xKhdfZbldQUOB6a9asFK3NuiJKKVvNaEmUieZIa7qlptKvX79oRyLE4WYD+mOSY9dgGoiCSY49j3nalRWd0IQQR6sR06Zhs9sbHwhkFe3g7yu/Zai3J2d36b7X/lhZlVo0D30wH06fA/yRPnnHlFbBY7slRvS8bpdDJ8a7dUKcy0ps7bKSWnustI7xwdSOCYHU9q2DdvvR97kt4K+mXUJsxH8/e0hoQ7eObTj5+NSghmBJaSU7CktthSWValdZpaqs9FPhq1aVVQcWRgCbrfFR4ggyAdOE/39AOJP6EyYxNpjdzeqHYhJo/9vP87qAVzD3+OEKEIUpcPIA44DxhxK42D8OrXX39h07VicnJfm11jEopew2m02HMphSRSaau+zsbMrLpd2QaFH6ADdgekp2Cm1bg1kx501gXZTiEkK0AL0vvRS7o9GF0wHwbd3IFnYS32cAJ59wRhNHJsRBiwemYqY23QZ8EcmT/+aqk8t/c9XJcrN6iEYMaFeUnJzc4qeNimZhPvA1JhG2MLRtC/A0ZsXJKZjqsvGhf+9vv4EgcAambUocJjkW/rcbk0ATh4EDID093Xf88cdXVldX2wF7TEyMBtBaK0mUieZOpluKFqIXpnJsJKZJKJhVcp7GrI4jzfmEEABUlZfzt6FDI3pOKxCI6PmEaIbKgC6YPp4/AzdhHkIJIcQLQE6dbb8D5rHnPfhYTPP+szFJr9uAdw/gOgFgR+grotIz02dodBeNLsidkvvrSJ//SOLN9H6ESTyuyZmS89v6xuzf40AhhBDREE6OXQccE9qWh0mOzQYWRCkuIUQzpi2LXZs3RzsMcQR57LHHyM7OpqCggJycHBISElizZg0xMTHccccd6ptvvgGgsLAwbtu2bY6dO3dSWlo6BPPh8WjQCVOxAaatXl/MAjgzMRUhJVGKSwjRPLxez7Z1mErU2izM7I43mzyiA6TR/YBjFeqgmu6lj0v/i1b6LICkwqR+y6Yva7Jp0N7x3lvR3Amgbfq23Mdzf4zwJYZg/ubHNTRAEmVCCNG8pGH6jd0I9AxtC5dyz8aUd0tDXiHEfrnmr3/F07p1RM+5v/3JxJFj9uzZlJWVAbBmzZ5FVD/++GPN9JI1a9a4au06MfR1NFKYVeZuAZZiegaJKPF6vcHf/OY3FR6PR+5/hIgSrXR3TPsXfJ18TTvt0KIjylzLFrS1atJrNUASZUIIEX3dgCsw1WMDQtsKMU+vZgOfYsq3hRDigHQ/80zi27WLdhiimZs7dy7z5s1jxYoVDBw4kPT0dF577TXi4+OZOHGi9eCDD9q01owePXrXwoULY1etWuWqrq5+edasWb+PduwRMgD4JxBeXaIE+Aq4C/glWkEJo2fPnsFnnnmmNNpxCHGk6Dq2ayeHzXEzplp2hb3c/mpwHx8l0jLTBivUOUB7FDuAz3Mm53wD0GdUH2dhUuHNQGp4fKWv8jZvpjeAojBncs474e2pE1KH2LTtHDTtUOywYft0w+QNC+te79ixx8ZXq+orFOoUy2a5lFZbtdLzcyfnzvNmekdgqnoB0Epf6s30Hg+Q6k7959xJcwMA3Sd0P0Fb+hKgm9a6wlLWgo0bNn7I7D2/0YwHM7oHHcGbgEQsvs3Jztmvaj9JlAkhRHR0xjTjH4FZuVJhGn2Gk2Of0wQrbwkhhBB1tWnThvj4eNxuN61btyYxMRG3201sbCxt27bF7XajtSY5OVnHx8drt9uN2+2uBAqiHXuEFGE+F5ViWhzcgPT+FEIcgdLGpx2ntPoGSAlvC8YFb8FUye6hy++6eJwu52xgWM1GU7c5IT0z/ZVsd/atZRVlbuCl2sdp9HOhf6wE3uk0qVOs2+d+R1t6aO3zWFgT0jLT/pnrzr2dSVgA3nHeQdWqejbQTqNR2hSnKa0KgHYKNUmjj691uQfD//iZn2cCAW+m9wnLssYQWoASBTZs96d1T1vEJIblTsotAug+rnv/oAr+Bx2aYqm429vdey37sSiCJMqEEOLwSca8EV0PnIP5I10BvINJkH0BVEctOiGEaMSuqkp2VJTt9/jNpUVNGI0QEeUD7sP0FpIpfs3UqlWrHO+8845r0qRJsoqoEPVQWr1IKEmm0W/asC1Ac6tWumvdsU63cyq6Jkn2pNLqQ42+DMUYjb7JW+ldHqwMTldxaqpGXwlkYE48TSkVtJT1C4CryvWExiTJFGoqmo81+goU9yvULem+9GXZZD9/zAPHpPjxfwC0AbRGz0Ax32bZOmqlLzbfAP9AczlwViiufynUZoCO+R39rnGuGwglz5RW72ub/rPS6kSN/ptC9cfH34AbmYTN8lkzMH3ItFb6RaXVGuBezGqi+ySJMiGEaFqJwCWYyrELASdQiZlOORuzEo7c7AkhjggfZq1i2vdzoh2GEJH2E6ZHaEWU4xCNWLdunX369OmxkigTYm/p49ITNPrs0MsV4dUtM8ZmfBhUwVygpslon1F9nIW68Kbw2CDBySiwqqw1Tpfz10AnFLdmPZP1DJAZmv6YARDrif396kmrq2ud50YApdXygApMrXWe64COGn0b8Lzf6b8WkyQDeC53Su49tcKfCpA9OfvP3nHeWJRJlCmtXsyeml2zgJk303t76J8+v/bfrbSqBFbZsA1TqIuAa9uPaX+Xq8qVyu7F0D7NnZx7F0B6Zvp3Gr2ksZ+lJMqEECLyYjGVYzcAFwAxQBWmYmw28G/M9A4hhDginZ96DB3jEvZ7/IkpHZswGiEOWWG0AxBCiEMVVMF0W3hWoeb78PasJ7I2ezO9W4Eu4W1FbYs6ECQ+9LK3HXshgN21x4I9x9CIwraFnQjSCkArfWoj5+kR3qCU+qzOqfa3kjd8DrfD5thSz/6YeHt8WtAKpocmZqJQNT+L7A3Zy73dvQEayYVJokwIISLDA5yHqRwbjinzDQLfAq8BbwHFUYtOCCEi6MpjejOgkzfaYQghouibb75xvv766+5LL720aurUqXHV1dXq8ssv92VmZlYotfeieDt27LDddttt8U8//XRZt27dahpuv/LKK+6srCz7o48+Wr569WrHtGnTYlevXu2wLIvTTz/d//DDD5enpKRY9cUwceLEuBNOOCEwYsSIqvC2++67r9VFF11UfcEFF1QDbN261fboo4/GLVu2zOlyufTVV1/tu+uuuyrri1GII5lSqrIm3aRw1dntrv3C0lbtdi+5wH/qOV+jyStLW9W1knM5KP6713lQVmh/dTh5pbWOa+zcDQjHXQ78q96YsEqVXVVoqyb8mu89rVeaEx+NLt8tiTIhhDh4buB8THLsCqAVYAGLMZVjbwLboxadEEIIIUQTycvLs7/xxhuen376yfnwww+XlZSUqHvvvbd1ZWWlevjhh/eaGpmSkmLl5eXZZ8yY4Z44cWLN/ilTpsTdfvvtlQC5ubm2Xr16BW666abKQCCgnnjiidjhw4cnzJ8/f1d9MSxcuDDGZrOBqdwH4Ouvv3Ydd9xxQYCSkhI1cODApFNPPdX/+OOPlxYXF9seeOCB+JKSEjV+/HiZaiuOKo4yR04wLlgFuIDz08elJ2RPzS5Oy0wbTK3m/gA5j+ds82Z6N2JWs4xRWo3Nnppd81C/+5ju7XDuUQFWrrVJPJVVlyUDWwE2Pr5xqzfTuwnoqpSKsbltD2ZNyioJH+ed4G2vgqo7ADa+q1U3dlfapLSPcifl+sCsTpn1ZNYGAG3TZeEm/5bN2iNuYAnQFfAEVfBveZPz1oR3dJrUKTamMuasnKk5eRljM6ygLagBZWFdnDYpbVLupFyfqlRXoWg0Sy6JMiGEODB2oB+mIf+1QGv2TI69hSwlL4Q4zD557DF+fP99AMI3skII0dTKy8vVjBkzinv27BkEqKqqKr3rrrtaT5w4sdxu37to4/rrr/c9//zznnCibPHixc7s7Gz79ddf7wMYNmxY9bBhw2oqXU455RR/SkpKu/z8fFunTp3qrSrbl1deecVtt9v1rFmzisPx2O12Ro0a1VoSZeJok/VMVlX6+PRXtNajgfZa6VXpmekrNXoQ5vPKHqs9KqX+pLX+J9DJslnLvOO8HyiUTdt0T0tbgzX6n8BCAG3p1eH0ki1oW+PN9BYC/8qZkvMQij+h+btGdw5UBZZ7M73vK5RDK30cFkMsZU0HFqW6Ut/Pq8xbrpU+FRisfGqdN9O7FGgfJNgD6ACggmp1OFKl1buhRNyKnCk5lylLTdE2fTHgsmP/Jm182ltKqzKtdXdVqc5BsR74IjTd9BPgYoXqhY+V3vHeDWiGYKZ57jNZJokyIYRonA3oj6kcuxZoF9q+BpMcew3Ijk5oQggBVWVllO3cGe0whBAtTGJiohVOkgGcddZZ/pKSEpWbm2vfsmWL7f3336+Z/vXwww+XX3/99b7MzMxWCxYscA4cONA/Y8YM969+9avqdu3aWQD5+fm2MWPGtJo7d25MUVGRLSYmRmutyc3NtR9MomzlypXO8vJydeWVV9Y0VSwvL1fbtm2zFRcXq4SEBHmyII4qHpdnTIWvIhX4FdBFozsDjwNXAcfWHps9OXtG2ri0GKXUI0qr7iju1+hwt7DNSqvl4bFBFXzOjn0AcD6KNkAbrXQKQM7knH+kj0uP0Ur/SWnVHXig1nk2KdQPAHMnzQ1kjM/4VZDgs2iuArqFvsBM/wQg54mcL72Z3ueBWzAzeFKBHQDZT2R/n56ZfrFGP4fmGIW6EyA0lXoXmm/C57EFbLdaTusjNKcDGWhSgfuBKZg2OQ2SRJkQQtSvdnLsakJPODDJsRcwc+LXRyc0IYSo30nDhnHmyJE1r2MTEw/4HOfOfp7S6qp69wWsA/6cKoQ4itWtYA2/bqj/V0pKijV06NCqGTNmePr27Rt455133C+//HLNNK3bb7+9tdvt1suXLy/s0KGDVVVVpWJjY9sFAoF6z6eU2iuG6urdrZeCwSBerzc4evToytpj7r///gq3e4+WTUIcFVZPWl0GDE2bkHayDuqOduyrs6dm52VMynii2ldtd2jHHm/kuVNzX+ozqs8/d7XZ1RtFe23TJfaAPT88DTIsb0reLmDo4EmDHdm+7HiAgCtQc7OQPTX7hT6j+vy99nmCOrglb0reHsUEWZOzCoBrjnngmJQqV9WJ9qA9NqiC+V63d2UOOeFhOmdKzt3Ab7pldmsDEOOOqUnIZ0/J/hI4rtv4bj2VpVIdylHlV/5f8rLy1jGbmnEbpm3YziT6pfvS+wR1MMURdCzbMG3D9rRJaTMtn6U8lqf+PyxIokwIIerqhUmO3QCEO1WvAV4CZgFroxSXEEJgBYNsX793jr6iqAiA1h064D3jjEO6RkXAj0LRI7FuW5DdWsfIB0whBBQVFdlWrVrlOOGEEwIAc+fOdSYkJFipqanB9PT04KBBg/x1j7npppt8N910U+uBAwdWOxwOLr300poP2z/88IPj6aefLu3QoYMFMGfOHKe1jwR927ZtrS1bttTM8SwoKLDl5ubWvD755JP9n3/+ecygQYP8cXFxUj0mWozcx3N/BH4Mv67dN6yuZdOX+WH3Kpn7MnfS3ABQb8/AAznPz0/9vAOYE369kY31DdOhBF39+0x/sjUN7DcmYWWTvbT2ptxJuUWNxSeJMiGE2J0c+z8gI7RtI/A0ZmrlgijFJYQQeygrKOBvF13U5NdJb5PMzItGNj5QCNGitWrVSt98882tH3roobLi4mLb2LFj4x944IGK+vqThV188cVVbrebMWPGxP/617/2OZ3Omn2nnnpq4KmnnopLTEzU+fn5tmnTpsXu61xDhw6tuuuuu1r37dvXn5ycbD3zzDOxHo+nJiF2yy23+F588cXYCy64oM19991XkZycbK1Zs8axYsUKx8svv1waoR+DEOIoI4kyIURLFU6OXcvu+fqb2J0cWwjIk0chomh89+4DXF27fjdp7twGS+NbMrtj79s4Wz3bhBCiqSQlJVmPPfZY2bRp02J9Pp8aP358+X333bfPJvkOh4NHHnmkbNGiRc4777xzjymRf//730vGjx/faty4ca3S0tKCr776aslzzz3nCVeYpaWlBa+55pqaY0aOHOnbvn277fXXX3fHxsbqhx56qHzevHnOnj17BsAk8hYsWFD41FNPxU6fPt1TWlpqS0tLC44YMcLXFD8PIcTRQe6mhBAtSSpwOSZBNiC0bQcwHXgdSY4J0axorSf58vKOG+v13vJETs5/ox3HjkK9AAAM1UlEQVRPc/PHH37AFbfPXrRCCNHkLrjgguoLLriguvGRu91+++2Vt99+e2Xd7R06dLBmzJixxxSx2pVfffr0CfTp06cs/FopxQMPPFDxwAMP1CTnzj333D1iSUpK0o899lj5gcQnhGjZJFEmhDjadQWGY5Jj/TFLARdiEmOzgc8AqVYRonmKQesuNng30+tdq4LBmybn5e27F8UR6omlX/Pd1nr7c+zBUVRCeujfIz+dieWKiXgsFf4D+rwrhBBCCHFUkUSZEOJolIJJjt3A7uRYETATkxz7HNiruawQR7LxnTsna7e7d7TjiCSldZtQiWe8hr7Y7Yszvd63tVJjpmZnF0c3usjaUlbMz7sKSHLH7nNcjK+iJlG2s7KCoBX5P2VtXB7auKRZvxBi304//XT/ww8/XNb4SCGEOLJIokyIFmZceno3LOuUqbm5H0Q7lghLAi7GVI79CvP3rRL4BHgN+ACQMglx1LJcrjOVZX0c7TgiqfY8aGUS3q1R6ial9bDfd+ly2mObN2+JVmxN5fMrR+NxOBvcX/LLL0ye8jIAn105SqZeCiGipkePHsEePXoEox2HEEJEmiTKhGhh7Fq3seC9TK93mbbbb56albU62jEdgjbApZjk2IWAE/BhKsZmA+8C0pNCtAgOrVcFlbo72nFE2O/QOqPW62K0XmtZ1s1TjsIkWVhhXh6v3X57vfusoHwmFUIIIYRoSpIoE6IlUqoQ6EswuHB8evon2Gz3Tc7KKoh2WPvJA5wHXA9cBsQAVcAXmOTYvwFZ7lu0OI/l5GwEno92HJGUmZ5+DZCBUmVAsYbRU7OzP4l2XE0tUFXFtvXrox2GEEIIIUSLJIkyIRoxzuu9Gmgf7TgiRWvdCfCglFKQoLW+SgWDF43zeh/zpKb+ddLcuc2xsb0bOB9TOTYciAOCwLeY5Ni/MKtXCiGOIgqCWutyrdRET7duTzfTv09NRinFr599tsH9Trf0ERNCCCGEiDRJlAnRCAVjgNOiHUekaKXqborREKO0nurbuNEDPBKFsOrjAi7AJMcuB+IBC1iMSY7NArZFLTohRJOztH7R5vePmLJly06ys6MdToO++//27jdEruoO4/jzu3dmNpNoNpIabBEkKbaKgqJiW30jFGptsbVUYxG0Vij1TWuhELB9USmC2RaktCRN1lJjlWhiFJGstZS2EhtjG4Nt0iZr2d1JpJvdNZGN+28md+69py/ubIzEZDP7Z87M7PcDww5n7sx92BczwzPnnDt0RK8Pzi1f3+jHd/1X33bbnF4XAAAA9aEoA2aSpt8zs+W+Y8wXC4JPp849LqmzNjQm5447s+90lUq7fGaTFEr6grJy7B5lV6+UpH2Snpa0XdKQn2gAGq2rVNruO8P5+Oexo3rqP3t9xwDmjXPu1N/x8fFgenzjxo0rqtXqGb+4AQDQTijK6rBuxw79o1Q66+NrLr5Yv7v//sYFQkOsP3Lkbd8Z5tOPV68edVIamEXOuUkz+0lHqbT5kWy2lg+BpJuUlWPfkrSqNn5Q0gZlBVm/n2gAcP4evfkrunLl3FbqL8nl2WQRTWW6NJOkSqUSnONQAADaAkVZHW5cvVqfuOCCM8Z7h4f12jvvaHmx6CEVUD+TLlSaPuGC4OH1AwMfeIpxvaT7lBVkn6yNHZT0G0lbJf3XUy4AmJVLL+zUZy662HcMYE6cczIzJ8kkpUuXLh1TdpVpSVqn7LNaZnb2X48BAGhhFGV1uPP6688YGxkb0zc2blRHLqdHbr/dQyqgPi4MhyxNr3qsVPJRRF2lrBi7V9Ka2lhJUpekpyQd8pAJAIBF7fRZY0mSqFqtSpLK5bJt3bq147TH93R3d/+t8QkBAGgcirI5mDx5Ug9s2aKRsTE9vnatrrvsMt+RgBk91td3TNKxBp5yuhy7R9LltbF3Jf1K2ab8fOEGAKBJpGmqOI5NkqrVqu3evZslEwCARYWibJaSNNUPt23ToaEhrbv1Vn392mt9RwKayWpJayV9W9KVtbH/6cNybLck9/FPBQAADeaccyZJQRDEYRimkgqSysoupCNJMjOuNg0AaHsUZbP0s5079edDh3TXDTfowVtu8R0HaAaXSbpD2eyxm2tj7yvbjP/3kv4ifxcMALCI1fZc8h0DaFruw7WVNjo6mgwPD8fKirIT3d3d9/tLBgBA41GUzcJvX39dT+/Zo8+vWaNH77jDdxzAp0slfVNZOXaTso1/R5WVY89L+oOk2Fs6AKg5fQ+mery97+/a/sB3z3lM6py+5FK9/Mtn9Uow94sCpkly6v5scwP1KJfLQRxnH9cHDhzoOHDgQIfnSAAAeENRVqe/9vaq69VXdfmqVdp0773Kh6HvSECjrZT0VWUb8n9RWTn2gaRnlJVjf5QUeUsHAPNoslKRTUye85iwdotU4c0PLem0q1xK2Q9c1dr9sp9EAAD4Q1FWh38PDur7zz6rzmJRT9x3n5YvWeI7EtAoF0m6XdnMsS8re+8oS+pRVo7tkDTlLR2wiMRJYtVqlXWE58HM5rzc+3PX3ahPbdt2XseG87y808xYMoqFFkhSR0fHaBAESyUtkbSxu7v7Ib+xWsOuXbuK+/fvL/jOgfYTRRFv/oBHFGV1ePCZZzQVRfraNddod3+/dvf3f+TxFUuX6rarr/aUDlgQgaTrJL2n7P2iIullSdsk7RTlGLAgnGpL7pyTc07lcjmIomyu0os7dqw8dSAlyseK49jGx8fn7Z+z6oor5uul6jY+Pu7t3Fg8NmzYsIICvn779u3jV3MAaEMUZXU4MZV1As/t3avn9u494/HPXnIJRRnazY8k/VrSa8r2HXtJ0pjPQMBiYJKrzSQyk5JCoTARhGFO0rLUuYcCaUSS0jTd7zVoE2NvL+D8VSqVuW+ut7i8pGzFNbDQhnwHABYjirI6vPHww+f84h3Owwa+QJPpVrb3GDPHgAaqfdI4y+67MAwrYRAUJC1zuVzP+r6+/nM9f7HK5/OVMAy5gAhw/lbU/q6TdFCSzKzkL05r6O7uvtN3BgDAwqEoqwN7kmGRoiQDGsCddkuTRFEUWZKmNjExEfb09HQy4+Pc3nrrrSWDg4N8rwFmIQiCNzZt2rTbdw4AAJoBXygBAGgCVptBJsmcc4mTEudcLo7jycHBwT9NH5dPkglfGZvZwMBAfmBgIO87BwAAAFobRRkAAM3AzEmSSbIgqHQUCpMWBMslvbu+VGKZz1mY2cvOuX/5zgG0smq1OuI7AwAAzYKiDACAJhDHscVJotQ5O378eP7QwYPF0RMn2Cx6Bps3b/6p7wwAAABoHxRlAADM0rxcWbH2GnEcWzWK5JzTkcOHO44cPtwhSWbG5vQAAABAg1CUAQBQP2e1pZLzwCRZoVA4XiwWE5MuSs0e+fnAwPrpc83TeQAAAADMgKIMi1o1iiyKIvOdA0BrWMjGqjajzJxk9uEFMAEAAAA0EEUZ2p5zLrtJStLUyuWyRXG2kunJLVtW+U0HoJXE1arF1eq8lutRFEmSXnzhhZWnBo3+HgAAAPCBogxtz7IlUtnaJikuFApTgVlBUtFJd08f58LwDX8pATS7fD4/pmyZ5LwKwjAnaZmTfiBpRJKCNN0/3+cBAAAAMDOKMiwWTpKZWRKGYSUIApNU7CqVtvsOBqA1hGF4ckFeNwgKkpalYfjKL/r6+hfiHAAAAADOD0UZ2tLpm/vEcawoiixJU3t/eDjf09PTOTkxEfrMB6B17Nq164J8Lrdg+4VVKpVgoV4bAAAAQH0oytCWTHK19VHmpNhJqdI0Xy6XhwcHB9/0mw5AS3DumMxeeG9kpCGnyyfJRENOBAAAAOCsKMrQrqZnf1gYhlMdhULFgqDTme3pKpXu8poMQEt47PDhXkl3+s4BAAAAoHEoytCW4iSxuHZly6NHjxaqUWRTk5MstwQAAAAAAGdFUYa24ZyTc9lEsmoUWbVWlPX29hZ7paLPbAAAAAAAoPlRlKFdODNzJpkkKxaLI2EYmqRVTrq7q1R6fvo4jxkBAAAAAEAToyhDW6pGkZ2MokDKNvYXBRkAAAAAAJjBjEXZ9FI2oGk5pyRJrFwuW1Rbbvnkli2rPKcCAAAAAAAt5lxFmTOzhgUB5sKkuFAoTAVmBUlFJ909/ZgLwz0eowEAAAAAgBYx04wyppOhFZjMkjAMK0EQmKRiV6m03XcoAAAAAADQWs4oypI0tSiKfGQBZiVJU3t/eDi/c+fOzqnJydB3HgAAAAAA0Jo+UpQFQVD1FQSYg3y5UhmaOnr0Td9BAAAAAABA6/pIUZbL5aZ8BQFmy6ROSXu6SqW1vrMAAAAAAIDWlZOkoeHhvO8gwGxNlssstwQAAAAAAHOWk6T+vr5if19f0XcYAAAAAAAAwJf/A6835sk2RkYNAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "fa93a4de",
   "metadata": {},
   "source": [
    "## 3. Detecting drift\n",
    "\n",
    "As you just did, if labels (e.g. `y_`) are available at test time, harmful drift can be detected by simply evaluating your ML model's performance. However, in deployment/production, labels often aren't available. What can we do to detect drift here?\n",
    "\n",
    "In the following sections you will use a number of different [Alibi Detect](https://docs.seldon.io/projects/alibi-detect/en/latest/) drift detectors to detect different types of drift. As you will discover, the general workflow is usually similar, with a preprocessing step, a statistical hypothesis test, and then a comparision of the resulting p-value to a certain threshold to decide if drift has occured.\n",
    "\n",
    "![drift_detection_pipeline.png](attachment:drift_detection_pipeline.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8700f882",
   "metadata": {
    "id": "a51888f2"
   },
   "source": [
    "### Detecting drift on the inputs\n",
    "\n",
    "The first type of drift we'll consider is input (covariate) drift. This is defined as $P(\\mathbf{X}) \\neq P(\\mathbf{X}_{ref})$, i.e. input drift is said to have occured when the distribution of the input data has changed.\n",
    "\n",
    "We can't apply a drift detector directly to the original data `x_`. Instead, we monitor the data in the embedding space. To start, we need a preprocessing function which will be applied to reference data (the data given to a detector at instantiation time) and test data (the data given to the detector's `.predict` method).\n",
    "\n",
    "This function should: \n",
    "- Accept a list of posts `x`.\n",
    "- Apply the classifier embedding model to transform `x` to the embedding space.\n",
    "- Return the embeddings as a ndarray.\n",
    "\n",
    "*Hint: the embedding model will return a `torch.Tensor`, but Alibi Detect expects `np.ndarray`'s, so you'll need to convert to this with `.cpu().numpy()`*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92ee2758",
   "metadata": {
    "id": "2bbf345d"
   },
   "outputs": [],
   "source": [
    "def preprocess_fn(x: list):\n",
    "    x = ???\n",
    "    return ???"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99b2e917",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    def preprocess_fn(x: list):\n",
    "        x = clf.embedding_model(x)\n",
    "        return x.cpu().numpy()\n",
    "\n",
    "</details>\n",
    "\n",
    "This is a good example of how you can often recycle components of your existing ML models to serve as the preprocessing stage of a drift detector. Now that you've defined the preprocessing stage, you can instantiate the actual drift detector.\n",
    "\n",
    "Import the `MMDDrift` detector from `alibi_detect.cd`, and set the following args and kwargs (see the [MMDDrift docs](https://docs.seldon.io/projects/alibi-detect/en/latest/cd/methods/mmddrift.html) for a full list of kwargs):\n",
    "\n",
    "- The first arg is expected to be a np.ndarray or list of reference data. This is the data we will compare to at predict time i.e. usually it will be from the same data distribution used to train the model. As a general rule, you shouldn't use the exact model training data as reference data, so instead, use the held-out `x_test` data.\n",
    "- Set the `backend` kwarg to `'pytorch'`, since our classifier is a PyTorch model (the default is `'tensorflow'`).\n",
    "- Set `p_val` to a decimal of your choice (the default is `0.05`). This is the p-value threshold below which we flag drift. \n",
    "- Pass your preprocess function defined above to `preprocess_fn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfe09d4e",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "db4aea0e",
    "outputId": "2d098717-e6c3-4168-cc6e-b88d933510ee"
   },
   "outputs": [],
   "source": [
    "from ??? import ???\n",
    "x_ref = ???\n",
    "dd = MMDDrift(x_ref, backend=???, p_val=???, preprocess_fn=???)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58ff7413",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    from alibi_detect.cd import MMDDrift\n",
    "    x_ref = x_test\n",
    "    dd = MMDDrift(x_ref, backend='pytorch', p_val=.05, preprocess_fn=preprocess_fn)\n",
    "\n",
    "</details>\n",
    "\n",
    "Now the detector is ready to be used! Pass the `x_nodrift` data to `.predict()` to check if this data has drifted from `x_ref`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9784fe85",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "80db1111",
    "outputId": "f7c99ed3-82ef-49ef-f807-97a384080541"
   },
   "outputs": [],
   "source": [
    "dd.??(??)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9177a685",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    dd.predict(x_nodrift)\n",
    "\n",
    "</details>\n",
    "\n",
    "The `predict` method returns a dictionary with the relevent predict info stored under `'data'`. The important field for us is `is_drift`, which equals `0` if no drift is detected (i.e. `'p_val'` $>$ `'threshold'`), and equals `1` if drift is detected (i.e. `'p_val'` $\\le$ `'threshold'`).\n",
    "\n",
    "Hopefully drift isn't detected above (`is_drift = 0`), meaning we have failed to reject $H_0: P(\\mathbf{X}) = P(\\mathbf{X}_{ref})$. Note that you might flag drift here sometimes, as a perfectly calibrated detector will still return the occasional false positive (we discuss this [later](#Important-note!)).\n",
    "\n",
    "Now, try the same for `x_drift`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28ad2f46",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "57005484",
    "outputId": "5f955430-ae9e-4878-b659-021f25a8c534"
   },
   "outputs": [],
   "source": [
    "dd.??(??)['data']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc31733b",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    dd.predict(x_drift)['data']\n",
    "\n",
    "</details>\n",
    "\n",
    "\n",
    "Hopefully, drift is detected! If it is, this means there is enough evidence to reject $H_0$ and flag drift (although sometimes it might be missed, unless the detector's [test power](https://en.wikipedia.org/wiki/Power_of_a_test) is 1)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe533d5b",
   "metadata": {},
   "source": [
    "### Detecting drift on the outputs (optional!)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8cc5f8f",
   "metadata": {},
   "source": [
    "The next type of drift we'll consider is output (prior) drift. This is defined as $P(\\mathbf{Y}) \\neq P(\\mathbf{Y}_{ref})$, i.e. output drift is said to have occured when the distribution of the output data (labels) has changed.\n",
    "\n",
    "If you have access to labels `y_` you can monitor output drift by passing them directly to the detector (you can try this!). When labels are not available, the ML model's outputs can be used as a proxy for the true labels, as is done below:\n",
    "\n",
    "*Note: You can either use the classifier's logits output, or its class predictions.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ca13c0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess_fn(x: list):\n",
    "    x = clf(x)\n",
    "    return x.detach().cpu().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61e38077",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd = MMDDrift(x_ref, backend='pytorch', p_val=.05, preprocess_fn=preprocess_fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67817bc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd.predict(x_nodrift)['data']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63a47a67",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd.predict(x_drift)['data']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dea12f23",
   "metadata": {
    "id": "7ad48493"
   },
   "source": [
    "### Detecting model drift"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32bb22ad",
   "metadata": {
    "id": "4a0e48b9"
   },
   "source": [
    "A downside (or upside!) of monitoring input/output drift is it doesn't tell you if the drift is likely to be harmful to your ML model. An alternative approach is to use a model uncertainty based detector such as the [ClassifierUncertaintyDrift](https://docs.seldon.io/projects/alibi-detect/en/latest/cd/methods/modeluncdrift.html) detector. This tests for a change in the number of instances falling into regions of the input space on which the model is uncertain in its predictions. \n",
    "\n",
    "Define your `ClassifierUncertaintyDrift` detector with the following args and kwargs:\n",
    "\n",
    "- Again, the first arg should be the reference data.\n",
    "- The second arg should be the ml model used to obtain uncertainties. Since our classifier provides a notion of uncertainty (the logits), this can be used directly.\n",
    "- Set the `backend` and `p_val` kwargs.\n",
    "- Check the [ClassifierUncertaintyDrift](https://docs.seldon.io/projects/alibi-detect/en/latest/cd/methods/modeluncdrift.html) docs to see what to set `preds_type` to.\n",
    "- Set `preprocess_batch_fn=clf.embedding_model`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15c950d8",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "0b17b8f3",
    "outputId": "6ea4cbcd-fc19-4194-ea98-eb2e37005d31"
   },
   "outputs": [],
   "source": [
    "from alibi_detect.cd import ???\n",
    "\n",
    "\n",
    "dd = ClassifierUncertaintyDrift(???, ???, \n",
    "                                preprocess_batch_fn=???, backend='pytorch', \n",
    "                                p_val=.05, preds_type=???)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc472513",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    from alibi_detect.cd import ClassifierUncertaintyDrift\n",
    "\n",
    "    dd = ClassifierUncertaintyDrift(x_ref, clf, \n",
    "                                    preprocess_batch_fn=clf.embedding_model, backend='pytorch', \n",
    "                                    p_val=.05, preds_type='logits')\n",
    "\n",
    "</details>\n",
    "\n",
    "Once instantiated, you should be able to make predictions on `x_nodrift` and `x_drift`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09582f60",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "26423a3f",
    "outputId": "ff80f477-cf5f-407f-ae20-882c7fc4876b"
   },
   "outputs": [],
   "source": [
    "dd.predict(x_nodrift)['data']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ced5426",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "005cff42",
    "outputId": "c0913742-0c3d-4761-849b-b9ff5b27ef46"
   },
   "outputs": [],
   "source": [
    "dd.predict(x_drift)['data']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "162d36a7",
   "metadata": {
    "id": "a6c7e40a"
   },
   "source": [
    "## 4. Accounting for context\n",
    "\n",
    "So far, **we have assumed that both the reference and test data are [i.i.d.](https://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables) samples from their underlying distributions**. However, in some *real world* scenarios this assumption doesn't hold. \n",
    "\n",
    "For example, consider the case where we have trained our classifier on a training set consisting of all twenty topics (to ensure the classifier sees a diverse range of posts during training). Now, imagine that we wish to use the classifier whilst the World Series and Stanley Cup are occuring, leading to a spike in news posts on baseball and hockey respectively. Furthermore, there is not too much news on Mac or Windows since there are no new releases or products planned anytime soon. Crucially, the distribution underlying each subpopulation (e.g. the distribution of *hockey* news itself) remains unchanged, therefore the classifier should still work well in this scenario. However, since the relative frequency of one or more subpopulations (i.e. news topics) has changed, if we assume the test instances are i.i.d (i.e. we use the standard `MMDDrift` detector), then we will detect drift.\n",
    "\n",
    "To ensure that drift isn't flagged when the distribution underlying each subpopulation is unchanged, we need a way to provide our drift detector with relevent context (*conditioning*) information.\n",
    "\n",
    "### Generate some example data\n",
    "\n",
    "We can use the `split_data` to generate some data matching the scenario described above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a99c046",
   "metadata": {
    "id": "43ee261b"
   },
   "outputs": [],
   "source": [
    "n_classes = len(classes)\n",
    "n_nochange_c = 1000 // n_classes  # equally subsample each class from 1000 instances\n",
    "\n",
    "n_change_c = [50] * n_classes  # 50 of each class (so 1000 in total), but then mod. below.\n",
    "n_change_c[4], n_change_c[5] = 5, 5  # few stories on Mac/Windows\n",
    "n_change_c[9], n_change_c[10] = 95, 95  # more stories on baseball/hockey\n",
    "\n",
    "(x_nochange, y_nochange), (x_change, y_change), (x_train, y_train) = split_data(dataset.data, dataset.target, \n",
    "                                                                      n_nochange_c, n_change_c, seed=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "144e7d50",
   "metadata": {},
   "source": [
    "Three sets of data have been generated:\n",
    "\n",
    "- `_nochange`: Data equally subsampled from all 20 subpopulations (topics) - *Used as reference data*.\n",
    "- `_change`: Data with the prevalence of some subpopulations/topics changed - *We'll monitor this for drift*.\n",
    "- `_train`: Held-out data with the same subpopulation prevalence as `_nochange` - *Used for training the classifier*.\n",
    "\n",
    "Visualise the `no_change` embeddings. You should now see all clusters for all 20 topics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cc9cc0d",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 607
    },
    "id": "6e7f68ad",
    "outputId": "dd0cec8f-eb4f-4215-b77f-9b1af0ce91a7"
   },
   "outputs": [],
   "source": [
    "emb_nochange = emb_model(x_nochange)\n",
    "plot_clusters(emb_nochange, y_nochange, classes, dr_model=umap_model, title='Reference data: clustered news topics')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64bf7c3b",
   "metadata": {
    "id": "7aa6b2c9"
   },
   "source": [
    "### Vanilla MMD detector"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dd9fb25",
   "metadata": {},
   "source": [
    "Before looking at context-aware drift detection, lets see how the standard `MMDDrift` detector performs here. To save a bit of time, the embeddings are pre-computed and passed straight to the detector (so that `preprocess_fn` doesn't have to be called for both detectors). \n",
    "\n",
    "For the reference data `emb_ref`, the `emb_nochange` embeddings computed previously are used. To save time, the ndarray `emb_change` is also computed here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83532fbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "emb_ref = emb_nochange.cpu().numpy()\n",
    "emb_change = emb_model(x_change).cpu().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c1b3cf1",
   "metadata": {},
   "source": [
    "Now init an `MMDDrift` detector, and use it to check for drift on `emb_change`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa510b7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd = MMDDrift(emb_ref, p_val=0.05, backend='pytorch', n_permutations=100)\n",
    "dd.predict(emb_change)['data']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59f93eec",
   "metadata": {
    "id": "520b3e80"
   },
   "source": [
    "You should find that drift is consistently flagged. This is expected since the vanilla MMD detector cannot take any external context into account and correctly detects that the reference and test data do not follow the same underlying distribution. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69b485ff",
   "metadata": {},
   "source": [
    "However, since the underlying distributions of the subpopulations are unchanged, a classifier trained on the full `x_train` should still work well on `x_change` (relative to `x_nochange`). This can be seen below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dfd02cc",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "f7f21ff6",
    "outputId": "cd98a551-6ec3-4e46-fdcb-040a99a13136"
   },
   "outputs": [],
   "source": [
    "n_train = 2000\n",
    "idx_train = np.random.choice(len(x_train), size=n_train, replace=False)\n",
    "x_train, y_train = [x_train[_] for _ in idx_train], y_train[idx_train]\n",
    "\n",
    "TRAIN_CLF = False  # Set to TRUE to train classifier, otherwise it will be loaded from disk\n",
    "filepath = 'classifier_full'\n",
    "if TRAIN_CLF:\n",
    "    # init model\n",
    "    clf_full = Classifier().to(device)\n",
    "\n",
    "    # Train model\n",
    "    train_model(clf_full, x_train, y_train, epochs=5, shuffle=True)\n",
    "    clf_full.eval()\n",
    "    \n",
    "    # Save model\n",
    "    torch.save(clf_full.state_dict(), filepath)\n",
    "else:\n",
    "    # Load model\n",
    "    clf_full = Classifier()\n",
    "    clf_full.load_state_dict(torch.load(filepath, map_location=device))\n",
    "    clf_full = clf_full.to(device)\n",
    "    clf_full.eval()\n",
    "\n",
    "_, _ = eval_model(clf_full, x_nochange, y_nochange)\n",
    "_, _ = eval_model(clf_full, x_change, y_change)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2df2e754",
   "metadata": {},
   "source": [
    "You should see a similar classification accuracy for the `_nochange` and `_change` cases above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70bc68e9",
   "metadata": {
    "id": "58c29960"
   },
   "source": [
    "### Context aware MMD detector\n",
    "\n",
    "We want to monitor whether the underlying distribution of each subpopulation has drifted, given the topic label of that subpopulation. In other words, **we wish to provide the posts' topics as context**. Since we don't have access to the true labels at test time, we instead **condition on the prediction probabilities of the classifier**. \n",
    "\n",
    "First, define a function `context` to:\n",
    "\n",
    "- Take the raw embeddings `emb`. \n",
    "- Compute the `logits` with the classifier (*Hint: We need the classifier to be reasonably accurate, so remember to use the classifier trained on all 20 newsgroups*)\n",
    "- Take the `logits`, and return an *ndarray* of prediction probabilities. (*Hint: you'll need the [softmax](https://pytorch.org/docs/stable/generated/torch.nn.Softmax.html) function for this.*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23eca6c3",
   "metadata": {
    "id": "b5b5f6a1"
   },
   "outputs": [],
   "source": [
    "def context(???):\n",
    "    \"\"\" Condition on classifier prediction probabilities. \"\"\"\n",
    "    logits = ???(???).detach()  # Important classifier trained on full dataset used here!\n",
    "    softmax_fn = torch.nn.Softmax(dim=-1)\n",
    "    return softmax_fn(???).cpu().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdb582a9",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    def context(emb):\n",
    "        \"\"\" Condition on classifier prediction probabilities. \"\"\"\n",
    "        logits = clf_full(emb).detach()  # Important classifier trained on full dataset used here!\n",
    "        softmax_fn = torch.nn.Softmax(dim=-1)\n",
    "        return softmax_fn(logits).cpu().numpy()\n",
    "\n",
    "</details>\n",
    "\n",
    "Now you're ready to define your [ContextMMDDrift](https://docs.seldon.io/projects/alibi-detect/en/latest/cd/methods/contextmmddrift.html) detector. Define it with the following args/kwargs:\n",
    "\n",
    "- arg 1: The reference data, as usual. You can use `emb_ref` for this.\n",
    "- arg 2: The context associated with the reference data. For this you can use your `context` function to transform `emb_ref`.\n",
    "- `p_val` and `backend` kwargs as usual.\n",
    "- Set `n_permutations=100` to match the `MMDDrift` detector's default (this usually defaults to 1000 for the `ContextMMDDrift`). This kwarg's controls how many permutations are used to compute the p-value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90a5cb45",
   "metadata": {
    "id": "60c8f1ca"
   },
   "outputs": [],
   "source": [
    "from alibi_detect.cd import ContextMMDDrift\n",
    "dd_cad = ContextMMDDrift(???, ???, p_val=.05, n_permutations=100, backend='pytorch')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44de12a8",
   "metadata": {},
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    from alibi_detect.cd import ContextMMDDrift\n",
    "    dd_cad = ContextMMDDrift(emb_ref, context(emb_ref), p_val=.05, n_permutations=100, backend='pytorch')\n",
    "\n",
    "</details>\n",
    "\n",
    "Now you're ready to use the detector! Try it out on the `emb_nochange` and `emb_change` data sets (*Hint: You must give the `ContextMMDDrift.predict()` method context associated with the test data*):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40a9d872",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd_cad.predict(???, ???)['data']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "201b3c4f",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "c79c5ea2",
    "outputId": "3d0be038-5983-4188-97fa-dab17ef4f89d"
   },
   "outputs": [],
   "source": [
    "dd_cad.predict(???, ???)['data']"
   ]
  },
  {
   "attachments": {
    "contextMMD_calibration.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEWCAYAAACwtjr+AAAgAElEQVR4Ae2dCbxdRX3Hv29/JGEJENkSwi4gyBK2qmArImCrSEXA1gooxQWsiFLXKoKtImqtFKmoUKsVBAoWrUipQKWKLUF2MLIIJIAQIATI9tZ+ficzMDm5y3nv5b537r2/fz4357xzZubMfGfO/Oc/M2cGLCZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiZgAiUj8M/A5yYhTtsAo0B3eNbVwHHr6LkHAguSsB4CXp/8PdHTu4E/nGggFfxvBvwceB74coX7U3XpDOB7U/XwMTy3UflSKwoqwzvUcjAF91Q2F03guf8E/M0E/Jfda77uKU18Y2VYmgi1aEQOL5guvdw7AvfXcH8j8PIa98dyS8pXL+6nEk+vSM7X5elJwFPABkERr8uw2yGsRuVLK7M7HjgReE2SyPcm5z6dRAKdk/gsP2riBJq5cTAXuGeciqaZ0z3xXHcIrUagLctzWZSNuoI+HiqjJcBFQH+VEnYv8CfJPWXcYmDvcO0y4PfA0tBtU61FqFbP/yTh6DTtNugDvgQ8AjwByPxeL+c+/tkV3Krl/iDwx/FGON4QWlj6U90S/x3iJ/c/CG7UxSS5HXgBOCZ0Z8ny+GhIk7hU6kbYtwq7WmmUpfHnwF+H5/0oPD/tlhODrwKPhZ/OdU0S4/Fh4EngceCEcC9/kAWlbsT4LHX7FQk7TXc+zO2B64Cng8X0r8BGeUfJ3yoH1wLPhPz8RHKvF/iX0MWn7qp9knsfAx4I96Qsj0zuRb4qJyq3vwNSK3bbpOvwv4Dzcl12BwC/BJ4N+V6r+zLNF3X9XVojzkkU6QD+PuTRc8CdwG7BQVoudSmmJ/X/xlCmVVbPAWKdIbe/AP4xlOXfAAcnHrcErgq8Zan/ZXJP8b88lH11q/4a2CO5n76Hulyr+7ta/uwS3tk/COVbjCuFpXgpfioXiq/iHUXxkCV0X8gj5Z94VpJ6aVL+qTzfASwLXezV8l/v/vzcQz4U4qfLql9uBZSfCwE9u5psCHw7vJ+PhmEE1VeSmN/Vyu/GoS7W+6/y/cPgTwfVwbcFLirDr0zulfpUGXEXMAdQAlWIq42tfBpQxRJF4KWAorwLWD+pzAQkSlpoI+h4T8e0kOsFVeFTfBSeKuPPp46TcxVIvWwx/teHsGILJn2pLwY+GV5aKdTUxE+fr+BV+QwBZ4f0SNnFSj4+vha7emlMeaThxTGgM4FfAS8DZoWK8azgMMZNbnoAVUrLgZkxoNwx/6wiYafpzgWXKe1DAhfFTcpayrCSKP+kDKUYxVx/7x8c6kVdGeKvl1B5rDRHeVuogFTJqhJQRbFFuCm+g6Eild/3BaUcK6SbQiNEykz5rMohjg9tFRSluClspUWKU2mpJHllUyvOqf9DgVuCIla8VAnH+KflUn7y5UXlUWVZ78DWwG+TRpPcqmyqElT+i40aeHIrUX58PfDeMzQIXxfuibm4HRX8fiQoaoUjyb8HadnJl/96+ZNvUKZhKT5SomqoqvFzboh3iEYWjx8Hdkq/GrWHxZu5Y700Kf9UF6mO0HtcK/+nhcaNutSj3AwcG/4Qg91DuVElr8bwW8K9/JjNlcA3gOnhPf4/4D3Bbb3y+x+hQaB3Wnnz2uBvr9B40Tukcq+GpNIXG6LBWTkPimjal6oXUK3JSiLLQK0hZYhEikcKqJKopauCK+0uSQta/sXS/VjI9VKqUlHrOYpaSGq5VhK1sNP4v6GGslEL+gJgdoWA4vPjLRWqgZyVl3/ZarGrlUY9I+URn6nworJRHigvoqji0n2J4rEimQSha7Jw1FqrJPln1Qs7n+5KYabX9LKptVdJ3l7jnioJWR1Rdg3pin/nj6owjggXxTcdX1OZVB5uHipnVcaxnMqLFE1UNmrlfjcX+DU1JpKk+TKWOKtClZJQvkSrJD62iLJJK9f3Az9L0q4Wb1SsuqyK7C9ChToclHp8lpS4yoBE8U8VuuKlxoAmv0jy70FadvLlP3h58ZDPn1rKRi3+L77oE2YEJagKW6J4pI1BWZOypCpJvTQp/9QQjlIv/1VOYr0mpZPWeTGMeFQjS41jSapsNClnVa5HRu+CGhCSWuVXDZKRKo3H84HY6AxBZZOWojKK19Y45gvfGjcn+Q+Zg1EeTsxZzeRSt5J+6vbRyy1L5k3hRX4z8P3gUVr2C0FRqRUZK8ZNY8AFj2pdqpJQi1Dmt34/rdHqlOmdj3+1R6krSS+oXkx12aQFsJIftabUiq0l+WenXQG1/NW7p3CUF1HSfNE1tcRVoUaRZaMXtojUC7teuvUiXQKoayBaDNXyWa3Jao0XxVXdrlGUBlk/0Sp9Z9JdoHKgLqj0OXm/CkcMlD51zSi8KGk+aQxLrfJYvnRUxRatjuin2jH/3DTOqR81hNTVpS4gNQbU0NEkjaKSxjmf/2KvCjlKvB/Trgoyiu6pNR8lDVeVmrqL5W+sUi9/aoWXL4OqY1Sm03jmOdcq3/XSlN6vl/+q06QYJH8WurBiWZJFIYWhd0TWpBq6aZkM3tAzZJFIkcdyJitHPRVR8unTdaVR74zKr7rP8qJw1UsQw9RR7mvmX5mUjSIbRSarWk0S9YEr8frF7jN1RSkj1MJUP3psXSpTdE0tc1kzsYWStr5CsJnlkrY61RqNItNarXb188s60k/hVStoysx8/GNY+aMyV/3EyhiZs+pqqDW9NH2Z82HFv/PPjuxknVVLo/zWC1vhqGBFSfMlXhvvsV7Y9eL2dyH+6k5Q5fmOXCs7jZde8u3SCwXPlfZvAqcAm4RyoO7eSuUpH6TKhLqUUv5pPilOsmxi+dJRXR1qLK1r+RowD5DVthNwenhAvfIhZ2mc8/mvSjllEe8rb2P3c0yL7kk5RUnDVT0kSz+WW1WqKbf03Yz+dayXP/XKUL4Mir/yOY1n+rx657XSJL9pfOrlv8YX1ehVF6TqutigVjg6Vxe/nqd6SePJaT7EeOoZsmykiGI507tSbRw7+tNRfpWH8pcX3fvbJEy5UX6pXq4qZVI2J4cCpwRqTCMOnFeKvFq06qpSH3maCeqLF1y1TpR4VUjVRAPxgq7MVKtQZnAUtbRUycg0ja0AvVjqRqokMq//KsRf/ZvVTG35VWs2dqGp1aACqOdJ1Pc6nkqxGrtaaSzyPBUeTYtWoVeBlVkfu4FWx3j8/080bOW1WqJq2SlvYgVaKUbqd5fFcGroV07HbCq5j9dU+Sh/1IKUaAJEHFwPl6oe1JLXIK/KlcZs1A0razyKOOpvlSlZ5CqD6iKKZSO6m+hRk0fUElYLV8pFVnIsb+py+tPwrqjB8+4KDxNXlWlVbB/MvZd6N1TuFbbKtcaDfhIqKg0aq+tM6dK4gsJOy46Un54tC1L5ovc2dq0pXmo4iou68ap1z9TLH71P4in+lURlUHmqOkDjDaov/jfpEankp9a1WmnK+6uX/xrT0mQnTcpQnSjlE0XlV1aH8nK/wCreS49q8Pxn+K5NSkb1vYYGqvHM+1WvkhrDcczmoOBAdaOsKZUrKTnlg8bOFa+qUiZlI6UhMJrNpS6PahMElBhB1ODrq3KFX+MhesnVMpHFEwtvJQDqx9YgtfrrNdsk37erPlVZTApD3TRyV+37FsFXf7sqd82suaLSA8M1vfwq0Koo1TrRC6w0S1QxfSeYp0eHa0UO1djVS6P6rNXalRmczjSJz1QeqMLUDBrNYlLaauVL9FfkONGwPxsGdqVsNJBZi7m6czQAr8pdlqXy+48KRFJlSB+gqqyp4pIVpckrRUXdvlIyavwovWpAqVKVqHUoK1yz4qTM9Lcq9nX9TqqSUflUw0bvhuKiCkyixpTGxpQ2lbvYcxBuZ4d/D93JUgDirDITReVY4wnqCVBLVwP+Cl+i1rh6FmQ9aJD6M7mxMYWrSQWKl8Z5pHhUwUr0TiivVC7FsFLZlLt6+aMuRHVVK88Vx7zondYHnv8W6hRVxHEQPu+2yN+10pT3XyT/9V6rl0ZKJ+2u1tiZ6i6VazUA1ditJupmlLIVK7HWLMCiXbXKF+WJJj+pC1aNAonqBPXOqHtWYaqe1PhPU0g6+NkUEXYkTWAcBKRspCRbQSpNPimaLjWqUiunqL8yu2vFNK1T3uu6FbVOI+fATKDJCciKVWtZ75m6g2TJVGulN3lSHX0TqE2g0cpGL5jW8ZKZVWkcQ/1+6p7RYLlMenXpRNFHnvIn/9XGSqJbH02gjAQ0sK3pxeoy1SC9xhirTc8uY/wdJxNoCgIa3NPYiwa81Weo8YxUmSgR6k+OoinMml4skTu516CdvsJWOArPYgImYAIm0IQEGmnZaJaELBMNfmsQUjPI1I2Qigbeo8SZJfpb7uReg6n6kFLhKDyLCZiACZhAExKIH641IuqajqoZF1H00ZamyuVF03ZPC9ZPXM5CftOZZPKra3nR+l76MX369Hk777xz/r7/NgETMAETqEHglltu0Uy9assk1fA5tluNVDZFY6Ivm/XTvHp90zGWfV/0NbR+7LzzzqPz5+fXrisaBbszARMwgfYk0NHRoSnxDZdGdqPpW5f0i1p9XFXry1x1m8XF5Mbqt+Gg/AATMAETMIHxE2iksrk5fPClAX5NENDHUvqIMZV0VVN9gaqP7SRyJ/dxgoDcaS0xiwmYgAmYQBMSaGQ3mr541ZpS+rJeM8kuDF/z6stX9XdJoei+vpDVV6r6EjV2oemrX30Vq69eFY7GdbSKrMUETMAETKAJCVRavK0JkwHz5s3zmE1T5pwjbQJjJzA4OMiiRYtYubLeguhjD7tVffT39zN79mx6erSU3UvS0dGh1e3TDQNfurkOzxpp2azDaDooEzABE3iJgBTN+uuvzzbbbENHR8u0mV9K4Do+Gx0d5emnn84U9LbbamRj8qWRYzaTnxo/0QRMoC0IyKLZZJNNrGgK5rYUsnhNpSVoZVMws+zMBEygXARs0YwtP6aal5XN2PLLrk3ABEzABMZBwMpmHNDsxQRMoFwE5oSxG7Xe18VP4RWR3//+9xx77LFsv/32mqTEG9/4Rn77W20jBV/96lfRoPzSpdpyqbJozOnAAw9c4+aee+7Jbrut3qPvhhtuyNLzrW9960U3t912W3btS1/6Unbt+OOPR+Mwe+yxBzvttBPvfOc7s7GZFz2U5MQTBEqSEY6GCZjA+AksevhhLlqiPRXXjZwws/7+Yhp0P/LIIznuuOO45BJ9kw633347TzzxRFbpX3zxxey7775cccUVnHCCNgStLM8//zwLFy5kzpw53HvvvWs5kuK59NJLOfHEE7N7CleKJZVzzjmHo446CsVJSu51r3sdd911F7291TYpTX1Pzrktm8nh7KeYgAm0GIHrr78+m0b83vdqp5TVIiUgS+WBBx7ghRde4HOf+xxSDrXk6KOP5gc/0L56ZG7f/nZtcvqSzJ07NxvYlxKTMvnpT3/K4Ycf/pKD5ExW3Yc+9CE233xzrr5auzqXR6xsypMXjokJmEATEZDloK6zSiJLR91rUjwLFizIrJ1K7nTtrW99a2b96PxHP/oRb3qTdsReU2S1XHbZZfzyl79k7733pq9Pi6tUF7n5zW+0m3N5xMqmPHnhmJiACbQIAVkzUjadnZ2ZMpGiqCaakjxz5sysK26XXXZh2rRpazmV9aMwFG7e8lnLMWQWUKXrU3nNymYq6fvZJmACTUvgFa94Bbfcoo/v15Q777yT++67j0MOOST76FRWjpTE8PAwGvzX79Of/vQano455hhOPvnkqopE3WL68v/aa6/l4IMPXsNvpT9uvfVWpLjKJJ4gUKbccFxMwASahoAG4T/xiU9wwQUXcNJJ2bZa3HHHHXzwgx/kjDPO4OMf1872q0WzxbTqgWaSVRJNNHj88cc59NBDeeyxxyo54cwzz+TJJ5+kq6v6psUa0zn33HOzsA477LCK4UzVRSubqSLv55qACawzArPnzqXIDLKiD1R49USD8VdeeSWnnnoqZ599djbNWVOZNV35/PPPX8O7lIksnI9+9KNrXI9/aOmdaveim1e96lXxdK3j6aefzllnncXy5cs54IAD0OSFMs1EU4RbZlEhL8S5VvnzBRNoWQKaIly2bqJmgF2J22QtxOkxm2YoIY6jCZiACTQ5ASubJs9AR98ETMAEmoGAlU0z5JLjaAImsBYBDYZbihOYal5WNsXzyi5NwARKQkBrjml/lqmuQEuCo240xEm8xG2qxLPRpoq8n2sCJjBuAtpxUlOJFy9ePO4w2s1j3KlzqtJtZTNV5P1cEzCBcRPQB45TtePkuCPd5h7djdbmBcDJNwETMIHJIGBlMxmU/QwTMAETaHMCVjZtXgCcfBMwAROYDAJWNpNB2c8wARMwgTYnYGXT5gXAyTcBEzCBySBgZTMZlP0MEzABE2hzAlY2bV4AnHwTMAETmAwCVjaTQdnPMAETMIE2J9BoZaPdexYA9wMfq8D6NOAe4A7gZ0C6icQwoJ2G9Luqgl9fMgETMAETaBICjVxBQNvJnQccAiwCbg5KQ8olyq3APsBy4H3AF4Fjws0VwJ7RoY8mYAImYALNS6CRls1+waJ5EBgALgGOyKG6PigaXf4VMDt333+agAmYgAm0AIFGKputgIUJI1k3ulZN3g1cndzU8qTzgxJ6S3LdpyZgAiZgAk1GoJHdaGNB8Y7QnfbaxJPGbx4FtgOuA+4EHkju6/Sk8PPqrzkw/tMETMAEykSgkZaNFMWcJLHqItO1vLwe+CTwZmBVcjO6VTfcDcBeyb14ekFQUvvMmjUrXvPRBEzABEygZAQaqWw0IWBHYFugFzi2wqwyKZBvBEXzZMJmJtAX/t4UeHWYtZY48akJmIAJmECzEGhkN9oQcApwDaCZaRcCdwNnhrEYTWc+B5gBXBaAPRIUzy5BCY0AUohfsLJpliLleJqACZjA2gQ61r7UnFfmzZs3On++5hNYTMAETMAEihLo6Oi4JQxHFPUyLneN7EYbV4TsyQRMwARMoPUIWNm0Xp46RSZgAiZQOgJWNqXLEkfIBEzABFqPgJVN6+WpU2QCJmACpSNgZVO6LHGETMAETKD1CFjZtF6eOkUmYAImUDoCVjalyxJHyARMwARaj4CVTevlqVNkAiZgAqUjYGVTuixxhEzABEyg9QhY2bRenjpFJmACJlA6AlY2pcsSR8gETMAEWo+AlU3r5alTZAImYAKlI2BlU7oscYRMwARMoPUIWNm0Xp46RSZgAiZQOgJWNqXLEkfIBEzABFqPgJVN6+WpU2QCJmACpSNgZVO6LHGETMAETKD1CFjZtF6eOkUmYAImUDoCVjalyxJHyARMwARaj4CVTevlqVNkAiZgAqUjYGVTuixxhEzABEyg9QhY2bRenjpFJmACJlA6AlY2pcsSR8gETMAEWo+AlU3r5alTZAImYAKlI2BlU7oscYRMwARMoPUIWNm0Xp46RSZgAiZQOgJWNqXLEkfIBEzABFqPgJVN6+WpU2QCJmACpSPQaGVzGLAAuB/4WIXUnwbcA9wB/AyYm7g5Drgv/HRuMQETMAETaFICjVQ2XcB5wOHArsDbwzFFdSuwD/BK4HLgi+HmxsBngP2B/cL5zNSjz03ABEzABJqHQCOVjZSELJoHgQHgEuCIHJrrgeXh2q+A2eH8UOBa4BlgSTiXlWQxARMwARNoQgKNVDZbAQsTJosAXasm7wauDjeL+j0JmK/f4sWLq4Xr6yZgAiZgAlNMoHuKnx8f/47QnfbaeKHg8QJAP2bNmjVa0I+dmYAJmIAJTDKBRlo2jwJzkvSoi0zX8vJ64JPAm4FV4WZRv/mw/LcJmIAJmEAJCTRS2dwM7AhsC/QCxwJX5RjsBXwjKJonk3vXAG8ANClAP53rmsUETMAETKAJCTSyG20IOCUoCc1MuxC4GzgzjLNI8ZwDzAAuC+weCYpHEwPOAqSwJPKjaxYTMAETMIEmJNDRhHGuGOV58+aNzp+vuQIWEzABEzCBogQ6OjpuCWPmRb2My10ju9HGFSF7MgETMAETaD0CVjatl6dOkQmYgAmUjoCVTemyxBEyARMwgdYjYGXTennqFJmACZhA6QhY2ZQuSxwhEzABE2g9AlY2rZenTpEJmIAJlI5AUWXzJqCo29Il0hEyARMwAROYWgJFFcgxYV8ZbQGw89RG2U83ARMwARNoNgJFlY0WytTSMg8A/wzcBGjF5fWbLcGOrwmYgAmYwOQTKKpsFLPnwgZn2pdmC+BI4NfAByY/2n6iCZiACZhAMxEoqmy06dmVwA1AT9g9Uztw7gF8uJkS7LiagAmYgAlMPoGiC3H+KfD3wM9zUdQum9r0zGICJmACJmACVQkUtWx+X0HRnB1C/VnV0H3DBEzABEzABMYwnfmQCrTUjWYxARMwARMwgboE6nWjvQ94P7A9cEcSmmah/SL526cmYAImYAImUJVAPWXzfeBq4PPAx5JQnvdmZgkNn5qACZiACdQkUE/ZjAIPASdXCGVjK5wKVHzJBEzABExgLQL1lI0smz8BtJObFE+6s6f+3m6tEH3BBEzABEzABHIE6ikbKRrJtuHogwmYgAmYgAmMmUA9ZbN3nRC1goDFBEzABEzABGoSqKdsvlzDt7rRXlfjvm+ZgAmYgAmYQEagnrL5I3MyARMwARMwgYkSqKdsZLlcB2i5mkpyRaWLvmYCJmACJmACKYF6yua1Qdlo87S8qBvNyiZPxX+bgAmYgAmsRaCesvlM8HHCWj59wQRMwARMwAQKEii6EOcmwNfC/jX65uYfAF2zmIAJmIAJmEBdAkWVjTZMWwy8FTgqnP+gbuh2YAImYAImYAJjWPVZO3OeBfwu/D4HbFaA4GHAAuD+3Npq0etBwVoaCkosXtdxGLgt/K5Kb/jcBEzABEyguQgUtWz+Ezg2KCf5ORq4pk5Su4DzAG1FsCvw9nBMvT0CHA9oWZy8rAD2DL8352/6bxMwARMwgeYhUG+CgFZ3jmuinQp8LyRNCucF4CM1krpfsGgeDG7UFaftpe9J/GiRT8lIOPpgAiZgAibQggTqWTbat2YDQEe5lXLST+e6Xku2AhYmDhYBulZU+oH5wK+At1TxdFJwM3/xYg0pWUzABEzABMpIoJ5lk8Z5JrAjICUQ5efxpAHHucCjYWVpfVh6J/BA7jkXAPoxa9YsWWAWEzABEzCBEhIoqmxOBD4IzA4D9gcAN9VZG02KYk6SZvnVtaIS3aob7gZgrwrKpmhYdmcCJmACJjCFBOp1o8WoSdHsCzwMaL00VfzPxptVjjcHS0jbE/SGCQZFZ5XJiuoL4W4KvDo31lPlkb5sAiZgAiZQRgJFlc1KQD+JlMBvgJeHv6sdNJ35lDBr7V7gUuBu4Ewgzi6TAtNYztuAb4T7Cm+XMBZzO3A98AUrm2qYfd0ETMAEyk+gaDeaFMJGwA+Ba4Elwcqpl8KfAPql8unkD1k/6l7Lyy+B3fMX/bcJmIAJmEBzEiiqbI4MyTsjWBobAj9tziQ71iZgAiZgApNNoKiyUby0a+drwnc3vwAGJjuyfp4JmIAJmEBzEig6ZqOur++ExTc1YH8R8KnmTLJjbQImYAImMNkEilo2fw7skUwS0IC91i3TGmkWEzABEzABE6hJoKhl81juY07NSIvfwdR8gG+agAmYgAmYQD3L5twwRrM0TEvWTDR9qX8I8H/GZwImYAImYAJFCNRTNlqbTKIN064M5zroi36LCZiACZiACRQiUE/ZaFJAFK0CsFP4Q3vUDMYbPpqACZiACZhALQL1lE30+4dhNpq2BOgIa54dBzRyIc74bB9NwARMwASanEBRZfNl4A1h100lWRbOxcC8Jk+/o28CJmACJjAJBIrORutJFI2i9VtA1ywmYAImYAImUJdAUctGEwS+lezUqe9u4uSBug+xAxMwARMwgfYmUFTZvBc4GfirgOtG4Ovtjc6pNwETMAETKEqgiLLpArTU/87AV4oGbHcmYAImYAImEAkUGbMZDuM1W0dPPpqACZiACZjAWAgUsWwUnnbO1MZnWjVgWfKAuAlacsmnJmACJmACJrAmgaLK5m/W9Oa/TMAETMAETKA4gXrKph/Q5IAdgDuBbwPa7tliAiZgAiZgAoUJ1Buz0XI1+wRFczigjzstJmACJmACJjAmAvUsm12B3UOIsmq80vOY8NqxCZiACZiACNSzbNLFNt195jJjAiZgAiYwLgL1LBvtzvlcCFkLcK4X/ta59rXZYFxPtScTMAETMIG2IlBP2eiDTosJmIAJmIAJTIhAvW60CQVuzyZgAiZgAiYgAlY2LgcmYAImYAINJ2Bl03DEfoAJmIAJmICVjcuACZiACZhAwwk0WtkcFhbxvB/4WIXUHAT8OqxKcFTuvradvi/8dG4xARMwARNoUgL1ZqNNJFmayXYecAiwCLgZuAq4Jwn0EeB44CPJNZ1uDHwmrF6gKdbavE1+l+Tc+U8TMAETMIEmINBIy2Y/QBbNg8AAcAlwRI7JQ8AdwEju+qHAtcAzQcHoXFaSxQRMwARMoAkJNFLZbAUsTJjIutG1IlLU70lhe+r5ixcvLhKu3ZiACZiACUwBgUZ2o01Gci4A9GPWrFnqbrOYgAmYgAmUkEAjLZtHgTlJmmcDulZEJuK3SPh2YwImYAImMIkEGqlsNCFgR2BboBc4NgzyF0neNcAbwg6h2iVU57pmMQETMAETaEICjVQ2WiX6lKAk7gUuDVtLnwnE7aT3DTPV3gZ8I9wXRk0MOCvMYJPSkh9ds5iACZiACTQhAa3e3BIyb9680fnz57dEWpwIEzABE5gsAh0dHfq0RJtkNlQaadk0NOIO3ARMwARMoHkIWNk0T145piZgAibQtASsbJo265on4nO22YaOjo4J/RSGxQRMoHkJNPt3Ns1Lvo1ivujhh7loySC+bTYAABXlSURBVOMTSvEJM7eYkH97NgETmFoCtmymlr+fbgImYAJtQcDKpi2y2Yk0ARMwgaklYGUztfz9dBMwARNoCwJWNm2RzU6kCZiACUwtASubqeXvp5uACZhAWxCwsmmLbHYiTcAEpprARD8BaPbp/576PNUl0M83ARNoCwIT/QSg2af/27Jpi2LuRJqACZjA1BKwspla/n66CZiACbQFASubtshmJ9IETMAEppaAlc3U8vfTTcAETKAtCFjZtEU2O5EmYAImMLUErGymln+hp7f7lMlCkBrsqN3zYKrTP9XPV/GaaBwaXERLH7ynPpc+i6Ddp0yWIYvaPQ+mOv1T/XyVwTLEoQzvwnjjYMtmvOTszwRMwARMoDABK5vCqOzQBEzABExgvASsbMZLzv5MwARMwAQKE7CyKYzKDk3ABEzABMZLwMpmvOTszwRMwARMoDABK5vCqOzQBKaOwESn3cq/xQSmkoCnPk8lfT/bBAoS8LTbgqDsrLQEbNmUNmscMRMwARNoHQK2bFonL5sqJaOjowwwytBo+OXOR0dhFP30DzZ/2xH8x/NP09EB3XTQ3dFBVzjqvK+jk2n6dXYxrVPnXfR1dNAhDxYTMIEpJ9BoZXMY8A9AF/At4Au5FPcB/wLMA54GjgEeAtTBfC+wILj/FfDenF//WTICK0aGWTw8yOKhQZ4aHuTp4UGeHR5m5y+dyc+WLWHV6CirRkfCTyqkuOx8zmf50jOLinsAZLav19HJ9M4u1u/sYqOubjbs7GbDrm426uxmo66u7G9d37SrJ/v1ddrYHxNkOzaBggQaqWykYM4DDtFKD8DNwFXAPUnc3g0sAXYAjgXODgpHTh4A9kzc+nSKCTw3PMRjQwM8NrQqOz4+NMCTQwOZgnlqaJBloyNrxbC/o5ON/mBflo+MZNbH9K7u7ChLpLejg56OjhctlWixyFLRv1jtyzb56J778/DvfoeeMCxrKGcJrRoZYfnoCCuy43D2vOUjw9k1HZ8bGebZ4SEeH1rO0uGhinFV5KWUouLJjt2rldAmrzswU56ymPptMa2Vz75gAvUINFLZ7AfcDzwYInEJcERO2ejvM8L9y4F/BNzvUS/XGnhfCuXhwVU8MrgyUyiPDq1CSkVK5oWR4TWePLOzm826e9m6u5+9+9dnVrAOZsUKuqub9Tq76Jj7Si5a8vgafsf6x8pFj/Gy7t6xeqvqfmB0hKXDwywdGcqUkCyx7BesMp0/OLiCJcNDmYLb/cJz+fELz2ThSQmqu2566LaLR3XfTe/oYnpnJ1Ky7sKrin9CN9QFq5KoRscw4TfK6vNwb2R0dHVF0rG6QlGD59aVL2RdLOt1drJeR1dm9epcedXp7tYJ5UkRz41UNlsBC5NIyLrZP/lbp6mbIWApsElwsy1wK/Ac8Cngxpxf/zlOAnpZn86UykoeHtRvVXZ8ZHAVS0aUDatFlsbm3T1s0d3HLtOmsWVPH1t292a/Lbp7M0US3Tbbsbejk1ndncyip2bUVaE9MzzI9vvM48PX/mi1xTS62nJaFroNHx4ZzhRSGpBaTBpDUheelND2n/gQlz+3OFPIUsazunrZuKubrjat5AZHR3h+ZDj7PTc8nDVk1Jh58TeanI8MM++qf+WK559CjYSB0dXjeCnveud7XvxNTntCnSWVJbPAO7uZ2fXST42nzV8s731kA4aVvftqAQKNVDYFHl/ViZrBW4dxHI3n/BB4RVA8qaeTAP1YvHhxen2Nc31joKmj45XZc+ey8CENJTWfaIxkowP2ySq6BwdW8tDgysxqSbu8ZnR2MbenjwPW2yA7zu3tZ+vuvsxqKU1l2NU1IUthvHmo9M/q7uX52+9ibk9/xQIg5a3xqGWJElo+OswydemNDPPM8BBbvvMYzlvy2Br+ZSFt0tWTjSXJSsyPJanii+NM3Rusj1rrZWqBSxEP6oeOGosbzZSBypwUQjY+N7L62h7fv4ATH1uwWrmMDLOyQpdrCkeV/wxZip1dqHwOPPVM1r2prlc1FHroyBS1+uqVR2oYdXZocHj1dbGNo4I6fuFP/pQbbrghm5CiZ6u7dcXoCCtHh7Nz5dWzI0MsGR7k90MD3LtqefZ3DENxO/Dem7jq+adfUkhBOcnKtdQn0Ehl86i2gEiiMBvQtVSiG1k9isuGQcEoj1cFh7eE8ZudgPmpZ+CC8GPWrFlpuVjDWTt8o6AKTy3FZ7IXZihrjasLSEplz0u+lVV0qtC26ennkBkz2bqnn7ndfVkFqhZ26bt8hocn1BV3wswt1igT6/IPsdM4Tj+dbFKl3jlh6915dmjwxQkUmkjx1PBANplC1qTGkWRlLq1SEb/mjhv57nNPZi9JjyrbMN7Vy+pzVbiqYHVUFDqTSlcKas57jufipU++mOzUNtA4mJSDFEamPILSyK4xisbD9rrsoqyizdxkymV0LWvuxcDDieITx+Y6uld3ue6YTdboZoOu1ZM2NEa2fmd3NlYmpaKfuiGVxlQ63rUTp02gK3bp/97Cnv0z0iDrnsuKemJokMdDV/KHz/48W7/vxKxb+cHBlS/6VzdqtFZlDbWzxfoilAonjVQ2mhCwI6DuMCkVTQD4s1wcNGHgOOAm4CjgutAgmQWog1xds9uFcOLYTy6I9vtTLVyNNTw1PJQNWqvlLMWiQXOJunA26OzKxjhU8C88+i+4/7r/ZuOu2l1G7UdyclMsy0W/HXrXq/nglaGVLQUkRaSJDSf+1Qd42999NigFWRMvKQa10jV2oXIRxzJGsrGMlx6z/cdP5YJna4+byVpYPWnjJWUmpSYrY2RgIFMEPR3dqxUdqyd3SCloQke0ODTdfLWC6cwUY2zEnHD0u7hV89mbSGRBzenpy36K9lv+9it86iOnZymIeaQu1jgD86HB1e3jtPtZ3c4batyyTbtL0+xupLJR5/8pwDVh6vOFwN3AmcFCkaL5NvDdMJFAykUKSXJQcDcIWQNK055Xj84GB+1yUAWy3rZb8+DAikS5DBJHVlRBSKHs2NvPzNCq0rRetXCjLPn5TVY0EUYTHPs7O9m8szcbL4jRffSi77PHV74c/yx8VPmR5fKeLbdj2bJlmb9YMmIR0d+qIGtViB1//nJOnYBlUTjCTeKwUh7FMTx1w2lSzaKVz2ep0djd7J4+Nn7tq7JJDem72STJXSfRbKSyUQR/En5pZD+d/CFb9G3J3/H03wD92k40QHrPqmXctWoZd69azoJVy9n/+qu4ccVzWfeIrJMde6exaVd31t8vC6ZWJdF2AJ3gNQioC00dUiMrV6IK0tI4Ahpf0k9d1RK9y9lnAoMDaLz0ld/5Oj94bnGmeLbu6WN2d19mFTYuRuUKudHKplypLVlsNM6iVtBdq5ZnykUK5neDK7POMFUL2/eux8HTZ/LF93+A07/+texDxDINEJcMp6NjAqUioPGnnXqnZT9NpvjoUcdw5He/ycKhVdl7rl4JKZ3tevvZvKu3VJM/GgHSyqYRVKuEKeWiacaa73/7qhe4e9WyrGtMzmVq79o3jYOmbcZufdOzc32jIjnt0h+y8TfOrxKqL5uACZSdgLrOnr7uRl49bcNsbO0Jfcc1sCKrDx4YXJl987NdTz879q6XjeuVPT3jiZ+VTRFqE5h22z9nK7Z70+Ece8anMiUTv2PZrKuHPfpmsFv/9Ey5bNvTv8Y4S5FoFXYzgfjrGV29vQwPDBR+XCkdTpBBKdPkSBUnUKL8V++EvlPTb//RURYNrcq62e4ZWM7dA8vZvKsns4Zk9bTS+I6VTZHiOoZpt5omGpd00SBh/J7ltlUvsHf/DPbqn5Edt+jRsnCTJGOIf6UYadrwRFYAaOS040rxrXhtHTCoGK4vNgeBCea/EtmIcqyZfBrj0U/fZD0wsJLfDqzg5yuW0reygx161msZa8fKZoKvirrG9DW+lnV5dGggW/JEEzw1FXSLrl526+7l3FcfzAv3PeiB/AmytncTaGUC+jh096ynY1rWYL1vYAXR2tmquzebzabZhc06bmtlM47Sq4+9Hh0cyMxfWTHxa+hNurp5Zd90VDD0ZXgsFMsfeMiKZhyc7cUE2pGAZpdu1dOX/bSSuiydBQMrstlsxz+2gCPX35RDZ8zM1udrJj5WNgVzS3PoFw6uymaSaAaZvl3QB2xbdvdlymWr7j5PLS3I0s5MwASKEdAkoT36Z2Tjup84/t3s/E//wNeWPMqNK5bylc22LxZISVxZ2dTICH3BrcUp9/7hd7n8+acyl/quZdfeadlXxVqCPlovNYLxLRMwAROYEAFNFHjyqqs5/99/kn2Hp6nUzSZWNkmOxfGXh4e0WOWqbA+U1bc72KtvRjYn3ktPJMB8agImMOkEdu2bPunPXBcPbHtlIwWj6cha10grImsxSy3foaXFdwkWzClveQcfmMhSHSWadrkuCo3DGAeBdi8D7Z7+cRSZVvPStspGixtKueinlXajgtm9b3q2vP463R54gtMuGzHlstUKcunT0+5loN3TX/oC2vgItpWy0S6UqxXMS5uE6QOqXfqnZUvte+2oxhc4P8EETKA9CbS8stHMsTknHcePX3g6+x5G2fyyrh72618/2yjMGx+1Z8F3qk3ABCaXQEsrG30A9f7H78u25BXWffpnZF/qamVWiwmYgAmYwOQRaGllo2nJH990aw56+S6cf0d+k8/Jg+wnmYAJmEC7E2j5DS72XW99Vi7M70bd7tnu9JuACZjA5BJoactmclH6aSZQYgKeelzizGmPqFnZtEc+O5XtTsBTj9u9BEx5+lu+G23KCTsCJmACJmAC2fbkxmACJmACJmACDSVgy6aheB24CZiACZiACFjZuByYgAmYgAk0nICVTcMR+wEmYAImYAJWNi4DJmACJmACDSdgZdNwxH6ACZiACawDAuFbKW0bPZ7fnG22WQeRGH8Q/s5m/Ozs0wRMwAQmj0CTfytly2byioqfZAImYAJtS8DKpm2z3gk3ARMwgckj0GhlcxiwALgf+FiFZPUBPwj3/xdIOxU/Hq7L/6EV/PqSCZiACZhAkxBopLLRpjHnAYcDuwJvD8cUzbuBJcAOwN8DZ4ebcn8s8ApACuvrgDehScn53ARMwASaiEAjlc1+wTJ5EBgALgGOyLHR398J1y4HDgY6gju5XwX8LoSj8CwmYAImYAJNSEAVe6PkqGCVnBge8BfA/sApyQPvCm4WhWsPBDdnAL8Cvheufxu4GpBCSuUkQD/Jy0OXXfhzjcOmwFNrXGmvP9o9/crtdmfg9Ld3HVDrHZgLzGp0ldjsU58vAPSrJ9qmc596jlr4frunX1nb7gyc/vauA6b8HWhkN5q2x5yTVOCzgfyWmakbKb4NgaeDu3p+k6B9agImYAImUGYCjVQ2NwM7AtsCvWHA/6ocDP19XLimbrfrgFFA1zVBQLPV5F/h/F/Or/80ARMwARMwgYzAG4HfAhqL+WRgcibw5nDeD1wWJgBImWyXcJN7+dPUZ81om4jEcZ2JhNHMfts9/cq7dmfg9DfzG7xu4t7uZWDdUHQoJmACJmACJmACJmACJmACJmACJmACJmACJmACRQlMZHmcos8os7t66T8NuAe4A/gZoPn1rST10h/T+tYwEaUVp8MXYXB0KAd3A9+PUFrkWC/9WwPXA7eG90Djyq0kFwJPAvqGsZLo28qvhXFy1QN7V3Lka7UJaDkbTSjQJAPNfru9wvI47wf+KQSj2W5al61VpEj6/wiYFhL8vjZMv5K+PvDz8NFwqymbImVAMztV0c4M5eBlrfIChCWt6tUB+i5PZV+iZbEeCuetcjgoKJBqykbKVR/IS+kcAGhNykmRRk59npQEJA+ZyPI4STBNe1ok/WrRLQ8p1AoN+vapVaRI+pXWs8IafCtbJeFJOoow+MuwZqHWJJSoFdwqUiT9+rRig5Bgfdf3WKskPqRDDalnaqRJS4T9S7DsVQdsBGxRw/06u9VKymYrYGFCRkvg6FoqqZshYCmwSeqgic/TtCkZldKfJk+LoKqF0ypSJP3qMtDHwv/RKonOpaMIg50A/X4RrDt1O7WKFEm/lsJ6R3g/fgJ8oFUSXzAdRRgVDGpszpp9uZqxpdauIwG9bOpCem280AZHNay+AhzfBmmtlUS98+pK+8Ng2aolvDvwbC1PLXRPq8//M/Bl4A+A7wK7ASMtlMZSJqWVLJt06RvBHsvyOKXMnDFGqkj6FeTrwwe2+rBWq2q3itRLv8ZqVKncEPrp1V+tlSpaadymHgPltSxepXswrKiuj66lfFpBiqRfFv2lIbE3AfqwXIuUtosUYdQuLMadTrXYtJ1BXB5HEwS0H04qJ+cmCMRCl7pp1vMi6d8rTKJolcolzasi6U/dS+m0kqJR2oowULdZ3NZDlay6nlulK7lI+tV1HK3bXcKYTSNXv0/L3GSdaxPKahME/jg3QcDLgI0zVyayPM44H1kqb/XS/1/AE8Bt4acWbitJvfSnaW1FZaP01WOgilXdiZoCf2dYgzDl0uzn9dKvGWgar1JjVO/BG5o9wbn4Xww8HixXWbGy5N4bfnKq/Nemlpq1p/xvtQZXDof/NAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETGB8Bfd/xj+Pzal8m0NwEWmkFgebOCcfeBEzABFqYgJVNC2eukzYhAvoK+zfAvwL3Apcn2zMoYL07Wp5eq+ZGuQ/YDHhTWLpdS/nrQ1pdy4vW5zoqufhCcn46cHPYb+Wz4fr0sICoPkbU1+HHJO59agKlJ2BlU/oscgSnkMDLga8DWtbkOUD7IUXRwo3/DhwZLuwPPBxWaPifsFeIlge6BPjr6KnAUV+0azkhLZe/JzAP0B4lWmZGy+HvEdZ4+2mBsOzEBEpDwMqmNFnhiJSQgNYN09Imku8Brwnn8aDN96KFkW7Gp0VgrwnLgchKya/RF/1XOkrZ6Cer6NfAzkH5aGmRQ8JePAeG7TEq+fc1EyglASubUmaLI1USAtpoKxVtthXXldOq2Vo1eAdgFvAW4Irg+NwwEUBL978nrCychqNz7acU3z8dtbusRGtXfT5YNbJsFP63Aa3OrP14pHQ+B3w6uPfBBJqCQCzsTRFZR9IEJpmA9qvXnieSPwN+nCgBLWIqZXRlWNhS4zpPB7dSSlrKXXJcOOYPGu9RF5lEiqsnnMsiehcwI/ytza60dfOWYZdVWVjneO/4QMeHpiGgJbktJmAClQksALQtxYVhleTzKzhTV5oG8+Oy9XKi3SAvA7T18nVh24u812+GMR8N+Gv8ZVlw8J9hjEhWk0QTB7TZnSwcKRmNFWkvmveF+z6YgAmYgAk0MYFae4I0cbIcdROYGgLuRpsa7n6qCZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZiACZjAmAj8PwixmATf3PWjAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "1011fd87",
   "metadata": {
    "id": "89f02c43"
   },
   "source": [
    "<details>\n",
    "<summary><b>CLICK FOR SOLUTION</b></summary>\n",
    "\n",
    "    dd_cad.predict(emb_nochange, context(emb_nochange))['data']\n",
    "    dd_cad.predict(emb_change, context(emb_change))['data']\n",
    "</details>\n",
    "\n",
    "Hopefully drift isn't detected! We don't expect the context-aware detector to consistently flag drift, as the change in frequency of news topics is permissible given the context provided.\n",
    "\n",
    "### Important note!\n",
    "\n",
    "We've been a bit vauge with the word \"*consistently*\" above. When run many times under *no-drift*, a perfectly calibrated detector should return a uniform distribution of p-values $\\mathcal{U}\\left[0,1\\right]$. Therefore, if you chose a p-value threshold of 0.05, you'd expect a false positive detection 5% of the time. In practice, to quantify how well calibrated a detector is, it must be run many times so that the p-value distribution under no-drift can obtained.\n",
    "\n",
    "Below is the p-value distribution under no-drift for the context aware detector. As you can see, it is relatively well calibrated. For homework, you can generate these plots yourself for the standard and context-aware detectors!\n",
    "\n",
    "![contextMMD_calibration.png](attachment:contextMMD_calibration.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d75ddaaa",
   "metadata": {
    "id": "33d1d5bb"
   },
   "source": [
    "## Homework!\n",
    "\n",
    "To learn more about drift detection, below are a number of topics you can explore in your own time.\n",
    "\n",
    "### Examining detector calibration and test power\n",
    "\n",
    "The performance of a drift detector can be quantified by measuring its:\n",
    "\n",
    "- **calibration**: The distribution of calculated p-values when there is no drift. For a perfectly calibrated detector, we expect the p-values to be uniformly distributed under $U[0,1]$ when there is no drift.\n",
    "\n",
    "- [**test power**](https://en.wikipedia.org/wiki/Power_of_a_test): probability of a \"true positive\" when drift is present, i.e. correctly rejecting the null hypothesis (with the null hypothesis being that there is no drift).\n",
    "\n",
    "To measure the above quantities, we must run many experiments, instantiating detectors and making predictions on many random splits of data. To measure the calibration of the standard and context drift detectors from [section 4](#4.-Accounting-for-context) we set up the same datasets used previously. Since we need to run multiple experiments on different splits, we increase the size of the data splits to 5000 instances.\n",
    "\n",
    "*Note: We are measuring* ***calibration*** *here, since we have set up a scenario where we don't expect drift (given the context!). You can measure* ***test power*** *in a similar way, by measuring the proportion of true positives in a scenario where you do expect drift (see this [Alibi Detect example](https://docs.seldon.io/projects/alibi-detect/en/stable/examples/cd_context_20newsgroup.html#Changing-the-subpopulation-distribution)).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65c153b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_classes = len(classes)\n",
    "n_nochange_c = 5000 // n_classes  # equally subsample each class from 5000 instances\n",
    "\n",
    "n_change_c = [250] * n_classes  # 250 of each class (so 5000 in total), but then mod. below.\n",
    "n_change_c[4], n_change_c[5] = 50, 50  # few stories on Mac/Windows\n",
    "n_change_c[9], n_change_c[10] = 450, 450  # more stories on baseball/hockey\n",
    "\n",
    "(x_nochange, y_nochange), (x_change, y_change), _ = split_data(dataset.data, dataset.target, \n",
    "                                                               n_nochange_c, n_change_c, seed=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5cb65a2",
   "metadata": {},
   "source": [
    "Before we set off our experiments, we compute all necessary embeddings and contexts so we don't have to run the embedding model on every loop. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27d94bf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "emb_nochange = emb_model(x_nochange).cpu().numpy()\n",
    "emb_change = emb_model(x_change).cpu().numpy()\n",
    "c_nochange = context(emb_nochange) \n",
    "c_change = context(emb_change)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b28d1120",
   "metadata": {},
   "source": [
    "We're now ready to run the experiments! On each loop we re-split the `_nochange` and `_change` datasets, init new detectors, and make predictions. If it is taking too long you can reduce `n_runs`, however too few runs will result in a noisy p-value distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38688f76",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "QulWdDXcGnsg",
    "outputId": "d7395981-c1d4-4207-ebdd-2a01877da91b"
   },
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "n_runs = 100\n",
    "n_drift = 1000\n",
    "p_vals_mmd = []\n",
    "p_vals_cad = []\n",
    "for run in tqdm(range(n_runs)):\n",
    "    rng = np.random.default_rng(run)\n",
    "    idx = rng.choice(len(x_nochange), size=n_drift, replace=False)  # assumes nochange and change are same len!\n",
    "\n",
    "    emb_ref = emb_nochange[idx]\n",
    "    emb_predict = emb_change[idx]\n",
    "    c_ref = c_nochange[idx]\n",
    "    c_predict = c_change[idx]\n",
    "    \n",
    "    dd_mmd = MMDDrift(emb_ref, p_val=0.05, n_permutations=100, backend='pytorch')\n",
    "    pred_mmd = dd_mmd.predict(emb_predict)\n",
    "    p_vals_mmd.append(pred_mmd['data']['p_val'])\n",
    "\n",
    "    dd_cad = ContextMMDDrift(emb_ref, c_ref, p_val=.05, n_permutations=100, backend='pytorch')\n",
    "    pred_cad = dd_cad.predict(emb_predict, c_predict)\n",
    "    p_vals_cad.append(pred_cad['data']['p_val'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5722b4d",
   "metadata": {},
   "source": [
    "The computed p-value distributions can be visualised with a histogram plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8580dc1f",
   "metadata": {
    "id": "c6b9abcf"
   },
   "outputs": [],
   "source": [
    "from workshop_utilities import plot_hist\n",
    "\n",
    "p_vals_mmd = np.array(p_vals_mmd)\n",
    "p_vals_cad = np.array(p_vals_cad)\n",
    "\n",
    "p_vals = [p_vals_mmd, p_vals_cad]\n",
    "title = 'p-value distribution for a change in subpopulation prevalence'\n",
    "plot_hist(p_vals, title, ylim=(0,0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cb34287",
   "metadata": {
    "id": "f715dbda"
   },
   "source": [
    "or we can use [Q-Q (Quantile-Quantile) plots](https://en.wikipedia.org/wiki/Q%E2%80%93Q_plot):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa58e10c",
   "metadata": {
    "id": "c5b40539"
   },
   "outputs": [],
   "source": [
    "from workshop_utilities import plot_qq\n",
    "plot_qq(p_vals_mmd, 'Q-Q plot MMD detector')\n",
    "plot_qq(p_vals_cad, 'Q-Q plot Context-Aware MMD detector')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84a3ae31",
   "metadata": {
    "id": "6734e435"
   },
   "source": [
    "A perfectly calibrated detector should have a Q-Q plot which closely follows the diagonal. Only the middle plot in the grid shows the detector's p-values. The other plots correspond to *n_runs* p-values actually sampled from $U[0,1]$ to contextualise how well the central plot follows the diagonal given the limited number of samples.\n",
    "\n",
    "\n",
    "Both plots should show that the context-aware MMD detector is well-calibrated, but the normal MMD isn't! The normal MMD detector should consistently flag drift due to the change in subpopulation prevalence. The context aware drift detector should behave as if there is no drift, since  the change in subpopulation prevalence is permissible given the context provided."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2a30bcf",
   "metadata": {},
   "source": [
    "### Speeding things up!\n",
    "\n",
    "The MMD detector and context-aware MMD detector can take some time to run on large datasets. As discussed [here](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/background.html#dimension-reduction), one solution is to add a dimension reduction step to a detector's `preprocess_fn`. For example, in this Alibi Detect [IMDB text example](https://docs.seldon.io/projects/alibi-detect/en/stable/examples/cd_text_imdb.html#MMD-PyTorch-detector), an Untrained Auto Encoder (UAE) is used to further reduce the 768-dimensional text embedding space to a 32-d space. Applying drift detectors in this lower dimensional space is a good way to speed things up. For an advanced homework task, you can try adding a dimension reduction step to the preprocessors used in this workshop."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b29f1ee",
   "metadata": {
    "id": "a64453f6"
   },
   "source": [
    "### Delving deeper into context-aware drift detection\n",
    "\n",
    "If the distribution of the subpopulations changes (instead of just the frequency), we would now expect the context-aware detector to flag drift. You can check this for yourself by adopting the dataset used in this Alibi Detect [newsgroups example](https://docs.seldon.io/projects/alibi-detect/en/stable/examples/cd_context_20newsgroup.html#Changing-the-subpopulation-distribution).\n",
    "\n",
    "We've only touched the surface of context-aware drift detection today. The detector can be applied to many different tasks by changing the context that is supplied. For more details, see the [ContextMMDDrift docs](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/contextmmddrift.html), the [newsgroup example](https://docs.seldon.io/projects/alibi-detect/en/stable/examples/cd_context_20newsgroup.html), and the [ecg example](https://docs.seldon.io/projects/alibi-detect/en/stable/examples/cd_context_ecg.html).\n",
    "\n",
    "\n",
    "### Interpretability of detections\n",
    "\n",
    "Sometimes you might want to interpret your drift detections in order to understand why drift is occuring. There are a number of things you can explore here:\n",
    "\n",
    "- **Feature-level attribution**: Some of the detectors listed [here](https://docs.seldon.io/projects/alibi-detect/en/stable/overview/algorithms.html) return a p-value for each feature, meaning you can attribute a drift detection to certain features. \n",
    "- **Instance-level attribution**: Some detectors, such as [ClassifierDrift](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/classifierdrift.html) return instance-level predictions, allowing you to attribute drift to certain instances.\n",
    "- **Detectors designed for interpretability**: The [SpotTheDiff](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods/spotthediffdrift.html) detector is specifically designed to give detections that are interpretable at the feature level when they occur.\n",
    "\n",
    "If you follow the \"Examples\" links at the bottom of each drift detector docs page, you can explore how some of the above capability can be used. We are also working on detection methods which take advantage of our ML interpretability library [Alibi Explain](https://docs.seldon.io/projects/alibi/en/latest/), so watch this space!\n",
    "\n",
    "### Online drift detection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb6ecfd0",
   "metadata": {},
   "source": [
    "A final thing we haven't covered today is \"online\" drift detection. The describes the scenario when test instances arrive sequentially instead of in batches. In such a scenario it is important for the detector to be well calibrated (we desire a certain Expected Run Time in the absense of drift), yet it flag drift quickly when it is present. Alibi Detect offers a number of state-of-the-art [online detectors](https://docs.seldon.io/projects/alibi-detect/en/stable/cd/methods.html#online) that are designed to meet these requirements. For more info, see [this paper](https://arxiv.org/abs/2108.00883) and a related YouTube video: [Protecting Your Machine Learning Against Drift: An Introduction](https://www.youtube.com/watch?v=tL5sEaQha5o).\n",
    "\n",
    "The syntax for the online detectors is consistent with the offline detectors used in this workshop:\n",
    "\n",
    "```python\n",
    "# Init detector\n",
    "from alibi_detect.cd import MMDDriftOnline\n",
    "detector = MMDDriftOnline(x_ref, ERT, W, backend='pytorch', preprocess_fn=encoder_fn)\n",
    "\n",
    "# Make a prediction on a test instance\n",
    "detector.predict(x)\n",
    "```\n",
    "\n",
    "For a good example, check out the [medical imaging example](https://docs.seldon.io/projects/alibi-detect/en/stable/examples/cd_online_camelyon.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "233a560d",
   "metadata": {},
   "source": [
    "## Follow-up questions\n",
    "\n",
    "Any follow-up questions or problems? Please feel free to get in touch at ashley.scillitoe@seldon.io "
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "name": "Copy of intro_to_drift_detection_master.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}