followup_ef_final.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "followup_ef_final.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tc0c_buJhhgt",
        "colab_type": "text"
      },
      "source": [
        "This notebook allows you to train and run the stochastic follow-up EF model presented here (https://ecs193-final-project.herokuapp.com/). In order to run this code, you will need have the cleaned data downloaded. This data is located here: https://drive.google.com/file/d/1ecNmerBZLMQlaWnoyNkE4bZrHPV8gkoL/view?usp=sharing"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-ZH7SP2IhdE5",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "ab2ad55a-eb0d-4a3f-aade-1dd57a2f330c"
      },
      "source": [
        "try:\n",
        "  # Use the %tensorflow_version magic if in colab.\n",
        "  %tensorflow_version 1.x\n",
        "except Exception:\n",
        "  pass"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "TensorFlow 1.x selected.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8282cRWykho8",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "145282f8-375e-42b1-f965-1cd8efa36d2b"
      },
      "source": [
        "import pandas as pd\n",
        "import sklearn\n",
        "import keras\n",
        "import tensorflow as tf\n",
        "import operator\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import warnings\n",
        "import math\n",
        "import itertools\n",
        "from keras import backend as K\n",
        "from time import time\n",
        "from keras.utils import plot_model\n",
        "from sklearn.model_selection import train_test_split,KFold\n",
        "from sklearn.utils import shuffle\n",
        "from sklearn.metrics import precision_recall_curve,average_precision_score,roc_auc_score,roc_curve,recall_score,precision_score,f1_score\n",
        "from datetime import datetime, date\n",
        "from sklearn.metrics import confusion_matrix\n",
        "from keras.layers.merge import concatenate\n",
        "from keras.layers import Input, Dense, Conv1D, MaxPooling1D, Softmax, Add, Flatten, Activation, Embedding, LSTM, Dropout, Dropout\n",
        "from keras.models import Model\n",
        "from keras.optimizers import Adam\n",
        "warnings.filterwarnings(\"ignore\")\n",
        "%matplotlib inline"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Using TensorFlow backend.\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jjxF4QQQl7Kq",
        "colab_type": "text"
      },
      "source": [
        "Here, we clean up the the .csv files downloaded. In particular, a lot of column values have no values, so we impute them with np.nan so they can be zero filled later."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nDJ6HQyPlQb4",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "df2 = pd.read_csv(\"combined_echo_data_new.csv\")\n",
        "df2 = df2.drop(\"Unnamed: 0\",axis=1)\n",
        "df2 = df2.dropna(axis=0,subset = [\"Pt. ID\"])\n",
        "exclude_mrns = []\n",
        "\n",
        "for ids in np.array(df2['Pt. ID'].unique()):\n",
        "  try:\n",
        "    ids1 = int(ids) \n",
        "  except ValueError:\n",
        "    exclude_mrns.append(ids)\n",
        "\n",
        "df2 = df2[~df2['Pt. ID'].isin(exclude_mrns)]\n",
        "df2['Pt. ID'] = df2['Pt. ID'].astype(int)\n",
        "\n",
        "vol, ef = [], []  # Created new cols\n",
        "\n",
        "for index, record in df2.iterrows():\n",
        "  volume_A4C = record[\"LA Vol A/L A4C i BSA\"]\n",
        "  volume_A2C = record[\"LA Vol A/L A2C i BSA\"]\n",
        "  volume_BP =  record[\"LA Vol A/L BP i BSA\"]\n",
        "  ef_A4C = record[\"LV EF, MOD A4C\"]\n",
        "  ef_A2C = record[\"LV EF, MOD A2C\"]\n",
        "  ef_BP = record[\"LV EF, MOD BP\"]\n",
        "  # First pick is the A4C channel\n",
        "  if not np.isnan(volume_A4C):\n",
        "    vol.append(volume_A4C)\n",
        "  # Next pick is from A2C channel\n",
        "  elif not np.isnan(volume_A2C):\n",
        "    vol.append(volume_A2C)\n",
        "  # Worst case, choose from BP channel\n",
        "  elif not np.isnan(volume_BP):\n",
        "    vol.append(volume_BP)\n",
        "  # If all channels are empty; set it empty\n",
        "  else:\n",
        "    vol.append(float(\"nan\"))\n",
        "  \n",
        "  # Do the same for EF\n",
        "  if not np.isnan(ef_A4C):\n",
        "    ef.append(ef_A4C) \n",
        "  elif not np.isnan(ef_A2C):\n",
        "    ef.append(ef_A2C) \n",
        "  elif not np.isnan(ef_BP):\n",
        "    ef.append(ef_BP) \n",
        "  else:\n",
        "    ef.append(float(\"nan\"))\n",
        "\n",
        "df2[\"LA Vol A/L Volume\"] = vol # New column for Volume added\n",
        "df2[\"Ejection_Fraction\"] = ef  # New column for EF added\n",
        "df2 = df2.drop([\"LA Vol A/L A4C i BSA\",\"LA Vol A/L A2C i BSA\",\n",
        "                \"LA Vol A/L BP i BSA\",\"LV EF, MOD A4C\",\n",
        "                \"LV EF, MOD A2C\",\"LV EF, MOD BP\"], axis=1) # Drop these obsolete cols\n",
        "df2 = df2.dropna(axis = 0, how=\"all\",subset=[\"Ejection_Fraction\"])\n",
        "ind = df2[df2['Ejection_Fraction'] <= 0].index.tolist() + df2[df2['IVS d, 2D'] > 5].index.tolist() \n",
        "ind = ind + df2[df2['LV d, 2D'] > 10].index.tolist() + df2[df2['LV s, 2D'] > 10].index.tolist() \n",
        "ind += df2[df2['RVSP (TR)'] > 125].index.tolist() + df2[df2['LA Vol A/L Volume'] > 200].index.tolist()\n",
        "df_echo = df2.drop(ind,axis=0)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "USkLHGXpoXub",
        "colab_type": "text"
      },
      "source": [
        "We also need to zero fill any np.nan values and normalize."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qpNIqwBPoUgs",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "zero_filled_values = np.nan_to_num(df_echo[df_echo.columns.to_list()[2:]].values)\n",
        "prior_values = df_echo[df_echo.columns.to_list()[0:2]].values\n",
        "all_values = np.append(prior_values, zero_filled_values, axis=1)\n",
        "\n",
        "df_echo_norm = pd.DataFrame(all_values,columns=df_echo.columns.to_list())\n",
        "\n",
        "# normalizing the values\n",
        "df_echo_norm['IVS d, 2D'] = df_echo_norm['IVS d, 2D'] / 5.0\n",
        "df_echo_norm['LV d, 2D'] = df_echo_norm['LV d, 2D'] / 10.0\n",
        "df_echo_norm['LV s, 2D'] = df_echo_norm['LV s, 2D'] / 10.0\n",
        "df_echo_norm['RVSP (TR)'] = df_echo_norm['RVSP (TR)'] / 125.0\n",
        "df_echo_norm['LA Vol A/L Volume'] = df_echo_norm['LA Vol A/L Volume'] / 200.0\n",
        "df_echo_norm['Ejection_Fraction'] = df_echo_norm['Ejection_Fraction'] / 100.0"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xdMZYxiHpdU8",
        "colab_type": "text"
      },
      "source": [
        "Now, we collect the follow-up EF samples using the data from the EF sample. We check to see if there is a follow-up EF sample 90-365 days following the initial EF date. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hSu2jaOuok0p",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "933f342b-67ac-4e95-e7e1-88df55b5c281"
      },
      "source": [
        "from tqdm import tqdm\n",
        "\n",
        "date_format = \"%Y-%m-%d %H:%M:%S\"\n",
        "\n",
        "x_values = None\n",
        "y_values = None\n",
        "\n",
        "groups = df_echo_norm.groupby(by=\"Pt. ID\")\n",
        "\n",
        "for grp_name, record in tqdm(groups):\n",
        "  sorted_record = record.sort_values(['Study Date'])\n",
        "  sorted_record['Study Date'] = pd.to_datetime(sorted_record['Study Date'])\n",
        "  for i in range(len(sorted_record)):\n",
        "    for j in range(i+1,len(sorted_record)):\n",
        "      mrn = sorted_record.iloc[i]['Pt. ID']\n",
        "      initial_echo_date = sorted_record.iloc[i]['Study Date']\n",
        "      end_echo_date = sorted_record.iloc[j]['Study Date']\n",
        "\n",
        "      delta = (end_echo_date - initial_echo_date).days\n",
        "\n",
        "      if delta > 90 and delta < 365:\n",
        "        input_indices = [0, 2, 3, 4, 5, 6, 7]\n",
        "        data_input = [sorted_record.iloc[i].values[input_indices]]\n",
        "        data_output = [[sorted_record.iloc[j]['Ejection_Fraction']]]\n",
        "        if x_values is None:\n",
        "          x_values =  data_input\n",
        "          y_values = data_output\n",
        "        else:\n",
        "          x_values = np.append(x_values, data_input, axis=0)\n",
        "          y_values = np.append(y_values, data_output)"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "100%|██████████| 35627/35627 [01:31<00:00, 389.13it/s]\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1VgTvCn9qU-7",
        "colab_type": "text"
      },
      "source": [
        "Now, we will undersample the data so there is an equal balance of data of data points where EF < 0.35 and EF >= 0.35."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ji3Ue0bNqUb5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def undersample(X, y, threshold):\n",
        "  assert (len(X) == len(y)), \"lengths of X and y are not equal\"\n",
        "  X_low = None\n",
        "  X_high = None\n",
        "  y_low = None\n",
        "  y_high = None\n",
        "\n",
        "  for i in range(len(y)):\n",
        "    if y[i] < threshold:\n",
        "      if X_low is None:\n",
        "        X_low = [X[i]]\n",
        "        y_low = y[i]\n",
        "      else:\n",
        "        X_low = np.append(X_low, [X[i]], axis=0)\n",
        "        y_low = np.append(y_low, y[i])\n",
        "    else:\n",
        "      if X_high is None:\n",
        "        X_high = [X[i]]\n",
        "        y_high = y[i]\n",
        "      else:\n",
        "        X_high = np.append(X_high, [X[i]],axis=0)\n",
        "        y_high = np.append(y_high, y[i])\n",
        "\n",
        "  if len(y_low) < len(y_high):\n",
        "    chosen = np.random.choice(len(y_high), len(y_low), replace=False)\n",
        "    X_high = X_high[chosen]\n",
        "    y_high = y_high[chosen]\n",
        "  else:\n",
        "    chosen = np.random.choice(len(y_low), len(y_high), replace=False)\n",
        "    X_low = X_low[chosen]\n",
        "    y_low = y_low[chosen]\n",
        "\n",
        "  X_new = np.append(X_low, X_high, axis=0)\n",
        "  y_new = np.append(y_low, y_high)\n",
        "  X_new, y_new = shuffle(X_new, y_new, random_state=0)\n",
        "\n",
        "  return X_new, y_new"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mYmiMbJQr1Dw",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x_values, y_values = undersample(x_values, y_values, 0.35)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Iw6rQ-tZsAda",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "all_training_data = np.append(x_values, np.expand_dims(y_values, axis=1), axis=1)\n",
        "df_final_training_data = pd.DataFrame(all_training_data,\n",
        "                                      columns=[\"Pt. ID\", \"IVS d, 2D\",\t\"LV d, 2D\",\t\n",
        "                                                \"LV s, 2D\",\t\"RVSP (TR)\",\t\"LA Vol A/L Volume\", \n",
        "                                                \"Ejection_Fraction\",\t\"Ejection_Fraction_Followup\"])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sVeWWxOMsZ3b",
        "colab_type": "text"
      },
      "source": [
        "We need to make sure we split by individual patients so same patients are not in both train and test."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SWZdlwAosCxf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "ids = df_final_training_data[\"Pt. ID\"].unique() # Get all the different patients"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mbLTxPRusppr",
        "colab_type": "text"
      },
      "source": [
        "We'll do uncertainty estimation here. First, we need to instantiate a stochastic model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OqIRinW5smeL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class StochasticModel(keras.Model):\n",
        "  def __init__(self):\n",
        "    super(StochasticModel, self).__init__()\n",
        "\n",
        "    self.layer_1 = Dense(100, activation='relu')\n",
        "    self.dropout_1 = Dropout(0.5)\n",
        "\n",
        "    self.layer_2 = Dense(50, activation='relu')\n",
        "    self.dropout_2 = Dropout(0.5)\n",
        "    \n",
        "    self.layer_3 = Dense(10, activation='relu')\n",
        "    self.dropout_3 = Dropout(0.5)\n",
        "\n",
        "    self.output_layer = Dense(1, activation='sigmoid')\n",
        "  \n",
        "  def call(self, x):\n",
        "    x = self.layer_1(x)\n",
        "    x = self.dropout_1(x, training=True)\n",
        "\n",
        "    x = self.layer_2(x)\n",
        "    x = self.dropout_2(x, training=True)\n",
        "\n",
        "    x = self.layer_3(x)\n",
        "    x = self.dropout_3(x, training=True)\n",
        "\n",
        "    output = self.output_layer(x)\n",
        "\n",
        "    return output"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lk66-vlhsunT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x_data = df_final_training_data[df_final_training_data.columns.to_list()[:-1]].values\n",
        "y_data = df_final_training_data[df_final_training_data.columns.to_list()[-1]].values"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Q4ThgmR0tv66",
        "colab_type": "text"
      },
      "source": [
        "We will split by patient IDs here and create the train and test sets ready for model input."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-CA24QTyteTk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "train_ids, test_ids = train_test_split(ids, test_size = 0.10)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "m5dZ8gBHtfxQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "df_train = df_final_training_data[df_final_training_data[\"Pt. ID\"].isin(train_ids)]\n",
        "df_test = df_final_training_data[df_final_training_data[\"Pt. ID\"].isin(test_ids)]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "L_IYQ_twtnp-",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x_train = df_train[df_final_training_data.columns.to_list()[1:-1]].values\n",
        "y_train = df_train[df_final_training_data.columns.to_list()[-1]].values"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ICy7SYMWtrpN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x_test = df_test[df_final_training_data.columns.to_list()[1:-1]].values\n",
        "y_test = df_test[df_final_training_data.columns.to_list()[-1]].values"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TIofJZgWtvdY",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 88
        },
        "outputId": "2363bc28-a885-4bbe-d3b7-53d6cc8b6a46"
      },
      "source": [
        "batch_size = 1000\n",
        "epochs = 50\n",
        "adam = Adam(lr = 0.002, beta_1 = 0.9, beta_2 = 0.999)"
      ],
      "execution_count": 22,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From /tensorflow-1.15.2/python3.6/tensorflow_core/python/ops/resource_variable_ops.py:1630: calling BaseResourceVariable.__init__ (from tensorflow.python.ops.resource_variable_ops) with constraint is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "If using Keras pass *_constraint arguments to layers.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r9SCjr1Kt3qu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def mse_loss(y_actual, y_pred):\n",
        "\treturn K.mean(K.square(y_pred - y_actual))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xqxH0RJZt6tN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "stochastic_model = StochasticModel()\n",
        "stochastic_model.compile(loss=mse_loss, optimizer=adam, metrics=[mse_loss])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ETr6mwCCt8wU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "stochastic_model.fit(x_train, y_train, batch_size = batch_size, validation_data=(x_test, y_test), epochs = 100)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y0Dz6oq4usTE",
        "colab_type": "text"
      },
      "source": [
        "Some helper functions."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6rhkEhzIurnE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def binarize_output(output,limit = 0.35):\n",
        "  processed_output = []\n",
        "  for elem in output:\n",
        "    if elem < limit:\n",
        "      processed_output.append(0)\n",
        "    else:\n",
        "      processed_output.append(1)\n",
        "  \n",
        "  return processed_output"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4D7fnb99uy52",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def basic_accuracy_metric(y_pred, y_truth):\n",
        "  count = 0\n",
        "  for i in range(len(y_pred)):\n",
        "    if y_pred[i] == y_truth[i]:\n",
        "      count = count + 1\n",
        "  return count * 100.00/ len(y_pred)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "P0q3nn83uH9h",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "y_test_binarized = binarize_output(y_test)\n",
        "stochastic_model_preds = binarize_output(stochastic_model.predict(x_test))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Fc8TrYwWvGH0",
        "colab_type": "text"
      },
      "source": [
        "Here, the run through every test point in the dataset and perform inference on in N (20 in this setting times). We use the mean as the prediction of the model and the variance as a measure of uncertainty."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AptFmOkju991",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "c9bf2e28-8f09-4d21-fade-40c014fb1284"
      },
      "source": [
        "from tqdm import tqdm\n",
        "\n",
        "correct_pred_uncertainties_negative = []\n",
        "correct_pred_uncertainties_positive = []\n",
        "\n",
        "incorrect_pred_uncertainties_negative = []\n",
        "incorrect_pred_uncertainties_positive = []\n",
        "\n",
        "correct_pred_mean_negative = []\n",
        "correct_pred_mean_positive = []\n",
        "\n",
        "incorrect_pred_mean_negative = []\n",
        "incorrect_pred_mean_positive = []\n",
        "\n",
        "all_predictions = []\n",
        "all_uncertainties = []\n",
        "\n",
        "num_samples_per_example = 100 # hyperparameter that can be changed\n",
        "\n",
        "for i in tqdm(range(len(x_test))): # iterating through every example\n",
        "  ground_truth = y_test_binarized[i]\n",
        "  sample_predictions = []\n",
        "  for j in range(num_samples_per_example):\n",
        "    sample_prediction = stochastic_model.predict(np.expand_dims(x_test[i], axis=0))[0][0]\n",
        "    sample_predictions.append(sample_prediction)\n",
        "  \n",
        "  sample_predictions = np.array(sample_predictions)\n",
        "\n",
        "  raw_prediction_mean = np.mean(sample_predictions)\n",
        "  raw_prediction_variance = np.var(sample_predictions) # measure of uncertainty, un-binarized\n",
        "\n",
        "  all_predictions.append(raw_prediction_mean)\n",
        "  all_uncertainties.append(raw_prediction_variance)\n",
        "\n",
        "  if raw_prediction_mean < 0.35:\n",
        "    prediction_mean = 0\n",
        "  else:\n",
        "    prediction_mean = 1\n",
        "\n",
        "  if prediction_mean == ground_truth:\n",
        "    if prediction_mean == 0:\n",
        "      correct_pred_uncertainties_negative.append(raw_prediction_variance)\n",
        "      correct_pred_mean_negative.append((raw_prediction_mean, y_test[i]))\n",
        "    else:\n",
        "      correct_pred_uncertainties_positive.append(raw_prediction_variance)\n",
        "      correct_pred_mean_positive.append((raw_prediction_mean, y_test[i]))\n",
        "  else:\n",
        "    if prediction_mean == 0:\n",
        "      incorrect_pred_uncertainties_negative.append(raw_prediction_variance)\n",
        "      incorrect_pred_mean_negative.append((raw_prediction_mean, y_test[i]))\n",
        "    else:\n",
        "      incorrect_pred_uncertainties_positive.append(raw_prediction_variance)\n",
        "      incorrect_pred_mean_positive.append((raw_prediction_mean, y_test[i]))"
      ],
      "execution_count": 38,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "100%|██████████| 314/314 [00:21<00:00, 14.79it/s]\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vUXZwAZvvWLc",
        "colab_type": "text"
      },
      "source": [
        "Next, we will produce the uncertainty estimate plots."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gCbUKWRKvcmw",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "correct_negative_errors = []\n",
        "correct_positive_errors = []\n",
        "\n",
        "incorrect_negative_errors = []\n",
        "incorrect_positive_errors = []\n",
        "\n",
        "for (prediction, ground_truth) in correct_pred_mean_negative:\n",
        "  correct_negative_errors.append(abs(prediction - ground_truth))\n",
        "\n",
        "for (prediction, ground_truth) in correct_pred_mean_positive:\n",
        "  correct_positive_errors.append(abs(prediction - ground_truth))\n",
        "\n",
        "for (prediction, ground_truth) in incorrect_pred_mean_negative:\n",
        "  incorrect_negative_errors.append(abs(prediction - ground_truth))\n",
        "\n",
        "for (prediction, ground_truth) in incorrect_pred_mean_positive:\n",
        "  incorrect_positive_errors.append(abs(prediction - ground_truth))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AXfvA0Xfv0Ls",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "outputId": "ce4a795d-0a51-472d-bcf0-ed18b85a0e09"
      },
      "source": [
        "plt.title(\"Uncertainty in negative class predictions\")\n",
        "plt.xlabel(\"Error\")\n",
        "plt.ylabel(\"Uncertainty\")\n",
        "plt.scatter(correct_negative_errors, correct_pred_uncertainties_negative, marker=\"s\", label='Correct', alpha=0.5)\n",
        "plt.scatter(incorrect_negative_errors, incorrect_pred_uncertainties_negative, marker=\"o\", label='Incorrect', alpha=0.5)\n",
        "plt.legend(loc='upper right');\n",
        "plt.show()"
      ],
      "execution_count": 40,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZwcdZ3/8dc7M8nMJOROBHORSAJCEBCS4ApyyRFYMe4umKDwQxaJIKCs/FzB64f8Fpe4ysoCq+QHGLlMIqKOLpcIWRblmCQEyMERSIBJCITc10wyw+f3R31n6Ol0z3TP9FHd/Xk+HvOY6m9VdX+qZro+9T2qSmaGc845l6lexQ7AOedcafHE4ZxzLiueOJxzzmXFE4dzzrmseOJwzjmXFU8czjnnsuKJw+WcpG9Luq1An/WgpPML8Vn5JGm7pI8U8PMWSPpyoT4vnySdIKkx4fUySSd0430+JenlnAZXpjxxxJQkkzQ+qewaSXcXIZasDjJm9kMzy2j5nm6TmZ1uZr/s7vrFkGp/mtk+ZvZ6sWIqJ2Y20cwWdLVc8nfMzP7HzA7Ka3BlwhOHS0sR/x9xBSWpqtgxuM75QaFEtVXPJV0p6V1Jb0u6IGF+naSfSHpD0hZJT0qqC/M+IemvkjZLej6xWh/Ohq+T9BdgJ3AX8Cng5tCccnNY7kZJb0naKmmRpE8lvEd7LULS2HBmd76kNyW9J+k7Yd5U4NvA9PDez0s6W9KipG39hqTfp9kP7Wfvkr4UtvPHkjZJWiXp9E724WpJ/1vSC2EfzZNUmzD/M5KWhP30V0mHJcw7UtJzkrZJ+nVY91/CvMGS/ihpfYjjj5JGhXnXpdmfJmm8pKMlrUs8eEr6O0kvhOlekq6S9JqkDZLmSxrSyTZOC9uwNawzNcUyB0h6LLzfe5LukTQoYf63JK0J2/qypE+H8imSFob3fkfSDWliaPtf/XZ4/9WSvpgwf46kn0l6QNIO4ERJIyT9JuzDVZK+lrB8XVhnk6TlwOQUf9eTw3RV+NzXQvyLJI2W9ERY/Pnwd5iuvZu8Dg7/X5sVNX99NinmWyT9V3jfZyQdEOZJ0r8r+l5ulfSipEPT/Y1Kkpn5Twx/AAPGJ5VdA9wdpk8AWoBrgd7AGUQH+sFh/i3AAmAkUAV8EqgJrzeE5XsBp4TXw8N6C4A3gYlAdXjvBcCXk2I5FxgalrkSWAfUpohzbNiW/wfUAYcDzcDBycuG1zXAxrb5oew54B/S7Kf22IAvAXuAi8I2XwKsBZRm3dXAs8AIYAiwArg4zPs48C5wdHiv88PyNUAf4A3g62H//D2wG/iXsO5Q4B+AvkB/4NfA71LFnOrvDbwGnJIw79fAVWH668DTwKgQy63Ar9Js3xRgS/gb9wp/+4+m2G/jwzI1wHDgCeCnYd5BwFvAiIS/5wFh+ingvDC9D/CJNHGcQPS/ekP4jOOBHcBBYf6cEOcxIc6+wCLg+2FffwR4HTgtLH898D/hbzYaWAo0Jv1dTw7T3wReDNshov+/oam+YyHOxjDdG1hJdGLTBzgJ2JYU84awj6uBe4C5Yd5pIf5B4TMPBj5c7GNKTo9PxQ7Af9L8YTJLHLuA6oT57wKfCF++XcDhKd73W8BdSWUPA+eH6QXAtUnz2w8yncS7qe3zSJ04RiUs+ywwI3nZhPk/A64L0xPDe9ek+dz22IgSx8qEeX3DZ++XZt3VwLkJr38E/Dwhhv+btPzLRAe944A1JCQk4ElC4kjxOUcAmzrbn3RMHP8C3BGm+xMdZPcPr1cAn05Y78NEybI6xefeCvx7V/stxbzPAc+F6fHh/+pkoHfSck8APwCGdfG/cQJR4uiXUDYf+F6YngPcmTDvaODNpPe4GvhFmH4dmJowbybpE8fLwLRMvmN0TByfIjoZ6pUw/1fANQkx35Yw7wzgpTB9EvAK4bvY2b4p1R9vqoqvVqKznkS9iQ4SbTaYWUvC651EZ37DgFqiM9dk+wNnh+r3ZkmbgWOJDkBt3uoquNDEsyI08WwGBobPTWddijjT+SXwBUkCzgPmm1lzVzElf46Z7QyTnX1Wurj2B65M2k+jiWonI4A1Fo4SQfs+k9RX0q2Kmgm3Eh1gBynztvt7gb+XVENUm1lsZm8kxPXbhJhWEP2v7JvifUaT+n+gA0n7SpobmqO2AncT/pZmthK4gijBvxuWGxFWvRA4EHhJUoOkz3TyMZvMbEfC6zeI9mObxP+5/YERSfv+2wnbOCJp+TdIL6N9kMII4C0zez/pc0YmvE75v2NmjwE3E9X635U0W9KAbsQQW5444utNorP1ROPo/EvS5j2gCTggxby3iGocgxJ++pnZ9QnLJN8yucNrRf0Z/wx8nqhpbBBRU4MyiC3ZXrdnNrOniZp+PgV8gaifpdDeIqr1JO6nvmb2K+BtYGRIbG1GJ0xfSdQ0crSZDSCqocAH+6fTW1Kb2XKiv/PpRNt/b1JcpyfFVWtma9JsQ6r/gWQ/DDF9LMR7bkKsmNm9ZnYs0QHdgFmh/FUzOwf4UCi7T1K/NJ8xOGneGKJmxPaPSYp7VdI29jezM8L8t+m4v8d0sm2Z7oNka4HR6jg4ZAxRTbNLZvYfZnYUcAhRcv1mN2KILU8c8TUP+K6kUaFD9GTgTOC+rlYMZ0l3ADeETsYqSX8TzmDvBs6UdFoorw2dgqM6ect3iNqZ2/QnanpYD1RL+j7Q3TOqd4Cx2nv01p1EZ217zOzJbr53T/w/4GJFndWS1E/S30rqT9S23wpcJqla0jSitu42/YmaCjcr6rj+P0nvnbw/U7mXqD/jOKI+jjY/B66TtD+ApOHh81O5HbhA0qfD/9BISR9NsVx/YDuwRdJIEg5ykg6SdFL432kK2/V+mHeupOHh/21zWOV90vuBpD7hxOMzSduV6Flgm6JO+brwf3qopLZO8PnA1YoGIYwCLu/kM28D/q+kCeHveJikoWFeZ3+HZ4hqEf8sqbeiASRnAnM7+SwAJE0O/ze9iZoZm+h8v5QcTxzxdS3wV6K2801E7e9fNLOlGa7/v4k6BRuIOptnEbW3vgVMI6r6ryc6I/smnf8v3AicpWgUy38Q9Yk8RNSO+wbRF6PL5q002g4eGyQtTii/CziUKNEVnJktJOpkv5lo/68k6kPBzHYTNSFdSHTAPBf4I1GnP8BPiQYCvEfUkf1Q0tsn789UfkXUn/KYmb2XtG498IikbeH9j06zDc8CFwD/TlQj/G+iWkOyHwBHhmX+C7g/YV4NUWf0e0RNMx8i6m8AmAosk7Q9xDXDzHal2Z51RPtxLVFH8sVm9lKauFuJEssRwKrw2bcRNYe2xftGmPcInddIbyBKNI8AW4mSaV2Ydw3wy9Ac9vmkGHYTJYrTw+f/J/C/0sWcZADRicemEOcG4N8yWK9kqGMzrXPxoGjo8LvAkWb2arHj6YqkZ4g61n9R7FjiJpyt321mndVqXQnxGoeLq0uAhrgmDUnHS9ovNFWdDxzG3jUL58pSdbEDcC6ZpNVEnbOfK3IonTmIqAmkH9Hw0LPM7O3ihuRcYXhTlXPOuax4U5VzzrmsVERT1bBhw2zs2LHFDsM550rGokWL3jOz4anmVUTiGDt2LAsXLix2GM45VzIkpb3Y2JuqnHPOZcUTh3POuax44nDOOZeViujjcM6Vtz179tDY2EhTU1OxQyk5tbW1jBo1it69k2/GnZ4nDudcyWtsbKR///6MHTuWjjctdp0xMzZs2EBjYyPjxo3LeD1PHHl0/YMrWL9t78dIDO9fw1WnH1yEiJwrT01NTZ40ukESQ4cOZf369Vmt54kjj9Zva2bU4L57lTdu2pliaedcT3jS6J7u7DfvHHfOOZeVvCYOSVMlvSxppaSrUsyvkTQvzH9G0tiEeVeH8pclnZZQ/nVJSyUtk3RFPuN3zrlsrFu3jhkzZnDAAQdw1FFHccYZZ/DKK68U5LPnzJnD2rVru14wB/KWOMLzlW8hehDKIcA5kg5JWuxComcRjyd62MyssO4hwAxgItHDYv6z7SlgRA/XmQIcDnxG0vh8bYNzzmXKzPi7v/s7TjjhBF577TUWLVrEv/7rv/LOO+90uW5LS0unrzNRyMSRzz6OKcBKM3sdQNJcoifPLU9YZhrRU7ggeiTqzeE5ztOAuWbWDKyStDK83yjgGTPbGd7zv4mexPajPG6Hc66M5GvQyuOPP07v3r25+OKL28sOP/xwzIxvfvObPPjgg0jiu9/9LtOnT2fBggV873vfY/Dgwbz00kvMnj27w+sVK1Zw1VVXsWDBApqbm7n00kv5yle+AsCsWbO4++676dWrF6effjqTJk1i4cKFfPGLX6Suro6nnniMOpqhdTdU9YG6QdC7Ll3oWctn4hhJx8eJNrL3Iy7blzGzFklbgKGh/OmkdUcCS4metzyU6NnHZwApb0IlaSYwE2DMmM6eZZ8/w/vXpOwIH96/pgjROOcgf4NWli5dylFHHbVX+f3338+SJUt4/vnnee+995g8eTLHHXccAIsXL2bp0qWMGzeOBQsWdHg9e/ZsBg4cSENDA83NzRxzzDGceuqpvPTSS/z+97/nmWeeoW/fvmzcuJEhQ4Zw88038+Mf/5hJh0+E7e9Cr6ooaVhr9HqfD+UseZTUqCozWyFpFtHzg3cAS4DWNMvOBmYDTJo0qSgPHfEht865J598knPOOYeqqir23Xdfjj/+eBoaGhgwYABTpkzpcP1E4utHHnmEF154gfvuuw+ALVu28Oqrr/Loo49ywQUX0LdvlPyGDBnS8QN3bY6Shqqi16qKOiV2bc5Z4shn5/gaYHTC61GhLOUykqqJHka/obN1zex2MzvKzI4jehh8YXqenHOuExMnTmTRokVZrdOvX7+0r82Mm266iSVLlrBkyRJWrVrFqaee2vWbtu7+IGm0UVVUniP5TBwNwARJ4yT1Iersrk9aph44P0yfBTxm0SMJ64EZYdTVOGAC8CyApA+F32OI+jfuzeM2OOdcRk466SSam5uZPXt2e9kLL7zAoEGDmDdvHq2traxfv54nnniCKVOmdPl+p512Gj/72c/Ys2cPAK+88go7duzglFNO4Re/+AU7d0ZNaxs3bgSgf//+bNu27YPmqUTWGpXnSN6aqkKfxWXAw0AVcIeZLZN0LbDQzOqB24G7Quf3RqLkQlhuPlFHegtwqVn7nvhN6OPYE8o352sbCsmvMneutEnit7/9LVdccQWzZs2itraWsWPH8tOf/pTt27dz+OGHI4kf/ehH7Lfffrz00kudvt+Xv/xlVq9ezZFHHomZMXz4cH73u98xdepUlixZwqRJk+jTpw9nnHEGP/zhD/nSl77ExRdfTF1tDU89+Gvq+vWLahrWCu+3wj5Dc7etlfDM8UmTJlncH+R05fwlaTvsfvL5I4oQkXOlY8WKFRx8cGYnWBVxkrZnV9SnkeGoqlT7T9IiM5uUavmS6hx3zrmeKpvk0JnedTkdfpvMbzninHMuK17jcFmriKq+cy4tTxwua37XX+cqmyeOmPCrzJ1zpcITR0zkoonHm5Ccc4XgneNlpK0JKfknVTJxzuXWPvvsU+wQOvXTn/60/aLBnvLE4ZxzMZeL267nMnF4U5XLmvfHuJK3bims+ANseQsGjoaDz4T9Ds3JWy9YsIBrrrmGYcOGtd8x9+6770YSDQ0NfP3rX2fHjh3U1NTw5z//md69e3PJJZewcOFCqqurueGGGzjxxBOZM2cO999/P9u3b6e1tZULLrigw+sHHniAyy+/nKVLl7Jnzx6uueYapk2bRmtrK9/61rd46KGH6NWrFxdddBFmxtq1aznxxBMZNmwYjz/+eI+20ROHy5r3l7iStm4p/PUmqB0EA0ZGV1j/9Sb45OU5Sx7PPfccy5YtY8SIERxzzDH85S9/YcqUKUyfPp158+YxefJktm7dSl1dHTfeeCOSePHFF3nppZc49dRT258auHjxYl544QWGDBnCnDlzOrz+9re/zUknncQdd9zB5s2bmTJlCieffDJ33nknq1evZsmSJVRXV7ffdv2GG27g8ccfZ9iwYT3ePk8czrnKsuIPUdKoGxS9bvu94g85SxxTpkxh1KhRABxxxBGsXr2agQMH8uEPf5jJkycDMGDAACC67frll18OwEc/+lH233//9sRxyimndLhteuLrRx55hPr6en784x8D0NTUxJtvvsmjjz7KxRdfTHV1dHjf67brOeCJo4x4E5JzGdjyVlTTSFQ7ICrPkZqaD75zVVVV3eqTgK5vu/6b3/yGgw46qHtB9oAnjjLiTUjOZWDg6Kh5qq2mAdC0NSrPo4MOOoi3336bhoYGJk+ezLZt26irq+NTn/oU99xzDyeddBKvvPIKb775JgcddBCLFy/u9P1OO+00brrpJm666SYk8dxzz/Hxj3+cU045hVtvvZUTTzyxQ1NV223XvanKpeTXczjXiYPPjPo0IKppNG2Fps1w5Hl5/dg+ffowb948Lr/8cnbt2kVdXR2PPvooX/3qV7nkkkv42Mc+RnV1NXPmzOlQY0nne9/7HldccQWHHXYY77//PuPGjeOPf/wjX/7yl3nllVc47LDD6N27NxdddBGXXXYZM2fOZOrUqYwYMaLHneN+W/Uy5Ldod5Umm9uqA3kdVVWK/LbqMeQ1AOdiZr9DKzpR9JQnjgLwmwI658qJXznunCsLldDsng/d2W95TRySpkp6WdJKSVelmF8jaV6Y/4yksQnzrg7lL0s6LaH8nyQtk7RU0q8k1eZzG5xz8VdbW8uGDRs8eWTJzNiwYQO1tdkdRvPWVCWpCrgFOAVoBBok1ZvZ8oTFLgQ2mdl4STOAWcB0SYcAM4CJwAjgUUkHAvsBXwMOMbNdkuaH5ebkaztKkV/P4SrNqFGjaGxsZP369cUOpeTU1ta2X6yYqXz2cUwBVprZ6wCS5gLTgMTEMQ24JkzfB9wsSaF8rpk1A6skrQzv92aIuU7SHqAvsDaP21CSvMPdVZrevXszbty4YodRMfKZOEYCiZdiNgJHp1vGzFokbQGGhvKnk9YdaWZPSfoxUQLZBTxiZo+k+nBJM4GZAGPGjOn51vSA1wCcc+WkpEZVSRpMVBsZB2wGfi3pXDO7O3lZM5sNzIboOo6CBprEawDOuXKSz87xNUDiNfyjQlnKZSRVAwOBDZ2sezKwyszWm9ke4H7gk3mJ3jnnXEr5TBwNwARJ4yT1IerErk9aph44P0yfBTxm0bCIemBGGHU1DpgAPEvURPUJSX1DX8ingRV53AbnnHNJ8tZUFfosLgMeBqqAO8xsmaRrgYVmVg/cDtwVOr83EiUXwnLziTrSW4BLzawVeEbSfcDiUP4coTnKOedcYfi9qpxzzu2ls3tV+ZXjzjnnsuKJwznnXFZKajiuc7nidyx2rvs8cTig8g6kfsdi57rPE4cD/EDqnMuc93E455zLiicO55xzWfGmqhJRaX0Qzrn48sRRIrwPIrf8jsXOdZ8nDgdU3oHUa2nOdZ8nDgf4gdQ5lznvHHfOOZcVTxzOOeey4k1VJaLS+iCcc/HliaNEeB9E7vjQZud6xhOHqzg+tNm5nvE+Duecc1nxxOGccy4reU0ckqZKelnSSklXpZhfI2lemP+MpLEJ864O5S9LOi2UHSRpScLPVklX5HMbnHPOdZS3Pg5JVcAtwClAI9Agqd7MlicsdiGwyczGS5oBzAKmSzoEmAFMBEYAj0o60MxeBo5IeP81wG/ztQ2VzDuQnXPp5LNzfAqw0sxeB5A0F5gGJCaOacA1Yfo+4GZJCuVzzawZWCVpZXi/pxLW/TTwmpm9kcdtqFjl3IHsQ5ud65l8Jo6RwFsJrxuBo9MtY2YtkrYAQ0P500nrjkxadwbwq3QfLmkmMBNgzJgx3QjflSuvMTnXMyU5HFdSH+CzwNXpljGz2cBsgEmTJlm+Y/KmHedcpchn4lgDjE54PSqUpVqmUVI1MBDYkMG6pwOLzeydXAfdXeXctOOcc4nymTgagAmSxhEd9GcAX0haph44n6jv4izgMTMzSfXAvZJuIOocnwA8m7DeOXTSTOXSS1UzenHNFjD42KiB7WUNqzey+r2dHDthWKFDdM7FXN4SR+izuAx4GKgC7jCzZZKuBRaaWT1wO3BX6PzeSJRcCMvNJ+pIbwEuNbNWAEn9iEZqfSVfsZezVDWj5Wu3guhQvnrDDtZva9qrxuQdyM65vPZxmNkDwANJZd9PmG4Czk6z7nXAdSnKdxB1oLs8Onb8cBo37eQnnz+i2KE452LGrxx3zjmXlZIcVRVHfm2Ac65SeOLIER9y65yrFJ44KkyqmlGrGdjeQ4e9tuScS8UTR4XxmpFzMbRuKaz4A2x5CwaOhoPPhP0OLXZUaXnnuHPOFdO6pfDXm2DXZhgwMvr915ui8pjyGkcZi+NtUOIYk3NFteIPUDsI6gZFr9t+r/hDbGsdnjjKWBxvgxLHmJwrqi1vRTWNRLUDovKY8qYq55wrpoGjoWlrx7KmrVF5TGVU45A01Mw25DuYSuPNNvDkyvXsbG5lx+4Wrpy/pL28kvaBq3AHnxn1aUBU02jaCk2b4cjzihtXJzJtqnpa0hLgF8CDZpb325RXAm+2gZ3NrQyo6w10vFdWJe0DV+H2OxQ+eXnHUVVHnhfb/g3IPHEcCJwM/CPwH+EGhHPM7JW8Reacc5Viv0NjnSiSZZQ4Qg3jT8CfJJ0I3A18VdLzwFVm9lSnb+CKIpvboBSq2Sw5ph27WwDo28fHaThXKjLu4wDOBc4D3gEuJ3qWxhHAr4Fx+QrQdV82B/xsm826m2iS5105f0nKz3XOxVemp3lPAXcBnzOzxoTyhZJ+nvuwSlsldHp7/4xzlSvTxPFdM5ufWCDpbDP7tZnNykNcJS3Tg6rfUdf3gXOlKNPEcRUwP6nsaqJmKtdN5VL76Il87YNKqPU5VyydJg5JpwNnACMl/UfCrAFEj3TtlKSpwI1Ej469zcyuT5pfA9wJHAVsAKab2eow72rgQqAV+JqZPRzKBwG3AYcCBvyjd867ZN6U5lz+dFXjWAssBD4LLEoo3wb8U2crSqoCbiF6Pngj0CCp3syWJyx2IbDJzMZLmgHMAqZLOoTo+eMTgRHAo5IODM8dvxF4yMzOktQH8J7VHPAmI+dcpjpNHGb2PPC8pHvNbE+W7z0FWGlmrwNImgtMAxITxzTgmjB9H3CzJIXyuWbWDKyStBKYImk5cBzwpRDfbmB3lnG5FLJtvvFE41zlyrSPY4qka4D9wzoiurzjI52sMxJIvEtXI3B0umXMrEXSFmBoKH86ad2RwC5gPfALSYcT1YK+bmY7MtyOgqiEg6r3EzhXuTJNHLcTNU0tIupzKJZq4EjgcjN7RtKNRB3330teUNJMYCbAmDFjChpkJRxUvfPZucqVaeLYYmYPZvnea4DE2zuOCmWplmmUVA0MJOokT7duI9BoZs+E8vuIEsdezGw2MBtg0qRJfm+tHIt753Ml1PqcK5ZME8fjkv4NuB9oP800s8WdrNMATJA0juigPwP4QtIy9cD5RBcYngU8ZmYmqR64V9INRJ3jE4BnzaxV0luSDjKzl4FP07HPpGL1tAZQbjWIUozZuVKRaeJo65uYlFBmwEnpVgh9FpcBDxMNx73DzJZJuhZYaGb1RE1gd4XO741EyYWw3HyipNACXBpGVEF0u5N7woiq14ELMtyGstbTGkDcaxApldhzmp0rF5ne5PDE7ry5mT0APJBU9v2E6Sbg7DTrXgdcl6J8CR0TmCugtppJw+qNLF/7wcNn+tZUcez44YULpO05zbWDOj6n+ZOXd5o8yq1m5VwxdHUB4Llmdrekb6Sab2Y35CcsF1dtNZPla7e2P0cDYOuubEdr91A3n9NckjUr52KmqxpHv/C7f74DcaWlb5/qDslix+4WGjftLFzncwk+p9m5ctHVBYC3ht8/KEw4rlQcO2FYh9eNm3byk88fUbgABo6OmqfaahoQ++c0O1cuMn0eRy3R7UEmArVt5Wb2j3mKy2Wpp8NPS274agk+p9m5cpHpqKq7gJeA04BrgS8CK/IVlMteTzt2S65juASf0+xcucg0cYw3s7MlTTOzX0q6F/iffAbm4ilWNZNuPKc5VvE7V6IyTRxtvaCbJR0KrAM+lJ+QXL7kYihqydVMkuQifh/S6ypdpoljtqTBwHeJrvbehxT3h3Lx5kNRc8P3o6t0mSaOP5vZJuAJ4CMA4VYizjnnKkyvDJf7TYqy+3IZiHPOudLQ1ZXjHyUagjtQ0t8nzBpAwrBcV1m8jd+5ytZVU9VBwGeAQcCZCeXbgIvyFZQrvHTJ4LX12zlg+D4dyhpWb2T4PrUpLwJ0zpW/rq4c/72kPwLfMrMfFiimshG3M/POhqKm6/BtWL2R4w/sWL587VZ27m7JW5xx50N6XaXrsnM8PAPjc4Anji4kJ4qG1Rvp16eavn2qO5ydF+vMvLNkdeX8JQWMpLR5c5yrdJmOqvqLpJuBeUD78727eJBTxUk+a2+7g2zB7xzrnHN5lGniaLt73bUJZZ0+yMm5ZHFrunPOdU9eH+TkMhOHA+qLa7Z0eDBTm007du9V1remivXbmvdqcuuqjd8vnHOuPGR6d9x9ifo4RpjZ6ZIOAf7GzG7Pa3QVIhYHVAO0d3F1Va+94hg7tB+Txw7xWoJzFSrTpqo5wC+A74TXrxD1d3SaOCRNBW4keub4bWZ2fdL8GuBO4ChgAzDdzFaHeVcT3cq9FfiamT0cylcTDQduBVrMLLaPke1bU8XWXXvaH3LUJo6jbz42amDa5FXQ52w452Iv08QxzMzmh4M5ZtYiqbWzFSRVAbcApwCNQIOkejNbnrDYhcAmMxsvaQYwC5geajQziC4+HAE8KulAM2v7zBPN7L1MN7JQkodpjh3ar73cz86dc+Ui08SxQ9JQogYNJH0C2NLFOlOAlWb2elhnLjANSEwc04BrwvR9wM2SFMrnmlkzsErSyvB+T2UYb1F4cnAprVva8bkhB5/pzw1xJS3TxPENorviHiDpL8Bw4Kwu1hkJJD4AuhE4Ot0yoRazBRgayp9OWrftAdMGPCLJgFvNbHaqD5c0E5gJMGbMmC5CdUpcDX8AABU/SURBVIVQkRfOrVsaPamwdlD0jPRdm6PXn7zck4crWZmOqlos6XiiW5AIeNnMinVxwrFmtkbSh4A/SXrJzJ5IXigklNkAkyZNskIHmY04HFALEUNF1shW/CFKGm3PRm/7veIPnjhcycp0VNWlwD1mtiy8HizpHDP7z05WWwOMTng9KpSlWqZRUjUwkKiTPO26Ztb2+11JvyVqwtorcZSSOBxQCxlDHIYfF8yWt6KaRqLaAVF5Z7x5y8VYprdVv8jMNre9CM/m6Oomhw3ABEnjJPUh6uyuT1qmHjg/TJ8FPGZmFspnSKoJz/2YADwrqZ+k/gCS+gGnAksz3AYXE23Dj5N/UiWTkjdwNDQlXR/TtDUqT6eteWvX5o7NW+v8X93FQ6aJoyp0WgPtI6b6dLaCmbUAlwEPAyuA+Wa2TNK1kj4bFrsdGBo6v78BXBXWXQbMJ+pIfwi4NIyo2hd4UtLzwLPAf5nZQxlug3OFd/CZ0LQ5Ovjb+9Hvps1ReTqJzVvqFf2uHRSVOxcDmXaOPwTMk3RreP2VUNYpM3sAeCCp7PsJ003A2WnWvQ64LqnsdeDwDGN2rvj2OzTqCE9sdjryvM6bnbrbvOVcgWSaOL5FlCwuCa//BNyWl4icKzf7HZpd/8TA0VHNpK0jHbpu3nKugDIdVfU+8LPw45zLp4PPjPo0IKppNG2NmreOPK+4cTkXZDqq6hiiC/X2D+sIMDP7SP5Cc+WqkMOPS3IEV3eat5wroEybqm4H/glYRHSPKOe6rZAH7FjcQLI7sm3ecq6AMk0cW8zswbxG4tIqybNm51zZyjRxPC7p34D7gfYjmD8BsDBK9qzZubjyCyx7JNPE0XaPqaPCb+FPAHTOlSK/f1iPdZo4JH0jTP4x/DZgPfCkma3KZ2DOOZcXfv+wHuuqxtE/Rdn+wHckXWNmc/MQk3M5E4cbSLqY8Qsse6zTxGFmP0hVLmkI8CjgiSOmsulQL+fO91KP3+WBX2DZY5n2cXRgZhsT713l8qs7Z83ZdKh757urKH6BZY91K3FIOhHYlONYXBp+1uxcDvkFlj3WVef4i4THxSYYAqwF/le+goqrnjbplHOTUE74EElXKH6BZY90VeP4TNJrAzaY2Y48xRNrnTXpZJIUvEmoEz5E0rmS0VXn+BuFCqTUeVLoIR8i6VzJ6FYfh4u/bDrUM102r01tPkTSuZLhiaNMZXMgz3TZvNaqfIikcyUj00fHOpdf3XnEqnOuKPJa45A0FbgRqAJuM7Prk+bXAHcS3QNrAzDdzFaHeVcDFxLdxv1rZvZwwnpVwEJgjZkld+DnTWdNOqmacLJZv7vKZqSWD5F0rmTkLXGEg/stwClAI9Agqd7MlicsdiGwyczGS5oBzAKmSzoEmAFMBEYAj0o60MzangXydWAFMCBf8afS2YH4+gdXdJkUunUgX7eU//nDHdRsX8um3vvyQv/jWFd3QPt7l1WnvA+RdK4k5LPGMQVYaWavA0iaC0wDEhPHNKInCwLcB9wcrkifBsw1s2ZglaSV4f2ekjQK+FvgOuAbxERezu7DEFXbuYfW/iMY3rqdv91+HwsHnMuGfSbslRyeXLmenc1Rbt2xu4Ur5y8BSrD24ZyLtXwmjpFA4pCYRj64Pftey5hZi6QtwNBQ/nTSum1Dbn4K/DOpb8DYTtJMYCbAmDFjurcFxRaGqO6qaqK3etFcHVWwJmx8nA37TNhr8Z3NrQyo693+uq0mkpxgutu85TcMdM5BiY2qkvQZ4F0zWyTphM6WNbPZwGyASZMmJV/9Xhrah6g2tRc1V+1D/+a3e/S23W3e8lqLcw7yO6pqDZA4lnJUKEu5jKRqYCBRJ3m6dY8BPitpNdGdeU+SdHc+go+FgaOjIakJalq3s63mw0UKyDnn8ps4GoAJksZJ6kPU2V2ftEw9cH6YPgt4zMwslM+QVCNpHDABeNbMrjazUWY2NrzfY2Z2bh63objCENW61m1g71PTspXalm28OuTE9kXamo8aN+1kx+4Wtu7aw9Zde+jbp6Qqk+Vt3VJ4/F/hd1+Nfq9bWuyInOuRvB1dQp/FZcDDRMNx7zCzZZKuBRaaWT1wO3BX6PzeSJQMCMvNJ+pIbwEuTRhRVZbS9Tsc0utvObDvAqq2rWV97315of9prNszEjbt3KtP4sr5S1I2QbkiCgMcGt55n40ttdS98Sp1S37AgiEzWFd3gA9ccCUpr6elZvYA8EBS2fcTppuAs9Osex3RyKl0770AWJCLOAstVZJoWL2R4f1rOHb88A7lyzeN4cKLbmh/fVon7+ud1zEUBjhsbG1iQN/eQA1VLdUc2/IUTw/+WGkOm3YVz9sziiBV5/TytVvbh9J2V6Znrp5gCihPAxycKyZPHBUoTk0jZXPlezpt9+BK4AMcXKnzxOH2VsAHKmU6NLhkE0x4TGld6x6wQdS0bqe2ZRtLPzSt2JE5121+k0PXUdsDlXZt7vhApSKPBGpLMMk/mdwjrKjCPbh2VfWn/+53aK4ewMIR56a8gNO5UuE1jpjo26ea9dub9jrT7k6/Q4/OzivogUoFq8XsdyjPHXAxj7R91h4g/J29X8mVIk8cRZCqc3rssL5MHjc4JwesHt34sIIeqFTIG0TGujnNuSx54iiCWB9E/IFKzrkueOJwHYXOXCCqaTRtjR6odOR5efk4HxrsXOnxxOE6KvADlfzaE+dKjycOt7cYPlAp1s17zlUYTxxlyM/OM1Mu+6lkr3FxJcsTRxnyg0VmymU/ldXjg11J8AsAnXPOZcVrHDnmzQauK/4/4kqdJ44cK/lmgwLep6pSlfz/iKt4njjcB9ruU1U7qON9qj55eaySh5+xO1dcnjgK4MmV61m/rZkr5y/pUB67A12J3KfKz9g7KpfRYa50eOIogJ3NrfTrU73XwS52B7oKuk9Vj8WoSS9WJx+uIuQ1cUiaCtxI9Mzx28zs+qT5NcCdwFHABmC6ma0O864GLgRaga+Z2cOSaoEngJoQ+31m9n/yFX/FNYn4faoyUyJNes7lS94Sh6Qq4BbgFKARaJBUb2bLExa7ENhkZuMlzQBmAdMlHQLMACYCI4BHJR0INAMnmdl2Sb2BJyU9aGZP52MbutMkkqrZYMfuFobvU5vz+HIuR/epKvuE28MmPW9acqUunzWOKcBKM3sdQNJcYBqQmDimAdeE6fuAmyUplM81s2ZglaSVwBQzewrYHpbvHX4sj9uQtVQHxivnL0mZgGInR/epKvs+iB426ZVF8nQVLZ+JYySQ+E1qBI5Ot4yZtUjaAgwN5U8nrTsS2msyi4DxwC1m9kyqD5c0E5gJMGbMmJ5uS+WI4X2qkhX9jN2b9FyFK7nOcTNrBY6QNAj4raRDzWyv55qa2WxgNsCkSZOKWivp6kBX9k07OVb0fVLgW887Fzf5TBxrgMRTsFGhLNUyjZKqgYFEneRdrmtmmyU9DkwFivtA7C50daAr+6adclPgW887Fzf5TBwNwARJ44gO+jOALyQtUw+cDzwFnAU8ZmYmqR64V9INRJ3jE4BnJQ0H9oSkUUfU8T4rXxtQ9CYR11GMhsCWQpOec/mSt8QR+iwuAx4mGo57h5ktk3QtsNDM6oHbgbtC5/dGouRCWG4+UUd6C3CpmbVK+jDwy9DP0QuYb2Z/zNc2FL1JJIcK2RyWl4Qb0yGw3szoKlFe+zjM7AHggaSy7ydMNwFnp1n3OuC6pLIXgI/nPtLiuf7BFTSs3sjytVs7lPftU83YYbkbiVXI5rC8HDBjelW7NzO6SlRyneOFlu8zyvXbmunXp5oBdb07lG/dtafH711W/Kp252LDE0cXCnFG2bemaq9EsWN3i/elJPIhsM7FhieOGDh2/PC9yho37fQ28kQ+BNa52PAnALrS0DYEtm4QbF0T/fZ7QzlXFF7jqBCJI5322/Uah217gsF73qF5nxGwrrpHB+CCjSyK4RBYH7LtKpEnjiIr1IGn/QC+bin89REYPghqPxY1+fRwWGsljyyKZXNinK536a5y2IYy5omjC/k+sBf8wNNwO7z3KrTujvoKhh0UDXON2cOaXDfF9HqXrJTDNpQ5TxxdiOUZZYaSm5D22/UaX2p8iN29BzJy2GDY0wRvPQOjJvuw1nIR0+tdslIO21DmPHEUWgGr4MlNSJ/Y9hR7agazZ08rSNA7PCNk3Ysw4ZS8xOAKrByudymHbShznjgKqchV8P7Nb7O+7wSGbXkxqm1U1wAGuzZFCczFUlaDD8rhepdy2IYy54mjkIpcBd9W82FqWrayunYiI3tviL6MVX3gIyf26PNjN7KozDpWsxp8UA7Xu5TDNpQ5TxyFVOQq+KtDTmTS2rvZ3qsO9j/mgy/k5At79L6x6geq9I7VcrjlezlsQ5nzxFFIRa6Cb9hnAgtHnMu+ax6JLqIrxy9kPmp1pVaDieH1Llkrh20oY544CqnAVfBUTUiNjGT4ARdzWpxqCbmU61pdpddgnEvBE0chFbgKXpQmpGKfnee6VudDQ53biyeOQivnKngczs5zXauLwdDQ2A0+cBXPE4fLnTicnWdZq+tyqGsMhobGavCBc3jicLkUg7NzIKtaXZdDXX1oqHN7yWvikDQVuJHomeO3mdn1SfNrgDuBo4ANwHQzWx3mXQ1cCLQCXzOzhyWNDsvvCxgw28xuzHXchbjba1k+qzoGZ+c550NDndtL3hKHpCrgFuAUoBFokFRvZssTFrsQ2GRm4yXNAGYB0yUdAswAJgIjgEclHQi0AFea2WJJ/YFFkv6U9J49Voi7vZblHWXL9ey8nPulnOuGfD7IaQqw0sxeN7PdwFxgWtIy04Bfhun7gE9LUiifa2bNZrYKWAlMMbO3zWwxgJltA1YASW0jrmj8YUvOVYR8NlWNBBIbtxuBo9MtY2YtkrYAQ0P500nrdkgQksYCHweeSfXhkmYCMwHGjBnTzU34wJOvvsf67U1cOX9Jh/KSblrKBz87d67slWTnuKR9gN8AV5jZ1lTLmNlsYDbApEmTrKefuXN3C/36VO/VvFTSTUvOh7o61w35TBxrgMRe0VGhLNUyjZKqgYFEneRp15XUmyhp3GNm9+cndFcpvLboXPbymTgagAmSxhEd9GcAX0haph44H3gKOAt4zMxMUj1wr6QbiDrHJwDPhv6P24EVZnZDvgJPdRa6Y3dLTs9C/UzXOVeq8pY4Qp/FZcDDRMNx7zCzZZKuBRaaWT1RErhL0kpgI1FyISw3H1hONJLqUjNrlXQscB7woqS2zoZvm9kDuYw91VnolfOXpBwFlcvPKLhi3x7EOVeS8trHEQ7oDySVfT9hugk4O8261wHXJZU9CSj3kVagONwexDlXkkqyc7wYyq5pKQ63B3HOlSRPHBmKRdNSLsXl9iDOuZKTzwsAXZwNHB1d2Z2o1G8P4pwrCE8clergM6PbgezaDPZ+9Ltpc1TunHOd8MRRqfz2IM65bvI+jkrmtwdxznWD1zicc85lxROHc865rHjicM45lxVPHM4557LiicM551xWZNbjR1XEnqT1wBtZrDIMeC9P4eSKx5gbHmNueIy5EacY9zez4almVETiyJakhWY2qdhxdMZjzA2PMTc8xtwohRjBm6qcc85lyROHc865rHjiSG12sQPIgMeYGx5jbniMuVEKMXofh3POuex4jcM551xWPHE455zLSsUlDklTJb0saaWkq1LMr5E0L8x/RtLYhHlXh/KXJZ0WtxgljZW0S9KS8PPzIsZ4nKTFkloknZU073xJr4af82MaY2vCfqwvYozfkLRc0guS/ixp/4R5cdmPncUYl/14saQXQxxPSjokYV5cvtcpYyzk9zpjZlYxP0AV8BrwEaAP8DxwSNIyXwV+HqZnAPPC9CFh+RpgXHifqpjFOBZYGpP9OBY4DLgTOCuhfAjwevg9OEwPjlOMYd72mOzHE4G+YfqShL91nPZjyhhjth8HJEx/FngoTMfpe50uxoJ8r7P5qbQaxxRgpZm9bma7gbnAtKRlpgG/DNP3AZ+WpFA+18yazWwVsDK8X5xiLJQuYzSz1Wb2AvB+0rqnAX8ys41mtgn4EzA1ZjEWSiYxPm5mO8PLp4FRYTpO+zFdjIWSSYyJz0nuB7SNCorN97qTGGOn0hLHSOCthNeNoSzlMmbWAmwBhma4brFjBBgn6TlJ/y3pU3mIL9MY87FuNnr6ObWSFkp6WtLnchtau2xjvBB4sJvrdldPYoQY7UdJl0p6DfgR8LVs1i1yjFCY73XG/AmA5eVtYIyZbZB0FPA7SROTzmRcZvY3szWSPgI8JulFM3utWMFIOheYBBxfrBi6kibG2OxHM7sFuEXSF4DvAnnrF+quNDHG7ntdaTWONcDohNejQlnKZSRVAwOBDRmuW9QYQ3V7A4CZLSJqUz2wSDHmY91s9OhzzGxN+P06sAD4eC6DCzKKUdLJwHeAz5pZczbrFjnGWO3HBHOBttpPrPZjgvYYC/i9zlyxO1kK+UNUw3qdqBOsrYNqYtIyl9Kx43l+mJ5Ix06018lPJ1pPYhzeFhNRJ9waYEgxYkxYdg57d46vIurQHRym4xbjYKAmTA8DXiWpI7OAf+uPEx0oJiSVx2Y/dhJjnPbjhITpM4GFYTpO3+t0MRbke53V9hTzw4uywXAG8Er4R/9OKLuW6EwJoBb4NVEn2bPARxLW/U5Y72Xg9LjFCPwDsAxYAiwGzixijJOJ2nF3ENXYliWs+48h9pXABXGLEfgk8GL4cr8IXFjEGB8F3gl/0yVAfQz3Y8oYY7Yfb0z4bjxOwkE7Rt/rlDEW8nud6Y/fcsQ551xWKq2PwznnXA954nDOOZcVTxzOOeey4onDOedcVjxxOOecy4pfOe5cjklqJRp+2maumV1frHicyzUfjutcjknabmb7dLFMlZm1pnud6XrOFYM3VTlXIJJWS5olaTFwdorX54TnMSyVNCthve2SfiLpeeBvirYBzgWeOJzLvbqEh+4skTQ9Yd4GMzvSzOYmvgaeAGYBJwFHAJMT7ibbD3jGzA43sycLthXOpeF9HM7l3i4zOyLNvHlpXk8GFpjZegBJ9wDHAb8DWoHf5CNQ57rDaxzOFdaOLl6n0uT9Gi5OPHE4Fw/PAsdLGiapCjgH+O8ix+RcSt5U5Vzu1UlakvD6ITO7qrMVzOxtSVcR3RVVwH+Z2e/zGaRz3eXDcZ1zzmXFm6qcc85lxROHc865rHjicM45lxVPHM4557LiicM551xWPHE455zLiicO55xzWfn/3C7Okz+IOrkAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QKIFIRfZv3Mm",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "outputId": "5e7beeca-0902-4a20-872f-8e5ec923c958"
      },
      "source": [
        "plt.title(\"Uncertainty in positive class predictions\")\n",
        "plt.xlabel(\"Error\")\n",
        "plt.ylabel(\"Uncertainty\")\n",
        "plt.scatter(correct_positive_errors, correct_pred_uncertainties_positive, marker=\"s\", label='Correct', alpha=0.5)\n",
        "plt.scatter(incorrect_positive_errors, incorrect_pred_uncertainties_positive, marker=\"o\", label='Incorrect', alpha=0.5)\n",
        "plt.legend(loc='upper right');\n",
        "plt.show()"
      ],
      "execution_count": 41,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO29fXxcdZX4/z6ZyfNTkyb0KS0tULCUJ6GtIg9SVlpAEXYBWxRWEUF8QN316y66sj/W16rgKisL7rrsiqyg22JFKbtUKtrKFhBboPQRKZTSpm3atGmek8nM5Pz+uHemdyYzyUwyM5kk5/165ZW5n/u59557587nfD7nnM/5iKpiGIZhGKlSMNoCGIZhGGMLUxyGYRhGWpjiMAzDMNLCFIdhGIaRFqY4DMMwjLQwxWEYhmGkhSkOI6uIyNdE5D9zdK01IvLxXFwrFUTkhyJy1yD7c/Zs3OtdIiKNubpethGR9SLyKffzx0Rk7TDPk1fvzVhAbB5H/iMiCsxV1Tc9ZXcDp6jqjTmWZT3wmKpmvMEbrXvKBSJyCc5za5jIMmSS4byL4/kdyyU24jBSQhzsfTGygoj4R1sGI3WsIRgHREwQIvJlETksIgdF5GbP/lIR+Z6IvCMibSKyQURK3X3vFZEXRKRVRF5ze6WR49aLyDdF5HmgG3gUuAh4UEQ6ReRBt979IrJPRNpF5GURuchzjrtF5DH382wRURH5uIjsFZEjIvJ37r7Lga8By9xzvyYi14vIy3H3+tci8mSS5+A1XXzCvc/visgxEXlbRK4Y5BnuEZGvisgOt/6PRaTEs/9WEXlTRFpEZLWITHfLRUT+2X3u7SKyVUTOcPc9IiL/KCLlwBpguntvnSIyPe7ZrBGRz8fJ9JqI/IX7+V0i8hv3+n8SkY8Mci+1rvwH3Hv5VZJ6d4rIWyLS4d73n3v2nSIiv3fflyMisnKo+03yfXxbRP7o1n1SRGrdfZF34RYR2Qv8zi3/pIjsdOV+RkRO9JzvMhF53ZXpQUA8+z4hIhs82/M9z+uQOGbBAe+YR87Ie1MgIl8X57dyWER+IiLVcTIPeH/d/YtEZJN7r4dE5L5k39GYR1XtL8//AMUZXnvL7sYZpgNcAoSAbwCFwJU4DX2Nu/8HwHpgBuAD3gcUu9tH3foFwGXudr173HpgLzAf8LvnXg98Kk6WG4HJbp0vA01ASQI5Z7v38h9AKXA2EADmxdd1t4uBlsh+t+xV4NokzykqG/AJIAjc6t7zZ4ADuObZBMfuAbYBM4Fa4HngH919lwJHgHNdmR4AnnP3LQVeBibhNGTzgGnuvkc857gEaBzkO/xL4HnPvtOBVvd65cA+4Gb3Gb/blef0JPfyv8BKoMb9zt6fSAbgemC6+90vA7o8sv838HfuvhLgwqHuN8n3sR84w72HXyR4F37i7isFrgbedM/pB74OvODWrwM6gOvce/ornHfe+31vcD9XAgdx3sUSd/s9id6xBO/NJ10ZTgIqgCeAR1N8f18EbnI/VwDvHe22I1t/NuIYPwSBb6hqUFWfBjqB08QxL30S+KKq7lfVsKq+oKoBnAb/aVV9WlX7VfU3wCYcRRLhEVXdrqohVQ0murCqPqaqR90638Np7E4bRNZ/UNUeVX0NeA3nB5jovAGcBvBGcHqROD/e/0nxmbyjqv+hqmHgv4BpwJRB6j+oqvtUtQX4JnCDW/4x4GFVfcWV6avA+SIyG+e5VwLvwlFKO1X1YIryefklcI6nh/0x4An3eh8C9qjqj91n/CpOI3x9/ElEZBpwBXC7qh5z34ffJ7qgqv5cVQ+43/1KYBewyN0dBE4Epqtqr6pu8JSnc7+Pquo2Ve0C7gI+IiI+z/67VbVLVXuA24Fvu+cMAd/yPJMrge2qusp9D7+P00FJxIeAJlX9nit7h6q+NIiMXj4G3Kequ1W1E+e7Xi6xprRk728QOEVE6lS1U1X/kOI1xxymOMYGYZxelpdCnBc1wlH3xxahG6fXU4fT63orwXlPBK4Xx0zVKiKtwIU4DWyEfUMJJyL/zzUvtLnnqHavmwzvDz4iZzL+C/ioiAhwE/C425imQvQ6qtrtfhzsWt57fQenN477/x3PuTpxRmYzVPV3wIM4o7rDIvKQiFSlKF8UVe3AGSksd4tuAH7qfj4ReE/c9/QxYGqCU80EWlT12FDXFJG/FJHNnnOewfHv7W9wRhR/FJHtIvJJV8507zf+mRYS+254958I3O+Rp8WVYQbOdxCtq063Ptm7OZPE73sqxHzX7mc/sR2OZO/vLcCpwOsislFEPjRMGfIeUxxjg704PW0vc4h9wZNxBOgFTk6wbx9Oj3CS569cVe/x1IkPu4vZFsef8TfAR3BMY5OANjz25zQYEOLn9tr6cHwrH8Xxs2SLmZ7Ps3BMW7j/vbb2chzT3H5Xxn9R1fNwzEunAl9JcO5Uwhf/G7hBRM7HUfbr3PJ9wO/jvqcKVf1MgnPsA2pFZNJgF3J78f8BfB6Y7H5v23C/N1VtUtVbVXU68GngX0XklDTuN0L8Mw3ivJMRvM9lH/DpuPssVdUXcExP0XO5HQnvueOfwUlJ9g31PcR8167MIeDQEMehqrtU9QbgBOBeYJX7row7THGMDVYCXxeRBtd59wHgKmDVUAeqaj/wMHCfOA5Zn4icLyLFwGPAVSKy1C0vEcfRPli45iFif5SVOD+sZsAvIn8PpN3j9px7tgyM3voJTi836DGZZIPPuc+4Fse+v9It/2/gZhE5x31u3wJeUtU9IrJQRN4jIoU4PoJeoD/BuQ8BkyOO1iQ8jdNofQNY6X534JjmThWRm0Sk0P1bKCLz4k/gmo3W4DT0NW7dixNcqxynEW0GECeYIurkFicwIfIeHHPr9qdxvxFuFJHTRaTMva9VrukwET8EvuqaJBGRahGJmOP+F5gvIn/hmo2+QOIRFzjPa5qIfElEikWkUkTe4+5L9o5F+G/gr0RkjohU4HzXK+NG8wkRkRtFpN793lrd4sGezZjFFMfY4BvAC8AGnB/xd4CPqeq2FI//f8BWYCPO8P9eoEBV9+E4JL+G04Dsw+k9DvZe3A9cJ07Uy78AzwC/Bt7AGQH1koJ5Kwk/d/8fFZFXPOWP4jRqjw3zvKnyM2AtsBvH1PGPAKr6LI59/hc4Pd+TOW5SqsLpuR/Duf+jwD/Fn1hVX8dplHa7ppjpCeoEcJyxH3BliZR3AEvcax7AMZXci+NLSsRNOD3714HDwJcSXGsH8D0ch+4h4EycgIAIC4GXRKQTWI3jI9ud6v16eBQnSKAJZxT1hWQVVfWX7n2tEJF2nBHQFe6+Izg+nXvca86Nk9d7ng6cQI+r3OvuAha7u5O9YxEedmV+Dngb532+Y5D783I5sN19ZvcDy13fzbjDJgAaeY84ocOHgXNVdVeWrrEHJ7Lm2WycfyIiWZwsaowuNuIwxgKfATZmS2kYhpEeNlvTyGvckYAA14yyKIZhuJipyjAMw0gLM1UZhmEYaTEhTFV1dXU6e/bs0RbDMAxjzPDyyy8fUdX6RPsmhOKYPXs2mzZtGm0xDMMwxgwiknSCsZmqDMMwjLQwxWEYhmGkhSkOwzAMIy0mhI/DMIzxTTAYpLGxkd7e3tEWZcxRUlJCQ0MDhYXxCbiTY4rDMIwxT2NjI5WVlcyePRsnca6RCqrK0aNHaWxsZM6cOSkfZ4rDyAn3rNlJc8fAZTTqK4u584oBSV4NIy16e3tNaQwDEWHy5Mk0NzendZwpDiMnNHcEaKgpG1DeeKw7QW3DSB9TGsNjOM/NnOOGYRhGWpjiMAzDyBBNTU0sX76ck08+mfPOO48rr7ySN954IyfXfuSRRzhw4MDQFTOAKQ7DMIwMoKr8+Z//OZdccglvvfUWL7/8Mt/+9rc5dGjIVWcJhUKDbqdCLhWH+TgMw5hQZCtQY926dRQWFnL77bdHy84++2xUla985SusWbMGEeHrX/86y5YtY/369dx1113U1NTw+uuv89BDD8Vs79y5kzvvvJP169cTCAT43Oc+x6c//WkA7r33Xh577DEKCgq44oorWLBgAZs2beJjH/sYpaWlvPjii5SWlg77XobCFIeRE+orixM6wusrk61+ahjZIVuBGtu2beO8884bUP7EE0+wefNmXnvtNY4cOcLChQu5+GJnGfhXXnmFbdu2MWfOHNavXx+z/dBDD1FdXc3GjRsJBAJccMEFLFmyhNdff50nn3ySl156ibKyMlpaWqitreXBBx/ku9/9LgsWLBjRfaSCKQ4jJ1jIrTFR2bBhAzfccAM+n48pU6bw/ve/n40bN1JVVcWiRYti5k94t9euXcuWLVtYtWoVAG1tbezatYtnn32Wm2++mbIyR/nV1tbm/J5McRiGYWSA+fPnRxv5VCkvL0+6rao88MADLF26NKbOM888M3whM4Q5xw3DMDLApZdeSiAQ4KGHHoqWbdmyhUmTJrFy5UrC4TDNzc0899xzLFq0aMjzLV26lH/7t38jGAwC8MYbb9DV1cVll13Gj3/8Y7q7HdNaS0sLAJWVlXR0dGThzgZiIw7DMIwMICL88pe/5Etf+hL33nsvJSUlzJ49m+9///t0dnZy9tlnIyJ85zvfYerUqbz++uuDnu9Tn/oUe/bs4dxzz0VVqa+v51e/+hWXX345mzdvZsGCBRQVFXHllVfyrW99i0984hPcfvvtOXGOT4g1xxcsWKC2kJNhjF927tzJvHmp+dEs/c1AEj0/EXlZVRN62m3EYRjGhGKiKodMYj4OwzAMIy1McRiGYRhpYYrDMAzDSAtTHIZhGEZamOIwDMMw0sIUh2EYRgaoqKgYbREG5fvf/3500uBIMcVhGIaR52Qi7XomFYfN4zAMY+LRtA12PgVt+6B6Jsy7CqaekZFTr1+/nrvvvpu6urpoxtzHHnsMEWHjxo188YtfpKuri+LiYn77299SWFjIZz7zGTZt2oTf7+e+++5j8eLFPPLIIzzxxBN0dnYSDoe5+eabY7affvpp7rjjDrZt20YwGOTuu+/m6quvJhwO87d/+7f8+te/pqCggFtvvRVV5cCBAyxevJi6ujrWrVs3ons0xWEYxsSiaRu88ACUTIKqGdDT6my/746MKY9XX32V7du3M336dC644AKef/55Fi1axLJly1i5ciULFy6kvb2d0tJS7r//fkSErVu38vrrr7NkyZLoqoGvvPIKW7Zsoba2lkceeSRm+2tf+xqXXnopDz/8MK2trSxatIgPfOAD/OQnP2HPnj1s3rwZv98fTbt+3333sW7dOurq6kZ8f6Y4DMOYWOx8ylEapZOc7cj/nU9lTHEsWrSIhoYGAM455xz27NlDdXU106ZNY+HChQBUVVUBTtr1O+64A4B3vetdnHjiiVHFcdlll8WkTfdur127ltWrV/Pd734XgN7eXvbu3cuzzz7L7bffjt/vNO/ZSLtuisMwjIlF2z5npOGlpMopzxDFxccXKPP5fMPyScDQadd/8YtfcNpppw1PyBFgznHDMCYW1TOhtz22rLfdKc8ip512GgcPHmTjxo0AdHR0EAqFuOiii/jpT38KOKnT9+7dm5IyWLp0KQ888ACRRLWvvvoq4IxK/v3f/z2qrLKRdt0Uh2EYE4t5V0Fvq+Pb0H7nf2+rU55FioqKWLlyJXfccQdnn302l112Gb29vXz2s5+lv7+fM888k2XLlvHII4/EjFiScddddxEMBjnrrLOYP38+d911F+CkY581axZnnXUWZ599Nj/72c8AuO2227j88stZvHjxiO/F0qobhjHmSSetOpDVqKqxSF6lVReRy4H7AR/wn6p6T9z+YuAnwHnAUWCZqu5x930VuAUIA19Q1Wfc8r8CPgUosBW4WVV7s3kfhmGMM6aeMaEVxUjJmuIQER/wA+AyoBHYKCKrVXWHp9otwDFVPUVElgP3AstE5HRgOTAfmA48KyKnAlOBLwCnq2qPiDzu1nskG/dgC74YhmEMJJsjjkXAm6q6G0BEVgBXA17FcTVwt/t5FfCgiIhbvkJVA8DbIvKme769rsylIhIEyoAD2bqB5o4ADTVlA8obj2Vm9qVhGJlDVXGaDyMdhuOuyKZzfAbgjW9rdMsS1lHVENAGTE52rKruB76Lo0AOAm2qujbRxUXkNhHZJCKbmpubM3A7hmHkKyUlJRw9enRYjeBERlU5evQoJSUlaR03puZxiEgNzmhkDtAK/FxEblTVx+LrqupDwEPgOMdzKqhhGDmloaGBxsZGrJOYPiUlJdHJiqmSTcWxH/AGRje4ZYnqNIqIH6jGcZInO/YDwNuq2gwgIk8A7wMGKA7DMCYOhYWFzJkzZ7TFmDBk01S1EZgrInNEpAjHib06rs5q4OPu5+uA36kz1lwNLBeRYhGZA8wF/ohjonqviJS5vpA/A3Zm8R4MwzCMOLI24lDVkIh8HngGJxz3YVXdLiLfADap6mrgR8CjrvO7BUe54NZ7HMeRHgI+p6ph4CURWQW84pa/imuOygb1lcUJHeH1lUNPzjEMwxiv2ARAwzAMYwCDTQC0lCOGYRhGWpjiMAzDMNLCFIdhGIaRFqY4DMMwjLQwxWEYhmGkhSkOwzAMIy1McRiGYRhpMaZyVRkTB0tpbxj5iykOIy+xlPaGkb+YqcowDMNIC1MchmEYRlqYqcpIm/HofxiP92QY2cIUh5E249H/MB7vyTCyhSkOI6ek2rNPJ6W9jRYMI7eY4jBySqo9+3QafBstGEZuMcVhZJzBRgCGYYx9THEYGcdGAIYxvjHFYaTNeFxSdzzek2FkC1McRtqMR4fzeLwnw8gWpjiMnJKNnr2NFgwjt5jiMHJKNnr2NlowjNxiisPIODYCMIzxjSkOI+PYCMAwxjeW5NAwDMNIC1MchmEYRlqY4jAMwzDSwhSHYRiGkRamOAzDMIy0MMVhGIZhpIUpDsMwDCMtTHEYhmEYaWETAHOArVBnGMZ4whRHDrD1KXKLKWrDyC6mOIxxhylqw8gupjjGEdbTNgwjF5jiGEdYT9swjFxgUVWGYRhGWmRVcYjI5SLyJxF5U0TuTLC/WERWuvtfEpHZnn1fdcv/JCJLPeWTRGSViLwuIjtF5Pxs3kMmiKxPEf9n61MYhjEWyZqpSkR8wA+Ay4BGYKOIrFbVHZ5qtwDHVPUUEVkO3AssE5HTgeXAfGA68KyInKqqYeB+4Neqep2IFAEDbTN5hvkXjpMLP4wtJGUY2SUlxSEik1X1aJrnXgS8qaq73XOsAK4GvIrjauBu9/Mq4EEREbd8haoGgLdF5E1gkYjsAC4GPgGgqn1AX5pyGaNILvwwpqgNI7ukOuL4g4hsBn4MrFFVTeGYGcA+z3Yj8J5kdVQ1JCJtwGS3/A9xx84AeoBm4McicjbwMvBFVe2Kv7iI3AbcBjBr1qwUxB37WE/bMIxckKriOBX4APBJ4F9E5HHgEVV9I2uSJcYPnAvcoaovicj9wJ3AXfEVVfUh4CGABQsWpKLoxjzW0zYMIxekpDjcEcZvgN+IyGLgMeCzIvIacKeqvpjgsP3ATM92g1uWqE6jiPiBauDoIMc2Ao2q+pJbvgpHcRgTHJvDYhi5I2UfB3AjcBNwCLgDWA2cA/wcmJPgsI3AXBGZg9PoLwc+GldnNfBx4EXgOuB3qqoishr4mYjch+Mcnwv8UVXDIrJPRE5T1T8Bf0asz8SYoEzEOSymLI3RIlVT1YvAo8A1qtroKd8kIj9MdIDrs/g88AzgAx5W1e0i8g1gk6quBn4EPOo6v1twlAtuvcdxlEII+JwbUQWO0vqpG1G1G7g5jfs1Rhnzw2SOiagsjfwgVcXxdVV93FsgIter6s9V9d5kB6nq08DTcWV/7/ncC1yf5NhvAt9MUL4ZWJCi3EaeYT1hwxj7pDoBMJEf4auZFMQwDMMYGww64hCRK4ArgRki8i+eXVU4JiTDyAlmzzeM/GEoU9UBYBPwYZw5ExE6gL/KllCGEc9Q9nzznRhG7hhUcajqa8BrIvIzVQ3mSCbDSJuJOOowZWmMFqk6xxeJyN3Aie4xgjO946RsCWYYxuBMRGVp5AepKo4f4ZimXgbCQ9Q1DMMwxjGpKo42VV2TVUkMwzCMMUGqimOdiPwT8AQQDW1R1VeyIpVhxDEce75FYhlGdkhVcUSy2non3ilwaWbFMYzEDKeht5nVhpEdUk1yuDjbghiGYRhjg6EmAN6oqo+JyF8n2q+q92VHrPzFzB+GYUx0hhpxlLv/K7MtyFjBzB9jg3vW7GTjnhZ2HGiPKS8r8jO7Lu9XGzaMvGaoCYD/7v7/h9yIYxiZobkjQHmRn6rSwpjy9h6bx2oYIyXV9ThKgFuA+UBJpFxVP5klufKCRGapjXta2HO0iwtPqR8lqYxUKSv2DVAUXX0hm1ltGCMk1aiqR4HXgaXAN4CPATuzJVS+kMgsteNAO90BmwOZTyRT8PWVxSyZPzWmvPFYt/miDGOEpKo4TlHV60XkalX9LxH5GfB/2RTMMFLFFLxh5JZUFUdkvN8qImcATcAJ2REpvykr8tPc2TvAGW7mD8MwJgqpKo6HRKQG+DrOOuEVwF1ZkyqPuXBuHY3HuvneR84ZbVGMQTAFbxjZI1XF8VtVPQY8B5wEICJzsiZVnrFh1xG6+46vW9XVF+LLj2+2uRt5jCl4w8geqS4d+4sEZasyKUg+EsmP1NzZO6C8oaYs4URAwzCM8c5QM8ffhROCWy0if+HZVYUnLHe8EhlNfPnxzQkn/RkjI1Oz8G1BI8PILUOZqk4DPgRMAq7ylHcAt2ZLKGNikKlZ+GYuNIzcMtTM8SdF5H+Av1XVb+VIJiMNMp47q2kb7HwK2vZB9UyYdxVMPSMDkhqGMV4Y0sehqmHgmhzIYgyDSK89/m9Y/pembfDCA9DTClUznP8vPOCUG4ZhuKQaVfW8iDwIrAS6IoUTZSGnCWND3/kUlEyC0knOduT/zqds1GEYRpRUFUckpvEbnrIJs5DThLGht+1zRhpeSqqccsMwDBdbyMk4TvVMxzwVGWkA9LY75VlgwozkDGOckWp23CnAt4DpqnqFiJwOnK+qP8qqdEZumXeV49MAZ6TR2w69rXDuTVm53IQZyRnGOCPVCYCPAM8A093tN4AvZUMgIz0ivfb4v2H12qeeAe+7wxlxtO93/r/vDvNvGIYRQ6o+jjpVfVxEvgqgqiERsdSjeUAqvfa0QnannmGKwjCMQUlVcXSJyGQchzgi8l6gLWtSGRllrC53a+u7G0Z+kqri+GucrLgni8jzQD1wXdakMgzGrsIzjPFOqlFVr4jI+3FSkAjwJ1W1xZtzjPXADcPIB1KNqvoc8FNV3e5u14jIDar6r1mVzoghX3rgpsCMCYWl4RlAqlFVt6pqa2TDXZvDkhxOUDKa5sQw8hlLw5OQVH0cPhERVY04x31AUfbEMtJh6/42vvz45gFlKJzZUM3WxjY27mkBoKTQx5kzqoGJMdHORkfGiLA0PAlJVXH8GlgpIv/ubn/aLTPygN5geIAJa8eBdhCio4EIY2lVvEzMLM8X854xRrE0PAlJVXH8LY6y+Iy7/RvgP4c6SEQuB+4HfMB/quo9cfuLgZ8A5wFHgWWqusfd91XgFiAMfEFVn/Ec5wM2AftV9UMp3oMxxrARgTHq5DgNz1gh1aiqfuDf3L+UcBv3HwCXAY3ARhFZrao7PNVuAY6p6ikishy4F1jmpjRZjrP64HTgWRE51U3xDvBFYCfOSoQThmQ98BK/bxSkMYwJQI7T8IwVUo2qugC4GzjRPUYAVdWTBjlsEfCmqu52z7ECuBrwKo6r3fOCs4b5gyIibvkKVQ0Ab4vIm+75XhSRBuCDwDdx5pdMGJL1wOP9G9nGkhMaE4ZIGh5vVNW5N01o/wakbqr6EfBXwMs4pqNUmAF4DYGNwHuS1XHTmLQBk93yP8QdGzE0fh/4G6AyRTmMDGMmJGNCYWl4BpCq4mhT1TVZlSQFRORDwGFVfVlELhmi7m3AbQCzZs3KgXSjR6IRQFgVdKATeKKNCmx0ZBiZJ1XFsU5E/gl4AojGNg6xAuB+wOtBanDLEtVpFBE/UI3jJE927IeBD4vIlUAJUCUij6nqjfEXV9WHgIcAFixYoKnc5FjFRgDJGc1nY6HAxnglVcURMTGd5/4Xhl4BcCMwV0Tm4DT6y4GPxtVZDXwceBEn99XvVFVFZDXwMxG5D8c5Phf4o6q+CHwVwB1x/L9ESsMw8gELBTbGK4MqDhGJOJ//x/2vQDOwQVXfHuxY12fxeZx1PHzAw6q6XUS+AWxS1dU4vpNHXed3C45ywa33OI4jPQR8zhNRZRjGMLFRkJEJhhpxJHJAnwj8nYjcraorBjtYVZ8Gno4r+3vP517g+iTHfhMncirZudcD6we7vmEYsdgoyMgEgyoOVf2HROUiUgs8CwyqOIzRZzR7mNa7NYzxSao+jhhUtcWdb2HkOaPZw7TerWGMT4alOERkMXAsw7IYw2Swnn2+MhFGIxYKbIxXhnKOb8VdLtZDLXAA+MtsCWWkx1js2Y9FmdNlvChAw4hnqBFHfAJBBY6qaleW5DHynIkwUhgp+fyMbBRkZIKhnOPv5EoQI79I1vht3d/GFWdMG1A+nkYKIyWfR1OjrbiM8cGwfBzG2GG4PcxkjV9kQaiRXttWCzSMsYspjnHOaPYwB7t2rjP6GoaROUxxjAPGot16KJlH4ifIZx+DYYwHTHGMMbLVKMafd+OeFnYcaKesyM+Fc+uGfd5kDCXrSPwE+exjMIzxgCmOMUa2GsX48+440E5VaSHtPcGYeiV+35gb3eSasTgCNIx0MMUxgUhnomBZsY/2niBdfaGYRvDSeSeYuWcI7PkY4x1THBOIwUYrW/e3seNA+4B9JYU+vveRc3Ih3sSmaVvs8qTzrrJV54y8xRSHAUBvMMy06tIB5QfbekZBmuFzz5qdUf+Ml7JiH7Mnl4+SVEPQtA1eeABKJkHVDOhpdbbfd4cpDyMvMcVh5CUjmX9SX1FCd19oQPnC2bUZlTFj7HzKURqlk5ztyP+dT5niMPISUxxjjGw5Xkv8vgGO8Ej5aDASP0GiKLDGY93563to2+eMNLyUVDnlhpGHmNN0NaMAACAASURBVOIYY6Ta+CVyhG/c08Keo11ceEr9gPpnNlRbCOtoUT3TMU9FRhoAve1OuWHkIaY4ximJHOF7jnTT3NE7QBlkKgVIPk682/BmM92BMF19oZjZ6nk1GXDeVY5PA5yRRm879LbCuTcN+5T5+F0Y4wdTHHlOJhuAC+fW0XisO2GU1D1rdo7YBJaPE++6A2GqSgsBYmTLq5HU1DMcR7g3qurcm0bk38jH78IYP5jiyHNy1QCMl15ovA+oy3WSlxXl+as+9YzcOsIt/NcYAXn+a8pPvKOArfvb6A2GAceRfGZDNTC6JoGkIalFfmbXDVRC44n4Z/7lxzcnVLwTGgv/NUaIKY5h4B0F7DjQHp3/0N4TjJZHer25sDUnyjPV2h0kGFJmTT7eaCaKmjJchuiBjyufgYX/GiPEFEeWyYWpKVGeqWC4n65AKEZZdPWFeKu5M2FK86EawHHVcMaTQg98XPkMLPzXGCGmOMYps2rLae8JsmT+1GhZpJEbTgOYSsMZ71+ImPFK/L5Ri2hKad7LOOyBD3rfJRb+a4wMUxx5zljKtJqqfyGXvfSUFNQ47IEPet9N4YyH/xoTC1Mcec6YNwONBSbaBLwshP8aEwtTHIOQzK7/VnNn9HNYNZoI0LtWxWiOCJKlRE9nol+uF3YaVbIwAS/vyXX4rzGuMMUxCMns+kDKqcYjpiZv2C4Qtftnwt4fb86KZIFNdO5U1/pOdWGnkZA3DvcUeuDDNhnafAljHGKKI0MM1gh+7yPnZNXen04jO9wGcLBRTCJSmUuSdqRSNhvhIXrgw1JkNl/CGKeY4sgQYyVcc6gG8JofbKClq48jnQH8BQXR8tKiApYvPDFpypJ4mjsC1FcW0x1wRlmH2nsJ9Suh/n46AkG+/PhmJ+nike7UTF9jsREeh9FahgGmOIw4Wrr6mFZdSncgTHHhccXR0Rsa5KjEeLPwrt3eFDV1RUKEdxxoH7BuRlJSaIRH0/SV6NrL9r9GX/lULpzrKRzj0VqGAaY4jFEkMgpZu70ppjzcrwMrpxAyO5qjvkTXDnc00N/ZAkw5Xjieo7WMCYMpjkEYS3MoomTAD7D3aDedgRDtvccb8HC/smbrQS6dd0LGRA31K0W+gmj22ggJl6sdgyGzu2oXM6/9EUfuiRKtZUwITHEMQibNGzlRQsP0A3jNLEc6A4TCiq9AKCoooLrMadQ7ekOc2VCd0Wfi9wl94f4BkVolhQlWHRyDIbNHK+ayvnY5C0vfsPkSxrjCFMcQpGo3H0ox5CS8dJjOWK+ZxV9QgGo//gIhlMhklCJJ05sXH1cKUypLAGLSosBA05LzHYSY2rOEsw48R03wHY4VTqFx6ge5Jc8b4abSk2HxtaMthmFkFFMcQ5Cq3XzwFA85iuXPQOqM0iIf3Z0hwv1Cv2rUKV5alN7a4/HPw6uAI8+uqy9EfUXJkOeKfgc1Z7KDM6Pl+RaxZhgTBVMc2SaXYaQZ8AMsXziLtdubaO3poysQZsak0ui+jXtauGfNzmGNnhIdk8paGUnngxT7ohMdI2R1kl7TNv7vqYcp7jzAscIpbKm82BlNuOcfk/4wwxgmpjiySdM2WPM30HUEyuug7jSocJ3L2Yjlz6AfIBhS+kL9MWa6vnA/T289CGTG9JassfWmfo+sLdIdCNMTDFPq+j/6Wh3ZvLPvszZJz62j3UHCldOpD3fywc5VbKq6kaMVc1Oe22IY44WsKg4RuRy4H/AB/6mq98TtLwZ+ApwHHAWWqeoed99XgVuAMPAFVX1GRGa69acACjykqvdn8x6GTaRB6joCpbUQ7IV9L8HM9zhKJBux/BlKXldW7KOvtZ9Qf39MeXmRn/Iif8r5ruJJ1V/kHYnsONAenVPS3htkUtnxCKzyIj8NNWUjM1ml4hdy6/T4eimUAgL+KgDmtqzjaMXcBCc1jPFN1hSHiPiAHwCXAY3ARhFZrao7PNVuAY6p6ikishy4F1gmIqcDy4H5wHTgWRE5FQgBX1bVV0SkEnhZRH4Td878INIgldc5SqPQteUf+RP4iqLmo4xPWhtG8rpEua6aOwL4C4STT6iIqTuSXFWD+Yu8z8FrmjrU3kuBSEz9k8J7uLD/ReZ0t1C59yQ2+M8HhtnjT8UvFK3TGy0K+CqoDBwc3jUNY4yTzRHHIuBNVd0NICIrgKsBbyN/NXC3+3kV8KCIiFu+QlUDwNsi8iawSFVfBA4CqGqHiOwEZsSdM6MM23YdaWzqToPGPzplviJnBOIxH+VDqpJk/gevX2Hv0W6C4X76wv109YWipqRUFdxQvor45Xgjczv2t/ZQXuQjEOynX5WGwG6W6VO0Ukarv566UDuXtK+AptOGZ/pLxS8UqeOhONxJR/G09K9nGOOAbCqOGYDXHtMIvCdZHVUNiUgbMNkt/0PcsTHdQhGZDbwbeCnRxUXkNuA2gFmzZg3zFgax5Tdtg3XfTu5QjTQ2FSdAwyJnpBHxdeRzfqUkBMP90RQkERMRpK7gmjsClBf5B0z2S2UEE1k3fefBdq4rfxXpnwThMmZOriAABAtCw/cZpeIXcuuUhoOgkygOd1IS6mDbCVenf73hRNhZhl0jzxiTznERqQB+AXxJVdsT1VHVh4CHABYsWDD8CQmJSMWh6m2Qyuuc0UZv67CVRq7zMNVXFkfnXoDjGAco9EuyQ7KC3ydR5RLq76e6r4lm6igsPC5HT0H58H1GqfiF3Dry1MP4Og7QXDiFLZVLaQrOgGPdqUdODSfCbiwmdzTGPdlUHPsBbxxog1uWqE6jiPiBahwnedJjRaQQR2n8VFWfyI7ox0nUYC85/GNqff0sfNcgCfde9eHriJ2wtqVyCeFXfdx5RfpyJDNprdl2MCsKJXJs5NxdfSHKi5zXpaxo5K/N3pYuJ3LLNX2BY6IqK/ZRVuSPKouKYj+nT3ec0cWFBajOZGa4gx5fabROrb8Xqk8avjCp+IWmnsFFt94X3Vw6nOsMZ4JmLjLs2ojGSJNsKo6NwFwRmYPT6C8HPhpXZzXwceBF4Drgd6qqIrIa+JmI3IfjHJ8L/NH1f/wI2Kmq95EDEjXYM44cYX+oNrZiooR702MnrPmBpgz7LnqD4Zys85HKnIuh8CqErkCYIl8Bc3mHj7CZmlATbQVTeT3QwAfqWqjkIB3F09jgP5+vfOTS4ydpmuT2wEs8pqVep7EbJVIeDQ5ngma210O3EY0xDLKmOFyfxeeBZ3DCcR9W1e0i8g1gk6quxlECj7rO7xYc5YJb73Ecp3cI+JyqhkXkQuAmYKuIRJay+5qqPp2t+wBYsXEvPX3HV+8rpoQKPcyKjSGWL3T9J8NMuDeSiWOTO3fxib4nWbSri47iaTSXnER9724qAwfZ318HTbfm1Y/fu+7G2u1NzPft4/L2p+ntr2JPXw3T+/bw8f7f8HrgNP4kU6kuOMjlVT+PdXzn4XrZKQc4DGeCZraTO9qaIcYwyKqPw23Qn44r+3vP517g+iTHfhP4ZlzZBiC3Rnagpy9MZcnxR/Vi+H0sCz5JT6ANtD/qUP1R9wXs8Exci0QQlRX7Ytam8DJcc9Lkzl0sOPAY21XoKJpKdfc7zDv8NAcrz6C1pIHSrsz2HL0KbmrPW5zV8Rw1wUMEKqZDk3/IayTKXXVO+P/o9VXRqmUUFISZ2n+MroJSpsoxDvsbaA2X0+vXBGtuhADX3tcMNIeorxzejPacMpwJmtlO7pjtEY0xLhmTzvHRZrdvNo+GrmKpbIT2/dFe747nQglDSjO5TneEuS3r6PVX0imAFFAZPEKfr4yKvmZaS2fR46t0zDkZ6jlGG+WmbfDCWqifBCVnOg1ZCgrK26jfs2YnACfoYZq1jmM9fYRVKZUuOqSEajrp7g+hQEuoNPdrbmTL5j+c0VK2R1hjMF29MfqY4hgmbxXM5pGik/nkNR77O5uT1h+MVGzk8T12X3sjzf56/AV9ABSFOgkUlFMc6jh+gmz0HDNg2oiE5rYVTGVSfyct+CkQ6JBSqqSbLqnAV+CkXC/t74LqHM7OzrbNfxgTNId1TKqMwXT1xuhjiiMB3oZ8a2MbXYEQvcEwBSKcXdjIpbzEDDnC4b4ToOmEhD/qGEdwXyja6CfyXfgO7+C60ItUBhyH8K7axdEcSBHie+zHCqdQGu4AKeT1g+2UhAuppJVASTXtPUHKi30j6jkmU2YfP7yV9qIT6OrzrNqn/dSE3uF3vambi8qKffwu9B6uDf6SSopo1zKatIbpvhb2MgvRfirporQ/zI+Ozk9sAizyp7ZeeSKSjSryxeafq0inPPQZGfmPKY4EeM0hDTVl/Nv6N/EVCCf37+EmXUMH5RzUyZSEk5tpvA3aoEnwmrZxScsKfGU1dBRNoTjUzoIDj7Fp+o00MiPhIc0dAQ7NWMKCA49RUesn4KuguGcGk9q3UnfyqZxSe8KIe47JTEIHD9dT2tNOVcXxqLLiUDuBsoaUc1hN7XmL64pepFIPEi6v4eihDqr6W2hkOi/5L+AU9jKFZvYymfW1V9LUPyuzJsDBRhXDsPlnPDNuriOdsjmiMcYlpjhSoECEcL9ymfyRVsrowGnEuqTC+XEPtzfqZs+d2dNIUCfTUjqb7qLJgOPDeK3yxti6bq9wyeFCDs1YwqbpNzK3ZR2VgYO0lZ3IS8UXcPPkntieIww+wz1NtlRe7Ci6kKuwvLOoU2nHXUVZ7ofKvmZKQ8co6u/n/tA1PNP/XgrDgpPz0nnuHyxNntbjUHtvdL3ytNKgDDaqGIbNP+NO+QyMenI9YdSYWJjiSIE5deVUlRYy72gblSfMBnFSb7T3BGN6o2n1PD3ZczsLqijrDzCtYwsHK8+iu7CGysBBtra38eXHNzO15y0uaVlBT0ElPQXldHccYYE4o5I/zLotesrGY92w+JyB18hgz7Wp9GTW1y7nQo9pbdsJVztZYlNxUO98CukPU9fzDqGCYnr8NVDQwpf8v6IxPJ2mklOiVbuDIeori2MawLJiX3Sk0Rs6HiJdX1mcehqUwUYV7/1s9m3+Q5mhXPk2vNlMV8C9R9ccuLJ5M6cX7OWWydvh4BYItEFJNUw9K+Y8qQYR5EzB2CTDcYUpjjQ4LCfQd+QIrer8IIPhfp7b+hY9vkpeTXeBI0/23MLOFsIFzoJJtT17CBcU0lE8jd5OZ3LfeztexFdWQ6G/ikLgQFs5vf4U0non6LlufKeFllX/wdoTbo6pOlRDsWHXEbr7Qu4s7xp+V/RhAMrEx4UViUONE9K2j3o9QnfYT0j9QD+dlFFb0MaHCl7mycr50arh/iLuvGJedCQBxIQ1//LVxgHLzqbEYKOKbNv8U1HmrnxdgXDULBcxB55duJ95ex4FrYLWPUABdB+DwvKknYLJnbuiI9P4+T05i1KzSYbjClMcafB84flcG/glNYXQLWWUazd1/jCbpl+b/hoVnuy5JYefo7cPAvip6D9KWOrYULuUEr+zaFFl4CAdRVNiDk8prXeCnnVLqIQZviMDGouhGoruvlBMgsKq0kL2tnSxvzVMt9srjpiLBlVC1TM5sWQT1E4GN136O4daKCicwoWlPfTMPa4IspYheKhIomza/FMxQw2SVHFuyzp6Ciqh8yD4S510/cFe6DgAU84cYM6KzPfp9VcSppDTOv4AK38PJy2Ghbdk5x6Hc8/GmMIURwISTVYDCJWfwi/4cy7lJUq7D7C3fzJrey7mpW2F9AYPsnFPCyV+H2c2VEfPM1gDGsmeO+X0iz3ZcxuYesX/x8KpZ9Dk9rQ7iqc5PU53ASF/gRDsamW/ry5GzgEmsQQ969L+LjrKprHhzeZogx+5R2+jn+wZlBX7QJ3Eg12e4wHqK0qSL6wUMVUc3AI9LdAfdjIHhwL4NUB78aykacqTmQAjijVtRjOSKBXnuytfz6r/oLbvUIw58NyDP2N/Qa2j7Iornfr+Ymc7gRM/Mt/H1x9keudWuqUQSmvg4GvwwgNM7VkCNWeSVWyS4bjzOZniSED8FxmzyFDHTPYUzWF/fw/lxT5mVZZT0NNJWaGfadVO4r2GmjI27DrCxj0tA16W6IuSRvbcXbWLWXDgMcAZaZxYHuTimYXwvltZOlhjl6BnXdrfwc7aa+l+Jzwgxbm30Y9/BolyVUUc00Oai7ymiqlnQKgP9m+E/iBUzeBg0RzKxMeu2sUJD0/2w/KasNJmtCKJUnW+Tz2DtSfcPOCZdxRPo7Szxfk+IwuEhQLHR05x54mMVhu6XiVUUExQ/eAvgUAHlEzirAPPxeRTywo2yTAv1t3JJKY4UsDbcC35599HPwdDyluHO2nvDVJQEJsJpdvNJpvUJJRGr/doxdyYCKr9vjp4Xwp5qBJcY33tcvwVc4GmwY9Nkdmht3nv3tUxc1AiYcQRhbvk8I8pDQfp8fUCvZQXn8iFp0xyzC01s+kJl7Cu6IJomvIIQ4WzZjwMNhfUz4PnvgP9ISifDBXT2drUxa+rFkVHmBG2NrYNeH921S5mXvsjUDENmnc4SkP7ofakGHNb5Nns76+jtKuVgkA7XQVlFPkLjiuakipqgu9k/56TmQZnnZ/RiD8jd5jiSJPeYJhp1aU0dwSiCxv5+pxw3VSY2vMWrFuT0o/F2zA2MsMJz610ypdOTT68/dGq1TQ0/Zaa4CE3nfvF/KblfaBOJFL5sSb2t/bQ3BGg0C/Mqi1P8yk4zA69zbWBX1IcOiFmDsqhiuuAc6K9rBlHjtBRMoVCbzTa5JOgsBiu+VcuAi4axvWTjUTuWbMz4Whk1M0CTdvgjTVwwulOqpquo9Dbzqbi6/FPP5OGuOpb97cNUIyNzCA48yYWTt4Owe7jUVWTT455l46niLnVabSPVIOq41cKB6C6Ad78LbOC7Rz70w/YUnkxTaUnR6+TUeWbqJM063znWZjDfExiiiNNSvxOOGhkYSOAcL8OWBc7EZM7dzGvZQXUn5TSj2VYjVzTNs7a9yi+shoCJdOpD3fywc5V7A5fzh7/HOorSuj2LNDUFQjT3hMcco2NRL37hb3P0+uvivpeIv/P6ngOuDZaL95HA2TVVDFSs0DW7NFeJ/Fkt5HuaWX67rdpTVD9zBnVfO9if4Iw1g8DH07tmpFGe+OPYPc6KK2FmpPgsLPa8tR3XcBUfwlLe9dmt9GONw2u+7Y5zMcwpjjS5MyGahpqyli7vSnqI3jrcGdKx0YjYtL8saTVkO18ip6CSgr9VZT1HaW2Zw+lwRZuD+zlh2W3cer88wGi8u867OS26u4LsXZ7U9LIqEQN5pYHOzmok2NncGsx06Q5pl68j6Y03OGsoRGJYkonxj8H8wGyZo9O4iROZi6a2vOWk1BypL3yqWfAVf98/Nm9/r9QXAXTznYCFCLkstE2h/mYxhTHMPHmoupXpTcU5mBbDyV+H43HuunqCw0Y7lcGDjoRMV5S+LGk1ZC17aOnoJzqvqNM69hCqKCYQ31lVPS3cEXnKta8Bnv8czjU3suhjl56g05k1OzQ21wQfJHpcoTSjtls6DkfGLx3fdb8Mzkr4vTsPBy7rnrTtmi9QX00acT4/2jVas7a92h0ImTpO7so3fwPbJl5E7dcl2IP3CUmH9n+tuhzKPH76A2Fo6sRxqfDH9FoJImT+FjhlITVz+p4zslCnEpHIxWFGun1Rxpt13To3HiOG+0J5jAfk/64QTDFMUwiuag27DrClKoSuvpCLJx9XCnUlhcxe3J5zMuyv7+OWn/cyzPMH8uGN5tp7ggMsOVf01ZCaX8rtT17CBUUEy4opkA7OcYkuqjgUl5iVempVJUWOiG1fSFuOLGDBQd+TW9lJQHfLIpD7VzSvgKaTuOeV33JG8p3u07P7qNwaLvTEBX4HMdtXKjn0Yq50cmKjce6j0eDpRHj39D025iJkFCML+Snoem3pGy6cfEq4x0H2plW7UzAbO8JUl7kjz6fwY7zktJoJImTeEvlkoQ/xJrgITY0Th00oeQ9a3biO7wjJrNA6Tu7qNn2j5x5/dcTjyDyodGe53l32vc7/32FcNFXcidDDsm2by3X4b6mONIkvufQ3NlLeZE/JuVFBG9iw3vW7GRLj5Pn6bktu5wfeH8XNQXd/HGaswDU1sa2mDQaJYU+zpxRzdb9A6NrugPhhFFbf+i5gEs6V1AabKHHNwl/f4Bi+tjGiXRSxinhwwPuKRLr7/VVBAtCsPMpmjuuiF7DO/ejqy9Ec0ctU3uWcP2RB5hdHIayGqg7zTF/9LTy3qPP86tjJ8dca+v+NtDjobTL9r/GMX895SXB4737JL3fmuAhAiXTY8oCvgpqug8MqDsY96zZGZNl1xsoMKm0KK1zpUWSSLrwq76ESwoHKqaj3a2DJpRs7ghwXejFAQr1WCfJTU/5kEp96hlw6hUDIsx4Yw3UnWJ+jjTJdbivKY40SWV+QyIia5DvrCphbss6agMH6Sibxq/9H4xmf/X2fIHonJCNe1pSlq+p9GS2zLyJmre+S0XfMTp9k9gmp9FCFZV0sb+/Libde4nfl3Bmek9B+YDGuzsQprU7SDDcT1+4nx0H2tlBPfO7a2isPIcLZ3vOUVLFRfX7ueia2KzA8c8r3NFAfaid5oDHbJKk93uscAr14c4YJ3txuJPmOFPPUGaByHogER9VJEIuEHQCHmaH3ubSjpeo7muicu9J0TT3GSHB/JE7r0hSt8nPxv/+hyETSqb6/cXIkA+p1Jt3wuyLYkc+Pa3mIB8DmOJIk/ghYaTnOtjysF68JhsgYU8znmPdfdHJdhH2t/bgL0gcyXXLdR+GppOivoO33upmYVGAkhBsmn4VSyqcCXuRxrWjY2DUU3QBpVg/N8FwfzQMOdLwtvRMoaKnHfA0XimaPiKO89JwCHTwdPBbKi/mg52rAGIa0i2VS1nqqTeSofns0NssCz7JMS3jQLiWkzpbmNf+COtrl1N/wunpp5YZCVPPSCmhZKKotSEXwMqHVOrmIB+zmOJIk/ghYWR9iGwsDzsiPL3KmtA7BMoajjc6Huori9nQcz6XtK8gWBCKMaEx7ypoDiW5wHGeLzyf0/r/x+ktpmn6iDjOp+xfG7MMb6JGran0ZDZVHXeyRxrSpmDidUuGYmbfbt4TeJ7i/gM0B+tZ27+IC4IvU1VbR7G/irKeIBfPnwo9rSwsfQMWXzuy2erDoKn0ZP4wREqQ+Ki14nAn/v4O5/tLeNI8yVSbD74WY1iY4hgD1JQVDUjrET8C8Wav9TZu9ZXXwLuucXrKQQbMzHZ65/Og6TRPYzLX05gM3VDu8c9hfc1yp3EdhunjaMVcXjthBkvjzFrx1FcW81rH8YmQAASHF5nyLtnL5e2r6JAKmnQyZeEObuh/kvpQgOa+d0HQXUURhp86P0fER611FE9jQ+1SFiZ6/vmUqTYffC3GsDDFMUIiYbne5WEh+w1JWbET7RS5ZryTPpJK27e/kUXnnA3vvmpAhFQkKstRIIlNF96GsqsvFJ34WOgriKnXVHoyLL52wPGZJJPRIR+p2ExxyQmU+6twVHINwc4WLqhqgillSXvBMTJ4e+4lM6EpnPHGdyhFFdnvzSzg3T+AfMpUmy++lnFA/Hsytectzup4jj+jGdadmfFRpSmOERIJyx10eViGbgAaj3UT7lcOtvVE95UUOnNCEmWBvfCU+phrep3O3lTazf76aK/S17GEhukDzR6DRV7Er3X+9NaDlLuzzCPmubJij3xDmEHyocdeX1mMb3+j82yCx02M5aVVvNF6iLbW3cdDW/u7KO3vcOaKeE+So577UMoybWWab36FfPC1jAMGdGheWOvMASqZn5V30xRHjhhJbzldu3pMeG0wGO1VjjQTauQeksWLp9KY5kMK6TuvmAclZw+0r/e08szuUzg0Y0lM5NvO2mvZEe9HGcWeezYmIZpfYRyRg3fTFEeajEaPOd1rJgrPzFQm1EEbprGUf2iwyXhxkW/AwGVxR7Hnno1JiOZXGEfk4N00xZEmo9FjTveayZIKJkttkTHyzQwyGEns603PhQZkqU3IWO25m19h/JODd9MUxzjBOyrZ4D8eXlteWuW8RIOktsgYY60xTWhfT9EsOJZ77uZXGN/k4N00xTFOiB2VnBMbXls6ZdDUFhkzs43lxjRdrOdu5Cs5eDdNcYxX0kltkclrjvHGNC1/kvXcjXwly++mKQ4js4zxxjQfor6GIh9Cmo2JjSkOwxhjjAXlZoxvCoauYhiGYRjHMcVhGIZhpIUpDsMwDCMtTHEYhmEYaWGKwzAMw0gLUdXRliHriEgzkG6ipjrgSBbEyTQmZ2YZC3KOBRnB5Mw0uZbzRFVNuKzphFAcw0FENqnqgtGWYyhMzswyFuQcCzKCyZlp8klOM1UZhmEYaWGKwzAMw0gLUxzJeWi0BUgRkzOzjAU5x4KMYHJmmryR03wchmEYRlrYiMMwDMNIC1MchmEYRlpMSMUhIpeLyJ9E5E0RuTPB/mIRWenuf0lEZnv2fdUt/5OILM1HOUVktoj0iMhm9++HoyjjxSLyioiEROS6uH0fF5Fd7t/HsyVjBuQMe57l6lGW869FZIeIbBGR34rIiZ59+fQ8B5Mzn57n7SKy1ZVlg4ic7tmXT7/1hHLm8rceg6pOqD/AB7wFnAQUAa8Bp8fV+SzwQ/fzcmCl+/l0t34xMMc9jy8P5ZwNbMuTZzkbOAv4CXCdp7wW2O3+r3E/1+SbnO6+zjx6NxcDZe7nz3i+83x7ngnlzMPnWeX5/GHg1+7nfPutJ5MzJ7/1+L+JOOJYBLypqrtVtQ9YAVwdV+dq4L/cz6uAPxMRcctXqGpAVd8G3nTPl29y5oohZVTVPaq6BeiPO3Yp8BtVbVHVY8BvgMvzUM5ckoqc61Q1sorTH4AG93O+Pc9kcuaSVORs92yWA5Foobz6rQ8i56gwERXHDGCfZ7vRLUtYR1VDQBswOcVj80FOgDki8qqI/F5ELhpFGbNxdsRvnwAAA0VJREFUbLqM9FolIrJJRP4gItdkVrQY0pXzFmDNMI8dCSORE/LseYrI50TkLeA7wBfSOTYP5ITc/NZjsBUAxycHgVmqelREzgN+JSLz43otRuqcqKr7ReQk4HcislVV3xpNgUTkRmAB8P7RlGMoksiZV89TVX8A/EBEPgp8Hciqf2i4JJFzVH7rE3HEsR+Y6dlucMsS1hERP1ANHE3x2FGX0x1eHwVQ1Zdx7KenjpKM2Tg2XUZ0LVXd7/7fDawH3p1J4TykJKeIfAD4O+DDqhpI59g8kDPvnqeHFUBkBJR3z9NDVM4c/tZjybVTZbT/cEZZu3EcXhFH1Py4Op8j1un8uPt5PrEOs91kz2E2EjnrI3LhONz2A7WjIaOn7iMMdI6/jePIrXE/Z1zGDMhZAxS7n+uAXcQ5LnP8nb8bp3GYG1eeV89zEDnz7XnO9Xy+Ctjkfs6333oyOXPyWx8gc7YvkI9/wJXAG+6L/Xdu2TdwekYAJcDPcRxifwRO8hz7d+5xfwKuyEc5gWuB7cBm4BXgqlGUcSGOzbYLZ9S23XPsJ13Z3wRuHuVnmVBO4H3AVvfHvBW4ZZTlfBY45H63m4HVefo8E8qZh8/zfs9vZR2eBjvPfusJ5czlb937ZylHDMMwjLSYiD4OwzAMYwSY4jAMwzDSwhSHYRiGkRamOAzDMIy0MMVhGIZhpIXNHDeMDCMiYZxQ0wgrVPWe0ZLHMDKNheMaRoYRkU5VrRiijk9Vw8m2Uz3OMEYDM1UZRo4QkT0icq+IvAJcn2D7BnfNhW0icq/nuE4R+Z6IvAacP2o3YBgupjgMI/OUehbW2Swiyzz7jqrquaq6wrsNPAfcC1wKnAMs9GSOLQdeUtWzVXVDzu7CMJJgPg7DyDw9qnpOkn0rk2wvBNarajOAiPwUuBj4FRAGfpENQQ1jONiIwzByS9cQ24noNb+GkU+Y4jCM/OCPwPtFpE5EfMANwO9HWSbDSIiZqgwj85SKyGbP9q9V9c7BDlDVgyJyJ07mUwH+V1WfzKaQhjFcLBzXMAzDSAszVRmGYRhpYYrDMAzDSAtTHIZhGEZamOIwDMMw0sIUh2EYhpEWpjgMwzCMtDDFYRiGYaTF/w+l9y70f6pGHwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZDmir9PTv6HG",
        "colab_type": "text"
      },
      "source": [
        "We finally need to save the model for future. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F1HPv85MwAVH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "stochastic_model.save_weights('stochastic_model_weights.tf')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8imDcLMgwFaD",
        "colab_type": "text"
      },
      "source": [
        "And some code to show you how to load this model and use it."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "V73yVPKAwJLK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def return_loaded_model():\n",
        "  loaded_model = StochasticModel()\n",
        "  loaded_model.compile(loss=mse_loss, optimizer=adam, metrics=[mse_loss])\n",
        "  loaded_model.predict(np.zeros((1, 6)))\n",
        "  loaded_model.load_weights('stochastic_model_weights.tf')\n",
        "  return loaded_model"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}