Clustering along lines

linear_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1147bdf50>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO19f6xb53nec0hdirqsSEeM43uRTuEtUjhhFrRYDVMq2m4FUlEpGqQR2q3B3VDATh1gyRAbKIj8QKtqRYQtxao/lnWg1xoLhq1GhzpY2k1S4wBNg/HGjjLElu3GrRqndtI0bmxLcmwrku1vf3x8eZ7z8juH5L2Xl+Th+wAHvPfwnMNDnu97vvd73h9f5JyDwWAwGPKJwqxvwGAwGAzTg5G8wWAw5BhG8gaDwZBjGMkbDAZDjmEkbzAYDDnGvlnfAOONb3yjazQas74Ng8FgWCh89atf/Z5z7ubQe3NF8o1GAxcuXJj1bRgMBsNCIYqiv017z+Qag8FgyDGM5A0GgyHHMJI3GAyGHMNI3mAwGHIMI3mDwWDIMYzkDQaDIccwkjcYDIYcw0jeYDAYcozlIfmtLeD4cf9qMBgMS4K5ynidKk6dAs6f93+fOzfbezEYDIY9wvKQ/MmTyVeDwWBYAiwPyR89aha8wWBYOiyPJm8wGAxTwk5cftN2Fy6PJW8wGAxTwk5cftN2F+bbkt+N4fXeey0qx2AwZOLkSaDd3p7LbyfnjoN8k7wMkadOZR8XGgzk3I9/fLxrGAyGpQW7/Ca1CeXco0enc2/5lmvGjagJzZfknBMngAcesKgcg8EAwBP4qVOeEjQxz2WktnNubraf+ImfcDNBr+dcu+1fp3G8wWBYWOju3m47B/jXUcfuFQBccCm8mh+5Zif6+zjzJb7+uDKQwWBYeOjuPm0NfdeRxv6z2HZkyWcNr+Miaxjm65slbzDkEt2uc/W6c51O3MUn6e67QUPbATIs+ZkTO287IvndIN7dmIfxcTYYGAwLhXrdU8DKin+t1yfrvvMo18yc2HmbmSYvmOQJyZDf7Sb380Axq2HdYDCMBd3l2ZIXwl+EybuR/E6Q9nSlBdTr/v9Oxw//m5vx8WkDgcFgmAuIHdZqDXdz7vrzbq8Zye8ErVbcChiawGV+t7ISHzPvLcNgWHIIgUs3111V3u9251uFNZLfLno952q1MMlrbG46F0X+lc+ft9ZgMOQMu9HNej3fxVutJKGH7LR5tN2M5LcLeZorK8OSi25ZaRZ/1jlp+wwGw9jYrbh17u5ZWvw8xlYYyW8XvV5Se+enqud3aSTP872QJydtnmgwGMbCuJHPWed2u841m84Vi8nuPgrzYtUbye8EvZ5zlYr/qZrNdE9Nt+ulnWYz2TrkeD1YZHl8DAbDriAt9kHbWNI9Aeeq1bg76gFk1P+zgpH8TiG6fK2W/lSFtGUwkJYV8tw4Nz+tw2DIKXgiri1tIfdm0//dbA7bciFNfl5tMyN5xjjkmhY8mxUKyZ4bGRRWVhYjk8JgWGCMsru00ipdlcma/2YFVdto3W6yi89aphFMneQB3AfgWQCP0b5DAD4P4K/7r28YdZ09IflxRDTdOhjjWObdbtJ7Ezo/dF7o3uZ1fmgwzAlC3ZVtrlBhMelGnY634CuV2IYLdfFQjuOk2bDTxF6Q/M8A+CeK5D8F4KP9vz8K4N+Pus5cWfJpoZPyhKvV5PuaoEeZF9rqT3Pla4fuvHh6DIY5Qa83ujvKcToqRrpfqEuFBoV5zXHcE7kGQEOR/JMA1vt/rwN4ctQ15kqTH5UEtbGRfD/LwyOinwThMrlrk0BfR+5D/AFm2RsMzrmkxS3dSfR0PdHWEAIvFJwrl/15ocm5ng3wufNkZ82K5C/T3xH/r867C8AFABcOHz489R9jbIxjifP7oSfPsg2bCzwb0NkXXCFJSF32hbw989jiDIY9gA5cq9WcW10dts1CFjx3TalMMm7i0zzaVTMn+f7/L4y6xp5a8tt9UmnnhSxsTfCFQmwycHodhwCE9PysmPp5bHEGw4TYTjPW2rl0D03yPBnmLqQnzYucq2hyTQg7tYDTWkSz6ed/hUJskW9uxlkW2qWv4+flOizxmDRjyDkm6Y5p5Cy208aGd6SKPcXK6LyFPu4WZkXyv6Mcr58adY2FsOQFIc2eTQmWXJrN5H5pgVqq4fvia2W1fCN+Qw4wSTMWJyswXBK42012NXmvUvFkzy6z7URTzyv2IrrmDwF8B8ANAN8CcCeAOoAv9EMoHwRwaNR15srxmgWOvmGHqLS+KPLWfKPhybpcjvdLxoV4ifiaTO4i52ivj4Zp8oacYFxClZgH3qQbcuYqJztJtxI3mA6LTItnWJSqI5YMtduQ1iHaebMZ/y2tiC1xTqXj/dVqfE1OwQtlX6TJNotiahgMKQgRaigJSScyhUheW/JiL3W7sYIaInZN+FmxDvMII/ndhm6BnOGqLfFmMynNNBreoi+V/N/y3vq6P69USg4S7LwVc0LkH54JGAwLilDAmraj5JUnuNINxN0l3YM1eMBb/mzhF4vZcRPzmOw0Ckby0wavGZaVZid/a/KWyBttfrD8wyl50oorlfEkHYNhzsAx6KGYdg4yk65RLsdNXxyv2soXFxh3JwmrrFR8l5JFutNi6RdxcmwkvxsY58kzoetl30V2WVvzpsTmZlKf37dv2NTQ81IxVeTa7NCdd9HQYCCMiivgsEdxaQlxR9FwN9O5g41GbOGzLq8t9bx0HSP53cA4Dk4eCPQasGzF87xTC4iAHwhE5imVXELGEXNFS0HbCTBeJFPFsFAY1cR0Irc+hxVJ+ZtlFyF2nhA3m7FTVuwnqf4tk2LW9bfTdeYVRvK7gUmJMRTMy7Hy0sJk/lkoxCZK1iaO3TQzaJIZxyKJjoaFwiibKBQ7wI5XJvm0OANtN/HGSiYrpNId8xaUlkXyBRiysbUFHD/u/z53Djh6dLzz7roL+N73gHe+Mz7/1lv9a6XiXx96CHjzm4F6HXj9daBajc8vlcLXvXoV+PrXgbU1oFgEfuzH/D0eOQK84x3Au98NnD8P3HOP/9x77/WvW1vxNU6e9J/53HPAqVPj/xYGw5g4eRJot/1rCEePxt3pnnt8k33xxficl1/2x8krAFy86I+9fNn/f+KE7wIAUCgAURQf+9JLwE03AWfOAK2W71o3bgAPPOC7wuXLfn/a/eUKaew/i20uLfmdDvlsNbOniKtcitW/uZn0MokAKftqtXjeyYU3tElTrw8vebMdi99g2AOEgsVEbpFKINwlAB+MxhPfSiV+v9EIJ6NzZE6erHjnsi35mRM7b3NJ8pOSYSirQhOttLJKxYuIIuNw7Je0YF66RgYEKZ1QLnv9Xq7TaMRpfTrM08jcMCMwwYaaoSb5Xi8ZkyDRMWkbkz9LMvrz2a7KW5cwkt9tZBHoOGXrOP2OSV2n64mwyJZ7oxFu6dXqeOECZsEb9hg80eSae2K7SGzBxkYycVyseF0VRCz5NNLnKBqt3bOjN08wkt9tZMVgpZGohAzUanE8vcw5sxyu/H69PhxPD/iZgPSEUc5UPQhZ9qxhBxi3/kta+ohuxjrgrFKJo47f8pb4OG3B8yZJ49LMeeDQy0PkBUbyu43tSCHcqrWOXq0mE59KpWEyLxbjjFnAuUOH/IAhWbMSP1atJuWaUZkeoWStPImVhl1HWqhj6JishCOx5HVMO0s1WRvbRqJesiXPA0ve7Rcj+XkAW/Ldrt8OHPBkLpY96/Kt1rA0w/NP6V2s40tPYcdrmukS6oV57wmGXQHbAkLI5XLYdkhbRZObmsQhiEspjchlO3QoGY0sjlZd5mmZbBYj+b3ApASphUo5X2qnbm7617U1N5i38rWFxLWcI/XrZWbABc9CPSCPXijDVMEpIGKHiBUtC3OI7cC2CDdBadYixehJq5ZieJ80bR2vwJPRcYu45gVG8nuBSc2GtJADXZkp5DGSc7WkI//zEoIhi5+JPlRQbVl6hmGASWwUburSfDY2kk1WXEM6+XscGUY3x7SCY6EaNKHEqmWAkfxeYKdSB8snsryNnrNKi5c5Mte7kW1lJbmSgljqLP3w3DlUB2eZeofBOZdto4Sigllq4ebDzbNQSBZa7XSSdsjKinPHjiUDw8RSZ/dStxt3gSiKP4snwHkoMrYTGMnPEuO2No7Y4dryIt3UaknLXUfaFItx7yoWkxUxOfnqwIG49/Dnplnyy9ZblhRZjzlrANDLILAzljdtR5TL4Yra1WrScg/l+nHBsSyLfZmarpH8LDGOjJOWrRFaobhYjFeckgzZAwfiGQBb9FqqqdeTlZvksyW8odOJ93ElJ+l9y9BbDEMIkaX2//OKljxpXF9PLqKmyZ8no5rg9++PA8XY0mdy12vsjEpZySuM5GeJSQqGaSeo7K9WkxUnuYXLe0LyHI/GcWnazOLcb16ZQWfo6nLGy2QeLSnGecRpWjw3O2kynU4szYTSPLRNwlY83xO7sHTo5jg5iHmGkfy8I030FMFTWq8eDDhjVtdg5f+1w5aDkZvNcPqg5IdrS55lpWXoPTnEqBj2kASjm6iuJyOWvDQrtux1uSW9FQp+sNDx8Z1OPEB0Okmil88X99IyEXoIRvKLAj0H1vnZuiV3Ot5s2tiIewNb3Fpj18HL0htDPU8+n2P75R7Tgp8Nc4tQmn9aFQwheSFwziCVEEmRZPbvj5tCt5ucPMr1ZXIpyyQ0Gr7ZieoY0vDlvniSydq8Rf4mYSS/KGB5Jo1EOUhZjl9dTSZChXLHeWYgPVRWVNBzbMD3SB4MmAmy1phddpNqTqHT/Dklgwt6sW0gES2SosHNQf5mBTFkL6ysJAPFeAEPad5M8pIQJXXk2ZLvdmN10ppXEkbyiwI9jw7Np7UMwx4pDp8Uj1UoAUpMq0ZjOMhYpxiWSumWfMgZu0zergVCaOwNKW+hyZ6QqpCzBGfpCBcpNKY3th+4oqTU3xOS5+akE8T53qxpDcNIflGhW7W0fA6B7HZ9z5GFvvUgwGGUIprqUEzuWaLPS9iDniFIj5Rr8L1lCb2GmSItQkZLNloxlKArOS6K4iAs7SLSljwTul5vVefmaZcU2y7igLVJYjqM5BcVWo8Piai6p7ElzyGVOn6tUIh7oRzPxzQacUgEa/1yTKUSjgRitlimvPI5R9rjSVsmmC13eeQyyZMgLG1ti4uoXk9a7zIwsAwTstTlc6WZVyrD7xvCMJLPA9KiWnTKIQ8ATPIshkoPYmctm06hOvYSbbO56Xvt5mby/tIGpGnPr828y0Sa/s7uGm07SJPiRbElY1X79mVCJ2O+lmd0/T22R0KJTBa8tT0YyecBaWSmLXmtsesa9jyHlt5dq8UVo+R/7tE81+aqUCF0u/F5Eooxzd5qQm0m0iJptJLHTUuOZX++RNjwYMEDhVjr1WpSm6/VkvKPNAm5Zsjxa2P25DCSzzN0r9CkzymJOnxBlwoMect0YRLu1aG69VwcTQaaUbLNTnq2sUIm0lwlWY9FW/LsguFmsbmZbDIcB18oxLHy4lSVzxL7Qa6/nfHZHnsSRvJ5RZo3jefCvBgJk32arMMyi1SYkoGh0xmuAauvxyGe+tppvdms8R1jFOnJY5HKjaHHwgMCh0Tq6NpQ8VMOuJLzhNT1LEL+l5WgJiFq3bytyXjMlOQBfBPARQBfy7oRZyQ/ObIWBgllzYolr9mAy/zpAiK8Vau+Tk5oP5cNDN2jLpUcyvA1syyIcX6etHFSHr8Qr1jboTw3eVQhTZ0ne7r4qUTian+/NEuOdZd7Ym1/EuiJqjUZj3kg+TeOc6yR/ITQvSmNDXSWqp7Ds5OV4+vFak9bg1aWKuTYORFd5XMkaofDKsZZcNwwwDgTnbRJnZBppZK05DltQit7etMhjzyBk/d5LNfr03BUr0BP+Mb5Pln7lx1G8nlFmh4fYgM+VnvjODlKiJh7fbOZXllKi7BcglDXl+UBhS17Qya2S2zymHnsbre9ll4sOnfwoH+V1AiuPSPVIyW8MkT8zaZ3zkoYZugxi2NXE3ooAzeUqGV2wHiYNck/BeD/AfgqgLsC798F4AKAC4cPH572b5FvjMsG2pLXmSmhQGXxkmmrfnU1uf6aZgTxCayuhitL7fS7GJxzYQWMSxexaqcfIYdByhjNdWWY/BuNeEbAKqBIPzJhy1IS9f2GCN0e/2SYNcm/uf/6JgCPAPiZtGPNkt8GJukNacdKFov0bBZvtUUu2bUi3rKJKBkvEkNXKMRr1TK58z2kzduXNGB6p1a7ECVPmJjYdbRsoZD0pesEKB77dSikVuf4uEkSn43Qd465ia4B8FsAfj3tfSP5bWDUvHeUycTXYHNMgp55kRHp2cViHG/PsfChRTzTFkLRIZcs2Mp9h8ol5xzjyhRplruuty6Jy0L0PNayslYsJpclaDbj8r+yhF+zGQ8I7L/nz7aol9lgZiQPoALgIP3dA3A87Xgj+W1gFImzRZxW+EzO5RWkuN6shFu0WslEJz1fD1Wokvm9XuInLeRSf7clc9JmjdP82EJjOxM4JxxxDhsnIOn67vITy7XFqtchlPJeKDXDKlnMBrMk+R/pSzSPAHgcwCeyjjeS3yHSGEITZZq5KLFusqYs92rx3vExLPrq6lSFQjrzSJ48O3nT5vejwjByiLTgJ36MobGdF+zQse6FQtKJKuoaPy79SNiZKkQvbhqeYOnKGhbeuPeYG7lm1GYkPyVkxaWH4unlVYp3cwFwfYz0aF5dKoriSB0Wh7Wg224PR/qkSUlLYsk7l4yKCSldoYmYDoDin5ZdK+yQbTbjtAd55WhcHsN1BeusBctMstl7GMkb0qFJVss63NM5W0YWFRFhl5cRkuNkCSAO35DMWT3IiA9ArxaxuZl7szA0Bme5I7RU02rFlnYUxeWBOT1CfOR8Pk/S5DFKjRrtSNXEHbITLLdtdjCSXzZMGnEjMXJ6yR1mA7medq7K+m+ix5fL2UsK8kDA+oBehUpYS2SfHDNGKNxQO1IZoexSIXLxh+/fH1eKrlT8/xsbcTilxLhzCQKdGsFJUKFiYob5gZH8smFSiSMrwkUzjWSwCqnrmjgs8sprqNgJawjMVKLXi/wjOkLO5v48DmuSl59YCDtNg2+1ktWjhcB5XyhZmSdXYrnrWnVC+voR5ewx5AZG8suGSefLbBqOOld6PlvebFYKq5TLyXm+xOHzYCBMI9a/xO0JO3FMHw80OdAD0iJfe71kNItUb9SuDa2eSZgj75ey/1xrRkfKyMDCjzXkE7eomfmGkbxhfIxKPxTTU4KotcTDKZNyLjtpeaFQXSyl3R5e2ER7+dLCKheM+NNulxUyXrtlYyNpleufglW0UikuK6Qdspub/rpScoiX7h3l+16idIWFg5G8YXyE2EebncwoOl5OFy0plz3jrK/HBC5sIf4AseS7Xb8xW3F9es6lZwdxVqLXDDDueJMW8SpqVbHoNfV6PZmCsLoaV4Ysl/2mJ0m8ceAT/5xplru+714vGWBlmD8YyRt2hpAgLOyga9HwsaEMWIm3F7Ae32z6fXJNMV05Okc+Oy0jaFxTc0qW/yQTjbRxqdcb/tm0WyP005bLw8eVSsmxlSNlORlq1Ng4qhaNYbYwkjfsHDrAWsw/Tq3ULMZ16lkQZqJnkhcdvtFIshVn+QjTsFN4kkIp8l12MZM2NAZqaUMTeloc+uZmknyLRb/dfrv/SUSC6XZHFwbln4zvZ5TWnjbDSPuJF0wpyyWM5A07B4vFodTLtPALDo9sNJKVsoTlJAwzjaWk8JkUUNMzCglDGdfM3GWRedSEIhSkpNMTuL6MkDvLLzxWihtEZJiDB2NVDPA/p/y9sZEuwaQR86TK1xwpZUsLI3nDzjHKcmbzMJToVKslq1cyg4kDt9n0zCUmqzAaaxMS9H3gwLAIzbOJPTA9Q185BD2mhCZFbMlLKd+0BbpkAZDQezJYyLgXWgiM7ylEzHOifBkmgJG8YXvYTu8NBVczI+m4eC0KC7NJaIlk0pbLw6xXrSZjAiWyRwd5T4l9xrVgtTqko1Cl1H5I3uF0A121OaTbiwQjP1toSV+5JyPm/CCL5AswGNJw6hRw/rx/HRcnTwL1OvDcc/7/dhs4fRpotYBmE3jrW4E77vDHCIpF4MQJ//eb3+xfb74Z+Mu/BK5di49529uAAwfi897zHuDWW/17APDCC+H7uOce4Phx4N57gSNH/HbvvX7f1tZkv4n6qu22v/XQpba2/H4A+JM/iY+9fNn/HPJVX34ZeOgh4Od+Lr6GXFu+7uuvA+fO+et85CNAtQqsryc/721vA+66C7jppvhnu/lmf52TJ4fv6dw54OjRbX99w6Igjf1nsZklP2fYrrk3KghcNA6WcMTU1CUNxURlp6uYxZxzL7IQ595rfwEL27vgeB3ljghZ+lq/D6UK8LVlqT6xzHkmID/NyspwzDsv6jFORI9hsQGTawx7Dk302vvIufmhOHsRrGu1OAtW6uVK8RXOruVVMOSYajVOBZVjo2i4rOI2cWez586i7d7f6I10bmr9XpcRlq+hs1Yl/FHcEKurPnaeI2tEh2fHL4dKLnDOmGFMGMkb9h7aZAzFEIa8llwhS87hLFi5RppXUo7Xnkkh+VIp2/xOA7Gj/Pmlir+/Xq09FBKZNokRFwGPcbwgl9wuj1lZG6/amJY6MGmEqWHxYCRv2HtoMzat+IlmP866YZZik7ZSGV7WiMm800lmzgq5syeTZxAhkzvlPh+utwe3eGez5x6ut92j3V5CQpExSitHo5J1G42YtHu99K/Ikxoeq+TnazTi9df5JzaJJr8wkjdMF6M0gBDLpOkKaRqHLF7Cq06nZQNJmQSuzMXx+DIQiI4vbNpsxn/rqJyeJ/Qj6AXL7uqwSAn6qdWSiUdy+/v3x9Gocp6cI8VA5dbk5+GoUa4eKeMnFyyTCQ3fmxUYyy+M5A3TxShTcVROf1oVTM1aQMxyImBrKSe0iVeS9wlziowj3kslC32203MPrvjXR7ue6O9s9gaTAV0Wn5UgTgdIs8rlNuQr8uIe8hNwrpf+DPZDy+REXBhSJWKcR2RYbBjJG6aL7XjzQrUANHvpkgdcT1cE7ErFr1LFzLmykp1QFUXee6mZVa5JGUQPrvRlmmJcQfNGperubvUSt8dRL/JxcjmOdZfjDx3yVr5Ez7AWr4lcvhLPHHQQUqhswbhlCQyLDyN5w3wjZL2LqMwMyKI2eym1Nc8ROpoJxWHLprV4LgNeU7Hk/34jmcD1cN1b972at+w/2+m5L1Xa7l2V3tDa5ULusloTTyDYf1wuJ8c8cciKysRVG7JKBhmWD0byhsVBr5esa9vpxN5I1s85xp5DJEulmCmFbUslz7SNhvvusU13A0V3rb7ubpQrfcs8ZtrvFyruarM1cKb2en0nabXn/r7hB4LHqy13BD33+aJn488X24NIm7NoJyYbosH3P35IruFJSK0W/wRsyctAECJ5qwppcC6b5C3j1TBfOHoUePvb4/+/+EXgtdeAW24BHnkEuHHDZ7heuQK85S0+nfMjHwFWV/3xb31rfO4rr/jX69d9yug3v4kDD/4p9uE17H/uO7h2DXgBVfzxzR8cnF95/SUcfOIhPPvhUzh/3ifVdjrAlavA088AV3EQv/vDZ3CxchS/8dpJnEMbv/HaSZzCSTxYbOOPcQK/99Rx7PvKFp54Is48/aEf8l8BAKJo+BaLReBTn4r3X74MvPiiz4z94Ad98u4dd8QZq3fc4b/6mTP+eNm/gwReQ16Rxv6z2MySX3JwNI1oHax9aAs9be06MXPF0pdlkBoN9/3I77uOWAa6XkkGqV9H0X0AXVcoOHcEPfeFUts9hvj9LbQGgTvFYqytV6suYdFrd4JIP0fQGyyMpSs3c6ilKEwcJpnmQDXH6nIDJtcYFgLMVNqLKSGPoVh2rW+ImC1aiRoMbqDont7XcFdB0TTFont5fcO9DC+AX0TT/fS+nnseXiv5PuLQzS20EpILE/PPln0W7BH0Bu8fgd93seLvr1drJyJH+ZUlGpZpRoXxm2N1uWEkb1gMhGLkdSVLNlV16KUOGNdhJ91uwuF6A0X38npj4Nx9sVhzN+C9pFdQcV9GfB0h+RsouE82uu4jt8dkLmNQvR6uAX8WfnB5tNJyZ9F271uLB4D9+5OpAuMQusGgYSRvWFykFX3hKly6gFmf8K80W+7uVi8myb4l/1p/E2J/ZtM7dP+h2hgw9KVycyDRXMGqu4T4vYfrbderJWUZUYrETyzRMsWiGzhp39/oJYKCeCZgJQgMO4GRvCE/YEknxJj9GPdLne5A/x4Y/4FaAc9gzb2w4s3op6IN9wMU3VPRhnO9nnt6n88qegWlgbRzA0X355td966Kt+R/Muq5cjmWZM40u+7ppt8v+36q2BvINp8vJqUcPTkxbd2wHcyU5AEcB/AkgEsAPpp1rJG8YSS0pNMPtxTLXP6/WvSvDxda7pONrrte68s6UuK47zl9Jepr7cXiQKr5/r6qc+22u1GiEgrkrJXB4yza7gPourNou62+tPNa0Q8iz6I+2Pd0s+2q1Vi2eb5YH1j1eqEtK0Fg2A5mRvIAigD+BsCPACgBeARAM+14I3nDxOj1BoR+tVh1rtNx12t1962St8IvVlruWZBnU9CXd26UygNL3gHuGlbc4xFl0xLJX0RzQOwv9B2yQvyXUfVROrTvIpruzwpt98ymv6f7Gx33csV/1lm03d0tP2Dd3eolrHez5g2TIovkpx0nfzuAS865bzjnrgO4H8B7p/yZhmXC0aPArf24+lvfDjzyCFauPIdnrr8JX6m34X73DP578zRu1OrAnXcOgslfetNbAADfxTq+fuOteCrawAuo4sP4NH5v5SO4UavjW7/4IbyICgDgRVTwa/h9vGffOdxRewA34QoAYAWv4TpWUMNVRACeKDRxHSsAgEIEHHv9HCr3/wFWrjyHf/HiH+Dj//hPcA5tnF89gU99/T3A+fP4ty/eg4frx/E7J3yQu6wKJas5GQw7Qhr778YG4JcA/D79/68AfDrteLPklwC7vdqUfq/Xc8+32u7+RieWaARkIt/Z7LkttBVBzUQAACAASURBVNyVqD8L6GvuV1AZyCxfK7fc8/DvP4/qQJr5ZKPrXoUPnbnRP1eib14qxSUWXkDVHUHP/eaaD595ZrMzKIPwfMvfy/Va3T1enbC2vcGggBnKNSNJHsBdAC4AuHD48OFp/xaGWWMCLSLB6ynnpXJ/PSDRUAhLr9YeRMy8hLK7jrh62JfRcl+qxDr7NcQ6uwPclypt92J/UHhd6fV/VWq6LbTc5f7gcBbtQZj+w/W4Jr3ETd7f6LiLaLqrxbgAvIVOGibFLEn+KIDz9P/HAHws7Xiz5JcAYzIYJxi128nzxilFn4ih16CkKCboGyi4G4C7gch9BptuCy13EU135VDDuVrNfffYpnuxWHOXyk33GWwOnL03UHDfQMNdhCf4ew913LUDNfdXpaZ7V6XnPr3pSxT/z2Nd93Cx5Wvg9GP5r9covbU/gE2kyduIYHCzJfl9AL4BYAOx4/UdaccbyRsEesGotPfkNW1NU+cCfN+XdO491HHPou5Oo+POoj2QZhwwkGOeL8Yk/Oq+0uDv5wrx/hv9V8mO/UHf8r9arLkttAZO2gdX4tmBiyLv2O10hsJpJuJt89Ia3AxJ3n82fh7AX8FH2Xwi61gjeYNgHAme86OE4FdW4goI8r7EtLer8SyAk2Jl+zV03TWU3GuA+78Hj7knai33d/Wme7V/AEszlzrdRP0bB7i/qzfduajtvvJ2XyT+VQrBlPj6J2rqg/VagLv5QxmWBjMl+Uk2I3nDpNAFvYA4EVb2n4+8tfu19fZgX6PhE2XLZb9+yF1R130vqruL/SzXB1faJOt4vf5llNyzqLvfXOu6M82uu4qKu9G3+F2l4h6vtnx4JYVkDtXDl8VJKpW4to7csFnjhm3CSN6QW4hawYmsvP5pve6GMk/b1Z77s0JcDbLddgOn6vOouXORXxDE9fxCIKfhZZ0PoOuOwO+7dsAT87VyzT3fag8iZAb6fmm/J2+RY3jkEcudV6IiPSlonKcVZTML3uCM5A0LjFFcxtJNo+EN504nriGztua5VFYCLJXc0AIfvZ5z/6bcdc+i7v71SjeRbSrXlwU/JGv1+4VVdxUV9w003ONR011CY+DI9Zp+IbbOe73kiiEynZCCN0L2tZpzvZ67u9WLk6UEWns3Ld5AMJI3LBTSln8dBTGWa7WkZd9uJ9dXZcteClhubDj3k1HPnY/iejeycPej3d6gHtoR9AbOWK3Js24/kGqkYLzsE6eB3KB2ELTb7mrT/3+1Scs+mSVvyICRvGGhwMSuuSwUGakN5WYzuWpgtxvL3zIIyAqDwrFH0BtINuciX+9dKk1+qdIezBbubvXc1WbLXS3QkoFRxb281hiMJN8vVBJLCrrV1fjvRiO+2XI5vHKIre1nmBBG8oa5wrgSTOh9ttaFD1WJmaHQeBk0KhVP7p2O59RqNV6ESmSYZ1F3R9BzrZZzdza9xf/h/V23hZbX3fsE/GilH0NfqMYae/8GLqLpnqrQurMHDiSJnqcVPN0Y9QPs5Icz5BpG8oa5wrgSTKgioygfWgXRfNnpeC49cMAvzVevx8fr5fW6XU/ovZqv+Q64QeVKKR0sJzxVabrnW17CGZQjoJW7rzZbw2GSHEFTq3ltSPZXq9klJ0PaFa8swsfoZRENSwMjecPMEDIux42B12Tc6/n9En0ohrKWwTud5D6Rv3lhj0olPr+lpO+7W71ByeDBB/dJegutweD08lrDOSCOh+d1aWVrNOIoGxl9ajW/nxd2Tas1HNKueNFX5+JjRIMykl86GMkbZoZJ5WUdElmpxOqGGLCa1NfW/HEi2+i4+c1NHzZ5Fm33T0u9oYGB+fTRbi+ZWSXaTz9L9u5Wz13q+BHoleiAE03+4XrbXWm2khfl0Em5ebnJcjkZTinHsPYUSpTq9eI1bOXeZETkdXANSwUjecOuY1z5d1KS5+QmXtiayVhyiRqNJCdWKvEiHK1WLIM3m25ouT6RxkVpkc+RImKJegr6y/Zv8LWo4F4s1twnG10H+BnAgIDX1pJkLR+g9Xghch7N9Pts7esfJGTZG5YORvKGXUdIVw+tU9rtehVB4tGzIEpFs+lla3aSsuGqY9e1VS5gXj0CHx750/u8JV8qxQTPfPloNzB66S/b7caZrM3mwMJPVMvkTFeJ0+QZAs8UeHoS0qB4FBIZiKcuWqPXP6pZ97mHkbxh1xHijlDhMBUCngktxYh1zoXHut2YI4VHpTyBHM/hlmpJ1yHZp9lML1YZ+rKi2csqUEMlCeTYzc2YsPm9VsuPTuVyXKCs0fDHlkrDoUKa5LVmJY6GUCQOy0Vm5ecaRvKGPUHIaSr8NM6apaGQyP5SrK5QCBciq9c9n7LhLHzGlvz6eqyGiHN2JAeSDi/33m7H4ZYDa5qnKVpWkYQovgATeGgU4n2FgnOHDvnXtbWkHqUtfHbo8ugaKuWZ9vDM4l9IGMkbpo60hMyQ31D/nxYmGTJom03n9u/33Lq5GZO0bBy0wteRwUGM63o9XM44cX99ouSIGrHkn2+RE7VYHCbxVivpAWapp1r1o47oUZypJdMKHoH0JkkCfMNpP4ReKTwNViZhoWEkb5g6QpF+OntfH+dc0uAUx6mQM1vnpVLMhSH/pCZznhFUKjFvNhpxElSW5NRuxzf3RK0V5kc9iowKh+QPYB2drX6ZYohXudFIFt5J0766XU/0UVwVM7M2hJVJyBWM5A1TRyhnh7XzUEi4c2GrvVaLQ8xDfsjQPsAPAlqGET51Lmm5pxm3ifsbJ8M0TfMZlbUqX7xU8l+Ipxghq53j77O0L9G05Pqsn7HnOm2FFcNCwkjesKdgS575RnhKImUox8hFUeyLDHGd3sSalyq9bOBybpCuKMlO4XHl6swvKSFE2nublp0qCK1aIj8MT2GiKA41YlIOzRJ4Pw8G2pIP6VSGhYaRvGFqGGeWr41Xtu61n5GTRjlE8uBBf9yxY56b5D3hRTmv08kO2Qw5h7dlzDJxpskhzWZM2PpDhIw5nl6mHBI3HxrdNjaSxM0DiU4TlllAyJI3aSZXMJI3TAVsietkJ+EVcY6Gas2IpS1lBiSiUHiKDVoxSsWgDXGjc8nom0wpJvD/xF+etfeQlznrZkbdGJO1hBjJ1mz6zxKNfn09JnWtZ4mDN/SQ9D2MGwZlmDsYyRumAm1MMkKBIVTHKyFjs3Ih+7ViIQ5Tngnw35zhr529+n6nJkOHkqZGBuFnXCfNw1woJH+0QsH/cJwsEDovRPKhTFrT6RcORvKGiaGN05CjMsv4E6udjdD19Zjz2HiVYyWYRPyRIY6TsHOx/mV/vT5876Ms+an9aNv9AP2jc3niEGGLtc8/lBRDK5f9/mPH4mmS/PChAmg6ocEknYWCkbwhE7rfs1qwQsUYs4y8tEEhFD0j6oXW6jl8kiUdDkfne2Djs9PZu99rJLY7yujwH/1jiAwjU5eQd1onSPExPIWSHzIkO2nyN7KfexjJGzLBsgv7EZngxbmZFsXHES3MQzKAsNSiZRqJkJEqvGLVy+eITq+LLIZqdc0FF6XpQqPi1dnJEdLXdVXLbjf5o4q3WRdDC02J5AfmB8qFzkKhmGbdzy2M5A3OuWwDk6PuxArX5KoNS5F4RXKRRE4toegiixK5xzlAPFvQdbvSZhE8A5mrQoyTWPJp2WHNZjLsSK8XK2FFjYb/YRuNeKBoNpPEzteS66TVyMmy5NMGLyP/mcNI3uCcS9a4YmjfWyicmo1L9vfpIorCHYVCPDjweaur/nOk3owEiAhncSCIzgHKygNaWJ7JCvmRLFYmZF6gNkunl03WkZXFS/QDkwesf9yQwyXNkTx1j7ZhFIzkDc65dJJPk181yXKGvSZtydcRxykXS+RVmCQARGetyj3pe9QDTej+cw1+CCzVaILnxAN2ZPCAwIlUnISwshIm6rQkqlBDWcgRNj8wkjc450bLNaFaViGS5e3AgZhjpP/zwKBDIiX6b3PTc8vmZlIqEolI5+1oKz83GNcZy5JJt5ssahZFXifj0VgWD9/cjB9SuRz/sPxAJXtMN4C0WcZcaWMG54zklw6jeCMUUMFSL3OKlkuYWxqNpE4fRT5ij6Xk1VV/TS4HzKqBTt7kAmOhkM3cGYyjpI4sHX9UmWLWz2Qrlfz5m5vD5zKBC+FLJttuPYhcPsTZw0h+yaB5Q/crfl/n3TQayeAKTbzCGXxcqzUcwCH5OSyxsMqg682zIcnROOMYiwvNG+PcfJpFzaU1JS4e8LWYWZMPTcE4Bn/fvuHIHK3vFwpxrP0495+230Izp4KZkDyA3wLwbQBf628/P+ocI/ndQRap6/c5Jh7wfZ3PGVXeXPqqSLwcQcO1tqS8r54Z8Ap2Av7scTgg936/UaO2Hj35f14mUCx5flBar+dFweUBsQOlVkvem0zlpOGk3bPAJJ+pYJYk/+uTnGMkPx2EAiX0e2LwsUEXSo7ijflE8m50WDcrCGlRO0zyaQrBqO+38EZh1pcYNWpzLGooOYFj70NJVIXC8LX5wbCTVpM5F+3PuudR/xt2BCP5nGLcfsIkrXNapE9zEpKWSoRD+L1Cwc/4JUxSjtHEzaqCDoXk6pHyHVgiWipMMh0JEaaunawdKvIji9O2UHDu9tv9wCAx9s7593mklrVnRfrRZQ9kABD9Lk27z/10a7aYJcl/E8CjAO4D8IaU4+4CcAHAhcOHD0/9x8gTQv0my08n/V8HbAgXcHikrNIkuTUhHx4bd+JIlUFAly7QUX9pS+6lhWLn3vDLsnTH1e15JJeHLA+k2Uw+IL3pDLTQxtqbNIZSKR4otBXBksyk38cwEaZG8gAeBPBYYHsvgFsAFAEUAHwSwH2jrmeW/GTIInRN/BwiqcMlQ2HX0of5fTEKQ1nycs21Nd/nNzeTUTtcgVKuLbwzjjy7dIZgyDue9eVDGlsobEo/wFotGfaknS96NZdQRI9M7eSzRzlXl+5hTh8zj64B0ADw2KjjjOR3jjSDSevg2tgTgmbrPGRdc76NLNCxsTHslNVZsuxYlWuwVJyWyRr6XkuBcSzfNDlEj+rtdjIBQaZnotHrOjXF4nBUjtSt2NiII3MKheT0Tz5TMuPSSiyHnERL94B3F7OSa9bp73sA3D/qHCP5yZHVN3S/F308FFnDhQ95eT427NgwrFSGfXWVSlwiWCz5tBLF8l7aoiOGMcAWM5MmSy6jUok5YqfVSkoyYrWvrAzHw0rDEeLnEE792Vn3rr3xZt1vC7Mi+f8G4GJfk/8ck37aZiSfxDjGTVbfCJ0/ymfHJM7yLa8nzdY6yz/6vXG+R1qpBcMYSAtH7PWG04P5hw5Z0vK+lATtdr31v7LiX0NOGJ1IIZtk2PJ00KJtpoqZyzXjbkbySYwrxY7bN5jQWS7l/XqWLoad9FcpByzH6QxZyYrlvJlJByLDhMgazUOrvoSmeNp5wqN2u+1Tmblh6Okhr2cr58pAw7G2HP0zSoIyjA0j+QXFpMZO1vu8LB4XMmTjr91ORtCsrTl36JD/W15FjpUoO46ck2sI8Ys1b/12htAWvCZ89qprLzxb4s7FA0EUxctzSVKDEL0s7yWNQBodSzzamx+SmsbNhDM454zkc4O0oAUdssyzaV3ud2Ul6V/rdJLh1BsbSYNNbzITD8mucn+SGKUX+TDMAEzyoSkVe8FlpNa1K2S/9q7zQruhqBtdpIhHe+1954Fo1DJkhiEYyecEaRKs9F0u56uX7xOHaqORtNa5n8q2sTFcClj6oAwQhUIyPj40Azdf2hxgVJROiFQ5+00PELz6y8qK/1u0u9XV2Psujc65ZHiW3IPMArhBpvkLbCo4EkbyOYP0A10dVupSyYId2prmmXjW2hOVipdqNNHLdUJhmCzRWt+cIbbzo6c9OM5mZStBGhk3GJFt2HKQMgvS0KrV4VhePTUNWQZmLYyEkfyCIysiRjs9ebFrHS4tM+NCwVvrehW4YnE4q1VH52lDKy1U2zAjhB7CdkfbUKqzWO2i26UtM8jWBA8W2lnL8bo6nl+gtccsx+2Swkh+gRCa2Wp5Rv4XyYQJVgwqHdLINWe0NMMrzMl5XO6EiV1LvFnhmYYZIPQQtjv6smzDhciEsOv1mMAlo01PC1mO4cYqVons1xY+NzCt3+sppMFIfpHAocehXBe2yFlrX18fXhCbE5mEtLUEw+HUzA082OjcmrTQbCP4OcV2H46OueXRXwia4+o5hEsvPC4avjTSlZV4VlCtDmuQ3MB0BJBZ8kMwkl8QcJx5uRxLodzmpQ9JgmGItLnwIGvzrMFvbCRnzZVKsrwvO25D1SLNybokCEXEMMlLA+S1Hmu1uOFxmeKQ86dSiRcb51h6XYZBh34aEjCSnxEmNTaEKFnWZB1dZr31engBD6kEy9Y8lxDma/I9aX3fuaRRNg5xm2G1BGDZhC10bmjSEHlxX5Fs1tbitSCFwLUDdpyiRxxTbw3OOZdN8gUYpoZTp4Dz5/3rODhxAqjVgHodqFSARgN49VX/3uuv+9fVVeDaNeAXfiF8jYsX/bnHj/vrHDjg98t1ajXg7Fng6NH4nNOn/TnlMvDkk8CRI8A73wl88YtAuw2cPDn63o8eBc6dS17XMOfY2vINZWtrvP1nzvhGdeUKcPAg0Gr5RlqtAh/8IFAs+uOuXQOuXgUKBeC114AXXgBeegl45RXf+D70IeC223wjO33aN7woAp5+GrjlFn/dN7zBX+vixeH7ACbvXMuMNPafxZZ3Sz4r+5ydp2KJ84xY5JTQ+3IML+0pBpFY8CLPhAwfHa1jssuSIE1jy8q606FVbHHzeo9s1e/bl2xY+nP1tFSHfYUcP2bJJwCTa+YDum3rmHPpQ9LmuY6MzJAlY1z8XBwZI9VgJUmp1Qqvoar7h8yCpeaMaPPWj3KOtAc8KutO/y/WhBA9lx/mTWrOa71dF0zScbydTvL+QiFoS95gjeTnBDpunZ2bEgLMTk6pE6MtbGnLYsmvriZDHgXSB6VevCZ1CVbQ4dBp1WhD32dJ+9RyYNRUVFslMoXUC4eLFq+dsGzVi34vjljd8NnyEGeR1NDRev6oaWgOG66R/AyQ1o44nl2MG37lkEZObmKHKVve2u/F7VuvBaHj7LloIJO9kPyovmARNQbnXNwoOYaeG65Y3XrFKV6FJmT1hywcHhTYguGqeFmLGOQ0a89IfgZIa0eS5S0ErduyJDhxIAGHD3c68Tk8K240hsMcQ9KqNr7S1lodBzk0iAw7gW60bElwLesoSg4CbOFLSGa97i32Y8diDbLTiWOHjx0bTuTQETvcuLPqb+QARvIzgG5H2qoWmUaMEZ7h8oLbUsdJ3ue+wRUh9cCgc0t0W+cl93LY5g2zQEgfD9Wz115+ncTB2+pq+DxZKd655OImemawJNl6RvJ7gKzIGV3LXfJGuBRHtxvnhrA/iSWZcjkZQcPviezJnxGScbSz1znT3g17AG5E3W6SiIXgQyVROftvfX2YwHu94YXI9bR4CRqukfweIESUuuSGZI+GjpV9UotdZrCc0SoyzerqcNSZvMcFAIFk0ENaEIJp74apIEsHZ62RtXwdWROqeS37jx1LZg4KoetZwRIsQGIkvwfgGaNASH7//jjs0bn0onos5Qhpi9VdKsXELzIL94eDB93Aouf1l6Vf7YSozZI3bAtZOrh2RvESgbxyjWj4IaLnrVKJO5Rek3YJLBQj+T0AZ2MzcbOMKFa6LrstK6ZxmyyVkmupSg0nljdDg0JakIERtWGqGJXpp9Hp+A6wtja8MIJIOqxxAsNyjd5EyxxlyeewMxjJ7wF0YhJr4BJsIO+Vy+HV0kKbkPz6+vDMt9mMjZy1tWH/lsGwZ9jOVFHO4Ww+55LWkUTWcDRBVodZWYkHA71GrZyvFzjOAYzk9wDcLkMSYCixr1hM7hNr/Pbb45wQyQiX47jssF7pyYjdsOfQ08tJGqGcK1l/1arfL+TPCxJLtMLa2nBGrXbYyjKEEsUgn6WTtnLUYYzk9wAcqhtqOzqHQ7bQMnuhIINKJZlpzvHzOTRMDIuCSS34kFTCZQrkGGnUEjkQKruatnGHqteT98naZ45gJL8HyJL5ej2vuYf8R6HEPtblo8hb8dz+dVmRHEqMhkXBpI1v3EGBp6kyU9CdJa1GjrbqJSZZZwtu9zvMIYzkp4BRpC7tSaJudNtjGUZHjUkJbl5ERy/ckYN2aVhGjNNwQ5l8HE7WaAxHHoSsJzleFhBP66zj1LyZ8w5nJD8FsGyiwTNDTfBco4bXVZC2q+PcedaZUznRYMiOqdcljLm4ku5g+/bFscxSQ0RHQzA4mSWrU815soiR/BSQRfIS816r+XwNLsfBpTW4AFmlMiw7ynvVamzJhwqRGQwLj1BMvThzeapbLCa1S3FKbW76V04w4evu3+/P5UQWdsaGOjJjWS15AL8M4HEArwO4Tb33MQCXADwJoD3O9RaJ5NOeucgzEveu10rlc1l21KG9ItnIYDDqcw2GhUaoYTNBc+cQCYfrePAMIIpiMpd6ITy1FmgrfoE71zRJ/u0AbgXw50zyAJoAHgGwH8AGgL8BUBx1vXki+VF+Gq2Py/8st4gxEpL8tOwoRosUJOOF7kcZGQZDLiGkG1rwQMfL62myLswkWyglXS+goFfEWgBMXa4JkPzHAHyM/j8P4Oio68wTyXPb0Bo4v6fbm5QV0LPCNCMlJLvIe3qxj7RrGQy5RigWXzod65yhUDPWTmUqLfvZUpNsRi7nukCa6CxI/tMA/iX9/wcAfmnUdeaJ5KUNhDRwLeUxya+ujhehNSo6J83hP+f+H4Nhb5C2BGAIWVaWWO16gYUFSx/fEckDeBDAY4HtvXTMtkkewF0ALgC4cPjw4T37Uca1iNOO44Ffy356tqeX2xunzaR9rlnyhqXDpFNhjSynLndOTsji8xYgnM3kmgB2mqjHkg23FZ49StvQy/WxI9VgMASQFVKp30/bx44vkWF0FmFa8smoqfycWVqzIPl3KMfrN+bN8bqTRL20XA0mey7zK+9xbsZ27sFgWBqErO9RHUUPBqyjiv7J02x9PP+f5YSdQ810mtE17wPwLQA/APBdAOfpvU/0o2qeBPDuca43T5q8Bs/yQoXI+P1Wy4c/al+PLtERirU34jcY3PY6grbEuf6NXgeTjw8tUciLPmiLTnfkOYAlQ02IrPBJGcTZUuf9bASMGuxDJD+HRoLBsHhgPX1UtExoms6aK4dg8nR9jmKbjeRHYJTeHjpWO955YBBrPmQ8ZH1u2j6DwTAhQvp6WmcM6f8ce687s5H84pG8tp41YY9LulkSn8FgmCHG7YwcNpcmzcyhJba0JL/TMEnnxo+iCuVcjDIgDAbDFKGt+XHIgK30Xi/WZaUWzpxiaUl+NyxpXmoyy9GfJvGYNW8wzAhZnS+tIwvJS7q5Ljw1p8gi+X3IMU6eTL5uB/fdB9y4ARSLwIkTwD33AA89BFy+DHz5y8nPunzZ/33iBHD8uH+9fBlotXZ2DwaDYRvQBHDvvcDHPw6cPu079kMP+U589ixw9Kg/5swZ4D3vAZ57zv/fbvvz5f1FRBr7z2Kbl+gahl6gZhyfC0s8ZsUbDHMCjn3WHZsxqubIHKaiI8OSL8x6kNktbG1563lra3eve+aMt8SbTW+V33GHH9zPnEk/5+RJf8zp07EhYDAYZozTp4F63b/qjr21FZMIAJw7F7beT50Czp/3U3omHNl/6tTefZ9xkcb+s9h2YslPU/sOFQybQwe7wWDIwqgaOOOQiE5t14Qwo6JmWAbH6zRJNxRhs0C1iwwGg3Oja+BMEo43KltyjzXapSD5aSItack0d4NhgbBblmAWkc+o5EEWyedGk58mjh5NSnRbW156M83dYFggcEfWTrytLeDIEb+NcuyJ0y3U8R94wEfmPPDA7t//dpHG/rPY5tGS32kpa4PBMIdIS0/faceekbMOyxonvxsQpzngjQBgd+LvDQbDDKE7MSe67KRjy2xhjhD5QWA+cNttt7kLFy7M+jYSEGlm0fMhDAbDFDAnBBFF0Vedc7eF3jNLfgTmcGA2GAzzgtBUf85gjlfCtBKqDAZDTpHlhJ0TGMkTOGnNCN9gMIyEDr3bDqZMNibXENgHswCzMIPBkAdMmWyM5Amsv1sEjcFg2BNMmWyWRq6ZdEa0G7Mwg8FgGIkpk83SWPIyI7p8GbjppplHPBkMBsOeYGlIXmZCly+b1m4wGJYHSyPXyIzozJm5j3gyGAyGXcPSWPICS24yGAzLhKWx5A0Gg2EZYSRvMBgMOYaRvMFgMOQYRvIGg8GQYxjJGwwGQ45hJG8wGAw5hpG8wWAw5BhztTJUFEX/AOBvZ30f28QbAXxv1jexh7Dvm18s03cF8vF93+Kcuzn0xlyR/CIjiqILactv5RH2ffOLZfquQP6/r8k1BoPBkGMYyRsMBkOOYSS/e7h31jewx7Dvm18s03cFcv59TZM3GAyGHMMseYPBYMgxjOQNBoMhxzCS3yGiKPrlKIoej6Lo9SiKblPvfSyKoktRFD0ZRVF7Vvc4LURR9FtRFH07iqKv9befn/U97TaiKDref36Xoij66KzvZ9qIouibURRd7D/PC7O+n91GFEX3RVH0bBRFj9G+Q1EUfT6Kor/uv75hlve42zCS3zkeA3ACwF/wziiKmgB+BcA7ABwH8HtRFBX3/vamjjPOuR/vb/9n1jezm+g/r/8E4N0AmgDe33+uecfP9p9nHmPH/yt8f2R8FMAXnHM/CuAL/f9zAyP5HcI595fOuScDb70XwP3OuR84554CcAnA7Xt7d4Yd4nYAl5xz33DOXQdwP/xzNSwonHN/AeB5tfu9AD7T//szAH5xT29qyjCSnx7eDOAZ+v9b/X15w4ejKHq0Pw3O1TQXy/MMGQ7An0VR9NUoiu6a9c3sEW5xzn2n//ffA7hlljez21i66NkIqQAAAadJREFUNV63gyiKHgSwFnjrE865/7XX97OXyPruAP4zgN+GJ4bfBvAfANyxd3dnmAJ+yjn37SiK3gTg81EUfb1v/S4FnHMuiqJcxZUbyY8B59y7tnHatwH8I/r/h/v7Fgrjfvcoiv4LgD+d8u3sNXLxDCeBc+7b/ddnoyj6LLxklXeS/24URevOue9EUbQO4NlZ39BuwuSa6eFzAH4liqL9URRtAPhRAA/P+J52Ff0OIXgfvBM6T/gKgB+NomgjiqISvCP9czO+p6khiqJKFEUH5W8Ax5C/ZxrC5wD8av/vXwWQq9m5WfI7RBRF7wPwHwHcDOB/R1H0Nedc2zn3eBRFfwTgCQCvAviQc+61Wd7rFPCpKIp+HF6u+SaAD872dnYXzrlXoyj6MIDzAIoA7nPOPT7j25ombgHw2SiKAM8N/8M5d262t7S7iKLoDwH8MwBvjKLoWwBOAvh3AP4oiqI74Uud//PZ3eHuw8oaGAwGQ45hco3BYDDkGEbyBoPBkGMYyRsMBkOOYSRvMBgMOYaRvMFgMOQYRvIGg8GQYxjJGwwGQ47x/wHvIeM/edfvBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.arange(-10,10,.01)  + np.random.normal(0,1,2000)\n",
    "y1 = x + np.random.normal(0,1,2000)\n",
    "y2 = -x + np.random.normal(0,1,2000)\n",
    "\n",
    "plt.scatter(x, y1, color='b', s=2)\n",
    "plt.scatter(x, y2, color='r', s=2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.empty((4000,2))\n",
    "X[:2000,0] = x\n",
    "X[:2000,1] = y1\n",
    "\n",
    "X[2000:,0] = x\n",
    "X[2000:,1] = y2\n",
    "\n",
    "np.random.shuffle(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectors = [[1,0],[0,1]]\n",
    "n_iterations = 10\n",
    "history = [vectors.copy()]\n",
    "assignments_history = []\n",
    "\n",
    "for iteration in range(n_iterations):\n",
    "    # E Step\n",
    "    distances = np.abs(X.dot(vectors))\n",
    "    assignments = np.argmax(distances, axis=1)\n",
    "    assignments_history.append(assignments)\n",
    "    \n",
    "    #M Step\n",
    "    #np.eig(A)[0]->eig_vals, np.eig(A)[1]->eig_vectors, we need 1st eigenvector\n",
    "    v1 = np.linalg.eig(np.cov(X[assignments==0].T))[1][0]\n",
    "    v1 = v1 / np.sqrt(np.sum(v1**2))\n",
    "    v2 = np.linalg.eig(np.cov(X[assignments==1].T))[1][0]\n",
    "    v2 = v2 /np.sqrt(np.sum(v2**2))\n",
    "    vectors = [v1, v2]\n",
    "    \n",
    "    history.append(vectors.copy())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "history = np.array(history)\n",
    "for i in range(len(history)):\n",
    "    fig = plt.figure(figsize=(10,6))\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_xlim(-2,2)\n",
    "    ax.set_ylim(-2,2)\n",
    "    ax.set_title(\"Clustering along lines\")\n",
    "    ax.plot([0,history[i,0,0]], [0,history[i,0,1]], 'b--')\n",
    "    ax.plot([0,history[i,1,0]], [0,history[i,1,1]], 'r--')\n",
    "    fig.savefig(\"plot_{}.png\".format(i))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "pngs = [Image.open(\"plot_{}.png\".format(i)) for i in range(10)]\n",
    "\n",
    "gif = Image.new(\"RGBA\", (pngs[0].width, pngs[0].height), (255,255,255))\n",
    "\n",
    "gif.save(fp=\"comb.gif\", format='GIF', append_images=pngs,\n",
    "         save_all=True, duration=1000, loop=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![test](https://s5.gifyu.com/images/combd74ef49de9a11965.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "assignments_history = np.array(assignments_history)\n",
    "for i in range(len(assignments_history)):\n",
    "    fig = plt.figure(figsize=(10,6))\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_title(\"Clustering along lines\")\n",
    "    ax.scatter(X[assignments_history[i]==0, 0], X[assignments_history[i]==0, 1], color='b', s=2)\n",
    "    ax.scatter(X[assignments_history[i]==1, 0], X[assignments_history[i]==1, 1], color='r', s=2)\n",
    "    fig.savefig(\"scatter_plot_{}.png\".format(i))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "pngs = [Image.open(\"scatter_plot_{}.png\".format(i)) for i in range(10)]\n",
    "\n",
    "gif = Image.new(\"RGBA\", (pngs[0].width, pngs[0].height), (255,255,255))\n",
    "\n",
    "gif.save(fp=\"comb_scatter.gif\", format='GIF', append_images=pngs,\n",
    "         save_all=True, duration=1000, loop=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![test](https://s5.gifyu.com/images/comb_scatter.gif)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}