anomalyFinal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anomaly Detection in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Anomaly detection is a technique used to identify unusual patterns that do not conform to expected behavior, called outliers. It has many applications in business, from intrusion detection (identifying strange patterns in network traffic that could signal a hack) to system health monitoring (spotting a malignant tumor in an MRI scan), and from fraud detection in credit card transactions to fault detection in operating environments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this article we will implement the anomaly detection algorithm and apply it to detect failing servers on a network.\n",
    "Before going into the original dataset let's first visualize how to proceed and carry out different functions. Let's start by importing few libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io as sio\n",
    "import math\n",
    "tic = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be using a dummy dataset to check whether all the functions are working properly and we would then use these functions on our original dataset. Let's import the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = sio.loadmat('D:/MLAssignments/AnomalyDetectionScratch/anomalyData.mat')\n",
    "X = dataset['X']\n",
    "Xval = dataset['Xval']\n",
    "yval = dataset['yval']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sio.loadmat() loads our dataset('anomalyData.mat') into the variable dataset. The variable 'X' contains the training dataset, 'Xval' the cross validation set and 'yval' the corresponding output for the 'Xval'. Let's see the array 'X' that we are going to use to fit it in a Gaussian model to detect anomalous examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(307, 2)\n"
     ]
    }
   ],
   "source": [
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see there are 307 training examples and each having 2 features. The features measure the throughput (mb/s) and latency (ms) of response of each server. While your servers were operating, you collected m = 307 examples of how they were behaving,and thus have an unlabeled dataset {x(1), . . . , x(m)}. You suspect that the vast majority of these examples are “normal” (non-anomalous) examples of the servers operating normally, but there might also be some examples of servers acting anomalously within this dataset. Now, let's visualize the dataset to have a clear picture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Throughput(mb/s)')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt8VPWd//HXJxCISMJFESK36K5K0Ba0WRXSdm29VAF1bfvoT7cXpO3Si3ar7e5vtVZt0a3d3a601W5d2oraum5/W7UKUmu03VpBqGBBMQFUDIJELl4giFxCPr8/zpnhzGRmcpLMTG7v5+Mxj5n5nu8558swOZ8536u5OyIiIu0p6e4CiIhI76CAISIisShgiIhILAoYIiISiwKGiIjEooAhIiKxKGCIiEgsChgiIhKLAoaIiMQysFAHNrPxwD3AGKAVWODuPzCzbwF/B+wIs37D3Zdk2P984AfAAOCn7v7d9s559NFHe1VVVX7+ASIi/cCqVat2uvuoOHkLFjCAFuDr7v6smZUDq8ysLtw2392/l21HMxsA/Ag4F9gCPGNmD7t7fa4TVlVVsXLlyjwVX0Sk7zOzTXHzFqxKyt2b3P3Z8HUz0ACMjbn76cBL7r7R3Q8A/w1cXJiSiohIHEVpwzCzKuBUYEWYdKWZPWdmd5rZiAy7jAU2R95vIX6wERGRAih4wDCzocD9wFXuvhv4MfAXwFSgCfj3TLtlSMs4ra6ZzTWzlWa2cseOHZmyiIhIHhQ0YJhZKUGwuNfdHwBw923ufsjdW4GfEFQ/pdsCjI+8HwdszXQOd1/g7jXuXjNqVKx2GxER6YSCBQwzM+BnQIO73xpJr4xkuwRYm2H3Z4ATzOw4MxsEXAo8XIhypq8HovVBREQyK2QvqVrg08DzZrY6TPsGcJmZTSWoYmoEvgBgZscSdJ+d4e4tZnYl8FuCbrV3uvsL+S7g/LoN7N53kBtmTcbMcHfmLa6noqyUq889Md+nExHp1QoWMNz9KTK3RbQZcxHm3wrMiLxfki1vPrg7u/cdZOHSRgBumDWZeYvrWbi0kTm1Vbg7wU2SiIhAYe8wejQz44ZZkwFYuLQxGTjm1FYl7zhEROSwfj01SDRoJChYiIhk1q8DRqLNImre4no1fIuIZNBvA0YiWCTaLF65ZQZzaqtYuLRRQUNEJIN+3YZRUVaa0maRqJ6qKCtVtZSISBrrS7+ka2pqvKOTD6b3hlLvKBHpT8xslbvXxMnbb6ukEtKDg4KFiEhm/T5giIhIPAoYIiISiwKGiIjEooAhIiKxKGCIiEgsChgiIhKLAoaIiMSigCEiIrEoYIiISCyFXKJ1vJn93swazOwFM/tqmP5vZrbOzJ4zswfNbHiW/RvN7HkzW21mHZvvQ0RE8q6QdxgtwNfdvRo4E7jCzCYDdcAp7v5eYANwbY5jfMjdp8ad50RERAqnYAHD3Zvc/dnwdTPQAIx198fcvSXMthwYV6gyiIhI/hSlDcPMqoBTgRVpmz4L/CbLbg48ZmarzGxujmPPNbOVZrZyx44d+SiuiIhkUPCAYWZDgfuBq9x9dyT9OoJqq3uz7Frr7qcBFxBUZ30wUyZ3X+DuNe5eM2rUqDyXXkREEgoaMMyslCBY3OvuD0TSZwOzgE96lgU53H1r+LwdeBA4vZBlFRGR3ArZS8qAnwEN7n5rJP184J+Ai9x9b5Z9jzSz8sRr4DxgbaHKKiIi7SvkHUYt8Gngw2HX2NVmNgO4HSgH6sK0OwDM7FgzWxLuOxp4yszWAH8CHnH3RwtYVhERaUfB1vR296eATMvXLcmQlqiCmhG+3ghMKVTZRESk4zTSW0REYlHAEBGRWBQwREQkFgUMERGJRQFDRERiUcAQEZFYFDBERCQWBQwREYlFAUNERGJRwBARkVgUMEREJBYFDBERiUUBQ0REYlHAEBGRWBQwREQklkKuuDfezH5vZg1m9oKZfTVMH2lmdWb2Yvg8Isv+s8M8L4ZLuoqISDcq5B1GC/B1d68GzgSuMLPJwDXAE+5+AvBE+D6FmY0EbgTOIFjL+8ZsgUVERIqjYAHD3Zvc/dnwdTPQAIwFLgbuDrPdDfxNht0/AtS5+5vu/hZQB5xfqLKKiEj7itKGYWZVwKnACmC0uzdBEFSAYzLsMhbYHHm/JUwTEZFuUvCAYWZDgfuBq9x9d9zdMqR5luPPNbOVZrZyx44dnS2miIi0o6ABw8xKCYLFve7+QJi8zcwqw+2VwPYMu24BxkfejwO2ZjqHuy9w9xp3rxk1alT+Ci8iIikK2UvKgJ8BDe5+a2TTw0Ci19Ns4KEMu/8WOM/MRoSN3eeFaSIi0k0KeYdRC3wa+LCZrQ4fM4DvAuea2YvAueF7zKzGzH4K4O5vAjcBz4SPeWGaiIh0E3PP2DTQK9XU1PjKlSu7uxgiIr2Gma1y95o4eTXSW0REYlHAEBGRWBQwREQkFgUMERGJRQFDRERiUcAQEZFYFDBERCQWBQwREYllYHsZzKwG+ABwLPAusBZ4XCOvRUT6l6x3GGZ2uZk9C1wLHAGsJ5go8P1AnZndbWYTilNMERHpbrnuMI4Eat393UwbzWwqcALwaiEKJiIiPUvWgOHuP8q1o7uvzn9xRKSvc3eCyawzv5eeq91GbzP7VzOrMLNSM3vCzHaa2aeKUTiRviJ9ks++NOlnR8yv28C8xfXJf7+7M29xPfPrNnRzySSOOL2kzgtXyptFsLDRicA/FrRUIn2ILpIBd2f3voMsXNqY/DzmLa5n4dJGdu872G+DaG/Sbi8poDR8ngHc5+5v6vZRJJ7oRRLghlmTkxfJObVV/ao6xsy4YdZkABYubUx+JnNqq7hh1uR+8zn0ZnECxiIzW0fQpfbLZjYK2FfYYon0DbpIpkp8HonPAeiXn0NvlatbbSWAu18DTANq3P0gsBe4uL0Dm9mdZrbdzNZG0n4ZWX2v0cwyNpyH254P82lFJOnVokEjob9eJBPVUFHR6jrp2XK1YdxpZsvN7LvAFMAA3P0dd389xrHvAs6PJrj7/3H3qe4+FbgfeCDH/h8K88ZaCUqkp9JFMhBts5hTW8Urt8xgTm1VSpuG9Gy5utVeYGZlwFnAJcD3zOxV4FHgUXfPOf7C3Z80s6pM2yz4afUJ4MOdK7ZI79Da2spNjzQkL5LXz6xOvof+dadhZlSUlaZUxyXuvCrKSvvN59CbdWhNbzM7DriA4M5hjLuf3k7+KmCxu5+Slv5B4NZsdw9m9grwFuDAf7r7ghznmAvMBZgwYcL7Nm3aFPvfI1JI8+s2sHvfQSrKBrJ7X0syWJQPHkjz/hYqykq5+twTYx2rL41d6Ev/lr4g72t6m9kYM7sIOAV4wN0vIpgipLMuA+7Lsb3W3U8jCE5XhAEmI3df4O417l4zatSoLhRJJH+ivaOiwWLh0kaa9wfv4waLvtYtNz04KFj0HnEmH/w8cAPwO4J2jNvMbJ6739mZE5rZQOCjwPuy5XH3reHzdjN7EDgdeLIz5xPpqs78Is5X7yh1y5WeJE632n8ETnX3NwDM7ChgGdCpgAGcA6xz9y2ZNprZkUCJuzeHr88D5nXyXCLtyhUQEtVKiYt84td93OqkTF1I2zt/VHuBpyPHEumqOFVSW4DmyPtmYHN7O5nZfcDTwElmtsXMPhduupS06igzO9bMloRvRwNPmdka4E/AI+7+aIxyinRYruqeroxMnl+3gW8//ALfXvRCSvqZtzzBvEUdq14yM66fWZ2Sdv3Mar7/+It9qqpKer6sdxhm9rXw5WvACjN7iKAR+mKCC3lO7n5ZlvTLM6RtJRhJjrtvJOjGK1JQ7VX3JNKgY9VK7s4f1m9n9ZZdAFw+fSLe6ty9/FW27d7PwmWNOM6NF54cq3ppft166uq3paTN/OEfGVw6gNWbd4HDDReqqkoKL2svKTO7MdeO7v7tgpSoC2pqanzlSo3zk/iidw0J6QHB3Tnu2iXJ7a/cMiPjxThxkXZ35i2qZ+GyxjZ5Jo0ZyuCBA1gTBpPo+SC1Abi1tRWAWbc9RX1TM5MryzmnejS/WL6JN/cepHrMUM44/ijuWrYp47HSjyeSSUd6SeUah9HjAoJIvrU3VUW2QXfpeb7/+AZ272tJXqy/OXMSK155g/qmw7W5kyvLWXRlLTc90pASMBLVS7v3HeT6mdWYGbc+tp4n1m3n3MmjOad6NO5OfVNz8ngjh5RyTvVomve3pJTt+pnVyTKv3ryLvz5xVLKtRd1ZpaviLtF6HTAxmt/d31vAcokURa6AkHidbdDd9TOr+cETL/GH9dvZ33KIhtf3JNsTFj/XxM49B1KOW9/UzF9c17Y5buZtT/FXE4dzz/LNLH95J2++c5C33z3A/pbgWIuurKWuPnVyhTf3HuS237/c9lg//CNgNLweBJap44fh7syv25DszltSUtLhxnvpmYr9IyBOL6l7CXpKPQ+0FqwkIkWWPlVFtA0DgjuNirKBzKmtonzwQG56pIFvzpgEQEXZQG5a3MCfX30r2VZRPWZoSvVQLpdPn0jroVbuWbGZhqZmXt+1j5OOGULD63uSecoGlmQNMumqxwyl4fU9Kfu/59ihQFC19XjDNuqbmlmx8Y3knYnaO3q3rvbg64w4AWOHuz9ckLOLdKP2pqoIqola+OaMSdy8ZB0LlzayYuMbnD1pFLv3tbBwWSNzplcxZfww7n761ZSLdcLl0yfi7vzymc3sazncXuitzjOb3k6+f2vvQd7aezBl330tqb/PjjqylD9942xO/84TvPHO4byTRg/NeO4tb+/n+WWb+NMrb1Lf1MzgMAAlqrUmV5ZTPniggkUv1F3jc9qdGsTMziYYmf0EsD9S4FwTB3YLNXpLZ2S6rYe21VGJxueERPq8xfVZ7yxOGn0k67e9k3xfVlrCvoOHA0H1mPJk9VF7TjpmCG+/e4htzcGf4ZRxFUwZN4x7lrfbyz2r/jrNel8Qp8NGHB1p9I4TMH4BTAJe4HCVlLv7Z2OXqEgUMCSfMv1BRm38zgU5g0W6k445kvXb30lJ+8wZ47lnRdsLfmkJROIKZQONfS2eDDgjjhjIW++2tNmvIyZXlrP4K++npCTWDEHSA8XtwZdLvueSmhLO1TTb3eeEjx4XLES6ItOa25nWsYg6/hu/SQkWk8YMZcrYiuwnyfB3nAgWZaWpf4oH01oL97U41WOG8vnaKqrHlLcJFp255Nc3NXPTIw2aVryX6o5p8+N8z5abWfa/GpFebn7d+uQfmrvT2toajphezyU/WpqSd9KYocyeNiElbeSQoCnw9IkjOHAo+x9rtGoq3b6DrUwaPZRJo4/MmmfHngN87byTKB3QdtvQwdn/lEcMCVZZHjwwNc/kynIqytSG0Rt119oicRq93w/MDqcc30/wO8nVrVZ6m2hbReIPan7dhmQPoj+/+jZTxlXwTONb1Dc1p7QvTB03jC1vv8u61/ewLq2BucSM2dMmcPfTOZeIyWnQAGPdtrYN1wCfOn0sv/jTa+zcc4Djv/GbNttLgN37s3dgnHnKMTxav4Odew4wYkgpJ44eyuRjh7FwaSNnHH+Uekn1Qt21tkicNoyJmdLdvcctPKE2DMkm2gXx+4+/yO53D9LqrazZspvVm99m5JBS3oz0UppcWc7AEnjutWYmjRnaJki8p3IoG7a/w/4cdxT5MrmyHFoPUb9tb6eP8Z5jy2lphYbXm5k6fjgPfGkaNz3SoHEYvVw+xmHkZaS3mQ119z25AkMiT4dKJ1Jk0S6I7o5hyWk7Zk+bwIGWQym9n4Dk+ynjhnHqhOFtAsbzTcX52g+woCxd/b34/NbD/76mt99N/iLVnUXvVuy1RXJVST1kZquBh4BV7v5OWKDjgQ8RLLH6E+BXBS2hSBelTxEe1V410qkThsfuBVUIiRsYJwgeh/zwczaTRg/NWr0FcNTQQbh7csS3pguRuLK2lLn72QRjL74AvGBmu8zsDeAXQCUw290VLKQoMvVi6oj2ejxls/zlnR3ep1ASQSJXsPj06eM48/ijsm4fOaSURVfWUlJS0udW8pPCy9lLyt2XuPsn3b3K3Ye5+1HuPt3db3b313PtK5Iv+biwZeqCGFVWWsKnzxiXfD9ySClTxlWwbts7TBo9tPOFL7KSkhLuenoTU8YNo7qyvM32N/ce5OZH1tHa2trptT6k/2q3W62ZPREnLUOeO81su5mtjaR9y8xeM7PV4WNGln3PN7P1ZvaSmV3T3rmk7+rKIkbRYyT2qR6T+eK/72Arv3lhOxBU+by59yCnTRjJ7DMn0PhG5xubi2315re5fNpETp0wnIamzCPIVzbuTN5xJbpiHnftkpQ5tVQtJZlkDRhmVmZmI4GjzWyEmY0MH1XAsTGOfRdwfob0+e4+NXwsSd9oZgOAHwEXAJOByzQOpP/Kx4Ut2gXxkb//AFPHD8uYL9Ht9JAHdxgLlzVy9/JX2dfS2qmBccU0/IigOXLNa7txnOsuOKnNYMCR4XiMpt0Hsg5MVLCQXHL9HXwBWEUwLciz4etVBI3gP2rvwO7+JPBmJ8p0OvCSu2909wPAfxOs8if9VD4ubFedc0Jyau8Hv1ybNd9bew9SVlqS0sU2KETwVGIw6sjS2OcttEljhjJl3DDefreF6spypoyrYPiQQXzn0Q2MHDIoJe8z153N5Mpy/mLU0JQpzqMKPVJYerdcjd4/cPfjgH9w9+MijynufnsXznmlmT0XVlmNyLB9LKlrhm8J06Sf6uiFLT391sfWc8l/LGPe4npaW1vbrLOdbl/6vBxAqx9+3vHOwTbb8+09x5a323Zy1JGlLLpiOqdOGA7AoAElPPCl6cFMuksbGT4kNbDdvGQdD18xnV9+YVq3jRSW3i3OSO9dZvaZ9ER3v6cT5/sxcBNBL8GbgH8H0uelyvSzMeu318zmAnMBJkyYkC2b9FJx1qxIXfnuxeQAvURaYiT36s1vZ1zYqCfaumsfb7xzkGFlA9i171DGPHsPHOLmJeuS/9a7lm1Krp0xubI8uazr4q+8P7nwU2I9jK+dd1K3jBSW3i1OwPiryOsy4GyCKqoOBwx3T65kb2Y/ARZnyLYFGB95Pw7YmuOYC4AFEIz07miZpGfJNA4gzoVtft0Gdr174PCgPAfHWb15V/LCWd/UnBIs3ju2gs1vvdtmHYqe4I13DjL8iIG8nWNG2iGDBnD306/yTONbnFM9OmXbOdXHAIcnGLx+ZjUrNgZLxiamA7n63BNTPm8N5pP2tBsw3P0r0fdmNgz4eWdOZmaV7t4Uvr0EWJsh2zPACWZ2HPAacCnwt505n/QuuVYQi17I0i9siZ5Udy3bxOXTJzJnelVyJDfAnOlVfHPmpDYr1z332m6AjBfmsoElbRYwKpbE1Oa5gkVicN7IIaUpiyIlNO8/xKIra5MLPyXuyNI7CxR7pLD0bnHuMNLtBU5oL5OZ3QecRdDLagtwI3CWmU0lqGJqJGhYx8yOBX7q7jPcvcXMrgR+CwwA7nT33JXO0uu1t4JYuuiFLddIboDysgFcePvSNukJiQvznOlVHDp0iHtWbO62YAFtpzbPJDGSO71xfuN3Lmiz7nj0M9EdhHRFnMkHF3G4DWEAUA38P3fvceMjNPlg79bVFcTSF5NJiE4sOPvMCTy0ZmubX+/VleWcM+kYmve3sHLjTta+nn0q8u703rHlHDxE1lX6EqsABhMLDkw2gEe3K2hIVF4mH4z4XuR1C7DJ3bd0qmQiOSTuFDrzi9jd2/R+unz6xGSbRmKluruXB3NHDRpgKWtXNDQ1s6N5Pzv3HGDQgOJfTEs4vJxlNiccfQRmJTS8vqvNtsunB5NKJz676Drk7XUWEIkrThvGH8xsDMH4CAdeLnippF/K1n22vYtbYr+7lm1i6vhhnDp+BI6ntGmUlw3kh797KblPpoWOEg3iuRZBKpQ4FWAv7nyXSaOzD51KfE4VZaUMGDBAvaAk79oNGGb2eeAG4HcEXV5vM7N57n5noQsn/UdHus+mS19MJj39q2f/JTc90lCMf0bBrcuwat/l0ycy7IhBlJSUpHxO6gUl+RanSuofgVPd/Q0AMzsKWAYoYEjedHUFsfSLI5DcPz0QXfIfS1m9ObVaZ/DAEvZ3Y0N3LkaOgUiAYVx1TtAPpb1eTwoW0hVxAsYWINrC1kzqSGyRvOjqL+JsF8f0QPTAl6Yz84d/pCGyKFJngsUAC6YKidOrqSuctl1/qyvLeSQyIA9Tu4QUXpyA8RqwwsweIvjuXgz8ycy+BuDutxawfNLPFOIXcTQQuTs3PdJAw+t7kkHk4tufSo7JSBhQAgPNci7Beshzr02RL5nGiTQ0NXPT4gaun1UNqF1CiiPOJJwvA7/m8F3xQ0ATUB4+RPKmqwslZRO9a0lv73jfxJFt8h9qhfKy4PdUYl2JwQPjz1k7Kcs06pmUlsCVZx3HiCFtJzWcMrYiGSxmT5vAK7fM4PJpQY+oP29+K3kXpnW5pRji9JL6djEKIpJrpHdnLojZlh9Nr/qqOGJgcuqQhMmV5QwdPJBZUyq4fmY1F96+NGX7iCGlOacUSV8DPJfqygq+/pFqvnr2Cbxn3uMpkx+eNnEkWNCK8a2LTsHMuPGik8Fg2BGDdFchRRVnAaUTzWyBmT1mZr9LPIpROOk/8rFQUlR7q/SlTivSQn1Tc3LW1sunT6S+qZmTxw7jhlmTKSkpSc7NlHBMeerU4ZkMHmBZF2wqKy3hmKGDOOmYITz32m6+9fBaLvqPp9l3sJWp44Ydnj12WSOnTRjBg1+uTblLuvHCk3VXIUUXpw3jf4A7gJ8CmafNFOmi9Ok9ss19FEd704ykN6xHq6i+//iLrH71barHDKWibCBmxrcffoHFzzelnGN9hu6t6fYf8pSG9YTE3cm+gwfYvucA1WOGsmbLruQKeVPCBZ6ivcRKSlJ/2+nOQrpDnKlBVrn7+4pUni7R1CC9X/r0Hq/cMqNTF8ds04wkFlHKlB9g3qL65MSFl0+fCA53Pb0JCNoyFl9Zyxm3/K7TU6RPrizn7EnHsGf/If786lus3pLavXf2tAnJqqdEuRQcpJDyMjVIuDwrwCIz+zLwILA/sd3dO7OankhWnR3pnYmZUT449et9eI6ltm0iiePfcGF4l7OskbuWbUpuT3RjNTPGDj8iJWBUjyln4449OXtUQRCAElVciQCVPvdVNFhEyyXSE+Rqw1gFrARmEwzeW8bhZVr1M17yKt8rwLW2tvJ4w7aUtFm3PdVum4iZJYNGVCJYzFtcz5otu7h82sRkGRteb6a6sqLdMj3059f4/uMbku/nLapvk+fbi17QanfSY2W9wwiXZxUpiq6O9I5KjLWILpwEJN9/c8akrMdz94wX8nmL67nxwpMzl9GDLq5RJXZ4WdeEt95toa5+G189+wRufmRdStVXYpLExF3NjReerLsL6XHizCX10QzJu4Dn3X17/osk/VW+5j6KBp/ywQNSusOeUz2am5esy1gtlQgW7V3IE+dI7heu7BecbyB1DduSDdjpBpcOoKSkhIojSpk6fjhTxw9LHhOCwKPustJTxekl9TlgGvD78P1ZwHLgxHASwoyr75nZncAsYLu7nxKm/RtwIXCAYEDgHHd/O8O+jQRTkBwCWuI2yEjvl6+R3lefeyKHDh1qs3BSYn3v9N5SiXNVHFHK6IrBjDxyUMp4kMXPN9HQ1JyxfMOOGJRy13HVOSfw0R8/jbuzJtKonZgkMDEWJH3+p0RVmIKF9FRxhq62AtXu/jF3/xgwmaDx+wzgn3LsdxdwflpaHXCKu78X2ABcm2P/D7n7VAUL6Qx35+Yl65LVUAmJ99fPrM54Yb7qnBOYcUplMPVGOMOtmbFzzwEmH1uRsX3h6nNPTLkTKikp4YEvTePUCcMzHj/BzNoELAUL6cniBIwqd4+2Hm4HTgx7SWUd6uruTwJvpqU95u6JSXGWA+M6WF6RWKLVUou/8v6UbedUj87YtTaxX3lZMPJ74dJGjrt2CQuXNTK5spzywQNzTrGekGhDuWvZppQG/LuWbepUA75ITxGnSuqPZraYYAAfwMeAJ83sSKBNdVIHfBb4ZZZtDjxmZg78p7sv6MJ5pJ+6+twTaW1tbbMWRvP+lqzjG9yd5v0tKe0eENyZnHH8URn3y5SmxYukL4oTMK4gCBK1BFPz3wPc78HPpA915qRmdh3Bcq/3ZslS6+5bzewYoM7M1oV3LJmONReYCzBhwoTOFEf6qMQv/Y4symRmXD+zmhUb32gzt1Smaqxc819Fj6/Fi6QvaLdKygO/cver3f2q8HWn76nNbDZBY/gnsx3H3beGz9sJBgyenqN8C9y9xt1rRo0a1dliSR+UravunNqqrL/0o11yo+rDNo3oV7a9+a8ylUekN4vTrbaZw1ObDwJKgXfcvf2RSm2PdT5BQ/lfu/veLHmOBErcvTl8fR4wr6PnEoGOd9UNgkzm2WsTc0tF8+Zr/iuR3iDOHUa5u1eEjzKC6qnb29vPzO4DngZOMrMtZva5cL9ygmqm1WZ2R5j3WDNLzJEwGnjKzNYAfwIecfdHO/WvE6FjXXUzzV47p7aK+qZmdu9radNgHQ0aCQoW0lfFacNI4e6/NrNrYuS7LEPyz7Lk3QrMCF9vBKZ0tFwi+dDREef5nP9KpKfr6EjvEqCG3GvSi/Rqcaux0ue/itOoLtKbxbnDuDDyugVoJFjXW6TPilONlc/5r0R6g3bXw+hNtB6GdIdsS8GK9AYdWQ8jzhKt48zsQTPbbmbbzOx+M9MIbZFQvua/Eunp4kwNshB4GDgWGAssCtNERKQfiRMwRrn7QndvCR93ARohJyLSz8QJGDvN7FNmNiB8fAp4o9AFExGRniVOwPgs8AngdaAJ+HiYJiIi/UjObrVmNgD4mLtfVKTyiIhID5XzDsPdD6ExFyIiQryBe0vN7HaCtSveSSS6+7MFK5WIiPQ4cQLG9PA5OmOsAx/Of3FERDpHAygLr92A4e6dWiRJRKRYci1kdfW5J3Z38fqMOJMPDiaY0rwqmt/dtUaFiHS76EJWQMokkHNqq3SnkUdxqqQeAnYBq4D9hS2OiEjHaCGr4okTMMa5+/kFL4mISCclgkaZ44xCAAANgklEQVQiWICmly+EOAP3lpnZezpzcDO7M5y0cG0kbaSZ1ZnZi+HziCz7zg7zvBiuAy4iklG2haz60mzcPUHWgGFma83sOeD9wLNmtt7MnjOz58P0OO4C0u9OrgGecPcTgCfC9+nnHgncCJwBnA7cmC2wiEj/lr6QVWJZ3YVLGxU08ixXldRYYGpXDu7uT5pZVVryxcBZ4eu7gf8F/iktz0eAOnd/E8DM6ggCz31dKY+I9D1ayKp4cgWMV9x9UwHOOdrdmwDcvcnMjsmQZyywOfJ+S5jWhpnNBeYCTJgwIc9FFZHeIO6yutI1uQLGMWb2tWwb3f3WApQnIdP/csb7SndfACyAYMW9ApZJRHowLWRVeLkavQcAQ4HyLI/O2mZmlQDh8/YMebYA4yPvxwFbu3BOERHpolx3GE0FGpz3MDAb+G74/FCGPL8FvhNp6D4PuLYAZRHJSdNNiByW6w6jy38VZnYf8DRwkpltMbPPEQSKc83sReDc8D1mVmNmPwUIG7tvAp4JH/MSDeAixTK/bkNKL5tEb5z5dRu6uWQi3SPXHcbZXT24u18W99juvhL4fOT9ncCdXS2DSGdougmRtrIGDP2il/5M002ItBVnpLdIvxQNGgkKFtKfKWCIZKHpJkRSKWCIZKDpJkTaijNbrUi/o+kmRNqyvvRLqaamxleuXNndxZA+ROMwpK8zs1XuXhMnr6qkRHLQdBMihylgiIhILAoYIiISiwKGiIjEooAhIiKxKGCIiEgsChgiIhKLAoaIiMSigCEiIrEUPWCY2Ulmtjry2G1mV6XlOcvMdkXy3FDscoqISKqizyXl7uuBqQBmNgB4DXgwQ9Y/uvusYpZNRESy6+4qqbOBl919UzeXQ0RE2tHdAeNS4L4s26aZ2Roz+42ZnVzMQomISFvdFjDMbBBwEfA/GTY/C0x09ynAbcCvcxxnrpmtNLOVO3bsKExhRUSkW+8wLgCedfdt6Rvcfbe77wlfLwFKzezoTAdx9wXuXuPuNaNGjSpsiUVE+rHuDBiXkaU6yszGWDiPtJmdTlDON4pYNhERSdMtK+6Z2RDgXOALkbQvArj7HcDHgS+ZWQvwLnCp96WVnkREeqFuCRjuvhc4Ki3tjsjr24Hbi10uERHJrrt7SYmISC+hgCEiIrEoYIiISCwKGCIiEosChoiIxKKAISIisShgiIhILAoYIiISiwKGiIjEooAhIiKxKGCIiEgsChgiIhKLAoaIiMSigCEiIrEoYIiISCwKGCIiEku3BQwzazSz581stZmtzLDdzOyHZvaSmT1nZqd1RzlFRCTQLSvuRXzI3Xdm2XYBcEL4OAP4cfgsIiLdoCdXSV0M3OOB5cBwM6vs7kKJiPRX3RkwHHjMzFaZ2dwM28cCmyPvt4RpKcxsrpmtNLOVO3bsKFBRRUSkOwNGrbufRlD1dIWZfTBtu2XYx9skuC9w9xp3rxk1alQhyikiInRjwHD3reHzduBB4PS0LFuA8ZH344CtxSmdiIik65aAYWZHmll54jVwHrA2LdvDwGfC3lJnArvcvanIRRURkVB39ZIaDTxoZoky/Je7P2pmXwRw9zuAJcAM4CVgLzCnm8oqIiJ0U8Bw943AlAzpd0ReO3BFMcslIpIP7k74gzjj+96qJ3erFRHpdebXbWDe4nqC37xBsJi3uJ75dRu6uWRdp4AhIpIn7s7ufQdZuLQxGTTmLa5n4dJGdu87mAwivVV3j/QWEekzzIwbZk0GYOHSRhYubQRgTm0VN8ya3OurpXSHISKSR9GgkdAXggUoYIiI5FWiGioq2qbRmylgiIjkSbTNYk5tFa/cMoM5tVUpbRq9mdowRETyxMyoKCtNabNIVE9VlJX2+mop6+0RL6qmpsZXrmyztIaISFH1pnEYZrbK3Wvi5FWVlIhInqUHh54aLDpKAUNERGJRwBARkVgUMEREJBYFDBERiaVP9ZIysx3Apjwd7mhgZ56OVSgqY36ojF3X08sHKmM2E9091nKlfSpg5JOZrYzb1ay7qIz5oTJ2XU8vH6iM+aAqKRERiUUBQ0REYlHAyG5BdxcgBpUxP1TGruvp5QOVscvUhiEiIrHoDkNERGLp9wHDzBrN7HkzW21mbWYutMAPzewlM3vOzE4rcvlOCsuWeOw2s6vS8pxlZrsieW4oQrnuNLPtZrY2kjbSzOrM7MXweUSWfWeHeV40s9lFLuO/mdm68P/yQTMbnmXfnN+LApfxW2b2WuT/c0aWfc83s/Xhd/OaIpbvl5GyNZrZ6iz7FuszHG9mvzezBjN7wcy+Gqb3mO9jjjL2qO9ju9y9Xz+ARuDoHNtnAL8BDDgTWNGNZR0AvE7QbzqafhawuMhl+SBwGrA2kvavwDXh62uAf8mw30hgY/g8Inw9oohlPA8YGL7+l0xljPO9KHAZvwX8Q4zvwsvA8cAgYA0wuRjlS9v+78AN3fwZVgKnha/LgQ3A5J70fcxRxh71fWzv0e/vMGK4GLjHA8uB4WZW2U1lORt42d3zNTix09z9SeDNtOSLgbvD13cDf5Nh148Ade7+pru/BdQB5xerjO7+mLu3hG+XA+MKce64snyOcZwOvOTuG939APDfBJ9/XuUqnwVTsH4CuC/f5+0Id29y92fD181AAzCWHvR9zFbGnvZ9bI8CBjjwmJmtMrO5GbaPBTZH3m8J07rDpWT/45xmZmvM7DdmdnIxCxUx2t2bIPgDAY7JkKcnfZ6fJbh7zKS970WhXRlWU9yZpSqlJ3yOHwC2ufuLWbYX/TM0syrgVGAFPfT7mFbGqJ78fQS04h5ArbtvNbNjgDozWxf+qkrINJF90buWmdkg4CLg2gybnyWoptoT1nf/GjihmOXrgJ7yeV4HtAD3ZsnS3veikH4M3ETwudxEUO3z2bQ8PeFzvIzcdxdF/QzNbChwP3CVu++2eGtQFPVzTC9jJL0nfx+T+v0dhrtvDZ+3Aw8S3OpHbQHGR96PA7YWp3QpLgCedfdt6Rvcfbe77wlfLwFKzezoYhcQ2Jaorguft2fI0+2fZ9iwOQv4pIcVxOlifC8Kxt23ufshd28FfpLl3N36OZrZQOCjwC+z5SnmZ2hmpQQX4nvd/YEwuUd9H7OUscd/H6P6dcAwsyPNrDzxmqABam1atoeBz1jgTGBX4ja3yLL+mjOzMWF9MmZ2OsH/6xtFLFvCw0Cil8ls4KEMeX4LnGdmI8KqlvPCtKIws/OBfwIucve9WfLE+V4UsozRNrJLspz7GeAEMzsuvPu8lODzL5ZzgHXuviXTxmJ+huF3/2dAg7vfGtnUY76P2crYG76PKbq71b07HwQ9TNaEjxeA68L0LwJfDF8b8COCHinPAzXdUM4hBAFgWCQtWsYrw/KvIWg4m16EMt0HNAEHCX6lfQ44CngCeDF8HhnmrQF+Gtn3s8BL4WNOkcv4EkGd9erwcUeY91hgSa7vRRHL+PPwu/YcwUWvMr2M4fsZBL1tXi5UGTOVL0y/K/H9i+Ttrs/w/QTVSM9F/l9n9KTvY44y9qjvY3sPjfQWEZFY+nWVlIiIxKeAISIisShgiIhILAoYIiISiwKGiIjEooAh/YqZ7elA3rPMbHohy5PhnN83sw/m4TiPZ5udVaSzFDBEsjsLKFrAMLORwJmenykffg58OQ/HEUlSwJB+z8wuNLMVZvbn8Jf56HCCuC8CV4drEHzAzEaZ2f1m9kz4qA33/1Y4SeD/mtlGM/v7yLE/E04iuMbMfm5m5Wb2SjhNBGZWEa51UAp8HHg0sm+jmX3HzJ42s5VmdpqZ/dbMXjazL4Z5Ks3sybCMa83sA+HuDxPMDiCSN5p8UASeIvhl72b2eeD/uvvXzewOYI+7fw/AzP4LmO/uT5nZBIIpJKrDY0wCPkSw1sF6M/sxcCJwHcHEcTvNbKS7N5vZ/wIzCSaJvBS4390PhgHoV2ll2+zu08xsPsHo6lqgjGDE7x3A3wK/dfd/NrMBBLMC4O5vmdlgMzvK3btjmhjpgxQwRIIJ534ZzuE0CHglS75zgMmRWVArEnP8AI+4+35gv5ltB0YDHwZ+5e47Adw9sa7ET4H/SxAw5gB/F6ZXAjvSzpmYH+p5YKgHayk0m9k+C1Znewa4M7xD+bW7R1e/204wxYQChuSFqqRE4Dbgdnd/D/AFgl/wmZQA09x9avgYG17AAfZH8h0i+DFmZJgq292XAlVm9tfAAHdPTCT3boZzJ47bmnaOVoKV2p4kWBXvNeDnZvaZSJ6y8JgieaGAIQLDCC64cHh2U4BmgiqmhMcIJnoEwMymtnPcJ4BPmNlRYf6RkW33EEzstzCS1gD8ZUcKbmYTge3u/hOC2VBPC9MNGEOwtKdIXihgSH8zxMy2RB5fI1hD+3/M7I/AzkjeRcAliUZv4O+BmrARu56gUTwrd38B+GfgD2a2BohOvX0vwRrS0SnrHyHomdURZwGrzezPwMeAH4Tp7wOW++HlP0W6TLPVinQDM/s4cLG7fzot/Slglru/3cXj/wB42N2f6MpxRKLU6C1SZGZ2G8EKijMybP46MAHoUsAA1ipYSL7pDkNERGJRG4aIiMSigCEiIrEoYIiISCwKGCIiEosChoiIxKKAISIisfx/TENpNI6rHlAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[:, 0], X[:, 1], marker = \"x\")\n",
    "plt.xlabel('Latency(ms)')\n",
    "plt.ylabel('Throughput(mb/s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian distribution\n",
    "To perform anomaly detection, you will first need to fit a model to the data’s distribution.\n",
    "Given a training set {x(1), ..., x(m)} (where x(i) ∈ R^n, here n = 2), you want to estimate the Gaussian distribution for each of the features. For each feature (i = 1 . . . n), you need to find parameters mean and variance(mu, sigma^2). For doing that let's write down the function that calculates the mean and variance of the array(or you can call it matrix) X.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimateGaussian(X):\n",
    "    n = np.size(X, 1)\n",
    "    m = np.size(X, 0)\n",
    "    mu = np.zeros((n, 1))\n",
    "    sigma2 = np.zeros((n, 1))\n",
    "    \n",
    "    mu = np.reshape((1/m)*np.sum(X, 0), (1, n))\n",
    "    sigma2 = np.reshape((1/m)*np.sum(np.power((X - mu),2), 0),(1, n))\n",
    "    \n",
    "    return mu, sigma2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu, sigma2 = estimateGaussian(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean:  [[14.11222578 14.99771051]]  variance:  [[1.83263141 1.70974533]]\n"
     ]
    }
   ],
   "source": [
    "print('mean: ',mu,' variance: ',sigma2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the mean and variance let's find out the probabilities of the dataset by writing down the multivariate Gaussian function. But before that let's look at what actually the function does."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multivariate Gaussian Distribution\n",
    "The multivariate gaussian is used to find the probabilities of each training example. It takes the mean(mu) and the co-variance matrix as it's parameters and returns the probabilities of the dataset passed. Let's see how to do that. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multivariateGaussian(X, mu, sigma2):\n",
    "     n = np.size(sigma2, 1)\n",
    "     m = np.size(sigma2, 0)\n",
    "     #print(m,n)\n",
    "     \n",
    "     if n == 1 or m == 1:        \n",
    "         sigma2 = np.diag(sigma2[0, :])\n",
    "     \n",
    "     X = X - mu\n",
    "     pi = math.pi\n",
    "     det = np.linalg.det(sigma2)\n",
    "     inv = np.linalg.inv(sigma2)\n",
    "     val = np.reshape((-0.5)*np.sum(np.multiply((X@inv),X), 1),(np.size(X, 0), 1))\n",
    "     \n",
    "     p = np.power(2*pi, -n/2)*np.power(det, -0.5)*np.exp(val)\n",
    "     \n",
    "     return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inside the function first we convert the sigma2 vector into a covarince matrix and then we simply apply the formula for the multivariate distribution to get the probability vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = multivariateGaussian(X, mu, sigma2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "p contains the proobabilities of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(307, 1)\n"
     ]
    }
   ],
   "source": [
    "print(p.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In an anomaly detection we use a threshold(epsilon) value of probability to determine whether an example is anomalous or not.\n",
    "\n",
    "P > epsilon (negative or normal case),\n",
    "P < epsilon (positive or anomalous case)\n",
    "\n",
    "To determine this threshold value epsilon we need a cross-validation set which we already have collected from our anomalyFile above. Let's see how we go about doing this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "pval = multivariateGaussian(Xval, mu, sigma2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pval contains the probabilities for the cross-validation vector that we are gonna use to calculate threshold. Now, we are going to use this pval vector to calculate the threshold value. Precisely, we need a function that can analyze some labeled datasets and give us the threshold on the basis of some underlying theory. Let's see this function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def selectThreshHold(yval, pval):\n",
    "    \n",
    "    F1 = 0\n",
    "    bestF1 = 0\n",
    "    bestEpsilon = 0\n",
    "    \n",
    "    stepsize = (np.max(pval) - np.min(pval))/1000\n",
    "        \n",
    "    epsVec = np.arange(np.min(pval), np.max(pval), stepsize)\n",
    "    noe = len(epsVec)\n",
    "    \n",
    "    for eps in range(noe):\n",
    "        epsilon = epsVec[eps]\n",
    "        pred = (pval < epsilon)\n",
    "        prec, rec = 0,0\n",
    "        tp,fp,fn = 0,0,0\n",
    "        \n",
    "        try:\n",
    "            for i in range(np.size(pval,0)):\n",
    "                if pred[i] == 1 and yval[i] == 1:\n",
    "                    tp+=1\n",
    "                elif pred[i] == 1 and yval[i] == 0:\n",
    "                    fp+=1\n",
    "                elif pred[i] == 0 and yval[i] == 1:\n",
    "                    fn+=1\n",
    "            prec = tp/(tp + fp)\n",
    "            rec = tp/(tp + fn)\n",
    "            F1 = 2*prec*rec/(prec + rec)\n",
    "            if F1 > bestF1:\n",
    "                bestF1 = F1\n",
    "                bestEpsilon = epsilon\n",
    "        except ZeroDivisionError:\n",
    "            print('Warning dividing by zero!!')          \n",
    "       \n",
    "    return bestF1, bestEpsilon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll implement F1-score principle to determine the best epsilon value. First, we use a variable called stepsize that takes a number of different epsilon values right from the minimum probability in pval to it's maximum. Then, we loop over all such values to find the optimum precision and recall. We are going to need a try except block because there can be cases where we divide by zero to calculate precision and recall. The epsilon value that gives the highest F1-score is taken as the best epsilon. The function returns the bestF1 and bestEpsilon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning dividing by zero!!\n",
      "Epsilon and F1 are: 8.990852779269493e-05 0.8750000000000001\n"
     ]
    }
   ],
   "source": [
    "F1, epsilon = selectThreshHold(yval, pval)\n",
    "print('Epsilon and F1 are:',epsilon, F1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find out the outliers we need to check the examples whose probabilties are less than the threshold value epsilon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "outl = (p < epsilon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to return the indices of the outliers to identify the faulty servers. This gives us a vector with binary entries where 1 means anomaly and 0 means normal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findIndices(binVec):\n",
    "    l = []\n",
    "    for i in range(len(binVec)):\n",
    "        if binVec[i] == 1:\n",
    "            l.append(i)\n",
    "    return l\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Number of outliers: 6\n",
      "\n",
      " [300, 301, 303, 304, 305, 306]\n"
     ]
    }
   ],
   "source": [
    "listOfOutliers = findIndices(outl)\n",
    "count_outliers = len(listOfOutliers)\n",
    "print('\\n\\nNumber of outliers:', count_outliers)\n",
    "print('\\n',listOfOutliers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, we have detected the faulty servers. One last thing in this dummy dataset is to visualize the anomalies graphically. Let's see how to plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VPXZ//H3nQUCEjZBiAoErdTEDdsUNWi1Lq0iarWt2k3AtrbW2tra+thFVLTa9mm17lZbAW0ff13UKotVpFoVEEUEFeJSkdUAgkDCTsj9++PMDDPJzGRIMpPlfF7XNdfMnDlzzp1hOPd8d3N3REQkvPLaOgAREWlbSgQiIiGnRCAiEnJKBCIiIadEICISckoEIiIhp0QgIhJySgQiIiGnRCAiEnIFbR1AJvr16+elpaVtHYaISIfy6quvrnP3/k3t1yESQWlpKfPmzWvrMEREOhQzW5bJfqoaEhEJuawlAjMbZGbPmlmVmS0ysx9Etl9nZqvMbEHkNipbMYiISNOyWTVUB1zp7vPNrBh41cxmRF671d1/m8Vzi4hIhrKWCNy9GqiOPK41syrggGydT0REmicnbQRmVgocDcyNbPqemb1uZg+YWZ9cxCAiIsllPRGYWQ/gEeAKd68B7gEOBoYTlBh+l+J9l5jZPDOb9+GHH2Y7TBGR0MpqIjCzQoIk8Bd3fxTA3de4+253rwfuB0Yke6+73+fuFe5e0b9/k91gRUSkmbLZa8iAPwFV7n5L3PaSuN3OBd7MVgxRvmMHPP00/POfsHEjWp5TRGSPbJYIRgJfB05u0FX0N2b2hpm9DnwG+GEWY+Dh2/9G7f6D8euug3vvxYcO5Ykf3MitM97J5mlFRDqMbPYaehGwJC9Nz9Y5G8WwbRujr7uMK06+lMFjLmD86HLu/MN0Lrzya6wadgR+6iEEBRcRkfDqEFNMNJc9/TQ9jj6SwWMuYOKspUyctRSA4Wd+iUuXz8bs820boIhIO9C5p5ioqcH69WP86PKEzccfV4bV1rZRUCIi7UvnTgSnnII//TS3TH4utqlw9y5W33k/fsYZbReXiEg70qkTgZeU8Mw5F3PBDy7g/jX/5v3SFTz3+DUsKOjDDfWl6j0kIkJnbyMw482vfpsVZcMZ9/4sbM4c9r/xGv7U5VCK9ylSQ7GICGAd4VdxRUWFt2Q9AndPuOg3fC4i0hmZ2avuXtHUfp26aiiq4UVfSUBEZI9QJAIREUlNiUBEJOSUCEREQk6JQEQk5JQIRERCTolARCTklAhEREJOiUBEJOSUCEREQk6JQEQk5JQIRERCTolARCTklAhEREJOiUBEJOSUCEREQk6JQEQk5JQIRERCTolARCTklAhEREJOiUBEJOSUCEREQk6JQEQk5JQIRERCTolARCTklAhEREJOiUBEJOSUCEREQk6JQEQk5JQIRERCTolARCTklAhEREJOiUBEJOSylgjMbJCZPWtmVWa2yMx+ENne18xmmNm7kfs+2YpBRESals0SQR1wpbuXAccCl5lZOXA1MNPdDwFmRp6LiEgbyVoicPdqd58feVwLVAEHAOcAkyO7TQY+n60YRESkaTlpIzCzUuBoYC4wwN2rIUgWwH4p3nOJmc0zs3kffvhhLsIUEQmlrCcCM+sBPAJc4e41mb7P3e9z9wp3r+jfv3/2AhQRCbmsJgIzKyRIAn9x90cjm9eYWUnk9RJgbTZjEBGR9LLZa8iAPwFV7n5L3EtPAGMij8cAj2crBhERaVpBFo89Evg68IaZLYhs+xnwK+BvZvYNYDnwpSzGICIiTchaInD3FwFL8fIp2TqviIjsHY0sFhEJOSUCEZGQUyIQEQk5JQIRkZBTIhARCTklAhGRkGuy+6iZVQAnAPsD24A3gWfc/aMsxyYiIjmQskRgZmPNbD7wU6Ab8DbBdBDHAzPMbLKZDc5NmCIiki3pSgT7ACPdfVuyF81sOHAIwehgERHpoFImAne/K90b3X1ButdFRKRjaLKx2Mx+Y2Y9zazQzGaa2Toz+1oughORDmDrVvyaa+DQQ+GQQ+Cqq/CNG9s6KtkLmfQa+mxkHYHRwEpgGPCTrEYl0oG4e9rnnZo7y0eewuJnX8H/8hd45BH8ww/5YMQJ3PbkoraOTjKUSSIojNyPAh5WbyGRPW6d8Q4Tpi6OXfzdnQlTF3PrjHfaOLLc8H//m6KP1nFW5XeZ8EERfsQRTDj3Sqq37mbg8zPClRQ7sExmH51iZm8RdB39rpn1B7ZnNyyR9s/dqdm+i4mzlgIwfnQ5E6YuZuKspYwbWYq7EyzL0XnZa6/R/7zRjDnhYCbOWhr7LD594smcn/dhp//7O4uUicDMSiIL0F9tZr8Gatx9t5ltJViAXiTUzIzxo8sBEi6C40aWMn50eTgugqWl2PTpjL+lPPb3A5y0vRobOrzt4pK9kq5q6AEze8nMfgUcRWRtAXff4u6rcxKdSDsXnwyiQpMEAM4+G1++nJljf0TXup0U7t7F1+dPpfY/s/Dzz2/r6CRDKROBu58BnAQ8B5wLvGRmj5rZJRpIJhKItgnEi28z6Oy8sJDbrr6bghde4I27vso7936di1fP59zPX8eE/6wIzefQ0aVtI3D37cC/IjfMbChwBnCnmQ109xHZD1Gkfaqvr+eGaVWxNoFrziyLPYfMSgYN2xE6WruCmeFDSvnPbZM58YQDMHdKe/fm01MX07OosEP9LWFmmWRsMxsIjAAceMXdV5tZF3ffme0AASoqKnzevHm5OJVIRm6d8Q4123fRs6iAmu11sSRQ3LWA2h119Cwq5IenDcvoGNGEES1dZPLe9qajJ7TOysxedfeKpvbLZEDZN4GXgfOALxJUEV2cqyQg0t7E9xaKTwITZy2ldkfwvKkLefwxolVJ0R5HNdt3dbgqlYYXfSWBjqXJEoGZvQ1Uuvv6yPN9gdnu/vEcxAeoRCDZ09xfsvEX7qi97S3UGscQSafVSgQEo4lr457XAiuaG5hIrjQ14rclg8HcPWlvoaZiiJeux1GoRytLzqUbR/CjyMNVwFwze5ygjeAcgqoikXarqfr3lgwGu+APc6jdvotPlfZJ2H74tU9xfsUgxp+VWZ2/u3P9lMRpGK6fsohe3Qqp2V7XKdoOpGNI12uoOHL/XuQW9Xj2whFpuUwv8s0ZDFZfX0/t9l0srq5lcXUtFx07iHnLNrK4upYtO3czcXbknGelTizRX/fn3j2LBSs2Mfa4IVx79mFc98SbTJq9jH49urBu807w9McRaS0Z9Rpqa2ojkL2Vaf27uzP0p9Njz9+/eVTC6/X19eTl5cUuwPX19QCMuGlmcLGOKC8p5pODe/HkorUJ26PdSvPyglrYW2e8w6ZtOxk/upzz7pnDghUbKSspprhrAZt31LG4upajDuzF8EG9mDxnz1IfYyuHcO1Zh8WeKxlIJlqz11CFmT1mZvPN7PXorXXCFMmOTEb8phsM5u5c8Ic5jL7jRW55+i0mTF3M7t27GX3Hi1x430ucecTAhPdN+d5ICvILEpIAQHHXAm6YVkV9fT319fVs3LqDSbOXcd49s3n00uM4dGAPqqpreXnpBhZX11JeUsyJw/rz8tINjf6m3z/zDtc98Sbn3j071o4RjTVeR/hxJ+1Lpr2GfgK8AdRHt7v7suyGtodKBLK3mioRuDsTpixm4uylDB/Um+GDegEwafYyxlWWUu/1THtjNes276Rv90I+2rordt+nWwEbttUlPW95STGLq/f0rSgbWEzV6lrKBvbgo627qNtdz/otu5r1N8WqjAhKCEGpYjYAj15aSV5eXmyQm9oTOqAZM/AJE7DXXoODDoIrr8QvuqhFpb/W7DX0obs/4e7vu/uy6K3ZkYlkWXwSGDeylPdvHsW4kaUJffZ//8y7vLZiA2MrhzB8UC8mzV7G3CXrGX5gL15bsYHJc5Yz+siS2AUciN1Hk0DZwGLe++XpCedeXF1L2cDi2POq1bX06V5A1erNrKnZwfotu+i3T5dm/V3RJFBUmMf8ZRs49+7ZLFixiQUrNnHm7S9wy9NvMfqOFzvsWIRQe/ZZtlzwFf4+8jx81Sq4/Xb8V79i+nfH52RK80ymob7WzP4IzAR2RDe6+6NZi0qkBcyMnkWFCSWAaDVRz6JgeY2a7btYsGITwwf1Zvzocl5+/6OEX/LRuv3rpyyiavXmRucoG1hMxZBeVP7q3wnb+3QvpGJILz6s3cG6LcGFe8PWxNJDdHtU1wLjmyOHcP+Ly9i5u+mL9/Zd9by+qmbPObsFiSYaZ3lJMdecWaZ2hA7Eb7qJp8b9mKsYRtXzKxk/+kT+8J0bOecX3+bVz30p650EMqka+jNwKLCIPVVD7u4XZy2qBlQ1JM2RbrBYsqqjeEtuOoPrpyxKaLCN17XA2FG35//OsP7dWb5xO9t31SfdvylFBXlsr6unS75llAzSWXLTGbHGaekg9t8fnzuXCQtqEr6TVXd+haKl72H9+jXrsK1ZNXSUu1e4+xh3Hxe55SwJiDRXw+6a8c+TNSbHO+hnT8aSQNeCvIR7ICEJfHy/7rzz4daMkkCX/MTnRQV5XHTMIHp1C0oqLU0CADdMq1K1UEczbBj28ssJ38mD16+gqEc3rHfvrJ8+k0Twkpml/h8j0o7Fjx6O3oLRw29z7l2zEvYtKylmzHGJM6z36VbAjrp69t2nkAsqDkh6jrfXbk15/q4FicX5nbsTX99eV88ryzYy++rPNLvtIF55SXFCW4h0EFddhV9xBQ/98gFwp2ztEn4/5bc8c8ZX8fz8pt/fQpm0ERwPjDGz9wnaCIygaujIrEYmsheSVQMB1GzbxcTZS3lt+caEnkHR3jwAww/sxfa63VRV11JVXZtw3LOOLAGDB19awYMv7f3MKjvqnGH7dWfp+m0pf+1Xra7l4J//q9H2wjxIV8jo2cWo2emxaqV+PbrEuqD2LCpQG0EH4mecwd+/diXH3DaBJR+txPbfn6dGf5Xv7Hs846Yuzvr8U5kkgtOb3kWk7TScTiLahXLxBzWxX/mT5yxnwYqNAAwf1CuWKMpKilmwclPC8aL19EUFeTw4dwX5Lfz/906KEsOhA3rw1prGDdFRTdU01ex0Pr5fd44Z2pcH566kIM8YWzkk0nU0Z3NCSiswM1adPIqqytMY/7lDsMJCPgeMy9G6DunmGurh7pvTdRWN7pOd0ESa1nA6iZ5FBcxYvCb2y3jS7GWUlRQnvGfBiuDCP7ZyCL8YdSgf+8VTCa9Hf7lvrwuuxK1QbZ9UuiSQqbfXbo1VTZ1++EDGjy5XQ3EHFZ0DK3rRN3K37Gm6EsHjZraAYG6hV919C4CZHQR8BjgfuB/4R9ajFEmh4ZxBUeUlxUz53khG3zmrUXVPVM+iQs6+a3Yuwkwr34JkE71PJV0Jok/3QnoWFSoJdHBtta5DujWLTyEYO/BtYJGZbTKz9cCfgRJgjLsrCUibS9YDaHF1UO+eKgkAPPTSMhZX1yb0BmoL0Yt/U0ngmKF9U76+YesuarfXJUypHU8Nx5JOU2sWTwemp9tHpKVausxhsjmDUinKN4b270GXgjwWrtxE14I8dtQ1r+9/Ln1qSG8mv7ScspJiKgb34qG5KxvtE73Yd6YlMCU3Mpl0bmYm25Ls84CZrTWzN+O2XWdmq8xsQeQ2au9Dls6kJYvDxO8/cdZSykuCKR/KG7QJlA0s5sgDegKwfbdTtbqWowf1ZuxxQ+iWSXeJduD1VTWUlxRTVV3L3+d/kHSfBSs3drolMCU30jUWFwHdgX5m1oeg7QKgJ7B/BseeBNwJPNhg+63u/tu9D1U6m5YsDhNlZhR3LYhN9nbWnbNYXF1Ln+6FbIjMDRTtJtq7WwEbt9VRVJjHpDkdZ7qsfj26sHBVDRcdOwggYSqMqD7dCznhY/3Iy8tr1joLEm7pSgTfBl4lmF5ifuTxqwSNx3c1dWB3fx74qBVilE4qWrcfnRBu6E+nx5LA3ly0fnjaMKZ8byTjRpbGLpIbtu7iqEgpIGrjtjr6dC9sNAK4pd1DmyuT85aVFLNu807KS4rp3b0rpxy6X8Lr0RLQhq272Lxzd8pR00oCkk66xuLb3H0o8GN3Hxp3O8rd72zBOb8XWdPggUhJQ0KspRetW2e8w7l3z+aGaVVcc2ZZwmsrN25rtH+0lBAvW91Dk+nbLZ8+3QvTnjff4J0JpzF8UG+qIt1gTzl0P6449RBmvrU2Yd8bp78VJMHK0lh/83TrLIgkk0kN6SYzu6jhRndvWOWTiXuAGwjWPr4B+B2QdN4iM7sEuARg8ODByXaRTiDVRSvZIjLRi1z84u6btu5kwYqNLFixkamvVyccZ/2WXa0yiVtr2q9XN95avTk2KjiZLvl53Dj9LY46oCcLVmykMD+Pmu27Yqui9evRhVGHDyA/Lz+o+nGo96CU03AK7vjqNlDJQJLLJBF8Ku5xEXAKQVXRXicCd18TfWxm9wNT0+x7H3AfBLOP7u25pP1pWOcfHQHc1EUr2gumZ1EBNdvruObMMm6YVkVx12AahaMO7MXClZsSVgfbd59C1m/Z1a6SAMBbqzfTqyifTdt3p9ynW5c8HnxpBWUlxRx1QE8WrtzEwsjo57KSYkaU9mHynOWMrRzCuMpSXluxgQUrNjFuZClA2im4lQQkmSYTgbtfHv/czHoBDzXnZGZW4u7Rn23nAm+m2186j2RdGqPTQKS7aMU3KEcbhOcuWR8bOby4upaxlUNiF8qo6CpgDaeLbivxcaRLAtEE1rd7YdIxENMuPx4zIy8vL+Xqa41GqEY+VyUBSaU5I2m2Aoc0tZOZPQzMAT5uZivN7BvAb8zsjciax58BftiM80sHk65LY/n+PRMWUYletKL93eMblKMNwfH3wwf1Yu6S9SnPHb34jqssZcyxbVfFuKPOKczgf1s0gX2UpC0DgimmgSbbVdpqhKp0TE2WCMxsCkGdPkA+UAb8ran3ufuXk2z+015FJ51Cw2kgmurSmOz5+NHlSReR2bFrd2xlrjHHDubxhR+wscF6wmUlxfToms+mbbubrJbJpnSTyHUtCLrBrkuxnvHYyiEYFmsTcBJLOcnaVUQylUkbQXyf/zpgmbs3HtYokkayi3mmF650I4e7FuTHqlMmvxQsJNOwgbiqujZYOnLzTtrjZbIwD4bt14M3PmhcFVQ2sJgRQ/swafYyxlYOYWzlkIQ2ATUGS2tosrDq7v8B3gZ6AX0JkoHIXmlul8aGI4eBhPsFKzdx1lGJ4xuTNRBHG5LbvrWgsV31sHLj9qSvHTO0L9eedRjjRpbSq1sXrj3rME4ctl+jdpVxI0vVGCzNlknV0DeB8cC/CUYX32FmE9z9gWwHJ51DS7o0xi9E37OogGMO2jeh19AxQ/fltRUbcvjXtL58Sz6+YWzlEHp16xIbLRz9jNQYLK0tk6qhnwBHu/t6ADPbF5gNKBFIRuIv5s3p0hh/4Ws4cvb6KYsSqknOvXtWbL2BqOZMLJcH5GIqunTjHAzjilODfhmZtKOINFcmiWAlEF95WQvs/Zp9Emot/RUb/774+17duiQkmEcvreTM21+INSADzZpdtEuB4Q474i7SRutXLR3SvzuLVm9J2DaushTHmTh7KZjq/SX7Muk+ugqYG5k59FrgJeC/ZvYjM/tRdsOTziR+RHCy583xw9OGNRqbULV6M+NGlvL+zaNis47Gy8+DLmm++XkG2+ucXQ1+qUefHTage7P6XTf08f32Yc3mxCqh8pJiJs5eirFn2UklAcm2TL7P7wH/ZM//g8eBaqA4chPJSEunnE4lvpQQXwUF8MkhjRdz2V0PPbsF8/2UDWz8Fa6PfNNTlSNW1ezMuNroyAOKOXTAPo22H3VAT95eu4V1m3dSVlLMkpvOYFxl6Z7F57sVcO1Zh2n9AMmJTEYWX5+LQKRza40ppxseL9lcRA2roHp22zNFdVR5STE9uhZw5pHFGBabphqCC/frq1KvagbExikcOmAf3lqzJeV+/Xp04Z+XHc/u3bs54vpnYmsgA3xiSF8wwwwevbQyaBA+qxwMirsWaPF5yalMeg0NA34MlMbv7+4nZy8s6Wz2dlBZOk2twBU9VpB86lhcXRs7z/VPLGLSnGWMG1nKL0YdyhfumZNw7BUbknfjbGjMsYN5+JXljbYXFeSxT9d81m/ZxbrNO7n+iUU88Xo12+vqKRtYzLTvH79nfqXKUq4ZXRZbZ1i9f6StZFI19HfgNeAXBD2IojeRvdIa8+QnTFcxpcEKXNt2NWqDiK8quvC+l5i79CPGHjeE4q4F3DCtigUrNzGguCvv/fJ0+sYtZtOUyS8tZ2fcAOWiguBv2F5Xz/otuzh0YI/Yfhu2BnMHTb18ZKwr6LiRpfTs1nixeSUBaQuZ9Bqqc/d7sh6JdHqZTjmdTjSZvLZ8IxNnLw161rCnp83vn3k3oV49WlXk7tRu3xWbyG3a5ccz4qZgxdW++wQNsl8ZMYg7n1sSe+/Xjzkw6drAyVz4qcE4zsIVm8Bo1IX1lZ+fQn5+fsLfoIu+tBcpSwRm1tfM+gJTzOy7ZlYS3RbZLpKxhoPK3r95VGxlsuYsmjJ8UK/E4+NMmr0s6bq80dk6p15+fGzd34N+9mTQUDuwB9O+fwJmxpYGkwG9saqWQwf0aDKWsoHF9OxWyHVnH85jl43k0UsrG+1z4/S3GpVWRNqLdFVDrwLzgDEEVUGz2bNc5bzshyadSapBZc2dGsEazBoUnYsn3S/taDKIF00C8UlqyU1nMG5kKQtWbqKoSz5jjks/a2nV6lpqtgVVSu7O6DteTHi9vKS42QlPJBdSVg1FlqkUaTWtMTVCrGQxO2hsjVYNQZAc0vU+qq+v58zbX0jYdubtLzDt+yekGflcwKYM2w1iayxEuoBOvXxPw3B5SXFsIR2R9iaTXkPnJdm8CXjD3dcmeU0kpZZOjRArWVSWMn/5RwmvzV++gRumVtGzW2Gj/vf19fWxpR6j1UHREcgjbprJyz87BTNLSFLROY0mzVnG8EG9OerAnkx7Y3XCSmhRC1ZuSij1XHNmWaxhGIKEoi6h0l5l0lj8DeA44NnI85MIRhcPi0w+16zVykSa64pTD+G6J95k4coaIGgorvd6Js9ZzsKVm4KG4wYlAzOjIC94fsxB+2JmHHPQvlSt3kx+JAE0TEp5eXmNBqj16lbIQ5GeQFFlJcWcOKyfVgeTDiuTRFAPlEXXGzazAQSL0B8DPE8zl60UaS4zo3f3rnsWa4mrHho+qDfjz0q+2M2cn57CuXfNYtLsZUyavQwIZvgEGvU2imo44V3tjt1s2Lqr0Syqxx60b2w/TQgnHU0miaA0ftF5YC0wzN0/MrPMKk9FWln0Ag0kJILHvluZ9sI7fHBvFjRY33jS7GVJRzcnu7BrYXjpjDJJBC+Y2VSCgWUAXwCeN7N9gI1Zi0wkA3s7LiHT3kbpRi/H76uqH+kMMhlZfBkwCRgOHA08CFzm7lvc/TNZjE0kpb0dl9Cwt1G8hskhYfTy1Aajl7c3LgQrCUhHl8mkcw78I3ITaRf2drGb+N5GDRd+b7jCWWvOiyTSEVhTA1zMrJY9U1B3AQqBLe7eeKL3LKmoqPB58zSGTRpLVa+fat/rpyyKtQmMH13e6HnDYw396fTY8/dvHqUkIB2Kmb3q7hVN7ZdJiSBhwnYz+zwwogWxibSavemhY2aNVjS79qzDYqWFhkmgpfMiiXQUmTQWJ3D3f5rZ1dkIRiTbMunn37D9Ib6bKGjpSOl89nZkcR5QQesv3SqSM02VIva2/UGko8ukRHBW3OM6YClwTlaiEWknNEJYwiSTNoJxuQhEpL3RCGEJiybHEZjZgWb2mJmtNbM1ZvaImR2Yi+BERCT7MhlQNhF4AtgfOACYEtkmIiKdQCaJoL+7T3T3ushtEtA/y3GJiEiOZJII1pnZ18wsP3L7GrA+24GJiEhuZJIILgbOB1YD1cAXI9tERKQTSNtryMzygS+4+9k5ikdERHIsbYnA3XejMQMiIp1aJgPKZpnZncBfgS3Rje4+P2tRiUjouTu2eTPU1MD+++NoLEe2ZJIIKiP3E+K2OXBy64cjIgJ3Pj6fkXfcwPCXZ2JFRXifPjx8wfdZc8JpSZcUlZbJZGSxFp8RkZxxd07+9f+wuNZ5+g//4qoLjuXBmyYx6ldX8lif/viph6hk0MoymXSuK8HylKXx+7v7hFTvERFpLlu+nLJ3F/DovU/xx1equWfhk8AASs7/Bt96/UnMzm/rEDudTLqPPk7QYFxH0EYQvYmItL4VK7CPfYyfn3d0wubTLjgNW768jYLq3DJpIzjQ3U/f2wOb2QPAaGCtux8e2daXoNG5lGAW0/PdfUOqY4hICJWX41VV/O6h5xM2z73rzxzzyU+iSqHWl0mJYLaZHdGMY08CGiaQq4GZ7n4IMDPyXEQkxvv04fnTv8wZP7mYG7ut4v0xBzFp2TRKn53OreWn09TyurL3UpYIzOxNoD6yzzgzWwLsAIxgTfsj0x3Y3Z83s9IGm88BToo8ngw8B/xPM+IWkU7KzJh/8RUcNOAAvvqvidhDv+HEE0/k9lv+DzvgQDUUZ0HKxevNbAMwPNUb3X1ZkwcPEsHUuKqhje7eO/4c7t4nxXsvAS4BGDx48CeXLWvydCLSicQvDJTsuTStNRavfz+Ti322uPt9wH0AFRUVKguKhIwWBsqddIlgPzP7UaoX3f2WZpxvjZmVuHu1mZUAa5txDJFW4e7Y7t2Qnw9m+sUpoZWusTgf6AEUp7g1xxPAmMjjMQRdU0Vy7s/3Ps6SoyvxoiLo1Qu//HJu/sc8bp3xTluHJpJz6UoE1S0ZNGZmDxM0DPczs5XAtcCvgL+Z2TeA5cCXmnt8keby5cv5/NUXc+OxX6HH+Lv5+Yj+vPnVb/Opf13O7Fv+pJKBhE66RNCi/wnu/uUUL53SkuOKtJTdfz/7jPka3U69hD/OWsofX/6AghGXMP+BSzj1YFMSkNBJVzWkC7Z0Tm+/jR13HONHl8c21eUXUFw5Anv33TYMTKRtpEwE7v5RLgMRyZnDD8efe44JUxfHNnWt28mWF2bjZWVtGJhI28hkigmRTsW/9S02H34UhW9r1a/dAAAK80lEQVTv5LvfvJifDO/DO+MuY+b+R7Dg7V2MP0RtBBIuSgQSOlZSwiP/+yDn/OlWyi8dhfXqxbAxY/j7J75Az6JCJQEJnZQji9uTiooKnzdvXluHIZ2MRq5KZ5fpyOJMJp0T6ZQ0clUkoEQgIhJySgQiIiGnRCAiEnJKBCIiIadEICISckoEIiIhp0QgIhJySgQiIiGnRCAiEnJKBCIiIadEICISckoEIiIhp0QgIhJySgQiIiGnRCAiEnJKBCIiIadEICISckoEIiIhp0QgIhJySgQiIiGnRCAiEnJKBCIiIadEICISckoEIiIhp0QgIhJySgQiIiGnRCAiEnJKBCIiIadEICISckoEIiIhp0QgIhJySgQiIiFX0BYnNbOlQC2wG6hz94q2iENERNooEUR8xt3XteH5RUQEVQ2JiIReWyUCB542s1fN7JI2ikFERGi7qqGR7v6Bme0HzDCzt9z9+fgdIgniEoDBgwe3RYwiIqHQJiUCd/8gcr8WeAwYkWSf+9y9wt0r+vfvn+sQRURCI+eJwMz2MbPi6GPgs8CbuY5DREQCbVE1NAB4zMyi5/8/d/9XG8QhIiK0QSJw9yXAUbk+r4hIs7jjzz6LzZkDJSXwpS/hPXoQ+THbKaj7qIhIKjt2sOS4k/lw7CV4TQ1Mm4YPG8Z9tz3CrTPeaevoWo0SgYhICn733dTv3MlxF97ChOMvwv/xDx778g844earqNm2E3dv6xBbRVuOLBYRadfskUc4+ObxXLTzQCbOWsrEWUuxwsNYuHsr48u6dprqIZUIRETSMDPGjy5P2FZcVEDnSAEBlQhERFL54hfxW27hhq0DYpvOXvwf1hTuw4CDD+40yUAlAhGRFPzSS3l3/TYu/ObZ/Pm/j/H+uw9w46zJfOOky5gwrUptBCIinZ117cq0X95H31dmc1H9SuyASnpMnMiI/6ygZ1Fhp2kjsI6Q0SoqKnzevHltHYaIhJS7J1z0Gz5vr8zs1UzWe1HVkIhIExpe9DtCEtgbSgQiIiGnRCAiEnJKBCIiIadEICISckoEIiIhp0QgIhJySgQiIiHXIQaUmdmHwLIsHLofsC4Lx20tiq/l2nuMiq9lFF96Q9y9yUXfO0QiyBYzm5fJqLu2ovharr3HqPhaRvG1DlUNiYiEnBKBiEjIhT0R3NfWATRB8bVce49R8bWM4msFoW4jEBERlQhEREIvFInAzJaa2RtmtsDMGi1sYIHbzey/Zva6mX0ih7F9PBJX9FZjZlc02OckM9sUt8/4LMf0gJmtNbM347b1NbMZZvZu5L5PiveOiezzrpmNyWF8/2tmb0X+/R4zs94p3pv2u5DlGK8zs1Vx/46jUrz3dDN7O/J9vDqH8f01LralZrYgxXuz/hma2SAze9bMqsxskZn9ILK9XXwP08TXrr6HGXP3Tn8DlgL90rw+CngSMOBYYG4bxZkPrCbo+xu//SRgag7j+DTwCeDNuG2/Aa6OPL4a+HWS9/UFlkTu+0Qe98lRfJ8FCiKPf50svky+C1mO8Trgxxl8B94DDgK6AAuB8lzE1+D13wHj2+ozBEqAT0QeFwPvAOXt5XuYJr529T3M9BaKEkEGzgEe9MBLQG8zK2mDOE4B3nP3bAyey5i7Pw981GDzOcDkyOPJwOeTvPVzwAx3/8jdNwAzgNNzEZ+7P+3udZGnLwEHtvZ590aKzzATI4D/uvsSd98J/D+Cz75VpYvPglVXzgcebu3zZsrdq919fuRxLVAFHEA7+R6miq+9fQ8zFZZE4MDTZvaqmV2S5PUDgBVxz1dGtuXahaT+z3ecmS00syfN7LBcBhUxwN2rIfhPAOyXZJ/28jleTFDCS6ap70K2fS9SbfBAimqN9vAZngCscfd3U7ye08/QzEqBo4G5tMPvYYP44rXn72GCsCxeP9LdPzCz/YAZZvZW5BdRVLJ153LancrMugBnAz9N8vJ8guqizZF65X8Ch+Qyvgy1h8/x50Ad8JcUuzT1Xcime4AbCD6TGwiqXy5usE+bf4bAl0lfGsjZZ2hmPYBHgCvcvcYyWyIyZ59hw/jitrfn72EjoSgRuPsHkfu1wGMExe94K4FBcc8PBD7ITXQxZwDz3X1NwxfcvcbdN0ceTwcKzaxfjuNbE60ui9yvTbJPm36OkUbB0cBXPVIR21AG34Wscfc17r7b3euB+1Ocu60/wwLgPOCvqfbJ1WdoZoUEF9m/uPujkc3t5nuYIr52/z1MptMnAjPbx8yKo48JGnPebLDbE8BFFjgW2BQtfuZQyl9hZjYwUm+LmY0g+Hdbn8PYIPiMor0vxgCPJ9nnKeCzZtYnUu3x2ci2rDOz04H/Ac52960p9snku5DNGOPbnc5Nce5XgEPMbGiklHghwWefK6cCb7n7ymQv5uozjHzf/wRUufstcS+1i+9hqvg6wvcwqbZurc72jaD3xcLIbRHw88j27wDfiTw24C6C3hpvABU5jrE7wYW9V9y2+Pi+F4l9IUEDVGWW43kYqAZ2Efy6+gawLzATeDdy3zeybwXwx7j3Xgz8N3Ibl8P4/ktQL7wgcrs3su/+wPR034UcxvhQ5Pv1OsEFraRhjJHnowh6obyXrRiTxRfZPin6vYvbN+efIXA8QXXO63H/pqPay/cwTXzt6nuY6U0ji0VEQq7TVw2JiEh6SgQiIiGnRCAiEnJKBCIiIadEICISckoE0imY2ea92PckM6vMZjxJzvl7M/t0KxznmVQzboo0lxKBhNFJQM4SgZn1BY711plC4CHgu61wHJEYJQLptMzsLDOba2avRX5JD4hMEPYd4IeRueBPMLP+ZvaImb0SuY2MvP+6yORwz5nZEjP7ftyxL4pMHrfQzB4ys2Izez8y7QBm1jMy53wh8EXgX3HvXWpmN5nZHDObZ2afMLOnzOw9M/tOZJ8SM3s+EuObZnZC5O1PEIxCF2k1YZl0TsLpRYJf4m5m3wSucvcrzexeYLO7/xbAzP4PuNXdXzSzwQTTEZRFjnEo8BmCOeffNrN7gGHAzwkmDltnZn3dvdbMngPOJJgU8ELgEXffFUks/2gQ2wp3P87MbiUYzTsSKCIYaXov8BXgKXf/pZnlE4w+x903mFlXM9vX3XM9zYh0UkoE0pkdCPw1MsdPF+D9FPudCpTHzWzZMzoXDDDN3XcAO8xsLTAAOBn4h7uvA3D36Lz+fwSuIkgE44BvRbaXAB82OGd0/qA3gB4ezGlfa2bbLVjV6hXggUiJ4p/uHr9a2FqCKQuUCKRVqGpIOrM7gDvd/Qjg2wS/uJPJA45z9+GR2wGRCzPAjrj9dhP8eDKSTGvs7rOAUjM7Ech39+hEYtuSnDt63PoG56gnWOHqeYJVxFYBD5nZRXH7FEWOKdIqlAikM+tFcCGFPTNWAtQSVPVEPU0wsR8AZja8iePOBM43s30j+/eNe+1BggndJsZtqwI+tjeBm9kQYK27308wy+UnItsNGEiw1KFIq1AikM6iu5mtjLv9iGCN4L+b2QvAurh9pwDnRhuLge8DFZHG38UEjckpufsi4JfAf8xsIRA/TfJfCNbJjZ9SfBpBT6W9cRKwwMxeA74A3BbZ/kngJd+zHKJIi2n2UZFWZGZfBM5x96832P4iMNrdN7bw+LcBT7j7zJYcRySeGotFWomZ3UGw0tyoJC9fCQwGWpQIgDeVBKS1qUQgIhJyaiMQEQk5JQIRkZBTIhARCTklAhGRkFMiEBEJOSUCEZGQ+/9Z12qRveqtqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[:, 0], X[:, 1], marker = \"x\")\n",
    "plt.xlabel('Latency(ms)')\n",
    "plt.ylabel('Throughput(mb/s)')\n",
    "plt.scatter(X[listOfOutliers,0], X[listOfOutliers, 1], facecolors = 'none', edgecolors = 'r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The red circles shows the faulty servers in the network. Congratulations! you've successfully written down the code to work on some original and pretty large dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's work this code out in some larger dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "newDataset = sio.loadmat('D:/MLAssignments/AnomalyDetectionScratch/anomalyDataTest.mat')\n",
    "\n",
    "Xtest = newDataset['X']\n",
    "Xvaltest = newDataset['Xval']\n",
    "yvaltest = newDataset['yval']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000, 11)\n",
      "(100, 11)\n",
      "(100, 1)\n"
     ]
    }
   ],
   "source": [
    "print(Xtest.shape)\n",
    "print(Xvaltest.shape)\n",
    "print(yvaltest.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll repeat the same steps as above but simply on a larger dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning dividing by zero!!\n",
      "\n",
      "Best epsilon and F1 are\n",
      " 1.3772288907613575e-18 0.6153846153846154\n"
     ]
    }
   ],
   "source": [
    "mutest, sigma2test = estimateGaussian(Xtest)\n",
    "\n",
    "ptest = multivariateGaussian(Xtest, mutest, sigma2test)\n",
    "\n",
    "pvaltest = multivariateGaussian(Xvaltest, mutest, sigma2test)\n",
    "\n",
    "\n",
    "F1test, epsilontest = selectThreshHold(yvaltest, pvaltest)\n",
    "print('\\nBest epsilon and F1 are\\n',epsilontest, F1test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we just need to check for the outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      " Outliers are:\n",
      " [9, 20, 21, 30, 39, 56, 62, 63, 69, 70, 77, 79, 86, 103, 130, 147, 154, 166, 175, 176, 198, 209, 212, 218, 222, 227, 229, 233, 244, 262, 266, 271, 276, 284, 285, 288, 289, 290, 297, 303, 307, 308, 320, 324, 338, 341, 342, 344, 350, 351, 353, 365, 369, 371, 378, 398, 407, 420, 421, 424, 429, 438, 452, 455, 456, 462, 478, 497, 518, 527, 530, 539, 541, 551, 574, 583, 587, 602, 613, 614, 628, 648, 674, 678, 682, 685, 700, 702, 705, 713, 721, 741, 750, 757, 758, 787, 831, 834, 836, 839, 846, 870, 885, 887, 890, 901, 911, 930, 939, 940, 943, 951, 952, 970, 975, 992, 996]\n",
      "\n",
      "\n",
      "Number of outliers are:  117\n"
     ]
    }
   ],
   "source": [
    "outliersTest = ptest < epsilontest\n",
    "listOfOl = findIndices(outliersTest)\n",
    "\n",
    "print('\\n\\n Outliers are:\\n',listOfOl)\n",
    "print('\\n\\nNumber of outliers are: ',len(listOfOl))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}