Performance_Matrix_From_Scratch_Without_Sklearn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing performance metrics in Pure Python without scikit-learn\n",
    "\n",
    "### All the .csv files used for this blog and JupyterNotebook are uploaded in [this Github Repo](https://github.com/rohan-paul/ML_Algo_From_Scratch/tree/main/Performance_Matrices_From_Scratch)\n",
    "\n",
    "Here I will compute following for our example dataset 5_a.csv\n",
    "\n",
    "\n",
    " - Confusion Matrix\n",
    "\n",
    " - F1 Score\n",
    "\n",
    " - AUC Score\n",
    "\n",
    "  Note 1: in this data we can see number of positive points >> number of negatives points\n",
    "\n",
    "  Note 3: we need to derive the class labels from given score\n",
    "  \n",
    "  ## $$y^{pred} = [0 \\text{ if } y_{score} < 0.5  \\text{ else }  1]$$\n",
    "\n",
    "- Compute Accuracy Score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10100, 2)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>proba</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.637387</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.635165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.766586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.724564</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.889199</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     y     proba\n",
       "0  1.0  0.637387\n",
       "1  1.0  0.635165\n",
       "2  1.0  0.766586\n",
       "3  1.0  0.724564\n",
       "4  1.0  0.889199"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "df_5_a = pd.read_csv('./5_a.csv')\n",
    "\n",
    "print(df_5_a.shape)\n",
    "df_5_a.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>proba</th>\n",
       "      <th>y_predicted</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.637387</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.635165</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.766586</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.724564</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.889199</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     y     proba  y_predicted\n",
       "0  1.0  0.637387          1.0\n",
       "1  1.0  0.635165          1.0\n",
       "2  1.0  0.766586          1.0\n",
       "3  1.0  0.724564          1.0\n",
       "4  1.0  0.889199          1.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_5_a['y_predicted'] = np.where(df_5_a['proba'] >= 0.5, float(1), float(0))\n",
    "df_5_a.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Checking to see if there's any 'proba' less than or equal to 0.5\n",
    "# And there is none. So all y_predicted will be classified as 1\n",
    "# df = df_5_a.loc[df_5_a['proba'] <= 0.5 ]\n",
    "# df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "actual_y_train_arr  [1. 1. 1. ... 1. 1. 1.]\n",
      "predicted_y_arr  [1. 1. 1. ... 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "# print(df_5_a.to_numpy())\n",
    "actual_y_train_arr = df_5_a.iloc[:, 0].values\n",
    "print('actual_y_train_arr ', actual_y_train_arr)\n",
    "predicted_y_arr = df_5_a.iloc[:, 2].values\n",
    "print('predicted_y_arr ', predicted_y_arr)\n",
    "y_proba = df_5_a.iloc[:, 1].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "#### In below implementation, for a binary class-label (1 and 0 ) I will have 'true-nagative' at the top left of the final confusion matrix, as the traversing of the unique_classes array will start from 0\n",
    "\n",
    "![img](https://i.imgur.com/MvapmPv.png[/img])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Explanations on the above calculate_confusion_matrix() funcion\n",
    "\n",
    "[In above implementation, for a binary class-label (1 and 0 ) I will have 'true-nagative' at the top left of the final confusion matrix, as the traversing of the unique_classes array will start from 0 ]\n",
    "\n",
    "#### 1.np.zeros()  - the first arg of np.zeros() is the shape which is a tuple of ints or simple int\n",
    "\n",
    "e.g., (2, 3) or 2.\n",
    "\n",
    "---\n",
    "\n",
    "#### 2. Now for each row, I need to compare the values between true_y_classes_array and predicted_y_classes_array\n",
    "\n",
    "So, I will implement this by comparing each element from true_y_classes_array with each of the unique array's elements and then the same for  predicted_y_classes_array\n",
    "\n",
    "It will give me a list of True / False. Lets see this example\n",
    "\n",
    "```python\n",
    "arr1 = [1, 1, 2, 6, 4]\n",
    "arr2 = np.unique(arr1) # [1 2 4 6]\n",
    "\n",
    "print(\"comparing \",  arr1 == arr2[0])\n",
    "# comparing [ True  True False False False]\n",
    "```\n",
    "\n",
    "That is, the above will compare the whole of arr1 with arr2[0] which is 1\n",
    "\n",
    "---\n",
    "\n",
    "### 3. An now use np.sum() to count number of True value\n",
    "\n",
    "```python\n",
    "print('np.sum of above True False numpy array ', np.sum(arr1 == arr2[0]))\n",
    "\n",
    "# np.sum of above True False numpy array  2\n",
    "```\n",
    "\n",
    "Note that above kind of operation will only work with `np.unique()` which returns a numpy array, and will not work with plain python list\n",
    "\n",
    "---\n",
    "\n",
    "### 4. So now for the confusion matrix, I need to fill up with matching counts of\n",
    "\n",
    "#### Note `unique_classes[0]` is 1 and `unique_classes[1]` = 0\n",
    "\n",
    "### For first row of my final confusion_matrix\n",
    "\n",
    "confusion_matrix[0,0] => i.e. i, j = 0, 0 => will have the Total 'True' count (i.e. `np.sum()`)  of following conditions\n",
    "\n",
    "`(true_y_classes_array == unique_classes[0]) & (predicted_y_classes_array == unique_classes[0])`\n",
    "\n",
    "Similarly for `confusion_matrix[0, 1]` => i.e. i, j = 0, 1  => will have the Total 'True' count (i.e. `np.sum()`)  of following conditions\n",
    "\n",
    "`(true_y_classes_array == unique_classes[0]) & (predicted_y_classes_array == unique_classes[1])`\n",
    "\n",
    "#### And for second row of my final confusion_matrix\n",
    "\n",
    "confusion_matrix[1,0] => will have the Total 'True' count (i.e. `np.sum()`)  of following conditions\n",
    "\n",
    "`(true_y_classes_array == unique_classes[1]) & (predicted_y_classes_array == unique_classes[0])`\n",
    "\n",
    "Similarly for confusion_matrix[1, 1] => will have the Total 'True' count (i.e. `np.sum()`)  of following conditions\n",
    "\n",
    "`(true_y_classes_array == unique_classes[1]) & (predicted_y_classes_array == unique_classes[1])`\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "---\n",
    "\n",
    "# Confusion matrix of binary class-labels placing 'true-positive' / TF at the top left of Matrix\n",
    "\n",
    "### If for a binary class-label I have to place 'true-positive' / TF at the top left of the final confusion matrix\n",
    "\n",
    "i.e. the top-left will have 1 instead of 0 in the final confusion matrix then only need to reverse the unique_classes variable as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1.]\n",
      " [0. 3.]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def calculate_confusion_matrix(true_y_classes_array, predicted_y_classes_array):\n",
    "\n",
    "  # extract all unique classes from the train y class\n",
    "  unique_classes = np.unique(true_y_classes_array)\n",
    "  # print('unique', unique_classes)\n",
    "\n",
    "  # initialize a matrix with zero values that will be the final confusion matrix\n",
    "  confusion_matrix = np.zeros((len(unique_classes), len(unique_classes)))\n",
    "\n",
    "  for i in range(len(unique_classes)):\n",
    "    for j in range(len(unique_classes)):\n",
    "      confusion_matrix[i, j] = np.sum((true_y_classes_array == unique_classes[i]) & (predicted_y_classes_array == unique_classes[j]))\n",
    "\n",
    "  return confusion_matrix\n",
    "\n",
    "# actual_y_class_list = [1, 3, 3, 2, 5, 5, 3, 2, 1, 4, 3, 2, 1, 1, 2]\n",
    "# predicted_y_class_list = [1, 2, 3, 4, 2, 3, 3, 2, 1, 2, 3, 1, 5, 1, 1]\n",
    "\n",
    "# Only binary class dataset\n",
    "actual_y_class_list = [1, 1, 1, 0, 0]\n",
    "predicted_y_class_list =   [1, 1, 1, 0, 1]\n",
    "\n",
    "print(calculate_confusion_matrix(actual_y_class_list, predicted_y_class_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[10000.   100.]\n",
      " [    0.     0.]]\n"
     ]
    }
   ],
   "source": [
    "def get_confusion_matrix(true_y_classes_array, predicted_y_classes_array):\n",
    "  \n",
    "  unique_classes = np.unique(true_y_classes_array)\n",
    "  # For a binary class the above will give me [0 1] numpy array\n",
    "  # so top-left of confusion matrix will start from 0 i.e. 'True Negative'\n",
    "\n",
    "  # But the challenge here asks that the top left will be 'True Positive'\n",
    "  # Hence I need to reverse the above numpy array\n",
    "  unique_classes = unique_classes[::-1]\n",
    "  # print('reversed unique', unique_classes) # will convert the above array to [1 0]\n",
    "\n",
    "  # initialize a matrix with zero values that will be the final confusion matrix\n",
    "  # For the binary class-label dataset, this confusion matrix will be a 2*2 square matrix\n",
    "  confusion_matrix = np.zeros((len(unique_classes), len(unique_classes)))\n",
    "\n",
    "  for i in range(len(unique_classes)):\n",
    "    for j in range(len(unique_classes)):\n",
    "      confusion_matrix[i, j] = np.sum((true_y_classes_array == unique_classes[j]) & (predicted_y_classes_array == unique_classes[i]))\n",
    "\n",
    "  return confusion_matrix\n",
    "\n",
    "confusion_matrix_5_a = get_confusion_matrix(actual_y_train_arr, predicted_y_arr)\n",
    "print(confusion_matrix_5_a)\n",
    "\n",
    "true_negative, false_positive, false_negative, true_positive = int(confusion_matrix_5_a[1][1]), int(confusion_matrix_5_a[0][1]), int(confusion_matrix_5_a[1][0]), int(confusion_matrix_5_a[0][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explanations and notes on above Confusion matrix function\n",
    "\n",
    "![img](https://i.imgur.com/1A3Izpg.png)\n",
    "\n",
    "#### Note `unique_classes[0]` is 1 and `unique_classes[1]` = 0\n",
    "\n",
    "### For first row of my final confusion_matrix\n",
    "\n",
    "`confusion_matrix[0,0]` => i.e. i, j = 0, 0 => will have the Total 'True' count (i.e. `np.sum()`) of following conditions\n",
    "\n",
    "`(true_y_classes_array == unique_classes[0]) & (predicted_y_classes_array == unique_classes[0])`\n",
    "\n",
    "Similarly for `confusion_matrix[0, 1]` => i.e. i, j = 0, 1 => will have the Total 'True' count (i.e. `np.sum()`)  of following conditions\n",
    "\n",
    "`(true_y_classes_array == unique_classes[1]) & (predicted_y_classes_array == unique_classes[0])`\n",
    "\n",
    "---\n",
    "\n",
    "### Now second row\n",
    "\n",
    "And for second row of my final confusion_matrix\n",
    "\n",
    "`confusion_matrix[1,0]`  => i.e. i, j = 1, 0 => will have the Total 'True' count (i.e. `np.sum()`)  of following conditions\n",
    "\n",
    "`(true_y_classes_array == unique_classes[0]) & (predicted_y_classes_array == unique_classes[1])`\n",
    "\n",
    "Similarly for `confusion_matrix[1, 1]`  => i.e. i, j = 1, 1  => will have the Total 'True' count (i.e. `np.sum()`) of following conditions\n",
    "\n",
    "\n",
    "`(true_y_classes_array == unique_classes[1]) & (predicted_y_classes_array == unique_classes[1])`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "# To check that the total num of elements of the original dataframe matches\n",
    "# with the counts captured in the confusion matrix\n",
    "# sum-all-the-elements-of-the confusion_matrix_5_a\n",
    "sum_all_elements_of_confusion_matrix = np.concatenate(confusion_matrix_5_a).sum()\n",
    "print(sum_all_elements_of_confusion_matrix == df_5_a.shape[0] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 1]\n",
      " [0 3]]\n"
     ]
    }
   ],
   "source": [
    "# Testing my custom confusion_matrix result with scikit-learn\n",
    "from sklearn.metrics import confusion_matrix\n",
    "sklearn_confustion_matrix = confusion_matrix(actual_y_class_list, predicted_y_class_list)\n",
    "print(sklearn_confustion_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 1 0 3\n",
      "0 100 0 10000\n"
     ]
    }
   ],
   "source": [
    "tn, fp, fn, tp = confusion_matrix(actual_y_class_list, predicted_y_class_list).ravel()\n",
    "print(tn, fp, fn, tp)\n",
    "print(true_negative, false_positive, false_negative, true_positive)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### From above we can see the values of the confution Matrix matches between scikit-learn and our custom-implementation\n",
    "\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "## F1 Score\n",
    "\n",
    "![img](https://i.imgur.com/ZPntYB0.jpg)\n",
    "\n",
    "![Imgur](https://imgur.com/qy5Fesd.jpg)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "My custom function's f1_score and accuracy_score  (0.9950248756218906, 0.9900990099009901)\n"
     ]
    }
   ],
   "source": [
    "# the below function will work only for\n",
    "# binary confusion matrix\n",
    "def get_f1_score_accuracy_score(binary_conf_matrix):\n",
    "    true_negative  = binary_conf_matrix[1][1]    \n",
    "    false_positive = binary_conf_matrix[0][1]\n",
    "    false_negative = binary_conf_matrix[1][0]\n",
    "    true_positive = binary_conf_matrix[0][0]\n",
    "\n",
    "    precision = true_positive / (true_positive + false_positive)\n",
    "    recall = true_positive/ (true_positive + false_negative)\n",
    "    \n",
    "    f1_score = (2 * (precision * recall)) / (precision + recall )\n",
    "    \n",
    "    sum_all_elements_of_confusion_matrix = np.concatenate(binary_conf_matrix).sum()\n",
    "    \n",
    "    accuracy_score = (true_positive + true_negative)/sum_all_elements_of_confusion_matrix\n",
    "    \n",
    "    return f1_score, accuracy_score\n",
    "\n",
    "\n",
    "print(\"My custom function's f1_score and accuracy_score \", get_f1_score_accuracy_score(confusion_matrix_5_a))  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now verifying the above F1-Score with that of sk-learn\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sklearn_f1_score  0.9950248756218906\n",
      "sklearn_accuracy_score  0.9900990099009901\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import f1_score\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "sklearn_f1_score = f1_score(actual_y_train_arr, predicted_y_arr)\n",
    "print('sklearn_f1_score ', sklearn_f1_score)\n",
    "\n",
    "sklearn_accuracy_score = accuracy_score(actual_y_train_arr, predicted_y_arr)\n",
    "print('sklearn_accuracy_score ', sklearn_accuracy_score)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From above we can see the values of1_score, accuracy_score matches between scikit-learn and our custom-implementation\n",
    "\n",
    "# Note on ROC Curve\n",
    "\n",
    "The ROC curve is plotted with TPR against the FPR where TPR is on y-axis and FPR is on the x-axis.\n",
    "\n",
    "In Machine Learning, performance measurement is an essential task. So when it comes to a classification problem, we can count on an AUC - ROC Curve. When we need to check or visualize the performance of the multi - class classification problem, we use AUC (Area Under The Curve) ROC (Receiver Operating Characteristics) curve. It is one of the most important evaluation metrics for checking any classification model’s performance. It is also written as AUROC (Area Under the Receiver Operating Characteristics)\n",
    "\n",
    "### What is AUC - ROC Curve?\n",
    "\n",
    "AUC - ROC curve is a performance measurement for classification problem at various thresholds settings. ROC is a probability curve and AUC represents degree or measure of separability. It tells how much model is capable of distinguishing between classes. Higher the AUC, better the model is at predicting 0s as 0s and 1s as 1s. By analogy, Higher the AUC, better the model is at distinguishing between patients with disease and no disease.\n",
    "The ROC curve is plotted with TPR against the FPR where TPR is on y-axis and FPR is on the x-axis.\n",
    "\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "# AUC-ROC Score Calculation\n",
    "\n",
    "The Receiver Operating Characetristic (ROC) curve is a graphical plot that allows us to assess the performance of binary classifiers. With imbalanced datasets, the Area Under the Curve (AUC) score is calculated from ROC and is a very useful metric in imbalanced datasets.\n",
    "\n",
    "TPR and FPR are defined as follows:\n",
    "\n",
    "- TPR = True Positives / All Positives\n",
    "- FPR = False Positives / All negatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "def get_single_tpr_fpr(df):\n",
    "\n",
    "    '''\n",
    "    Note, this implementation is only for binaly class labels (0 and 1)\n",
    "    :param df: the dataframe should have 'y' and 'y_predicted' as its labels\n",
    "    :return: a list containing tpr and fpr\n",
    "    '''\n",
    "\n",
    "    tp = ((df['y'] == 1.0 ) & (df['y_predicted'] == 1)).sum()\n",
    "    fp = ((df['y'] == 0.0 ) & (df['y_predicted'] == 1)).sum()\n",
    "    tn = ((df['y'] == 0.0 ) & (df['y_predicted'] == 0)).sum()\n",
    "    fn = ((df['y'] == 1.0 ) & (df['y_predicted'] == 0)).sum()\n",
    "\n",
    "    tpr = tp / (tp + fn )\n",
    "    fpr = fp / (fp + tn)\n",
    "\n",
    "    return [tpr, fpr]\n",
    "\n",
    "\n",
    "# While computing AUC score you need to calculate \"TP,\"FP\" at every threshold by using actual \"y\" and predicted \"y_pred\".\n",
    "\n",
    "def calculate_all_thresholds_tpr_fpr_arr(df_original):\n",
    "\n",
    "    '''\n",
    "\n",
    "    :param df_original: the original dataframe, which should have a 'proba' label\n",
    "    :return: two arrays, tpr_arr_for_all_thresholds, fpr_arr_for_all_thresholds\n",
    "    '''\n",
    "\n",
    "    tpr_arr_for_all_thresholds = []\n",
    "    fpr_arr_for_all_thresholds = []\n",
    "\n",
    "    sorted_df = df_original.sort_values(by=['proba'], ascending=False)\n",
    "\n",
    "    unique_probability_thresholds = sorted_df['proba'].unique()\n",
    "\n",
    "    for threshold in tqdm(unique_probability_thresholds):\n",
    "        sorted_df['y_predicted'] = np.where(sorted_df['proba'] >= threshold, 1, 0)\n",
    "        tpr_fpr_arr = get_single_tpr_fpr(sorted_df)\n",
    "        tpr_arr_for_all_thresholds.append(tpr_fpr_arr[0])\n",
    "        fpr_arr_for_all_thresholds.append(tpr_fpr_arr[1])\n",
    "\n",
    "    return tpr_arr_for_all_thresholds, fpr_arr_for_all_thresholds\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now finally the AUC Score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10100/10100 [00:15<00:00, 664.72it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "My Custom function ROC-AUC Score for 5_a.csv:  0.48829900000000004\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "all_tpr_together_5_a, all_fpr_together_5_a = calculate_all_thresholds_tpr_fpr_arr(df_5_a)\n",
    "auc_score_5_a = np.trapz(all_tpr_together_5_a, all_fpr_together_5_a)\n",
    "print('My Custom function ROC-AUC Score for 5_a.csv: ', auc_score_5_a)\n",
    "# My Custom function ROC-AUC Score for 5_a.csv:  0.48829900000000004"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What exactly is np.trapz() function\n",
    "\n",
    "numpy.trapz() function integrate along the given axis using the composite trapezoidal rule. Say, I have one list of 100 numbers as height for Y axis, and as length for X axis: 1 to 100 with a constant step of 5. I need to calculate the Area that it is included by the curve of the (x,y) points, and the X axis. This is where I can use trapz() function\n",
    "\n",
    "```py\n",
    "import numpy as np\n",
    "\n",
    "# The y values.  A numpy array is used here,\n",
    "# but a python list could also be used.\n",
    "y = np.array([5, 20, 4, 18, 19, 18, 7, 4])\n",
    "\n",
    "# Compute the area using the composite trapezoidal rule.\n",
    "area = np.trapz(y, dx=5)\n",
    "print(\"area =\", area) # area = 452.5\n",
    "\n",
    "```\n",
    "\n",
    "### What is Composite Trapezoidal Rule\n",
    "\n",
    "The trapezoidal rule gives us a technique to approximate the integral on a given interval [a, b], but we cannot reduce the error because the error depends on the width of the interval over which we are integrating.\n",
    "\n",
    "By dividing the interval [a, b] into many smaller intervals, and applying the trapezoidal rule to each, this allows us to find a better approximation the integral.\n",
    "\n",
    "An animation that shows what the trapezoidal rule is and how the error in approximation decreases as the step size decreases\n",
    "\n",
    "\n",
    "![Img](https://upload.wikimedia.org/wikipedia/commons/thumb/1/10/WikiTrap.gif/220px-WikiTrap.gif)\n",
    "\n",
    "[Source](https://en.wikipedia.org/wiki/Trapezoidal_rule)\n",
    "\n",
    "Read more [here](https://secure.math.ubc.ca/~pwalls/math-python/integration/trapezoid-rule/) and [here](https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/13Integration/comptrap/complete.html) on Composite Trapezoidal Rule\n",
    "\n",
    "When passed two arguments (y and x both), trapz approximates\n",
    "\n",
    "![Imgur](https://imgur.com/Cxr9X3o.png)\n",
    "\n",
    "And when you do not pass x, then trapz assumes $x(t) = t$ meaning $dx/dt = 1$\n",
    "\n",
    "namely, the area between a curve and the x axis. \n",
    "\n",
    "An alternative for numpy `trapz()` function is `simps()` function from scipy. `simps()` function is used to get the integration of y(x) using samples along the axis and composite simpson’s rule by using scipy.integrate.simps() method.\n",
    "\n",
    "```py\n",
    "\n",
    "import numpy as np\n",
    "from scipy.integrate import simps\n",
    "\n",
    "y = np.array([5, 20, 4, 18, 20, 18, 3, 4])\n",
    "output = np.trapz(y)\n",
    "\n",
    "print(\"trapz() Output : \",output)\n",
    "\n",
    "# simps\n",
    "\n",
    "arr = np.array([5, 20, 4, 18, 20, 18, 3, 4])\n",
    "\n",
    "output = simps(arr)\n",
    "print(\"simps() output =\", output)\n",
    "\n",
    "# trapz() Output :  87.5\n",
    "# simps() output = 88.66666666666666\n",
    "\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sk-learn roc_auc_score for 5_a.csv:  0.48829900000000004\n"
     ]
    }
   ],
   "source": [
    "# Checking result with scikit-learn\n",
    "sklearn_roc_auc_score = roc_auc_score(actual_y_train_arr, y_proba)\n",
    "print('sk-learn roc_auc_score for 5_a.csv: ', sklearn_roc_auc_score)\n",
    "# sk-learn roc_auc_score:  0.48829900000000004"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## As we can see above the roc_auc_score matches between my Custom function and that of sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'AUC=0.4883')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxH0lEQVR4nO3de5yN5d7H8c9vGMaZSEKaDsQ4lcah2pVSe6eTnupply12xpnoIBt5CCkqJI3DmKGyq4lKKaSUsJXDlAghUShyyPk0ZuZ6/riXZoyZMbZZa83M+r5fr/Va9+Fas363GfOd677udd3mnENEREJXWLALEBGR4FIQiIiEOAWBiEiIUxCIiIQ4BYGISIhTEIiIhDgFgYhIiFMQSIFmZl+a2V4zK55pW4dM7Zqb2bYM62ZmPc1stZkdNrNtZjbdzOqf5fubmY0wsz2+xwgzs1y8brKZOTO7PMO2SDOb7TueHWb2qpkVzbD/ZjP71swOmNkmM+uUYd9NZva9me3z1THDzKqdzbFI6FIQSIFlZpHA9YAD7j7Ll48BegE9gfOAWsAHwB1n+XU6AfcADYEGwF1A55xeYGZ/AS7LYtc4YCdwIXAlcCPQzfeacGAGMBEoB/wdGGVmDX2vXQv8zTlXHqgK/AiMP8tjkRClIJCCrC2wBHgNaJfbF5lZTaA78JBz7gvn3HHn3BHn3JvOueFnWUM7YKRzbptz7ldgJPDPHN67KDAWeDSL3ZcA05xzx5xzO4BPgLq+fecBZYGpzrMc+AGIAnDO/e6c+y3D10oFLkckFxQEUpC1Bd70Pf5mZhfk8nUtgG3OuWXZNTCzvr7TLFk+MjStC6zMsL6S9F/eWXkcWOicW5XFvpeBB82spO+0Tku8MMA59zvwNvCImRUxs2uAi4H/ZKi5hq+2o0Bv4IUc6hD5k4JACiTf6ZWL8f6C/gb4CWidy5dXBLbn1MA5N9w5Vz67R4ampYH9Gdb3A6WzGicws4vwThsNzOZtF+KFyAFgG5CEd7rqpLd9rz0OLAKeds5tzVDzFl9tlYABwLqcjlHkJAWBFFTtgE+dc7t962+RfnooBQjP1D4cOOFb3oN3Hj4vHMI7ZXNSWeCQy3o2x5eBIc65/Zl3mFkY3l//7wOl8H6ZVwBG+PbXBhLxekHF8AKjj5mdNqbhnPsDeB34MONgs0h2FARS4JhZCeAB4Ebf1TU78E65NPQNnm4BIjO97BLgF9/y50B1M4vO4T36m9mh7B4Zmq7BGyg+qaFvW1ZaAC9mqBngazNrjTcGUAN41TdmsQeYAtzua1cP2OCcm+ucS3POrQdm4Z0+ykpRoDKnhpRIlhQEUhDdgzcYGoV3dc2VQB280yVtgXfwzqU38V3eWQsvKBIBnHM/4l2h87bvstJiZhZhZg+aWV9fm+ecc6Wze2So5Q3gCTOrZmZVgSfxBq+zUgsvKE7WDN5VRjN8PZvNQFczK2pm5fF6OCfHElYANX2XkJqZXQbceXK/md1rZleYWZiZnQ+MAlb4egciOXPO6aFHgXrgnUIZmcX2B4AdeH8Nt8f7y/wAsBHoC4RlaGt4l4+uAY4Av+IFSN2zrMXwBmX/8D1eACzD/kPA9dm81gGXZ1i/EvgS2AvsBqYBF2Q6vtXAQbwxhBEnjwnvKqTNwGHfv0EicHGwv1d6FIyHOacb04iIhDKdGhIRCXEKAhGREKcgEBEJcQoCEZEQV+A+bFKpUiUXGRkZ7DJERAqUb775Zrdz7vys9hW4IIiMjCQpKSnYZYiIFChm9kt2+3RqSEQkxCkIRERCnIJARCTEKQhEREKcgkBEJMT5LQh8N+feaWars9lvZvaKmW00s1Vm1shftYiISPb82SN4Dbgth/0tgZq+Ryd0o20RkWwlJyfjr0lC/RYEzrmFeNPyZqcV8IbzLAHKm1le3TVKRKRwSE5macuWNKpQgX/fe69f3iKYYwTVgK0Z1rf5tp3GzDqZWZKZJe3atSsgxYmIBNvhw4d5ol07rvnkE9YcOULswoV+6RUUiMFi51yccy7aORd9/vlZfkJaRKRQ+eKLL2jQoAGjExMxoA8wf948zCzP3yuYQfArcFGG9eq+bSIiIWvfvn107NiRFi1asGnTJhpceCFLgRHnnUeJq67yy3sGMwhmAm19Vw81A/Y757YHsR4RkaD68MMPiYqKIj4+nmLFijF06FCSLriAaIAHH/Tb+/pt0jkzextoDlQys23AICAcwDk3AZgN3I53P9kjwCP+qkVEJD/buXMnPXv25J133gGgWbNmJCQkEBUVBVOneo1KlfLb+/stCJxzD51hvwO6++v9RUTyO+ccb775Jr169eKPP/6gZMmSPPfcc/To0YMi/frB7Nnw889e43bt/FZHgZuGWkSkMNi6dStdunRh9uzZANxyyy3ExcVxySWXeA1eeglOXiFUujRUy/KiyjxRIK4aEhEpLNLS0hg/fjx169Zl9uzZlCtXjoSEBD799NP0EHjvvfQQWLUKtmyB8uX9VpN6BCIiAbJhwwY6duzIwoULAbjnnnuIjY2latWqpzZ844305Xr1wA+XjGakHoGIiJ+lpKTwwgsv0LBhQxYuXEjlypWZPn0677///ukhsH49zJzpLb/zjt9DANQjEBHxq5UrV9K+fXu+/fZbANq2bcuoUaOoWLFi1i/49NP05YYNA1ChgkBExC+OHz/Os88+y/Dhw0lJSaFGjRpMnDiR2267DZKS4Lrr4ODB01946JD33KYNXHFFQGpVEIiI5LGvv/6amJgYfvjhBwC6d+/O888/T5kyZbwGc+Z4p4ByctNNfq4ynYJARCSPHDp0iAEDBvDKK6/gnKNWrVrEx8dz/fXXew1++gnuvBN++cVbf+QRGDr09C8UEQHZnTryAwWBiEge+Oyzz+jUqRM///wzRYoUoU+fPgwcOJCIiIj0RosWwbp13nKRInD33X79fEBuKQhERM7B3r17efLJJ5kyZQoAV155JQkJCTRq1AieffbUD4YlJ3vPt98OiYlw8lRRkCkIRET+SzNmzKBbt27s2LGD4sWLM2jQIHr37k14eLjXYPp02L//1BeZQatW+SYEAMxftz7zl+joaJeUlBTsMkQkFMXFQe/e7EhO5tETJ3g3LQ2Aa81ICA+ndlimj2YdO+Y9L1iQfilo0aJ+nUAuO2b2jXMuOqt96hGIiOSS++gjph48yGPAXqAUMBzo5hxhJ0/7ZFajBlx1Vb7qAWSmIBARyYVffvmFzkuWMNe3/tdbbmHi2LFEXnxxzi8sVswbGM7HFAQiEtpSU3PcnZaWxrjx4+nbvz+HDx+mAjC6d2/avvCCX24bGQwKAhEJXaNGwVNPge9cf2brgRhgsW/9fmAsUKVly4DMARQomnRORELXvHnpIRAW9ufjhBnPm9EQLwQuAN4zY3pYGFUiIwM2B1CgKAhEJHTNmeM9f/yxd4ooNZUVSUk0vfJK+jvHceCRRx7hhz/+4N60NK/N5s0B/dRvICgIRCR0HDvmXdd/8nFyUrdSpTh27Bj9+/encePGrFixgosvvpi5c+cyefJkKlSoENy6/UxjBCISGhYvhltuSb+2P4P/bN1Khy5dWL9+PWZGz549GTZsGKVLlw5CoYGnIBCRwik1Ff74I319wYL0EChbFoCDztGvRAli27YFoHbt2iQkJHDttdcGutqg0qkhESl8nIPGjaFy5fTH0097+x59FPbvZ+60adSrUIHYnTspWrQoTz/9NCtWrAi5EAD1CESksElJga1bYcUKb71SpfR9JUrwxw038Hi7drzhuy9wo0aNSEhI4Morrwx8rfmEgkBECg/noEmT9BAID4ddu/7c/e6779K9e3d27txJ8eLFGTx4ME8++SRFi4b2r8LQPnoRKVicg23bvL/6s9t/MgSqVoV77gFg+/bt9OjRg/fffx+A66+/nvj4eGrVqhWAovM/BYGIFBxPPw3PP5+7tr/+inOO16ZM4YknnmDfvn2ULl2aESNG0KVLF8IyzxQawhQEIlIwHDoE06alr0dGZt/2jjvYvHkznTp1Yt68eQC0bNmSCRMmUKNGDf/WWQApCEQk/zt2DC67DHbu9NZnzfLu8pWF1NRUYmNj6VevHkeOHOG8885jzJgx/OMf/yg0k8TlNQWBiORPKSne/X2dg99/Tw+B5s3hmmuyfMkPP/xATEwMX3/9NQAPPPAAY8eOpXLlygEqumBSEIhI/vTQQ/Duu6duq1YN5s8/remJEyd44YUXGDJkCMnJyVx44YWMGzeOe3yDxZIzBYGI5B/HjsGaNd7yyat/AOrV857btTvtJd988w3t27dn1apVAMTExPDSSy9Rvnx5PxdbePg1CMzsNmAMUASId84Nz7S/BvA6UN7Xpq9zbrY/axKRfOxvf4OFC0/d9v336UGQwdGjR3nmmWcYOXIkqampXHLJJUyaNIkWLVoEqNjCw29BYGZFgFjgVmAbsNzMZjrn1mZoNgCY5pwbb2ZRwGwg0l81iUg+tnNneghccIF3GuiKK6BOndOaLly4kA4dOvDjjz9iZjz++OMMHTqUUkG4KXxh4M8eQRNgo3NuE4CZJQKtgIxB4ICyvuVywG9+rEdE8qvDh+Hyy9PXV6yACy88rdmBAwfo27cv48ePByAqKoqEhASaNWsWqEoLJX8GQTVga4b1bUDTTG2eAT41s0eBUsAtWX0hM+sEdAJ0DbBIYXLggPdLf+dOOHjQ29a9O1SpclrT2bNn06VLF7Zu3UrRokXp378//fv3p3jx4gEuuvAJ9mDxQ8BrzrmRZnYNMNXM6jnnTrmBqHMuDogDiI6OdkGoU0T84eab4Ztv0tcvvhheffWUJrt37+bxxx/n3//+NwDR0dEkJCTQoEGDQFZaqPkzCH4FLsqwXt23LaMY4DYA59zXZhYBVAJ2+rEuEQmWtDRYssQ7FQSwcaP3XLu2N1X0ww//2dQ5x/Tp0+nRowe7du0iIiKCoUOH8thjj4X8JHF5zZ//msuBmmZ2CV4APAi0ztRmC9ACeM3M6gARwC5EpHCKjYWePU/f/p//nHIf4N9++41u3brx4YcfAnDjjTcSHx/P5RnHESTP+C0InHMpZtYDmIt3aehk59waMxsCJDnnZgJPApPM7HG8geN/Oud06keksNqwIX35Ft+QYNOmf4aAc46EhAR69+7N/v37KVOmDC+99BIdOnTQJHF+5Nf+le8zAbMzbRuYYXktcJ0/axCRfGTuXO958GAYOPCUXZs2baJjx4588cUXANxxxx1MmDCB6tWrB7rKkKOIFRH/O3IEPv4YTl7nn+HS0NTUVEaPHk29evX44osvqFSpEm+99RYfffSRQiBANOIiIv7Xty+MHZu+fumlAKxZs4aYmBiWLl0KwEMPPcSYMWM4//zzg1FlyFIQiIh/7N7tTRCXlgbLlqVv796d5MaNGT5kCM8++ywnTpygWrVqjB8/nrvuuit49YYwBYGI+EfbtjBnzqnbZsxgebVqtL/uOlavXg1A586dGTFiBOXKlQtCkQIKAhHxl5MhcNllcPXVHDnvPAZ+8QWjY2NJS0vjsssuY9KkSdx0003BrVMUBCLiBxkvE/3gA77cvZsOHTrw008/ERYWRu/evRk8eDAlS5YMXo3yJwWBiJydLVvg88+9O4dlx/eJ4f1An1deIW7SJADq1atHQkICTZo0CUChklsKAhE5O61bw+LFZ2z2MdClWDF+nTSJ8PBwBgwYQN++fSlWrJj/a5SzoiAQkdxzLj0ErrnGmyMok13HjtFr6VLe3rQJkpNp0qQJCQkJ1Mvi5jKSPygIRCT3vvoqffmttyAy8s9V5xyJiYn07NmT3bt3U6JECYYNG0bPnj0pUqRI4GuVXFMQiEjubdmSvpwhBLZt20bXrl35+OOPAbj55puZNGkSl/o+OCb5m6aYEJHcO9kjuOEGANLS0pg4cSJRUVF8/PHHlC1blkmTJjFv3jyFQAGiHoGIZG3XLkhMhOPH07f98IP3XLUqGzdupGPHjnz55ZcA3H333YwbN45q1aoFvlY5JwoCEcnas8/CK6+ctjkFePnQIf6vfn2OHTvG+eefz9ixY3nggQcws8DXKedMQSAiWZs3z3uuWBH++U8Avt+1i5jPP2e5byygTZs2jB49mkqVKgWpSMkLGiMQkaytXes9v/wyx4cNY1CpUjR66y2W//or1atXZ9asWUydOlUhUAioRyAi6dasgffe82YMjYiAY8dYUq4cMY0asdYXDF27dmX48OGULVs2yMVKXlEQiEi6Rx/1po4GDgP/FxbGy61a4ZyjZs2axMfHc4PviiEpPHRqSERg/34YPvzPq4I+b9KE+uXLMzotDTOjT58+rFy5UiFQSKlHICIwdSr068c+4Ckg3ncjmQYNGpCQkEB0dHQwqxM/U49AJNS99hpMm8aHQBQQDxQrVoyhQ4eSlJSkEAgB6hGIhLLffuP3Rx6hJzDNt+maa64hPj6eqKioYFYmAaQgEAlFc+bgZs/mzRUr6AX8AZQsVoznhwyhe+/emiQuxCgIRELQlrZt6bJ7NyfvKHxLmTLErVzJJZdcEtS6JDg0RiASKpYtI61rV8Zffz11fSFQPiKCya1b8+ny5QqBEKYegUiI2NCrFx2WLGGRb/1/wsKIXbuWCxUAIU89ApHCLCmJlA4deKFxYxr6QqByRATTO3bkva++UggIoB6BSKG28sknab9wId/61tsBI+fPp2KzZsEsS/IZ9QhECqFjx44xYMAAon0hUKNkST55/HFe+/xzhYCcRj0CkULmq6++IiYmhnXr1gHQA3juk08oc/31wS1M8i2/9gjM7DYzW29mG82sbzZtHjCztWa2xsze8mc9IoXZoUOH6NWlC3+57jrWrVvHFaVKsQgYC5T5y1+CXZ7kY37rEZhZESAWuBXYBiw3s5nOubUZ2tQE+gHXOef2mlllf9UjUph99tlndOrUiZ9//pkiQB9g4OHDRADUrw+6c5jkwJ+nhpoAG51zmwDMLBFoBazN0KYjEOuc2wvgnNvpx3pECp29e/fy5JNPMmXKFACuLFKEhNRUGpUo4c0hBHDddcErUAoEf54aqgZszbC+zbcto1pALTNbbGZLzOy2rL6QmXUysyQzS9q1a5efyhUpWN5//32ioqKYMmUKxYsX57lu3ViWmkojgF694IEHvIduJi9nEOzB4qJATaA5UB1YaGb1nXP7MjZyzsUBcQDR0dEuwDWKBMa6dTBgABw5kmOzHceP02PNGt77/XcArqtQgfh69aj97bfpjQYN8melUsj4Mwh+BS7KsF7dty2jbcBS59wJYLOZbcALhuV+rEskf3r9de82kdlwwBvA48BeoBQwHOi2dy9hixalN7zvPu82kyK55M8gWA7UNLNL8ALgQaB1pjYfAA8BU8ysEt6pok1+rEkk/4mNhRkz4McfvfXmzeGpp05p8svOnXQeO5a5vr/6/9aoERMffZSLK2e6viIsDHSFkJwlvwWBcy7FzHoAc4EiwGTn3BozGwIkOedm+vb91czWAqnAU865Pf6qSSRfGjAA9u1LX+/UCW6/HYC0tDTGjRtH3759OXz4MBUqVGD06NG0bdsW05VAkkf8OkbgnJsNzM60bWCGZQc84XuIhAbnoHt3WL3aWz9wwHt+7z246CLw3RFs3bp1dOjQgcWLFwNw//33M3bsWKpUqRKMqqUQC/ZgsUjo+flnGD/+1G2VKsFdd0F4OCdOnODFF19k8ODBJCcnU6VKFWJjY7n33nuDUq4UfgoCEX/ZuhU6d4Y//jh1+7Fj6csLFnjPtWtDeDgrVqygffv2fPfddwA88sgjjBw5kgoVKgSmZglJCgIRf5k9G+bMyX7/LbfADTcA3iRxg/v148UXXyQ1NZXIyEji4uK49dZbA1SshDIFgYg/pKZCly7e8pVXnn4qyAwaNgTgP//5DzExMWzYsAEzo2fPngwbNozSpUsHtmYJWQoCEX/wzfwJQJs2kMXUzwcPHqRf797ExsYCULt2bRISErj22msDVaUIoPsRiOS9LVv+POVDWBg8cfpFcZ988gn16tUjNjaWokWLMmDAAFasWKEQkKA46yAwszAz+4c/ihEpFD77LH2A+H/+55SZP/fs2UO7du1o2bIlW7Zs4eqrryYpKYmhQ4cSoU8DS5BkGwRmVtbM+pnZq2b2V/M8ivfJ3wcCV6JIAXHsGNx2G/Tv7603awaJiQA453j33XeJiorijTfeICIighEjRrBkyRIa+sYKRIIlpzGCqXhTmnwNdAD6Awbc45z7zv+liRQwq1bB3Lnp6w88AEWLsn37drp3786MGTMAuP7664mPj6dWrVpBKlTkVDkFwaXOufoAZhYPbAdqOOeO5fAakdA0ahS89JK3XL48rFyJu+giXpsyhSeeeIJ9+/ZRunRpXnjhBTp37kxYmIbnJP/IKQhOnFxwzqWa2TaFgEg2Jk+G7du95f/9XzanptLpr39l3rx5ALRs2ZIJEyZQo0aNIBYpkrWcgqChmR3AOx0EUCLDunPOlfV7dSL53euvw//9H/z2GwCpM2fy6k8/0b9ePY4cOULFihV5+eWX+cc//qFJ4iTfyjYInHNFAlmISIGUmOhNJQGsLV+eDs8+y9fLlgHw97//nVdeeYXKmaeKFslncrpqKMLMHvNdNdTJzPThM5EsnACevecerjpyhK+XLePCCy/kgw8+IDExUSEgBUJOv9xfx/sZXwTcDtQFegWiKJGC4pv9+2kPrPrgAwA6dOjAiy++SPny5YNZlshZySkIojJcNZQALAtMSSL52PbtcNNNHP3tN545fpyXkpNJAy6pUoVJ//43LVq0CHaFImctt1cNpWigSwRYtoyF69fTAfgR78qJx4sXZ+j8+ZSqXTvIxYn8d3IKgit9VwmB9/Ouq4Yk9GzcCDffDLt2ccA5+p44wcl5RKNq1yZh7Fia3XADFCsW1DJFzkVOQbDSOXdVwCoRyY+WLoWtW5kNdAa24f2nefqvf6XfzJkUL148uPWJ5IGcgsAFrAqR/MY5aN6c3QsX8hjwpm9zdKNGTI6Lo/7VVwexOJG8lVMQVDazbG8q75wb5Yd6RIIjNfWUVbdvH9MWLuRRYBcQER7Os88/T69evShaVFdSS+GS0090EaA06Z8sFimcevWCV175c/U3oCsw07d+4403Eh8fz+WXXx6M6kT8Lqcg2O6cGxKwSkQCITUVkpNP3ea7r7ADEszo7Rz7gTLAS02b0uGLLzRJnBRqOf10qycghcuePVC9OpQseerjxx/5CWjRtCkdfSFw5513snbrVjotWaIQkEIvpx6BPhkjBd+BA97AL0BSEuzY4S377gaW6hxjUlMZkJrK0aVLqVSpEq+88goPPvigJomTkJHtnzrOuT8CWYhInuvSBcqV8+4PUL68d/cwgGuvhaNHWb18Odc2bMiTKSkcdY7WrVuzdu1aHnroIYWAhBRd/iAF18GD3u0hs/Pll+nLZX2ffwwLI7lVK54fPJhhw4Zx4sQJqlWrxvjx47nrrrv8Wq5IfqUgkILpo4+8G8NnuuwzS99/D/XqAbB8+XLat2/P6tWrAejcuTMjRoygXLly/qxWJF9TEEjBkpwMu3bB/PnpIVCpUvbta9eGWrU4cuQIAwcOZPTo0aSlpXHZZZcxadIkbrrppsDULZKPKQik4DhxAq64An7+OX3boEHwzDM5vmz+/Pl06NCBTZs2ERYWRu/evRk8eDAlS5b0a7kiBYWCQAqOvXvTQ6BqVShTBm6/Pdvm+/fvp0+fPsTFxQFQv359EhISaNy4cQCKFSk4/HqBtJndZmbrzWyjmfXNod19ZubMLNqf9Ughcf758OuvsG4dNGmSZZOPPvqIqKgo4uLiCA8PZ/DgwSQlJSkERLLgtx6BmRUBYoFb8SZtXG5mM51zazO1K4N357Ol/qpFQseuXbvo1asXb7/9NgBNmzYlISGBunXrBrkykfzLnz2CJsBG59wm51wykAi0yqLdUGAEkMN1gCLATN/sP2lpp+1yzvHWW29Rp04d3n77bUqUKMGoUaNYvHixQkDkDPw5RlAN2JphfRvQNGMDM2sEXOScm2VmT2X3hcysE9AJoEaNGn4oVfKt337zpoYA+O477zkl5ZQm27Zto2vXrnz88ccA3HzzzUyaNIlLL700gIWKFFxBGyw2szBgFPDPM7V1zsUBcQDR0dG6T0KoWL4cmjZNnyLipEcfBSAtLY1Jkybx1FNPcfDgQcqVK8fIkSNp3769Phkschb8GQS/AhdlWK/u23ZSGaAe8KXvP20VYKaZ3e2cS/JjXVJQzJqVHgK+D4RRujTcdx8//vgjHTt2ZMGCBQC0atWKcePGUbVq1SAVK1Jw+XOMYDlQ08wuMbNiwIOkT/GOc26/c66Scy7SORcJLAEUAuKZOhUGD/aWW7f2Ph38/fekLFrES/Pm0aBBAxYsWMD555/PO++8w4wZMxQCIv8lv/UInHMpZtYDmIt3k5vJzrk1ZjYESHLOzcz5K0hI+v132LTJ++TwSY88AsCqVauIiYkhKcn7W6FNmza8/PLLVKxYMRiVihQafh0jcM7NBmZn2jYwm7bN/VmLFACHDsHll3vPJw0ZwvHrr+e5QYN47rnnSElJoXr16kycOJHbc/gwmYjknj5ZLMG3YQNs3+49ToZAs2ZQpgxLatYkplEj1q71Pn7StWtXhg8fTtmTs4mKyDlTEEhwff89NGhw6rZLL+XwvHkMGDCAMa1b45yjZs2axMfHc8MNNwSnTpFCTEEgwfXLL+nLvl/ynzduTMf69dm8eTNFihShd+/eDBo0iBIlSgSpSJHCTTdjleA6+Snhu+5i34cf0qFmTW4ZOZLNmzfTsGFDli5dyvDhwxUCIn6kHoEEj3OwbBkAH/7+O12joti+fTvFihVj4MCB9OnTh/Dw8CAXKVL4KQgkeJYs4fdhw+gJTPMFwjXXXENCQgJ16tQJbm0iIUSnhiQonHNMffNNooBpQMmICMaMGcOiRYsUAiIBph6BBNyWLVvo0qULc+bMAeDWSpWIW76cyMjI4BYmEqLUI5CASUtLY9y4cdStU4c5c+ZQvlgxpgBz69dXCIgEkXoEEhAbNmygQ4cOLFq0CID/AWKTk7kQsr3LmIgEhoJA/ColJYWRI0cyaNAgjh8/zgWlShF7+DD3ATzwAJQsCR06BLtMkZCmIBC/WblyJe3bt+fbb78FoN0FFzDq9985D+DWW+Gdd4Jan4h4NEYgee7YsWMM6N+f6Kuv5ttvv6XGeefxyWOP8VpYmBcCV18NzzwT5CpF5CT1CCRPffXVV8TExLBu3ToM6AE898cflHn55fRGiYneLKMiki8oCOTs7d0LM2bA8eN/bjp07Bj9Z87k1QULcM5xRdmyxB84wF8A2rdPf23t2nDZZQEvWUSypyCQszdsGIwc+efqp0An4Be8OxD1AQYeOEAEQMeOEBcXjCpFJJcUBHJmO3fCu+/CiRPeuu8S0L3h4Txx6aW8tn49AFdWrMjkm27iqkqVvHbFi0P37sGoWETOgoJAzmzAAJg06ZRN7wPdS5Vix/r1FC9enEGDBtG7d29NEidSACkIJHsnTsCbb8KSJd56lSrsuOMOenz5Je/99BPs28d1111HfHw8tWvXDm6tIvJfUxBI9mbP/vPG8Q54/cEHeeL119m7dy+lS5dm+PDhdO3albAwXYUsUpApCCRr334LU6cC8DPQuVYtPvVdAvq3v/2NiRMncvHFFwetPBHJOwoCyVqbNqT98AOxQL+iRTm8YQMVKlTg5Zdf5uGHH8bMgl2hiOQRBYF4VwVNmQJHj/65ad0vv9ABWAyQksL999/Pq6++ygUXXBCsKkXETxQEAqNHw/DhAJwAXgQGA8lAlQsuIHbcOO69994gFigi/qQgEIiPB+DbcuWIKVGC73bsAKB9ixa8NH06FSpUCGZ1IuJnutwj1O3cydHdu+kHNDl4kO927CAyMpLPPvuMhHnzFAIiIUBBEMp27uQ/bdtyJTAcSHOOXr168f3333PLLbcEuTgRCRSdGgpRBw8epN/ddxO7dCkAdYoXJ2H+fK655pogVyYigaYeQajZu5dPWremXrVqxC5dSlFgQEQEK77+WiEgEqLUIwghe/bs4Yk77uANXy/gaiABaDhqFFx1VVBrE5Hg8WuPwMxuM7P1ZrbRzPpmsf8JM1trZqvM7HMz00dV/cA5x/Tp04m64greWLqUCOAFM5aMGkXDKVOgXbtglygiQeS3HoGZFQFigVuBbcByM5vpnFubodkKINo5d8TMugIvAH/3V02haPv27XTr1o0PPvgAgBuASUCtIUPg8ceDWZqI5BP+7BE0ATY65zY555KBRKBVxgbOufnOuSO+1SVAdT/WE1Lc0aNMvuMO6kRG8sEHH1AmPJzxlSszH6h13XXw2GPBLlFE8gl/BkE1YGuG9W2+bdmJAeZktcPMOplZkpkl7dq1Kw9LLJw2b97MX6+9lpjZs9mfnExLYM2JE3TZudP7hj/xBJQuHeQqRSS/yBeDxWbWBogGbsxqv3MuDogDiI6OdgEsrUBJ3b2bV++7j/6LF3MkNZWKwBigdXx8+iRxFSvCHXcEs0wRyWf8GQS/AhdlWK/u23YKM7sFeBq40Tl3PPN+yZ21a9cSc/fdLPnpJ8AbaHkFqNypE8TEBLU2Ecnf/BkEy4GaZnYJXgA8CLTO2MDMrgImArc553b6sZaC6/XXYdasbHcnp6UxYt06nv3hB5LT0qgKjAsPp9XkyRAeDrfdFrhaRaRA8lsQOOdSzKwHMBcoAkx2zq0xsyFAknNuJt5El6WB6b5TF1ucc3f7q6YCqWdPOHAgy11JeAMrq3zrHfEuuyrfpw+0aROY+kSkwPPrGIFzbjYwO9O2gRmWNaFNRs7BoEGwalX6tkOHvOcpU6BkSQCOJiczaNo0Rn78MWnOcekFFzCpUydurlcPIiLg1luDULyIFFT5YrBYfLZuhaFDT99eoYL3F37RoixYsIAOHTqwceNGwsLCeOLxxxkyZAilSpUKfL0iUigoCPKTY8fSl2fMSF9u0IADR47wr3/9iwkTJgAQFRXF5MmTadq0aYCLFJHCRkGQn6SkeM9RUXDPPX9unjVrFl26dGHbtm0ULVqUp59+mn79+lG8ePHg1CkihYqCIJiWLIEhQyA52Vs/OR5QpAgAu3fv5rHHHuPNN98EoHHjxiQkJFC/fv1gVCsihZSCIJgmTIA5p3+Y2tWowTuJiTz66KPs3r2bEiVKMHToUB577DGK+EJCRCSvKAiCZe9e7zMCAG3bwsMPA/Drnj10mzqVmQ89BEDz5s2ZNGkSl19+ebAqFZFCTjemCZYvvkhfbt8e16IFkzZvJqpTJ2bOmkXZsmWZOHEin3/+uUJARPxKPYK8tHUrdO4Mf/xx5rZ79njP9erxU/XqdGzRgvnz5wNw5513Mn78eKpX12SsIuJ/CoK8NGtWluf8s5MKjDn/fAbUr8/Ro0epVKkSY8eO5e9//3v6JHEiIn6mIDhXTz4Jvr/kOTlF9lVXwbhxOb5s9U8/EfP88yzzvbZ169aMGTOGSpUq+bNaEZHTKAjOxbFjMGrU6dvbtIFmzbJ8SXJyMs8//zzDhg3jxIkTVKtWjQkTJnDnnXf6uVgRkawpCM7FG2+kL3/zjfdcogTUrp1l82XLlhETE8Pq1asB6Ny5MyNGjKBcuXL+rlREJFsKgnNxMgjCw6FRo2ybHTlyhIEDBzJ69GjS0tK47LLLiI+Pp3nz5oGpU0QkB7p89L8xaBDUqgVJSd767NnZNp0/fz7169dn5MiRAPTu3ZtVq1YpBEQk31CP4L8RG5t++WdEhDc3UCb79++nT58+xMXFAVC/fn0SEhJo3LhxICsVETkj9QjOxsCBEBmZ/jmBr76C7duhatVTmn300UdERUURFxdHeHg4Q4YMISkpSSEgIvmSegRnIy4Ofv/dW46MhOhob3zAZ9euXfTs2ZPExEQAmjZtSkJCAnXr1g1CsSIiuaMeQW589ZU3JrDTd1vlr7+GH374MwScc7z11lvUqVOHxMRESpYsyejRo1m8eLFCQETyPfUIcmPWLPjxR285MtK7QqhYMQC2bt1K165dmeW7wXyLFi2Ii4vj0ksvDVKxIiJnRz2Cs9GtG6xfD8WKkZaWxsSJE6lbty6zZs2iXLlyxMfH89lnnykERKRAUY/gbFStCsWK8eOPP9KxY0cWLFgAQKtWrRg3bhxVMw0ai4gUBOoRnMny5fDccwCkpKby4osv0qBBAxYsWEDlypWZNm0aM2bMUAiISIGlHsGZ+CaFWwXEvPkmSRs2APDwww8zevRoKlasGMTiRETOnYIgJ3v2cPxf/2IY8HxYGCkbNnDRRRcxceJEWrZsGezqRETyhE4N5WBJYiKNgKFASloa3bp1Y/Xq1QoBESlU1CPIwuGHHmJAYiJjAAfUjIggfu5cbrjhhmCXJiKS5xQEJzkHqanM+/xzOr3zDpuBIkBvYFC/fpRQCIhIIaUgAHCOfTfdxJMLFjDZt6khkDBvHlffdBOE6QyaiBRe+g3nHB+88w5RvhAoBgwzY/lf/qIQEJGQENI9gt937ODR+vWZvns3ANcACWvXUqdOneAWJiISQCH5565zjqlTpxJVty7Td++mFPBK0aIsuv9+hYCIhBy/9gjM7DZgDN64a7xzbnim/cWBN4CrgT3A351zP/uzpi1bttC5c2c++eQTAG4F4syIPHHCn28rIpJv+a1HYGZFgFigJRAFPGRmmW/lFQPsdc5dDowGRvirnrS9e4kdPpy6UVF88sknlC9Xjil33slcINI5f72tiEi+588eQRNgo3NuE4CZJQKtgLUZ2rQCnvEtvwu8ambmXN7+Zk5NTeXWK65g/q5dANwLxO7fT5WPP/Ya9O2bl28nIlKg+DMIqgFbM6xvA5pm18Y5l2Jm+4GKwO6MjcysE9AJoEaNGmddSJEiRWhWpQprd+8mtnRp7itePH1n6dJw++1n/TVFRAqLAnHVkHMuDogDiI6O/q96CwOXLaP3kSOcd955eVqbiEhB588g+BW4KMN6dd+2rNpsM7OiQDm8QeM8FxERQUREhD++tIhIgebPy0eXAzXN7BIzKwY8CMzM1GYm0M63fD/wRV6PD4iISM781iPwnfPvAczFu3x0snNujZkNAZKcczOBBGCqmW0E/sALCxERCSC/jhE452YDszNtG5hh+Rjwv/6sQUREchaSnywWEZF0CgIRkRCnIBARCXEKAhGREGcF7WpNM9sF/PJfvrwSmT61HAJ0zKFBxxwazuWYL3bOnZ/VjgIXBOfCzJKcc9HBriOQdMyhQcccGvx1zDo1JCIS4hQEIiIhLtSCIC7YBQSBjjk06JhDg1+OOaTGCERE5HSh1iMQEZFMFAQiIiGuUAaBmd1mZuvNbKOZnXYfSjMrbmbv+PYvNbPIIJSZp3JxzE+Y2VozW2Vmn5vZxcGoMy+d6ZgztLvPzJyZFfhLDXNzzGb2gO97vcbM3gp0jXktFz/bNcxsvpmt8P18F+hbDprZZDPbaWars9lvZvaK799jlZk1Ouc3dc4VqgfelNc/AZcCxYCVQFSmNt2ACb7lB4F3gl13AI75JqCkb7lrKByzr10ZYCGwBIgOdt0B+D7XBFYAFXzrlYNddwCOOQ7o6luOAn4Odt3neMw3AI2A1dnsvx2YAxjQDFh6ru9ZGHsETYCNzrlNzrlkIBFolalNK+B13/K7QAszswDWmNfOeMzOufnOuSO+1SV4d4wryHLzfQYYCowAjgWyOD/JzTF3BGKdc3sBnHM7A1xjXsvNMTugrG+5HPBbAOvLc865hXj3Z8lOK+AN51kClDezC8/lPQtjEFQDtmZY3+bblmUb51wKsB+oGJDq/CM3x5xRDN5fFAXZGY/Z12W+yDk3K5CF+VFuvs+1gFpmttjMlpjZbQGrzj9yc8zPAG3MbBve/U8eDUxpQXO2/9/PqEDcvF7yjpm1AaKBG4Ndiz+ZWRgwCvhnkEsJtKJ4p4ea4/X6FppZfefcvmAW5WcPAa8550aa2TV4dz2s55xLC3ZhBUVh7BH8ClyUYb26b1uWbcysKF53ck9AqvOP3BwzZnYL8DRwt3PueIBq85czHXMZoB7wpZn9jHcudWYBHzDOzfd5GzDTOXfCObcZ2IAXDAVVbo45BpgG4Jz7GojAm5ytsMrV//ezURiDYDlQ08wuMbNieIPBMzO1mQm08y3fD3zhfKMwBdQZj9nMrgIm4oVAQT9vDGc4ZufcfudcJedcpHMuEm9c5G7nXFJwys0TufnZ/gCvN4CZVcI7VbQpgDXmtdwc8xagBYCZ1cELgl0BrTKwZgJtfVcPNQP2O+e2n8sXLHSnhpxzKWbWA5iLd8XBZOfcGjMbAiQ552YCCXjdx414gzIPBq/ic5fLY34RKA1M942Lb3HO3R20os9RLo+5UMnlMc8F/mpma4FU4CnnXIHt7ebymJ8EJpnZ43gDx/8syH/YmdnbeGFeyTfuMQgIB3DOTcAbB7kd2AgcAR455/cswP9eIiKSBwrjqSERETkLCgIRkRCnIBARCXEKAhGREKcgEBEJcQoCkVwys1Qz+y7DI9LMmpvZft/6D2Y2yNc24/Z1ZvZSsOsXyU6h+xyBiB8ddc5dmXGDbwrzRc65O82sFPCdmX3k231yewlghZnNcM4tDmzJImemHoFIHnHOHQa+AS7PtP0o8B3nODGYiL8oCERyr0SG00IzMu80s4p4cxqtybS9At58PwsDU6bI2dGpIZHcO+3UkM/1ZrYCSAOG+6ZAaO7bvhIvBF52zu0IWKUiZ0FBIHLuFjnn7sxuu5ldAiwxs2nOue8CXJvIGenUkIif+aaDHg78K9i1iGRFQSASGBOAG3xXGYnkK5p9VEQkxKlHICIS4hQEIiIhTkEgIhLiFAQiIiFOQSAiEuIUBCIiIU5BICIS4v4fhnO3Pzcjw+0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting of ROC-AUC Curve\n",
    "\n",
    "plt.plot(all_tpr_together_5_a, all_fpr_together_5_a, 'r', lw=2)\n",
    "plt.plot([0, 1], [0, 1], 'k-', lw=2)\n",
    "plt.xlabel('FPR')\n",
    "plt.ylabel('TPR')\n",
    "plt.title('AUC={}'.format(round(auc_score_5_a, 4)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Annotating threshold value in the RUC-AOC curve\n",
    "\n",
    "When we only plot the TPR and the FPR against each other we'll loose the threshold information. However, we can easily add them to the plot. In the below example with some randomly generated X and Y data, I annotated every 5th value but this should be enough the see the relationship (high confidence - bottom left, low confidence - top right)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAFzCAYAAAB2A95GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABxdUlEQVR4nO3dd1iUx/bA8e9QVBQVewEMGitNULDEhhWVBI09am7UGG8STXJjfrYkliTmqtEUW3oxN4m9IEnsvRfs2LuCvYPSFub3x8oGpEhblnI+z8MT9t155z27G+Ew78wcpbVGCCGEEELkLitLByCEEEIIURhJEiaEEEIIYQGShAkhhBBCWIAkYUIIIYQQFiBJmBBCCCGEBUgSJoQQQghhATaWDiCzypcvr11cXCwdhhBCCCHEU+3bt++W1rpCas/luyTMxcWFkJAQS4chhBBCCPFUSqmLaT0ntyOFEEIIISxAkjAhhBBCCAuQJEwIIYQQwgLy3Zyw1MTFxREWFkZ0dLSlQymQihUrhpOTE7a2tpYORQghhCgwCkQSFhYWRsmSJXFxcUEpZelwChStNbdv3yYsLIzq1atbOhwhhBCiwCgQtyOjo6MpV66cJGBmoJSiXLlyMsoohBBC5LACkYQBkoCZkby3QgghRM4rMEmYpVlbW+Pl5YW7uzsvvPAC9+7dMz139OhR2rRpQ506dahVqxaffPIJWmvT8ytXrsTHxwdXV1e8vb157733LPAKhBBCCJGbCmUSFnQgnGaTN1B99N80m7yBoAPh2e7Tzs6OgwcPEhoaStmyZZk9ezYAUVFRBAYGMnr0aE6ePMmhQ4fYsWMHX3/9NQChoaEMGzaM33//nWPHjhESEkLNmjWzHU9SBoMhR/sTQgghRPaZLQlTSv2slLqhlApN43mllJqhlDqjlDqslGpgrliSCjoQzpilRwi/F4UGwu9FMWbpkRxJxBI1bdqU8HBjf3PnzqVZs2Z06NABgOLFizNr1iwmT54MwGeffcYHH3xA3bp1AeOI2htvvJGiz8jISAYOHIiHhweenp4sWbIEAHt7e1ObxYsXM2DAAAAGDBjA66+/TuPGjRk5ciQuLi7JRudq1arF9evXuXnzJt27d8fX1xdfX1+2b9+eY++DEEIIkReZYzAmK8y5OnIOMAv4XxrPdwJqPf5qDHzz+L/Z8tGfRzl25UGazx+4dI/Y+IRkx6Li4hm5+DDz9lxK9RzXqqUY/4Jbhq4fHx/P+vXrefXVVwHjrciGDRsma/Pss88SGRnJgwcPCA0NzdDtx08++YTSpUtz5MgRAO7evfvUc8LCwtixYwfW1tbEx8ezbNkyBg4cyO7du3nmmWeoVKkSffv25d1336V58+ZcunQJf39/jh8/nqHXKoQQQuQ3iYMxUXHxwD+DMQBdvR1zNRazjYRprbcAd9Jp0gX4nzbaBTgopaqYK55ETyZgTzueUVFRUXh5eVG5cmWuX79O+/bts9Xfk9atW8fQoUNNj8uUKfPUc3r27Im1tTUAvXv3ZsGCBQDMnz+f3r17m/odNmwYXl5eBAYG8uDBAyIjI3M0diGEECKvmLr6JFFx8bSyOkRZjIM2UXHxTF19MtdjseQ+YY7A5SSPwx4fu/pkQ6XUEGAIQLVq1dLt9GkjVs0mbyD8XlTKYBzsWPDvpk+LOU2Jc8IePXqEv78/s2fP5u2338bV1ZUtW7Yka3vu3Dns7e0pVaoUbm5u7Nu3j/r162fpuklXLj65jUSJEiVM3zdt2pQzZ85w8+ZNgoKC+PDDDwFISEhg165dFCtWLEvXF0IIIfIT2/vn+cH2d9pb72eWoQvTDMZBiSup5Abmli8m5mutv9da+2itfSpUqJCtvkb418HO1jrZMTtba0b418lWv4mKFy/OjBkz+PzzzzEYDPTr149t27axbt06wDhi9vbbbzNy5EhjPCNG8N///pdTp04BxqTo22+/TdFv+/btTZP94Z/bkZUqVeL48eMkJCSwbNmyNONSSvHiiy8yfPhw6tWrR7ly5QDo0KEDM2fONLU7ePBg9t4AIYQQIi+KiYR1E1hTdCRNrY7x37iXmG7obnq6qoNdrodkySQsHHBO8tjp8TGz6urtyKRuHjg62KEwjoBN6uaRo/eBvb298fT0ZN68edjZ2bF8+XImTpxInTp18PDwwNfXl2HDhgHg6enJV199xUsvvUS9evVwd3fn3LlzKfr88MMPuXv3Lu7u7tSvX5+NGzcCMHnyZJ5//nmee+45qlRJ/25u7969+f333023IgFmzJhBSEgInp6euLq6ppoACiGEEPmW1nBoAczygW1fcs35eTolfMn38S8Q9/iGYE4OxmSGSrpfVY53rpQL8JfW2j2V5wKAYUBnjBPyZ2itGz2tTx8fHx0SEpLs2PHjx6lXr16OxCxSJ++xEEKIfOfKAVg5Ci7vhqoNoPNUcPJh4d7LjFxyGDAOxozwr2O2SflKqX1aa5/UnjPbnDCl1DzADyivlAoDxgO2AFrrb4EVGBOwM8AjYKC5YhFCCCFEwRd0IJypq08Sc+8a40os4YX49agS5aHLbKjfF6zy1iwsc66OfElrXUVrbau1dtJa/6S1/vZxAsbjVZFDtdbPaq09tNYhT+tTCCGEECI1QQfCGbv0AP4RS9lQ9D06GTbyS0IAf7X6C7z7mxKwoAPhjAv+ZwvTM/u30btdY6o4u5j28EzNkiVLUEqReDfu9u3btG7dGnt7e9MUo8yy5OpIIYQQQogcsWnFApaqH6llG87meE8+NrzMWe1IseDzBB2PMLXbdvoW0QbjtlQ6IZ47a7+hYu+JVHNyYt68kQQGBuLq6pqs74iICKZPn07jxv9sZ1qsWDE++eQTQkNDCQ1NdV/6p8pb43JCCCGEEJlx5zzM78dXcR9hi4FXY9/jlbhRnNXGOV7RhgSu3o82fSUmYACxV09h41AFW4fKXIs00KdPH5YvX57iEmPHjmXUqFHJtnMqUaIEzZs3z9YWTzISJoQQQoj8J/YhbP0CdswEKxu+senPl5HtiTVOPzdxdLDj77dbmB4n3S/UEHEbm1LGra+qOtjh5OjE7t27k52/f/9+Ll++TEBAAFOnTs3RlyAjYUIIIYTIP7SGI4thli9snQZuXeGtEKoEvI+1bfJRqdS2nsjMfqEJCQkMHz6czz//PMdfBshImNlMmDABe3t7Hjx4QMuWLWnXrp1F4/Hz82PatGn4+KS6SlYIIYTI+64eNm45cWkHVKkPPX6Gak0A6OptbDLx72PcioylvH0RPgxwTbH1ROLjqatPcq5kOWyi7pj2C520KgxHx3/aR0REEBoaip+fHwDXrl0jMDCQ4ODgHPl9KkmYmX388ceWDkEIIYTI3x7eho0TYd8csCsDL8x4vOIx+YhWV29HKpQsSr8fd/N1v4Y0ql421e66ejvS1dsRg6EltWt/TX2HWGJjY5k/fz5z5841tStdujS3bt0yPc7pAQ25HZmDPv30U2rXrk3z5s05edJYCHTAgAEsXrwYMCZkvr6+uLu7M2TIEBI3yt27dy+enp54eXkxYsQI3N2Ne9tGR0czcOBAPDw88Pb2Nu2SP2fOHLp160bHjh2pVauWqQQSwBtvvIGPjw9ubm6MHz8+N1++EEIIkbPiDbD7e5jpDft+hUb/hrf2Q8NXUiRgWWFjY8OsWbPw9/enXr169OrVCzc3N8aNG0dwcPBTz3dxcWH48OHMmTMHJycnjh07lrnrZzXwPGvlaLh2JGf7rOwBndLeOwRg3759zJ8/n4MHD2IwGGjQoAENGzZM1mbYsGGMGzcOgJdffpm//vqLF154gYEDB/LDDz/QtGlTRo8ebWo/e/ZslFIcOXKEEydO0KFDB1ONyYMHD3LgwAGKFi1KnTp1eOutt3B2dubTTz+lbNmyxMfH07ZtWw4fPoynp2fOvh9CCCGEuZ3bDKtGw41jUL0VdJoCFXO+ckvnzp3p3LlzsmNp3cXatGlTsscXLlzI1rVlJCyHbN26lRdffJHixYtTqlQpAgMDU7TZuHEjjRs3xsPDgw0bNnD06FHu3btHREQETZs2BaBv376m9tu2baN///4A1K1bl2eeecaUhLVt25bSpUtTrFgxXF1duXjxIgALFy6kQYMGeHt7c/To0Uxn5UIIIYRF3bsEC/8F/wuE2Ejo/Tv8a7lZEjBLK3gjYU8ZsbKU6Oho3nzzTUJCQnB2dmbChAlER0dnub+iRYuavre2tsZgMHD+/HmmTZvG3r17KVOmDAMGDMjWNYQQQohcE/sItn8F26cDClp/CM8NA1s7S0dmNjISlkNatmxJUFAQUVFRRERE8OeffyZ7PjEZKl++PJGRkaZ5Yg4ODpQsWdK0L8n8+fNN57Ro0YI//vgDgFOnTnHp0iXq1Em7yvuDBw8oUaIEpUuX5vr166xcuTJHX6MQQgiR47SGo8uMW05sngJ1A+CtEGg1okAnYFAQR8IspEGDBvTu3Zv69etTsWJFfH19kz3v4ODAa6+9hru7O5UrV072/E8//cRrr72GlZUVrVq1onTp0gC8+eabvPHGG3h4eGBjY8OcOXOSjYA9qX79+nh7e1O3bl2cnZ1p1qyZeV6sEEKIXJdYnPrKvSiqOtgxwr9Oiu0XcrOfHOnrWqhxy4mL26CSB3T7Hlyy/rsr6EA4E/82TsN58499qW5RkZeoxBV6+YWPj49OLJ6Z6Pjx49Srl3/vFUdGRmJvbw/A5MmTuXr1KtOnT7dwVMnl9/dYCCHys6AD4YxZeoSouHjTMTtba9P+VrndT7b7enQHNn4KIT9DsdLQZiw0HJCtFY85+dpyklJqn9Y61T0tZCQsD/j777+ZNGkSBoOBZ555hjlz5lg6JCGEEHnI1NUnkyUXAFFx8YwNCuXMjcgM9/Prjgs50k9W+1I6Hs/ry2h++VuKGiI5VKk7252HEH2nNKw9k6nrZzSeqatP5tnRMBkJExki77EQQlhO9dF/k9Zva2srleF+4hPS/p2fmX6y0lcjdYxx1r9Sz+oSuxJc+Tj+X5zUz2TqmlmJRwHnJwfk2HUyS0bChBBCiHysqoOdqeh0Uo4Odmwf3SbD/SQtXp2dfjLV1/0wWDMWji6F0s7Q4VeauHZhhcpc0pfVeKo65N3J/QVmdWR+G9HLT+S9FUIIyxrhX4ditsl/ZadVdPpp/WS0eHW2+4qLgs2fwUwfOLkCWo2GoXuMBbdzOAHLUDx5UIEYCStWrBi3b9+mXLlyKDN8sIWZ1prbt29TrFixpzcWQghhFl29HYmJi2fUUmNFGMcsrmpMWrw6u6sj0+zLqyocC4Y1Hxg3XnXtAh0mgkO1TF8jR+LJo/PBoIAkYU5OToSFhXHz5k1Lh1IgFStWDCcnJ0uHIYQQuSYjWy/kxS0jMiKxeLVZ+rpxHP73BpzfDBVd4ZU/oXrLHLlWluLJ4wpEEmZra0v16tUtHYYQQogC4MmtDsLvRTHm8QhU4i/4jLTJqWslthv/51HT46xeLzWrVq3inXfeIT4+nsGDByerYZzUkiVL6NGjB3v37sXHx4e1a9cyevRoYmNjKWJjzdRedWkTuwaKloROU8FnEFgXiDTDbArE6kghhBAip6Q1wVspKFXMFoAH0XGk9uszaZuMyGg/abXLyoT6pOLj46lduzZr167FyckJX19f5s2bh6ura7J2ERERBAQEEBsby6xZs/Dx8eHAgQNUqlCeqtfWETp3LP4/XiF87jvGckMlymU5poJGVkcKIYQQGXQllQQMjNV1Xnw86jRnx4WntsmIjPaTVru0Ys2oPXv2ULNmTWrUqAFAnz59WL58eYokbOzYsYwaNYqpU6eajnmXi4G/+sK1w7i5NSXKOoKY9pPTrewikpMkTAghhEgive0gJgS6AbD22PWntsmIjPaTVrvsbr8QHh6Os7Oz6bGTk5OplnGi/fv3c/nyZQICAoxJWORNWDIYjiyCklWh+08sOaFp0NAgCVgmFZgtKoQQQoickJGtDnJqO4SM9mOp7RcSEhIYPnw4n3/+OcRFw/3LsPBl4+rHliPgrRCOWtVj1OjRfPfdd2aNpSCSkTAhhBB5Xm4Wne7q7UjIxTv8vusSANZK0b2hY4o2kP3tEDLaj7m2X3B0dOTy5cumx2FhYTg6/tNnREQEoaGh+DVrDFF3ufYgjsALdgQvnodPmy6EhYXx4osv8r///Y9nn302W7EURjIxXwghRJ6W20Wnn2wTdW4f99Z/T9niNrwz9PUMrx7cs2cPQ4YMAYx7Lk6YMIEXX3wxU/Gam8FgoHbt2qxfvx5HR0d8fX2ZO3cubm6Pb4XePAmrRsPZDVC+Dn6/xTJt9o/4+Phw7949WrVqxfjx4+nWrZtlX0geJhPzhRBC5FtpFa/+6M+jma53+NGfR5/aV9I2OiGeO2u/oWLviVRzcmLevJEEBgamunpw+vTpNG7c2HTM3d2dkJAQbGxsuHr1KvXr1+eFF17Axibv/Oq1sbFh1qxZ+Pv7Ex8fz6BBg3Bzc2Pc+yPxKXKeQKsNYFsCOk4G38GwuL3p3FmzZnHmzBk+/vhjPv74YwDWrFlDxYoVLfVy8h0ZCRNCCJGnpVe82txiwo9zb9tcKvX+BAUMKX0YgDFjxiRr95///If27dszdepUpk2bho9P8oGP8+fP06RJE8LDw/NUEpZCQgIc/APWfwQPb0GDl6HNOLCvYOnI8i0ZCRNCCJFvpbVasWLJosx9rXEqZ6St7w+7uRERk25fSdsYIm5jU6qCKQ4nxwysHkxi9+7dDBo0iIsXL/Lbb7/l7QTs8h5YORKuHACnRtBvEVT1tnRUBVoe/r9BCCGEMK4MTG0e1/ud61GzYslM9fV+53pP7SutNiP86xARGp6sv8TVg3PmzEn1eo0bN+bo0aMcP36cV155hU6dOuW9WrwR12DdBDg0D+wrQ7cfwKOnWYpsi+RkiwohhBB5WldvRyZ186B4EeMWDY4OdlmalJ+0L0cHO1QafSVtY1OyHDZRd0xt0lw96OeHi4sLu3btIjAwkCenzdSrVw97e3tCQ0Oz9iaYgyEGtn0FMxtC6BJo/i68FQKevSQByyUyJ0wIIUS+MGrxYTafusmu99vm2jWfunrwCX5+fqY5YefPn8fZ2RkbGxsuXrxI06ZNOXz4MOXLl8+1+NN0arVx1eOdc1C7E/h/CuVkiwlzkDlhQgghRBakuXpw3Dh8fHwIDAxM89xt27YxefJkbG1tsbKy4uuvv7Z8AnbrDKweA6fXQLla0G8J1Gpn2ZgKMRkJE0IIkS9YYiSswIh+AFumwq5vwKYY+I2CRv8GmyKWjqzAk5EwIYQQojBKSIDD840T7yOvg1d/aDsOSlaydGQCScKEEEKIgil8H6wYCeEh4OgDfeaBU0NLRyWSkCRMCCGEKEgib8C6j+Dg71CiInT9Bjz7gJVsiJDXSBImhBBCFASGWNjzHWz+DOKi4Lm3oeUIKFbK0pGJNEgSJoQQwmyCDoQzdfVJrtyLoqqDHSP862Rpf6+gA+H8efgKj2LjaTZ5Q5b7KbDOrIOVo+H2aajZzljrsXwtS0clnkKSMCGEEGYRdCA82c7z4feiGLXkMNcfRNG2XsYnhq8/fp0v1p4mxpBg6mfM0iMAkojdOQerP4CTK6BsDXhpAdT2l81W8wnZokIIIYRZNJu8IdWajznF0cGO7aPbmK3/PC0mErZ+DjtngZUttBoBTd4Em6KWjkw8QbaoEEIIkeuupJOAzXwp44Wh35p3INP9F1haw5FFsHYcRFw1TrhvNwFKVbF0ZCILJAkTQghhFlUd7FIdCXN0sOOF+lUz3M/klSdS7aeqg1224st3rhyElaPg8i6o4gW9/gfOjSwdlcgGWa8qhBDCLEb416GYbfJfM3a21ozwr5PpfuxsrbPdT7718BYEvw3f+8HtMxA4E17bKAlYASAjYUIIIcyiq7cjO8/eYkFIGADWStG9oWOmJ9Mnts+JVZb5Snwc7P0RNk6CuIfGOV+tRoKdg6UjEzlERsKEEEKYRdCBcJYfumJ6HHk2hKmvdqKKswuTJ09O87wlS5aglCJxEdaePXuYMDCAh/OHU2rl+/xfnfsFPwE7uxG+bQ6rRoNjA3h9O3T8ryRgBYyMhAkhhDCLqatPEh1n3FZCJ8RzZ+03VOw9kWpOTsybN5LAwEBcXV2TnRMREcH06dNp3Lix6Zi7uzshISHY2Nhw9epV6tevzwsvvICNTQH8FXb3gnHLiRN/gcMz0Gcu1OksW04UUDISJoQQwiySrl6MvXoKG4cq2DpU5lqkgT59+rB8+fIU54wdO5ZRo0ZRrFgx07HixYubEq7o6GhUQUxIYh/Chk9hViM4uwHajIWhe6BugCRgBZgkYUIIIcwi6epFQ8RtbEpVMB13cnIiPDw8Wfv9+/dz+fJlAgICUvS1e/du3Nzc8PDw4Ntvvy04o2BaQ+gSmOULWz6Dei/AsBBo+X9gW+zp54t8TZIwIYQQZpGZ1ZEJCQkMHz6czz//PNW+GjduzNGjR9m7dy+TJk0iOjraLDHnqmtHYE4ALB4ExcvCwJXQ4ycoXcDnuwkTsyZhSqmOSqmTSqkzSqnRqTxfTSm1USl1QCl1WCnV2ZzxCCGEyD1dvR2Z8IJxzpdNyXLYRN1hUjcPuno7EhYWhqPjP8lGREQEoaGh+Pn54eLiwq5duwgMDOTJCin16tXD3t6e0NDQXH0tOerRHfhrOHzXEm4ch+e/hCGb4ZnnLB2ZyGVmS8KUUtbAbKAT4Aq8pJRyfaLZh8BCrbU30Af42lzxCCGEyH0BnsZNWScMCqSM4Tb1HWKJjY1l/vz5BAYGmtqVLl2aW7duceHCBS5cuECTJk0IDg7Gx8eH8+fPYzAYALh48SInTpzAxcXFEi8ne+INsOcHmOEN++aA72vw9n7wGQRW1k89XRQ85ryp3gg4o7U+B6CUmg90AY4laaOBUo+/Lw1cQQghRIFjbWPDrFmz8Pf3Jz4+nkGDBuHm5sa4cePw8fFJlpA9adu2bUyePBlbW1usrKz4+uuvKV++fC5GnwPObzXudn/jKFRvCR2nQKUnxyVEYWO2At5KqR5AR6314MePXwYaa62HJWlTBVgDlAFKAO201vvS61cKeAshRP4RER2Hx4Q1fBhQj8Etalg6nNx37xKsGQvHgqB0NfCfCPUCZcVjIZKXC3i/BMzRWn+ulGoK/KaUctdaJyRtpJQaAgwBqFatmgXCFEIIITIhLgq2T4dtXwIK/N6HZm+DbSGrdynSZc4kLBxwTvLY6fGxpF4FOgJorXcqpYoB5YEbSRtprb8HvgfjSJi5AhZCCCGyRWs4ttw4+nX/Eri9CO0/AQfnp58rCh1zJmF7gVpKqeoYk68+QN8n2lwC2gJzlFL1gGLATTPGJIQQQpjH9aPGeV8XtkJFN3jlL6jewtJRiTzMbEmY1tqglBoGrAasgZ+11keVUh8DIVrrYOA94Ael1LsYJ+kP0OaapCaEECJDgg6E51ix7L8PG9dbTfz7OL9sv1AwC28/ugObJsHen6BoSeg8DRoOBGtLz/gReZ3ZJuabi0zMF0II8wk6EM6YpUeIios3HbOztTbt75XZvkYvPWyqH5mdvvKkhHjjVhMbJkL0PeNWE60/MG68KsRjeXlivhBCiDxk6uqTyRIwgKi4eCYEH01x/GmmrDyRLAFL7Gvq6pP5Pwm7uANWjjTuev9MM+g0BSp7WDoqkc9IEiaEEMIkadHtpO5FxTFm6RGzXiNfuB8Ga8cZ6z2WcoIevxgn38uWEyILJAkTQghhUtXBjvBUkqTKpYoRNLRZpvrqOns71x6krPGYtLB3vhEXDTtnwtYvjLchW46E5v+BIiUsHZnIxyQJE0IIYTLCv06qc8JGd6pL5dLFMtXX6E51U+0rtQLeeZbWcOJvWP0+3LsI9V6ADhOhjIulIxMFgCRhQgghTBLnao1dHkpEtAHHbKyOTDwnp1Za5rqbJ41bTpzbCBXqwb+WQw0/S0clChBJwoQQQiTT1duR87ceMn39abaPbpPtvvJN0pUo6h5sngJ7vjfebuw4BXxfBWtbS0cmChhJwoQQQgiAhAQ4+Dus+wge3YaGr0CbsVAinxULF/mGJGFCCCHEpd3GLSeuHgTnJtB/CVT1snRUooCTJEwIIUTh9eAqrBsPhxdAyarQ7Ufw6CFbTohcIUmYEEKIwscQA7u+hs1TISEOWrwHzYdDUXtLRyYKEUnChBBCFB5aw6nVsHoM3DkHdQLAfyKUrWHpyEQhJEmYEEKIwuHWaVg1Bs6shXK1jPO+arazdFSiEJMkTAghRMEW/QC2fAa7vgHb4tDhU2g0BGyKWDoyUchJEiaEEKJgSkiAQ/Ng3QR4eAO8+0Pb8WBf0dKRCQFIEiaEEKIgCtsHK0dA+D5w8oW+88GxoaWjEiIZScKEEEIUHBHXYf1HcPAPsK8EXb8Fz95gZWXpyIRIQZIwIYQQ+Z8hFvZ8B5umgCEamr0DLUdA0ZKWjkyINEkSJoQQIn87vQ5WjYbbp6GWP/j/F8rXtHRUQjyVJGFCCCHyp9tnYfUHcGollH0W+i6E2v6WjkqIDJMkTAghRP4SEwlbp8HO2WBdBNp/DI3fkC0nRL4jSZgQQoj8QWs4vNBY6zHiKtTvC+3GQ8nKlo5MiCyRJEwIIUTed+UArBwFl3dDVW/o9Rs4+1o6KiGyRZIwIYQQeVfkTdjwMez/DUqUhy6zjSNgsuWEKAAkCRNCCJH3xMfBnh9g02SIewhNh0KrkVCstKUjEyLHSBImhBAibzm7AVaOhlsn4dm20HEyVKht6aiEyHEyniuEECKZoAPh/Lz9PADNJm8g6EB47lz4znmY3w9+exHiY+Cl+dB/iSRgosCSkTAhhBAmQQfCGbP0CFFx8QCE34tizNIjAHT1djTPRWMfwtYvYMdMsLKBtuOgyVCwLWae6wmRR0gSJoQQwmTq6pOmBCxRVFw8U1efzPkkTGsIXQJrx8GDcPDoBe0/glJVc/Y6QuRRkoQJIYQwCb8XlerxK2kcz7Krh41bTlzaAVXqQ4+foVqTnL2GEHmcJGFCCCGINSTw+dqTaT5f1cEuZy708DZsnAj75oBdGXhhOni/DFbWOdO/EPmIJGFCCFHInbsZyTvzD3Ik/D5Na5TlwOV7RMclmJ63s7VmhH+d7F0k3gAhP8PGTyEmAhr9G/xGGRMxIQopScKEEKKQ0lqzYO9lPvrzGEVtrfi2f0M6ulcm6EA4U1ef5Mq9KKo62DHCv0725oOd32K89XjjGFRvBZ2mQMV6OfdChMinJAkTQohC6N6jWEYvOcKqo9doVrMcX/TyolIp42rErt6OOTMJ/94lWPMhHFsODtWg9+9Q93lQKvt9C1EASBImhBCFzI6ztxi+4BC3H8YwplNdXmtRAyurHEyMYh/Bjhmw7UtAQesP4blhYJtD88qEKCAkCRNCiEIi1pDAF2tP8d2Ws1QvX4IfX2mGu2MOlgHS2jjqteZDuH8Z3LpBh0+gtFPOXUOIAkSSMCGEKASSTr5/qVE1xj5fj+JFcvBXwPWjxnlfF7ZCJQ948TtwaZZz/QtRAEkSJoQQBVhak+9zzKM7sGkS7P3RWFw74AtoOEC2nBAiAyQJE0KIAureo1jGLD3CylDj5PvPe3pRuXQOlQJKiDfu9bVhIkTfA59XofX7ULxszvQvRCEgSZgQQhRAZp18f2G78dbj9SPg0gI6TobK7jnTtxCFiCRhQghRgCSbfF8uhyff3w8z1nkMXQKlnaHnr+DaRbacECKLJAkTQogCwmyT7+OiYcdM2PYF6ARoNRqavQNFime/byEKMUnChBAinzPb5Hut4cRfsPoDuHfROOrVYaJx41UhRLZJEiaEEPlY0sn3zz1r3Pk+Rybf3zgBq0bBuU1QoR78KxhqtMp+v0IIE0nChBAinzLL5Puoe7BpMuz5HoraQ6ep4DMIrOXXhRA5Tf5VCSFEPpC0qHaV0sWoV6UUG07eyLnJ9wnxcOA3WP+xce+vhgOgzVgoUS5H4hdCpCRJmBBC5HFBB8IZs/QIUXHxAFy5H82V+9E0rVGWnwb4Zn/y/aVdsHIkXD0E1ZpCpylQpX4ORC6ESI8kYUIIkcdNXX3SlIAldelOVPYSsAdXYO14OLIQSlaF7j+Be3fZckKIXCJJmBBC5HFX7kVl6vhTGWJg5yzY8jkkGKDF/0GL4VCkRDaiFEJkliRhQgiRx1V1sCM8lYSrqoNd5jrSGk6tglVj4O55qPu8ccuJstVzKFIhRGZYWToAIYQQ6RvhX4eiNsl/XNvZWjPCv07GO7l5Cn7vDvP6gHUReHkZ9PlDEjAhLEhGwoQQIo/r6u3I1ftRTFl1EgBHBztG+Nehq7fj00+Ovg+bP4Pd34JtcfCfBI1eA2tbM0cthHgasyZhSqmOwHTAGvhRaz05lTa9gAmABg5prfuaMyYhhMiP2tStxJRVJ/m6XwM6e1R5+gkJCXBoLqybAA9vQYOXoc04sK9g9liFEBljttuRSilrYDbQCXAFXlJKuT7RphYwBmimtXYD/mOueIQQoqBYtWoVderUoWbNmkyenOJvWwgLgZ/asWTSENTIs4Q0ng2BM00J2KVLl7C3t2fatGm5HLkQIilzzglrBJzRWp/TWscC84EuT7R5DZittb4LoLW+YcZ4hBAi34uPj2fo0KGsXLmSY8eOMW/ePI4dO2Z8MuI6LHsDfmxLxPXLTD/7DI0bN4YKyeeODR8+nE6dOlkgeiFEUuZMwhyBy0kehz0+llRtoLZSartSatfj25cpKKWGKKVClFIhN2/eNFO4QgiR950KPUDNmjWpUaMGRYoUoU+fPixfthS2z4CZDSF0MTR/l7HX2zPq42kUK5a8jmRQUBDVq1fHzc3NQq9ACJHI0qsjbYBagB/wEvCDUsrhyUZa6++11j5aa58KFWQ+gxCi8Lp9/RrOzs6mx0429whfPQPWjgWXZvDmLvaXDeTylesEBAQkOzcyMpIpU6Ywfvz43A5bCJEKc07MDweckzx2enwsqTBgt9Y6DjivlDqFMSnba8a4hBAi/7t91rjf144/gZLQbzHUak9CQgLDu7dhzpw5KU6ZMGEC7777Lvb29rkerhAiJXMmYXuBWkqp6hiTrz7AkysfgzCOgP2ilCqP8fbkOTPGJIQQ+VqFcqXZd3AjzA4Cm2KEVWqHo09jqNUegIiICEJDQ/Hz8wPg2rVrBAYGEhwczO7du1m8eDEjR47k3r17WFlZUaxYMYYNG2a5FyREIWa2JExrbVBKDQNWY9yi4met9VGl1MdAiNY6+PFzHZRSx4B4YITW+ra5YhJCiHxLJ/Ci1VZeC1/EzLPnOV+pL449JjG/zfPMfeufFZKlS5fm1q1bpsd+fn5MmzYNHx8ftm7dajo+YcIE7O3tJQETwoLMOidMa71Ca11ba/2s1vrTx8fGPU7A0EbDtdauWmsPrfV8c8YjhBD5Uvh+qgW9yJdFvsFgX5VZ07/E//M91GvkR69evXBzc2PcuHEEBwdbOlIhRCYorbWlY8gUHx8fHRISYukwhBDC/CJvwvqP4MDvGOzKM/pBN9r0fofOnhnYKV8IkScopfZprX1Se87SqyOFEEI8KT4Ods6GmQ3g0Hx4bhjnX9rC4vhWoOTHthAFhdSOFEKIvOTMeuOqx1snoWY76DgZytci4VqEpSMTQuQwScKEECIvuHMeVn8AJ/+GsjXgpQVQ2x+UsnRkQggzkSRMCCEsKSYStn0BO2aBlQ20mwBN3gSbopaOTAhhZpKECSGEJWgNRxbD2nEQcQU8+xgTsFJVLB2ZECKXSBImhBC57eohWDkKLu2EKl7Q61dwbmTpqIQQuUySMCGEyC0Pb8GGT2Dfr1C8HATOBK/+YCUrHoUojCQJE0IIc4s3QMhPsPFTiH1onPPVaiTYOVg6MiGEBUkSJoQQ5nRuE6wcDTePQ43Wxi0nKta1dFRCiDxAkjAhhDCHuxdhzQdw/E9weAb6zIU6nWXLCSGEiUxEEEKInBT7CDb+F2Y3wnBqHd/Z9KPOtY9oFmRH0MErWe52w4nrALz5x36aTd5A0IHwnIpYCGEhMhImhBA5QWs4ugzWjIUHYVx2DODlSwFciHMAIPxeFGOWHgGgq3fmaj8GHQjnq3WnTY+z05cQIu+QAt5CCJFd10KNW05c3AaVPaDTZzSbF034vagUTa0UlCleJFPd330US0IqP6odHezYPrpNVqMWQuSC9Ap4y0iYEEJk1aM7xhWPIT9DMQd4/kto8ApYWXPl3t+pnpKgoZNH5Uxd5vddl1I9fiWVJE8IkX9IEiaEEJkVb4B9vxgTsOgH4Psa+I2G4mVNTao62KU6EuboYMfErh6ZutzGEzdT7auqg13mYxdC5BmZnpivlLJSSvUzRzBCCJHnXdgG37eCFf8Hldzh9a3Q+bNkCRjACP862NlaJztmZ2vNCP86mb5kTvYlhMg70hwJU0qVAoYCjkAwsBYYBrwHHAL+yI0AhRAiT7h3GdaONU6+L+0Mvf4H9QLT3HIiccL8+8uO8Cg2HkcHO0b418nSRPrEc6auPsmVe1FUzUZfQoi8I82J+Uqp5cBdYCfQFqgIKOAdrfXB3ArwSTIxXwiRq+KiYPsM2PYloKH5u/Dc21CkeIZOH7HoENvP3GLHmLbmjVMIkSdldWJ+Da21x+MOfgSuAtW01tFmiFEIIfIWreF4MKz+EO5fAteu0OETcKhm6ciEEAVEeklYXOI3Wut4pVSYJGBCiELh+jFYNQrOb4GKbvDKX1C9haWjEkIUMOklYfWVUg8w3oIEsEvyWGutS5k9OiGEyE1Rd2HjJNj7IxQtCZ2nQcOBYC0LyYUQOS/Nnyxaa+u0nhNCiAIlIR72/w/WfwzR94yJV5sPU6x4FEKInJTe6shiwOtATeAw8LPW2pBbgQkhRK64tAtWjIBrh+GZZtBpinHXeyGEMLP0xth/xTgvbCvQGXAD3smNoIQQwuweXIG14+DIIijlCD1+BrduaW458aSgA+FP3TIi6EA4fx+5yqPYeJpN3iDbSgghkkkvCXNNsjryJ2BP7oQkhBBmFBcNO2fB1i8gwQAtR0Lz/0CREhnuIuhAOGOWHiEqLh5IvaB2RtoIIQq3jK6ONKgM/nUohBB5ktZwciWsHgN3L0C9F6DDRCjjkumupq4+aUquEkXFxTNqyWEW7bsMQMiFu8QYElK0mbr6pCRhQggg/STM6/FqSDCuiJTVkUKI/OnmSVg1Gs5ugAp14eUgeLZ1lrtLq3B2jCGBmLgE0/eZOVcIUfikl4Qd0lp751okQgiR06Lvw6YpsOc7sC0BHaeA76tgbZutbtMrzr34jecAaDZ5gxTdFkKkK70C3qnXMxJCiLwuIQH2/wYzG8Kur8G7P7y9H5q8nu0EDDJWUFuKbgshnia9kbCKSqnhaT2ptf7CDPEIIUT2XN4DK0fClQPg3Bj6LYaqXjl6icQ5Xf+36BCGBJ1qcW4pui2EeJr0kjBrwJ5/dswXQoi8K+IarJsAh+ZBySrQ7Ufw6JHhLScyq6u3IzM3nKZulVLM7tsgzTaSdAkh0pJeEnZVa/1xrkUihBBZYYiBXd/AlqkQHwvNh0OL96CovaUjE0KIdKWXhMkImBAibzu12rjq8c45qNMZ/D+FsjUsHZUQQmRIeklY21yLQgghMuPWGeN+X6fXQLla0H8J1Gxn6aiEECJT0ivgfSc3AxFCiKeKfmC87bjrG7C1gw6fQqMhYFPE0pEJIUSmpTcSJoQQeUNCAhyeb5x4H3nduOVE2/FgX9HSkQkhRJZJEiaEyNvC98GKkRAeAk6+8NI8cGyYollGCmpnVEaLc1+8/YizNx9y8JIU5xZCZJ4kYUKIvCnyBqz7CA7+DvaVoOu34NkbrFLuMZ2TxbIzU5zbkKCzfT0hROElSZgQIm8xxBrLDG3+DOKioNk70HIEFC2Z5ilpFdT+IOgIh8PuZ+ryC0IuPbWvtNpIcW4hRGZIEiaEyDtOrzNuOXH7NNTqAP6ToHzNp56WVlHshzHxLAq5nKkQHsbEp3k8sa+02khxbiFEZkgSJoSwvNtnYfUHcGollH0W+i6E2v4ZPj29gtrbR7fJVChpFd5O2pcU5xZC5IT0CngLIYR5xUQaVzx+3QQubIX2H8ObuzKVgEHOFsuW4txCiNwiI2FCiNynNRxZBGvHQcRVqN8X2o2HkpWz1F3iPKyxQaFExBhSLaid2b7SWx0pxbmFEDlBkjAhRO66cgBWjoLLu6GqN/T6DZx90z0lJ7efyIiMFN6W4txCiOyS25FCiNzx8BYEvw3ftzbWeuwyGwZvyFACNmbpEcLvRaH5ZzuIoAPhKdpExBggjTbmtmrVKurUqUPNmjWZPHlymu2WLFmCUoqQkBAALly4gJ2dHV5eXnh5efH666/nVshCCAuTkTAhhHnFx8HeH2HjJIh7CE2HQquRUKx0hk5Pa/uJ/yw4yMglhwGINSSkOC/plhH37t3jzJkznD59mitXrtC9e3dcXFyy/dISxcfHM3ToUNauXYuTkxO+vr4EBgbi6uqarF1ERATTp0+ncePGyY4/++yzHDx4MMfiEULkD5KECSHM5+xG45YTN0/As22h42SoUDtTXaS37cOgZtUB+Hbz2VSfP/THf7Ed64/BYEge1tmzfP3115mKIz179uyhZs2a1KhRA4A+ffqwfPnyFEnY2LFjGTVqFFOnTs2xawsh8i+5HSmEyHl3L8D8fvBbVzBEw0vzof+STCdgkPa2D44OdozuVJfRnerimEYbfftCsgTM3t6e5557ji5duqC1znQsaQkPD8fZ2dn02MnJifDw5LdC9+/fz+XLlwkICEhx/vnz5/H29qZVq1Zs3bo1x+ISQuRtkoQJIXJO7EPYMBFmNTKOgrUdB2/uhjqdQKksdZmdLSO+X/AX9evXB8DKyorIyEh27NhBx44dqVevHuPGjSM0NDRHE7LUJCQkMHz4cD7//PMUz1WpUoVLly5x4MABvvjiC/r27cuDBw/MGo8QIm8waxKmlOqolDqplDqjlBqdTrvuSimtlPIxZzxCCDPRGo4shlm+sGUquHaBt0KgxXtgWyxbXXf1dmRSNw+K2hh/XDk62DGpm0eKLSMmdfPA0cEOlaRNvxZ12b59O35+fiQkJFCuXDn69+9PuXLlOHnyJJ988gkeHh64uroyfvx4jh49mqUYHR0duXz5n535w8LCcHT8J76IiAhCQ0Px8/PDxcWFXbt2ERgYSEhICEWLFqVcuXIANGzYkGeffZZTp05lKQ4hRD6jtTbLF2ANnAVqAEWAQ4BrKu1KAluAXYDP0/pt2LChFkLkIVcOaf1TR63Hl9L62xZaX9xplsv0/WGn7v719iydGxkZqVu1aqUB/eeff+q4uDi9Zs0a/dprr+ly5cppwPTl6uqqx48fr48ePZrh/uPi4nT16tX1uXPndExMjPb09NShoaFptm/VqpXeu3ev1lrrGzduaIPBoLXW+uzZs7pq1ar69u3bWXqdQoi8BwjRaeQ05hwJawSc0Vqf01rHAvOBLqm0+wSYAkSbMRYhRE57eBv+ehe+bwW3TsIL0+G1jVCtiaUjS6FEiRKsXbuWAwcOEBAQgI2NDe3bt+f777/n6tWrrF69msGDB1O2bFmOHTvGRx99hJubG25ubnz00UccO3Ys3f5tbGyYNWsW/v7+1KtXj169euHm5sa4ceMIDg5O99wtW7bg6emJl5cXPXr04Ntvv6Vs2bI5+fKFEHmU0maaC6GU6gF01FoPfvz4ZaCx1npYkjYNgA+01t2VUpuA/9Nah6TS1xBgCEC1atUaXrx40SwxCyEyIN4AIT/Dxk8hJgIaDQG/UWBXxqyX7ffjLmLiElj8xnNmu0ZcXBwbNmxg0aJFLF26lLt375qec3Nzo1evXvTs2ZN69eqZLQYhRMGilNqntU51upXFJuYrpayAL4D3ntZWa/291tpHa+1ToUIF8wcnhEjd+S3wXQtYOQKq1Ic3tkOnyWZPwHKLra0t/v7+/Pjjj1y/fp1Vq1YxaNAgypQpw9GjRxk/fjyurq54eHjwySefcOLECUuHLITIx8yZhIUDzkkeOz0+lqgk4A5sUkpdAJoAwTI5X4g86N4lWPgK/PoCxEZC79/hX8uhYsEdEUpMyH766SeuXbvGypUrGThwIA4ODoSGhjJu3Djq1auHp6cnEydO5OTJk5YOWQiRz5jzdqQNcApoizH52gv01VqnuvwovduRSfn4+OjEch9CCDOLi4Lt02Hbl4AyrnZ8bhjYpr4vlznlxu3IjIiNjWX9+vUsXLiQZcuWcf/+fdNznp6epluWtWtnfk80IUTBk97tSLMlYY8v3Bn4CuNKyZ+11p8qpT7GuFIg+Im2m5AkTIi8QWs4thzWjIX7l8C9O7T/GEo75filMlKcO+hAOKOWHCbGkIBjLhTwzqjY2FjWrVvHwoULCQoKSpaQ1a9f35SQ1apVy4JRCiEsyWJJmDlIEiaEmV0/CitHwYWtUMkDOk0Bl2ZmuVRi4e2ktSHtbK2T7QOWkTZ5QUxMTLKELOmGq15eXqaErGbNmhaMUgiR2yQJE0I83aM7sGkS7P0JipWCNmOh4QCwsn7qqVnVbPIGwlOpDWlrrXCtUgqAY1cfEBef8ueUo4Md20e3MVts2RETE8PatWtZuHAhy5cvT5aQeXt707NnT0nIhCgkJAkTQqQtIR72zTGWG4q+Bz6vQuv3obj596qqPvpv0voJ5FfHuBJ608mbqT6vgPOTU9ZhzGtiYmJYs2aNKSGLiIgwPdegQQNTQvbss89aMEohhLlIEiaESN3FHbByJFw7Ai4toONkqOyea5dPayQs6ShXRtrkF9HR0aaELDg4OFlC1rBhQ1NCVqNGDQtGKYTISXlynzAhhAXdD4PFg+CXThB1D3r+Cq/8masJGGSvOHfSNvlFsWLFCAwM5Pfff+fGjRsEBQXRt29f7O3t2bdvH6NHj+bZZ5/Fx8eHzz77jPPnz1s6ZCGEGclImBCFSVw07JwJW78AnQDN/gPN3oEixS0WUtCBcN5dcBANaa58zMgKyvwsKiqK1atXs3DhQv78808iIyNNz/n4+Jgm9bu4uFguSCFElshImBCFndZw/C8eftkQNkxkRbQ73aynE+TwL4smYBnV1duR7aPbcH5yANtHt6HY9SPUqVOHmjVrMnny5DTPW7JkCUopkv7hdvjwYZo2bYqbmxseHh5ER1u+bK2dnR1du3Zl7ty53Lhxg6VLl9KnTx9KlChBSEgII0eOpHr16jRq1Ihp06YhpduEKBhkJEyIgu7GCVg1Cs5t4pR2Znzcv9iZ4AZAURsrhrWuiV+dihYJbdPJG8zaeIYYQ4Lp2NO2n4iPj6d27dqsXbsWJycnfH19mTdvHq6ursnaRUREEBAQQGxsLLNmzcLHxweDwUCDBg347bffqF+/Prdv38bBwQFra/OtAM2OqKgoVq5cyaJFi/jzzz95+PCh6blGjRrRq1cvevTowTPPPGPBKIUQ6ZGJ+UIURlH3YPMU2P0dFC3JF4aezI5sSTx5M+FIKr1J9zt37mTChAmsXr0agEmTJgEwZsyYZO3+85//0L59e6ZOncq0adPw8fFhxYoVzJ07l99//928L8AMHj16lCwhe/Tokem5xo0bmxKyatWqWTBKIcST0kvCbHI7GCGEmSXEw4HfYf3H8Og2+AyE1h8y85NdaW4H8eO/LFOydfD/Uv+D6koqqyEThYeH4+z8T1laJycndu/enazN/v37uXz5MgEBAUydOtV0/NSpUyil8Pf35+bNm/Tp04eRI0dm81XkjuLFi9O9e3e6d+/Oo0ePWLFiBYsWLeKvv/5i9+7d7N69m/fee48mTZqYErKk75MQIu+RJEyIguTSbuOWE1cPQrXnjLvdV/EEoKqDXZpbPbRzrZTLgf5z7dRiquqQ9dqUCQkJDB8+nDlz5qR4zmAwsG3bNvbu3Uvx4sVp27YtDRs2pG3btlm+niUUL16cHj160KNHDx4+fJgsIdu1axe7du1i+PDhNG3a1JSQOTnlfMkpIUT2yMR8IQqCB1dh6RD4uQNE3oDuP8HAFaYEDPLmVg9ZicnR0ZHLly+bHoeFheHo+M/8sYiICEJDQ/Hz88PFxYVdu3YRGBhISEgITk5OtGzZkvLly1O8eHE6d+7M/v37c/6F5aISJUrQs2dPFi5cyM2bN1mwYAE9evTAzs6OnTt38u677+Ls7EyzZs346quvCAsLs3TIQojHJAkTIj8zxBi3m5jZkPjQZcyx7onrzYk0+7ssQQevJGva1duR7g3/SVaslaJ7Q0eLbvXQ1duRSd08cHSwQ2EcGXtaTUhfX19Onz7N+fPniY2NZf78+QQGBpqeL126NLdu3eLChQtcuHCBJk2aEBwcjI+PD/7+/hw5coRHjx5hMBjYvHlzign9+VmJEiXo1asXixYt4saNG8yfP5/u3btTrFgxduzYYUrImjdvzvTp0wkPD7d0yEIUanI7Uoj8SGs4tQpWjYG757lSpR2vhAVyOq48AI/uRTFm6RGAZIWwl+z755duvNYs2ReOzzNlLZ6IZeb6NjY2zJo1C39/f+Lj4xk0aBBubm6MGzcOHx+fZAnZk8qUKcPw4cPx9fVFKUXnzp0JCMj7pY+ywt7ent69e9O7d28iIyP566+/WLRoEStWrGD79u1s376d//znPzRv3pyePXvSo0cPqlataumwhShUZHWkEPnNrdOwajScWQfl60CnyTRbRKpzq4oXseZ5zyoA/HX4Ko9i41O0yY/lf0TWRURE8Pfff7Nw4UJWrFhBTEwMAEopmjVrRq9evejevbskZELkENmiQoiCIPrB4y0nvgXbEtB6DPgOBmvbdAthVyldDICr91PflDS/FMIWOS8iIoK//vqLhQsXsnLlymQJWfPmzU0JWZUqVSwcqRD5lyRhQuRnCQlwaC6s+wge3oQGL0ObcWBfwdSksBXCFjnvwYMHyRKy2NhYwJiQtWjRwpSQVa5c2cKRCpG/SNkiIfKrsBD4qR0sHwplXGDIRgicmSwBA+Mqw2K2yf85F+RC2CLnlSpVir59+xIUFMTNmzf5/fffCQwMxNbWli1btjBs2DCqVq2Kn58fX3/9NdeuXbN0yELkezISJkReFHEd1k0wjoDZV4YOn4BHT1AqzVPm7r7I+8tCgcJbCFvkvPv37/Pnn3+ycOFCVq9ebRohs7KyomXLlvTq1Ytu3bpRqZJl9poTIq+T25FC5BeGWOOcr82fQXwMNB3KX6X6MmlD2FMTp4wkYUJkx/379wkODjYlZHFxcYAxIWvVqpUpIatY0TK1SIXIi+R2pBD5wem18E1TWDsWXJrBm7sIKvcaI/48R/i9KDTGFZCjlx5mcchlYg0Jpq/FIZf5+M9jpq7O7N9G73aNqeLswuTJk9O85JIlS1BKkfQPm8OHD9O0aVPc3Nzw8PAgOjr1Cf2i8CldujQvv/wyf/75Jzdu3ODXX38lICAAa2trNm7cyBtvvEGVKlVo27Yt3333HTdv3rR0yELkaTISJoSl3T5r3O/r9GooVxM6ToZa7YG0J9OnRyfEc+WHf1Ox90SqOTnxaNFI5s2bl2JT0oiICAICAoiNjWXWrFn4+PhgMBho0KABv/32G/Xr1+f27ds4ODhgbZ33i34Ly7l37x7Lly9n4cKFrF27NtkIWevWrenVqxcvvvgiFSpUeEpPQhQ8MhImRF4UEwFrx8PsxnBxB3SYCG/sNCVgkH4h6xH+dUxfScVePYWNQxVsHSpzLdJAnz59WL58eYrzx44dy6hRoyhWrJjp2Jo1a/D09KR+/foAlCtXThIw8VQODg688sor/P3331y/fp1ffvmFTp06YWVlxfr16/n3v/9NlSpVaN++PT/88AO3bt2ydMhC5AmShAmR2xIS4NB8mOkD278Cz17w1j547i2wKZKsaVqFrB0d7BjauqbpyzFJO0PEbWxKVTCd7+TklKI8zf79+7l8+XKK3eJPnTqFUgp/f38aNGjAZ599lgMvWBQmZcqUYcCAAaxYsYLr16/z888/07FjR5RSrFu3jiFDhlC5cmU6dOjAjz/+yO3bty0dshAWI0mYELkpfD/87A/L/g2lHWHweuj6NZRMfWVZRreVyMz2EwkJCQwfPpzPP/88xXMGg4Ft27bxxx9/sG3bNpYtW8b69esz+yqFAKBs2bIMHDiQlStXcv36dX766Sf8/f1RSrF27Vpee+01KlWqhL+/Pz/99JMkZKLQkSRMiNwQeROWD4Mf2sDdC9Dla3h1HTilOk3AJLHottXjnSmsFKkW3U5aCNumZDlsou6YCmGHhYXh6PhP+4iICEJDQ/Hz88PFxYVdu3YRGBhISEgITk5OtGzZkvLly1O8eHE6d+7M/v37c/rdEIVQ2bJlGTRoEKtWreLatWv8+OOP+Pv7A8bb4IMHD6Zy5cp07NiRn3/+mTt37lg4YiHMTybmC2FO8XGw53vYNBnioqDJ69ByJBQrlaHTgw6EM2bpEaLi/qn5aGdrbUqwUmMwGKhduzbr16/H0dERX19f5s6di5ubW6rt/fz8mDZtGj4+Pty9e5e2bduybds2ihQpQseOHXn33XcLbJFrYXm3bt0iKCiIRYsWsX79euLjjf+v29jY0K5dO3r16kXXrl0pU6aMhSMVImtkYr4QlnB2A3zTDFa/D86N4M2dxsn3GUzAAKauPpksAQOIiotn6uqTaZ5jY2PDrFmz8Pf3p169evTq1Qs3NzfGjRtHcHBwutcrU6YMw4cPx9fXFy8vLxo0aCAJmDCr8uXLM3jwYFavXs21a9f4/vvvad++PVprVq1axaBBg6hYsSKdO3dmzpw53L1719IhC5FjZCRMiJx25zys+RBO/AVla4D/JKjtn+5u92lJqzC3FN0WBd3NmzdZtmwZixYtYsOGDSQkJABga2tL+/bt6dWrF126dMHBwcGygQrxFLJjvhC5IfYhbP0CdswEKxtoNQKavAk2RbPcpRTdFgJu3LhhSsg2btyYLCHr0KEDvXr1IjAwUBIykSdJEiaEOWkNoUtgzViIuAKefaDdBChVJdtdZ2VOmBAF2Y0bN1i6dCmLFi1i06ZNyRIyf39/evbsSZcuXShdurSFIxXCSJIwIczl6iFYOQou7YQqXtB5qnH+Vw4KOhDOJ38d4/bDWCrYF+WDgHqSgAkBXL9+nWXLlrFw4UI2b95sSsiKFCliSsgCAwMlIRMWJUmYEDnt4W3Y8AnsmwPFy0G78eDVH6zMs9Zlw4nrDJoTwvKhzajv7GCWawiRn12/fp2lS5eaErLE322Jq3wTE7JSpTK+MEaInCCrI4XIKfEG2P0dzPSGA78Z53y9tQ8a/MtsCZgQ4ukqVarEG2+8wcaNG7ly5QqzZ8+mVatWxMXFERwczMsvv0zFihXp2rUrf/zxBw8ePLB0yEJIEiZEhp3bDN+1gJUjoWoDeH07dPwv2DlYOjIhRBKVK1fmzTffZNOmTYSHhzNr1ixatmxJbGwsy5cvp3///lSsWJEXX3yRuXPnEhERYemQRSElSZgQT3P3Iix4Gf4XaFwB2WcuvLwMKta1dGRCiKeoUqUKQ4cOZfPmzYSFhTFz5kxatGhBbGwsQUFB9OvXjwoVKtCtWzfmzZsnCZnIVZKECZGW2Eew8b8wuxGcWQdtxsLQPVA3IEt7fgkhLKtq1aoMGzaMLVu2EBYWxowZM2jevDmxsbEsW7aMvn37UrFiRbp37878+fOJjIy0dMiigJMkTIgnaQ1Hl8EsX9g8Beo+D8NCoOX/gW0xS0cnhMgBVatW5a233mLr1q1cvnyZ6dOn06xZM6Kjo1m6dCkvvfQSFSpUoEePHixYsEASMmEWsjpSiKSuhRq3nLi4DSp7QKfP4JnnstRV0IFwpq4+yZV7UVR1sGOEf50sbS0hW1QIkXvCwsJYsmQJCxcuZMeOHabjdnZ2BAQE0LNnTwICAihRooQFoxT5iWxRIcTTPLoDGz+FkJ+hmAO0HQsNXgEr6yx1l1ObrMpmrUJYzuXLl00J2c6dO03HExOyl19+mcDAQAtGKPIDScKESEtCPOz7BTZMhOgH4DsYWo8BuzLZ6jatckPl7Ysw4yXvDPfz9rwD3IqMTXFcyhYJkbsuX77M4sWLWbhwIbt27TId37JlCy1atLBgZCKvkyRMiNRc2Ga89Xg9FKq3hI5ToJJrjnSdVuHtnCIFvIWwnEuXLrF48WIuXLjAp59+SsmSJS0dksjD0kvCbHI7GCEs7n6Ysc7j0aVQuhr0+h/UC8yxFY8PouMoZmud7BZiovL2RZjVt0GG+xo2d3+qI2FVHeyyFaMQIuuqVavG8OHDLR2GKAAkCROFR1wU7JgJW78wPvZ7H5q9DbY5l9Dsu3iHd+YfJCouHhsrhSHhn/EwO1trPgxwpUmNchnu78MA11TnhI3wr5NjMQshhLAM2aJCFHxaw7Fg435fGz+FOh1h2F7wG5VjCZghPoEv156i57fGybtL3mjKtJ71cXSwQ2Gcw5WVyfRdvR2Z1M0j2/0IYWn37t3j66+/Nj3etGkTzz//fKptBw8ezLFjx8wSR0hICG+//Xa6bZ6MNatmzJhBvXr16NevX7b7Suqrr77i0aNH2TrH3t4+y9c/ePAgK1asyPL5OcXFxYVbt26lOP7f//7X9P2FCxdwd3fP8jViYmJo164dXl5eLFiwIMv9pEXmhImC7cZx47yv85uhoht0mgLV059Em5GtJZK2qViqKHa21ly4/YgXvR35qIsbpYrZmvNVCZHvXLhwgeeff57Q0FDAmIRNmzaNv/76y8KRpfRkrGnRWqO1xiqNurF169Zl3bp1ODk5Zei6BoMBG5un36BycXEhJCSE8uXLZ6jf1M6xt7fP8t5nc+bMISQkhFmzZmX4nIy+tsxI631I+toy+lmmZdeuXXz44YesW7cuw+c8+VqlgLcofKLuGpOvb5rB1UPQeRr8e0uGErAxS48Qfi8KDYTfi2LM0iMEHQhPs831BzFcuP2Il5tU48veXmkmYKtWraJOnTrUrFmTyZMnpxnDkiVLUEqR+MfGH3/8gZeXl+nLysqKgwcPZvYdEcKiRo8ezdmzZ/Hy8mLEiBEAREZG0qNHD+rWrUu/fv1IHBTw8/MjJCSE+Ph4BgwYgLu7Ox4eHnz55Zcp+h0wYACvv/46Pj4+1K5d25TURUdHM3DgQDw8PPD29mbjxo1A8hG4CRMmMGjQIPz8/KhRowYzZsxIM9ZEFy5coE6dOvzrX//C3d2dy5cvM3XqVHx9ffH09GT8+PEAvP7665w7d45OnTrx5Zdf8vDhQwYNGkSjRo3w9vZm+fLlgDGhCQwMpE2bNrRt25arV6/SsmVLvLy8cHd3Z+vWrcmuP2PGDK5cuULr1q1p3bo1APPmzcPDwwN3d3dGjRqV4j1K7RyADz74gPr169OkSROuX78OwM2bN+nevTu+vr74+vqyffv2ZH3FxsYybtw4FixYYBodmjBhAi+//DJNmzalVq1a/PDDD6b3ukWLFgQGBuLqmnLR0xtvvIGPjw9ubm6m9w2MydX48eNp0KABHh4enDhxAoDbt2/ToUMH3NzcGDx4MKkNIo0ePZqoqCi8vLxMI5Dx8fG89tpruLm50aFDB6KijCvXz549S8eOHWnYsCEtWrQwXSfRjRs36N+/P3v37sXLy4uzZ8/i4uLCyJEj8fDwoFGjRpw5cwb45//Dxo0bM3LkyBRxpSkxk88vXw0bNtRCpCneoPXen7WeUl3rCQ5a/zVc64e3M3z6c5PW62dG/ZXiq86HK/Qbv4foN34P0XU+XJFqm+cmrU+zX4PBoGvUqKHPnj2rY2JitKenpz569GiKdg8ePNAtWrTQjRs31nv37k3x/OHDh3WNGjUy/HqEyCvOnz+v3dzcTI83btyoS5UqpS9fvqzj4+N1kyZN9NatW7XWWrdq1Urv3btXh4SE6Hbt2pnOuXv3bop+X3nlFe3v76/j4+P1qVOntKOjo46KitLTpk3TAwcO1Fprffz4ce3s7KyjoqL0xo0bdUBAgNZa6/Hjx+umTZvq6OhoffPmTV22bFkdGxubItYnX4dSSu/cuVNrrfXq1av1a6+9phMSEnR8fLwOCAjQmzdv1lpr/cwzz+ibN29qrbUeM2aM/u2330yvo1atWjoyMlL/8ssv2tHRUd++bfw5NW3aND1x4kSttfHnxoMHD1LEkLTf8PBw7ezsrG/cuKHj4uJ069at9bJly9I9R2utAR0cHKy11nrEiBH6k08+0Vpr/dJLL5k+h4sXL+q6deum6OuXX37RQ4cONT0eP3689vT01I8ePdI3b97UTk5OOjw8XG/cuFEXL15cnzt3LtX3MvE1GwwG3apVK33o0CFTrDNmzNBaaz179mz96quvaq21fuutt/RHH32ktdb6r7/+0kCy15SoRIkSpu/Pnz+vra2t9YEDB7TWWvfs2dP0ObRp00afOnVKa631rl27dOvWrVP0lfT/l8TYEj+fX3/91fTcK6+8ogMCArTBYEjRBxCi08hpZGK+KDgu7YIVI+DaYXimmfHWY2WPTHVxJZW9vQCi4xI4fT3S9H1mzgXYs2cPNWvWpEaNGgD06dOH5cuXp/jrcOzYsYwaNYqpU6em2s+8efPo06fPU1+HEPlBo0aNTLfqvLy8uHDhAs2bNzc9X6NGDc6dO8dbb71FQEAAHTp0SLWfXr16YWVlRa1atahRowYnTpxg27ZtvPXWW4DxtuAzzzzDqVOnUpwbEBBA0aJFKVq0KBUrVjSNCKXnmWeeoUmTJgCsWbOGNWvW4O1t3P8vMjKS06dP07Jly2TnrFmzhuDgYKZNmwYYR+ouXboEQPv27SlbtiwAvr6+DBo0iLi4OLp27YqXl1e6sezduxc/Pz8qVKgAQL9+/diyZQtdu3ZN97wiRYqYRgQbNmzI2rVrAVi3bl2y+XgPHjwgMjLyqXPIunTpgp2dHXZ2drRu3Zo9e/bg4OBAo0aNqF69eqrnLFy4kO+//x6DwcDVq1c5duwYnp6eAHTr1s0U29KlSwHjnmyJ3wcEBFCmTMb2c6xevbrpfWzYsCEXLlwgMjKSHTt20LNnT1O7mJiYDPX30ksvmf777rvvmo737NkTa+vMbfAtSZjI/x5cgbXj4MgiKOUEPX4BtxeztOVEVQe7VDdZdXSwY+3wVkDaG7Gmt21EeHg4zs7OpsdOTk7s3r07WZv9+/dz+fJlAgIC0kzCFixYYLqNIUR+V7RoUdP31tbWGAyGZM+XKVOGQ4cOsXr1ar799lsWLlzIzz//nKIf9cS/9ScfZyeG1CQtWaS1ZsyYMfz73/9O9xytNUuWLKFOneQrm3fv3p2sv5YtW7Jlyxb+/vtvBgwYwPDhw/nXv/6V0ZeTYba2tqb3KenrTkhIYNeuXRQrlrk6uWl9BmmVdzp//jzTpk1j7969lClThgEDBhAdHW16PvFzyehnkp4nP+OoqCgSEhJwcHDI0tSOpK816fdZKWUlc8JE/hUXDVumwUwf4+rHliNh2B5w75blPb9G+NehqE3yfxZPbgkxwr8OdrbW6bbJrISEBIYPH87nn3+eZpvdu3dTvHjxbK30EcJSSpYsSURERKbOuXXrFgkJCXTv3p2JEyeyf//+VNstWrSIhIQEzp49y7lz56hTpw4tWrTgjz/+AODUqVNcunQpRQKUE7H6+/vz888/myaCh4eHc+PGjVTbzZw50zSP6cCBA6n2d/HiRSpVqsRrr73G4MGDU33NSeNr1KgRmzdv5tatW8THxzNv3jxatWqV5dfUoUMHZs6caXqcWpKSWl/Lly8nOjqa27dvs2nTJnx9fdO9zoMHDyhRogSlS5fm+vXrrFy58qmxtWzZkrlz5wKwcuVK7t69m2o7W1tb4uLi0u2rVKlSVK9enUWLFgHGJPnQoUNPjQEwrZJcsGABTZs2zdA5aTFrEqaU6qiUOqmUOqOUGp3K88OVUseUUoeVUuuVUs+YMx5RQGgNJ1bA141hwydQs40x+WrzARTJXlHdrt6OvNW2pulxaltCZGXbCEdHRy5fvmx6HBYWhqPjP+0jIiIIDQ3Fz88PFxcXdu3aRWBgIElXAs+fP980DC5EflOuXDmaNWuGu7t7isnuaQkPD8fPzw8vLy/69+/PpEmTUm1XrVo1GjVqRKdOnfj2228pVqwYb775JgkJCXh4eNC7d2/mzJmTbEQkp2Lt0KEDffv2pWnTpnh4eNCjR49Uk52xY8cSFxeHp6cnbm5ujB07NtX+Nm3aRP369fH29mbBggW88847KdoMGTKEjh070rp1a6pUqcLkyZNp3bo19evXp2HDhnTp0iXdc9IzY8YMQkJC8PT0xNXVlW+//TZFm9atW3Ps2LFk2zZ4enrSunVrmjRpwtixY6latWq610l8jXXr1qVv3740a9Ys3fYA48ePZ8uWLbi5ubF06VKqVauWarshQ4bg6en51K1B/vjjD3766Sfq16+Pm5tbhu8y3L17F09PT6ZPn57qYpHMMNsWFUopa+AU0B4IA/YCL2mtjyVp0xrYrbV+pJR6A/DTWvdOr1/ZoqKQu3kSVo2GsxugQl3oOBmeTf+HSqKMbD0BMGvjaaatNs4dcUynXWYYDAZq167N+vXrcXR0xNfXl7lz5+Lm5pZqez8/P6ZNm4aPj3FVc0JCAs7OzmzdutU0r0wIYVyV9vzzz9OjRw9Lh1JoTZgwAXt7e/7v//7P0qGYXVa2B7FU2aJGwBmt9bnHQcwHugCmJExrvTFJ+11AfzPGI/Kz6PuwaQrs+Q5sSxjrPPq+CtYZ248rcVuJxJ3nw+9FMXrpYSJjDHT2qGJqt+LIVWasP2N6nLhFBZCtRMzGxoZZs2bh7+9PfHw8gwYNws3NjXHjxuHj40NgYGC652/ZsgVnZ2dJwIQQogAx50hYD6Cj1nrw48cvA4211sPSaD8LuKa1nphevzISVsgkJMDBP2D9R/DwFjR8BdqMhRIZ/ysE0p5Mn1GODnZsH90my+cLIYQonPL8Zq1Kqf6AD5DqkjCl1BClVIhSKuTmzZu5G5ywnMt74Mc2EDwMytaAIZvghemZTsAg/e0jPgp0M31l5XwhRM7JSkmenDRhwgTTNhIZkV5ZnKTll5KW2HnuuedM5yZONM+oJ0sGZTbejEqrJFBa5syZw7BhqY6xZLpEUnrlniZNmkTNmjWpU6cOq1evTreft99+O8W1Fy5ciKurK25ubvTt2zdTcZmDOW9HhgPOSR47PT6WjFKqHfAB0EprneomHVrr74HvwTgSlvOhijwl4hqsmwCH5kHJKtDtR/DokaUVj3HxCcxYf5q0/qdxdLDjledcTI+/33Iu09tPCCFyzldffUX//v0pXry42a5hjhI6qfnxxx9TPb5jxw7gnyQsM8nAwYMHCQkJoXPnzhk+J3Fj0LTKK+U1X3/9darlno4dO8b8+fM5evQoV65coV27dpw6dSrVvblCQkJSrJ48ffo0kyZNYvv27ZQpUybVVay5zZyfyF6gllKqulKqCNAHCE7aQCnlDXwHBGqtLf9uCMsyxMC2r2BmQwhdAs2Hw7AQ8OyZpQTs4u2H9Px2JzM3nKGRSxmK2aa/9QSYZ/sJIURKqZWsSau8TqKPP/4YX19f3N3dGTJkSKpla9IqY/RkeaA7d+7QtWtXPD09adKkCYcPHzb1cejQoRQleCIjI2nbtq2plE7SlXQGg4F+/fpRr149evToYRrJSyy/9KTE0ZnRo0ezdetWvLy8+PLLL2nZsmWyLSGaN2+ebNuE1EoGgTE5ebL0UkbLKz18+JCAgADq16+Pu7t7siLVM2fOTFE6KL33LdH58+dNq0U//PDDFM8n+uKLL3B3d8fd3Z2vvvoKSFnuKanly5fTp08fihYtSvXq1alZsyZ79uxJ0W98fDwjRozgs88+S3b8hx9+YOjQoaZNXitWrJhmbLkmra30c+IL6IxxheRZ4IPHxz7GmHQBrAOuAwcffwU/rU8pW5Q/LNsfpp+btF67PC7ns2x/WPonnFyl9XQvrceX0npuH61vn83G9dbpMUsOadexK7XH+FX6z0PhmYop07ELITItvZI1qZWiSXqO1lr379/fVHYnqbTKGD1ZHmjYsGF6woQJWmut169fr+vXr6+1TrsET1xcnL5//77WWuubN2/qZ599VickJOjz589rQG/btk1rrfXAgQP11KlTtdb/lF968nUlltV5siTOnDlz9DvvvKO11vrkyZM6td93qZUMSqv0UkbKKy1evFgPHjzY1N+9e/dM8aZWOiit9y1pXC+88IL+9ddftdZaz5o1K1kZoUQhISHa3d1dR0ZG6oiICO3q6qr379+f4r1KaujQoaaSQ1prPWjQIL1o0aIU7b766iv9xRdfJHuvtda6S5cuesSIEfq5557TjRs31itXrkxxrjlgqbJFWusVwIonjo1L8n07c15fWEZqKxHTXGF46wysHgOn10C5WtB/CdTM3P8WKa8Xzdw9l6lRvgS/DW6M4+NbiV29HTO0wjGj7YQQWZdeyZq0bNy4kc8++4xHjx5x584d3NzceOGFF1K0S62MESQvD7Rt2zaWLFkCQJs2bbh9+zYPHjwAUi/BExAQwPvvv8+WLVuwsrIiPDzcVOLI2dnZtM9V//79mTFjRpa2a+jZsyeffPIJU6dO5eeff2bAgAEZOi+t0ksZKa/UokUL3nvvPUaNGsXzzz9PixYtTP2mVjoovfct0fbt201tXn755VSLim/bto0XX3zRtMt8t27d2Lp1qym+rLpy5QqLFi1i06ZNKZ4zGAycPn2aTZs2ERYWRsuWLTly5AgODg7ZumZ2SNkikeOmrj5pSogSRcXFM3Z5KJfuGIfpixgi8b38E17h8zBYFWVP9f9wqGpvEi7awsXTmbreD1vPpbgeQLQh3pSACSHyjqeVrElNdHQ0b775JiEhITg7OzNhwoQ0z8lsCZ2MnP/HH39w8+ZN9u3bh62tLS4uLqbrZ6dsUlLFixenffv2LF++nIULF7Jv374MnZdW6aWMllfav38/K1as4MMPP6Rt27aMGzcuWb9ZKR2U1fcgPU/b9BqMlQjOnDlDzZrGTbcfPXpEzZo1OXPmDE5OTjRu3BhbW1uqV69O7dq1OX369FN39zen/DFLT+Qraa0kjIg28OXaE1xc/yPddnSlYdhvLI5rTrOH03j5eCOmrT/PF2tPZforIjr1Hw5X76X/Q10IYRnplaxJq7xOYsJTvnx5IiMjWbx4cZr9p1bG6ElJyxpt2rSJ8uXLU6pUKSD1Ejz379+nYsWK2NrasnHjRi5evGjq69KlS+zcuROAuXPnJitCnp7UXuvgwYN5++238fX1TbVAdVbKP0Ha5ZWuXLlC8eLF6d+/PyNGjEizNFSi9N63RM2aNWP+/PkAprap9RMUFMSjR494+PAhy5YtSzYKl5rAwEDmz59PTEwM58+f5/Tp0zRq1ChZm4CAAK5du8aFCxe4cOECxYsX58wZ496PXbt2NY2Q3bp1i1OnTll870UZCRM5Lq0i2G1LXebHCotQ4SFoR1/iO06hh2MDsrvPdYvPNnAllYRLVjQKkTclLVmT9FYe/FNep2rVqmzc+M9+3g4ODrz22mu4u7tTuXLldEcvEssYPXjwwFTG6EkTJkxg0KBBeHp6Urx4cX799VfTc4kleG7dumUqwdOvXz9eeOEFPDw88PHxoW7duqb2derUYfbs2QwaNAhXV1feeOONDL0Pnp6eWFtbU79+fQYMGMC7775Lw4YNKVWqFAMHDkz1nNatWzN58mS8vLwYM2ZMhq4DxvJKx48fN9U6tLe35/fff+fMmTOMGDECKysrbG1t+eabb9LtJ733LdH06dPp27cvU6ZMSbWEEkCDBg0YMGCAKYkaPHjwU29Furm50atXL1xdXbGxsWH27NmmlZGdO3fmxx9/TLdckr+/P2vWrMHV1RVra2umTp1KuXLl0r2muZlts1Zzkc1a874n52iV5z5jiiygu9UmsK8E7T4Cz96QQ8uln7weGFc0Pq2eoxCi4MnvZYyuXLmCn58fJ06cyDdbSoj0WapskSikEhOf//55iMCYv/iP7TKKq1h47h1oOQKKljTL9TJSF1IIIfKq//3vf3zwwQd88cUXkoAVEjISJjItI4Wwd6yaT5VdH1GdK2xVDYlu8wntWzRLo0chhBCiYMrzZYtE/pF46y/8XhSaf7afCDrwuBjC7bNc/bYrz+36NzohgYGxI3g56j3eXhPxT5tcsGrVKurUqUPNmjWZPHlymu2WLFmCUsq0oeIff/yBl5eX6cvKyirZ5olCCCFETpGRMJEpaRXCtlfR/J/dn/RN+JNYbcMMw4v8Et+JuCR3vHOrCHZ8fDy1a9dm7dq1ODk54evry7x583B1dU3WLiIigoCAAGJjY5k1axY+Psn/UDly5Ahdu3bl7NmzZo9ZCCFEwSQjYSLHpNx+QtPFahvrirzHgIRlHCnTjtYxn/N9/AvJErDUzzWPPXv2ULNmTWrUqEGRIkXo06dPshIjicaOHcuoUaNSXTkFMG/ePPr06WPucIUQQhRSkoSJTEm67YO7OsfiIh8xvcjX3LUqB6+uo+E7CyjikPoS4dzaMiI8PBxn539qxzs5OREenvxW6P79+7l8+TIBAQFp9rNgwQJeeukls8UphBCicJMkTGTKCP86VLV9yH9tfiC4yFieUdd4P+F1Tr4QBM6+pjZ5uQh2QkICw4cP5/PPP0+zze7duylevDju7u65GJkQQojCRLaoEBkXH0fXmGACin4KcY/4Kb4Ti0r05c2ODZKtjrT0lhFPK20RERFBaGgofn5+AFy7do3AwECCg4NN88Lmz58vo2BCCCHMSibmi4w5uxFWjYabJ+DZtvxgP4SpIZpTn3aydGQpGAwGateuzfr163F0dMTX15e5c+fi5uaWans/Pz+mTZtmSsASEhJwdnZm69atFi9pIYQQIn+Tifki6+5egPn94LeuYIiGl+ZD/yXcsXOxcGBps7GxYdasWfj7+1OvXj169eqFm5sb48aNIzg4+Knnb9myBWdnZ0nAhBBCmJWMhInUxT6EbV/C9hlgZQMt34MmQ8HWuJJwyqoT/LT1fJ4cCRNCCCHyCilbJDJOawhdAmvHwYNw8OgF7T+CUmkXRRVCCCFE5kkSJv5x9TCsHAWXdkCV+tDjZ6jWxNJRCSGEEAWSJGECHt6GjRNh3xywKwMvTAfvl8HK+qmnCiGEECJrJAkrzOINEPIzbPwUYiKg0b/Bb5QxERNCCCGEWUkSVlid32K89XjjGFRvBZ2mQMV6lo5KCCGEKDQkCSts7l2CNR/CseXgUA16/w51nwelLB2ZEEIIUahIElZYxD6C7dNh+1eAgtYfwHNvgW3u1HMUQgghRHKShBV0WhtHvdZ8CPcvg1s3aP8xODg//VwhhBBCmI0kYQXZ9aPGeV8XtkIld3jxW3BpbumohBBCCIEkYQXTozuwaRLs/RGKlYaAz6HBALCWj1sIIYTIK+S3spkEHQhn6uqTXLkXRVUHO0b416Grt6N5L5oQb9zra8NEiL4HPq9C6/eheNkcvUzQgXB+23mR2PgEmk3ekDuvTQghhChgJAkzg6AD4YxZeoSouHgAwu9FMWrJYW5HxtDetbJZrlk0fBdltoylyK1Qoh2f427LT4gr7wpRQNSjHLvO2mPX+Gz1SWIMCYDxtY1ZegRAEjEhhBAiE6SAtxk0m7yB8HtRuXKtKtxmjO1cAq13Eq7L8WlcP1YkNAZyd8sJRwc7to9uk6vXFEIIIfI6KeCdy66kk4B93rN+jlzDKj6GWmfnUPfMDyidwLGab3Lq2UG0tbGjbY5cIXXvLTqU6vH0XrMQQgghUpIkzAyqOtilOhLm6GBH94ZO2etcazjxN6x+H+5dhHqB0GEirmWewTV7PWfIF2tPpfraqjrIfmNCCCFEZlhZOoCCaIR/Hexskxe/trO1ZoR/nex1fOME/NYVFvQD2+Lwr2Do/RuUeSZ7/WaC2V6bEEIIUcjISFgWPG3lY1dvR0Iu3uH3XZcAsFaK7g0dsz5xPeoebJ4Cu7+DovbQaSr4DLLIlhOJryHXV34KIYQQBYyMhGVS4srH8HtRaP5ZHRh0IDxZmyX7/nkcrzW/LlyOo8uz1KxZk8mTJ6fZ/5IlS1BKERISAgnxxO3+iVeaO+PxymfU+xEmxbwCjYdYdM+vrt6ObB/dhvOTA9g+uo0kYEIIIUQWyEhYJk1dfdK09USiqLh4PgwK5djVBwDM3X0pWRudEM+1lbPxfG0qOyf2wtfXl8DAQFxdk8/iioiIYPr06TRu3BiuhcIP77FozV5irEtz5OA2HpWuhaurKy+98iouLi5mf61CCCGEMB8ZCcuktFYBRsYY+G3nRX7beZHIGEOy52KvnsLGoQp3rMpQpEgR+vTpw/Lly1P0MXbsWEa99RrFIi9B8DCIvIFqNISHlXwwVHAjKiqKIkWKUKpUKbO8NiGEEELkHknCMighQfPt5rOktauao4Mdxz/pyPFPOuL4xEpBQ8RtbEpVMK0gdHJyIjw8PFmb/Xt2cfnABgJOjYKHt8CrP7wVQo/hUylRogRVqlShWrVq/N///R9ly+bsDvhCCCGEyH2ShGXAtfvRvPzzbiavPIGnUymK2SZ/255cHZjaCkJrK5X6CkKtSTj+N8P7tuNz7/PwbGtwbACNBkOREuzZswdra2uuXLnC+fPn+fzzzzl37pxZXqcQQgghco/MCXuKVaHXGL30MDFxCUzp7kEvH2eWH7zy1NWR8M8KwipVqmJ9fafpeFhYGI6OjnDrNKwaTUToWkKvxeC3uBwEbefatWsEBgYSHBzM3Llz6dixI7a2tlSsWJFmzZoREhJCjRo1LPJ+CCGEECJnSNmiJJJuPVG5dDFcyhVn57k7eDiWZnofL2pUsM9SvwaDgdq1a7N+/XocHR3x9WnI3KFNcbu+xLjfl98YaPQaWNsC4Ofnx7Rp0/Dx8WHKlCmcOHGCX375hYcPH+Lr68v8+fPx9PTMyZcuhBBCCDOQskUZ8GTR7av3o7l6P5q2dSvyTf+GFLHJ+p1bGxsbZs2ahb+/P/FR9xnkZsDt6kLGHa+NT7e3CGz6cprnDh06lIEDB+Lm5obWmoEDB0oCJoQQQhQAkoQ9ltrWEwAnrkVkKwFL1NmzIp1HOEH4dXBqBJ0/4+Oq3qm23bRpk+l7e3t7Fi1alO3rCyGEECJvkSTssbS2nsh2YeqI67BuAhyaC/aV4cXvwbMXKJW9foUQQgiRr0kS9lhaRbezXJjaEAu7v4XNn0F8DDR/F1q8B0VLZjNSIYQQQhQEskXFYzlamPr0WvimKawdCy7N4M1d0G6CJGBCCCGEMJGRsMcSi27P3X2JBA1WiswX3b59Fla/D6dWQbma0G8x1GpvvqCFEEIIkW9JEvZYYtHthMc7diRoWLIvHJ9nyj49EYuJgC3TYNfXYF0UOkyERv8GmyLmD1wIIYQQ+ZIkYY+lVZh76uqTaSdhWsPhhbB2HEReA69+0HY8lKyUCxELIYQQIj+TJOyxTK+ODN8PK0dB2B5wbAh9/gCnVPdiE0IIIYRIQZKwxzK8OjLyJqz/CA78DiUqQJevof5LYCVrHIQQQgiRcZI5PPbU1ZHxcbDza5jZEA7Nh+eGwVv7wLufJGBCCCGEyDQZCXsscd7Xx38d487DWCqULMoHnesZj5/dACtHw62TULMddJwM5WtZOGIhhBBC5GdmHcJRSnVUSp1USp1RSo1O5fmiSqkFj5/frZRyMWc8mVX84WWY3w9+exES4uClBcZtJyQBE0IIIUQ2mW0kTCllDcwG2gNhwF6lVLDW+liSZq8Cd7XWNZVSfYApQG9zxZSepAW87YjmX1ELabX2bww2Nti0mwBN3gSbopYITQghhBAFkDlHwhoBZ7TW57TWscB8oMsTbboAvz7+fjHQVinLFFVM3KKiidUxNhT9P96yCeKvhMZ0t5llLDkkCZgQQgghcpA554Q5ApeTPA4DGqfVRmttUErdB8oBt5I2UkoNAYYAVKtWzSzBJm5FcUM7cEWXY2js2+zXtVH3zXI5IYQQQhRy+WJZn9b6e621j9bap0KFCma5RuJWFOd0VbrHfsR+XTvZcSGEEEKInGTOJCwccE7y2OnxsVTbKKVsgNLAbTPGlKYcLeAthBBCCPEU5kzC9gK1lFLVlVJFgD5A8BNtgoFXHn/fA9igtdZmjClNXb0dmdTNA0cHOxTg6GDHpG4emSvgLYQQQgiRQWZLwrTWBmAYsBo4DizUWh9VSn2slAp83OwnoJxS6gwwHEixjUVu6urtyPbRbTg/OYDto9sQPHMsFStWxN3dPd3z9u7di42NDYsXLzYdGzlyJG5ubtSrV4+3334bC+WWQgghhMijzDonTGu9QmtdW2v9rNb608fHxmmtgx9/H6217qm1rqm1bqS1PmfOeDJrwIABrFq1Kt028fHxjBo1ig4dOpiO7dixg+3bt3P48GFCQ0PZu3cvmzdvNne4QgghhMhH8sXEfEtp2bIlZcuWTbfNzJkz6d69OxUrVjQdU0oRHR1NbGwsMTExxMXFUalSJXOHK4QQQoh8RJKwbAgPD2fZsmW88cYbyY43bdqU1q1bU6VKFapUqYK/vz/16tWzUJRCCCGEyIskCcuG//znP0yZMgWrJwp4nzlzhuPHjxMWFkZ4eDgbNmxg69atFopSCCGEEHmRFPDOhpCQEPr06QPArVu3WLFiBTY2Npw+fZomTZpgb28PQKdOndi5cyctWrSwZLhCCCGEyENkJCwbzp8/z4ULF7hw4QI9evTg66+/pmvXrlSrVo3NmzdjMBiIi4tj8+bNcjtSCCGEEMnISFg6XnrpJTZt2sStW7dwcnLio48+Ii4uDoDXX389zfN69OjBhg0b8PDwQClFx44deeGFF3IrbCGEEELkAyq/7V/l4+OjQ0JCLB2GEEIIIcRTKaX2aa19UntObkcKIYQQQliAJGFCCCGEEBaQ725HKqVuAhfNfJnywC0zX0NknnwueY98JnmTfC55j3wmeVNufC7PaK0rpPZEvkvCcoNSKiSt+7fCcuRzyXvkM8mb5HPJe+QzyZss/bnI7UghhBBCCAuQJEwIIYQQwgIkCUvd95YOQKRKPpe8Rz6TvEk+l7xHPpO8yaKfi8wJE0IIIYSwABkJE0IIIYSwgEKdhCmlOiqlTiqlziilRqfyfFGl1ILHz+9WSrlYIMxCJwOfy3Cl1DGl1GGl1Hql1DOWiLMwedpnkqRdd6WUVkrJKjAzy8hnopTq9fjfylGl1NzcjrEwysDPr2pKqY1KqQOPf4Z1tkSchYlS6mel1A2lVGgazyul1IzHn9lhpVSD3Iqt0CZhSilrYDbQCXAFXlJKuT7R7FXgrta6JvAlMCV3oyx8Mvi5HAB8tNaewGLgs9yNsnDJ4GeCUqok8A6wO3cjLHwy8pkopWoBY4BmWms34D+5HWdhk8F/Kx8CC7XW3kAf4OvcjbJQmgN0TOf5TkCtx19DgG9yISagECdhQCPgjNb6nNY6FpgPdHmiTRfg18ffLwbaKqVULsZYGD31c9Fab9RaP3r8cBfglMsxFjYZ+bcC8AnGP1SiczO4Qiojn8lrwGyt9V0ArfWNXI6xMMrI56KBUo+/Lw1cycX4CiWt9RbgTjpNugD/00a7AAelVJXciK0wJ2GOwOUkj8MeH0u1jdbaANwHyuVKdIVXRj6XpF4FVpo1IvHUz+Tx8L2z1vrv3AysEMvIv5PaQG2l1Hal1C6lVHojASJnZORzmQD0V0qFASuAt3InNJGOzP7eyTE2uXERIcxBKdUf8AFaWTqWwkwpZQV8AQywcCgiORuMt1f8MI4Wb1FKeWit71kyKMFLwByt9edKqabAb0opd611gqUDE7mvMI+EhQPOSR47PT6WahullA3GoePbuRJd4ZWRzwWlVDvgAyBQax2TS7EVVk/7TEoC7sAmpdQFoAkQLJPzzSoj/07CgGCtdZzW+jxwCmNSJswnI5/Lq8BCAK31TqAYxvqFwnIy9HvHHApzErYXqKWUqq6UKoJxgmTwE22CgVcef98D2KBlYzVze+rnopTyBr7DmIDJPBfzS/cz0Vrf11qX11q7aK1dMM7TC9Rah1gm3EIhIz+/gjCOgqGUKo/x9uS5XIyxMMrI53IJaAuglKqHMQm7matRiicFA/96vEqyCXBfa301Ny5caG9Haq0NSqlhwGrAGvhZa31UKfUxEKK1DgZ+wjhUfAbjpL4+lou4cMjg5zIVsAcWPV4ncUlrHWixoAu4DH4mIhdl8DNZDXRQSh0D4oERWmsZyTejDH4u7wE/KKXexThJf4D8cW9eSql5GP8gKf94Lt54wBZAa/0txrl5nYEzwCNgYK7FJp+9EEIIIUTuK8y3I4UQQgghLEaSMCGEEEIIC5AkTAghhBDCAiQJE0IIIYSwAEnChBBCCCEsQJIwIUShoJSKV0odTPLlopTyU0rdf/z4uFJq/OO2SY+fUEpNs3T8QoiCp9DuEyaEKHSitNZeSQ8opVyArVrr55VSJYCDSqk/Hz+deNwOOKCUWqa13p67IQshCjIZCRNCCEBr/RDYB9R84ngUcJBcKugrhCg8JAkTQhQWdkluRS578kmlVDmMdS+PPnG8DMaai1tyJ0whRGEhtyOFEIVFituRj7VQSh0AEoDJj8vM+D0+fghjAvaV1vparkUqhCgUJAkTQhR2W7XWz6d1XClVHdillFqotT6Yy7EJIQowuR0phBDp0FqfByYDoywdixCiYJEkTAghnu5boOXj1ZRCCJEjlNba0jEIIYQQQhQ6MhImhBBCCGEBkoQJIYQQQliAJGFCCCGEEBYgSZgQQgghhAVIEiaEEEIIYQGShAkhhBBCWIAkYUIIIYQQFiBJmBBCCCGEBfw/YAeD+rTR8n8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn import metrics\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "x = np.random.randint(40, 400, 100).reshape(-1, 1)\n",
    "y = np.random.randint(0, 2, 100)\n",
    "\n",
    "model = LogisticRegression()\n",
    "model.fit(x, y)\n",
    "probs = model.predict_proba(x)\n",
    "fpr, tpr, thresholds = metrics.roc_curve(y, probs[:,1])\n",
    "\n",
    "# %%\n",
    "plt.subplots(figsize=(10, 6))\n",
    "plt.plot(fpr, tpr, 'o-', label=\"ROC curve\")\n",
    "plt.plot(np.linspace(0,1,10), np.linspace(0,1,10), label=\"diagonal\")\n",
    "for x, y, txt in zip(fpr[::5], tpr[::5], thresholds[::5]):\n",
    "    plt.annotate(np.round(txt,2), (x, y-0.04))\n",
    "rnd_idx = 27\n",
    "plt.annotate('this point refers to the tpr and the fpr\\n at a probability threshold of {}'.format(np.round(thresholds[rnd_idx], 2)),\n",
    "             xy=(fpr[rnd_idx], tpr[rnd_idx]), xytext=(fpr[rnd_idx]+0.2, tpr[rnd_idx]-0.25),\n",
    "             arrowprops=dict(facecolor='black', lw=2, arrowstyle='->'),)\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.xlabel(\"FPR\")\n",
    "plt.ylabel(\"TPR\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## After plotting ROC-Curve - next question is how the threshold relates back to the values of the variable (x) for identification of the cut off.\n",
    "\n",
    "Simple ans is we can not.\n",
    "\n",
    "X was our input matrix on which we performed the prediction. The thresholds are only related to the prediction from the classifier (\"probabilities\" values)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "\n",
    "# C. Computing the best threshold\n",
    "\n",
    "### Compute the best threshold (similarly to ROC curve computation) of probability which gives lowest values of metric A for the given data 5_c.csv\n",
    "\n",
    "we will be predicting label of a data points like this:\n",
    "\n",
    "## $$ y^{pred} = [0 \\text{ if } y_{score} < \\text{ threshold }  \\text{ else } 1 $$\n",
    "\n",
    "$ A = 500 \\times \\text{number of false negative} + 100 \\times \\text{numebr of false positive}$\n",
    "\n",
    "Note 1: in this data we can see number of negative points > number of positive points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "df_5_c.shape  (2852, 2)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>prob</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.458521</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0.505037</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0.418652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>0.412057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0.375579</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   y      prob\n",
       "0  0  0.458521\n",
       "1  0  0.505037\n",
       "2  0  0.418652\n",
       "3  0  0.412057\n",
       "4  0  0.375579"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_5_c = pd.read_csv('./5_c.csv')\n",
    "print('df_5_c.shape ', df_5_c.shape)\n",
    "df_5_c.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "actual_y_train_arr_5_c  [0 0 0 ... 1 1 1]\n",
      "y_proba_5_c  [0.45852068 0.50503693 0.41865174 ... 0.65916054 0.45626546 0.65916054]\n"
     ]
    }
   ],
   "source": [
    "actual_y_train_arr_5_c = df_5_c.iloc[:, 0].values\n",
    "print('actual_y_train_arr_5_c ', actual_y_train_arr_5_c)\n",
    "\n",
    "y_proba_5_c = df_5_c.iloc[:, 1].values\n",
    "print('y_proba_5_c ', y_proba_5_c)\n",
    "\n",
    "unique_probability_thresholds_5_c = np.unique(df_5_c['prob'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# First I will modify the above get_tpr_fpr() function to return the value of A\n",
    "\n",
    "def get_A_metric(y_actual, y_probabilities, threshold):\n",
    "    tp = 0\n",
    "    fp = 0\n",
    "    tn = 0\n",
    "    fn = 0\n",
    "\n",
    "    min_a = float('inf')\n",
    "\n",
    "    for i in range(len(y_probabilities)):\n",
    "        if y_probabilities[i] >= threshold :\n",
    "            if y_actual[i] == 1:\n",
    "                tp += 1\n",
    "            else:\n",
    "                fp += 1\n",
    "\n",
    "        if y_probabilities[i] < threshold:\n",
    "            if y_actual[i] == 0:\n",
    "                tn += 1\n",
    "            else:\n",
    "                fn += 1\n",
    "\n",
    "    A = (500 * fn) + (100 * fp)\n",
    "\n",
    "    return A\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Now a function to traverse the entire unique probability thresholds array\n",
    "# and return the minimum value of A and also the corresponding threshold\n",
    "\n",
    "def get_minimized_a(y_actual, y_probabilities, total_threshold_arr ):\n",
    "  min_a = float('inf')\n",
    "  min_t = 0\n",
    "\n",
    "  for threshold in total_threshold_arr:\n",
    "    a = get_A_metric(y_actual, y_probabilities, threshold)\n",
    "    if a <= min_a:\n",
    "        min_a = min(a, min_a)\n",
    "        min_t = threshold\n",
    "\n",
    "  return min_a, min_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(141000, 0.2300390278970873)\n"
     ]
    }
   ],
   "source": [
    "print(get_minimized_a(actual_y_train_arr_5_c, y_proba_5_c, unique_probability_thresholds_5_c ))\n",
    "# (141000, 0.2300390278970873)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that thats the minimum value of **A** which is 141000\n",
    "\n",
    "and the corresponding threshold is 0.2300390278970873"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>pred</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>101.0</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>120.0</td>\n",
       "      <td>100.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>131.0</td>\n",
       "      <td>113.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>164.0</td>\n",
       "      <td>125.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>154.0</td>\n",
       "      <td>152.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       y   pred\n",
       "0  101.0  100.0\n",
       "1  120.0  100.0\n",
       "2  131.0  113.0\n",
       "3  164.0  125.0\n",
       "4  154.0  152.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from sklearn.metrics import r2_score\n",
    "\n",
    "df_5_d = pd.read_csv('./5_d.csv')\n",
    "df_5_d.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "actual_y_train_arr_5d  [101. 120. 131. ... 106. 105.  81.]\n",
      "predicted_y_arr_5d  [100. 100. 113. ...  93. 101. 104.]\n"
     ]
    }
   ],
   "source": [
    "# print(df_5_d.to_numpy())\n",
    "actual_y_train_arr_5d = df_5_d.iloc[:, 0].values\n",
    "print('actual_y_train_arr_5d ', actual_y_train_arr_5d)\n",
    "predicted_y_arr_5d = df_5_d.iloc[:, 1].values\n",
    "print('predicted_y_arr_5d ', predicted_y_arr_5d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coefficient of determination or $R^2$ calculation\n",
    "\n",
    "[See Definitions of $R^2$ error](https://en.wikipedia.org/wiki/Coefficient_of_determination#Definitions)\n",
    "\n",
    "Coefficient of determination also called as $R^2$ score is used to evaluate the performance of a linear regression model. It is the amount of the variation in the output dependent attribute which is predictable from the input independent variable(s). It is used to check how well-observed results are reproduced by the model, depending on the ratio of total deviation of results described by the model.\n",
    "\n",
    "Assume R2 = 0.68\n",
    "It can be referred that 68% of the changeability of the dependent output attribute can be explained by the model while the remaining 32 % of the variability is still unaccounted for.\n",
    "R2 indicates the proportion of data points which lie within the line created by the regression equation. A higher value of R2 is desirable as it indicates better results.\n",
    "\n",
    "R² is calculated by taking one minus the sum of squares of residuals divided by the total sum of squares.\n",
    "\n",
    "![Imgur](https://imgur.com/X0zWQho.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9563582786990937"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Checking the r-square metric with sklearn\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "sklearn_r2 = r2_score(actual_y_train_arr_5d, predicted_y_arr_5d)\n",
    "sklearn_r2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9563582786990937\n"
     ]
    }
   ],
   "source": [
    "def calculate_r2_score(y_train, y_predicted):\n",
    "    y_train_bar = y_train.mean()\n",
    "    # y_train_bar = np.mean(y_train)\n",
    "\n",
    "    sum_squared_residual = ((y_train - y_predicted)**2).sum()\n",
    "    sum_squared_total = ((y_train - y_train_bar)**2).sum()\n",
    "\n",
    "    return 1 - (sum_squared_residual/sum_squared_total)\n",
    "\n",
    "print(calculate_r2_score(actual_y_train_arr_5d, predicted_y_arr_5d))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## So the r-square measures matches between my custom function and sklearn's inbuilt funcion\n",
    "\n",
    "---\n",
    "\n",
    "# Mean Squared Error (MSE) / Mean Squared Deviation (MSD)\n",
    "\n",
    "The Mean Squared Error measures the average of the errors squared. It basically calculates the difference between the estimated and the actual value, squares these results and then computes their average.\n",
    "\n",
    "Because the errors are squared, MSE can only assume non-negative values. Due to the intrinsic randomness and noise associated with most processes, MSE is usually positive and not zero.\n",
    "\n",
    "![Imgur](https://imgur.com/GWH0ap5.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "177.16569974554707\n"
     ]
    }
   ],
   "source": [
    "def calculate_mse(y_actual, y_predicted):\n",
    "    mse = np.mean((y_actual - y_predicted)**2)\n",
    "    return mse\n",
    "\n",
    "print(calculate_mse(actual_y_train_arr_5d, predicted_y_arr_5d))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "177.16569974554707"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Checking the Mean Square Error metric with sklearn\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "sklearn_mse = mean_squared_error(actual_y_train_arr_5d, predicted_y_arr_5d)\n",
    "sklearn_mse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## So the Mean Squared Error measures matches between my custom function and sklearn's inbuilt funcion\n",
    "\n",
    "---\n",
    "\n",
    "# Mean Absolute Percentage Error (MAPE)\n",
    "\n",
    "The Mean Absolute Percentage Error measures the error between actual and forecasted values as a percentage. It achieves so by calculating it similarly to MAE, but also dividing it by the actual value, expressing the result as a percentage.\n",
    "\n",
    "By expressing the error as a percentage, we can have a better understanding of how off our predictions are in relative terms. For instance, if we were to predict next year’s spending, an MAE error of $50 could be both a relatively good or bad approximation.\n",
    "\n",
    "![Imgur](https://imgur.com/zCFJTE6.png)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12.912029940096867\n"
     ]
    }
   ],
   "source": [
    "# The below will return the % value i.e 12.9 means 12.9%\n",
    "def calculate_mean_absolute_percentage_error(y_actual, y_predicted):\n",
    "    mape = np.mean((np.abs(y_actual - y_predicted)) / np.mean(y_actual)) * 100\n",
    "    return mape\n",
    "\n",
    "print(calculate_mean_absolute_percentage_error(actual_y_train_arr_5d, predicted_y_arr_5d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When one of the actual data-point is zero\n",
    "\n",
    "https://en.wikipedia.org/wiki/Mean_absolute_percentage_error\n",
    "\n",
    "Problems can occur when calculating the MAPE value with a series of small denominators. A singularity problem of the form 'one divided by zero' and/or the creation of very large changes in the Absolute Percentage Error, caused by a small deviation in error, can occur.\n",
    "\n",
    "As an alternative, each actual value (At) of the series in the original formula can be replaced by the average of all actual values (Āt) of that series. This alternative is still being used for measuring the performance of models that forecast spot electricity prices.[2]\n",
    "\n",
    "Note that this is equivalent to dividing the sum of absolute differences by the sum of actual values, and is sometimes referred to as WAPE (weighted absolute percentage error).\n",
    "\n",
    "So in that case the formulate becomes\n",
    "\n",
    "![Imgur](https://imgur.com/KTg47Gk.png)\n",
    "\n",
    "The derivation is as follows.\n",
    "\n",
    "![Imgur](https://imgur.com/Lm45BRD.png)\n",
    "\n",
    "Basically,\n",
    "\n",
    "```python\n",
    "\n",
    "mean(actual_value) = sum(actual_value) / n\n",
    "\n",
    "hence n * mean(actual_value = sum(a)\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12.91202994009687\n"
     ]
    }
   ],
   "source": [
    "def calculate_mean_absolute_percentage_error_2(y_actual, y_predicted):\n",
    "    mape = ((np.sum(np.abs(y_actual - y_predicted))) / np.sum(y_actual)) * 100\n",
    "    return mape\n",
    "\n",
    "print(calculate_mean_absolute_percentage_error_2(actual_y_train_arr_5d, predicted_y_arr_5d))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}