[2] Unsupervised k-NN.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.neighbors import NearestNeighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Unnamed: 0', 'Flow.Duration', 'Total.Fwd.Packets',\n",
       "       'Total.Backward.Packets', 'Total.Length.of.Fwd.Packets',\n",
       "       'Total.Length.of.Bwd.Packets', 'Fwd.Packet.Length.Max',\n",
       "       'Fwd.Packet.Length.Mean', 'Fwd.Packet.Length.Std', 'Flow.Bytes.s',\n",
       "       'Flow.IAT.Mean', 'Flow.IAT.Std', 'Fwd.IAT.Mean', 'Fwd.IAT.Std',\n",
       "       'Fwd.IAT.Max', 'Bwd.IAT.Mean', 'Bwd.IAT.Std', 'Fwd.PSH.Flags',\n",
       "       'Fwd.Packets.s', 'Bwd.Packets.s', 'Max.Packet.Length',\n",
       "       'Packet.Length.Mean', 'Packet.Length.Std', 'FIN.Flag.Count',\n",
       "       'SYN.Flag.Count', 'ACK.Flag.Count', 'URG.Flag.Count', 'ECE.Flag.Count',\n",
       "       'Down.Up.Ratio', 'Subflow.Fwd.Packets', 'Subflow.Fwd.Bytes',\n",
       "       'Subflow.Bwd.Packets', 'Subflow.Bwd.Bytes', 'Init_Win_bytes_forward',\n",
       "       'Init_Win_bytes_backward', 'min_seg_size_forward', 'Active.Mean',\n",
       "       'Active.Std', 'Idle.Mean', 'Idle.Std', 'Protocol:Unknown',\n",
       "       'Protocol.TCP', 'Protocol.UDP', 'ProtocolName.FreqEnc',\n",
       "       'Source.IP.Value', 'Destination.IP.Value', 'Source.Port.Value',\n",
       "       'Destination.Port.Value'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset = pd.read_csv('final_dataset.csv')\n",
    "df = dataset.copy()\n",
    "df.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Remove dummy column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.drop(columns=['Unnamed: 0'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Resampling - create reference dataset and dataset for testing phase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(250000, 47)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "REFERENCE_SIZE = 250_000\n",
    "df_reference = df.sample(n=REFERENCE_SIZE)\n",
    "df_reference.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "reference_dataset_indeces = []\n",
    "\n",
    "for index, _ in df_reference.iterrows():\n",
    "    reference_dataset_indeces.append(index)\n",
    "\n",
    "with open('[k-NN] reference indeces.csv', 'w', newline='') as f:\n",
    "    write = csv.writer(f)\n",
    "    write.writerow(['index'])\n",
    "    write.writerow([idx] for idx in sorted(reference_dataset_indeces))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_testing = df[~df.isin(df_reference)].dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entire dataset has 3577296 samples.\n",
      "Reference dataset has 250000 samples.\n",
      "Testing dataset has 3327296 samples.\n"
     ]
    }
   ],
   "source": [
    "print(f\"Entire dataset has {df.shape[0]} samples.\")\n",
    "print(f\"Reference dataset has {df_reference.shape[0]} samples.\")\n",
    "print(f\"Testing dataset has {df_testing.shape[0]} samples.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. k-NN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 5\n",
    "neighbors = NearestNeighbors(n_neighbors=K)\n",
    "neighbors.fit(df_reference.values)\n",
    "distances, _ = neighbors.kneighbors(df_reference.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances_nei = distances.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x304930dc0>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkJUlEQVR4nO3deXgd1X038O8PO06ahBQS+00pJogk0NZQkriKQ97wUt42C0sCb9O8LU5TkpY3NItDCNlMAbM62JitYBds1hgTb5iAwIu8yXhfZEuWJVmyZVmrbUmWrd3az/vHnStG13eZO3fWc7+f5/Hje+eOZs5svzlzzplzRCkFIiLSyzl+J4CIiJzH4E5EpCEGdyIiDTG4ExFpiMGdiEhDDO5ERBoam2oGEXkZwDcBNCulrojzuwD4LwA3AOgB8AOl1L5Uyx0/frzKyclJO8FERNls7969J5VSE1LNlzK4A3gVwFwACxP8fj2AS41/XwLwnPF/Ujk5OSgsLLSweiIiihKRWivzpSyWUUptBnAqySw3A1ioInYCOE9ELrCWTCIicoMTZe4XAqg3fW8wphERkU88rVAVkdtFpFBECltaWrxcNRFRVnEiuDcCuMj0faIx7SxKqQVKqVylVO6ECSnrA4iIyCYngnsegFsl4ioA7Uqp4w4sl4iIbLLSFHIxgGsBjBeRBgD3A/gAACilngewCpFmkFWINIX8N7cSS0RE1qQM7kqpqSl+VwB+6liKiIgoY3xDlSgLlR1rx766034ng1xk5SUmItLMjc9sBQDUzLrR55SQW5hzJyLSEIM7EZGGGNyJiDTE4E5EpCEGdyIiDTG4E5FWuvsG0dU36HcyfMemkESklcvvzwfAZp7MuRMRaYjBnYhIQwzuREQaYnAn24aHFd7Y24DBoWG/k0JEMRjcybalhfX41fL9eHV7jd9JIaIYDO5k26nu/lH/E1FwMLgTEWmIwZ2ISEMM7kREGmJwJyKK0dbTjxlvl6JvcMjvpNjG4E5EFGP2mkos3FGLt4uO+Z0U2xjciYhiDA+ryP9K+ZwS+xjciYg0xOBORKQhBnciIg0xuBMRaYjBnYhIQwzuREQaYnAPmctnrMH9b5f6nQwiCjgG95Dp7h/C73fU+p0MIgo4BnciogTC+woTgzsRhVxbTz+W7qlzdJkiji7OF6EP7rWt3fjlsv0Y4FBvRFnpzqXF+O2KA6g80el3UgLFUnAXketEpFJEqkRkepzfPyUiBSJSJCIlInKD80mN71fL92PFvgYU1bV5tUpf1J/qQW1rt9/JIAqck119AID+QWbwzMammkFExgCYB+BrABoA7BGRPKVUuWm2ewEsU0o9JyKTAKwCkONCerPW/3qswO8knEWFuFMlIt1ZyblPAVCllKpWSvUDWALg5ph5FICPGZ//FEB4+8mktOlQPknhUtLQxsxFClaC+4UA6k3fG4xpZg8A+J6INCCSa/9ZvAWJyO0iUigihS0tLTaSSxRf5YnOkW5aSW9rSk/gprnb8MbeBr+TEmhOVahOBfCqUmoigBsAvCYiZy1bKbVAKZWrlMqdMGGCQ6umbFdc34ZvPL0Z8zdXu76urr5BNHX0ur4eSuzoyUjdU1VLl88pCTYrwb0RwEWm7xONaWa3AVgGAEqpHQA+BGC8EwkkSqXhdA8AoLSx3fV13fTsVnzpdxtcXw9RpqwE9z0ALhWRS0RkHIBbAOTFzFMH4O8BQET+CpHgHphyl4oTHegdCO9YiBQc1SfZYonCIWVwV0oNApgGIB/AQURaxZSJyEMicpMx2y8B/FBE9gNYDOAHKiC1Ha1dfbju6S2YvqLE76QQEXkmZVNIAFBKrUKkotQ8bYbpczmArzibNGd090Vy7HvrTvucEiIKm2BkUe0J/RuqREQAoBzsCUaH5r0M7kQ+GBwaxoaDTX4nQwsCDSKxCxjciXzwzMYq3Pb7QmyqbPY7KaQpBneyLczlkX5rOBVpvtna1e9zSkhXDO42PLfpCP7jtUK/kxEYfCwmXzBzkZSl1jI02uw1FX4ngUx0eoI41d2PyhOd+PJnPuF3UgJLh8pOLzDnTvrQ4KKfumAnpr6w0+9kkAYY3D3y6OqDuOze1X4nQ0saZdxR2cQBJ8gZLJbxyPz33O/UKtuFKeNeZ1SoErkl9Dn3VOWtTr7YQOSUwlq+MR0GYY4foQ/uUakqWdiig4isC3+80Ca4B8nru2qRM30lx3Qk8pBOraacwODugsfzKwFEBnYgCoOd1a3Imb4SrcZg02HCppHxMbiTbcwo6ePFLZEK/311bf4mhBzD4E4ZY86J/MDMRXIM7hpYXliPtp7s7aMkIOPC2BLelPuHeQlrtA/uIb7uLalq7sSv3yjBz5cU+50U3wkfIYhGaB/co3S97nsHIi1yWjrPrgirOdmNsmPODxr9dnHjyKDURABQ19qDGo4vGyh8Q1Vj1z6+CQBQM+tGR5f78yXFmHDuB/GvV13s6HIpvK6ZUwDA+XON7MuanLuuiuvbfFlvvCcFv2n6cEZkC4N7iDWc7sG9b5X6nQwyzFxZjpzpK/1Ohq+6+gbxX+sPY3BIjxf4wlxnp32xTIiPTUodZ6y9JDU8rHDOOczXOql3YAgnu/ow8fwPj0x7YctRH1PkjzsWF+G9Qy0j3x9bU4GFO2qRM/7DuPnzF3qaFievdR3q6LIm5+7lsQraDaVfk1xUkPx40V5cPbvA72Q4Lt1mpXn7j6H9zMDI9+6+IQDwtOsNDeKwK7ImuLthY0UTupN0McCTTl8FlS2pZwoVnq3xBC2jlg5tgrsfZWP//mohfrOixPsVE5GrdLjVhT64+102VtfK9t5+C3Oll2687P88zG8meyH0wZ3CZcnuOuRMX4m+wSHHl+33jZ7e5+bbwjzO1oQ+uPPmnVqQLoYn1h0CALT3DKSY030PvVOOhTtq/E4GkStCH9yjghTAgibdUajaewbwyLvlGHCxlc2BxnYsK6x3bflWvLztKGa8XZZ0nqaOXsx/74hrRQAsWgi2+94qxV+EdGB7bYI7OefR1Qfx4tajeGf/MdfWcdvvC/GbN9yvjP76U+/h+feO2P77Hy3ai0dXV+BIC/tNCZt5BVWOLKcvpCOqMbjTWQaGIrnJoeFw5CqTVeIdaurCrNUVtpfd2Rtp6soctkUB2k1zjBHRshWDe5apa+3BqW49+35nyVxwhP1YxBbzdvQOuNIIwE2WgruIXCcilSJSJSLTE8zzTyJSLiJlIvIHZ5Npnx85rqBl8swn6jVzCnDNY/69WRm0fUOjZXp4on9/uqff887l4l3rSimUNmbe7fWVD6zF917clfFyvJQyuIvIGADzAFwPYBKAqSIyKWaeSwHcDeArSqnLAdzpfFLtiR5uPwZyiLfKIHQs5dTA3ToG6ukrSnDV7zb4nYy48vYfwz/P3+HKsp2+PB5ZeRBfnLne2YUmkiTxC3fU4pvPbsW2qpMZr2ZPzemMl+ElKzn3KQCqlFLVSql+AEsA3Bwzzw8BzFNKnQYApVSzs8nMXBAfE4eHlSdFJGtKT2C7Ayd3rGhZdxD3rV1L9tTjREev38mI647FRdh19JTfyQiVihMdAIBaH142XFN6HIt313m+3igrwf1CAOY2aw3GNLPLAFwmIttEZKeIXBdvQSJyu4gUikhhS4tufXOk75mNhzH54XVocjmY/GxxEb7r5iMl26GSIZNi0LVlJ/Cd57ZbXkbQnxx/tGgf7n7zgG/rd6rL37EALgVwLYCJADaLyF8rpdrMMymlFgBYAAC5ubkBPzT2WT051x9sAgA0d/Thkx/7kJtJ0lrQL3Ky5iev78PgsMLgsMIHxiTOMKT73ka2spJzbwRwken7RGOaWQOAPKXUgFLqKIBDiAT7rOblSeh25rmjd8CRiik3cYBs8krf4BDuXFKExrYzficlISvBfQ+AS0XkEhEZB+AWAHkx87yFSK4dIjIekWKaaueSqZewtZmeubIcVz6wFt98dqsDaQ/XthPFU1DRgreKj+HBvORvOPspZXBXSg0CmAYgH8BBAMuUUmUi8pCI3GTMlg+gVUTKARQA+LVSqtWtRAeJzpnFaBjOxhGGdHT0ZDdypq9EUZ17rT546w4OS2XuSqlVAFbFTJth+qwA3GX8o5DT+YaVzTYbw+H9sSi2VDW4lFIYGlYYO4bvW6aLe0xDWw634Oa5W31b/x2Li3Dbq3s8W58bpVxe9dPvd0530c5aKKVGtjdoJYYLNlfjs/esRluPnm9Vu0n7AbL9kOr6UCqSO3brQvrV8v1o6vD27UCzPBc7HEvGyQeOa+YU4N2fXe3gEoPp3rdK0drVj8qmTtfX9a8v7cL1V1yA737pU5b/5o29DQCA5s4+nPfhcUnnTXY5pTuIiA4tcrQP7r7mRCyeH4EtBglYLs5MKQWlgHPOcW/nNZzuSVmBvKb0OBZsrsYYF9PhttJjzrWCSra7thw+iS2HT6YV3K1IvucjvwbticQLoS6WGRgaxvF2iy8Aheja6+kfxIEG/5odurmrnLrIZq2uwKf/cxUGXexz3izRDfiuZfuxr67Nt1fT23sGMH1FCXr6nelSIlNBi6GBzTh5INTB/Z4/Hgh0O1O77lhcjG/N3YqOXv9HKwqqV7bXAAAGQ9ItsVue2XgYS/bU4w+7/HvN3Yp+F/pED8KRD0IaEgl1cN9Y4XwXBkv31Lmea051QhTXR3KBfQPBGyQg01jqRk4qyBeY28JS3DA/gwFTYgUhNx6ENKQS6uBu5tS+/u2KA/iWyy1NgvgSU0tnHzotPCn843PbLS+zq28Q3Q71QGlJkpNgysz1ePCd9F84Cd6RCqf2M3wK9Zo2wT2V6pZufO7Bta4t/yev7027O98gxfgvzlyPa+dsSjlfcX2b5WVecX8+Lr8/336ikklz3zV39uGVbTUZrDAEWTWbgrJlAboctJA1wR1wN/ew6sAJ15btldaY7ofTbT7mh3iPxyv2NmB9eVPc+aMv8ljxo0X70Oxjk9IwcuKpNCg3m7DLquCum2MuVSanW574zIbDKD/WYWleL55Wfrl8P/7fwsK4v60pS+8mbGdgk94Bd4djK2loQ870ldhTc8qRG3AYyo/D7kz/EIY9rvwPZXA/1d2fxiAXwct9OpWiY+3BaSm066i1roRinw7s6veoCWQq8eLit//ber2EHVsORwZeKah4f0ycTHrENL+wc2YgGE0qneZnFBgcGsZfzViDB2zU+WQilMF98sPrMPnhdZ6us39w2Padd15BlcOp0Vff4BCWFdYHstLZqvLj1p5i/BK7a833hV8s3R/IgaDtng7RTbvvrdKzimV7+gdx95slrlf2DgxFEr+ssD7FnM4KZXCP8rK/icvuXY3/WLTX2swxJ+Kc/ErnE5Rs9eGNi3hy7SH85o0SrE1QZh5P0G8EzZ29aQdMNwKO1cx9bwZNcBMdCbtHKJ0HklSnwa7q0U+Xr++sw+Ld9XEzX1bXG+RTL9TB3esXWNbFCTjJzoFUJ0ii3906YV7YPLqL/XXlTfjcg2ttlxG7kc6WzkgFZlevteIBc5FCUPsDmTJzA376elH8HxPsQzdbdoXV4abOkTFRzeyWSEXrK+xkDoJ5po0W6uBuRSYBqLtv0FaFWirRNKVKm9MVXU+sG/0E8btVB9F+ZsC1itkg8uuijA6pqD0L11vuI+ttLfprT23GdU9veX9VFtZlvoY2pdFSSgfaB/dMXH5/Pq5wq522Zvx8PA3wk7ErUuU0+weH8eA7ZWjvsVa0c1YmwuUderLLn+alQe+iwWna9AqZbRd4LCdz+Sv2RQZzCHJ5YjxW9oFOzf4SbctbxY14ZVsNegeG8ei3/xoAMCe/AqWN8St6g1qc5ZXoef5WUSOGlcK3J0/0N0EO0Sa4Z+JM/xCe21SFaX+X9WN6Y2g4MvJNuvy6D3gVrJ1aT83Jbnzio8n7Jc9UtFWXuXXXvALn+nZJJpN293bOu0xEb2rRtd65tBgAGNx1Mq+gCnMLqvDxj6R/0e33sWteN5gf+UOWcceOI6nb2judS023ffm1j2/CpAs+ltE6d1afyujvz5LmLqk40YHmjj5cc9kEx5JQ1dzl2LKsih66vsEhV3qt9BvL3IGRZmqJXoxJd0DhlCMxuRA2Sxs7tGhPb94zdy0txotbqhPOGyss3T+b28H/ZkVJ2n+/teqkk8lBc8foMRFSnZ/XPb0Ft768O611pCriq23ttvV3qVi5mS/aWYerZ2/MbEUBxOBuwT/YfOMw0Wk10lrGXnIS8ro9vZnTbc1FgDeLGvHIyoOOLtdpQWxjnypFr+2sHfXdr4FGEon3NGS3WMzqYD7NndYqeWeuLEdn7wBOj7xjE7zjH5W1xTLNHb1Yuqce0/7us34n5awT16vTJZCBKc00ZXdV4Oic7avbjuLav/gfI99tt/8O3mlhO01ON0F9YctRvLDl6Mj37r4htPX0pxzf1Q9Zm3O/Y0kRnlh3KGELgkzYCZoNp3vOKvfzKnBl0i+JW9xK0uCwPmWr5iKH3oEhPPBOOf5p/g4fU+T8jSG988C5lVtd7Y7qVnz+oeRdofjVu2rWBvcz/ZFy9iEXsym/WLrf0nxdvYO4enYB/vOPB0ZNt5Oy4vq2jEdL0tni3d727+GmgaFhvGoMNxjNUJhfuluyp37kPE+Hk6fPo6sPpj3OQbrs5gO8ztR43eRUm+Ae3W1zNx525GSqau7MeBlWHwl7jNf/3zPeoMvkFPg/87Zl9JJIOhd2YO4hgUmIt1aVHk85zyvbj6acx03z37NeIZ6uMIw34CdtgnvU42sPjfqe6vC/sLk6bjOorz652cFU2dfc0et6zifs7LQ4SVeiG64Tub+ak9043JR+ZiJeCVNP/9CoohGvy88TrS9VIE70d03GYCnmLr6T7fN0MjbBK4x0VtZWqEatPJA69+On3TXOtWm2cqE/vtZei5sgVsJ5oac/876Hrn18EwCgZtaNtpfhx/6P172B07npaO+Yy/daK06z22+N0zwdOzgB7XLusTrS6Dp1TanzQ+Xl7T921rT3Ow6LfJCY6X5y8zHainR3wWfvWZ3xOhtO99j+26DUbziZDKsNAn7wanpt3cMkkyeyY21nsHi3//3YaB3clVL4zvMJWg/EOXj76tocT8OiHbUp5/GrtUpA4lJcdiufWiy2Vza7enaBrXVZ4WZzUzdPmyW76/CNp5IXTVaeyLxeyk9u7b//OWtjIN7P0Dq4H7P4AoOT7lpWjO40Wyik8xamU870D6G2NZJjTfWyVaa8fHP2y49uGPnsdV8lbnnk3fKUOcEn1x1K+nu6pr95AJU26gHSsa68Cd96dmvaI5y1Gy8QvVHYkNH6dS9z1zq4ez0gLQC8afSoaEU0xx4dEzOal46mus/h/i7Me+Nbc7c6uuxkvHxz1jyAy5Pr/Htj94SDGYsXtx7F3W8eSD2jIdNy70R/3dE7gNtNA49nevP/xdJiHGhsR3f/4Kh1/uV9q/HMhsMJ/67+dKSbic4AlGsHmdbBPRkn79pv7kucg8ikv5OvzEre30Um22CloyY/mpo5WYqxsSK9wRmaOxMH5HTHRb3K9AQRNomOwbI99SmHP3Ti+A0MqbOfRBw8L55/7wjWWGhGmg4r9TZnBoZsj3pmh6XgLiLXiUiliFSJyPQk8/2jiCgRyXUuicH31PrEj8Txgns6QbO1K8k4sQF6s9TpG0GPjZdvMrGzuhVTZiYOyNP+kGCYvBT8riTvODOAoyfjd8rlpXj7wfzCVaq6CSd346zVFfjRon0OLjFSb9Pa1ZdyOx5+t9zR9SaTMriLyBgA8wBcD2ASgKkiMinOfOcC+DmAXU4n0o6ak93YWNGc8PcgxEUrlW3355V5kJLM/WFXHdp7BrCyJHGOKFHPf2alxyJdKN/7lvViCDtiB6wuO+Z8NxSx0m0e50Sx4vzN1fjfRlNLp8W7oZvboyf9Wxe6lk73Rup0Q4a/eWQ9lsepB6huef+8b+rwrh7QSjv3KQCqlFLVACAiSwDcDCD2FvQwgNkAfu1oCm261qUT2gkJX/QIW/2fKcE1rT344sz1CbtNBoC/nbMp5SKjA2O7XV2ybI/73RBMX1GCPxk3ZuT74TT7LC9pbMfnLzov4e9eZ1CsnJ+FtdZ6mHxlW03KecqOWRsrwW6QHnThJNt+5OzumL/5rHf1W2ZWimUuBGC+EhqMaSNEZDKAi5RSSV+lFJHbRaRQRApbWvwZrLa6JXKBFbnQ7DHbPLNxdCuYZIHdD8ku+YEh9++kS/bUWwpiUU41m3Q6k/B4fiVypq90tOhtWWHqm2v0Rh/rSIszA3t4WUTih4wrVEXkHABPAvhlqnmVUguUUrlKqdwJE+yN4rJwR42tv4vK9G791SffO2ta/an0Kk2DlkEPQhGVVasPHEerCwMsB+FCT1aMGI9bHVzGBvG5ARsEpkSz0c/cYqVYphHARabvE41pUecCuALAJuPx6M8A5InITUqpQjjsvxOMBelVsyg/hgNz2+ZD8Z+iHl1VgX/50sWOrCNRrnRwaBiT7s+3NMxZW08/fvz6PnwuSVGFFduqTuIhH4J5qntouh2+pfP2tRNiD6FTTwi/Xr4f+WXJW+EU1jo8tKBPvCx6tZJz3wPgUhG5RETGAbgFQF70R6VUu1JqvFIqRymVA2AnAFcCOxBpThRPEAeeSGTRzlpsjzNUWnQLvL5o73s7fqVtV98g7lpW7Oq615U3pQzsJQ1tAN4v9mk8ndlwev/yoj91/k4/ISXK0Lh1Jbi13FSBHUj/6dgvh5qCk/lLmXNXSg2KyDQA+QDGAHhZKVUmIg8BKFRK5SVfgrPaPQ58bpi1ugIA8Jd/du6o6dGWBm5U9ACwNQhwOi9l2WFlW/c3tOPKiee5mg6d7K9vc2W5bhSH6SbV+xBeZkEtlbkrpVYppS5TSn1GKTXTmDYjXmBXSl3rVq7dCY+tqbTcXIuc40i79bSbumW+SifZGTjDjozfbE6wn3sHhq3MRknsqzttqUmwE7LuDdX1B5tw39ulficj6/xwoXP3e6tBO9p3TlD884KdSX+PfaT3695kvhE/neQFPV1lMthNKm09A5aaBDtBm+CeTi+CdoonKDPbj7Ta/+OY+hSrPT92hazvkZe2+jtqUpR5DIGn17/fx8twiOq1ErHSmdzcjcFqHWSXPoN1BOwR3IqKDLpMzaQP8rDKZH+FUdCKlV7fFdMzZQax3o1jmbB7b5Od1akzGfWnerAhzWapQaRNzj3bXD27AAMaPIGkE8Be2RaMnK1XfrvigK2OpjJuOeZCBt1OP/tWpXMPtPL04UVg96I76qwM7gHLENmmS3/lVmXX1gIHj3eMtKzykhu9gQZlOMuglCxtqnT/BqJNcNclYKfTosLL7kPdsvXw2e39zxKQ8onvvpC8QtQNr26v8XydYWlTHmZe5Mu0Ce66SKd7hSccHn3HD0usdOBlZLf8znVlVCnsoXeT9MxphdWKU7f7+19qof+ZICwzqLQJ7un0DBeQjGBcbr3ApAPumexi5QW602m+s5KsS2rdaNNa5u2iRhyyWANvd/BlomwQppvoFx5eh9yLz/c7GWnzIgJpE9zfLGoEiqy9Kr+m7ITLqbHPj3Ffg27X0VN4ct0h5Iz/iN9JIRMvuk22wmof8tkmdMF9xtiFmHROrd/JcM2Htp2DL44LfxNHR1UA3wMgJwA1zu/E6O/Py/4EfzuOlapu+vOiKcCkZ1xdhzZl7roISGaIsliygcIpPEKXc39o8Fa/k0BElJGXJufiUy6vgzl3IiINMbgTEWmIwZ2IyGNevGvD4E5EpCEGdyIiDTG4ExF5zIu35BnciYg0xOBORKQhBnciIq+xtQwREdnB4E5E5DEvuvxlcCci8lg6gwvZxeBORKQhBnciIg0xuBMRaYjBnYjIY6xQJSIiWxjciYg0xOBOROSxwPTnLiLXiUiliFSJyPQ4v98lIuUiUiIiG0TkYueTSkREVqUM7iIyBsA8ANcDmARgqohMipmtCECuUupKAG8AeMzphBIRkXVWcu5TAFQppaqVUv0AlgC42TyDUqpAKdVjfN0JYKKzySQi0kdQ+nO/EEC96XuDMS2R2wCsjveDiNwuIoUiUtjS0mI9lURElBZHK1RF5HsAcgHMife7UmqBUipXKZU7YcIEJ1dNRBQaNa3drq/DSnBvBHCR6ftEY9ooIvJVAPcAuEkp1edM8oiI9DMnv9L1dVgJ7nsAXCoil4jIOAC3AMgzzyAiXwAwH5HA3ux8MomI9NF+ZsD1daQM7kqpQQDTAOQDOAhgmVKqTEQeEpGbjNnmAPgogOUiUiwieQkWR0REHhhrZSal1CoAq2KmzTB9/qrD6SIiogzwDVUiIg0xuBMRaYjBnYhIQwzuREQaYnAnItIQgzsRkYYY3ImINMTgTkSkIQZ3IiINMbgTEWmIwZ2ISEMM7kREGmJwJyLSEIM7EZGGGNyJiDTE4E5EpCEGdyIiDTG4ExFpiMGdiEhDDO5ERBpicCci0hCDOxGRhhjciYg0xOBORKQhBnciIg0xuBMRaYjBnYhIQwzuREQaYnAnItIQgzsRkYYY3ImINMTgTkSkIUvBXUSuE5FKEakSkelxfv+giCw1ft8lIjmOp5SIiCxLGdxFZAyAeQCuBzAJwFQRmRQz220ATiulPgvgKQCznU4oERFZZyXnPgVAlVKqWinVD2AJgJtj5rkZwO+Nz28A+HsREeeSSURE6bAS3C8EUG/63mBMizuPUmoQQDuAT8QuSERuF5FCESlsaWmxleC8aV+x9XdEREGx8o6rXV/HWNfXYKKUWgBgAQDk5uYqO8u4cuJ5qJl1o6PpIiLSjZWceyOAi0zfJxrT4s4jImMB/CmAVicSSERE6bMS3PcAuFRELhGRcQBuAZAXM08egO8bn78DYKNSylbOnIiIMpeyWEYpNSgi0wDkAxgD4GWlVJmIPASgUCmVB+AlAK+JSBWAU4jcAIiIyCeWytyVUqsArIqZNsP0uRfA/3U2aUREZBffUCUi0hCDOxGRhhjciYg0xOBORKQh8avFooi0AKi1+efjAZx0MDlhwG3ODtzm7JDJNl+slJqQaibfgnsmRKRQKZXrdzq8xG3ODtzm7ODFNrNYhohIQwzuREQaCmtwX+B3AnzAbc4O3Obs4Po2h7LMnYiIkgtrzp2IiJIIXXBPNZ5r0IlIjYgcEJFiESk0pn1cRNaJyGHj//ON6SIizxjbWiIik03L+b4x/2ER+b5p+t8Yy68y/tbzEbFE5GURaRaRUtM017cx0Tp83OYHRKTRONbFInKD6be7jfRXisg3TNPjnt9Gr6y7jOlLjR5afR2/WEQuEpECESkXkTIR+bkxXdtjnWSbg3eslVKh+YdIr5RHAHwawDgA+wFM8jtdaW5DDYDxMdMeAzDd+DwdwGzj8w0AVgMQAFcB2GVM/ziAauP/843P5xu/7TbmFeNvr/dhG68BMBlAqZfbmGgdPm7zAwB+FWfeSca5+0EAlxjn9Jhk5zeAZQBuMT4/D+DHxuefAHje+HwLgKUebvMFACYbn88FcMjYNm2PdZJtDtyx9vSid2DHfhlAvun73QDu9jtdaW5DDc4O7pUALjCdPJXG5/kApsbOB2AqgPmm6fONaRcAqDBNHzWfx9uZg9GBzvVtTLQOH7c50QU/6rxFpDvtLyc6v43AdhLAWGP6yHzRvzU+jzXmE5+O+dsAvpYNxzrONgfuWIetWMbKeK5BpwCsFZG9InK7Me2TSqnjxucTAD5pfE60vcmmN8SZHgRebGOidfhpmlEE8bKp6CDdbf4EgDYVGZ/YPH3UslSS8YvdZhQRfAHALmTJsY7ZZiBgxzpswV0HVyulJgO4HsBPReQa848qclvWugmTF9sYkP34HIDPAPg8gOMAnvA1NS4RkY8CWAHgTqVUh/k3XY91nG0O3LEOW3C3Mp5roCmlGo3/mwH8EcAUAE0icgEAGP83G7Mn2t5k0yfGmR4EXmxjonX4QinVpJQaUkoNA3gBkWMNpL/NrQDOk8j4xObpo5YlPoxfLCIfQCTIva6UetOYrPWxjrfNQTzWYQvuVsZzDSwR+YiInBv9DODrAEoxegza7yNSjgdj+q1GK4OrALQbj6L5AL4uIucbj39fR6Rc7jiADhG5ymhVcKtpWX7zYhsTrcMX0eBj+AdEjjUQSectRuuHSwBcikjFYdzz28iZFiAyPjFw9v7zZfxiY/+/BOCgUupJ00/aHutE2xzIY+1HJUSGFRg3IFJDfQTAPX6nJ820fxqRWvH9AMqi6Uek3GwDgMMA1gP4uDFdAMwztvUAgFzTsv4dQJXx799M03ONE+sIgLnwoXINwGJEHk0HECkzvM2LbUy0Dh+3+TVjm0qMC/MC0/z3GOmvhKlFU6Lz2zh3dhv7YjmADxrTP2R8rzJ+/7SH23w1IsUhJQCKjX836Hysk2xz4I4131AlItJQ2IpliIjIAgZ3IiINMbgTEWmIwZ2ISEMM7kREGmJwJyLSEIM7EZGGGNyJiDT0/wGX+PAmbJ01MQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Including self-sample distances - zeros\n",
    "plt.plot(distances_nei.mean(axis=1))\n",
    "plt.plot([0,250_000], [np.percentile(distances_nei, 99.5), np.percentile(distances_nei, 99.5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x30497f8b0>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXjElEQVR4nO3de3BcZ33G8e9vd3W1ZUm25AuWg+zGCZiQkkQ4pg3XJiRO22RoudgzDAFC0ilNhxbaGTMwKU07nVIKMzAEiGlTJkAJKS1gBjOmTUMvTByshDiJ7TpRHBtLcWzZlm+yLnv59Y9z1l7JkrWSd3V2j57PzM6e855Xe95XZ/346D03c3dERCReElE3QERESk/hLiISQwp3EZEYUriLiMSQwl1EJIYU7iIiMZSaqoKZPQj8DnDE3a+aYLkBXwRuBc4CH3T3p6b63La2Nu/s7Jx2g0VE5rInn3zyqLu3T1VvynAHvgF8GXhokuXrgdXh63rgq+H7RXV2dtLd3V3E6kVEJM/MDhRTb8phGXf/b+D4RarcDjzkge1Ai5ktK66ZIiJSDqUYc18OHCyY7w3LREQkIrN6QNXM7jazbjPr7u/vn81Vi4jMKaUI9z5gRcF8R1h2AXff7O5d7t7V3j7l8QAREZmhUoT7FuADFlgHnHT3QyX4XBERmaFiToX8DvA2oM3MeoG/AGoA3P1rwFaC0yB7CE6F/FC5GisiIsWZMtzdfeMUyx34o5K1SERELpmuUBURKZMn9h3jCz/dSyabm/V1K9xFRMqk+8AAX/rPHrIRPBRJ4S4iEkMKdxGRGFK4i4jEkMJdRCSGFO4iImVm2KyvU+EuIhJDCncRkRhSuIuIxJDCXUSkTDyCi5fyFO4iImVms388VeEuIhJHCncRkRhSuIuIxJDCXUSkTCI8nqpwFxEptwiOpyrcRUTiSOEuIhJDCncRkRhSuIuIlEmEx1MV7iIi5WYRXKKqcBcRiSGFu4hIDCncRURiSOEuIlImukJVRCTGdIWqiIiUhMJdRCSGFO4iIjGkcBcRKROP8BpVhbuISJlV7DNUzewWM9trZj1mtmmC5ZeZ2WNm9ksze8bMbi19U0VEpFhThruZJYH7gfXAGmCjma0ZV+3TwCPufg2wAfhKqRsqIiLFK2bPfS3Q4+773H0UeBi4fVwdBxaE083Ay6VrooiITFcx4b4cOFgw3xuWFfoM8H4z6wW2An880QeZ2d1m1m1m3f39/TNorohI9YjDFaobgW+4ewdwK/BNM7vgs919s7t3uXtXe3t7iVYtIlLZKvWWv33AioL5jrCs0J3AIwDu/jhQD7SVooEiIjJ9xYT7DmC1ma00s1qCA6ZbxtX5FfBbAGb2WoJw17iLiEhEpgx3d88A9wDbgD0EZ8XsMrP7zOy2sNongLvMbCfwHeCD7lGONomIzG2pYiq5+1aCA6WFZfcWTO8GfrO0TRMRqW56hqqIiJSUwl1EJIYU7iIiMaRwFxGJIYW7iEi5RHjSoMJdRKSMorjdLyjcRURiSeEuIhJDCncRkRhSuIuIlImuUBURiamIjqcq3EVE4kjhLiISQwp3EZEYUriLiJRJHJ6hKiIiE4ji+amgcBcRiSWFu4hIDCncRURiSOEuIlImHuE1qgp3EZEy0hWqIiJSMgp3EZEYUriLiMSQwl1EpEx0haqISEzpGaoiIlIyCncRkRhSuIuIxJDCXUSkTPQMVRGRmLKIrlFVuIuIxFBR4W5mt5jZXjPrMbNNk9R5r5ntNrNdZvbPpW2miIhMR2qqCmaWBO4HbgJ6gR1mtsXddxfUWQ18EvhNdx8ws8XlarCIiEytmD33tUCPu+9z91HgYeD2cXXuAu539wEAdz9S2maKiMh0FBPuy4GDBfO9YVmhK4ArzOznZrbdzG6Z6IPM7G4z6zaz7v7+/pm1WESkSsTh9gMpYDXwNmAj8HUzaxlfyd03u3uXu3e1t7eXaNUiIpXJ3UlEdNpKMavtA1YUzHeEZYV6gS3unnb3l4DnCcJeRGTOyrmTiOjmMsWE+w5gtZmtNLNaYAOwZVydHxDstWNmbQTDNPtK10wRkeqTcyo33N09A9wDbAP2AI+4+y4zu8/MbgurbQOOmdlu4DHgz939WLkaLSJSDXLukd0VcspTIQHcfSuwdVzZvQXTDnw8fImICMEB1YrdcxcRkZkJxtyjWbfCXUSkTCr9gKqIiMxAzsEU7iIi8eIalhERiZ9cTgdURURix9Geu4hI7GjMXUQkhqK8iEnhLiJSJrqISUQkhnQRk4hIDFX0jcNERGRmNOYuIhJDrtsPiIjEjy5iEhGJIQ3LiIjEkKM9dxGR2Kn0B2SLiMgM6FRIEZEYCsbcFe4iIrES7LlHs26Fu4hImbg7EWW7wl1EpFz0DFURkRjSRUwiIjGki5hERGJI93MXEYkhRxcxiYjEji5iEhGJIV3EJCISQ7qISUQkhvSwDhGRGNIDskVEYiiXo7LH3M3sFjPba2Y9ZrbpIvV+38zczLpK10QRkepU0XvuZpYE7gfWA2uAjWa2ZoJ6TcDHgCdK3UgRkWrkDhbRrcOK2XNfC/S4+z53HwUeBm6foN5fAZ8FhkvYPhGRqpWr8CcxLQcOFsz3hmXnmNm1wAp3//HFPsjM7jazbjPr7u/vn3ZjRUSqSVWf525mCeALwCemquvum929y9272tvbL3XVIiIVrdIfkN0HrCiY7wjL8pqAq4Cfmdl+YB2wRQdVRWSuy+Uq+IAqsANYbWYrzawW2ABsyS9095Pu3ubune7eCWwHbnP37rK0WESkSoxkctSlohl0n3Kt7p4B7gG2AXuAR9x9l5ndZ2a3lbuBIiLVKgj3ZCTrThVTyd23AlvHld07Sd23XXqzRESq33A6S31Nhe65i4jI9Ll7GO7R7Lkr3EVEyiCddXJO5Y65i4jI9I1ksgDacxcRiZPhdA6AOoW7iEh8DKeDPXcNy4iIxIiGZUREYig/LFOvPXcRkfjI77lrzF1EJEa05y4iEkP5A6pRjbkXdfuBivKTTfDKs1G3QkTkot4wOMLDtWdYvbUZasdF7dLXw/q/Lev6tecuIlIGOQ/eo7qfe/XtuZf5fzsRkVLYtv0An/7Bc/xiw29Rv6B+1tevPXcRkTI4PZwBoKm+JpL1K9xFRMrg9HCaVMJ0y18RkTg5M5Jhfn2qeh+QLSIiFzo1lKapPrrDmgp3EZEyOH42zcLG2sjWr3AXESmD44MjLJyncBcRiZWBwTQL59VFtn6Fu4hIGRwbHGHRfO25i4jExtnRDMPpHK0acxcRiY8jp0YAaG/SsIyISGwcPjUMwNIIbjuQp3AXESmxV8JwX7JAe+4iIrHRd2IIgFe1NETWBoW7iEiJvXxiiJbGGubV6QpVEZHYePnEMK9qjm6vHRTuIiIl1ztwNtIhGVC4i4iUVCab46Wjg1y+eH6k7VC4i4iU0MGBIdJZV7iLiMTJvv4zAKxsa4y0HUWFu5ndYmZ7zazHzDZNsPzjZrbbzJ4xs0fN7NWlb6qISOXbc+gUAKuXNEXajinD3cySwP3AemANsNHM1oyr9kugy92vBr4H/F2pGyoiUg2ePDDA5YvnsyCiZ6fmFbPnvhbocfd97j4KPAzcXljB3R9z97Ph7Hago7TNFBGpfNmc8+SBAd7Y2Rp1U4oK9+XAwYL53rBsMncCP5logZndbWbdZtbd399ffCtFRKrArpdPcmo4w/UrF0XdlNIeUDWz9wNdwOcmWu7um929y9272tvbS7lqEZHI/dfeYKf1htVtEbcEirk2tg9YUTDfEZaNYWY3Ap8C3uruI6VpnohI9di2+xWuuayFtvnR3TAsr5g99x3AajNbaWa1wAZgS2EFM7sGeAC4zd2PlL6ZIiKVrXfgLM/1neLm1y2NuilAEeHu7hngHmAbsAd4xN13mdl9ZnZbWO1zwHzgX8zsaTPbMsnHiYjE0o92HgJg/VWVEe5F3bLM3bcCW8eV3VswfWOJ2yUiUjWyOefbTxxgbedCXr1oXtTNAXSFqojIJfuPPYfpHRjiwzd0Rt2UcxTuIiKX6J9+/hLLmuu58bVLom7KOQp3EZFL8POeo2zfd5w7b1hJKlk5kVo5LRERqTKjmRx/+aNddLQ28P51lXVLreieASUiUuW+9OgLPH/4DP/wgS7qa5JRN2cM7bmLiMzAU78a4Cs/6+Hd13Vw45rKGWvPU7iLiEzTyyeGuOfbT/Gqlgbu/d3xN8mtDAp3EZFp6B04y/s2P87p4Qxfe/91kd/adzIacxcRKdLB42fZ+PXtnBxK862PXM9Vy5ujbtKkFO4iIkU4ePwsGzZv5/Rwmm9/5Hqu7miJukkXpXAXEZnCU78a4KPfeoqhdJZ/vmtdRe+x52nMXURkErmc8+D/vsT7HnicmpTx8N3VEeygPXcRkQk923uSv/7xbp546Tg3vnYxn3/PG2hurMyDpxNRuIuIFDh4/Cx//9O9/PDpl1k4r5a/edfr2bh2BWYWddOmReEuIgIMDI7y5cd6+ObjB0gk4J63X84fvHUVTRV6quNUFO4iMqcNjmR46PEDfOVnPQyOZHjPdSv405uuYGlzfdRNuyQKdxGZk/pODPHQ4/t5+BcHOTmU5u1XtrNp/Wu5cmlT1E0rCYW7iMwZ6WyOR/cc4bs7fsV/Pd8PwM2vW8pdb1nFtZe1Rty60lK4i0isuTu7D53i+0/18f1f9nFscJQlC+r46NsuZ8PaFXS0NkbdxLJQuItI7Lg7PUfO8ONnD/GjnS/zYv8gNUnjHa9ZzHu7VvDWK9or6sEa5aBwF5FYyGRzdB8Y4NE9h/n33YfZf+wsZrBu5SI+fMNKbr1qGa3zaqNu5qxRuItI1Tpyepj/ef4oj+09wn8/38+p4Qw1SeNNv9bGnW9exc1rlrB4QXWf9TJTCncRqRpHz4zQvX+A7fuO8fOeo7xw5AwAbfPruPl1S3nHaxbz5ivamV+naNNvQEQq0mgmx95XTrOz9wTP9J6g+8AA+/oHAaivSfDGzoX83rUdvHl1G2uWLSCRqK4rSMtN4S4ikcvmnH39Z9jZe5Jnek+ws/ckew6dYjSTA6ClsYZrL2vlvV0reGNnK1ctb6YuVVnPLK00CncRmVXD6Swv9p/h+cOn2XPoNDsPnuC5vpMMjmYBmFeb5KrlzXzwNzq5uqOZX+9ooaO1oeru7RI1hbuIlJy7c/TMKPuPDfJS/yA9/WfY1z/Ivv4z7D82SM6DerWpBGuWLeDd13Xw+o4Wfr2jmVXt80lqiOWSKdxFZNry4d13Yoi+gSFePjFE34khegeC94PHz3JmJHOufm0qQeeiRq5Y0sRvX72MK5c2ccWSJla2zaMm5uebR0XhLiIXGE5n6T89ci60+8LQLpweCcfD8+bXpVje0sDy1gauX7mQVy9qpLNtHqva5tHR2qi98VmmcBeJOXdncDTLibOjnBxKc3IozamhNCfOBtPHB0fpPz3C4dPDHDk1wpHTI5wcSl/wOYvm1bK8tYErlzbxjtcsZnlrw7kw72hpZEFDSuPiFUThLlIF3J2hdPZcOOeDeXxQn1selueXZ/KD3BOoTSZob6pj8YI6VrXPY92qRSxZUMfipnqWNNcHAd7SQEOtzk6pJgp3kTJKZ3MMp7MMp4P3kUyWodEcZ0czDI5mGBzJMjiS4cxIhlPDmTCoR8cEdf6Vzk4e0AmD5oaa86/GWi5b2EhzQ4rmhhpaGmppbqhhQbi8pfF83cbapPa4Y6iocDezW4AvAkngH9z9b8ctrwMeAq4DjgHvc/f9pW2qyKVxd0azOYbTOUbygZvJjgnfoXQwPzLBssL6I/myTFA+NBpMnytPZxnO5MheZI95PDNoqkvR3Hg+jJc1N7BgXBgHYX0+qJsba5hfm9JFPDLGlOFuZkngfuAmoBfYYWZb3H13QbU7gQF3v9zMNgCfBd5XjgbLpXN3cn7+PeeOh++5CZYVLs/mglcm/54Ny9zJ5nLn5jNj6uUuqH9BeW78ZwfLMlknnc2RPvcerGM0myMzvjznjGZyjIbzo5nglc6eL5+pVMKor0lSX5OgLhW8B/PBdEtDDfU1Sery5WGdhoI6dfnpVIJ5dSkaa5PMr0sxL3zNr0vpoKOUTDF77muBHnffB2BmDwO3A4XhfjvwmXD6e8CXzczcvfjdliI9suMgm/9nHxAEEMCYlTgXlE1Uz8/V8wvLJmn1tD9ngnpMWW/ydUy2fLI+Txbc1SCVMJIJI5UwalIJapIJasLpVMKC+WSCmmQw3ViboiZppJIJalMJ6sLltanz77VJo64mSV1qbDAHQVwY2PkAP1+m0/Wk2hQT7suBgwXzvcD1k9Vx94yZnQQWAUcLK5nZ3cDdAJdddtmMGtw6r5YrlxQ8BsvGvOXXM0FZcfXOf54V1Ltg8bgyu7DMxv/EZG2YfH0TtXVMU6doQzJhmEHCjET4buG0EZZdpE6+LJgP6uUDN3hPkExAMpEYE8bJhJFK2sTliQTJpI0pT4xbnl+niMzcrB5QdffNwGaArq6uGe1D3rRmCTetWVLSdomIxE0xf2v2ASsK5jvCsgnrmFkKaCY4sCoiIhEoJtx3AKvNbKWZ1QIbgC3j6mwB7gin3w38ZznG20VEpDhTDsuEY+j3ANsIToV80N13mdl9QLe7bwH+EfimmfUAxwn+AxARkYgUNebu7luBrePK7i2YHgbeU9qmiYjITOn8LhGRGFK4i4jEkMJdRCSGFO4iIjFkUZ2xaGb9wIEZ/ngb465+nQPU57lBfZ4bLqXPr3b39qkqRRbul8LMut29K+p2zCb1eW5Qn+eG2eizhmVERGJI4S4iEkPVGu6bo25ABNTnuUF9nhvK3ueqHHMXEZGLq9Y9dxERuYiqC3czu8XM9ppZj5ltiro902Vm+83sWTN72sy6w7KFZvbvZvZC+N4alpuZfSns6zNmdm3B59wR1n/BzO4oKL8u/Pye8Gdn/akXZvagmR0xs+cKysrex8nWEWGfP2NmfeG2ftrMbi1Y9smw/XvN7OaC8gm/3+FdWZ8Iy78b3qEVM6sL53vC5Z2z1GXMbIWZPWZmu81sl5l9LCyP7ba+SJ8rb1u7e9W8CO5K+SKwCqgFdgJrom7XNPuwH2gbV/Z3wKZwehPw2XD6VuAnBA9lWgc8EZYvBPaF763hdGu47BdhXQt/dn0EfXwLcC3w3Gz2cbJ1RNjnzwB/NkHdNeF3tw5YGX6nkxf7fgOPABvC6a8BfxhOfxT4Wji9AfjuLPZ5GXBtON0EPB/2Lbbb+iJ9rrhtPav/6Evwi30TsK1g/pPAJ6Nu1zT7sJ8Lw30vsKzgy7M3nH4A2Di+HrAReKCg/IGwbBnwfwXlY+rNcj87GRt0Ze/jZOuIsM+T/YMf870luJ32myb7fofBdhRIheXn6uV/NpxOhfUsom3+Q+CmubCtJ+hzxW3rahuWmeh5rssjastMOfBTM3vSgmfKAixx90Ph9CtA/jmCk/X3YuW9E5RXgtno42TriNI94RDEgwVDB9Pt8yLghLtnxpWP+axwef75xbMqHCK4BniCObKtx/UZKmxbV1u4x8EN7n4tsB74IzN7S+FCD/5bjvUpTLPRxwr5PX4V+DXgDcAh4PORtqZMzGw+8K/An7j7qcJlcd3WE/S54rZ1tYV7Mc9zrWju3he+HwG+D6wFDpvZMoDw/UhYfbL+Xqy8Y4LySjAbfZxsHZFw98PunnX3HPB1gm0N0+/zMaDFgucTF5aP+SyL4PnFZlZDEHLfdvd/C4tjva0n6nMlbutqC/dinudascxsnpk15aeBdwLPMfYZtHcQjOMRln8gPMtgHXAy/FN0G/BOM2sN//x7J8G43CHglJmtC88q+EDBZ0VtNvo42ToikQ+f0LsItjUE7dwQnv2wElhNcOBwwu93uGf6GMHzieHC318kzy8Of///COxx9y8ULIrttp6szxW5raM4CHGJBzBuJThC/SLwqajbM822ryI4Kr4T2JVvP8G42aPAC8B/AAvDcgPuD/v6LNBV8FkfBnrC14cKyrvCL9aLwJeJ4OAa8B2CP03TBGOGd85GHydbR4R9/mbYp2fCf5jLCup/Kmz/XgrOaJrs+x1+d34R/i7+BagLy+vD+Z5w+apZ7PMNBMMhzwBPh69b47ytL9LnitvWukJVRCSGqm1YRkREiqBwFxGJIYW7iEgMKdxFRGJI4S4iEkMKdxGRGFK4i4jEkMJdRCSG/h/MJN6zRfGTcAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "distances_mean_non_zeros = distances_nei.mean(axis=1)[distances_nei.mean(axis=1) != 0]\n",
    "plt.plot(np.sort(distances_mean_non_zeros))\n",
    "plt.plot([0,250_000], [np.percentile(distances_nei, 99.5), np.percentile(distances_nei, 99.5)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Distribution of distances' value is that the vast majority is slightly close to zero. Ploting sorted data indicates that anomalous threshold can be taken as the 99.5 centile (empirical adjustment). Finally, each newly arrived sample in the network traffic with the distance higher than threshold, could be interpretted as anomalous sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Separation of anomalous samples from test dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "anomalous_distance_threshold = np.percentile(distances_nei, 99.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "anomalous_df = []\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "for index, row in df_testing.iterrows():\n",
    "    distance, _ = neighbors.kneighbors([np.array(row.tolist())])\n",
    "    \n",
    "    if distance.mean() > anomalous_distance_threshold:\n",
    "        anomalous_df.append(index)\n",
    "\n",
    "stop = time.time() - start\n",
    "avg_time = stop/df_testing.shape[0]\n",
    "print(f\"Avg processing time of single sample {avg_time:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discovered 21310 samples to be suspected as anomalous. <br>\n",
    "Avg processing time of single sample 0.00251"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('[k-NN] anonmalous indeces.csv', 'w', newline='') as f:\n",
    "    write = csv.writer(f) \n",
    "    write.writerow(['index'])\n",
    "    write.writerows([element] for element in anomalous_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anomalous_full_df = df.iloc[anomalous_df,:]\n",
    "anomalous_full_df.to_csv('anomalous_dataset.csv', index=False, header=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.5 64-bit ('3.10.5')",
   "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.10.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "16ce62d89139d6b459a4fa414ee333754219537208081da9003b7fab6a6c11db"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}