examples/trustscore_mnist.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trust Scores applied to MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to know when a machine learning classifier's predictions can be trusted. Relying on the classifier's (uncalibrated) prediction probabilities is not optimal and can be improved upon. *Trust scores* measure the agreement between the classifier and a modified nearest neighbor classifier on the test set. The trust score is the ratio between the distance of the test instance to the nearest class different from the predicted class and the distance to the predicted class. Higher scores correspond to more trustworthy predictions. A score of 1 would mean that the distance to the predicted class is the same as to another class.\n",
    "\n",
    "The original paper on which the algorithm is based is called [To Trust Or Not To Trust A Classifier](https://arxiv.org/abs/1805.11783). Our implementation borrows heavily from https://github.com/google/TrustScore, as does the example notebook.\n",
    "\n",
    "Trust scores work best for low to medium dimensional feature spaces. This notebook illustrates how you can **apply trust scores to high dimensional** data like images by adding an additional pre-processing step in the form of an [auto-encoder](https://en.wikipedia.org/wiki/Autoencoder) to reduce the dimensionality. Other dimension reduction techniques like PCA can be used as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "tf.logging.set_verbosity(tf.logging.ERROR)  # suppress deprecation messages\n",
    "from tensorflow.keras import backend as K\n",
    "from tensorflow.keras.layers import Conv2D, Dense, Dropout, Flatten, MaxPooling2D, Input, UpSampling2D\n",
    "from tensorflow.keras.models import Model, load_model\n",
    "from tensorflow.keras.utils import to_categorical\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.model_selection import StratifiedShuffleSplit\n",
    "from alibi.confidence import TrustScore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train shape: (60000, 28, 28) y_train shape: (60000,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAM3ElEQVR4nO3dXahc9bnH8d/vpCmI6UXiS9ik0bTBC8tBEo1BSCxbQktOvIjFIM1FyYHi7kWUFkuo2It4WaQv1JvALkrTkmMJpGoQscmJxVDU4o5Es2NIjCGaxLxYIjQRJMY+vdjLso0za8ZZa2ZN8nw/sJmZ9cya9bDMz7VmvczfESEAV77/aroBAINB2IEkCDuQBGEHkiDsQBJfGeTCbHPoH+iziHCr6ZW27LZX2j5o+7Dth6t8FoD+cq/n2W3PkHRI0nckHZf0mqS1EfFWyTxs2YE+68eWfamkwxFxJCIuSPqTpNUVPg9AH1UJ+zxJx6a9Pl5M+xzbY7YnbE9UWBaAivp+gC4ixiWNS+zGA02qsmU/IWn+tNdfL6YBGEJVwv6apJtsf8P2VyV9X9L2etoCULeed+Mj4qLtByT9RdIMSU9GxP7aOgNQq55PvfW0ML6zA33Xl4tqAFw+CDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJ9Dw+uyTZPirpnKRPJV2MiCV1NAWgfpXCXrgrIv5Rw+cA6CN244EkqoY9JO2wvcf2WKs32B6zPWF7ouKyAFTgiOh9ZnteRJywfb2knZIejIjdJe/vfWEAuhIRbjW90pY9Ik4Uj2ckPS1paZXPA9A/PYfd9tW2v/bZc0nflTRZV2MA6lXlaPxcSU/b/uxz/i8iXqilKwC1q/Sd/UsvjO/sQN/15Ts7gMsHYQeSIOxAEoQdSIKwA0nUcSNMCmvWrGlbu//++0vnff/990vrH3/8cWl9y5YtpfVTp061rR0+fLh0XuTBlh1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkuCuty4dOXKkbW3BggWDa6SFc+fOta3t379/gJ0Ml+PHj7etPfbYY6XzTkxcvr+ixl1vQHKEHUiCsANJEHYgCcIOJEHYgSQIO5AE97N3qeye9VtuuaV03gMHDpTWb7755tL6rbfeWlofHR1tW7vjjjtK5z127Fhpff78+aX1Ki5evFha/+CDD0rrIyMjPS/7vffeK61fzufZ22HLDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJcD/7FWD27Nlta4sWLSqdd8+ePaX122+/vaeeutHp9/IPHTpUWu90/cKcOXPa1tavX18676ZNm0rrw6zn+9ltP2n7jO3JadPm2N5p++3isf2/NgBDoZvd+N9LWnnJtIcl7YqImyTtKl4DGGIdwx4RuyWdvWTyakmbi+ebJd1Tc18AatbrtfFzI+Jk8fyUpLnt3mh7TNJYj8sBUJPKN8JERJQdeIuIcUnjEgfogCb1eurttO0RSSoez9TXEoB+6DXs2yWtK56vk/RsPe0A6JeO59ltPyVpVNK1kk5L2ijpGUlbJd0g6V1J90XEpQfxWn0Wu/Ho2r333lta37p1a2l9cnKybe2uu+4qnffs2Y7/nIdWu/PsHb+zR8TaNqUVlToCMFBcLgskQdiBJAg7kARhB5Ig7EAS3OKKxlx//fWl9X379lWaf82aNW1r27ZtK533csaQzUByhB1IgrADSRB2IAnCDiRB2IEkCDuQBEM2ozGdfs75uuuuK61/+OGHpfWDBw9+6Z6uZGzZgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJ7mdHXy1btqxt7cUXXyydd+bMmaX10dHR0vru3btL61cq7mcHkiPsQBKEHUiCsANJEHYgCcIOJEHYgSS4nx19tWrVqra1TufRd+3aVVp/5ZVXeuopq45bdttP2j5je3LatEdtn7C9t/hr/18UwFDoZjf+95JWtpj+m4hYVPw9X29bAOrWMewRsVvS2QH0AqCPqhyge8D2m8Vu/ux2b7I9ZnvC9kSFZQGoqNewb5K0UNIiSScl/ardGyNiPCKWRMSSHpcFoAY9hT0iTkfEpxHxL0m/k7S03rYA1K2nsNsemfbye5Im270XwHDoeJ7d9lOSRiVda/u4pI2SRm0vkhSSjkr6UR97xBC76qqrSusrV7Y6kTPlwoULpfNu3LixtP7JJ5+U1vF5HcMeEWtbTH6iD70A6CMulwWSIOxAEoQdSIKwA0kQdiAJbnFFJRs2bCitL168uG3thRdeKJ335Zdf7qkntMaWHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSYMhmlLr77rtL688880xp/aOPPmpbK7v9VZJeffXV0jpaY8hmIDnCDiRB2IEkCDuQBGEHkiDsQBKEHUiC+9mTu+aaa0rrjz/+eGl9xowZpfXnn28/5ifn0QeLLTuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJMH97Fe4TufBO53rvu2220rr77zzTmm97J71TvOiNz3fz257vu2/2n7L9n7bPy6mz7G90/bbxePsupsGUJ9uduMvSvppRHxL0h2S1tv+lqSHJe2KiJsk7SpeAxhSHcMeEScj4vXi+TlJByTNk7Ra0ubibZsl3dOvJgFU96Wujbe9QNJiSX+XNDciThalU5LmtplnTNJY7y0CqEPXR+Ntz5K0TdJPIuKf02sxdZSv5cG3iBiPiCURsaRSpwAq6SrstmdqKuhbIuLPxeTTtkeK+oikM/1pEUAdOu7G27akJyQdiIhfTyttl7RO0i+Kx2f70iEqWbhwYWm906m1Th566KHSOqfXhkc339mXSfqBpH229xbTHtFUyLfa/qGkdyXd158WAdShY9gj4m+SWp6kl7Si3nYA9AuXywJJEHYgCcIOJEHYgSQIO5AEPyV9Bbjxxhvb1nbs2FHpszds2FBaf+655yp9PgaHLTuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJMF59ivA2Fj7X/264YYbKn32Sy+9VFof5E+Roxq27EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBOfZLwPLly8vrT/44IMD6gSXM7bsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5BEN+Ozz5f0B0lzJYWk8Yj4re1HJd0v6YPirY9ExPP9ajSzO++8s7Q+a9asnj+70/jp58+f7/mzMVy6uajmoqSfRsTrtr8maY/tnUXtNxHxy/61B6Au3YzPflLSyeL5OdsHJM3rd2MA6vWlvrPbXiBpsaS/F5MesP2m7Sdtz24zz5jtCdsTlToFUEnXYbc9S9I2ST+JiH9K2iRpoaRFmtry/6rVfBExHhFLImJJDf0C6FFXYbc9U1NB3xIRf5akiDgdEZ9GxL8k/U7S0v61CaCqjmG3bUlPSDoQEb+eNn1k2tu+J2my/vYA1KWbo/HLJP1A0j7be4tpj0haa3uRpk7HHZX0o750iEreeOON0vqKFStK62fPnq2zHTSom6Pxf5PkFiXOqQOXEa6gA5Ig7EAShB1IgrADSRB2IAnCDiThQQ65a5vxfYE+i4hWp8rZsgNZEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoMesvkfkt6d9vraYtowGtbehrUvid56VWdvN7YrDPSimi8s3J4Y1t+mG9behrUvid56Naje2I0HkiDsQBJNh3284eWXGdbehrUvid56NZDeGv3ODmBwmt6yAxgQwg4k0UjYba+0fdD2YdsPN9FDO7aP2t5ne2/T49MVY+idsT05bdoc2zttv108thxjr6HeHrV9olh3e22vaqi3+bb/avst2/tt/7iY3ui6K+lrIOtt4N/Zbc+QdEjSdyQdl/SapLUR8dZAG2nD9lFJSyKi8QswbH9b0nlJf4iI/y6mPSbpbET8ovgf5eyI+NmQ9PaopPNND+NdjFY0Mn2YcUn3SPpfNbjuSvq6TwNYb01s2ZdKOhwRRyLigqQ/SVrdQB9DLyJ2S7p0SJbVkjYXzzdr6h/LwLXpbShExMmIeL14fk7SZ8OMN7ruSvoaiCbCPk/SsWmvj2u4xnsPSTts77E91nQzLcyNiJPF81OS5jbZTAsdh/EepEuGGR+addfL8OdVcYDui5ZHxK2S/kfS+mJ3dSjF1HewYTp32tUw3oPSYpjx/2hy3fU6/HlVTYT9hKT5015/vZg2FCLiRPF4RtLTGr6hqE9/NoJu8Xim4X7+Y5iG8W41zLiGYN01Ofx5E2F/TdJNtr9h+6uSvi9pewN9fIHtq4sDJ7J9taTvaviGot4uaV3xfJ2kZxvs5XOGZRjvdsOMq+F11/jw5xEx8D9JqzR1RP4dST9vooc2fX1T0hvF3/6me5P0lKZ26z7R1LGNH0q6RtIuSW9L+n9Jc4aotz9K2ifpTU0Fa6Sh3pZrahf9TUl7i79VTa+7kr4Gst64XBZIggN0QBKEHUiCsANJEHYgCcIOJEHYgSQIO5DEvwEvYRv57rmVLgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "(x_train, y_train), (x_test, y_test) = tf.keras.datasets.mnist.load_data()\n",
    "print('x_train shape:', x_train.shape, 'y_train shape:', y_train.shape)\n",
    "plt.gray()\n",
    "plt.imshow(x_test[0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare data: scale, reshape and categorize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train shape: (60000, 28, 28, 1) x_test shape: (10000, 28, 28, 1)\n",
      "y_train shape: (60000, 10) y_test shape: (10000, 10)\n"
     ]
    }
   ],
   "source": [
    "x_train = x_train.astype('float32') / 255\n",
    "x_test = x_test.astype('float32') / 255\n",
    "x_train = np.reshape(x_train, x_train.shape + (1,))\n",
    "x_test = np.reshape(x_test, x_test.shape + (1,))\n",
    "print('x_train shape:', x_train.shape, 'x_test shape:', x_test.shape)\n",
    "y_train = to_categorical(y_train)\n",
    "y_test = to_categorical(y_test)\n",
    "print('y_train shape:', y_train.shape, 'y_test shape:', y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, xmax = -.5, .5\n",
    "x_train = ((x_train - x_train.min()) / (x_train.max() - x_train.min())) * (xmax - xmin) + xmin\n",
    "x_test = ((x_test - x_test.min()) / (x_test.max() - x_test.min())) * (xmax - xmin) + xmin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define and train model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example we are not interested in optimizing model performance so a simple softmax classifier will do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sc_model():\n",
    "    x_in = Input(shape=(28, 28, 1))\n",
    "    x = Flatten()(x_in)\n",
    "    x_out = Dense(10, activation='softmax')(x)\n",
    "    sc = Model(inputs=x_in, outputs=x_out)\n",
    "    sc.compile(loss='categorical_crossentropy', optimizer='sgd', metrics=['accuracy'])\n",
    "    return sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"model\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_1 (InputLayer)         [(None, 28, 28, 1)]       0         \n",
      "_________________________________________________________________\n",
      "flatten (Flatten)            (None, 784)               0         \n",
      "_________________________________________________________________\n",
      "dense (Dense)                (None, 10)                7850      \n",
      "=================================================================\n",
      "Total params: 7,850\n",
      "Trainable params: 7,850\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "Epoch 1/5\n",
      "60000/60000 [==============================] - 1s 12us/sample - loss: 1.2706 - acc: 0.6963\n",
      "Epoch 2/5\n",
      "60000/60000 [==============================] - 1s 9us/sample - loss: 0.7030 - acc: 0.8422\n",
      "Epoch 3/5\n",
      "60000/60000 [==============================] - 1s 9us/sample - loss: 0.5762 - acc: 0.8618\n",
      "Epoch 4/5\n",
      "60000/60000 [==============================] - 1s 9us/sample - loss: 0.5155 - acc: 0.8706\n",
      "Epoch 5/5\n",
      "60000/60000 [==============================] - 1s 9us/sample - loss: 0.4787 - acc: 0.8759\n"
     ]
    }
   ],
   "source": [
    "sc = sc_model()\n",
    "sc.summary()\n",
    "sc.fit(x_train, y_train, batch_size=128, epochs=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate the model on the test set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy:  0.8862\n"
     ]
    }
   ],
   "source": [
    "score = sc.evaluate(x_test, y_test, verbose=0)\n",
    "print('Test accuracy: ', score[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define and train auto-encoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ae_model():\n",
    "    # encoder\n",
    "    x_in = Input(shape=(28, 28, 1))\n",
    "    x = Conv2D(16, (3, 3), activation='relu', padding='same')(x_in)\n",
    "    x = MaxPooling2D((2, 2), padding='same')(x)\n",
    "    x = Conv2D(8, (3, 3), activation='relu', padding='same')(x)\n",
    "    x = MaxPooling2D((2, 2), padding='same')(x)\n",
    "    x = Conv2D(4, (3, 3), activation=None, padding='same')(x)\n",
    "    encoded = MaxPooling2D((2, 2), padding='same')(x)\n",
    "    encoder = Model(x_in, encoded)\n",
    "\n",
    "    # decoder\n",
    "    dec_in = Input(shape=(4, 4, 4))\n",
    "    x = Conv2D(4, (3, 3), activation='relu', padding='same')(dec_in)\n",
    "    x = UpSampling2D((2, 2))(x)\n",
    "    x = Conv2D(8, (3, 3), activation='relu', padding='same')(x)\n",
    "    x = UpSampling2D((2, 2))(x)\n",
    "    x = Conv2D(16, (3, 3), activation='relu')(x)\n",
    "    x = UpSampling2D((2, 2))(x)\n",
    "    decoded = Conv2D(1, (3, 3), activation=None, padding='same')(x)\n",
    "    decoder = Model(dec_in, decoded)\n",
    "    \n",
    "    # autoencoder = encoder + decoder\n",
    "    x_out = decoder(encoder(x_in))\n",
    "    autoencoder = Model(x_in, x_out)\n",
    "    autoencoder.compile(optimizer='adam', loss='mse')\n",
    "    \n",
    "    return autoencoder, encoder, decoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"model_3\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_2 (InputLayer)         [(None, 28, 28, 1)]       0         \n",
      "_________________________________________________________________\n",
      "model_1 (Model)              (None, 4, 4, 4)           1612      \n",
      "_________________________________________________________________\n",
      "model_2 (Model)              (None, 28, 28, 1)         1757      \n",
      "=================================================================\n",
      "Total params: 3,369\n",
      "Trainable params: 3,369\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "Train on 60000 samples, validate on 10000 samples\n",
      "Epoch 1/8\n",
      "60000/60000 [==============================] - 29s 477us/sample - loss: 0.0606 - val_loss: 0.0399\n",
      "Epoch 2/8\n",
      "60000/60000 [==============================] - 34s 572us/sample - loss: 0.0341 - val_loss: 0.0301\n",
      "Epoch 3/8\n",
      "60000/60000 [==============================] - 43s 715us/sample - loss: 0.0288 - val_loss: 0.0272\n",
      "Epoch 4/8\n",
      "60000/60000 [==============================] - 48s 806us/sample - loss: 0.0265 - val_loss: 0.0253\n",
      "Epoch 5/8\n",
      "60000/60000 [==============================] - 41s 680us/sample - loss: 0.0249 - val_loss: 0.0239\n",
      "Epoch 6/8\n",
      "60000/60000 [==============================] - 39s 649us/sample - loss: 0.0237 - val_loss: 0.0230\n",
      "Epoch 7/8\n",
      "60000/60000 [==============================] - 33s 545us/sample - loss: 0.0229 - val_loss: 0.0222\n",
      "Epoch 8/8\n",
      "60000/60000 [==============================] - 29s 484us/sample - loss: 0.0224 - val_loss: 0.0217\n"
     ]
    }
   ],
   "source": [
    "ae, enc, dec = ae_model()\n",
    "ae.summary()\n",
    "ae.fit(x_train, x_train, batch_size=128, epochs=8, validation_data=(x_test, x_test))\n",
    "ae.save('mnist_ae.h5')\n",
    "enc.save('mnist_enc.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ae = load_model('mnist_ae.h5')\n",
    "enc = load_model('mnist_enc.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Trust Scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize trust scores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = TrustScore()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key is to **fit and calculate the trust scores on the encoded instances**. The encoded data still needs to be reshaped from (60000, 4, 4, 4) to (60000, 64) to comply with the k-d tree format. This is handled internally:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Reshaping data from (60000, 4, 4, 4) to (60000, 64) so k-d trees can be built.\n"
     ]
    }
   ],
   "source": [
    "x_train_enc = enc.predict(x_train)\n",
    "ts.fit(x_train_enc, y_train, classes=10)  # 10 classes present in MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the trust scores and closest not predicted classes of the predictions on the test set, using the distance to the 5th nearest neighbor in each class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Reshaping data from (1000, 4, 4, 4) to (1000, 64) so k-d trees can be queried.\n"
     ]
    }
   ],
   "source": [
    "n_samples = 1000 # calculate the trust scores for the first 1000 predictions on the test set\n",
    "x_test_enc = enc.predict(x_test[:n_samples])\n",
    "y_pred = sc.predict(x_test[:n_samples])\n",
    "score, closest_class = ts.score(x_test_enc[:n_samples], y_pred, k=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's inspect which predictions have low and high trust scores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5\n",
    "\n",
    "# lowest and highest trust scores\n",
    "idx_min, idx_max = np.argsort(score)[:n], np.argsort(score)[-n:]\n",
    "score_min, score_max = score[idx_min], score[idx_max]\n",
    "closest_min, closest_max = closest_class[idx_min], closest_class[idx_max]\n",
    "pred_min, pred_max = y_pred[idx_min], y_pred[idx_max]\n",
    "imgs_min, imgs_max = x_test[idx_min], x_test[idx_max]\n",
    "label_min, label_max = np.argmax(y_test[idx_min], axis=1), np.argmax(y_test[idx_max], axis=1)\n",
    "\n",
    "# model confidence percentiles\n",
    "max_proba = y_pred.max(axis=1)\n",
    "\n",
    "# low score high confidence examples\n",
    "idx_low = np.where((max_proba>0.80) & (max_proba<0.9) & (score<1))[0][:n]\n",
    "score_low = score[idx_low]\n",
    "closest_low = closest_class[idx_low]\n",
    "pred_low = y_pred[idx_low]\n",
    "imgs_low = x_test[idx_low]\n",
    "label_low = np.argmax(y_test[idx_low], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Low Trust Scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image below makes clear that the low trust scores correspond to misclassified images. Because the trust scores are significantly below 1, they correctly identified that the images belong to another class than the predicted class, and identified that class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAETCAYAAABN+bF9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3debgcVZn48e8bIpuIsijbAGETBIwgy7BKoiKjLAIiKCqgI/xQQRwGBxCVBBFBAo6OoiC7IqLgoCI7JApBUDYFISwDAUTCGtaEROD8/jh1odN09e3uu3TXvd/P89ST3HrrVJ3qrrer++3qOpFSQpIkSZIkSdU1ptsdkCRJkiRJ0sBY4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIrriQJPRCwbESkiJrTRZlJE3D6E3WpbRNweEZNq/p4ZEYcMcJ3TIuL7A+7cMIqItSLi0Yh4c7f7MhQi4pcR8Z/d7sdwMDebrtPc7DHmZr9tzM0eZW6OHOZm03Wamz0mIo6PiP/pdj+Gg7nZdJ3mZo8ZyHmz3wJPRJxZJMNpDWLHFbGLOtn4KLAJcFIrC0bEPhHxfIPQrsDhg9qr1vrz1oi4LCL+ERHzIuKhiPhBi0l0DHBSSumZIe7j5yPi/oh4MSJuioit+1l+14i4PCIej4jnIuKGiNipwXIHRcSMiJgbEX8v9nuJmkWOAo7o9guKuTkglc3NWsWblYeL53rZFpr0ZG7Wtd0qIl5q9IbK3BwVKp2bEfHdiLixOPZnttG0J3MzIraJiOsi4ski72bUf5AonofUYFq0ZjFzs/rMzSE0FOfN4kN7o9z8W81i3wb2jojVB2tfOmFuDoi5OYSG4rxZLLdkRHwvXvusfW9E7F6zSMfnzVav4HkI2D0i3ljTqbHAXsCD7W60l0XEGwZrXSmlx1NKcwa4jqdSSs8NVp/a8Arwv8COwNuBfYD3AT9u1igiVgZ2Bs4Yys5FxB7Ad8nJvSFwHXBJRKzSpNk2wNXA9kWbi4H/rU3UiNiTfLL7JvAO8jH+oWJbAKSUbgPuAz45iLvUKXOzAxXPzVpnALe2smCP52Zf26WAs4GrGsTMzR5lbi5gDHAW+ThuSY/n5vPA94D3AOsCRwOTI+LzdcvNAVaonVJKL/YFzc3uMDcXMNJys69t6XmT/KG9Ni/HAc8Bv+hbIKX0OHA58LkB7cTgMDc7YG4OnaE6bxbP/xXAWsDuwNrkz9r39y0zkPNmqwWevwL3FB3osz3wIjCtdsGIGBMRX4t8xce8iLgtIj5ct8wmRQXsxYi4BfjX+g1GxLoR8bvIV1o8FhHnRsTyre5YRIwrqr17RsS1xbZmRMQHapaZUCzzoYj4U0TMB7YrYjvW9PH+iPhmRCxc0/ZtEfHrojL3QER8pkEfFrhkLiLeHBE/jIhHivXeGRF7RL5U8AzgjTXV9UlFmwUumYuIpSLirIiYXWz7yohYrya+T0Q8HxHvi3wJ3wsRMTUiVmv1sQNIKT2ZUvpRSummlNIDKaWryNXh/r5R2AO4PaX06gtxTZ92jIi7i32fGgP7tuBg4MyU0o9TSnemlA4EHqHJCSqldFBK6diU0p9SSvemlCYDN5FfIPpsAVyfUvpJSmlmSulq8gtO/TH6G+DjA+j/YDE3R1lu1qzvIGBx4IQWm/RsbtY4jXyS/2ODmLlpbvZ8bqaUDkwp/Q9wdxvNejY3i/cAP08p/S2ldH9K6afAZbz+vUBKKc2qnRqsztw0N83NxobkvFl8aK/Nya3I7xtOr1vU3DQ3zc3Ghuq8+WngrcCHU0rXFu9rr00p/bludR3lZjv34DkNqD2oPkM+SFLdcgcBXwYOBd5JvgrkVxGxAUDky+l/R65IbQwcBkypXUFErAD8Abgd2BR4P7AE8OuIaPe+Qd8mV9E2IFfKfh0RK9UtcxzwVWAd4IaI2A44B/g+sF6xr7uRq3d9zgTWLPq2M7m6PK6sExER5CtGtiE/qeuSD5r55Grgl1jwG7ApjdfEmeQXqA+TH5s5wKURsVjNMouQL7P7DLA58BbgRzV96Xsx2qesvw36vyL5m4Df97Po1sCNDeYvAhxJ3vfNgYXIx0UU69+6SMpm01eKZRcGNiJ/41DrcvKHwHa8CZhd8/e1wAYRsVmxrVWAncjPXa0/AZvWPe7dYm6OstyMiA3Jz+Ne5KvtWtHTuRn5m43lyN90NGJumps9n5sd6uncrFW89mzB698LLFZ8+Ph7RFxULFfP3DQ3652JuTmU5816+wKXppQeqpv/J2CliFijxfUMJXPT3GxFT+dmrZLz5s7AdOB/ImJWRNwR+X5P9Vd2dXbeTCk1nchP8EXAUsBc8qVEywPzgFX64jXLPwx8vW4d04CfFv/fD3gaWKIm/kly4k4o/j4KuKpuHUsVy2xa/D2JXLkr6/e4YvkjauaNIVcHjy7+nlAs85G6tn8AvlY3b2fyJVdB/slSArasia8KvAxMqpk3Ezik+P+25A9i7yjp7z7A8w3mTwO+X/x/rWK776mJvxl4BvhszXoSsHbNMp8onq8o/l4JmAHs0sLzfy45qRPwW2Cxfpa/FZjcYN/KHq/3F38vRn4BazYtXSy7Yv3jUMz/OnBXf/tUs/wXyJeqrtpg/nzgn8V2zu577GqWGV/E1mh1e4M9YW6OytwE3gjc1ffY1DxWy1Y1N8lvzh4FVmt2DGFumps9nJt1fTgEmNnisj2bmzXL/b14PF5ucKxuDuxN/mCzNXA++X3DWuamuWludve8WbN83/P84QaxJYvY+8xNcxNzczjOmzPIV6adTi4ifQSYBUypW66j8+ZYWpRSmh0R/0uu0j0NTEspPVgUxACIiCWLB2N6XfNryfdKgHzvhL+mlGpv8FR/WeFGwHui8U2g1iBXs1r16rpTSq9ExA3kamat+grgRuRq2aE188aQD4rlyfvwSm0/UkoPRMQ/mvRjQ+CRlNKdbfS9Xt92a/fpmYi4jQX3aV5K6a6av/8BLEx+0XoqpfQwuXrciv8AJpNfZL4F/Dfw/5osvxj5gK1X9nitC1yZUpoL3NtinwYsIj4CHA/skVJ6oGb+NsDXgM8DN5AT/bvkx+DrNauYW/zb9W8izc1Rl5vfA65NKV3QZh97MjcjYhHgPPKbk/ubLGdumpud6MZ5s109mZt1tiZ/s70ZcFxE3J9S+knRrz9S8/hGxHXkN98HAl+sWYe5+XrmprnZtlbPm3X2Jf+05HcNYubm65mb5uZAlZ43yc/zY8C+KaWXgZsiYhngOxHx5VRUd+gwN1su8BROJ//O83kWfEPditT/Iq8aQ34BajTk26NtbrcVLzTY/mTglw2Wfbzm/+3s01Cr7ctLJbF2Lzck5d/szgJmRMRTwDURcXR6/eWdfZ4gJ3Z/fVxA5BsdX9JPd45JKR1TbONl8mWptZYr+tpUROxG/uZ/r5TSb+vCRwPnppROLf6+LfLN3k6NiKNSSn2P7dLFv4/TG8zN0ZOb7wNWjoi9i7/73vXMiojjUkpHlLTr1dxcgfxm4oyIOKOYNyZvOl4CPpRSuhxzsz/mZvuG5LzZgV7Nzdc68dqHyNsiYjnyt9o/KVn25Yi4kfwtcC1z8/XMzcbMzcE5b/b1Z2HyVXY/rjlX1jI3X8/cbMzcHJzz5iPAP4viTp87yffIWpbXnv+OcrPdAs9V5MvjlwUurA+mlJ4tqmRbsuDd3LcC7ij+fyewT0S8MaXUd6BvVreqm8k32HogpfTPNvtYbzPyyEl9v0vclHz5cDM3A+uklBpW+CJiBvng3ZT8e8a++0Gs2GSdtwArRMQ7Sqqq88m/E2zmzmK7m5Mv6+urYr+TIb6LeKEvYRdpsswtvL5i3de20ePV91jcSL68u5mnAFJK8yPiJvJliLUvitsCTa9qiDz83FnA3imlRsfB4uRkrvUyr32I7rM+8HBKaShOAJ0wNxk1ufkB8jckfTYhvxmaQL45YZlezc2HyY9Trc8XbXYhX3oM5qa5+Xq9lpud6tXcLDOGJu8DimNmPPCXupC5aW7WMjcZ8vNmn53Jx/nrhiAvrE/+6fNt/fR5uJibmJv96NXcLFN/3pwO7BkRY1JKfffSfDv5581P1CzX2Xmz/jdb9ROv/83jm4Alm8S/BDxLvuPz28m/b3wZeFcRX4JchTqPfEOpbcnJWPubyBXJly39inyDp9XJN5c6BXhTau83kQ+Rb1i1Nvly/heBf0kL/iZy2bq225Ff6I4qHth1inV8u2aZS8gvhJuTD5SryZXmst9EjiFf6nZHsf7Vin3fuYhvUfRlW/IL2uKp7jeRxd8Xkg/SrcmJ9ptiHxdLJb+trN9PWrvPxw7kav/6xWO5fdH3P/ZzvOxIPjDH1szbp3g8/1TzeE0l3y0/mq2vyXb2IL9IfZb8LcZ3i8d/1ZplzgbOrvn7Y0U/DiJf+tg3LV2zzCTy8fuxmufoXuCCBnlxWid9H6wJc3NU5maDx7PhY1Wl3GywjtcdQ5ib5mYFcpP808ENgBPJl6tvUEwLVzE3yT+z2oF8Nc5awL+Tj9Vja5Y5snh+Vi/6eXrR903NTXMTc3NIcrPBOkqPIeBK4Ip+2l7Vbr/NTXPT3Oz4vLlyMe9/imNmO/I9e45vkBdtnzfbTrgWEnIM+T4JDxUPyG19B1XNMv9KrlrOI3/DsyM1CVcssxa58jmb/Puzu4oHYeE2E+4T5Crei8U6Plh2INa1/wBwDbmS9iy54ndATXy54mCfW+zrZ8l3YW+YcMXfbwF+TH7BeZGcfLvXxH9YHKypbz28PuGWIl+B0ve4XAmsV3dw95dwfY/NPk0ev/eTXyCeLrZzN/nu70v1c7wsVDwe29f3iXwn9nuK5/33wJqdJFvNej9fPMbzyMOd198Eaxr5t7u1f6cGU+0yY8lvVu+peW5Pqt1vYFHyjcY2G0j/Bzphbo7K3GzweJQ+VlXJzQbtX3cMYW6amxXITcrPM+OqmJvkD1F/I/+04JniGPw8MKZmme8ADxTrfIw8HOzmdes1N83NvuPL3ByE3GzQvuExRC5avFL73DRY5i7gY+amuWluDjw3aeG8WSy3WXHMzAXuJxf6Fq6Jd3ze7LvL9YgTEePID9YmKaVGw6hpiETE/sBHU0rvK/7eh/yisURXOzZIIuIL5FEIPtDtvlSRudk95qaaMTe7x9xUM+Zm94yC3NyePOjI+NT4/jxqwtzsnlGQmx2fN9u9B4/Uih8DS0fEm1NKz3S7M0Pgn+TL76SqMTel3mRuSr1ppOfmG4FPW9xRBY303Oz4vGmBR4Mu5TuCH9PtfgyVlNIp3e6D1AlzU+pN5qbUm0ZBbv6i232QOjEKcrPj8+aI/YmWJEmSJEnSaDEc49RLkiRJkiRpCFng6TERMTMiDhngOiZERIqIZQerX9JoZ25KvcnclHqTuSn1JnNzZLPAM4giYlxxoG/c7b50IiKOiIjpEfFCRPjbPY0YVc9NgIjYLiL+GBFzIuLpiLi6232SBqrKuVnz5rbR9NFu908aiCrnJrz6AbY+L4/tdr+kgapybhZ9Py0i7ouIucW/34qIxbrdt5HEmyyr1iLAr4BpwFe62xVJfSJiZ+AM4AhgH3Jx/t3d7JMkrgNWqJv3RfKoF5cMf3ck1TkK+GHN3893qyOSAFgHWAj4HHAP8A7gFGAZYL8u9mtE8QqeYRQRa0TEryNiVnGVzM0RsUODRZeIiJ9GxPPFsofUrefNEXFKRDwWEc9FxO8Ho4qbUvp6SukE4JaBrkuqkl7OzYhYCPge8F8ppZNSSnellO5MKZ0zkPVKVdDLuZlSmp9SmlU7AbsB56aU/CCpEa2Xc7PGc3U5al5qxOvl3EwpXZpS2ieldFlK6b6U0u+AbwIfGch6tSALPMNrCfK3etsC7wIuAH4VEevULXcwcCf5G/ojgWMiYleAiAjgd8BKwA7AhsAfgKsjov6bxFcVl/JNGtS9kUaOXs7NjYCVgfnFSXpWRFweERu2v5tS5fRybtYvPwFYi/xtpDTSVSE3D4mIJyPi1si3IVi4nR2UKqoKuVlrSWB2m23UTErJaZAmYByQgI3baHM98NWav2cCV9QtcypwbfH/95IvMV2sbplbyd/wA0wo+rFsTXwGcECLfdotHxrdf0ydnAZjqnJuAh8r2jxQ5OZGwOnAM8AK3X5snZwGMlU5Nxv062fALd1+TJ2cBmOqem6SP7xOBMYDnwWeAE7t9uPq5DTQqeq5Wbe+VYvcPLjbj+tImrwHzzCKiDeSK6Q7kH+3/wZgUeCvdYv+scHfuxb/3whYHHg8F1dftSiwRtm2U0r1VVtJhR7Pzb4rLb+ZUjq/6O9+wPuBvYDj+mkvVVaP52ZtP5cptndwq22kKuv13EwpnVjz518j4lngvIg4NKX0ZH/tparq9dys6edywKXAFcB3Wm2n/lngGV5TgH8DDiHfWGoOcDbQziWjY4BHga0bxJ4daAelUaqXc/OR4t87+maklF6KiHuAVQawXqkKejk3a+0FvAx4byyNFlXJzT43FP+uCVjg0UjW87kZEcsDVwO3A59KxeU8GhwWeIbXVsDZKaULACKirwp6d91ymzX4+87i/zcDywGvpJTuG8K+SqNJL+fmTcA8YG3g2qJ/Y4r+XTaI25F6US/nZq3PAr9MKT0zROuXek1VcrPPBsW/jzRdSqq+ns7N4h4+U4G/AR9PKb00mOuXBZ6h8vaIqD9YZ5ATa5eI+DXwT/Llc4s2aL9ZRBwOnE/+feNewCeK2JXAdODXEfFfxXqXJ1dqr0wpXdOoQxExA/h+Sun7ZZ2OiFWApcm/7SQi+k6G9yZHHtDIULncTCk9GxE/AiZHxN/Jv5s+AFgK+EkrOy1VQOVys2a5rYB1cYhXjUyVy82I2Jz8YXUq+X51m5B/AvKblNKDrey0VAFVzM0VgWnAP4AvAcvW/ATs8ZTSy/3ss1pggWdoNLpE+53k3+afBlxDvlv4f9M44U4k3xTuCOAF4Ot9995IKaWI+BBwNPBj4G3kS+imky+/K7M2sGw//T4K2Lvm777h0ieSk1Gquqrm5peB+cBZ5N9E3wxMTCn5TaRGiqrmJsC+wJ0ppektLCtVTRVzcx6wB/mD7SLkQQp+DHy7SRupaqqYmx8gjza5FlBfbF2N/CWmBij8yZskSZIkSVK1jel/EUmSJEmSJPUyCzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxI6LAExGTIiL1M40bwu0fGxE3DtX6e0VEbBgR10bEixHx94g4vI22i0fEHcVzsX5dbPOImBoRT0fE7Ii4PCLePVjbVneYl8Ojk9yIiFkNnotJdcv8MCJuioh5ETGjn/WtGxFzIuKJAe6OhoG5OTyGMDebnjMj4v0R8dtiXS9ExK0R8akh2EUNMnNzeHT6njIidoiIGyJibpF7l9TFV4uIi4vz4WMRcUJEvKEmvnJEnBsRd0fEKxHxo8HeNw0Nc3N4DGFuNn1PGxFLRMRZEXFbRLwUEZcO1j51w4go8ABTgBVqpruAE+rmPVTfKCIWHsY+DrvB3L+IWBq4ApgJbAQcAnw1Ir7Q4iq+B9zfYL1LAZcA/wf8K7AVMBu4PCIWG6RtqzvMywZ6KC+/woLPxZQGy5wO/KyfPiwK/ByY2nLH1W3mZgNVyM1WzpnFvJuBXYHx5Dw+MyI+MuAd01AzNxvohdyMiI8BZwGnAe8CtgTOrom/AbgUeAOwBbAX8CngmJrVLAY8BnyTnKOqDnOzgSrkZo1m72nHAnPIn1cv66D7vSWlNOIm4HZgUoP5PwfOB74G/AN4qJg/CzigbtnrgSk1f+9RrHcu8BT5w8wywP5Aqps+VtKv1YCLyG/GXgDuAHatia9c9PFJ8kF2E7BVTfxA4D5gPnA3sHdNbNFi2/sBvynWf3QReyf5pPM88CjwU+CtbT6m/1H06w01844G7m+h7R7AreSkS8D6NbGtinkr1Mx7R+1yA9m2U+9M5mXv5GWjx7bJsl8FZjSJ/7CY9gee6PZx5tT+ZG5WJzdp4ZxZ0u43wDndPtac2pvMzd7ITWDh4rH9VJNldgFeAparmffZov+LN1j+SuBH3T7GnDqbzM3q5Gbd8k3f0xbLnApc2u1jbCDTSLmCpx3bAasD2wIfbKVBRKwKnAP8iPxGahtyckCuGH4f+AuvVXAvLFnVKUAA7yEnwiHAs8U2lgSuAZYHdiJ/6/atYnki4uPk6vG3gfWLvpweEdvWbeMo4FfF+n8cESsDfwD+TK6GbgcsC1xQs3//VlxauFmTh2FzYFpK6Z818y4DxkXECmWNImI14L+BPYF5DRb5G/lF6LMRsXBxNcBnyd9O3jOQbatSzMthzMvCERHxRETcEhH/VXsZeasiYlfy435wu21VGeZmb+VmK+fMRpYs2mnkMDeHLzc3A5YDUuSfPD4SEZdExDvr1vuXlNKjdetdnPwFp0YPc7O3cnPUGdvtDnTBs8B+dQdPf1Yi/5zt/JTSrGLebX3BiHgBeKkmVmZV4LSUUl/b+2piewNvAXZOKT1dzLu3Jn5I0bbv97onRsSmwKHky9n6/CSldGZN374NXJdS+lrNvH2ARyJifErpr+Rq613kinGZ5clV5VqP1sQeqW8QEWOBc8nV3TsiYp36ZVJKsyNiIvC/wKRi9v8B26aU+gpCbW9blWNeDlNeFk4gn4ifIl9KfgywCnBAk20tICJWAU4CPphSmhsRrTZVtZibPZSbLZ4zFxARuxXr+lKT/qp6zM3hy83Vi3+/Qb7K4O/AQcDvI2LtlNLjRdtH69rVrlejh7nZW7k56ozGK3j+2mbCQX6zdQ1wV0T8MiL2i4hlOtj2fwNHR8T0iDgqIjaoiW0I3FSTcPXWAabXzbsWWLduXv0NuDYCto2I5/smXkvmNQBSStemlNZJKf2l3R3qx2TyTzZ+ULZARLyJ/JvIq4FNyZef3wn8uvhmUqODeTl8eUlK6fiU0rSU0l+LE/mXgM9HxBKttI9czfkZcEJK6ZbB7p96irnZQ7nZ7jmzKAadCeyfUrp1sPurrjI3hy83+z4vHZlSujCldCPw7+Sfsew5yNtS9Zmb5mZXjcYCzwsN5r1CcXlajVcviS6SdCLwIfJvGT8H3BMR72hnwymlk8gH+k/IyfKniDisnXU0Wm3d3/X7N4Z8Gd8GddNaLFiN7c8s8iVwtZariTXyPuCDke9G/hKvVWVvjYjTiv/vRa7K7ptSuiml9Efg4+QXmR0GsG1Vi3k5fHnZyA3kx3qNFpdfhHwTu2/V5PcPgGWKv/dqY9vqbeZmb+VmK+dM4NXizm+BL6eUTm9jm6oGc3P4crPvyoE7Xu1sSvPJV8+tMoD1amQyN3srN0ed0VjgaeRx8u8ZAYiIN5IPylellF5JKU1PKR1JrlTOBj5ahOcDC7WyoZTSgymlH6WUdiPfRX+/InQL8O6IeEtJ0xnkD1S1tqLmgC5xM7Ae+QZV99ZNz7fS58IfgQmx4L0AtgVmppTKLjXfk/y7475E36WYvxvw9eL/iwOvpJRqXzxeIb+Y9B2fnWxb1Wde9m+wcqPvG55W34TOI//2uvZE/k3g6eL/v2lj26oec7N/Q5WbrZwziYj3kW+0eWhK6YdtbE/VZm72r5Pc/BP5Bspr980obkOwGvBAzXrfFRFvrVvvHPK9UzS6mZv9G6rcHHUs8GRXA3tHxFYRsT75cuZXRcTWEXF4RGxc3HdiF3KS9h3wM4E1ImJ8RCwbJUPGRcT3I+IDEbFaRLybfND2reNs4DngwojYMiJWj4hdImLrIn488O/FJXtrRcTB5KT/dj/79t2irz+LiE2K9X4gIk7r62ex3zMiotlN4M4mv4E8LSLWi4jdyTdWrR2+dY9iPcsCpJTuSynd3jfx2g0g700pPVz8/zJguYj4bkSsHfmmWGeSq8O/b3XbGpHMyyHIy4h4T0R8MSLeVezzx4v+/CLV3Byy2J8NyFcLLBIRGxTT2JTdXpff/wBeLv4uu/xXI4O52b3c7PecGfmGmBcB3wEuiIjli2nZfvZd1WduDs372SfJQzB/MyLeHxFrk69aXYR8r0mA35GvGvhpkcPbke+hdVJKaU6x3ug7lwJLkK963SAa3KNSI4652b3cbPqetmaZ9YpllgbeVMTH97PvvSn1wFBegz3Rz9B1DeYvBfwSeAZ4iDwixatD15G/qb4MeAx4kTxs3H/UtH8j+dK0p2k+dN3J5Bf/F4t1ncOCwymOIw+t9wz5zdqNwJY18fqh6/apifUNXbdDg+2uU9O/OeQK7XeBhYr4vxVtN+vncX03+beZLwIPA1+pi/cN47d8Sft1aDCUK/kO838s9vsp8tCRm7azbafen8zLhvkw7HkJ/Cv5G4+nyTe7u5M8bOSide2u5/XDcjbLb4dJr+hkbr5uu72em03PmcXz1ih3mw4N69R7k7n5uu127f0s+QPjd4r9fRa4Chhf12514JKib48DJ7LgkM99+2ZuVnwyN1+33V7PzX7f05KvlK2Pv9jtY62TKYodkiRJkiRJUkX5Ey1JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKs4CjyRJkiRJUsVVvsATEftExPPd7sdwiYgzI+KibvcDICLGRUSKiI273Rf1HnOze8xNNWNudo+5qWbMze4xN9WMudk95mb7er7AExG7RsTVEfF0RLwQEbdFxDcj4m3d7lsjg5UQHsydiYj1IuL8iLivePwmdbtPI5W5aW62KyKWjIjvRcQ/ImJeRNwbEbt3u18jjblpbrYjIqYVj1v99Ldu922kMTfNzXYURYVGublot/s20pib5mY7ImLfiLgmImYXx8zUiNiq2/3q09MFnoj4JvBL4FZgB2Bd4CBgNeBzXezaiBMRb+h2HwbJ4sBM4KvA/d3tyshlbg6fkZKbxX5cAawF7A6sDeyDeTqozM3hM1JyE9gVWKFmGgc8B/yii30acczN4TOCchNgDgvm5woppRe726WRxdwcPiMoNycA5wHvBf4VuAu4LCLW6manXpVS6skJ2BRIwMEl8bcU/+4DPF8X+3/AvcD84t99G8TvBl4EngAuA8bWxD8N3LblUigAABwTSURBVFHE7wb+AxjTX3tgUtHn2mlCSf/HAF8DHgLmAbcBH66J169nWjH/TOAi8gvPw8Bs4Axg8Zq2AfwX8H/A3GLdn6yJjyvW+XHg6mKZA0r6uTBwDPBA0c/7gC/WrWfj4u+FgNPIH9jmAvcU/ah97N4JXAU8CzwP/AWYWMTeAHwP+EexrYeAYwdwDN0OTOr2sTzSJnPT3OwkN4H9ij4u3O1jeKRO5qa52UluNuj/J4CXgJW7fUyPlMncNDc7yc1Gx4OTuWludj83G/Q/gFnAgd0+plNKPV3g+W7xhLyhn+UWSDhgF+CfwAHA24EDi793LOIbk9+4fAJYFXhXkVBji/i+wCPAbuTK7Y7FE3ZAf+2BJcjVvCuA5Yup4YeZos2zwJ5FP48CXgY2KOKbFAfzdsV6lq5JuGeAHwPvAD4APA0cXrPub5Irif9W7MOewAvA9nWJMrNmP/+lpJ/nAn8HPgKsDkwE9ipJuDcU+7FJEdu96Nu/16zvNuCnwDrAmsXztXkR+88iyd4DrAJsAXy6pu0kILVxDFngMTfNzR7JTeBi4CfAycXzdkfRrulx5GRumpvDd94s2kwDLur28TySJsxNc7Oz8+Y+xeP4QNHvi4ANu308j6QJc9PcHJzz5iJFHz7ZTrshO6673YEmD9TFwF9aWK4+4aYDp9ctcyZwbfH/XYsD9k0l63sQ+FTdvC8Bd7TY/kxaeGNEroZ+vW7eNOCnqcHBXLf+h4CFaub9GLiy+P8bydXMreva/Tdwcd26/7OfPq5VLPdvJfGGfaxb5ti+vhV/PwvsXbLs98jV1iiJHwDMaOMYssAzBJO5aW42iPebm8AM8rdQpwMbkU/is4Ap3T6mR8pkbpqbDeLtnjffXvTvw622cWrpcTU3zc36eCvnzc2BvYENgK2B88k/2Vqr28f0SJnMTXOzQbyt82bR5nhygWrJbh/TKaWevgdPdNjuHeSkq3Ut+feUkKudDwD3R8Q5EbF3RLwJICLeCqwMnBwRz/dN5INmjf7atyoilgRW7KefzdyRUnq55u9/AH03AVsXWBS4tG4fPlezD31u7Gc7GwKvAFNb6BMAEbF/RNwYEY8X2/0PcnW0z4nAqcWNzI6IiHVqYmeST2J3R8QPImL7iHj1GE0pfT+lVLu8usPcLGdulhsDPEa+hPmmlNIFwNeBz0VEp8eUFmRuljM3W9P3rfLv2mij/pmb5czNEimlP6aUzkop3ZpSugbYg/xzmANb3Qf1y9wsZ2621peDyD+n2zWl9Gyr7YZSLxd47gbWiIiFB2l9CSCl9BzwbvLlXA8ChwMzImJFXns89ic/8X3T+sB6LbQftH72458N2vT1ve/fHVlwH9YjX15X64UO+9hQROxBrtyeSb7UbwPgJPLvKnNHU5pEflG4kHxJ3F8j4jNF7GZylfbwYj/OAq6oTTr1BHOznLlZ7hHg7ro3C3eSb4y+bAe7pdczN8uZm/33ZWHy1QJnpJRe6nSf1JC5Wc7cbFFx/ryRfMWDBoe5Wc7c7L8vXwKOBj6UUvpTp/s02Hr5g/PPyJd/HdAoGBFvKWl3J7Bl3bytyPd7ACCl9FJK6eqU0uHA+GI7O6SUHiVXJ9dIKd1bP/XXvgjPJ9/8qVRR3ftHP/2cX/zbdF0N3EG+YdSqDfbhgTbXdSv5GJnY4vJbATcUlc+bi8esvopLSumelNL3Ukrbk2+S9dma2HMppfNTSp8DtiffnXzNNvutoWVuZuZme7k5HViz7gT6dvLl5k+0sR6VMzczc7Oz8+bO5GLraR20VXPmZmZuDuA9bXG163jyFyYaHOZmZm62mZsRcTDwDfI9h65tp+1QG9vtDpRJKd0QEd8Gjo+IfwEuIP+2bTXg38l3K5/coOnxwC8j4ibgcvKNnz5B/i0jEbED+SD4A/AU+WB6EzlRAY4E/icinib/LvMN5ArqSimlb7XQfibwwYhYG3gSeCalVF8B7evnURFxD3AT8Eny72vfXcQfI/+2cbuImAm8mFJ6poXH7bmImAJMKU4EfyDfjGsz4JWU0in9raNmXXdHxC/Il7gdBNwM/AswLqX0kwZN7gb2iYgPkp+fjwHbkO+8TkQsBkwhD0U4E1iOIkmL+MHkk9at5KrxnuTfUP69iB9AvvlY6WVzRQW+77LDRYHlI2ID8u9m7y1rp9aZm+YmHeQm8EPyG6jvRsT3yd+eTAZOSim18k2S+mFumpt0lpt99gOuSind1+r+qjXmprlJZ+9pjwSuJ48StCTwRfIHfYfuHiTmprlJZ7n5ZfJNpj9J/qnX8kVobiuP35BLPXAjoGYT8FHyzaCeIV/idTv5AX1rEd+H1w9btz/5Cf8ndcPWkZ/gqeRkmFus79N17T9OPrheJB8s1wIfa6U98FZyoj8HLQ9bN598t++d65b5LPmyvJepG7aubrlJwO01fwf597l91dXHyb/l3LaIj6Ofm1XVrGsR4Nvkm3TNI//294BG6yFfGnda8Zg9Xfz/68DMmvjPyMk2j1xVPoXihlTk3/7fXDx2zwK/B7ao28/UT3/7+lQ/Tev2sTzSJnPT3Kzbz9RCnzcDriuen/vJoyA4bPogT+amuVm3n6mFPq9OvgfC7t0+fkfyZG6am3X7mfrp73d4bdjox8jDZG/e7eN4JE7mprlZt5+pn/7OpPHnzTO7fSynlPLdoyVJkiRJklRdvXwPHkmSJEmSJLXAAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFdfxMOkR4d2ZNdo9kVJ6a7c7Uc/clMxNqUeZm1JvMjelHpRSinbbeAWP1LkHut0BSQ2Zm1JvMjel3mRuSiOEBR5JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKs4CjyRJkiRJUsVZ4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqbmy3OyBJg2Hq1KmlsQkTJpTGpk2bVhqbOHHiAHokSZIkScPHK3gkSZIkSZIqzgKPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnMOkSxoRmg2F3mm7TodXlyRJkqTh5hU8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKs5h0gfRD3/4w9LYfvvtVxr7xS9+URr7zGc+03D+3LlzW++YpI44TLokSZKqbplllimNHXTQQQ3nf+1rXyttk1IqjUVEaezcc88tjZ1xxhmlsSuuuKI0pgV5BY8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWco2i1aYUVViiNbbfddqWxZnca/+hHP1oaO+WUUxrOnzp1amkbaaSaNGnSiN6eJEl6zW677VYae/e7310aK3v/DDBz5syBdEnqWSuuuGJp7JJLLimNrb/++g3nP/TQQ6VtzjrrrNLY+PHjS2Mf/OAHS2OrrrpqacxRtFrnFTySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqzmHS2/TII4+UxmbNmlUaW2WVVTra3le+8pWG82+88cbSNs8991xH25J63TbbbNPtLkgaRJtuumlp7LzzziuNjRs3rjSWUmo4/5xzzilt86lPfao0JmloffjDHy6NNRuKefHFFy+NzZ8/vzQ2adKklvolVc3WW29dGisbCh3g4osvbjj/sMMOK23zt7/9rfWO1bj55ptLY5tssklprNnw6s2GgB+NvIJHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqzgKPJEmSJElSxTlM+iC67rrrSmPNhoJtZuLEiQ3n77DDDqVtzj333I62JfW6CRMmdLsLkhpYcsklS2OnnnpqaazZ8Mhjx5a/RXnllVda61iNPfbYozQ2efLk0ti9997b9rak0WrNNdcsja233noN5++1116lbRZbbLHSWEqpNLbVVluVxi688MLSWKemT59eGrv77rtLYzfccEPD+bNmzRpwn6Razc6bJ598csP5nQ6F3qmFFlqoNLbRRhuVxhwmfUFewSNJkiRJklRxFngkSZIkSZIqzgKPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniHCZ9EF1++eWlsQMOOKA01mwo2DKbbLJJacxh0lVlkyZN6nYXJLVpl112KY195CMf6WidZcMHQ/Ohy8uGXl9iiSVK2xx66KGlsX333bc0Jo1Giy++eGnsN7/5TWlsnXXWGYruNPTe97532LYFsNNOO3XU7rnnnms4//HHHy9tc+yxx5bGTj311I76oZFh6tSppbGJEyeWxq699tqh6E7bXnzxxdLYFVdcMYw9qTav4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxDpM+iJoNk95suNctt9yy7W3tsccepbGTTz65NHbXXXe1vS1Jkjp1wQUXlMbmzJlTGvvyl79cGms2hPA3vvGNhvO/8pWvlLYZP358aUzSgjbccMPS2GAPhX7RRReVxo477rjSWEppUPvRn4gojTXry6qrrtpwfrNh3pdZZpnWO6ZR5bHHHusoNtiafbZdb731SmMPP/xwaazZZ2ktyCt4JEmSJEmSKs4CjyRJkiRJUsVZ4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZzDpA+To48+ujR2ySWXtL2+5ZZbrjR2/vnnl8be+c53tr0taSSbPHlyt7swIJMmTRrR21M1nHXWWR3FhsKVV17ZcH6zYdJXXHHF0thKK61UGms2pKs0Uk2fPr00NmZM+XfHv/71rxvO33HHHUvb3HbbbR31oyquu+66hvPPPffcYe6J1J4ll1yyNHb22WeXxsaOLS8/nHzyyQPqkzKv4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxDpM+TJoN5Th79uzS2FJLLdX2tt785jeXxpoNaffss8+2vS1Jg2PChAmlsSOPPLKjdkNhm222KY1NnDhxGHsiDZ5mw6Qvv/zypTGHSZcWtPjii5fGVl999Ybz77vvvtI2Rx999ID7JGnwTZo0qTQ2bty40tj9999fGjvvvPMG0CP18QoeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFecw6cNkzpw5pbETTzyxNPaNb3yj7W2ttNJKpbFmQxz/9re/bXtbklrXbEjzqVOnDl9HBqCTfXD4dEkaHXbaaafS2Hrrrddw/uWXX17aZu7cuQPuk6Tmxo5tXBI47rjjStt88YtfLI29/PLLpbE999yzNDZz5szSmFrnFTySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqLlJKnTWM6KyhXmfhhRcujV111VUN52+++ealbSKiNHbjjTeWxj70oQ+Vxp588snS2Ch2U0pp4253ot5Izs1OX6861SyXyvTSUOjTpk0rjU2ePLk0duSRR5bGmu1fmU4exwEyN0exbbbZpuH8q6++urTN888/XxrbdNNNS2N33XVXaextb3tbR+tsNoRspy6++OKG83/6058O+rb6YW6OcCuttFJp7KGHHmo4/9BDDy1tc/zxxw+4T2qJuTmKffWrX204v9l7xZdeeqk0dvDBB5fGfvCDH7TeMZFSavtNtFfwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqzgKPJEmSJElSxY3tdgcE8+fPL43Nmzev4fxmo9KMGVNet9t44/Ib5K+wwgqlMUfRklrTbASqodBspKyJEyd2tM5mI2V1MoqW1Otmz55dGttuu+1KY0cddVRpbIsttiiNrbjiiq11rM7pp59eGrvwwgtLY81G0JQGU9n71ma233770pijaEmDo9lnwC9+8Yttr+/MM88sjTlSVnd5BY8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKc5j0Hvfggw82nJ9SKm3zyiuvlMaatdtxxx1LY7fffntpTBoukydPLo0NxfDkkyZNKo2VDU8+3MOIdzoUulRlK6ywQmlsl112aXt9K6+8cmnsO9/5Ttvr68+MGTNKY1OmTCmNnXfeeaWxOXPmDKhP0mB46aWXSmOzZ89uOP9d73pXaZstt9yyNDZ9+vTWOyaNAuuvv35p7LLLLiuNveUtb2k4v9lQ6Pvvv3/L/dLw8goeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFRfNhs1u2jCis4Zqy0477dRw/q9+9avSNhFRGmv2fF977bWlsfe///2lsWZDYo5wN6WUNu52J+qN1tzs9LWsU8M5THqzodDL+jEQU6dOLY11sn/NXpOGiLnZQ8aNG1caO/vss0tjiy66aGls+eWXL42ttNJKLfVrMFx//fWlsWOPPbY09uc//7k0NmvWrAH1qceZm6PYz3/+84bzd99999I2zXJlm222KY29+OKLrXdMYG5WxhJLLFEaazas+S677FIau+aaaxrO33XXXUvbPPXUU6UxDZ6UUttvor2CR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKs4CjyRJkiRJUsU5THqPGzt2bMP5t9xyS2mbddddtzTW6fN92GGHlcamTJnS0TpHAIeU7CGDPbT3cJs8eXJpbNKkScPXEQZ/yHmHSc9Ga26WDY0M8NGPfnQYe1LupJNOKo1961vfKo01G9L8lVdeGVCfRihzcxRbfvnlG86//PLLS9usv/76pbFmry177rln6x0TmJuVscEGG5TGbrrpptLY/PnzS2MTJ05sOP/6669vvWMaEg6TLkmSJEmSNApZ4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIprPAa3esZLL73UcP6pp55a2ubEE08c9H7ssMMOpbFRPEy6ekjZEI9QjSHUjzzyyNLYtGnTOoo10+wx6VSzod41eh144IGlsbPOOqs0tuyyy5bGttxyy9LYvvvu23D+HXfcUdqmWR8lDY4XXnih4fzZs2d3tL5HH310IN2RetbCCy9cGmv2GbCZo446qjTmcOgji1fwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqzgKPJEmSJElSxVngkSRJkiRJqjiHSa+ou+66a1i3N378+NLYqquuWhp74IEHhqI7UluaDaE+adKk0lizocuH01AMaT4UOh2yXSPb448/Xhq75JJLSmOLLbZYaaxZ3pa58MIL224jafCccMIJDedvvfXWpW3uvvvu0tixxx474D5JvWjnnXcujW244YYdrbPT4dVVPV7BI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIcJr2iLr300tLYn//859LYxhtv3NH2llxyydLYgQceWBo75JBDOtqeNFyaDbfcbNjvsiHUJ0yYMLAO9YBm+z158uSO2knt+sIXvlAaGzduXNvrO+eccwbQG6m7Pv3pT5fGbr755tLYX/7yl0HvS7P3hBdccEFp7H3ve1/D+XPnzi1tc8QRR5TGHn300dKY1OtWX3310tgxxxzT0TqnTJlSGps9e3ZH61T1eAWPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkinOY9BHod7/7XWms02HSm9l8880HfZ1SL2g27HdZrNkw6c1iZcOuD5Vm+zZx4sTh64hUYrHFFuuo3YMPPthw/pNPPjmQ7khd9bnPfa40Nn369NLY4YcfXhpbeumlS2Of/OQnS2Nf+MIXSmMrr7xyaWzevHkN5zcbErrZsOtSlTXLo9VWW62jda6yyiqlseOOO66jdQ6na665pjR24YUXDmNPqs0reCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFVcpJQ6axjRWUMNudVXX700dvfddw/69g477LDS2JQpUwZ9ez3kppTS4I87P0DmpmRuVsUiiyxSGrv66qtLYxtvXP707r///g3nn3HGGa13TEPF3OzQ5ZdfXhp7z3veUxq74447SmNLLLFEaWzNNddsrWN17rnnntLYEUcc0XD++eef39G2NKjMzWF21VVXlcYmTJgwfB3pIdddd11pbOuttx7GnvSOlFK028YreCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFXc2G53QIPvoYceKo0dfPDBpbGy4SsBFl988dJYs6FsJUkqc+CBB5bGNttss9LYM888UxpzOHSNROecc05pbPz48aWxDTbYoDQ2b9680tiMGTNKY5/5zGdKY7fddltp7IUXXiiNSWrNE088URr7+c9/Xhprlu+nnHJKw/nNXj9uvfXW0lgz6623XmnM14jB4RU8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKs8AjSZIkSZJUcZFS6qxhRGcNpZHjppTSxt3uRD1zUzI3q+Kaa64pjW2xxRalsZNOOqk01mxkLnWduTkEVlpppdLYqquuWhqbM2dOaazTEXJUWeam1INSStFuG6/gkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHFju90BSZI0Oj399NOlsa9+9aulsRNOOGEouiNV0sMPP9xRTJI08ngFjyRJkiRJUsVZ4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIpzmHRJktQVO+64Y7e7IEmSNGJ4BY8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKs4CjyRJkiRJUsVZ4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVdzYAbR9AnhgsDoiVdCq3e5ACXNTo525KfUmc1PqTeam1Hs6ystIKQ12RyRJkiRJkjSM/ImWJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLF/X8jbCG0nLBDgwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(imgs_min[i].reshape(28, 28))\n",
    "    plt.title('Model prediction: {} (p={:.2f}) \\n Label: {} \\n Trust score: {:.3f}' \\\n",
    "              '\\n Closest other class: {}'.format(pred_min[i].argmax(),pred_min[i].max(),\n",
    "                                                  label_min[i], score_min[i], closest_min[i]), fontsize=14)\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### High Trust Scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The high trust scores on the other hand all are very clear 1's:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAETCAYAAABN+bF9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZxcVZ338c8vhLCqCEECyCKIBBQIIhiUKD4MIgqI4AIIsrgy4jgqosijIiMuiAuKOCOiyODoM8iioo4Ksgg4gCyy77uyy5YQCEnO88c5DTeVut1d3V1ddTuf9+t1X0ndc++pU1X3W8uv7xIpJSRJkiRJktRck3o9AEmSJEmSJI2OBR5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWq4vijwRMTUiEgRsW0H6xweEdd0cVgdi4hrIuLwyu07IuLgUfZ5bkQcO+rBjaOI2CAi7o+IF/R6LN0QEadExCd6PY7xYDYH7dNs9hmzOeQ6ZrNPmc2Jw2wO2qfZ7DNmc8h1zGafMpv1hizwRMSJJQwntGn7amk7cyR3vgTYEjhuOAtGxH4RMbtN027AoWM6qmGKiGMi4i8R8VRE3NHBql8CjkspPdaloQEQEf8cEbeX8V0WEbOGsc6HI+L6iJgbETdGxHvaLPP8iPh2RPw9Ip6OiFsi4p2VRY4ADuv1G4rZHBWz2UXdyGZEvKM85kcjYk5EXBkR+7Z0Yzabz2x20QizuVfJ25MRcV9EnBwR01qW+WhE3FDye09EfDciVqwsYjabz2x2UTeyGRFLR8TnIuLW0u9fI+JNLd2YzeYzm13Upe+0+5VtunVatrLYiLM53D147gbeGRErVAY2GXgPcFend9rPImLpseorpfRgSunJUfbxj5TSE2M1pg5NAn4MnDTcFSJiLWBX4EfdGlS5n3cBx5DDvTlwEfDbiFh7kHUOBL5KDszLgc8D342InSvLLA38AdgAeCewIbAfcPvAMimlq4HbgL3H9EGNjNkcAbPZPd3KJvAw8EVgJrAp+XGcEBFvHljAbPaG2VzERMvma4H/JD+ml5dxbgz8pLLMXsBRwJHARuRt/M3lvgCz2StmcxFLXDbJn5kHAv9S2v4dOD0iNh9YwGz2htlcxETL5nC+0wI8CaxenVJKTw00jiqbKaVBJ+BE4EzgcmD/yvy3kn/0/hg4szJ/EvBZckifBq4G3trS55bAZcBTwBXAW4AEbFtZZmPg18ATwAPAT4FplfbDgWsGGfe6pc+9gAvKfd0AvLGyzLZlmTcDlwDzgJ1K286VMd5O/uIypbLui4BfAHOBO4EDgGuAwyvL3AEcXLn9AuB7wL2l3+uBd1XGUZ0OL+ucCxxb6eOF5Tl/pNz3WcDLK+37AbOB7cp45gDnAC8Z6rUe5Lk8GLijg2WvaJk3MKadgZvKYz8HWG8UY7oYOL5l3s3AlwdZ5yLgmy3zvg5cULn9AXKYpgxx/5+rrteLCbNpNpegbNasd3lrv5hNs2k2xzKbBwN3tszbH5hduX0scF7LMl9o3dYwm2bTbI53Nv8OfLRlmVOBk82m2cRsdiubw/m9uR+VrA7S14iy2ck5eE4oG9WAA8hVs9Sy3EeBTwKfAjYBTgdOi4gZAGWX3V+Tf0S/Cvg0cHS1g4hYHTi/bDBbAf8ErAj8IiI6PW/QUcC3gRnkPTN+ERFrtizzVeD/AtOBiyNiB3IF/Fhy5e0A4O3k6t2AE4GXlrHtSq4ur1s3iIgI4DfA68lvwBsDHyeH/CLgX1m0knd0+544EXg1+Q1vq7LO/0TEcpVlliHvZncAsDWwErlqPzCWdctuYPvVjXcUZgF/aTN/GXIFc/8ypqXI20WUMc2KiNlDTJ8py04BtgB+33IfvwdeM8jYliGHvWousFWlkr4rcCHwnbK763WRj79trbRfUtZbjt4zm2ZzOJqezWdFth15D7vzW5rNptlsdSJmc6TZvBBYPSJ2LrmbCuxBfl0GXADMiIiZ5b7WBnZpWQbMptlc3ImYzW5ms+6zdZuWeWbTbLY6EbPZ7e+0y0XEnZEPaz4zKnvWVYwsmx1UVF9YBrcBMI1cLV17oL2y/N+Az7X0cS6lWkzeQ+JRYMVK+95UKqrkXZrObunjhWWZrTqsqB5WmTeJXNH7YktFdfeWdc8HPtsyb1dyVTCAl5X1XltpXwdYQE1FFdgeWAhsVDPe/WhTyaNSUS3PfQJeV2l/AfAY8L5KPwnYsLLMu8vrFeX2muTq8tu6UFG9EvhCm8dW93z9U7m9HPkNbLBp5bLsGq3PQ6XKeeMgY/sScD+5oh/kN/z7Sl+rl2VuIIfyh+RQ716WObqlr03Leut3WlUdqwmzaTaXoGxWntPZwDPknL63TV9m02yei9kck2yWZXYDHifnLpG/3C7XssyHyT8gBpY5aeC5M5tm02z2JpvAf5H33tiwbC/bk3+oP202zabZ7Onvza2BfckFwVnAz8nZ3GAssjmZYUopPRIRp5OrdI8C56aU7ioFMSCfnLY8GRe2rH4Bebc0yMdnX5VSqp7g6c8ty28BvC7anwRqfXI1a7ie7TultDAiLiZXM6taK4BbkKtln6rMm0TeKKaRH8PC6jhSSndGxN8HGcfmwL0ppes7GHurgfutPqbHIuJqFn1MT6eUbqzc/jswhfym9Y+U0t/I1eNuWI7Fq5ZQ/3xtDJyVUpoL3NKlMQ34N/LrdxE5cPeTdz88pIwP8uv8APD+lNIC4LKIWAX4ZkR8MpW0kT98ID/enjKbZnOYmp5NyLtQzyD/hW074BsRcUdK6ezKMmZzcWbTbI5IRGwMfIec0d+R/+L7NeA/yH9JJiJeTz5U4p/Ju7O/lHzOgi+QvwgPMJuLM5tmc0SGk03yXi7HA9eRfyTeSt4b5oCW7szm4sym2RypIb/TppT+TOX5jYiLyEWrj5DPmTVgRNkcdoGn+GEZ4GwW/dAejjT0Is+aRN6trt0l3+7v8H6HY06b+/8CcEqbZR+s/L+Tx9Rt1bHMr2nrdHfDkXiIHOx2ap+vyGck/+0QfX8ppfSlch8LgNVa2lcjV0jb33kO9QER8cGy7L3kCv8TPPe63gs8U4o7A64HlgemVpZbufxb3R56yWyazaE0PZuklBby3AfzlRGxEfAZoFrgMZuLM5vtmc0hskne/f6SlNLXyu2rImIO8KeI+ExK6R7yiVx/mlL6QVnm6sgnSf1BRByRUhp4bs3m4sxme2ZzDLKZUnoQ2DXylXlWIf8A/wr5sKUqs7k4s9me2Ryj77Qt6yyIiL+Q956qGlE2Oy3wnE3eBXcqcEabwT1eqmSvZdEv3NuQq8eQfyzvFxErpJQGNvSZLV1dTr6C0Z0ppWc6HGOrmcAf4dnjErci7wY1mMuB6SmlthW+iLiBvPFuRa7ODRxzvsYgfV5BPlZ2o5qq6jzycYKDub7c79aU806UKvYmdPks4h24gsUr1lD/fA08F38h/2V+MP8ASCnNi4jLyLshVt8UtyefPG5QZZu6p4xjD/IunwN7CVwI7BURkyrzXkbebe6hSjevAP6WUurGB8BImE3M5hCans12JpGPda4ym2azymwyqmwuT/6CWzVwe9IQy0TLPLNpNqvMJl3PJuU+ngL+Fvn8H7sD/92yntk0m1Vmk/H9Tlu2mU2Bv7Y0jSybgx2/lSrHRFZuPw94/iDt/0o+HnRP8g/jI8hvOJuV9hXJVaj/Rz6h1PY8t+vgwDGRa5APkzmNfIKn9cgnl/o+8LzU2TGRd5NPWLUheZfhp4AXp0WPiZzasu4O5ONZjyhP7PTSx1GVZX5LPmP71uQN5Y/kSnPdMZGTyLtiXVf6f0l57LuW9teUsWxPfkNbPrUcE1lun0HeSGeRg/bL8hiXSzXHVrY+ToZ5TCR5N+sZwDfIVf8ZZaq9whT5zOUPAZMr8/Yrz+cllefrHOAqWo7RH+5EPhv8POB95F0JjynP/zqVZU4CTqrcfhmwD7k6uhXwM/Kll9etLLMWefv9TtlmdiCH82ttcnHCSMY+VhNm02wuWdk8rGxr65V+P1HG/iGzaTYxm93K5sA4Dizb1GuBS4HLKsscTt5+96i8RrcAp5pNs4nZ7GU2X00+T8965Tk+m7z3zkpm02xiNruVzeF8p/18eX3WK+P8YRn7VmORzY4DN4xATuK5y9bNKxvlri3rvJpctXyaXKnamUrgyjIbkCufA5dnu5H8o3tKh4F7N7mK91TpY8e6DbFl/TcCfyLvufE4ueJ3UKV9tbKxzy2P9X0Mfdm6lcjHwj5YxnMd8M5K+/fKxpoY5WXrhgjcwHOz3xCv/blludZp3UHWWao8H29pHRP5TOw3l9f9POClIwlbpd9/Ls/x0+RLDLaeBOtc8rG7A7c3Ild8nySfKOwMKicHqyw3s2wzc8mXLDyCRS9ZuGxZf+Zoxj/aCbNpNpegbAJfLmOcS/7rykXAni3LmE2zObB9mc0xyGaZ9xHg2vK63ku+6suLK+2TyV9WB/J5N3Ac8EKzaTYxm73M5uvL6/JUeS1OAtZo6cNsms2B7ctsjkE2Gd532m+SL3v/NLnI+Dtg67HK5sBZrieciFiX/ON8y5RSu8uoqUsi4kPAO1JK25Xb+5HfNFbs6cDGSER8GHhrSumNvR5LE5nN3jGbGozZ7B2zqcGYzd4xmxqM2ewds1mv03PwSMNxPLByRLwgpfRYrwfTBc+Q/2oiNY3ZlPqT2ZT6k9mU+pPZrGGBR2Mu5StQfanX4+iWlNL3ez0GaSTMptSfzKbUn8ym1J/MZr0Je4iWJEmSJEnSkmI8rlMvSZIkSZKkLrLA02ci4o6IOHiUfWwbESkipo7VuKQlndmU+pPZlPqT2ZT6k9mc2CzwjKGIWLds6K/q9VhGIiIOi4gLI2JORHjsniYMsyn1J7Mp9SezKfUns6mhWOBR1TLAacC3ej0QSYswm1J/MptSfzKbUn8ym11mgWccRcT6EfGLiLivVC0vj4id2iy6YkScHBGzy7IHt/Tzgoj4fkQ8EBFPRMR5Y1HFTSl9LqX0deCK0fYlNYnZlPqT2ZT6k9mU+pPZlAWe8bUi8Ftge2Az4FTgtIiY3rLcx4HrgVcCnwe+FBG7AUREAL8G1gR2AjYHzgf+GBGr191x2ZXv8DF9NNLEYTal/mQ2pf5kNqX+ZDaXcJN7PYAlSUrpr8BfK7OOjIidgbcDX6zMvzildGT5/00RsSU5hKcBbwBmAKumlOaWZT5b+tkHOKrm7m8EHhqbRyJNLGZT6k9mU+pPZlPqT2ZTFnjGUUSsQK6Q7gSsDiwNLAtc1bLon9vc3q38fwtgeeDBXFx91rLA+nX3nVJqrdpKKsym1J/MptSfzKbUn8ymLPCMr6OBNwEHAzcDTwInAVM66GMScD8wq03b46MdoLSEMptSfzKbUn8ym1J/MptLOAs842sb4KSU0qkAETFQBb2pZbmZbW5fX/5/ObAasDCldFsXxyotScym1J/MptSfzKbUn8zmEs4CT3e8LCLmt8y7gRyst0XEL4BnyLvPLdtm/ZkRcSjwc2Bb4D3Au0vbWcCFwC8i4pDS7zRypfaslNKf2g0oIm4Ajk0pHVs36IhYG1gZWLfcnlGabkkpzR7sAUsNYTal/mQ2pf5kNqX+ZDbVlgWe7vhJm3mbkE9cdQLwJ+AR4Fu0D9w3gE2Bw4A5wOdSSj8HSCmliHgz+SRZxwMvIu9CdyF597s6GwJThxj3EcC+ldsDl697A3DuEOtKTWA2pf5kNqX+ZDal/mQ21VaklHo9BkmSJEmSJI3CpF4PQJIkSZIkSaNjgUeSJEmSJKnhLPBIkiRJkiQ1nAUeSZIkSZKkhrPAI0mSJEmS1HAWeCRJkiRJkhqu8QWeiDg8ItIQ07pdvP+vRMRfutV/P4iI6TXP67bDXH9SRPyxrLNTS9t9bfo9vE0f+0fEVRHxVEQ8GBHHj82jU7eYze7rcjY3iohfRcTDEfF4RFwYEdu1LPOSiPhNRDwZEQ9ExNcjYukxfIjqArPZfSPNZkT8b5t1TmxZxs/NCcpsdt9IsjnIOikiPlJZ7qSIuK1k7oGIOC0iXtamv50i4uKImBsRj0TEb7v0cDVGzOb4KN9LPxkRN0XE0xHx94j4whDrLB8R3yvfV2eX3K3essz3IuKy0ucNbfoY1ffpfjO51wMYA0cD/165fS5wZpk/4MHWlSJiSkppXneH1jtdenzbAjdWbv9jmOt9Bpg7RPuPKrdnVxsj4hDgX4BPApcAKwDrD/O+1Ttms40mZDMiAvgtcHXp+yngI8CZEfGylNLdpZDzP8BdwGuAacBJwEJyVtW/zGYbfZTN7wFHVG63+/z0c3NiMptt9EE2bwZWb5m3B/l1Oa0y73+B7wP3AFOBfwN+HxHrpZQWAkTEHsB3gUPJr+9kYLORPgiNG7PZRhce37HkbH4CuAZYmZylwXwX2A54B/AY8B3gjIiYmVJKleV+CLwK2HqQvrZlZN+n+0tKaUJN5I3h8Dbzfwb8HPgs8Hfg7jL/PuCglmX/Fzi6cvtdpd+55Bf6HGAV4ENAapn2qBnXS8hvBI8Ac4DrgN0q7WuVMT4MPAlcBmxTaf8IcBswD7gJ2LfStmy57w8Avyz9f7G0bUL+ETYbuB84GVi1w+d0eun/FSN4PV4L3AGsUfrYqaV9see/pX1V8o/LWb3etpxGN5nN5mQTeHGZt3Vl3vLV5YC3AfOB1SrLvK88xuV7vb05dbQtmM0+yWbr81izjJ+bS8hkNvsnm236+RPwyyGW2arc1zrl9pTyGu3T623LaXST2exKNjcr97t+B+tMJX8X3b0y76VlnK9vs/z/BW5oM39M3hf6ZWr8IVod2gFYD9ge2HE4K0TEOsBPyFXbjYDXk4MB8GNypfGv5Mr+6sAZNV19HwjgdeQQHAw8Xu7j+eQPimnALsCmwJfL8kTEnuQK8VHAK8pYfhgR27fcxxHkvyRsAhwfEWsB5wOXAluUxz8VOLXy+N5UdkGbOYyn4zcRcX9E/Cki3jrUwhHxQnLA38vgFdDDIuKhiLgiIg6JRQ/x2BF4Glg/Im6IiHsi4pTyumjiMJv9lc17gVuBfSNihYiYDHywLHtxWWZr4K8ppfsr6/2OXAjyr5ETh9kcx2wW+5bPxGvKbvkrtFnGz02ZzfHP5sD9TAe2Kc9D3TLPA/Yjf5b+rcyeCawGpIi4MiLujYjfRsQmndy/+p7ZHFk23wpcD7w9Im4v0wkRscog62wFLAX8fmBGSukWcpHqNYOsV2fE7wv9ZCIcotWJx4EPpJSe6WCdNcnnKvp5Sum+Mu/qgcaImAPMr7TVWQc4IaU0sO5tlbZ9gZWAXVNKj5Z5t1TaDy7rDuwa+I2I2Ar4FPCHynL/mVI6sTK2o4CLUkqfrczbD7g3IjZNKV1FrrTeyOCHUD0GfAz4M/nwi92A0yPinSmlnw+y3g/Iz9vZEbFszTJfJ78h/IMcxC8BawMHlfb1yBXjT5Kryk+Q31jOiohNUkpPDXL/ag6z2UfZTCktiIj/Q/4C8UTp+0Fgh5TSwC7I08h/pam6v9KmicFsjm82f0z+QXgv+cv3V4CXAztXlvFzU2A2xzubVR8gF20WO3dORHycnLcVgGuB7VJK80vzeuXffytjuAf4KHBeRGxY+XxVs5nNkWVzPWBD8ufd3sDS5M+70yPi9ansatNiGjA3pfREy/z76ey76Fi8L/SNJa3Ac1WHYYP8JepPwI0R8XvyBn5qSunhDvv5FnBMROwCnA2cllK6srRtDlxWCVur6cA3WuZdQD4Gv6r15FtbALMiYjaLW5/8fFxQ+q+VUrq3jH/ApRGxGnAIeTfExUTE+8m7Ce45RN9fq9y8KiKeBE6MiE+nlGaT3+ymAAemlM4vfe9J3tVxB+AXg/WvxjCbz+l5NiNiKeA/yF8+P0zeZfZD5A+7LVJKDww2Lk0oZvM5Xc1mWe97lZtXR8Td5B9/G6eUrivL+LkpMJtVXc/mgIhYBngPcFxKaUGbRX4I/Ir8g/0Q4JSImJVSeprnLm7z+ZTSGaW/95LPZbcXcMxQ969GMJvPGXY2yflYBtgrpXQX5IsFkPdc2rT82xWjfV/oN0vaIVpz2sxbSNk1reLZXZ1LQN8AvJl8HOOBwM0RsVEnd5xSOo68kf8nsDFwSUR8upM+2nXbcrv18U0i/wV+Rsu0AYtWYkfi4tJPne3Ih2k8GRHzee4EkL+IiLOH6Dd47mSQ95Z/rxtYoLzZPUD+i6UmBrPZX9l8E/BPwJ4ppT+nlC4j/8Uykb/YQv6xuFpL36tV2jQxmM3xy2bdOgyxnp+bSyaz2Zts7ga8EDihXWNK6dGU0s0ppXOBt5cx7lKa22VzHnmvPbM5cZjNkWXzXmDOQHGnuL78W5eP+4DlyiGRVasx+u+iI/nM7gtLWoGnnQepnBk/8rHui7yYKaWFKaULU0qfJ1cpHyGfqRvyX7aXGs4dpZTuSin9e0rp7cCR5B9MAFcAr4yIlWpWvYF8QtSqbah8QNS4nLxr9+0ppVtapnZV1k7M4LkPqnYOJv+IHAj5VmX+R8gnYh2sX3gulBeWfzccWKA8Ty8C7uxsyGoYszkyY5HN5cu/CwdWKrvGLuS5z40/A5tFxKqVvrcnn7Sva39lUV8wmyMzVDbr1mGI9fzc1ACzOTKdZPP9wB9SSsPJ0sAP+mXKv5eQTwhbzeZk8l61ZnNiM5tDuxBYISLWrMx7Wfm3Lh+XAAvI3z8BiIj1yYd7XdTBfbczks/svmCBB/5IPqHhNhHxCuDEamNEzIqIQyPiVRGxNvnKMavz3MZ+B/lEhptGxNSImNLuTiLi2Ih4Y0S8JCJeSd4QB/o4iXyM/BkR8dqIWC8i3hYRs0r714D3RsQHImKDcnzvO8gnwRrMMWWs/xURW5Z+3xj5hFVTyri2iXwSxtqTokbEeyPiXRExvUyfJn/AfaeyzCL9pJTuSSldMzBVHutdKaXbyzqvi4h/iYjNyvOyZxnzf6dy4taUj9v8HXBsRMwsr9GPgdupnFBLE5LZ7FE2ybsJzwF+FBGbRMSGETEw5t+UZX5N/qvjySXDO5DPBXJcSunJIR6/ms1sdiGbEbFRRBwWEVtExLoRsTP5ZOgXk3ff93NTQzGbXchmZf4G5MsoH9+m340j4uCIeGVErB0R25AP7Xiccq6esifdCcCREfFPEbEh+RLPywA/HeLxq9nM5hDZJH+/vAb4cfmM24KctT9RzkfU2k9K6SHy3krfjIg3lHV+TP7MPL/yvGwQETPI5+VZJiJmlGlyaR/yfaFRUh9cymssJ4a4bF2b+S8ETiGfXOlu8l+wn71sHfkM4b8j79r8FPmScR+rrL8Cebe0Rxn8snX/Qf4x9FTp6ycsennhdckfBI+Rf1j9BXhtpb31snX7VdoGLlu3U5v7nV4Z35Pk6uwxwFKl/U1l3ZmDPKfvI+8i92QZ3yWtj3OoftqNEXh16etR8km3ridfvm7ZlnVXIr8RPkK+rN8ZwLq93tacOpvM5mL327fZLPNnAmeRT+T6OPkvITu0LLMe+Yvrk+S/Tn0DWLrX25pTZ5PZXOx+e5LNkqfzS+aeAm4mn2Bypco6fm4uQZPZXOx+e/q5CXyVfPLWxT7nyHvhDDy3T5PPq3MSsEHLcssA3yzLPU4+T8qmvd7WnDqbzOZi9zvqbJblXgycznOXW/8xMLXSvlg/5L3O/5382TmnrL96S7//y+KXmk/AtNI+5PtCk6YoD0qSJEmSJEkN5SFakiRJkiRJDWeBR5IkSZIkqeEs8EiSJEmSJDWcBR5JkiRJkqSGa3yBJyL2i4jZvR7HeImIEyPizF6PA6BcxjVFxKt6PRb1H7PZO2ZTgzGbvWM2NRiz2TtmU4Mxm71jNjvX9wWeiNgtIv4YEY9GxJyIuDoijoyIF/V6bO2MVSDcmEcuIj4aETdExNyIuCcivhsRK/Z6XBON2TSbnYiIc8vz1jpd2+uxTTRm02x2ys/N8WE2zWYnIuIdEfGXyvZyZUTs2+txTURm02x2ot+/0/Z1gScijgROAa4EdgI2Bj4KvAQ4sIdDm3AiYulej2EsRMRewFHAkcBGwHuANwPH9HJcE43ZHD8TJZvAbsDqlWld4Angv3s4pgnHbI6fiZJNPzfHh9kcPxMlm8DDwBeBmcCmwI+AEyLizT0d1QRjNsfPBMpmf3+nTSn15QRsBSTg4zXtK5V/9wNmt7R9ELgFmFf+fX+b9puAp4CHgN8Bkyvt+wPXlfabgI8Bk4ZaHzi8jLk6bVsz/knAZ4G7gaeBq4G3Vtpb+zm3zD8ROJP8xvM34BHyG/7ylXUDOAS4FZhb+t670r5u6XNP4I9lmYNqxjkF+BJwZxnnbcC/tPTzqnJ7KeAE4PbS581lHNXnbhPgbOBxYDbwV+ANpW1p4NvA38t93Q18pcPt5ljgvJZ5XwCu6fU2PVEms2k2R5LNNuN/NzAfWKvX2/REmcym2RxJNvFz02yazb7MZs1juBz4cq+36YkymU2zORbZpM++0/Z8AIM8UceUF2TpIZZbJHDA24BngIOAlwEfKbd3Lu2vKi/Au4F1gM1KoCaX9vcD9wJvJ1dudwbuG9ggB1sfWBH4f8AfgGllmlIz7o+VjW6vMs4jgAXAjNK+ZdmYdyj9rFwJ3GPA8eS/tL0ReBQ4tNL3kcCNwJvKY9gLmAO8pSUod1Qe54trxvlT4B5gd2A94A3Ae2oCt3R5HFuWtneWsb230t/VwMnAdOCl5fXaurR9ooTsdcDawGuA/SvrHg6kIbaHPcrzM7PcXhu4Ajiq19v0RJkwm2ZzBNlsM/5zgTN7vT1PpAmzaTb93OzLCbNpNkf5uUn+Mb1deew79nqbnigTZtNsTsDvtD0fwCBP1G+Avw5judbAXQj8sGWZE4ELyv93Kxvs82r6uwvYp2XevwLXDXP9E4fzApOroZ9rs3GcXP6/yMbc0v/dwFKVeccDZ5X/r0CuZs5qWe9bwG9a+v7EEGPcoCz3ppr2tmNsWeYrA2Mrtx8H9q1Z9tvkamvUtB8E3DCM5/bD5Gr6M2V8J9X16dT5ZDbNZpv2YWWzsvzLyvjeOtx1nIb1vJpNs9na7udmH0xm02y2afgf+MwAABLNSURBVB9uNl9ALkA8Q96T471DreM0/Mlsms027Y3/TtvP5+CJEa63ETl0VReQj6eEXO28E7g9In4SEftGxPMAImJVYC3gPyJi9sBE3mjWH2r94YqI5wNrDDHOwVyXUlpQuf13YOAkYBsDywL/0/IYDqw8hgF/GeJ+NgcWAucMY0wARMSHygnhHiz3+zFydXTAN4AflBOZHRYR0yttJwIzgJvKCR7fEhHPbqMppWNTStXl293/68m7Iv4z8EryG+S25N3NNTbMZj2zOTwDf7n6dQfraGhms57ZrL9/Pze7z2zWM5uDe6L0syVwGPCNiNhuuI9BQzKb9czm8PTdd9p+LvDcBKwfEVPGqL8EkFJ6gvwF5p3k6umhwA0RsQbPPR8fIr/wA9MrgJcPY/0xG+cQnmmzzsDYB/7dmUUfw8vJu9dVzRnhGNuKiHeRK7cnknf1mwEcRz6uMg80pcPJbwpnkHeJuyoiDihtl5OrtIeWx/Fj4A/V0A3DF4GfppR+kFK6OqV0OvAZ4JCImDyax6dnmc16ZnPosUwB9gV+lFKaP9LHpLbMZj2zWc/Pze4zm/XM5iBSSgtTSreklK5MKX2dfDLgz4zmsWkRZrOe2Rx6LH35nbafCzz/Rd7966B2jRGxUs161wOvbZm3DfkkVgCklOanlP6YUjqUfFb6FYCdUkr3k6uT65c300WmodYvzfPIJ3+qlVJ6vNzPYOOcV/4dtK82riOfMGqdNo/hzg77upK8jbxhmMtvA1xcKp+Xl+estYpLSunmlNK3U0pvIZ8k632VtidSSj9PKR0IvAX4P+RjJ4drefKxpVULGHmFXoszm5nZ7CybA3YFppb+NbbMZmY2/dzsN2YzM5sj+9ysmgQsM8o+9ByzmZnNCfSdtm//MpNSujgijgK+FhEvBk4ln3zpJcB7yWcrb7f78NeAUyLiMuD35BM/vZu8yzERsRN5Izgf+Ad5Y3oeOagAnwe+ExGPko/LXJpcQV0zpfTlYax/B7BjRGxIvrzhYyml1growDiPiIibgcuAvYFZ5b4AHiAf27hDRNwBPJVSemwYz9sTEXE0cHRERBnniuRLLC5MKX1/qD4qfd0UEf9N3sXto+Qz978YWDel9J9tVrkJ2C8idiS/PnsAryefeZ2IWA44mvzXhzuA1SghLe0fJ+/idiW5arwX+RjKe0r7QeSTjw2229yvgI9HxF9Kvy8F/o18nGrfVFabzGyaTUaWzQEfAM5OKd023Mer4TGbZhM/N/uS2TSbjCCbEXFY6e82clHnzcA+5BP6agyYTbPJRPxOm/rgRECDTcA7yCeDeoy8i9c15LN2r1ra92Pxy9Z9iPyCP0PLZevIL/A55DDMLf3t37L+nuSN6ynyxnIBsMdw1gdWJQf9CRj2Zevmkc/2vWvLMu8j75a3gJbL1rUsdziVy5mS/+r2EZ6rrj5IPpZz+9K+LkOcrKrS1zLAUeSTdD1NvhTeQe36Ie8ad0J5zh4t//8ccEel/b/IYXuaXFX+PvD80v7+8rw/QQ7aecBrWh5nGmK8k8lvmjeX1+du8m57L+z1tjzRJrNpNlseZxrGmNcjH2f9zl5vvxN5Mptms+VxpiHG6+fmOE1m02y2PM40xHi/XMnlP4CLgD17vR1PxMlsms2Wx5mGMea+/U4bZYCSJEmSJElqqH4+B48kSZIkSZKGwQKPJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ034sukR4RnZ9aS7qGU0qq9HkQrsymZTalPmU2pP5lNqQ+llKLTddyDRxq5O3s9AEltmU2pP5lNqT+ZTWmCsMAjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkNZ4FHkiRJkiSp4SzwSJIkSZIkNZwFHkmSJEmSpIazwCNJkiRJktRwFngkSZIkSZIazgKPJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ1ngUeSJEmSJKnhLPBIkiRJkiQ1nAUeSZIkSZKkhrPAI0mSJEmS1HAWeCRJkiRJkhrOAo8kSZIkSVLDWeCRJEmSJElqOAs8kiRJkiRJDWeBR5IkSZIkqeEs8EiSJEmSJDWcBR5JkiRJkqSGs8AjSZIkSZLUcJN7PQCNzKc+9anatp122qm2bf/9969tu+WWW0Y1Jkkjt/vuu9e2HXLIIbVtr371q7sxHGnCOfTQQ2vbjjzyyNq2tddeu7btnnvuGdWYpCXJqquu2nb+xz72sdp17r///tq2Y445ZtRjkjQym222WW3bscceW9v24Q9/uLbtqquuGtWYlLkHjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkNZ4FHkiRJkiSp4SzwSJIkSZIkNZxX0Wqot7/97bVtW2yxRW3bGmusUdvmVbSk7poxY0Zt2/HHH1/bFhHdGI40Ia2yyipt5x944IG166SUujUcScUOO+zQdv5gV7g77rjjujUcSaPwox/9qLZt8803r217z3veU9t28MEHj2pMytyDR5IkSZIkqeEs8EiSJEmSJDWcBR5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsN5mfQ+t9Zaa3U0X1L/2nrrrWvbVlpppdq2M844oxvDkSak5Zdfvu38Nddcc5xHIqnq5JNPbjt/sMstS5pYzjnnnF4PYcJzDx5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkN52XS+9xqq63Wdv6LXvSicR6JpNFadtllR7Te6aefPsYjkVR100031bbNnj17HEciqWrGjBm9HoK0xJo2bdqI2gYzebLlh25zDx5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkN53XKJGmc7L333iNa79JLLx3jkUhLnkmT6v+mteGGG9a2rbjiirVtjz766KjGJGlwW265Za+HIC2x1lxzzdq21VdffUR93njjjSMdjobJPXgkSZIkSZIazgKPJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ1ngUeSJEmSJKnhLPBIkiRJkiQ1nJdJb6iIqG1LKY3jSCRVbbrpprVtm2yySW3bJZdcUtt2yy23jGpMkmDhwoW1bRdddFFt28MPP9yN4UiSNCHdddddtW0PPfTQOI5kyeQePJIkSZIkSQ1ngUeSJEmSJKnhLPBIkiRJkiQ1nAUeSZIkSZKkhrPAI0mSJEmS1HAWeCRJkiRJkhrOy6Q3lJdCl/rTLrvsUts2eXL9W+68efNq2+bPnz+qMUka3N/+9rfatrlz547jSCRJ6g877rjjiNY777zzatu8THr3uQePJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ1ngUeSJEmSJKnhLPBIkiRJkiQ1nAUeSZIkSZKkhvMy6X1ujTXW6PUQJHVgpJk96aSTxngkkiRJ0sh85jOfGdF6Z5999hiPRJ1wDx5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkN52XS+9wuu+wypv0tXLhwTPuTtKidd965tm3BggW1bb/85S+7MRxpibPFFlu0nR8RtesM1iZJ0pJo6aWXrm1LKdW23Xrrrd0YjobJPXgkSZIkSZIazgKPJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ1ngUeSJEmSJKnhLPBIkiRJkiQ1nJdJ7wOTJ9e/DLNmzeq4v/POO6+27YILLui4P0mLmzlzZtv5a6yxRu06g2XzgQceGPWYJMFll13Wdv5gl3QdrE2SpIms7jvtUkstVbvO/Pnza9v8vdlb7sEjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkN51W0+sCkSfV1tg022KDj/k477bTRDEfSMHz6059uOz8iate5/vrruzUcScUee+zR6yFIamP69Olt5w/2uXnttdd2aziSio022qjtfL/TNpN78EiSJEmSJDWcBR5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4L5PeBz74wQ/2egiSOrTKKqv0egiS2vjVr37Vdv5XvvKVcR6JpKr58+d3vM6FF17YhZFIqjrggAM6XufUU0/twkg0FtyDR5IkSZIkqeEs8EiSJEmSJDWcBR5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsN5mfQ+EBFj2t8tt9wypv1JGr4HHnigtu24444bx5FIS6Y5c+a0nT/YZ+1Yfw5LWlzd99OUUu06U6ZM6dZwpCXKeuutV9u2xRZbdNzfmWeeOZrhqIvcg0eSJEmSJKnhLPBIkiRJkiQ1nAUeSZIkSZKkhrPAI0mSJEmS1HAWeCRJkiRJkhrOAo8kSZIkSVLDeZn0Ceiss87q9RCkCWHNNdesbZs+fXrb+XfeeWftOtdee+2oxyRpZAa7FPNgbZLGRt1nakTUrrPvvvvWtr3//e8f9ZikJcWUKVNq25ZbbrlxHIm6zT14JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkNZ4FHkiRJkiSp4SzwSJIkSZIkNZyXSe8D66yzTsfr/O53v6ttW7BgwWiGI6mYNm1abdvUqVPbzv/qV7/areFI6pK11lqrtm355ZevbXvyySe7MRxpQnr00Ufbzk8pjfNIpCXPXnvt1fE68+fPr22bN2/eaIajLnIPHkmSJEmSpIazwCNJkiRJktRwFngkSZIkSZIazgKPJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ3nZdLHyWCXQt9nn3067u+Xv/xlbdvChQs77k/S4vbcc8/atrrLus6ZM6dbw5HUJTNnzqxtW3nllWvbvEy6NHx+Pkq9s8IKK3S8zmCfcXPnzh3NcNRF7sEjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkNZ4FHkiRJkiSp4bxM+jjZZpttatumTp1a23brrbe2nf+zn/1s1GOSNLh99tmnti0i2s6/9NJLuzUcSaNQl1lJkia6LbfcsuN1TjnllNq2O++8czTDURe5B48kSZIkSVLDWeCRJEmSJElqOAs8kiRJkiRJDWeBR5IkSZIkqeEs8EiSJEmSJDWcBR5JkiRJkqSG8zLp42SnnXYa0XqPPPJIR/MljZ1lllmmtu22225rO/+KK67o1nAkDcO8efPazn/44Ydr11l55ZVr26ZNm1bbds899wx/YJIk9chmm23W8TqXXXZZbdv8+fNHMxx1kXvwSJIkSZIkNZwFHkmSJEmSpIazwCNJkiRJktRwFngkSZIkSZIazgKPJEmSJElSw1ngkSRJkiRJajgvkz5O7rvvvhGtd+21147xSCSNhfPPP7/t/AULFozzSCRV3X///W3nn3POObXr7L777t0ajqQhXHTRRbVts2bNqm3be++9a9tOPvnkUY1JkprKPXgkSZIkSZIazgKPJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ1ngUeSJEmSJKnhLPBIkiRJkiQ1nJdJHyff+973atv22Wef2rYbb7yxG8ORNEpvectb2s6fPLn+bXX+/PndGo6kLrnvvvt6PQRpQtt3331r2771rW/Vti211FLdGI40IQ32m3KllVZqO//666/v1nDURe7BI0mSJEmS1HAWeCRJkiRJkhrOAo8kSZIkSVLDWeCRJEmSJElqOAs8kiRJkiRJDRcppZGtGDGyFaWJ47KU0qt6PYhWZlMym1KfMptSfzKbUh9KKUWn67gHjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkNZ4FHkiRJkiSp4SzwSJIkSZIkNZwFHkmSJEmSpIazwCNJkiRJktRwFngkSZIkSZIazgKPJEmSJElSw1ngkSRJkiRJajgLPJIkSZIkSQ1ngUeSJEmSJKnhLPBIkiRJkiQ1nAUeSZIkSZKkhrPAI0mSJEmS1HAWeCRJkiRJkhrOAo8kSZIkSVLDWeCRJEmSJElqOAs8kiRJkiRJDWeBR5IkSZIkqeEs8EiSJEmSJDWcBR5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGs4CjyRJkiRJUsNZ4JEkSZIkSWo4CzySJEmSJEkNN3kU6z4E3DlWA5EaaJ1eD6CG2dSSzmxK/clsSv3JbEr9Z0S5jJTSWA9EkiRJkiRJ48hDtCRJkiRJkhrOAo8kSZIkSVLDWeCRJEmSJElqOAs8kiRJkiRJDWeBR5IkSZIkqeEs8EiSJEmSJDWcBR5JkiRJkqSGs8AjSZIkSZLUcBZ4JEmSJEmSGu7/A+rz5kO5PZEKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(imgs_max[i].reshape(28, 28))\n",
    "    plt.title('Model prediction: {} (p={:.2f}) \\n Label: {} \\n Trust score: {:.3f}' \\\n",
    "              '\\n Closest other class: {}'.format(pred_max[i].argmax(),pred_max[i].max(),\n",
    "                                                  label_max[i], score_max[i], closest_max[i]), fontsize=14)\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### High model confidence, low trust score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where trust scores really matter is when the predicted model confidence is relatively high (e.g. $p\\in[0.8, 0.9]$) but the corresponding trust score is low, this can indicate samples for which the model is overconfident.The trust score provides a diagnostic for finding these examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAETCAYAAABN+bF9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd7wU1fnH8e8DCEgRVLAX7GILih1UEgR7YosVFaPYYqLkp0YliSWixhY1JBoVgyWxYIsaBbsRa8RYETUKWMAEC01pwvn9cc7VYdkzd3fv3t2dvZ/36zUvuPNMObM7z87sszNzzDknAAAAAAAAZFerajcAAAAAAAAATUOBBwAAAAAAIOMo8AAAAAAAAGQcBR4AAAAAAICMo8ADAAAAAACQcRR4AAAAAAAAMq4mCjxm1s3MnJn1K2Kec83szWZsVtHM7E0zOzfx92QzO62Jy3zKzEY0uXEVZGbLm9l/zWy9arelOZjZT83sgWq3oxLIzdRlZjE3Nwi52aXabWkOZnapmf2h2u2oBHIzdZnkZo0xs9Fm9n/VbkclkJupyyQ3awy52eg85GaNagG5+ZKZHVDKvI0WeMxsVEiGkXlivwuxB0tZeQuwjaQ/FTKhmQ02szl5QvtLOqusrSq8PS4ybNPI7GdLesg5934zt/EAM5tgZvPDv/sVMM9uZva8mc02s8/M7O9mtmEiPiqyzV8lFnODpN5mtlNzbFehyM0myWxuJoWTlU/Ce92tgFkulPQn59zMZm7XSWY2yczmmdn4QnLFzA4zs1fN7Gsz+9TMbjWzVRLxpyK5+VZiMZdIOsrM1m2O7SoUudkkmc7NyD56QgGzZjk3f2xmL5vZDDP7Kkx7VM5izpc0rNon4uRmk5Cbzdu+sudmmGY5M7vazKaG8+X/mNlBiUnIzezLem72N7PnwnezT8P73aaAWbOem419j/2tpIvNrOgLcgqd4SNJB5lZx0Sj2kg6UtKHxa60lpnZMuValnNuunPu6yYu4wvn3OxytakId0haNWe4VdIHkl6OzWRmHSQdK2mpD+hyMrMdQhv/KqlX+He0mW2XMs86kv4u6RlJW0raVdKykh5KTHaKlt7uDyTd2TCBc26+pL9J+nn5tqhk5GYJMp6bSX+R9GohE5rZmpL2DfM0GzM7WNJV8gfeLSU9J+lhM1srZZ4+km6RdJOkTUM7N5HP6wb7a8m87CFptpbMzemSHpF0Ytk2qHTkZgnqJDeHaMl99aa0iesgNz+XdIGk7SVtIb8dI81sz4YJnHNvyB9LB5Vxk0pFbpaA3Gw+zZWb4f1/VNIGkg6StJGkwZImNUxDblYHuemZ2ffkv4c9Kr/vHyzph5IubmS+rOdmId9jH5LUWdIeRTfcOZc6SBol6UFJr0g6OjH+R/IfEDdJejAxvpWkX8sn6XxJb0j6Uc4yt5E0XtI8Sf+WtJckJ6lfYppNJP1D/gT+f5Juk7RKIn6upDdT2t0jLPMwSePCuiZKGpiYpl+YZk9JL0laIGnvENsn0cZJkoZLapuYdyX5YsFcSVMk/UTSm5LOTUwzWdJpib+7SLpG0rSw3Lfld+SGdiSHc8M8T0kakVjG8uE1/zKs+zFJmybigyXNkdQ/tOcrSU9KWqex97qR/aCDpBmSzm5kugMlfSHJ8rzOe8t/GZ0XXtveTWjPHZIezRn3mKTbGmnbIkmtE+O+H9rWLTJPnxDfMWf8zvL7d4emvK5NfE9GidxssbkpX4x8XNIP0vbhxPSnSfp3zriGNu0j6d2w7U9KWrcJ++WLkq7PGfeepIsaaduUnHFHS5qTMs/hkr6RtGbO+CMlfVytvCQ3W3ZuhnYcWOQ8dZWbYZpXcpcr6TeSxpGb5Ca5WZnclHScfPGmbSPrJzfJzadUhdyUL5zk5tg+YZ2d6zg3C/oeK+lGSbcW2+5iLvkZGXaqBj+Rr5q5nOlOkXS6pF9K2lzSvZLuMbNekmRmneQT6QNJW0s6U9JlyQWY2aqS/im/w2wrf6VFJ0l/L+EypUskXS1fHXs0LGP1nGl+J+lXkjaW9KKZ7SZfSRshX3n7iXxx4MLEPKMkrR/atq/8l4oesUaYmclX4naRf5M3kfQL+SR/TtKpkr7Wd78qXJZ/SRolaTv5D7xtwzxjzGzZxDTt5C+z+4mkHSR1lXRtoi09wqWOg2PtzeMgSR3ld7Q0O0ka78JemeMy+f1ia/n3/8FwxY/MbC0zm9PIcG1iWTvI/1KfNFbSjilt+5ekhZKONbPWZtZZ0lGS/uWc+ywyzxBJbznnnssZ/7KkNqEd1UZutrDcNLMt5d/HIyUtTps2YSflv/qunaRz5Ld9B0mt5fcLC+vaqYDcPDtM21ZSby2dm48oPTeflbSqme1jXjdJh2jJq+tyDZE0xjn3Uc74lyStbrXxDDBys4XlZnCV+VuA/2VmJxTw+tdNboZp+stfKfDPnPBLkrbNed2rhdwkN1tCbu4bpvuD+dtEJph/pkzu1SPkJrmZa5Qqk5vt5AsxSXMltZfPi5is52ah32Nfkn8vi1NERXV5+Rd8A0mryFdL12qIJ6b/RNJvcpbxlEL1Sb6aPENSp0R8kBIVVfn7QR/PWcbyYZpti6yoDkuMayVf0bsgp6J6QM68/5T065xx+8pXBU3ShmG+Pon42vJXh+StqEoaIP9FrGekvYOV5xcxJSqq4bV3knZOxLtIminp2MRynKSNEtMcHt4vC3+vLl9d3q+ICuZzku4tYLr7JN2UM67hdT48Ma5T2A8a2t1G/gMsbVgpMf8CSUfmrOdISfMbad9Okj6V/+V/sXzVfKXItF3kP9BOicS/kHRMsVXVcg0iN1tkbsoXWt9peG0Sr1VjV/C8Kum8PNsWe712DX8vq8Zzc4Uw7Wq5r0MY/xtJ7zTSvv0lzZIvwjr5A9+ykWkb3ucf5YktF2L9yU1yUxU+bsr/otxX/iT//+R/1fxVvedmeE3nhGnmKc+xUf72LSdpPXKT3BS52ey5GV6XefI/zvaWdID8OfBl5Ca5WQu5KWlgaO8g+e+Cq4fXxUk6tI5zs6DvsfK3qy2W1KaYfCrkAUaSJOfcl2Z2r3yVboakp5xzH4aCmCT/IK/wYjybM/s4+cvSJKmnpNedc8kHPD2fM31vSTtb/odArSdfzSrUt8t2zi02sxflq5lJuRXA3vKV7F8mxrWS3ylWkd+Gxcl2OOemmNnUlHZsKWmac+7tItqeq2G9yW2aaWZvaMltmu+ceyfx91RJbeU/tL5wzn0iXz0uiJltKl9p3KuAyZeV9N9ILNnuOcl2O+e+kfSfQttUCvMPtxop6Wb5SzA7y3+432lmP3DO5V4JMUj+fb8lssi58ttbVeRmi8vNq+Uvpb67yDYuq6V/JZHir9cmkh5zzs1V8+fmJpL+IP9AubHyvypdKunP8ge8XEPkLz3+R57Y3PAvufkdcrNCx03n3G8Tf75qZq0lDZN/Rk1MPeTmbPkvzp3kL9m/wswmO+ceT0xDbi6N3CQ3S1JgbraSv+1oiHNukaTxZraipN+b2ekufIMUuUluLqliuemce8R8D2B/lC/izZffp3dS+hXqWc/NQs2VL/a1ly/8FaTgAk9wo/z9eHPkq1fFcI1P8q1W8ifu+bp8ixUPmuKrnL9bSTpP0ug8005P/L+YbWpuybZ8E4kV/RTu4Dj5e1zHFDDtZ/KJXRTzD6ua0MhktzrnGno8+FTSyjnxlcP4mJ9K+so5d0ZivYPkt21H+QND0hBJdzvnvogsbwUtuT9UE7nZcnKzv6Q17bteahrOej41s98554ZF5kvLzejrZb63gIcbadOFzrkLwzoWqfjcPEvSS865S8Pfr5vvue4ZMzvbOfdxoj1t5W+tvD4UhnOtEP4lN79DbubXnMfNBi9KWs7MVnbOxd6HzOdm+IGk4aT5VTPrKd+jZrLAQ24ujdzMj9wsT25Ok7QwFHcavC3/XM1u+u79JzeXRm7mV/bcdM5dYWa/ly+EfCl/VdRF8rfXxWQ9Nwv9HruCpHk5hcpGFVvgeVz+kqJu8rfiLME5NytUyfpoyYN6X3335f1tSYPNrKNzrmFH3z5nUa/IP/NlinNuYZFtzLW9pCekb+9L3FbSXY3M84qkjZ1zeSt8ZjZRfufdVv7WpYYCxWopy/y3/P14PSNV1QXy9wmmeTusdweFe9tDFXtzNdNTxM2svaQjJF2d5wqXfP4tf4lcPtsrJKv5J+RvJn81jeSrvr0aWfasxP+fl78M8dLEuAEK70dEB/lETWr4e4kPIzPbVtL35O9VXYr553u0l99XagG5qRaTmwPlfyFpsI38yVA/+Qe/xfxbS/+aJMVfr4bX4mU1nptfSJJzboGZjZfPxeQJywBJaVccFZyb8pcvd1O8p77N5C+JfaORNlcKuakWk5v59JL/lXFGyjT1kpu5bW+XM24zSZ+kfJmuNHJT5KbqOzeflXSYmbVKnMNvKP/4geSzJ8lNcjOp4rkZriabGtZ1qPyP72nfsbKem4V+j91MpXzXjN275b6792uUlrznsbOk5VLip8p/ET9U/kPk/LBR3wvxTvJVyTvkHyg1QD4Zk/dEriZ/SeE98g94Wlf+4VLXKTxRW4XfE/mR/AOrNpLv5myepDXckvdEdsuZdzf5Lwjnhxd247CMSxLTPCz/BWIH+R3lCflKc+yeyFbyb+aEsPx1wrbvG+I7hrYMkP9A6+By7okMf98nv5PuJJ9o94dtXNZF7q3M3U4V8Qwe+duUFklaq7Fpw/Sbh+lXzLP+CWH7Ng3v/38ldSxkuXnWs6N85fjM8P6cFd6z7RLTXKTEvbXyvQ0tlv81YANJW8lflfRhbjsk3SDp3ZT1D5b0filtL9cgcrNF52ZsGSnT7SN/QtcmMW5weD1fSrxeT0p6XYme8IrcLw+WP4E4Vv4y36vC6792YpqbJd2cpx0nhn2qj/xD0cfnWf5jyul5ICd+rnLuqSc3yU1VIDflc2xI2P715HNgpqSr6jk35W9z2TXEe8o/32ShpBPy5MVIcpPcJDcrlptryu+/fwj7zG6SPpZ0KblJbqoGcjNMd3pYx6byz8pa0NDeOs7NRr/HJt6XXxfd7mITroCEbKXvuq1bEHbKfXPm2U6+GjVf0mvhTfo24cI0G8hXPhu6Z3tH/gOqbZEJd7h8NWxeWMYesR0xZ/6Bkp6Rr3LPkq/4nZyIrxx29rlhW49V493WdZV0vfwHzjz55DsoEb8m7KxOTey2rpGEa3htBhfw/j8t6aEik+F5ST/Ns/4fyifZ/PD+b1NKsiWWe6D8B8cC+Q+h/fPsm5Nzxh0i/2DlOeF9eEDSJjnTdA7xM1LWPVbSmU1pf1MHkZstOjcLea1ypmsdXo+9ctsk30vCe+F9f1rS+k3cN08Kr/F8+XzLfUDdU/L31SfH/UzSW+F9nSbfs8QaOdOsK1+kPShl3e9IOoTcJDcrnZuSdpf/VXG2/KX4b8j39JL6cERlPDflf0x5L7y2X4R959CcZbSX/0K9PblJbpKblcnNMM32YZ+ZK98N9/lashtucpPcbNi/qnJOK1+0mhHW80Lytavz3Gzse+zqIbZGsW1ueMp13TGzHvIfZNs45/J1o4ZmYma7y1c4N3HOLTKzfvIV1O4u3h15ZpjZZvKXhG7onJtZ7fZkDblZPWZ2gqQfO+f6h78Hyx/QO1W1YWViZnvJX+66hcv/fB6kIDerpwXk5k/le70bWO22ZBG5WT3kJtKQm9XTAnLzUkldnHPHFTtvsc/gARrlnBtjZn+UtIakKdVuTzNYTb5rO4o7yJrrJa1gZl3qdP/tKOloijvIoHrPzYXyv2gCWUNuArWp3nPzf5IuK2VGCjxoFs65q6vdhubinHuk2m0ASuF8TxoXVrsdzcU5d2e12wCUogXk5nXVbgNQCnITqE0tIDcvbXyq/Or2Fi0AAAAAAICWoqh+6gEAAAAAAFB7KPDUGDObbGanNXEZ/czMmVm3crULaOnITaA2kZtAbSI3gdpEbtY3CjxlZGY9wo6+dbXbUqzQ9pFm9oGZzQ3/XmRmy1a7bUBTZTk3k8ysvZm9Vg/bAkjZzk0za2Vm95vZh2Y2z8ymmdmtZrZ6tdsGNBW5CdSmLOdmEue0zYcCDxpsLKm1pBMlbSr/RP0j5bs7B1AbLpP0cbUbAeBbT0g6SNJGkg6QtK6ke6vaIgASuQnUOs5pmwkFngoys/XM7O9m9qmZfWVmr5jZ3nkm7RR+aZgTpj0tZzldzOw6M/ufmc02s6ebWvl0zo1xzg12zo11zn3gnPuHpOHyB0WgrtVybiaW/SNJ35fUpEtqgSyp5dx0zi12zl3pnHvBOTfFOfecpIslbWNm7ZuybKDWkZtAbarl3Ewsm3PaZkSBp7I6SXpY0gBJ35N0t6R7zGzjnOl+IeltSVtJOkfShWa2vySZmUn6h6TVJe0taUtJ/5T0hJmtGltxuPzt3CLbu5ykL4ucB8iims5NM1tD0jWSDpM0t9iNAzKspnMzZ/oVJB0u6UXn3LxC5wMyitwEalNN5ybntBXgnGMo0yCphyQnaesi5nlB0q8Sf0+W9GjONDdIGhf+/wNJcyQtmzPNq5LOCP/vF9rRLRGfKOnkItq1tqTPJP2i2q8rA0NThyznpvytk09L+r9St4WBoVaHLOdmYrrfSfoqzP98chkMDFkdyE0GhtocspybnNNWZmgjVIyZdZSvkO4taVVJy0hqL+n1nEmfz/P3/uH/vSV1kDTdF1e/1V7SerF1O+dyq7Zp7VxZ0hhJj0r6faHzAVlV47l5tqQFkq5oZDqg7tR4bja4VNJI+R9GzpF0q5nt4cLZK1CPyE2gNtV4bnJOWwEUeCrrMkm7y99v+J6kryXdLKltEctoJem/knbKE5vV1Aaa2SryD6Z7U9IRHATRQtRybvYPy1yYc5B9wczucM4d3oRlA7WulnNTkuSc+0z+itd3zextSR9J6ivpmaYuG6hh5CZQm2o5NzmnrQAKPJXVV9LNzrm7Jd89nHwV9N2c6bbP8/fb4f+vSFpZ0mLn3AflbFy4p/JJSW9JOtQ59005lw/UsFrOzaMldUz8vZqksfLPE3i2jOsBalEt52Y+Dc82bNfM6wGqjdwEalMt5ybntBVAgad5bGhmucWRifKJtZ+Z/V3SQvnL5/I9zX97MztL0l3y9zceKb/jS9Jj8gnwdzM7Iyx3FflK7WPOuby/SpjZREkjnHMjIvHVJD0laaqkUyV1S1RWpzvnFjWyzUAWZC43nXOTcqafE/77vnOO7iVRLzKXm2a2g/zDKcdJmiF/Av1b+WcbjGt8k4FMIDeB2pS53OSctjIo8DSPv+YZt7n808pHyl8a+qWkK5U/4a6QtIWkYfIPh/uNc+4uSXLOOTPbU9IFkq6XtJL8JXTPyl9+F7ORpG4p8YGSNgjDhzmxdeQPikDWZTE3gZYgi7k5V9KBks6X/0Vymvzz6w529NSD+kFuArUpi7mJCjAesQIAAAAAAJBtrRqfBAAAAAAAALWMAg8AAAAAAEDGUeABAAAAAADIOAo8AAAAAAAAGUeBBwAAAAAAIOMo8AAAAAAAAGRcXRR4zOxcM3ONDD2acf0Xm9nLzbX8WmFmW5rZODObZ2Yfm9lZjUy/e8r7sU9iuvPM7Hkz+9rM5uVZzmpm9oiZTTOz+Wb2oZldZWadm2M7UT7kZmUUm5thnm3N7Akzm2lmn5nZn8ysQ840O5jZk2Y2w8y+DHm4VSK+q5k9YGafmtlXZvaqmR3RHNuI8iEvK6O5jplh2r3N7EUzmxty8+E8y3rRzOaY2VQzu8DM6uKcr56Rm5VR4jGzZzjefW5ms8zsWTPrn4i3j7xfgxPTdDKzm8zsDTP7xszGNNMmoszIzcpojtxMTHe0mb0elj3dzK5PxOrqfLZeDvaXSVo1Mbwj6fKccR/lzmRmbSvYxoor5/aZ2QqSHpU0WVJvSadJ+pWZ/TRltie05HuwqqQrJH0ZltVgGUl3ShoRWc4iSXdJ2kvShpKOlbS3pD+WtjWoIHIzj2rnppmtHeZ5S9I28rnVW9J1iWmWl/SwpPclbSepr3zuPmJmy4bJ+kp6RdL+kraQdKOkUWZ2QLm2D82CvMyj2nmpAo+ZZnaIpJskjZT0PUl9JN2ciG8j6QFJf5fUS9IgSQdLOr8sG4fmRG7mUe3cNDOTPx5KUr8w33hJD5rZmjmTH6El3687ErE2kr6WdLWksU3cFFQWuZlHVnLTzM6Q9FtJF0naVFJ/SQ8lFlVf57POubobJL0p6dw842+XLxT8WtJUSR+F8Z9KOjln2hckXZb4++Cw3LmSvpD0pKQVJZ0gyeUMh0TatY6kB+VP1r6SNEHS/on4mqGNn8sfAMZL6puI/0zSB5IWSHpX0lGJWPuw7uMk3R+Wf0GIbS5pjKQ5kv4r6VZJ3Yt8TYeGdi2TGHeBpElFLKOVfNJeHYkPkjSvwGWdUcy6GWpjIDdrIzcl/Ty8zpYYt42kxZLWDH/3De1eNTFNzzBus5Rl3y/pr9Xe1xiK2ofIyxrIyzzLWOqYKalteP2PSJnvCknP54z7saTZkpat9v7GUNR+RG7WQG5KWiO0aYfEuA5h3N457d67wHbcIGlMtfcxhtIGcjNTudld0jxJOxXZnsyez9bLFTzF2E3SupIGSNqjkBnCr91/lXSt/BecXeSTQ/K/oo2Q9Jq+q+DeF1nUdZJM0s7yiXCapFlhHctJekbSKpJ+KF89vChMLzM7VL56fImkzUJbbjSzATnrOF/SPWH514fq5T8l/Uu+qrmbpG6S7k5sX8Nl4dunvAw7SHrKObcwMW6spB5mtmrKfEkDJa2txFUCpTCzNSTtK+nppiwHNYfcrFxutpM034UjWDA3bFOf8Pdb8icIx5pZWzNrL3/13PuS3ktpz3JhPtQH8rK2jpnbS1pZkguXkE8zs4fNbPPENO3kT2aT5krqJH/FD+oDuVm53Jwmf+w7ysw6mlkbScfLfwl/MWfaa8PtHy+Z2THhCgO0LORmbeXmHpLmS1rPzCaGW79Gh9c8TWbPZ9tUuwFVMEvScTk7T2NWl/8l7S7n3Kdh3BsNQTP7StI3iVjM2pJGOuca5v0gETtKUldJ+zrnZoRx/0nETwvzXhv+vsLMtpX0Sy15u9MtzrlRibZdIuk559yvE+MGS5pmZls4516Xr7a+I38CGLOKfFU56b+J2LSUeRscJ+kF51zucgpiZvdK2l2+gnxPWB7qB7lZudx8TNLFZjZU/qRhOUkXhtiqkuSc+9LMvi/pXknnhtj7kgY45+bna4iZHShpR0mnprQX2UJe1tYxc93w72/lf+38WNIpkp42s42cc9PlT4hPMLOD5E+wV5P0qzBfocUl1D5ys0K56ZxbZGY/kP9SPVv+atfpknYLOSf5xwn8StJT8ldGDJA/vi4v/6UZLQe5WVu5ua78d8fT5a9Qmi1fpHrMzDZ3zuV7/mumz2db4hU8rxeZcJKvSD4j6Z1Q8TvOzFYsYd1XSrrA/MOfzjezXonYlpLGJxIu18aSns0ZN07SJjnjch/A1VvSAPMPWpxjZnP0XTKvJ0nOuXHOuY2dc68Vu0GFMrNVJO0j6frGpk1xkvzrtL98ZfnSMjQNtYPcrFBuOuf+LekYSWfJH2ynyp9IzJA/OMr8Q8xvlH8uyLbyt2y9Lenv4WqeJYRi0ChJJzjnXi1ne1FV5GVtHTMbztvOcc7d55x7WT6XF0g6LLTvfknD5G8BmS+ft/eH+RY3V5tRceRmhXLTzFpL+rN8QbWP/HPpHpR0r5mtFNa90Dk33Dn3rHPu3865S+SvjjijnG1BJpCbNZSb8sfNtpJOdM496px7QdKhknrIX22Uu8zMn8+2xALPV3nGLVa4PC1hmYb/hCT9vqQ95e9lPFHSe2bWs5gVO+f+JL+j3yKfLC+Z2ZnFLCPfYnP+zt2+VvJVzV45wwZashrbmE/lLwtPWjkRa8zR8l8k72hswhjn3DTn3ETn3L2Sfirp54nkRfaRmxXMTefcKOfcSvL3L68o/wtjV333a8+R8r+YDHHOjXfOPS9/QNxY/iHn3woHwwckne6cu7GItqP2kZe1dcxs+AVzQsMI59wC+avr1kqMu1hSlzCum757CGXy11xkG7lZudzcXdKukg51zj3vnBsvf4Wdkz9WxrwoqbuZdSmifcg+crO2cjPfcfNzSf9T4rgp1c/5bEss8OQzXYnLls2so/xO+S3n3OJQlT9HvlL5pfxDCyX/y1nrQlbknPvQOXetc+5AScP13W1G/5a0lZl1jcw6Ud89G6NBXyV21ohX5J8WPsk595+cYU4hbQ6el9TPzJZJjBsgabJzLvVS83D/8bGS/uacy/ehV4qGfbeun04PcrMAJeemJDnnPg3rO1T+suInQ6iDpMXOLfGcnsXyB81vjx3mu6J8UNIvnXPXFNFuZBd52bjmOma+JOkbSRslpm8j/2DNKckJnTc1XH5+qHwRqKRbpJEZ5GbjSsnNDuHfb6+AC8fGxUr/LtVL/taUYtqH+kRuNq65crPhqqTkcbOrpJWUOG7W0/ksBR7vCfmHM/U1s83kL8v6lpntZGZnmdnWZraWpP3kk7Rhh58s/+CmLcysm0W6jDOzEWY20MzWMbOt5HfahmXcLH9P4H1m1sfM1jWz/cxspxC/VNIx4ZK9DczsF/JJf0kj23ZVaOvfzGybsNyBZjayoZ1huyeaWdrDF2+W/2I30sw2NX9v/y+UuK/YzA4Oy+mWM29/+fsf896eZWZrh0sI1/J/Wq8wdAjxfc3siLDeHmb2Q/n7mp92zn3cyPYj28jNZspNMzvFzLY0s43M7BT5nndOTxyMx0pa2cyuCtNsLv/6f6XwgHPzD957UNLvJd1tZquEIfczAPWFvKzSMTP86jhS0nAz29XMNpL0R/kHK98WlruMmf3CzDYLw/nyzxH4WU7BFvWH3Gye3HxG/tj3FzPbPBwTG9rzUJhnfzP7SVjm+mZ2gnxPSn9yzi1KLHtT8+e8K0jqbP58d4tGth3ZR25WKTedfwbQWEkjzGz78PrfJGmSpEfCcuvrfNbVQFde5R7USNd1ecYvL2m0pJmSPpL/5ezbruvknxI+Vv5Srnny3cYNTczfUf7StM7Wq/QAACAASURBVBlK77ruz/K/oM0Ly/qrpJUT8R7yXevNlN9ZX5bUJxHP7bpucCIW7Z5R/paKhvZ9LV+hvUpS6xDfPcy7fSOv61byVdB5kj6RdHZOvKEbv1Vyxt8h6ZWU5d6upbv/+7Y98k8/fzG8Ll/LP6RruKQu1d7XGIobyM2l1lu13Ax5+YX8Mzr+ne+1Cbn3fNjuL+QfzrxtzvuWL3cnVntfYyh8IC+XWm+tHzPbyZ+E/k/+qrvHJW2RiLeRL8LOlL/N61lJu1Z7P2MofiA3l1pvNY+Z28sfA78Iefec/INcG+L7SHpV/gv0HPkej05qaFtiuk+19DFzXrX3NYbiBnJzqfXWbG6GabrKF9W+lO+K/T5JPXLet7o5n7WwUQAAAAAAAMgobtECAAAAAADIOAo8AAAAAAAAGUeBBwAAAAAAIOMo8AAAAAAAAGRc5gs8ZjbYzOY0PmV9MLNRZvZgtdshSea7LHdmtnW124LaQ25WD7mJNORm9ZCbSENuVg+5iTTkZvWQm8Wr+QKPme1vZk+Y2Qwz+8rM3jCz4Wa2UrXblk+5EoKduTRmNsTMnjGzL8M+86SZ9a12u+oRuUluFsPMljGz35jZ+2Y2z8xeM7Pdq92uekRukpvFMLMfm9nLif3lVTM7qtrtqkfkJrlZDHKzcshNcrNUZnZoeA1roiAm1XiBx8yGSxot6VVJe0vaRNIpktaRdGIVm1Z3zGyZarehTPpJukPSDyRtJ+kdSWPNbINqNqrekJuVU0e5eYH8vvFz+f3lWkn3mtmWVW1VnSE3K6eOcvNz+fzcXtIWkv4iaaSZ7VnVVtUZcrNyyE0Ug9ysnDrKTUmSma0r6VJJz1S7LUtwztXkIGlbSU7SLyLxruHfwZLm5MSOl/QfSQvCv0PyxN+VNE/SZ5LGSmqTiB8taUKIvytpqKRWjc0v6dzQ5uTQL9L+VpJ+LekjSfMlvSHpR4l47nKeCuNHSXpQ/oPnE0lfyn/gd0jMa5LOkPS+pLlh2YMS8R5hmYdKeiJMc3KknW0lXShpSmjnB5J+nrOcrcPfrSWNlDQpLPO90I7ka7e5pMclzZI0R9Jrkr4fYstIulrS1LCujyRd3MT9yCR9Kuln1d6n62UgN8nNUnIzzHtKzri7Jd1a7X26XgZyk9wsJTcj2/CKpIuqvU/Xy0BukpvkZm0O5Ca5WWpuhuW8KOmohter2vvzt22rdgNSXrSrwhuyTCPTLZFwkvaTtFDSyZI2lPSz8Pc+Ib61pG8kHS5pbUnfCwnVJsSHSJom6UD5yu0+8gWCkxubX1In+atHHpW0ShjaRto9NOx0h4V2ni9pkaReIb5N2Jl3C8tZIZFwMyVdL6mnpIGSZkg6K7Hs4fJXruwetuEwSV9J2isnUSYntnONSDtvk/SxpAMkrSvp+5KOjCTcMmE7tgmxg0Lbjkks7w1Jt0raWNL64f3aIcT+LyTZzpLWkrSjpKMT854ryRW5H7ULbRhUzHwM5KbIzbLmpvwvkcfnjLtV0uRq79P1MojcJDebeNyUP2HvH7Z9j2rv0/UyiNwkN8nNmhxEbpKbJeampEsk3ZF4vSjwFPCiPSTptQKmG6wlE+5ZSTfmTDNK0rjw//3DDts5srwPJR2RM+5USRMKnL+gN1i+GvqbnHFPKfyanbsz5yz/I0mtE+Oul/RY+H9H+WrmTjnzXSnpoZxl/18jbdwgTLd7JJ63jTnTXNzQtvD3LElHRaa9Wr7aapH4yZImFrkfXSr/gbFctffpehnITXIzT7zR3JT0N0lvS9pI/helAZK+ljS/2vt0vQzkJrmZJ17QcVNSF/kvOQvlfy0+prF5GAofyE1yM0+c3KyBgdwkN/PECzmnHShfuGq4wqug96NSQy0/g8dKnK+nfNIljZO/n1Ly1c4pkiaZ2V/N7Cgz6yxJZtZd0pqS/mxmcxoG+Z1mvcbmL5SZLSdptUbamWaCc25R4u+pkhoeAraJpPaSxuRsw4mJbWjwciPr2VLSYklPFtAmSZKZnRAeCDc9rHeofHW0wRWSbggPMhtmZhsnYqMk9ZL0rpn90cz2MrNv91Hn3AjnXHL6xtpyivzljfs752YVOh8aRW7GkZtxp8j/0jNB/nLmEfKX+y4udBvQKHIzjtxMNzssZxtJwyRdYWb9C90GNIrcjCM305GbzYvcjCM386+7e1jGUc65GYW2uZJqucDzrqT1zKxtmZbnJMk5N1vSVvKXc30o6SxJE81sNX33epwg/8Y3DJtJ2rSA+cvWzkYszDNPQ9sb/t1HS27DpvLVxqSvSmxjXmZ2sHzldpT8pX69JP1J/r5K31DnzpX/ULhP/pK4183sJyH2inyV9qywHTdJejSZdEW05VT5B9Pt6Zx7qdRtQl7kZhy5GeGcm+6c21f+V5+15S+bnSN/nzXKg9yMIzdTOOcWO+f+45x71Tl3ufwDR89uyrZhCeRmHLmZgtxsduRmHLmZ36aSVpX0uJl9Y2bfSDpS0p7h742auIlNVssFnr/JfxE4OV/QzLpG5ntbUp+ccX3lfzWWJDnnvnHOPeGcO0v+qfQdJe3tnPuvfHVyvfBhusTQ2PwhvED+4U9R4WqSqY20c0H4N3VZeUyQf2DU2nm2YUqRy3pVfh/5foHT95X0Yqh8vhJes9wqrpxz7znnrnbO7SX/kKxjE7HZzrm7nHMnStpLvjes9YtptJn9QtJv5e8BHVfMvCgIuemRm0XmZljOPOfcJ/L3kR8g6e/FLgNR5KZHbpaQmzlayT/DDuVBbnrkJrlZa8hNj9wsPDf/Jf8Q52Rh6375nrR6yT/8uaraVLsBMc65F83sEkmXmtka8r2tfCz/gKZj5J9Wfl6eWS+VNNrMxkt6RP7BT4fL38soM9tbfif4p6Qv5HemzvKJKknnSPqDmc2Qvy9zGfkK6urOuYsKmH+ypD1C9e5zSTOdc7kV0IZ2nm9m70kaL2mQpJ3CuiTpf/L3Nu5mZpMlzXPOzSzgdZttZpdJuszMLLSzk3wXi4udc9c1tozEst41szvlL3E7Rf7J/WtI6uGcuyXPLO9KGmxme8i/P4dI2kX+yesys2UlXSb/68NkSSsrJGmI/0L+gWOvyleND5O/h/LjED9Z/uFjaZfNnS7/0K9B8pferRJCcwt5/dA4cpPcVGm5uZ2k1cMyVpd/iF0r+YfUoQzITXJTpeXmsLC8D+S/OO4p6Qj5h4aiDMhNclPkZk0iN8lNFZmbzrmvJL2ZHBfexzbOuTfzzVNxrgYeBJQ2SPqx/MOgZspf4vWm/Bf47iE+WEt3W3eC/Bu+UDnd1sm/wU/KJ8PcsLyjc+Y/VH7nmie/s4yTdEgh80vqLp/os6WCu61bIP+0731zpjlW/rK8Rcrpti5nunMlvZn42+Q//Buqq9Pl7+Uc4Ap8WFViWe3kv4B9Epb1vr57wvsSy5G/NG5keM1mhP//RqGXnBD/m3yyzZevKl+n8ABk+SfKvxJeu1mSnpa0Y852ukbaOzm0KXcYVe19ud4GcpPczNlO10h7d9F33YF+JulmSatVez+ux4HcJDdzttM10t6L5LuZnSv/ReI5SYdWez+ux4HcJDdzttM10l5ys0IDuUlu5mynK3L/Wer1quZgoVEAAAAAAADIqFp+Bg8AAAAAAAAKQIEHAAAAAAAg4yjwAAAAAAAAZBwFHgAAAAAAgIxrUjfpZsYTmtGiOees2m3IRV4C+sw5173ajchFbgLkJlCjyE2gNhWdm1zBAwCoN1Oq3QAAeZGbQG0iN4HaVHRuUuABAAAAAADIOAo8AAAAAAAAGUeBBwAAAAAAIOMo8AAAAAAAAGQcBR4AAAAAAICMo8ADAAAAAACQcRR4AAAAAAAAMo4CDwAAAAAAQMZR4AEAAAAAAMg4CjwAAAAAAAAZR4EHAAAAAAAg4yjwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyDgKPAAAAAAAABlHgQcAAAAAACDjKPAAAAAAAABkHAUeAAAAAACAjKPAAwAAAAAAkHEUeAAAAAAAADKOAg8AAAAAAEDGUeABAAAAAADIuDbVbgAAAKh9bdu2jcb69+8fjR122GHR2IorrhiNrbbaatHYG2+8EY1NmjQp7/jf/OY30XkAAADqAVfwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxtGLFgBk2Prrrx+NvfPOO9FYq1bx+v7ixYujsb322isaGzNmTDSG2tGjR49obOjQodHYoEGDorGuXbtGY2YWjU2fPj0a69SpUzS23HLLRWOxXrvSet4aPXp0NAaU21lnnVX0PG+++WY09sADDzSlOQCQ1+233x6NHXTQQdHYOeecE4399re/bVKb0Diu4AEAAAAAAMg4CjwAAAAAAAAZR4EHAAAAAAAg4yjwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxdJMOACXYZZddip7n6aefLvu6zjjjjGjMOReNpXWFnjZfWgy1Zf311887/plnnonO07p162gsrUvzf/3rX9HY3XffHY3df//90Vj79u2jsTQffPBB3vFt27YtaXlAuQ0fPjzv+LTP12+++SYamzdvXpPbVE1mFo2ldbe8cOHCktY3cODAaOyWW26JxkaPHl3S+oBKWXPNNaOxY445Jhq74YYb8o5POw6nfV6ddNJJ0dhDDz0UjY0fPz4aQ+G4ggcAAAAAACDjKPAAAAAAAABkHAUeAAAAAACAjKPAAwAAAAAAkHEUeAAAAAAAADKOAg8AAAAAAEDG0U16Cdq0ib9sPXr0iMYGDRoUjXXq1KkpTVrKPffcE41NmDAhGpsxY0ZZ2wFkWVpXjn379s07/pNPPonO07Nnz2isa9eu0dhf/vKXaGyttdaKxtKk5Xpad+50YZkd/fr1yzv+3nvvjc6T1m35mDFjmtqksoltW5qxY8eWvyFAhaSde5b7HLLS0rpJv+yyyyrYEql79+7RGN2koxb06tUrGnvyySejsc6dO0djvXv3zjt++vTphTcsIS2P9tprr2iMc8zy4AoeAAAAAACAjKPAAwAAAAAAkHEUeAAAAAAAADKOAg8AAAAAAEDGUeABAAAAAADIOAo8AAAAAAAAGdeiu0n/3ve+F43ttttu0djee+8djfXp06dJbSqXoUOHRmPvvvtuNHbWWWdFY/fdd1+T2gRUS1p3jcOGDYvG0j4H3nrrrbzjL7nkkug8aV3BpnVdXWpX6GluueWWaOzUU08t+/pQeTfccEO1m9Akq6++ejR2zz33FL28zz77rCnNAcrmpJNOyjt+wIABZV9X2vGvVs5ZARTumGOOicbSukJPM27cuLzjBw0aVNLyUF1cwQMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyDgKPAAAAAAAABlHgQcAAAAAACDj6r6b9OOOOy4aGzx4cDS23XbbRWOffvppNHbjjTdGY8OHD4/G5syZE43FrLTSStHY/vvvH42dc8450dhtt90WjY0ePToaO/LII6MxoBLS9sHTTz89GuvZs2c09tprr0VjAwcOzDs+rSvmK6+8MhrbeeedozHnXDSWJq278xEjRpS0TKCcunbtGo3dfvvt0ViHDh2isYMPPrhJbQKa27XXXlvU+KbYddddo7GxY8eWfX2lmDJlSjT23nvvlbTMt956Kxr7/PPPo7F77rmnpPUBlbL55puXNN+XX34ZjV1zzTV5xzdHN+mbbbZZ2ZeJJXEFDwAAAAAAQMZR4AEAAAAAAMg4CjwAAAAAAAAZR4EHAAAAAAAg4yjwAAAAAAAAZBwFHgAAAAAAgIyri27Sd9lll2jswgsvjMbatWsXjZ188snR2MiRI6OxBQsWRGPlltYd84QJE6KxtO4h07pO3m+//aKx3r17R2Pjx4+PxoBi9OjRIxobPnx4NLbqqquWtL5YV+hSPP/S2njEEUeU1I40kydPjsZuueWWsq8PKFZaV+hXX311NLbjjjtGY++++240dtdddxXWMKAFWGONNSq6vkWLFkVjsXPytGPV+++/3+Q2AVnz05/+NBrbeeedozHnXDQ2ffr0aGzWrFmFNSzBzIqeR5KmTp1a0nwoHFfwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyLi66Cb9gQceiMY6duwYjZ133nnR2DXXXNOkNtWyG264IRo78sgjo7Ftt902GuvevXuT2oSWpW3bttHYWmutFY3de++90VhaV+gLFiyIxi6//PJoLNYVuhT/bDnjjDOi83Tp0iUaK9W0adOisZkzZ5a0zLR8TutmEy3XcsstF42lHaPTukJPy9tjjz22sIYBLUDnzp2jsaFDh5Z9fV988UU0Nnjw4GjsH//4R9nbAtSjH/7wh9FYWlfoH330UTS2//77F92O++67LxrbZJNNorG0Nvbu3bvodqA4XMEDAAAAAACQcRR4AAAAAAAAMo4CDwAAAAAAQMZR4AEAAAAAAMg4CjwAAAAAAAAZR4EHAAAAAAAg4+qim/S07iEXL14cjc2ePbs5mlPzFi5cGI3Nnz+/gi1BS5XWFfrEiRPLvr60rtB//etfl7TMnXfeOe/44447rqTlleqCCy4oab6rrroqGttpp52iscMOOywaa473DrWjS5cu0dj9998fjfXp0ycaS+tKdebMmdHYVlttFY2tsMIKecdPmjQpOk+nTp2isbT9esaMGdEYUCl9+/aNxjbaaKOyr699+/bR2EEHHVR07KmnnorOc9NNN0Vjaef4QK3bddddo7F+/fqVtMwnn3wyGivlHK05jnFpy+zVq1c0Nnfu3Ghs+vTp0dgXX3xRWMPqCFfwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyLi66Cb9hBNOiMbeeeedaOyNN95ojubUvPXWWy8a6927dzQ2derUaGzcuHFNahNalquvvjoaM7Oyr+/ZZ58t+zJjSm1/q1bxevubb74ZjW244YbR2JlnnhmNpXXBmdb1bI8ePaIxukmvb2uuuWY0ltZNc6m6d+8ejV155ZVFLy8tN9O6a//qq6+isWeeeSYaS+v6edq0adHY6NGjo7H58+dHY2i5Hn744WjsggsuiMbOO++8ktbXoUOHaGzQoEFFL++II46IxrbbbrtobOHChdFY2mdEWpfKX3/9dTT2zTffRGNAsYYMGRKNtW7duqRljho1Khpr27ZtNLbjjjvmHX/SSSeV1I40e+yxR0mxNLNnz47Gzj///Gjs97//fUnrq3VcwQMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyDgKPAAAAAAAABlXF71o/fnPf652EzIlrYeDtJ4R/vSnP0Vjc+bMaVKbUH922WWXaCytx5203mxKtfXWW0djaT3rTJkyJRrbZJNN8o4vtf1pPVf17NkzGiu1B4C09f3vf/+LxtJ6H0F9S+vNbdddd43GVl999WhsxRVXjMZ22GGHaGzBggXR2Oabb553fFrPVV26dCmpHaX2BpL2uXPzzTdHY8OGDYvGLrroomgMLdett94ajaX1CnX44YdHY7HjX3NI62koTam9/9x+++3R2PDhw6OxCRMmlLQ+1LcVVlghGhs4cGDZ15d2jOjWrVs0tsUWW5S9LeV22223RWMffPBBNJbWM2W94goeAAAAAACAjKPAAwAAAAAAkHEUeAAAAAAAADKOAg8AAAAAAEDGUeABAAAAAADIOAo8AAAAAAAAGWdN6ZLYzMrfnzHK4swzz4zG0rp5vOqqq6Kx008/PRpbtGhRYQ2rM865eF+3VVIreZnWRfADDzxQwZaU7sMPP4zGYt1NdujQoaR1pXWb3Bxdx8+cOTMau+uuu6Kx448/vuxtaQbjnXNbV7sRuWolN+tdu3bt8o6fP39+dJ42bdpEY3369InGDjzwwGisd+/e0djyyy8fjW244YbRWNqx9uSTT47GrrvuumiswsjNOvejH/0oGttuu+3yjt9xxx2j8+y8885NblO5pB2n77zzzmjspZdeisYuv/zyJrWpjMjNCrv22mujsSFDhkRjlT5fjGnVKn6dyOLFi6OxVVZZJRqbPn16k9pUp4rOTa7gAQAAAAAAyDgKPAAAAAAAABlHgQcAAAAAACDjKPAAAAAAAABkHAUeAAAAAACAjKPAAwAAAAAAkHF0k55hAwYMiMbGjBkTjU2cODEaGzhwYDT2ySefFNawFoRu0uPqoZv0SnZFWeluL0eMGBGNnXrqqWVfX4XR3StqWseOHaOxjz76KBrr0qVLNDZlypRobN111y2sYc2P3MRS2rZtG421a9cuGjvjjDOisV69ekVje+65Z2ENK5N58+ZFY7/85S/zjk87RjcTcrPCunbtGo1deOGF0diuu+4ajaXl0gsvvBCNde/ePe/4XXbZJTpP2nnrU089FY3tvvvu0djChQujsRaMbtIBAAAAAABaGgo8AAAAAAAAGUeBBwAAAAAAIOMo8AAAAAAAAGQcBR4AAAAAAICMo8ADAAAAAACQcXSTngFnnnlm3vE///nPo/O89tpr0diQIUOisY8//rjwhoFu0lNssskm0Vha16ZHHHFEczSnJK1axWvgixcvrol1zZgxIxq75ZZborE66Ao9Dd29IrO22mqraOzee++NxlZZZZVoLNbVbVq3uc2E3ERFLLvsstHY8ssvH42tscYa0didd94Zja255pqFNaxArVu3LuvyCkBuZkSnTp2isbRzyVmzZkVj/fr1yzv+sccei86T1k366NGjo7FDDjkkGkNedJMOAAAAAADQ0lDgAQAAAAAAyDgKPAAAAAAAABlHgQcAAAAAACDjKPAAAAAAAABkHAUeAAAAAACAjGtT7Qa0JCuttFI0dskll0Rjhx12WN7xY8eOjc6z7777RmOLFi2KxoBymTBhQjQ2ZMiQaOztt9+Oxg4//PBobPr06dFYWpft3bt3j8bSuid3rrw9d06aNCkau+uuu6KxESNGRGMffvhhk9oEoPLSzhVWWGGFaGzu3LnR2Pvvv9+kNgFZk5YPabF11103Gps2bVo0Vu5u0oGYOXPmlDRf27Zto7G082tkD1fwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyDi6Sa+gtK6O+/TpE429/vrrecffcMMN0Xk6d+4cjXXt2jUamzx5cjRWS2JdxaZ1U/naa681V3NQpIULF0Zjv/vd76Kx22+/PRpL6ya9Z8+e0djxxx8fjR1zzDHRWLn985//jMaGDRsWjaW9lkAtWH311aOxTz75pIItqR3du3ePxk444YSSlnnmmWdGY2mfj8i+Y489NhobMmRINPbWW29FYz/5yU+a1KZK6Nu3bzTWvn37aOzEE0+Mxvr37x+NpZ1bl+r9998v+zKBfNK+/x199NEVbAmaG1fwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyDhzzpU+s1npM2dYWjeJl19+eTQ2ePDgaKx169ZNadJSPvvss2isVat4XW/OnDnRWJs2bZrUpnymTp0ajY0ZMyYai3Vxueyyy0bnaY7uLZ1zVvaFNlFLzcvmsGjRomisKZ+d+ZjFd6Wzzz47GkvrVr4FG++c27rajchVz7m5++67R2PHH398NLbffvs1R3MqJu24eMABB0RjaTm92WabRWMjRoyIxk455ZRorIaQmyXq1q1bNPbSSy9FY2uvvXY0NmvWrGhs6NCh0diECROisTS77bZb3vG77LJLdJ60Y21aN+lt27YtvGHN7Ouvv47GYvk+ZcqU5mpODLlZ51ZaaaVobNq0aUUvL+375B133BGNHXLIIUWvq4UrOje5ggcAAAAAACDjKPAAAAAAAABkHAUeAAAAAACAjKPAAwAAAAAAkHEUeAAAAAAAADKOAg8AAAAAAEDGlb/f6xZg++23j8aOOeaYaGz48OHR2GuvvdakNhUjrVvMtK5u07rD69ixYzT29NNPR2PDhg2LxpZbbrlobOHChXnHr7jiitF5gGI988wz0dhOO+1U1nWl5deFF14YjaXlSVp+AcXq2rVrNDZ69OhobOrUqdHYySefHI2NGTOmsIaVwY9//ONoLO24su+++0Zj66yzTjRmZtHYH//4x2jstNNOi8ZQ39K6BF911VVLWmba8WPkyJElLbMUace/xYsXV6wdkjRv3rxobPLkydFY2rn1BRdcEI1VoTt0YCnOFd8bfaVzE4XjCh4AAAAAAICMo8ADAAAAAACQcRR4AAAAAAAAMo4CDwAAAAAAQMZR4AEAAAAAAMg4CjwAAAAAAAAZRzfpJXjyySejsY022iga+/LLL6Oxzz//vEltKpdHH320out78MEHo7GZM2dGYxdffHHe8Z07d25ym4AGV1xxRTS21VZb5R3foUOHktaV1t1kWveVPXv2LGl9QLHSuvZeZpllorENNtggGrvqqqua1KZipLW/lC5iGzN79uxo7LbbbovGhg4dGo198803TWoTsuu+++6Lxl588cVoLO0Y0a1btya1qdq+/vrraCztnPuaa66Jxl599dVo7OGHHy6sYUAL9/LLL1e7CS0aV/AAAAAAAABkHAUeAAAAAACAjKPAAwAAAAAAkHEUeAAAAAAAADKOAg8AAAAAAEDGUeABAAAAAADIOGtK16BmVv5+RYEMcc7F+92tEvKyMo4//vi8488444zoPGuvvXY0ltaF86RJk6KxwYMHR2Pjxo2LxurceOfc1tVuRK56zs0DDzwwGttiiy2isd69e0dju+++WKtrfQAAAotJREFUe5PalOu5556LxtJyLM3YsWOjsUceeSQamz59eknrqwPkZoWtueaa0Vha1+vrrbdeSeu7//77o7G07txLMXHixGjs8ccfL+u6WgBys861a9cuGvvDH/6Qd/whhxwSnadjx47R2A9+8INo7Omnn47GkFfRuckVPAAAAAAAABlHgQcAAAAAACDjKPAAAAAAAABkHAUeAAAAAACAjKPAAwAAAAAAkHH0ogU0Ab1oIVfXrl2jsXvuuaek2C233BKNzZw5s7CGtSz0BgLUJnITqE3kJpbyq1/9Khrr379/NDZkyJBo7D//+U+T2tQC0YsWAAAAAABAS0OBBwAAAAAAIOMo8AAAAAAAAGQcBR4AAAAAAICMo8ADAAAAAACQcRR4AAAAAAAAMo5u0oEmoJt0oCbR3StQm8hNoDaRm0Btopt0AAAAAACAloYCDwAAAAAAQMZR4AEAAAAAAMg4CjwAAAAAAAAZR4EHAAAAAAAg4yjwAAAAAAAAZBwFHgAAAAAAgIyjwAMAAAAAAJBxFHgAAAAAAAAyjgIPAAAAAABAxlHgAQAAAAAAyDgKPAAAAAAAABlHgQcAAAAAACDjKPAAAAAAAABkHAUeAAAAAACAjKPAAwAAAAAAkHEUeAAAAAAAADKOAg8AAAAAAEDGUeABAAAAAADIOAo8AAAAAAAAGdemifN/JmlKORoCZNDa1W5ABHmJlo7cBGoTuQnUJnITqE1F56Y555qjIQAAAAAAAKgQbtECAAAAAADIOAo8AAAAAAAAGUeBBwAAAAAAIOMo8AAAAAAAAGQcBR4AAAAAAICMo8ADAAAAAACQcRR4AAAAAAAAMo4CDwAAAAAAQMZR4AEAAAAAAMi4/wcyA/PTbDniXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    ax = plt.subplot(1, n, i+1)\n",
    "    plt.imshow(imgs_low[i].reshape(28, 28))\n",
    "    plt.title('Model prediction: {} (p={:.2f}) \\n Label: {} \\n Trust score: {:.3f}' \\\n",
    "              '\\n Closest other class: {}'.format(pred_low[i].argmax(),pred_low[i].max(),\n",
    "                                                  label_low[i], score_low[i], closest_low[i]), fontsize=14)\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see several examples of an over-confident model predicting the wrong class, the low trust score, however, reveals that this is happening and the predictions should not be trusted despite the high model confidence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following section we will see that on average trust scores outperform the model confidence for identifying correctly classified samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of Trust Scores with model prediction probabilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s compare the prediction probabilities from the classifier with the trust scores for each prediction by checking whether trust scores are better than the model’s prediction probabilities at identifying correctly classified examples.\n",
    "\n",
    "First we need to set up a couple of helper functions.\n",
    "\n",
    "* Define a function that handles model training and predictions:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_sc(X_train, y_train, X_test):\n",
    "    clf = sc_model()\n",
    "    clf.fit(X_train, y_train, batch_size=128, epochs=5, verbose=0)\n",
    "    y_pred_proba = clf.predict(X_test)\n",
    "    y_pred = np.argmax(y_pred_proba, axis=1)\n",
    "    probas = y_pred_proba[range(len(y_pred)), y_pred]  # probabilities of predicted class\n",
    "    return y_pred, probas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Define the function that generates the precision plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_precision_curve(plot_title, \n",
    "                         percentiles, \n",
    "                         labels, \n",
    "                         final_tp, \n",
    "                         final_stderr, \n",
    "                         final_misclassification,\n",
    "                         colors = ['blue', 'darkorange', 'brown', 'red', 'purple']):\n",
    "    \n",
    "    plt.title(plot_title, fontsize=18)\n",
    "    colors = colors + list(cm.rainbow(np.linspace(0, 1, len(final_tp))))\n",
    "    plt.xlabel(\"Percentile\", fontsize=14)\n",
    "    plt.ylabel(\"Precision\", fontsize=14)\n",
    "    \n",
    "    for i, label in enumerate(labels):\n",
    "        ls = \"--\" if (\"Model\" in label) else \"-\"\n",
    "        plt.plot(percentiles, final_tp[i], ls, c=colors[i], label=label)\n",
    "        plt.fill_between(percentiles, \n",
    "                         final_tp[i] - final_stderr[i],\n",
    "                         final_tp[i] + final_stderr[i],\n",
    "                         color=colors[i],\n",
    "                         alpha=.1)\n",
    "    \n",
    "    if 0. in percentiles:\n",
    "        plt.legend(loc=\"lower right\", fontsize=14)\n",
    "    else:\n",
    "        plt.legend(loc=\"upper left\", fontsize=14)\n",
    "    model_acc = 100 * (1 - final_misclassification)\n",
    "    plt.axvline(x=model_acc, linestyle=\"dotted\", color=\"black\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The function below trains the model on a number of folds, makes predictions, calculates the trust scores, and generates the precision curves to compare the trust scores with the model prediction probabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_precision_plt(X, y, nfolds, percentiles, run_model, test_size=.2, \n",
    "                      plt_title=\"\", plt_names=[], predict_correct=True, classes=10):\n",
    "    \n",
    "    def stderr(L):\n",
    "        return np.std(L) / np.sqrt(len(L))\n",
    "    \n",
    "    all_tp = [[[] for p in percentiles] for _ in plt_names]\n",
    "    misclassifications = []\n",
    "    mult = 1 if predict_correct else -1\n",
    "    \n",
    "    folds = StratifiedShuffleSplit(n_splits=nfolds, test_size=test_size, random_state=0)\n",
    "    for train_idx, test_idx in folds.split(X, y):\n",
    "        # create train and test folds, train model and make predictions\n",
    "        X_train, y_train = X[train_idx, :], y[train_idx, :]\n",
    "        X_test, y_test = X[test_idx, :], y[test_idx, :]\n",
    "        y_pred, probas = run_sc(X_train, y_train, X_test)\n",
    "        # target points are the correctly classified points\n",
    "        y_test_class = np.argmax(y_test, axis=1)\n",
    "        target_points = (np.where(y_pred == y_test_class)[0] if predict_correct else \n",
    "                         np.where(y_pred != y_test_class)[0])\n",
    "        final_curves = [probas]\n",
    "        # calculate trust scores\n",
    "        ts = TrustScore()\n",
    "        ts.fit(enc.predict(X_train), y_train, classes=classes)\n",
    "        scores, _ = ts.score(enc.predict(X_test), y_pred, k=5)\n",
    "        final_curves.append(scores)  # contains prediction probabilities and trust scores\n",
    "        # check where prediction probabilities and trust scores are above a certain percentage level\n",
    "        for p, perc in enumerate(percentiles):\n",
    "            high_proba = [np.where(mult * curve >= np.percentile(mult * curve, perc))[0] for curve in final_curves]\n",
    "            if 0 in map(len, high_proba):\n",
    "                continue\n",
    "            # calculate fraction of values above percentage level that are correctly (or incorrectly) classified\n",
    "            tp = [len(np.intersect1d(hp, target_points)) / (1. * len(hp)) for hp in high_proba]\n",
    "            for i in range(len(plt_names)):\n",
    "                all_tp[i][p].append(tp[i])  # for each percentile, store fraction of values above cutoff value\n",
    "        misclassifications.append(len(target_points) / (1. * len(X_test)))\n",
    "    \n",
    "    # average over folds for each percentile\n",
    "    final_tp = [[] for _ in plt_names]\n",
    "    final_stderr = [[] for _ in plt_names]\n",
    "    for p, perc in enumerate(percentiles):\n",
    "        for i in range(len(plt_names)):\n",
    "            final_tp[i].append(np.mean(all_tp[i][p]))\n",
    "            final_stderr[i].append(stderr(all_tp[i][p]))\n",
    "\n",
    "    for i in range(len(all_tp)):\n",
    "        final_tp[i] = np.array(final_tp[i])\n",
    "        final_stderr[i] = np.array(final_stderr[i])\n",
    "\n",
    "    final_misclassification = np.mean(misclassifications)\n",
    "    \n",
    "    # create plot\n",
    "    plot_precision_curve(plt_title, percentiles, plt_names, final_tp, final_stderr, final_misclassification)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detect correctly classified examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The x-axis on the plot below shows the percentiles for the model prediction probabilities of the predicted class for each instance and for the trust scores. The y-axis represents the precision for each percentile. For each percentile level, we take the test examples whose trust score is above that percentile level and plot the percentage of those points that were correctly classified by the classifier. We do the same with the classifier’s own model confidence (i.e. softmax probabilities). For example, at percentile level 80, we take the top 20% scoring test examples based on the trust score and plot the percentage of those points that were correctly classified. We also plot the top 20% scoring test examples based on model probabilities and plot the percentage of those that were correctly classified. The vertical dotted line is the error of the classifier. The plots are an average over 2 folds of the dataset with 20% of the data kept for the test set.\n",
    "\n",
    "The *Trust Score* and *Model Confidence* curves then show that the model precision is typically higher when using the trust scores to rank the predictions compared to the model prediction probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = x_train\n",
    "y = y_train\n",
    "percentiles = [0 + 0.5 * i for i in range(200)]\n",
    "nfolds = 2\n",
    "plt_names = ['Model Confidence', 'Trust Score']\n",
    "plt_title = 'MNIST -- Softmax Classifier -- Predict Correct'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Reshaping data from (48000, 4, 4, 4) to (48000, 64) so k-d trees can be built.\n",
      "Reshaping data from (12000, 4, 4, 4) to (12000, 64) so k-d trees can be queried.\n",
      "Reshaping data from (48000, 4, 4, 4) to (48000, 64) so k-d trees can be built.\n",
      "Reshaping data from (12000, 4, 4, 4) to (12000, 64) so k-d trees can be queried.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEeCAYAAAA9/ZrRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydd3hUZdbAfyedQKgJvUQUCwo2sK+gIhbWAu6uFcWGjXXdVdfu2kFhd2U/0bUsKDZ0EQuKBV2xIEoRUEQRRJSEUEJPgbT3++PcIcNkkswkk5lMOL/nmWfuvPe99565M/eee8p7XnHOYRiGYRjxREKsBTAMwzCMcDHlZRiGYcQdprwMwzCMuMOUl2EYhhF3mPIyDMMw4g5TXoZhGEbcYcprD0VEficii0WkWESciAyMtUyNARFZJSKzYi0HVC+LiFwrIj+IyE7vt8sWkXt8y1EXtAkiIiMCr4tgbUbsaBDlJSIDvR/Zichj1fRpLyIlXp9ZAetmee0rRSQlyLa+C7VfkGPeFNC3lYjcKSKLRGSLiBSIyM8i8oaIXOH1GeEnb22vWUSJUGSv4373BV4GtgKjgOHA9945vEdEWkfoKzQKRKSNiNwtIvO881giIjki8pqIDBMRibWMoSIiJwATgB+Aq9HfbkNMhYoAAfcM36tARBaIyJ9EJDHWMtaH+lxbIrKviDzuPbAUeg+cP4rIUyLSvyHkjTZ+D2CHhLpNUkMKBOwALhCRG51zOwPWDQcEKKth+72Aa4DxdTm4iLQE5gE9ganARKDE+3wc8CfgGeBTTx5/7gD2D9K+ri6yhEsYsteFgehvf4Nz7mu/Y14D/A14FthSx303KkTkCOBNoD3wFvAisA3oApwOvAZcBzweKxlrYD8gsIrAyd77Zc65Tb5GEXkAGAMEXmfxxsvADPTe0BkYATwKHAiMjJ1YADwPTEGvw3AZSB2uLRG5HHgCvZe+DCxC75n7AucAV4rIgc65pXWQqTGRjZ6fVeh3rJWGVl6vA+cDZwGvBqy7FP2TnlTNtsXAz8CdIjLRObe9Dse/EuiF3qSrKEAR6QjgnFsJrAxYdwWwv3PuhTocNxKEJHsd8W27qcZecY53jqYDacAA59znAV3uF5FTgDZRFy4Egjzwgffb+Ssu73MZNT8Iho1n7aQ654oiud9a+Nr/mhORJ4DvgStE5C7nXNCHRxFJBhKdczsaSjDnXDlQ3lD7D0REBgFPAUuBU5xzawLW3wb8MYLHE6C5c66gmvUZdbwPNwzOuYi/0KcMB9wELAbeDVh/hLf+bKAAmBWwfpbXfqbX7/6A9fd47f2CHdOv7d9e20F1+A6z9PRE/vyEePywZQeOB2ai7sBi4Gvg8oA+LshrFfpEGGzdPQHnvDf6JJwHFAEfAft5fYZ5xyz29jkyiIznohbQr6iVkA+8AfQN6DcEqAD+E9DeAliGWsAdazkff/dkvjqMc7gqyP9xMPAK+oBTjD45f4AqxMDtDwT+C+R6328t8DEwxK9Pmnc+l3nncAvwLTC2OlnQJ9Ngv49vve/3yQ7YRyvgYWCFJ88G9Am+Z0C/Ed72g4C7gJ+AUmBElP7vAwm4fv3WTfXWHR3wXQ8E/gHkoEploN82g7zfaAtqtXxT3f8AfVD8wTs/K4Ab0IdrF7DPEYFtXnsK8FfUYihCr7/5wChv/bPV/Hb31HJOFqDXQO8wzmNzYLT3+/n+f5OBHtWc7xGo52Gp1993vc/y/n8+z88m/O6HqGV8jSdjEXq//hg4oRq5zvH2ucXrvwz4l3fufOc16H+7uldDW16g7q5/iEgX51yu13YZsB54u6YNnXNvicjnwJ9FZIJzbm2Yx/7Je79URG5x+nQaL4Qlu4icgVq6a9Gb9nbgPOAZEenpnLvD6zocVTJDgT+jyqMAVQYtA9pBL3p/nvP6PwRkATcC74vIXcAjqItjInA58KSILHW7WzyjgI3oE+VaYG/UHTRbRA5zzi0HcM69IyKPor/9TOfcFG/7x1GL9PQQ/g/noC6e52rpVxsjgLboTSAHdTleAXwkIic45z4DEJF2wP+8bf4N/AJkAv2AI4F3vHUT0GtgMnrzTfK+04k1yLAB/e1GAr+h0p1drRtbRFoBXwDd0d/kO6ATcC3wlYj0c879ErDZOCAZeBp1ry6rQaYGx7MG9vE+5gesfhF9mPA9pOR524xEz/+XwINAIepufUJE9nbO3ey3/xuAf6IP2bcD6ehD9/oQ5UsB3keVwQfAC6iy7INeZ48BTxLateW/372Aw4DPXIguQc/6fB84FlU4f0f/V9cAg73fOydgsxuAdujvvRZY7beuBfAJMBsNo7T3W/c86lWbCkwCUoELgZkiMsw595afXA+i53Ypeq7z0Ov+HOBuNGzzkNfnKeAzb9OaQzQN/RTlnZidwO3eumao9h3nfa7W8vKWj/H29W+/9fcQmuXVBn3Cd96JmArcgsaMEmr5DrOIreUVsuxAInqj3AJ09mtPQf945UCvIOcvO2A/QdsD1k0HxK/9eq99G9DNrz0Lz08fsJ/mQfZ9gPcfeTygPQV9gt2KPgEO9441LoTzl+H1/SbM874qyP8xmMwd0JvQDL82n6fgD7UcY5P/dmHK8myw/2Ww3w6NFRcDBwf07eH9Xs/6tY3wtl8GpMfg/+67fu9GFX4W0Be9qTpgTpDvOgtICthPJ+9/91KQY4z3roWe3ufWqGJb6v+dga7ofalWywu1uBzwUJDjJQSRucq1Vc35OMPr/68wzuGV3jaPBLQP8dqfD3K+NwHtg+xrlrf+gSDrhnrrRga0J6HX68949wgqvWz/A9IC+otfP588I0L9vg2eKu+c24i6iUZ4TcNQV8bEELf/AnUrXe5lyYVz7M3A4ajbZCuq6cegmv0nERkczv6iSZiyH473dO38/OLOuRLUGkpA446R4F/O+7d5+J6S3nLO7Xpqc85tQG+EvQK+VyHoE7WItBSRTNSqWIZaJ/59S1A3o6BW5ePoxXFbCHK29N63hfi9qsUnsyd3C8/CKge+CpB5q/d+mpdwUx1bgQNF5KD6ylYdnsVyIfpUmysimb4XesP+EnWHBvKEi26MK5B70f/DetQaugy9f5wdpO+jrqpH4neoFfAf/+/sfe/p6LUwyOs7GLW0Jvh/Z6fWyYshynshsBm4L3CFc64ixH0Eoy7/36Gom3F0gBzvoC7Ns0Qk8J4/2TlXk5U5LkjbRahn542A89saPcfZVF73F3rvt7mAeKTzqP1rBScabkNQs/IdETkO/TPOdeFlx9yGPomMRm/iIePdRG8FbvVuOkcDf0B/gNdF5GDn3Ipw9lkdnpumWUCzL405K6C92Dm3lRoIQ/a9vE2+C7IbX1vP0L5FrawM+LzZe/85SN/N6FP+LkTkUOB+9EmreUD/Kvtwzv0kIn9Bn8CLgfOdc6UhyOm76DNC6FsjIrI36n46Bb1AdxPRT9ZPRGQy+qB2oYjMAz4EXgn4v9+Aul2+FZGVaKxgOjC9njc8f7JQr8dgqk+lD3asH0M9gIhkoVa/j3Ln3AYRaYY+oPqz1TlXHMJun0Jjhg5Vsj+6gOSUWmQ9wHv/sIZjdPDefdfED0H6hHp/6gUsCrwxR4C6/H/3AtZ4D76BfAccglq1/sqqpt97g3MuWGbkAZ5cNbn1Onj77oX+lotr6FsnoqW83kcD2H8DTkB9sCHjnPtBRCahGUdH1rpB9fvZiMbZ3haR1aiP9TzggbruM4DxwCUBbT7FEnhjfo5Ka7RWoiB7qFSXbVVd+64xVCLSHbUEtqEKbBl6g3JoEkiLavZxhvfeDE0fr/Vhwzm3XUR+AfYXkWYh3jirCi/SwpO5uSfjt+hTZwX6ULVbnMo5d4mIjAVOQ2NTNwJ3iMgNzrnHvD5vig4mPh0YgFoClwOficggz+KsL77z/iFqvYdKOFbXPHZ/OPkFfeo+F31g9edS1OVZG8udczUpHn+Cyer73hfjxcCCEPgA1hhZ4r0f2sDHqen3rm6doA9EF9Sw7RK/ZV8CRkSJivJyzpV7T6S3oU/PL9dhN39DT9bDqD+2vnzpvXeJwL58PIIGbP3xJRWcHNC+hroTKLvvYjwwSN/eAX1qIuJ/sACGogrqTOfcx/4rPMuySmq4iPwRjSWNQV3Oz4pIX+dcdTcmf6ahAfLh6BN9XTgJHW90mXNutxuyN7aqCs65JejFO9YblPoVMMZLOnJen03of+UFz8U3Bo2fnIVaHvVlAxoDbRmGMgiXC9nd0+B7QHifqv/3YF6BhmC5954fwvf2XRP7o1mz/vQmNH5EH5BSXfChDT7Curaccz+LyELgWBHZ3zkXzDoMZCVwqoi0DmIx9UYfGgOTXurCcnSc2ZeumrR6P35EH+QOBubW0C/se080y0P9G/VnX+2cCzsO4cVyxqNPqqeHso2IHC3Vj2j3+dAjNrjPObfUOfdhwGuH9wpsr/G4Ycr+NZrccan/+C8v++hm9I/xZghfwfdHbBtC37rgs852q2ghIldSOfbMv/1gYCzqVrsDtTRbAs8H8d0H4xH0Jv6IiBwdrIOIDBaR8+og82ACYnQi0jZQLu8m8jMaW0kTkcTA39VTaAu9jxE595778UXgCBH5XbA+ItI+WHsYx5gd8J+e7bXnBfm/h/KwEQleRR+C7vXcl7shWrUm1fs4E1W414lIul+frtRsVfjzIppcdWeQY/n/Z+pybd3ivU+RIOM6vf/SDSLiU7RvoPf0WwP6nYZacG9FyC092TvO6GArRaSD38eXvPeHJHi1JN85Cvv8RMttiHPuVzTjpj48jKYKh1oS5UL0hv4OqvU3onGA01H35VJCTByJASHL7lm2o9Ckhnki8hTq2joXOArNhFpe9RBV8Fl0D4vIi2jW1hLPkogE76KuiOdFy4ZtRtN6T0eHBuz6P4pIc7SawTbgIu+iWygit6DptrdQzcXjwzm3VkR+iyruz0XkDSrdlp2BU9HszZrc2J/jDT/wXH05aOxgOOpC7OPX92I0tf911LVZij5snQK86pwr9hRXnoi8hSqs9VRWktmMxr4ixR3o+X1VRF5Ff98S1NV3OjpGZ0QEjxdznHM5opVinkFLnj2PujOz0N/qbNQKWeWc2yw6xGMc8IXnHUpHy24tJzSX3XjUrX2naKmmD9Dr5kDUxe1LDgn72nLOzRRN+38CWCYi/hU29kHj/3sDvsSfZ9GwxS3ef/VTr9+1aHzq9hC+T60456Z6YZxRInIYGs7IR7M0j/aO2dPrO1dEHkav169F5BX0etoLTa45AvUQLEXvWdeKiG/s43rn3P+ojlDTEsN5UcOAwyB9a0yVD9L/z1T6UGtLlT8IjQnNRt10Jd4JWogq0pY1yDWL2KbKhy07eqOcid6cd3h9Lw/S7x6qT4n/K+p+KCX4IOXsgP7Z/v2CnMNVAW3Howphu/cHfcf7rrv1RRVzBTqey397QS+WUuDIEM9lW9Tt7Eu7L0GV0FTUhenfd1WQ/2Nf4D1UuWz3ZP0NAWnrqFJ7DlVchd7vsBiNe6V6fVJQpet7INnpHXMifsMZapBlt2PW9puiN+O7UEVb7Mn/PZoAc6RfvxEEGYAbxf/7QEK/Z1T7//Xrcyz6MLfe+73XoBb8jVRN2b4Kjb/WdZByGvqg8B163W1B44HXhnJthfB990MV2I/ow98OT94ngUMD+voGKa/0vvd6NDmoRzXne0Q1x5xFwLUbpM9wNNvYd79Zhbrqzw3S93z0XrbduzZ+QGPIKX59Tke9SDsIYZCyL8feMAzDMOIGmxLFMAzDiDtMeRmGYRhxhykvwzAMI+4w5WUYhmHEHVFLlY82mZmZLjs7O9ZiGIZhxBULFizId84FlrNrdDRZ5ZWdnc38+fNjLUZEWb1a695269YtxpIYhtFU8cqqNXqarPJqigwfrlM4zZo1K7aCGIZhxBhTXnHEnXdWqUBjGIaxR2LKK44YNGhQ7Z0MwzD2ACzbMI5YuXIlK1fGw2wOhmEYDYtZXnHEZZddBljMyzAMI2rKS0QmAr9FKwVXmf7cK40/Hi3OWIQWjPzaW3cJlVMOPOCcey46Ujcu7r333liLYBiG0SiIpuX1LPAYOhdMME5Dp4zuhc6T9ARwpIj4KoL3QysNLxCRt1zwqa6bNAMGDIi1CIZhGI2CaM7n9ak3x0x1nAVMdlrm/ksRaS0indDS/TOdzjyLiMxE52Gqy2zMcc2yZcsA2G+//WIsiWEYAGxaBj9MAVdee98Is3o1/PRT1fajj4bUzK7Qd2TUZYomjSnm1QVY7fc5x2urrr0K3sRtIwG6d+/eMFLGkKuuugqwmJdRDUX5ULxBl4vzoWgdFK7V951bQRL0RQKIQPFG2LG5sj2UV0U5FORCyTZwFbr/XftO1P3utk1ikP2E0Ib3npDot+y1JyTu3se/X0ICSBIkBLwkCRKT/dYlwrZfYMO3gPPakv3ek6u0FZcm8+2SJAp3JNP7wCQ6tC2ifPFEEtlJRYU3IbDvzW8O5WCzTlW3XqosVE9XB12qzE0MshBK2hxKiimv+ME59xTwFEC/fv2a3ERlDz30UKxFMMKltBh2bILynXqT3bkVdm7Wd+egolQ/g97A/W+2CYlQsl1fvnU4VUzF+VC0QZVVcT4UrFElFQxJgOTmejzngApVPKmtIbUVIPrZVXjrXPWfRaBZFiS30H1n9oG01n59nN92FWqR+G+/2zr/dv9+Zbv39/UlyP6ogAr/fXp9Ksp1PxUBy66s8rykZEC7AyExRX+H0hJ9ryinvKyMRCmjvLSM9esqSJBSXHk5PRNKSU4pI311KS63nC2tB/PgF/9ku+tKWZlQUgplpXDFZSV06ez44stEpryaTGmZUFYKpaVQVgYP3reDrl0c/52WzIQnUigqFnbsqNRYn88qoH2W418TUnjs8VSSkhztsxxpzRwJAq//t4ju3XZXgv4kpTSpW3tQGtM3zAX86x519dpyUdehf/usqEnViDjmmGNiLULTxadIdt0ond4Eizd6Fsw6T+lsg5KtsO1X2P4r7Nii2yQk6lN62Q4o9ZRN8UYoWtsw8iZnQFobaNYO0tpC617QZl9I90rSpbbV5Ywu0KwDJDer+rgvojLv9pgf8My3m9ngW/brLwmesq3FVKhifgR5tgx6rFrW7Wqvvs05p+I5x/dLYerrCWzeVMp99zhcUnOGX5LM7NmQkQEtWsCmTY62beGN1ytwFY6RVyVQUABt2jj+dH05xx9XTimOEudIT0rikUtSEal8NnDehNnOQc/94cJLdl/nHJSXt6C8HP54A1wzShXazp2wfTsUFEDbzBa4RBh4IrRuCxs3CuvWCatXw5Ah0CU7g5Ztaj7lTZ3GpLzeAkaJyBQ0YWOrcy5PRN4HHhIR3081GLgtVkLGkiVLlgBw0EFVkjUNH2U7Yd18yJ2tFokrh4oy74ndWy4tgtICKClQRbNzK2zPgbKi0I+T1haad4bUliAput+yHZCUCs3a6h2q7QGQ0V37JqaokktpoUonra3nKktUBSSJnnylfpZCqVpMqW0AB+Xe7PHNMrV9l5WWWLvyiDUiu/QlQGGhkJQEqan6ee1aWLxYPxcXw/r1cOKJ0K2bxnXefRdatqx8lZVpbCcjA2bPhtdf121++glWrdLjvP8+dO0KDz8M//iHKqaNG/V4PXvCRZdCWhoceqi+FxVBYSF07y4ceyzsLEnEOXj6GWjtGZetWiWS1MB3zQrPmCwvhy5d4KST9JwUF+v3TkjQc7CnE81U+ZdRCypTRHLQDMJkAOfcv4EZaJr8CjRV/lJv3SYRuR+Y5+3qPl/yxp7GqFGjAIt5AXon2bkVtq6ELT9B/reQ8ynkfQXlO7RPYmplLGXXTT5RFUxSc1Ukyc3VDdblOEjL9OIqXtwG1K3WLAuad1SllNJa3W3J6Wq1BFMau+I2VL43IYqL1Tpo3RqSk/Um+7//wdatsGGDvtavh9NPVwX0449w8smqOBISIClJ9/HII3DFFbBkCRx/fNXjTJoEAwfCe+/BH/9Ydf1rr0HfvjBzJowfD23bqrI67DA9xtat+vPsuy8MH64yZ2fDb3+rSrFVK0hPh7vvrrSIKioqDTdfW2pqdJ8NEhIqz5NPubdure9lZXq+ExOjJ09jRVywaGIToF+/fq6pVZWfN0/1d//+/WMsSYQp2Q7bczURYPtq2PwjFOZBaaGuKy3wlgvUOiopgLJCtXZ8SAK07Q0d+0HnY6DjEdC6p8aOxIvp4FCXl1QG+XcpF2n81ksM+eknmDYNnn8evvtOb/J//jPcfjvk5akS8adlS7juOrjgAlUijz5aab2UlamiGTQI9t5bLaWZM+HII/UnSE6GDh301ayZ3si3b9f9bNmi7wB9+kC7drq+rMwLfaFWVHKyHispSZehUhmJ/dQ1IiILnHP9Yi1HbZjyMqKDq4BNP6pLb918zfIqyNFEg9KC3ftKksZrkptDUjPvla6fk9Mrl1NbQcu9oFU2tNkH0juotZVgj6XhsmCBuu2aN4d169QaSk+H0aOhpASOOgp++AH699fldu2gUycYPFi3//prdcu1b6/t6emqVHzKItkzUn0KJDGx0nrwKR2zJhoH8aK8GlPMy6iFRYsWAXDIIYfEWJJacE5TqNfOgzVfwJo5sH6hxpcAEtM0uaBVT+h0DDTvAOntoXknjRG16aWZYJJQ+YgsnjsPt7tbzgiJoiKYOhUWLYKcHI3tFBTAjBm6/o47NEbko00bOOIIHUskAn/9q7rajj5aFVOg5VKfeV9NaRl1wSyvOGLgwIFAI4t5bfsF1i2E/G8g/zvY9L3GoHzJDwnJ0HZ/yOwLWX0g61B17SWne3EoU0QNgXMahzrmGLWc/vIXmDgRUlIq3XEdO8KECbq8erXGV3wuva5dtW9Skr7MzbbnYJaXEXEeffTRWIugd8UtK+Dn9+CHlyDvy8p1zTtB671h33PUgmp7IHQ+Cpq3VyWWGGREpRFxFi2Cq66CuXM11jRwIPzud/p+7LGqrFJTVTmlpaly2nvvWEttGOFhyiuOiJm7sGwHLH8dVr0Hv/5PY1UAGT3gsD9DxyPV1ZfaUpVUakuNPSWmeANrjUhRXg7LlmmM6ptvoEcPuPpqtYx+9ztVXKtWaUzqgQc0s65VK+jcWdPK09Ji/Q0MIzLYnSWOiHq2YeE6WPQ4LH5CKzmkZECno6HPFdDpSGjZE9IztX23zD2jIbj4Ys34KyzUzykpcMop+hKB/HzYZx845xy46CI48MDKTDvDaGqY8oojbr75ZqCBY17OwaYfYP44LThaVgRdfgMHjlbFlZQGKS01ZmWZfQ1GRYVm8H3wAdxwA2zbptl8w4ZpirjvlZGhSsw5TbjwZfKlpHhD1gyjiWLKK4547LHHGvYAm5bBJzfDzzMAgX3Oht4XQ+aBflUdTFk1FMuWaVLFggWwcCFs2qSW06GH6kDbMWNUgZlSMgxTXnFFg5WF2rEVFj0GXz2k/qc+V8ABF0Kz9jreKq21uQQjSGGhWlTffqtjp66+WqtCzJ0L//wn9OqlyRVHHqkVKnr1UgvLlJZhVGLKK4744osvgAgW6C3aCIv+T2NaRevVPXjcQ1oOKa21V5PPgiZ1Yft2VVKlXiXxnj21ZNKNN2odvsJCfU7o0gWOO04TLA49VK2uHj00GzApycZAGUZ1mPKKI26//XYgAjEvVwEr3oIPr9aK6Z2OhoH/hPaH6lQX6e21/p8RFhUVOpj3xRe10KyPs8/W2nslJfDpp1oV/KyzKpWWbxxVcjWlEg3DqIoprzjiySefrP9Oykth9t0wfyxkdIMzpkLWwRrpb95Ri9XaHTQstm9Xt15xsaaq9+0LF16o9f0SEtSiSk/XJIr58ysL2hqGUXdMecUR++23X/12sOVnmHGBDize+0w49n5ISNFUd3MRhkVhIUyfDi+/rFXPP/5YldKECTqmqoU9AxhGg2LKK4745JNPABgwYED4Gy/7L3xwpc4RNfAf0PNMbc/oai7CMNi2TUstvfyy1gvs0AHOP19dhl26WGKFYUQLU15xxN/+9jcgzJhXaRF8dhss/BdkHgQn/B+06KiVL1p0spJNIeKrjr5xI7z5psaszjxTswJbttSSS2ZpGUb0MOUVR0ycODG8DXK/gPcu0VqE+58P/W+DlOZaxT0pSGlwA9DswOXLdZ6qZcs0yeK77+DDD9Wq+vRTrQWYYnrfMGKGKa84omfPnqF3XjUT3jxbZ/09ZRJ0Okqnm2/WzsZs1UBpKYwaBU89VdnWrp1OnLh+vc5V1b276X3DiDWmvOKIDz/8EIBBgwbV3HHtfHjjTM0mPO05SGunmYSpLaMgZfyxfTuMHQunngqZmfq+1146LcgBB2i9wObNNaXdMIzGgV2OccQDDzwA1KK8dm6F6b/XWYZPm6zTlLTobJmEATgHr74KU6aoG3DTJlViI0fqTMFDhphb0DAaM6a84ojnn3++5g7OwYyLYHsOnP4itOiiisvqEQKQm6tuv5ISuPtutbY6ddIyTJdcolPaZ2SYhWUY8YBdpnFEt27dau4w9xFY+TYccTt0OHyPV1wlJVqKaepU+Owz2LBBx2O1bavLd9wBV1yhn9PTTWkZRjxhl2sc8d577wFw6qmnVl2Z8xl8cRf0GAwHXLTHK66vv4bf/x5WrtQ41tFHw0EH6eDhbt10huHkZFVahmHEH1FVXiJyKjAeSASecc6NCVjfA5gIZAGbgIucczneukeAIUACMBP4k3PORVH8mDNmjJ6uKsqrKB/eOQ/SO8BxD2ol+ORmMZAwtvz6K+Tk6DxXoDMIP/kkDB1aOXjY6gcaRtMgaspLRBKBCcDJQA4wT0Tecs4t9es2DpjsnHtORE4ERgPDReQY4Figr9fvc2AAMCta8jcGpkyZEnzFe5dA0QatU9isHTRrG13BYkxRkU55P26cWlcvvaTJFtOnawUMcwcaRtMjmpf1EcAK59xKABGZApwF+Cuv3sBfvOWPgTe8ZQekASmAAMnAuijI3Kjo2LFj1cYfXtHJI/vfCm33g+ad95hxXJs3wyOPqHW1ebNWvRg5UqcR6dRJK1+YlWUYTZNoKq8uwGq/zznAkQF9FgPDUNfiUCBDRNo55+aIyDWLfWgAACAASURBVMdAHqq8HnPOfR94ABEZCYwE6N69e+S/QYyZPn06AGeccYY27NwGH98AbfeHAy7QtPg9pE6hL9X94Yd1wsZLL9W4VocOWqrJ6gsaRtOmsTlUbgIeE5ERwKdALlAuIvsABwBdvX4zReQ3zrnP/Dd2zj0FPAXQr1+/JhcP+/vf/w74Ka8590HRWjjx/6BZ5h4xCLmkBBYu1MrtRx8Nb7wBhx+u8a30dFNahrGnEE3llQv453p39dp24Zxbg1peiEgL4Bzn3BYRuRL40jlX4K17Fzga2E15NXWmTp1a+WHLzzoL8j5nQ1ZfjXU1cXJyNIPwm29gxgzNGjz5ZLW0DMPYs4jmc+o8oJeI7CUiKcB5wFv+HUQkU2RXwOY2NPMQ4FdggIgkiUgymqxRxW3Y1MnMzCQzM1M/fHIjkACH/cWbi6vploOoqIAxY7RU0+LFcP/90LMn9Ohhissw9lSiZnk558pEZBTwPpoqP9E5952I3AfMd869BQwERouIQ92G13mbTwVOBL5Fkzfec85Nj5bsjYVp06YBMKx/Jqx4HQ69XsdzpbWJsWQNR0UFXHwxvPiixrZuuUVT4Tt0MBehYezJSFMdKtWvXz83f/78WIsRUQYOHAg4Zl22DQrXwLB3oXUvSM2ItWgNQmkprFsHd92lsxSPHKlKq00byyI0jIZCRBY45/rFWo7aaGwJG0YNvPnmm/Dts/DVDXD8I5qkkdIi1mI1CAsWwJo10Ls3/O1vmqjRvTukpcVaMsMwGgPmeIkjWqU6Wi15GDL7wF5DIL19kzJBnNPK7i+9BCedBLffrq5BU1yGYQRiyitecI5XRl/MK7Pz4IhbNc6V1HTu5jt26Lito46CCy9U1+Bjj2nMq0cPU1yGYeyOuQ3jhYJcnnjpHUhpw7ntD9cMwzinrAy2bIGdO+GLL+C88zT9fcwYnU8rPV0/J9tUZIZhBGDKK16Y/w9mXFYBZz6jca44r6RRWAh5ebBxoyZhHH00jB6tJZ4SErQ2YdeuprgMwwiOuQ3jgaJ8WPIf0vc9hfSOB8a11VVSAqtXw88/w733wgknQH6+Kqxzz62sS5idbYrLMIzqMeUVDyx+Ekq28cKq3rzw2oy4ne6kqAh++UVdhTfcAC+8ABdcoO7B4mJN2MjOtoK6hmHUjimvxo5z8MOLkNWXZ16fzTMvTou1RHVi82adb+uDD2DwYPjf/zS2dd99OmVJSoomZqQ03UIhhmFEEIt5NXbyvoRN38Ox9zNz2jnQumesJQoL59QtuHGjzmK8dKnGuJ56Sgvqbt+ulpZVzDAMIxxMeTV2lkyChGTocQrJGR0gJX4SNcrLYf16WLJEU9779IEbb4TbblNFVVCglTPaN63haoZhRAFTXo2Z0h2w/DXocTKktuTZl6dBQhIjRoyItWQ1Ul6uimnDBrW0rrpKC+jOmqXjtZzT9c2bm+IyDKNumPJqzKx8C3Zsgl5DIak5z05+AaBRK6/yck3KKCuDzz6D66/XubaeflozCYuKtE+7dtC2rSkuwzDqhimvxsySSVq/sNNR0KwNs2bNirVEtZKfry7C55+HBx+Evn1h0iTIytL4Vtu26iq0xAzDMOqDhcgbKwVr4NePYJ+hkJgGSY0/PX7bNs0qTE/XebdOOw1ee02TMQoK9L19e1NchmHUH7O8GitLnoWKUtj7TEhtDZLA008/DcCVV14ZW9mCUFwMX3+tcayWLWHcOF12TpVa27Zar9AwDCMSmOXVGKkoh+8mQtYhOl+XN+3JK6+8wiuvvBJj4apSXAw//aQDj//wB413tWih7kN/i8swDCNSmOXVGPn1I9jyk87ZlZi8q3r8hx9+GGPBdqe8XK2q9evhkUfUVfj00zrouKxMkzO6dIGMpjlXpmEYMcSUV2Nk8ZOQ3AK6nwxp7WItTVDKyrRiRnk5TJ+uCRrXXAOnn65FdxMSdA6u9PRYS2oYRlPElFdjo3gTrHwb9j8PktMhpfmuVY8//jgA1157baykA9QduGaNxrNmz4abb4bjj4dbb1XFlZ4OHTtqarxhGEZDYDGvxsbyaVBRAj2HaKwrofL5Yvr06UyfPj2GwikbNugcXM2awW9+o7GuSZOgtFQHIXfqZIrLMIyGxZRXY+PHV6F5J8jsA6m7p+e9++67vPvuuzESTNm2rTI5IzdXFdjNN+v0JRUVanFZjULDMBoau800JnZshtWzYK8hkJi6K1GjsbBzp87Fdcst8NFHaoGBug8LCzWr0ObgMgwjGkRVeYnIqSKyTERWiMitQdb3EJGPROQbEZklIl391nUXkQ9E5HsRWSoi2dGUPSosf13HdmUP9sZ27V47afz48YwfPz4molVU6MzH48fD55/rrMeHHKKKa/t2LffUsmVMRDMMYw8kaspLRBKBCcBpQG/gfBHpHdBtHDDZOdcXuA8Y7bduMjDWOXcAcASwvuGljjKr3tNyUJl9NVkjgI8++oiPPvooBoLBpk3wxhuaCn/JJTrrcWmpKq727bX8k2EYRrSIZrbhEcAK59xKABGZApwFLPXr0xv4i7f8MfCG17c3kOScmwngnCuIltBRwznI+VTrGCYkqtswgLfeeisGgukg5PXr4dlnoX9/uOceTZUvKdEJJJs1/spVhmE0MaLpNuwCrPb7nOO1+bMYGOYtDwUyRKQdsC+wRUSmichCERnrWXK7ISIjRWS+iMzf4AvIxAublkHROlVeqS0bTbn18nJ1F6anwyuvqOWVnFw5ANkUl2EYsaCxJWzcBAwQkYXAACAXKEctxN946/sDPYERgRs7555yzvVzzvXLijc/1i8z9b1jfx2gHIRx48Yxbty4qInkHKxbp0qrvFwrZWRlacmnzEwbgGwYRuyIptswF+jm97mr17YL59waPMtLRFoA5zjntohIDrDIz+X4BnAU8J9oCB4VVn8M6R0ho0e1WYZz5syJqkhbtsDkyXD77arIRozQrMKWLTVBwzAMI1ZEU3nNA3qJyF6o0joPuMC/g4hkApuccxXAbcBEv21bi0iWc24DcCIwP2qSNzSuQuNdXX6jiRoS3CB+7bXXoiZSYaHOfHzffXDUUTB8uLoKU1I0Jb6ReDUNw9hDiZrb0DlXBowC3ge+B151zn0nIveJyJlet4HAMhH5EegAPOhtW466DD8SkW8BAZ6OluwNzoYlsGMjdDoSUmKfb15SAt9+WzkL8hNPaGZhYqLGuWwQsmEYsSaqtQ2dczOAGQFtd/stTwWmVrPtTKBvgwoYK371xbuOqHFg8pgxYwC49dYqQ+QihnOwdi2MHavxrjfeUBfhjh2Qna0V4w3DMGKN3YoaA79+DBndoGU3SKx+muFFixY1qBjOQX6+psbfeCOccAIcfLCO5erSxapnGIbReDDlFWsqyiH3M+gxWLMMawgmTZkypeHEqFBLKydH6xNmZKiltX27ZhbanFyGYTQmLHoRa9Z9DSXboOORu2ZMjgUbN6rVNXKkWl1gmYWGYTReTHnFml+92ZE7HxG0qoY/999/P/fff3/ERSgq0oklr7sOFi2CU09VSwy09JNlFhqG0dios9tQRFoToPycc5vqLdGexupZ0HIvaN4ZEmoOKi1btizihy8vh1Wr1OJavBj++U9VXgUFOiDZ5uUyDKMxEpbyEpEewL/RlHb/zAIBHGC3unCoKIc1X+jEk8nNazVxXnjhhYiLsHkzPPwwfP01/PvfcMYZmiqflKRp8oZhGI2RcC2vSUBr4HJgDaqwjLqybgGUFmiKfAziXTt3aqxr2DDo3FkVV2lpZcFdG89lGEZjJVzldQRwlHNuSUMIs8fx6//0vWO/WuNdAHffrUPi7rvvvogcft06TX/v319fpaWq0Lp3h9TaxTEMw4gZ4SqvnwG7rUWKnE+gZXZI8S6A1atX19onVIqK4N571VN5770a+9qxQxVXWuOawNkwDKMK4TqG/gSMFpF9GkKYPQpfvKvjESHFuwAmTZrEpEmT6n/oCq1b+NxzqrBA0+K7drUpTgzDiA/CtbzeRC2vZSKyEyjzX+mci31hvnhh/UJvfFe/qMe7Nm+GO+/UhIzbbtPMwg4doHnzqIphGIZRZ8JVXqMaRIo9EV+8q0P/kOJdALfddhsAo0ePrvNhS0rg8cdh4UL4xz8qFVabNnXepWEYRtQJS3k5555rKEH2OHI+0XqGGV1DincBbNy4sV6HLC+HlSth/Hg49lj4/e/V6ure3QYiG4YRX4Q9SFlEUoELgd5oqvx3wMvOuZ0Rlq3p4io03tX9pJDjXQBPPfVUvQ67fr2O35o2TWdB3rFDXYc2I7JhGPFGWAkbItIbWA78AzgSnc34UeBHETkg8uI1UTZ8Azu3qMswSvGu7dvhq680IWOffSrLPmVlReXwhmEYESXcbMPxwEKgu3PuN8653wDdgcWoEjNC4ZeP9L1T6PEugJtuuombbrop7MOVlcG778I558DkyZptWFQEnTrZ/FyGYcQn4d66jgX6O+e2+Rqcc9tE5A7gy4hK1pTJmQUtukBG95DjXQDFxcV1OtzatZpV2LmzVtPwZReau9AwjHglXOW1Ay0PFUgrb51RG64CcmdDtxPCincBTJgwIezDlZbCI49oosYLL2ih3YwMaB3sVzQMw4gTwnUbTgeeFpFjRSTRex0HPAm8FXnxmiD5S2DnZujQD1IafobHxYu14O5ZZ8Hxx2tbZqZlFxqGEd/UpcLGcuAz1NLaAXwC/AjcEFnRmii/ePN3dToSElNq7hvADTfcwA03hH6ai4vhu+/UXXjnnRrn6tDBpjkxDCP+CXec1xbgLBHpBezvNX/vnFsRccmaKjmfaC3DluHFu8KlvBzy8mDgQPj8cx2cnJwMLWI3WbNhGEbEqFOumXNuOWqBGeHgnMa7uhwXdrwL4NFHQ0/oLCiAt99Wd6GIKq8uXcxdaBhG06BW5SUi/wJuc84VesvV4py7vpZ9nYqm2ycCzzjnxgSs7wFMBLKATcBFzrkcv/UtgaXAG865+CtVtfF72LGxshhvA+EcTJkCf/qTDkQeOlTjXDbNiWEYTYVQLK8+QLLfcnXUODGliCQCE4CTgRxgnoi85Zxb6tdtHDDZOfeciJwIjAaG+62/H/g0BJkbJzmz9L1DaPN3BXLdddcBtWcdbtyo05wccICWgCottexCwzCaFrUqL+fcCcGW68ARwArn3EoAEZkCnIVaUj56A3/xlj8G3vCtEJHDgQ7Ae0C/esgRO3I+h7R20Co77GQNgGYhzlfyr39pvGvCBB2gnJlpg5ENw2ha1PuW5s3tleOcq22cVxfAfzbFHLTElD+LgWGoa3EokCEi7YDNwN+Bi4BBNcgyEhgJ0L179zC+RZTImwPtD4Ok9DoFn8aNG1fjeufg11/hiSfguOOgX7/K+oWGYRhNiXBrGz4kIpd4yyIiM9E0+TwROSoC8twEDBCRhcAAIBcoB64FZvjHv4LhnHvKOdfPOdcvq7EV7du+Bratgg6HNdj4ro0bYflyTcz48581Nb59e0uNNwyj6RGu5XUhcK63fBpwCFqc90I0PlWTWzEX6Ob3uavXtgvn3BrU8kJEWgDnOOe2iMjRwG9E5FqgBZAiIgXOuVvDlD925Hqhug6HQ1LdMidGjhwJBK8uv3075OdrnOuddzTOBdDSpgc1DKMJEu4g5Q6ouw/gdOBV59xc4P+AQ2vZdh7QS0T2EpEU4DwCqnKISKaI+GS6Dc08xDl3oXOuu3MuG7XOJseV4gLI/RwS06Bd7zrFuwDatWtHu3btqrTv3KmTS95xBxQWatuOHTog2VLjDcNoioSrvDYCPbzlwYBXHp0koMbbpHOuDJ2J+X3ge1TxfSci94nImV63gcAyEfkRVZQPhilf4yV3NmT1hZSWIOGedmX06NFVZlEuK4OcHHj4YZgxQ+fsKiqCtm11+hPDMIymSLhuw9eAlzzl0hZVRKDuw1qrbDjnZgAzAtru9lueCkytZR/PAs+GI3TMKSmA/G+hz5URn79r3Tr47DNVXDffDN26qUILYqAZhmE0GcJVXn8BfkHn8Pqrc85zUtEJeCKSgjUp1nwBrlzjXYlpdd7NpZdeCsCkSZMAtbBWr9a6hXvtBVdfrfUMu3e3JA3DMJo24dY2LENT1gPb/xkxiZoiOZ8BAu0PrXO8C6Bbt8p8F+d0nq5//ENdhW+8ofUMMzNtni7DMJo+oZSHOgxY5Jyr8JarxTn3dcQka0rkfg5t9oNm7SCh7ibRfffdt2u5oEAzCm+/HU49Ffr2VaurTZtICGwYhtG4CcXymg90BNZ7y47gyRkOrVlo+FNRBuvmwd5nRayeoXNqbaWm6sSSp56qWYbt2pm70DCMPYNQlNdewAa/ZSMcNnwDpYU6ODm5fv68iy66CIAnn3yB11/X4rsTJ1ZaWzamyzCMPYVQahv+EmzZCJHc2fqedQgk1D3eBbDffvtRUaGp8U88oZZX27ZqdXXqZPULDcPYcwjrdicio4AtzrkXAtovAlo65x6PpHBNgrw5Wow3ozsk1E+73HXXXeTlwYsvwsqV8PTTOkC5RQt1HxqGYewphDta9gZ2L67rYxXw53pL0xRZO08HJyfXrRivP0VFWr/w8cehd2+NdfmqxlslDcMw9iTCNQW6ouO8Asnx1hn+7NgCW1ZAz99Ccv0GJzsH5557Hjk58PPPU/jPf3R25ObNbZJJwzD2PMJVXmvRahqrAtoPA/IjIVCTYs0cfc86uF7juwC2bYN99z2E3r3h4ovhlFO0GG+nThGQ0zAMI84IV3m9BPxLRAqBWV7bCcCjwIsRlKtpkDcHEGh3UL2U144dOk/XyJG37optlZRo7UKrX2gYxp5IuDGvvwGz0ZqGRd7rXeAL4K7IitYEyPsKWvWE9LoPTi4vh9xcTYsfOBA2eIMWduyAxjZlmWEYRrQIS3k550qdc+cD+wEXeK/9nXPnOedKG0LAuMU5WLfAS9ao++DkLVu0DNT48VBefg63334OO3dqCSgrA2UYxp5KnXK3nXPLRWQbsME5VxFhmZoGW1fCjo0a70qqm5YpK9MJJu+5R5d///ujadNG0+OzsyMqrWEYRlwR7jivZHSOrWuAZsC+wEoReRj4xcZ5+ZH7hb7XI1ljyxaYOhU+/BDuvReuuOImtm3TSSbT6l6c3jAMI+6pS8zrDOAiYKdf+1xgRIRkahqs/VKnP2m7HyQkh715WZmO6fr8czjuOLjsMi28m5EBrVs3gLyGYRhxRLhuw/OBy5xzn4iIv7twCWqFGT7y5kLmgTq+qw4jiLdsgYQEeOYZTc4AuPLKM0lPh7fffivCwhqGYcQX4SqvzgQfpJxUh301XcpLIP8b2P/COs2cXF4OS5eqa7B5c02HLyiAk08+yZI0DMMwCF/hfAccT9VByn8AFkRCoCbBuq9VgbU/GBLDL3+xbRs88ADMnw/z5mnB3aQkuOWWP5EQrqPXMAyjCRKu8roXeEFEuqFzd/1eRPZHU+aHRFq4uCXPq6yReUjYyRoVFTor8syZcOutanVt2wbdu2OKyzAMwyPccV7TUStrMFCBJnD0As5wzn0YefHilDVfQrMsyOgSdiV5n9WVnQ0jR2qSRqtWOqbrtNNO47TTTmsYmQ3DMOKIkO+sIpKEKq2vnHMD6nIwETkVGI9abc8458YErO8BTASygE3ARc65HBE5BHgCaAmUAw86516piwxRYd38Og1Odk5T41eu1MrxKSk6V1dmpq4/44wzGkBYwzCM+CNk5eWcKxORacD+wMZwDyQiicAE4GS0Cv08EXnLObfUr9s4YLJz7jkROREYDQxHy1Bd7A2O7gwsEJH3nXNbwpWjwSnepAOU9z477MHJBQWwbBn06AGnn66Kq21bSPYy7a+99toGENgwDCP+CDeKshjYp47HOgJY4Zxb6ZwrAaYAZwX06Q38z1v+2LfeOfejc265t7wGWI9aZ40PX7yrfV9ICj1Zo6JC6xb+8Y/w8cf6OSkJ2rRpIDkNwzDimHCV1z3A30XkbBHpJiJt/V+1bNuF3SeyzPHa/FkMDPOWhwIZItLOv4OIHAGkAD+FKXt0WONVks/sE1ayRmEh/PyzKqzkZK0a36ULJPrV8x00aBCDBg2KvMyGYRhxRrjZhu9479MA59cu3ue6lU6v5CbgMREZAXwK5KIxLj2ISCfgeeCSYDUVRWQkMBKge/fu9RSljuR9Ba33gbS2IKE9GzgHX36pMyM/8giceaZWjA+cZPLcc89tAIENwzDij3CV1wn1OFYu0M3vc1evbReeS3AYgIi0AM7xxbVEpCWqPO9wzn0Z7ADOuaeApwD69evngvVpUHyV5LufGFayRnExPPqoDko++WR1Gfrm7fLnyiuvjKCwhmEY8UtIyktE0oGxwNlAMvAhcL1zLpzZk+cBvURkL1RpnYeOD/M/TiawybOqbkMzDxGRFOB1NJljahjHjC6bV8DOzV4l+dBniZw9G959F667rrKiRnL45RANwzD2GEKNed2LFt59B3gZzRh8IpwDOefKgFHoRJbfA686574TkftE5Eyv20BgmYj8CHRAK9iDji07HhghIou81yHhHD8qrJmt72FUki8uhrFjoUULuPpqKC2tPklj4MCBDBw4MDKyGoZhxDGhug2HAZc756YAiMiLwGwRSXTOlde8aSXOuRnAjIC2u/2WpwJVLCvn3AvAC6EeJ2bkfanp8W32DbmS/LJl8NVXcOWVanUlJKjlFYwRI0ZETlbDMIw4JlTl1Q34zPfBOTdXRMrQQr2rq91qT2NteJXki4s1zjVnjmYZ7twJ3bpV39+Ul2EYhhKq2zARKAloK8MqyVdStgPyl0BmX0gJLVljxQpVWq1ba4wrLa16qwugtLSU0tLSCAlsGIYRv4SqfAQtyOs/AWUa8LSIFPkanHNnVtlyT2HdQqgohfaH6CSUtVBSAtdco+O4XnpJra7s7JoNtpNPPhmAWbNmRUZmwzCMOCVU5fVckLbGH4OKJmu+0PfMviEla3zwgWYZ3nUXFBVp/cK0WnTeFVdcEQFBDcMw4p+QlJdz7tKGFiTuWfsVpHeEjK6QUPNY7fJyuP9+aN8eLrpIkzTa1lafBLjooosiJKxhGEZ8YzNERYq1cyGrT0gzJ7/9NsydC6NG6efMzNDm6ioqKqKoqKj2joZhGE0cU16RoGgDbPtFp0GppZJ8eTk8/TR07Ajnn69KK1g1jWCcfvrpnH766REQ2DAMI76xbMFIsMarJJ91aK3xru3bYcwY2LhRFVmHDqHPkHzNNdfUU1DDMIymgSmvSLDmC5BEaHdgjYOTy8ogN1dnRm7eXDMLW7UK/TBWmNcwDEMxt2EkWDtXK8k3a1djrvtrr8EJJ8BPP2lqfMeOIY1l3sXWrVvZunVrBAQ2DMOIb8zyqi+uQivJZ59SY7yrrEzdhamp6ips06b21PhAzjpL5+60cV6GYezpmPKqL5uWQck2yDoEkqvXRu+8A4sWwUMPhZ4aH8j1119fD0ENwzCaDqa86kuuNzg56xBICJ6sUVGhVlf79jrRZLt2dZvyZNiwYbV3MgzD2AOwmFd9yfsCUjKgzT6QEPxZYPZsnSn5mmtUaYWTpOFPfn4++fnhTKFmGIbRNDHLq76snQuZfXTm5CDZF85B587w/PNw8MHQsmXdJ5r83e9+B1jMyzAMw5RXfSgthI1Loe9InQYlCEVFmqxx4ok6xqu6iSZD4cYbb6z7xoZhGE0IU171IW+uZhtmHVbt4OSLL4YePbQUVIsWmm1YV84444y6b2wYhtGEsJhXfdhVSb5PUOW1cCFMm6auw5ISTdSoD2vXrmXt2rX124lhGEYTwCyv+pD3JWR0hxYdglaSHz0a0tPhwgs11hXuuK5AzjvvPMBiXoZhGKa86opzmqzR8UhN1ghg40Z4800tvtuiRf2tLoBbb721/jsxDMNoApjyqivbV0PRemh/cNDKGs8/r67CoUPV6qpPrMvHqaeeWv+dGIZhNAEs5lVX1vgNTg4S7zrwQPjTn2DffSNjdQGsXr2a1atXR2ZnhmEYcUxUlZeInCoiy0RkhYhU8YGJSA8R+UhEvhGRWSLS1W/dJSKy3HtdEk25g7JmDiSmQtv9q1SSLy2F7t1VeTVvHhmrC2D48OEMHz48MjszDMOIY6LmNhSRRGACcDKQA8wTkbecc0v9uo0DJjvnnhORE4HRwHARaQv8DegHOGCBt+3maMlfhbwvoV1vSGlZZXDyiy9qnKtvX8jKitwh77zzzsjtzDAMI46JpuV1BLDCObfSOVcCTAHOCujTG/ift/yx3/pTgJnOuU2ewpoJxC4AVF4CGxbrzMkByRolJXDTTfDMM6rT0mueWDksBg0axKBBgyK3Q8MwjDglmsqrC+AfsMnx2vxZDPiqzw4FMkSkXYjbIiIjRWS+iMzfsGFDxASvwvrFUL5T411Ju+e/v/GGZhoOG6Y1DEOdJTkUVq5cycqVKyO3Q8MwjDilsWUb3gQ8JiIjgE+BXKA81I2dc08BTwH069fPNYSAAKyZre9Zh1ZJ1vjPf9RVeOyxmmUYSS677DLAxnk1BioqKsjPz2fLli2Ul4f8FzWMmJKYmEjr1q3JzMwkIZJP1jEgmsorF+jm97mr17YL59waPMtLRFoA5zjntohILjAwYNtZDSlsjeR9Cc2yoEXn3SrJr1wJH30EI0dCs2b1H5QcyL333hvZHRp1JicnBxEhOzub5ORkJJwpsQ0jBjjnKC0tZd26deTk5NC9e/dYi1Qvoql65wG9RGQvEUkBzgPe8u8gIpki4pPpNmCit/w+MFhE2ohIG2Cw1xYb1s6FrIMhZfdivHPnQkYGnHNO5NLj/RkwYAADBgyI/I6NsCksLKRLly6kpKSY4jLiAhEhJSWFLl26UFhYGGtx6k3UlJdzrgwYhSqd74FXnXPfich9InKm120gsExEfgQ6AA96224C7kcV4Dzgkb2jdQAAHY1JREFUPq8t+hTlw9afVXn5DU4uK4NDD4Wvv4Zu3TRFPtIsW7aMZcuWRX7HRp2Id7eLsWfSVP63UY15OedmADMC2u72W54KTK1m24lUWmKxwzc4uf0hkFQ5gMs/j6JFC0hqgDN71VVXARbzMgzDaGwJG42fvDkgidDuoF2Dk8vL4YIL1GX4739D+/YNc+iHHnqoYXZsGIYRZzQN+zGa5H0FbXpBWqtdleRnzYIFC+D44yExUZM1GoJjjjmGY445pmF2bhgRIj8/HxEJy0Nwzz33cNBBBzWcUGGydu1aBg8eTPPmzXfFNLOzsxk3blyN27Vo0YJnn302ChIaprzCwVXA2vka7/JmTi4vh/Hj1VV49tnQunVkx3b5s2TJEpYsWdIwOzf2CEaMGIGIcPnll1dZd8sttyAi/Pa3v42BZKExa9Ysfvvb35KZmUmzZs3Yf//9+eMf/8iqVasiepxx48axZs0aFi1aRF5eHgDz5s3j2muvjehxjLpjyiscNv4Apdu9ZA01r+bNg3fe0Tm7mjVT12FDMWrUKEaNGtVwBzD2CLp168arr766W8ZZWVkZkydPbtTp008++SQnnXQS7dq147///S/ff/89//nPf6ioqOCBBx6I6LFWrFjB4YcfTq9evejYsSMAWVlZpEeyZI5RL0x5hcOuwclaSb6iAp54QsdzXX65loKKVBHeYIwdO5axY8c23AGMPYK+ffvSq1cvXn311V1t77zzDmlpaQwcOHC3vhUVFdx///1069aN1NRU+vTpw5tvvrlbn3nz5nH44YeTlpbGoYceyldffVXlmEuXLmXIkCFkZGTQvn17zj///LBmBc/JyeH666/nuuuu47nnnuOEE04gOzubY489lgkTJuzmzps2bRp9+vQhNTWVbt268eCDD+JcZc2C7OxsHnjgAa666ipatmxJ165dd7uusrOzefPNN5k8eTIiwogRI3a1+x9nxYoVDBw4kLS0NPbbbz/efvvtKnLn5uZy3nnn0aZNG9q0acOQIUNYvnz5rvU+d+mUKVPYe++9ycjI4OyzzyY/P3+3/Tz33HO7vlOHDh245JLK2uRbt25l5MiRtG/fnoyMDAYMGMD8+fNDPrfxiimvcMibo4V4W/eEhGQKC+Gvf4Vp07SaRkOM7fKnf//+9O/fv2EPYtSLgQOrvh5/XNcVFQVf7wuR5OcHX//KK7p+9eqq6+rK5ZdfzsSJlcm7EydO5NJLL60yZm38+PGMHTuWhx9+mG+//ZahQ4cybNgwFi1aBEBBQQFDhgyhZ8+ezJ8/nzFjxnDTTTftto+8vDyOP/54DjroIObOncuHH35IQUEBZ511FhUVFSHJ+9///peSkpJqJ2Rt3bo1AAsWLOD3v/89w4YN49tvv2XMmDGMHj2axx57bLf+//znP+nTpw9ff/01t9xyC3/961+ZM2cOoMp40KBB/OEPfyAvL4/x48dXOV5FRQVDhw6loqKCOXPmMHHiRO655x527ty5q09RUREnnHACaWlpfPLJJ8yZM4dOnToxaNAgioqKdvVbtWoVr7zyCq+//joffPABCxcu5I477ti1/sknn+Sqq67i0ksv5ZtvvmHGjBm74oPOOYYMGUJubi5vv/02Cxcu5Pjjj+fEE0/c5e5ssjjnmuTr8MMPdxFn0oHOvXSsc1tXu4oK51ascG7VKud+/dW5n35yrqIi8of0Z+HChW7hwoUNexAjJJYuXRq0fcCAqq8JE3RdYWHw9ZMm6foNG4KvnzJF1//6a9V14XLJJZe4IUOGuE2bNrm0tDT3448/ury8PJeSkuJ++eWXXet9dO7c2d17773/3969x+lY5g8c/3zHzJgDMwiDIccMsg7D6qQSemVThlQsRYvarE1tDsm21KbdvNYsKa9KCT/bDwmFH1psaxBCLDGOI6shtExY49R8f3/c90zPM8YYzDPPYb7v1+t5mfvw3Nf3nvvxfOe67uu+rnzneLf27t1bVVXfffddjY+P11OnTuVtnzFjhgL6+eefq6rqH/7wB23fvr3XMY4fP66Arl+/XlVVR48erTfffPNl4x44cKDGxcVd8fx69eql99xzj9e60aNHa2JiYt5y7dq1tWfPnl77NGjQQF999dW85c6dO2vfvn299qldu7b+5S9/UVXVzz77TMPCwvTAgQN521etWqWATnUv6JQpU7RBgwaa4/HFcPHiRa1UqZLOnj07L7ayZctqVlZW3j5jxozR+vXr5y0nJibqCy+8UOD5rlixQmNjY/XMmTNe65s3b65jx44t8D2ql//8qqoCGzUAvsOv9LKu8kV1/jQcT4dmv4bIWPbtc3oXpqZCcjJUr37JzCjF7rnnngPsOa9AVtiliYkpfHvlyoVvr1Wr8O1Xo2LFinTr1o0PPviAChUq0K5du0vud508eZJDhw5xxx13eK1v27Ytixc7j2ump6fTrFkzypX7abSZ2267zWv/TZs2kZaW5rVPrn379tGmTZsrxqtatKFK09PT6dy58yXxvvLKK5w8eZI4d8DRZs2aee1To0YNjh49WqQycstJTEz0+p3dcsstXg8Ab9q0if3791M+343wM2fOsG/fvrzl2rVrEx8fX2AsR48eJTMzkw4dOhQYx6ZNmzhz5gxV8s29dPbsWa8yQpElr6I6vM7pbVg1GcpEMXEifPcd1KwJERG+7aiRa8KECb4vxJQa/fr1o2/fvpQrV44//vGPV/XeqxkSKycnh86dOxfYzTwhIaFIx2jYsGFeMq1Ro0aRy/bkGXNERMQl24rahFlUOTk5tGjRglmzZl2yrVKlSsUSS05ODgkJCaxateqSbXHFPTJ4gLF7XkWV6XbWqNyMk/+NZMYMuO8+54HkqlV9X+sCaNGiBS1atPB9QaZU6NChA5GRkXz//fd07dr1ku1xcXHUqFGDNWvWeK1fvXo1TZo0AaBx48Zs27bNq+fiunXrvPZPTk5m+/bt1K5dmwYNGni98tdKLufhhx8mMjKS119/vcDtWVlZefEUFG/NmjWLXFZRNG7cmMzMTA4e/Gmmpi+//NIr6SQnJ7N3714qV658yXl7Jq/CVK1alcTERFasWFHg9uTkZI4cOUJYWNglZVT11WgJAcKSV1EdWgvx9SGmMh9MK0NWFgwY4NS6fDGOYUE2bNjAhg0bSqYwE/JEhK1bt7J//37KXqab7LBhwxg3bhwzZ85k9+7djBo1ilWrVuV1yujVqxfh4eH069eP7du3s2zZMl577TWvYwwaNIgffviBHj16sH79ejIyMli+fDlPPfUUp06dKlKstWrVYvz48bz11lv07duXf/7znxw4cIC1a9fyzDPPMGzYMACGDBnCypUrefnll9m9ezcffvghqampDB8+/Dp+U5fq2LEjjRo1ok+fPmzZsoW1a9fyu9/9jnCPceF69+5NQkICKSkprFy5kv3795OWlsaQIUO8ehxeye9//3smTJjA+PHj2b17N1u2bCE1NTUvjjvuuIOUlBSWLFnC/v37Wbt2LaNHjy6wNhZKLHkVhSoc+RKqNufHMuV46y1o2tR5Va5cMrUucL5Icv+TGlMcypcvX2jz0uDBgxk2bBjDhw+nadOmzJ8/n7lz59K8eXPAGVFi0aJF7Nmzh+TkZIYOHcrYsWO9jpFbewsLC6NTp07cfPPNDBo0iLJly142aRbkN7/5DcuWLePYsWN0796dpKSkvG7sL730EuDURObMmcPcuXNp2rQpI0aMYMSIEcX+fGRYWBjz588nJyeHW265hT59+vDSSy95nU9MTAxpaWnUq1ePRx55hEaNGtG3b19OnDhBxYoVi1zWwIEDmTRpEu+99x5NmzalU6dObN++HXD+AFm8eDHt27fnySefJCkpiUcffZRdu3Zdc/NqsJCi3ggNNq1bt9Zie9bh+G6YmgS3v8LpRs8z55NyxMRAq1ZQr57vRtTIL3d0jUAaRqe0Sk9Pp3Hjxv4Ow5hrUtjnV0Q2qWrrEg7pqlmHjaLIdKvfVZP5/kQk99wDFy5ApUoll7jAkpYxxuSyZsOiOPQFRJTjX5lNePmVCH74wWlJLOnOPF988QVffPFFyRZqjDEByJJXURxeD1WaMWFyVWbOEs6dg/h438zZVZiRI0cycuTIki3UGGMCkDUbXsn5U3A8nXONnmb2vFi6dnVqXO5oNCXq3XffLflCjTEmAFnyupLD60FzSNt7O9nZwkMPOdOf+HIA3stJSkoq+UKNMSYAWbPhlbgPJ7+z6C6qVVOaNXM6avjDypUrWblypX8KN8aYAGI1rys5tBaNr094uXh69BCionw3U/KVjB49GrCxDY0xxpJXYdyHky9Uu4fXXs4mJyqOSpVK7qHk/DynsDDGmNLMkldhTuyBsyc4cLYVZWNiyP7RGRncX+rVq+e/wo0xJoCU6D0vEekkIrtEZK+IXDKrnIjcKCKfi8hmEdkqIve76yNEZLqIbBORdBF5sUQCdh9O7ja0C1OmRxMd7Yxl6C/Lly9n+fLl/gvAGGMCRIklLxEpA0wCfgE0AX4pIk3y7fYS8JGqtgR6Au4ctDwClFXVnwGtgF+LSB2fB33oC85qHDuONOGWW8v4pXu8pzFjxjBmzBj/BmGClogU+sodJ9AXzp49i4iwaNGiK+771Vdf8cADD5CQkEBUVBS1a9fmkUceITMz02fxmeBTks2GbYC9qpoBICKzgBRgh8c+CuSOWxEPHPJYHysi4UA0cB446fOID69n86Gf0zjpIvUbRPq1yRBgxowZ/g3ABDXPaeEXLVrEk08+6bUu+jI9kS5cuHDJnFO+cujQITp06EBKSgpLly6lYsWKfPPNNyxcuJDTp0/7rNzz588TGRnps+Ob4leSzYaJwEGP5W/ddZ5eBh4TkW+BxcAz7vqPgf8Ch4F/A+NU9Xj+AkTkKRHZKCIbjx07dn3Rnj+F/iedz7a3pcsDF4mKKvkRNfKrVasWtWrV8m8QJmhVq1Yt71XBbUbwXBcfH8/OnTsREebMmcPdd99NVFQU06dP55133qFy5cpex1u6dCkikpdU/vOf/9CrVy+qVKlCVFQUDRo04O233wagTp06ADz44IOICI0aNSowxrS0NLKzs5kyZQotW7akTp06tGvXjtTUVK/nHA8ePEjPnj254YYbiImJoVWrVqxevTpv+5tvvkm9evWIjIykYcOGTJ8+PW9bbi1w8uTJdOnShdjY2LzJOLdt20anTp0oV64cCQkJPPbYY1z3d4nxiUDrsPFLYJqqporIbcAMEWmKU2v7EagBVARWicjy3FpcLlWdDEwGZ1T564rk8DqEHNZ+cxsvDg/ze5MhOF8WAJ06dfJzJKZAnz8HR7eUbJlVW8A9xT/D9ogRI0hNTaV58+aULVuWBQsWFOk9e/bsYcmSJVSuXJmMjIy8SSI3bNjAjTfeyIwZM+jYsaPXvFeeqlWrxrlz5/j000/p1q1bgTM2nzx5kjvvvJM6deqwYMECEhIS2LJlC7kzZMycOZOhQ4fyxhtv0L59exYtWkS/fv2oUaMG9957b95xRo0axeuvv86ECRMoU6YMBw8e5K677uK3v/0tb7zxBtnZ2YwYMYLu3buTlpZ2Lb9G40MlmbwyAc9qQ013naf+QCcAVV0rIlFAZaAXsFRVLwBHRWQN0BrIwFfch5N/Nbgm1WuEExXls5KKLHcWWUtexteef/75AmdXLsyBAwdo3bo1rVs7s2nk1rYAqlSpAkCFChWoVq3aZY/Rrl07hgwZwqOPPkpcXBxt2rShXbt29O7dO6/VYfr06WRlZfHJJ5/k1SAbNGiQd4xx48bRv39/nn766bxz+fLLLxk7dqxX8nr88ce97vMNHz6c22+/nVdffTVv3bRp06hevTpbt26lWbNmV/X7ML5VkslrA3CTiNTFSVo9cZKSp38DHYBpItIYiAKOuevb49TEYoFbgeL/c9NDTuZaLsbU56520ZwtG04gNIfPmjXL3yGYwvigBuQvuQnoagwaNIiePXuybt067r33Xrp06ULbtm2v+jjjxo3jhRde4B//+Afr1q3jnXfeYcyYMSxZsoQ777yTzZs306pVq7zEld/OnTt5/vnnvda1bduWP/3pT17r8p/jpk2bWLVqFeXKlbvkmPv27bPkFWBK7J6Xql4Efgt8BqTj9CrcLiJ/FJEu7m5DgCdF5F/ATOAJddoCJgHlRGQ7ThKcqqpbfRgs5w9sYOuxNpy5EEt8vP8eTPaUe2/CGF+LjY31Wg4LCyP/xLUXLlzwWk5JSeHAgQM8++yzHD58mPvuu4+BAwdeU/lVqlShR48ejB8/np07d1K9enVee+21azpWrvxNkPnPMScnh65du7Jlyxav1549e7xqbCYwlOg9L1VdjNMRw3PdKI+fdwB3FPC+0zjd5UvG8V1EyQn+/nVbHu4eS/nyJVZyoRYuXAg4N72NKUlVqlQhKyuLs2fPEuW2oW/Zcun9vapVq/LEE0/wxBNP0LFjRwYMGMCkSZMIDw9HRPjxxx+vuuyoqCjq1q2b1zGkZcuWzJ8/n6ysrAJrX40aNWLNmjX07t07b93q1atp0iT/kznekpOTWbp0KXXr1qVMmTJXHacpWTYwbwGObHbud5Wv9zOiYiMDoskQIDU1ldTUVH+HYUqh22+/ncjISF588UX27t3L7Nmzee+997z2GTlyJAsWLGDv3r1s376dTz75hKSkJMLCwggPD6dmzZosX76c7777Lq8jR35z586lb9++LF68mD179rBz507+/Oc/s2LFCrp16wZAnz59KF++PF27dmXNmjVkZGQwf/58Vq1yBhUYNmwYU6ZMYfLkyezZs4e//vWvzJkzh+HDhxd6jrk1xl69erFhwwYyMjL4+9//Tv/+/Tl//nwx/BZNsVLVkHy1atVKr9W//tJfT7war8vn7tNTJ3Ou+TjF7dixY3rs2DF/h2FUdceOHf4O4brMmTNHnf/+3tLT0xXQbdu2XbLto48+0vr162tUVJTef//9OnXqVAX01KlTqqo6atQobdSokUZHR2ulSpX0wQcf1F27duW9/+OPP9Z69eppeHi4JiUlFRjX7t27dcCAAdqwYUONiYnR+Ph4bdmypb755pte++3fv1+7d++ucXFxGhMTo61atdLVq1fnbZ84caLWrVtXIyIi9KabbtKpU6fmbcvOzlZAFy5cWOD5p6SkaHx8vEZHR2tSUpIOHjxYL168WPgvNMgU9vkFNmoAfIdf6SWq19ejPFC1bt1aN27ceE3vzRjdjMysBBJ/NZXaP6uJtSCY/NLT02ncuLG/wzDmmhT2+RWRTap69T12Spg1G+aTc/YUdcp/zcmyzYksHxtQiWvevHnMmzfP32EYY4zfBdpDyn4Xpue5kDycphdbUr5CADzc5WHixIkAPPTQQ36OxBhj/MuSV37RN3Dx56O5uDuT6NgA6anh+vTTT/0dgjHGBARLXpcRGRNFZFQAtRkC8fHx/g7BGGMCgt3zKkCZ8DKUrxQgD3d5mD17NrNnz/Z3GMYVqp2dTGgLlc+t1bwKEBkdSWR0YDUZAnkjdPfo0cPPkZiIiAiys7OJ8fc8OcZcpezs7BKb4saXLHkFkcWLF195J1MiqlatSmZmJomJiURHRxc4+rkxgURVyc7OJjMzk4SEBH+Hc90seQUR+ys/cMTFOXOmHjp06JIx/owJVBERESQkJOR9foOZJa8g8re//Q2Axx57zM+RGHASWCh8CRgTjCx5BZH3338fsORljDGWvILIsmXL/B2CMcYEBEteQSQUeggZY0xxsOe8gsi0adOYNm2av8Mwxhi/s+QVRCx5GWOMI2SnRBGRY8CB6zhEZeD7YgonWNg5lw52zqXHtZx3bVWt4otgilPIJq/rJSIbg2FOm+Jk51w62DmXHqF83tZsaIwxJuhY8jLGGBN0LHld3mR/B+AHds6lg51z6RGy5233vIwxxgQdq3kZY4wJOpa8jDHGBB1LXvmISCcR2SUie0VkhL/j8QURqSUin4vIDhHZLiLPuusricgyEdnj/lvR37EWNxEpIyKbRWSRu1xXRNa713u2iATeLKTXSUQqiMjHIrJTRNJF5LZQv9Yi8jv3s/21iMwUkahQu9Yi8oGIHBWRrz3WFXhdxTHRPfetIpLsv8iLhyUvDyJSBpgE/AJoAvxSRJr4NyqfuAgMUdUmwK3AIPc8RwArVPUmYIW7HGqeBdI9lscC41W1AXAC6O+XqHzrDWCpqjYCmuOcf8heaxFJBAYDrVW1KVAG6EnoXetpQKd86y53XX8B3OS+ngLeLqEYfcaSl7c2wF5VzVDV88AsIMXPMRU7VT2sql+5P5/C+TJLxDnX6e5u04Gu/onQN0SkJtAZeN9dFqA98LG7SyieczxwFzAFQFXPq2oWIX6tcQYdjxaRcCAGOEyIXWtVTQOO51t9ueuaAvyPOtYBFUSkeslE6huWvLwlAgc9lr9114UsEakDtATWAwmqetjd9B0Q/HOFe5sADAdy3OUbgCxVveguh+L1rgscA6a6zaXvi0gsIXytVTUTGAf8Gydp/QBsIvSvNVz+uobcd5slr1JMRMoBc4HnVPWk5zZ1nqEImecoROQB4KiqbvJ3LCUsHEgG3lbVlsB/yddEGILXuiJOTaMuUAOI5dLmtZAXatc1P0te3jKBWh7LNd11IUdEInAS14eqOs9dfSS3KcH996i/4vOBO4AuIvINTnNwe5x7QRXcpiUIzev9LfCtqq53lz/GSWahfK07AvtV9ZiqXgDm4Vz/UL/WcPnrGnLfbZa8vG0AbnJ7JUXi3ORd4OeYip17r2cKkK6qf/XYtADo6/7cF/i0pGPzFVV9UVVrqmodnOv6D1XtDXwOPOzuFlLnDKCq3wEHRSTJXdUB2EEIX2uc5sJbRSTG/aznnnNIX2vX5a7rAqCP2+vwVuAHj+bFoGQjbOQjIvfj3BspA3ygqq/5OaRiJyJtgVXANn66/zMS577XR8CNONPJPKqq+W8IBz0RaQcMVdUHRKQeTk2sErAZeExVz/kzvuImIi1wOqlEAhnAr3D+cA3Zay0irwA9cHrWbgYG4NzjCZlrLSIzgXY4054cAUYDn1DAdXWT+Fs4zadngF+p6kZ/xF1cLHkZY4wJOtZsaIwxJuhY8jLGGBN0LHkZY4wJOpa8jDHGBB1LXsYYY4KOJS9jgpiI1BERFZHWBS0bE6oseZlSRUSmuV/uKiIXRCRDRMa54/0FNBH5p4i8lW/1QaA6sMUPIRnjN+FX3sWYkLMceByIAO7EeYA3Fhh4tQcSkUh3BgK/UNUfcQZgNaZUsZqXKY3Oqep3qnpQVf8X+BB36ggRaSIi/ycip9yJ/maKSLXcN7o1t0Ui8oKIfIszdiAiEikifxKRAyJyzq3RDfZ4X1GP+6yIZIrICRGZKiIxuduBu3HmXsutOdYpSjPhlco2JhhZ8jIGsoEIdyDTNOBrnLndOgLlgE9FxPP/yt1AM5yhdjq466YDfYDngcY4Ex1mQd4AqUU57p1AU3d7D6AbzuSZuP+uBabiNBNWx3uKiwJdRdnGBBVrNjSlmoi0AXrhzDo7EPiXqr7gsb0PzoR/rYEv3dVngX654+KJyE04g/3+QlWXuvtkeBRT1OOeBJ52mwLTRWQOTnL8s6r+ICLngTPuYLu5x7nSKRa1bGOCiv3lZUqjTiJyWkTO4tRm0oBngFbAXe620yJymp9qN/U93v91vgFdW+IMcPz5Zcor6nF3uIkr1yGg6rWc4DWUbUxQsZqXKY3SgKeAC8Ahd84n3Ga0/wOGFvCeIx4///cqyyvqcS/k26Zc/x+YRS3bmKBiycuURmdUdW8B678CHgUO5Ca0ItqCkyTuAZYWsP1aj5vfeZypeq5GcZVtTECxZkNjfjIJiAdmi8gtIlJPRDqKyGQRKX+5N6nqbpw5lN4Xke7uZKZ3isjj13PcAnwDtHF7GFYuYoeL4irbmIBiycsYl6oewpkuPgenBrUd58v/nPsqTB/gf4GJwE5gGk7SuN7jehqHU/vaARzDmXDQl+dkTMCyySiNMcYEHat5GWOMCTqWvIwxxgQdS17GGGOCjiUvY4wxQceSlzHGmKBjycsYY0zQseRljDEm6FjyMsYYE3T+HyIOHoJhM4eQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "run_precision_plt(X, y, nfolds, percentiles, run_sc, plt_title=plt_title, \n",
    "                  plt_names=plt_names, predict_correct=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}