Gerrymandering solution (HackMIT Puzzle Hunt 2018)

Gerrymandering.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>voters_by_block</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>party_A</th>\n",
       "      <td>[0, 9342, 147565, 81158, 32798, 54066, 15916, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>party_B</th>\n",
       "      <td>[6924, 21879, 60357, 48212, 40474, 37607, 5134...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                           voters_by_block\n",
       "party_A  [0, 9342, 147565, 81158, 32798, 54066, 15916, ...\n",
       "party_B  [6924, 21879, 60357, 48212, 40474, 37607, 5134..."
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.read_json('voters.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pa, pb = pd.read_json('voters.json').voters_by_block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pa = np.array(pa).reshape((10, 10))\n",
    "pb = np.array(pb).reshape((10, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[     0   9342 147565  81158  32798  54066  15916  11472     30  34549]\n",
      " [ 22546   9545  91396  54423  23600   3516  14781  56783  51528  21106]\n",
      " [ 10335 110916   3850      0  28175  11319  32015   3613  43891 107716]\n",
      " [  2537  20001   1551  11495  36812  21660  34372  22045  98262  79401]\n",
      " [ 23489   6086  36712  15739    963  24644   2812  54262  91882   8103]\n",
      " [   155  40104   1796  15550 124579  90032 216559   4500  21988  44536]\n",
      " [  8700  15032      0     84  14453      0  37881  25023  28917  39379]\n",
      " [ 16071      0  99449  75566   4729  53554 169628      0   4300  24676]\n",
      " [  6646  12693  27422  51045  27725   1811  34193  92718  30292  40070]\n",
      " [ 34455   6966  61631 113609   4743  22771  35197  83484   8078   4218]]\n",
      "[[  6924  21879  60357  48212  40474  37607   5134  15904 121612  77183]\n",
      " [ 20320  14360 185957  53870  75933   1531  29939  47286   3189  52423]\n",
      " [ 29587 331906  20971 122268      0  59639 141474  23432 141645  97232]\n",
      " [  4356  34300   3955  11297 107940  73371  65330  25874 133599  68691]\n",
      " [ 35143  11873 160269  21626   1523  10544   3012  65759  58141   7980]\n",
      " [  3419  15332  13295  11002 124911  26349 103241   5950  28931  91693]\n",
      " [ 18034  14960   3887    117  33855 162086  81009  29056  37306  52465]\n",
      " [ 93901  74963 105269  88581  85313  95589  99885 198307  10834  43071]\n",
      " [ 52401  16590  25088  91794  72238 127455 125767 264146  71503  24481]\n",
      " [ 34769   8756  84916 181061   9347  25309  52771  68561   7159  10667]]\n"
     ]
    }
   ],
   "source": [
    "print(pa)\n",
    "print(pb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10cc2ba90>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD3CAYAAAA9vL6wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFu1JREFUeJzt3X20XWV94PHvCXkDmoQgBKTFUqT8wK6FtIKJEELoQJGgtaVrWtZSUXGoLvEFhzWiCNW2qHXG6hoKXVgUQUpmbBFHzEyE4UUmRCFLBqyM+gsvFco7pHkhgQTuPWf+2Cf0rhjvOefm7J2zz/l+WHvl3H2fs3/PJTe/+7vPfp79NFqtFpKk+pi2uzsgSeqNiVuSasbELUk1Y+KWpJoxcUtSzUwv9eIzf7WSKStzZ+1VRRgAHnn7ayuLde3/3K+6WM3HK4t1xIx9K4nzvzdlJXEAto2/XFms9S9urizWm/Y/orJYqx6/tbGr13j5uYe7zjkz9jt0l+PtLlbcklQzpVbcklSp5vju7kElTNyShsf42O7uQSVM3JKGRqvV3N1dqISJW9LwaJq4JalerLglqWa8OSlJNTMiFXfX87gjwjnfkgZaa3ys66POJq24I+JQ4IvAMcBYO3n/GPhoZq6toH+S1D1vTgLwFeATmXn39hMRsQj4GnB8mR2TpJ45VALA7IlJGyAz7yqxP5I0dc3x7o8a61Rx/ygirgK+C2wE5gDLgH8qu2OS1LMRqbg7Je4PAH8ALAbmApuAFcC3Su6XJPWuTzcdI2IGcBVwCDALuAT4CXA10ALuB87NzGZEfAo4HRgDzsvMNRFx2K62nax/kybuzGxRJGkTtaTB17+bk+8A1mXmOyNiX+C+9nFRZn4vIq4A3hYRjwAnAguBg4FvAsdSTOqYcls65FzncUsaGq1W38au/xG4vv26QVEhvwG4o31uJfB7QAI3t4vcRyNiekTs34e2Jm5JI6JPY9yZuRkgIuZQJPCLgC+0ky7A88A8iiHkdRPeuv18YxfbTspFNZKGR7PZ/dFBRBwM3A5cm5nLgYlvmgNsoLjvN2cn53e17aRM3JKGR6vZ/TGJiDgAuBm4IDOvap++NyKWtl+fBqwCVgOnRsS0iHgNMC0zn+tD20k5VCJpePRv788LgfnAxRFxcfvcR4BLI2Im8FPg+swcj4hVwA8oCuFz223PB66cattOnWu0WuXt5+tmwbvGzYJ3jZsF77q6bRa89a5vdJ1zZi/6k9puFlxqxX3YPgeVeflXPPNixyGhvrnp26+qLNYPZ22tLNaGbVsqi/Xf1j1USZxX7Tmnc6M+OWP+UZXF2rJPdav+7n7hkcpi9YULcCSpZnzIlCTVjIlbkuqlVeG9ht3JxC1peDjGLUk141CJJNWMFbck1YwVtyTVjBW3JNXMWL13b++WiVvS8LDilqSacYwbIuJ2iv3WJmoArcw8rrReSdJUWHED8HHgSuAPKbbukaTBZcUNmXl3RFwLHJWZbhgsabBZcRcy879U0RFJ2mXOKpGkmilxY5hBYuKWNDwc45akmjFxS1LNeHNSkmpmvLr9OHenUhP3c1s3lnn5V0xrNLhlfjW7US96elUlcQDGmtV9Ex71qt+oLNacWXtWEufgPfevJA7ANU/dXVmsUxZUtzHx8pmvqSxWXzhUUh9VJW1JA87ELUk14xi3JNVLq+k8bkmqF4dKJKlmnFUiSTVjxS1JNWPilqSa8SFTklQzI1JxT+v1DRGx41ZmkjQYmq3ujxr7pRV3RLwVuAx4GfhkZn6j/amVwO9W0DdJ6k2fZ5VExELg85m5NCIWUGzlOB/YAzgrMx+KiHOA91Fs73hJZq6IiP2A5cCewBPAezLzhV7aTtavySruTwJHAwuB90XEu9rnG1P5HyBJZWs1m10fnUTEx4CvALPbp/4zcF1mLgEuAo6IiAOBDwPHA6cCn2uPSvwZsDwzTwDupcihXbft1LfJEvdLmbk+M9cBbwM+GBEnAfX+HUPS8OrvUMlDwBkTPj4e+LWIuAV4O/A94I3A6szclpkbgQeBo4DFwHfb71sJnNxj20lNlrh/HhFfjIi9M/P59hdwOeATnSQNplaz+6ODzPwmxVDxdocA6zPzZOBR4AJgLjDxMajPA/N2OL+zc53aTmqyxH028E+0K+zM/BfgJOAfOl1UknaLcm9OrgNubL/+DnAMsAmYM6HNHGDDDud3dq5T20n90puTmTkGXL3DuaeB8zpdVJJ2i7FSl7zfCSwDrgWWAP8PWAN8JiJmA7OAI4H7gdXttlcDpwGremw7qZ6nA0rSwOrjUMlOnA+cFRHfB94MfDYznwIupUi2t1HMwNsKXAKcGRGrgTcBl/XStlNHXIAjaXj0eX52Zv4cWNR+/Qhwyk7aXEkxTXDiuacpkvuU207GxC1paHQzzW8YmLglDY+ar4jslolb0vAwce+6TdsmXbXZN3+05bFK4gCMV7jz+vw9f6WyWPP2qGbndYDNM+d0btQHMxp7VBIH4JC5B1QW60Pbqvu+uGtmdbXdsf24iBspSFK9uOekJNWNiVuSasZZJZJUM1bcklQzJm5JqpfWuEMlklQvVtySVC9OB9yJiNgTaGbmtpL6I0lTZ+KGiHgd8FlgPXAdxf5r4xHxkcxcUUH/JKl7ozHE3bHivgK4mGLLnuuBw4GtFPuimbglDZTW2Ghk7k6Je1pm3gHcEREnZeYzABExVn7XJKlHo5G3OybujIivAH+ame8GiIiPA0+V3TFJ6pU3JwvnAG/NzIk/xx6j2H5HkgaLFTe0E/a3dzj396X2SJKmyIpbkurGiluS6qU1ItMmTNyShkbLiluSasbELUn1YsUtSTVj4u6D1847qMzLv+JfNj9bSRyAxQteV1msB154srJYdz7zk8pixfyDK4mzZt3aSuIAvP3AhZXF2rStut3rL1r/g8pifagP12iNN/pwlcFnxS1paFhxS1LNtJpW3JJUK1bcklQzrZYVtyTVihW3JNVM01klklQv3pyUpJrpd+KOiIXA5zNzaUQcDfwNMA5sA87KzKcj4hzgfcAYcElmroiI/YDlwJ7AE8B7MvOFXtpO1q9pPXwBC3r8miWpUq1W90cnEfExig3SZ7dP/VfgQ5m5FLgBuCAiDgQ+DBwPnAp8LiJmAX8GLM/ME4B7gff10rZT335pxR0Rh+9w6usRcRZAZla3JE2SutTnivsh4Azg2vbHZ2bm9uXM0yk2Tn8jsDoztwHbIuJB4ChgMfDZdtuV7dcP9dD2S5N1bLKhkluAFyhK9wYQwJeBFvC7nb9mSapWP6cDZuY3I+KQCR8/CRARxwEfBJZQVM4bJ7zteWAeMHfC+Z2d69R2UpMNlRwD/AT4XGaeBNyXmSdlpklb0kAaH290fUxFRPwJcAVwemY+C2wC5kxoMgfYsMP5nZ3r1HZSvzRxZ+YzwB8Dp0fEhZ2/JEnavVqtRtdHryLiHRSV9tLMfLh9eg1wQkTMjoh5wJHA/cBqYFm7zWnAqh7bTmrSm5OZOZaZ51EMl3R9I1OSdodWs9H10YuI2AO4lKIiviEivhcRf56ZT7XPrwJuAz6ZmVuBS4AzI2I18Cbgsl7adupPo9XN7dUpOnLBGyvZcrnKx7oeu+9hlcWq8rGuT29eX1msqh7r+sDGxyuJA9U+1nXZtlmVxXrv83dVFmvj5od2eYD6p7+5rOucc+QD/6u2k76dxy1paLgAR5JqZrw5GiO6Jm5JQ6PEkd+BYuKWNDSaPtZVkurF53FLUs04VNIHD2yoZjrWvzvgqEriANz+zI8ri7Vg730qi/Xw64+oLNb3HzuwkjhnPLa8kjgAcw8+qbJYd/7K/pXFWrrvkZXF6geHSiSpZpxVIkk1MyIjJSZuScPDoRJJqhlnlUhSzYzIJu8mbknDo4UVtyTVyphDJZJUL1bcO4iIacCrgSczc1SGkiTVyKgkpklnq0fEV9t/LgTWUmxJf39ELKqgb5LUkxaNro8667TM6Dfaf34GOC0zFwInA58vtVeSNAXNHo4663Z96HhmPgCQme4/KWkgjdPo+qizTmPc8yLiHmDviHgvcB3w18AjpfdMkno0IjuXTZ64M/MNETELeD3wAsVvGD8GvlpB3ySpJ82aV9Ld6jirJDO3AWsmnLqivO5I0tT5kClJqpm633Tslolb0tBoNhwqkaRaGd/dHaiIiVvS0HBWiSTVjLNK+mDOzD3LvPwr1qx/gL1mzKokVrPCbaSf2ry+sliLHqzu63p6y88qiXPi0e+vJA7Aa+e+urJYs6fNqCzWi62xymL1g7NKaqSqpC1psDlUIkk143RASaqZcStuSaoXK25Jqpl+Je6ImAFcAxxCMT38HGAMuJriHuj9wLmZ2YyITwGntz9/XmauiYjDum07lf75eFZJQ6PV6P7oYBkwPTOPA/6CYk+CLwIXZeYJQAN4W0T8DnAisBA4E7i8/f5e2vbMxC1paPRxI4W1wPT2lo1zgZeBNwB3tD+/kmJTmcXAzZnZysxH2+/Zv8e2PXOoRNLQ6OOS980UwyQ/A/YD3gIsycztU8WfB+ZRJPV1E963/Xyjh7bP9to5K25JQ6PZ6P7o4KPATZl5OMV+BNcAMyd8fg6wAdjUfr3j+WYPbXvWU+KOiP0iYkQm3Eiqmz4OlawHNrZf/yswA7g3Ipa2z50GrAJWA6dGxLSIeA0wLTOf67FtzyYdKomI9wAHAyuA5cBWYK+I+EBm3jKVgJJUlj5OB/wScFVErKKotC8EfghcGREzgZ8C12fmeLvNDygK4XPb7z+/h7Y96zTG/QFgKXAj8PuZuTYiDgK+DZi4JQ2Ufj2rJDM3A3+8k0+duJO2nwY+vcO5td22nYpOQyUvZ+YWikH0h9uBn2B0nuUiqUb6OMY90DpV3DdGxLcpJpCviIibgDcDt5XeM0nq0ahspDBpxZ2Zf0UxkbwBPAosAC7NzI9X0DdJ6kmTVtdHnXWzy/sd/NtEckkaWD6rRJJqpt51dPdM3JKGhhW3JNXMWGM0am4Tt6ShMRpp28QtaYg4VNIHe8+cXeblX1Hlbugr5y+uLNaPZlW3m/fFz66qLNY7D1pUSZzLLziokjgAcz9yQ2WxZuxRXb1V5e71/VD3aX7dsuKWNDRGI22buCUNEYdKJKlmxkek5jZxSxoaVtySVDMtK25JqhcrbkmqGacDSlLNjEba7vA87oiYW1VHJGlXjdHq+qizTluXPRUR762kJ5K0i1o9/FdnnRL3j4DfjojbIuIXNr6UpEHS7OGos05j3C9m5gcj4hjgExFxGXAr8HBmXlp+9ySpe3WvpLvVKXE3ADLzh8AfRcQ8YAkQZXdMknpV90q6W50S99UTP8jMjcB32ockDZTxlhU3mXlNVR2RpF3lPG5JqhnHuCWpZhzjlqSacahEkmrGoRJJqhlnlUhSzThU0gcbtm4p8/Kv2G+v6p6F9bHGY5XF+q3GAZXFOmLer1UW67on764kzoOXVLdO7B/2re6JEJ+bVt334L577F1ZrH7w5qQk1Uy/x7gjYgFwD3AKMEaxKLEF3A+cm5nNiPgUcHr78+dl5pqIOKzbtlPpV6eHTElSbTRpdX10EhEzgC8DL7ZPfRG4KDNPoHgcyNsi4neAE4GFwJnA5VNo2zMTt6Sh0Wq1uj668AXgCuCJ9sdvAO5ov14JnAwsBm7OzFZmPgpMj4j9e2zbMxO3pKExTqvrYzIR8W7g2cy8acLpRmZuf+PzwDxgLrBxQpvt53tp2zPHuCUNjT7OKjkbaEXEycDRwNeBBRM+PwfYAGxqv97xfLOHtj2z4pY0NPo1VJKZSzLzxMxcCtwHnAWsjIil7SanAauA1cCpETEtIl4DTMvM54B7e2jbMytuSUOj5Hnc5wNXRsRM4KfA9Zk5HhGrgB9QFMLnTqFtz0zckoZGGUve21X3dr8wYT8zPw18eodza7ttOxU9Je72T489MvPFjo0lqWIueQci4nDgs8BLwKUUA/TTI+ITmfmNCvonSV1zyXvhSuAvKaasrABeT3EX9BbAxC1poIxK4u40q2R6Zt4C3ACsy8zHM3ML8HL5XZOk3vR5Ac7A6lRx/zwi/nu73eaI+AzFBPInS++ZJPVoVCruTon7XcAyYC2wGfgo8ALF5HRJGihupABk5hhw44RT55fbHUmauvHWaDzY1XnckoZG3ceuu2XiljQ0HOOWpJpxjFuSaqbpUIkk1YsVtyTVjLNK+uC35x9a5uVf8aMN/1xJHIAl+xxcWaxzpm2uLNaxG6vbOfzQea+uJE6jkiiFddOre7T92g2PVxZr21i9Fkk7VCJJNeNQiSTVjBW3JNWMFbck1cx4a3x3d6ESJm5JQ8Ml75JUMy55l6SaseKWpJoZlVklXa8aiIgq1zNIUs9aPfxXZ512eX8tcDlwJHBQRNwDPAz8x8x8qoL+SVLXRmXJe6eK+3Lgw5n568AJwO3AXwNfLbtjktSrUdksuFPinpeZawEy8y7g+My8B5hfes8kqUfNVqvro8463Zx8OCKuAFYCbwF+GBGnA1tK75kk9ajulXS3OlXc7wF+DPwesAb4T8A64MyS+yVJPWvS6vqos067vL9EMc490V3ldUeSpm5UKm7ncUsaGqMyq8TELWlo1P2mY7dM3JKGhkMlklQz/VoRGRHTgL8FXg9sA/5DZj7Yl4v3QXUb5UlSyfq4AOcPgNmZ+Sbg4xQLDweGiVvS0OjjApzFwHfhlcWHx5Td916UOlSy6vFbfTBVTby4uzugrp2zuzswwMZeerxfOWcusHHCx+MRMT0zx/p0/V1ixS1Jv2gTMGfCx9MGJWmDiVuSdmY1sAwgIhZRrCAfGM4qkaRf9C3glIj4PtCgePzHwGiMyrxHSRoWDpVIUs2YuCWpZkzcklQzA3NzcncsMY2IhcDnM3NpiTFmAFcBhwCzgEsy88aSYu0BXAkE0ALen5n3lxGrHW8BcA9wSmb+rMQ4/5diehbAP2dmaTeKIuITwO8DM4G/zcxStumLiHcD725/OBs4GjgwMzf0Oc4M4BqK779x4Jyy/q4iYhbwNeBQir+vczPzgTJijbpBqrgrXWIaER8DvkLxj6ZM7wDWZeYJwJuBy0qM9VaAzDweuAj4TFmB2gnhy5S8diciZgONzFzaPspM2kuB44DjgROBg8uKlZlXb/+aKH74fbjfSbttGTA9M48D/oISvyco1gZtzsxFwIco93t9pA1S4q56ielDwBklxwD4R+Di9usGUNok/sz8H8Cftj/8daCMRLDdF4ArgCdKjAHFb2B7RcTNEXFbe05tWU6lmK/7LeA7wIoSYwEQEccAv5WZf1dSiLXA9PZvtHOBl0uKA/A6im0OycwEjiwx1kgbpMS90yWmZQXLzG9S7jfx9jibM/P5iJgDXE9RCZcZbywirgH+BriujBjtX/Ofzcybyrj+Dl6g+CFxKvB+4LoSvy/2oygY/v2EWGU/tuFC4M9LvP5mimGSn1EMo11aYqz7gLdERKP9A/ZX28N36rNBStwDvcR0V0TEwcDtwLWZubzseJn5LuBw4MqI2LuEEGdTLE74HsXY7Ncj4sAS4kBRMf59ZrYycy3FnqevLinWOuCmzHypXTFuBfYvKRYRsQ8QmXl7WTGAj1J8TYdT/PZyTXv4qQxXUfw7XgX8IXBPZo6XFGukDVLiHuglplMVEQcANwMXZOZVJcd6Z/vmGhSVarN99FVmLsnME9vjs/cBZ2XmU/2O03Y27fsdEXEQxW9mT5YU607gze2K8SBgb4pkXpYlwK0lXh9gPf/2m+y/AjOAsqrgY4FbM3MxxRDhwyXFGXkDM6uEAV9iugsuBOYDF0fE9rHu0zKzjJt6NwBfi4j/Q/EP9LyS4lTpq8DVEXEnxUyZs8v6TSwzV0TEEmANRVFzbskVY1B+cvsScFVErKKYKXNhZm4pKdYDwF9GxCcp7q+8t6Q4I88l75JUM4M0VCJJ6oKJW5JqxsQtSTVj4pakmjFxS1LNmLglqWZM3JJUM/8f30k+WO7+5s8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cabef60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(pa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10cdb68d0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD3CAYAAAA9vL6wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAF0tJREFUeJzt3XmQXWWZx/HvzdYdoBOWCBGLdRgeQWWTIUASCBYadpVxZigLQTaxQJkoM8gqLixDOaIiIE7YlRSULIopQQoCGCIYRRAo4QkQFgeIkChZyNZ975k/zknsacO953Tf8/Z97/19UqfS9/Tb53mTdJ5+7nve97yVJEkQEZF4jBjuDoiISDFK3CIikVHiFhGJjBK3iEhklLhFRCIzqsyLjx27XZApK5t2bRwiDAAf6tk2WKxl1dXBYq2qrQ0Wq5rUgsR5fulrQeIA9IwZGyzWheP3CRbr0RHvBIs165W7KkO9Ru/ihblzzugJOw453nBRxS0iEplSK24RkaBq1eHuQRBK3CLSPqp9w92DIJS4RaRtJIHunww3JW4RaR81JW4Rkbio4hYRiYxuToqIRKZDKu7c87jNTHO+RaSlJdW+3EfM6lbcZrYjcDmwN9CXJe+ngS+5+4IA/RMRyU83JwG4FjjH3X+z7oSZ7QvcAEwus2MiIoVpqASA7v5JG8DdHyuxPyIig1er5j8i1qji/oOZXQ/cCywFeoDDgKfK7piISGEdUnE3StynAZ8ApgDjgGXAbOCukvslIlJc5Dcd86qbuN09IU3SStQi0vp0c1JEJC5JEvfYdV5K3CLSPjTGLSISGQ2ViIhEpkkVt5mNBGYCBiTA54HVwI3Z62eA0929ZmYXAocDfcAMd59vZjsNtW29/mkZu4i0j2pv/qO+IwHcfTJwPnAx6Sry8919KlABPm5mewEHApOAY4Crsq8fUttGnVPiFpH2UavlP+pw958Cn8tebge8DXwYeDg7dw9wMOlU6fvcPXH3V4FRZvaeJrStq9ShkvFdG5V5+fXWNP7p2TQ/3DzcXevjF48MFmuHri2Dxbr9jd8GibPtuK2CxAHYYkxPsFhnvvlQsFiTJliwWE3RxJuT7t5nZjcBnwQ+BXw0myINsBwYT7q+ZUm/L1t3vjLEtnWp4haR9tGkinsddz8e2Jl0vHtsv0/1kFbhy7KPB56vDbFtXUrcItI+mpS4zewzZnZO9nIlaXL9nZlNy84dCswF5gHTzWyEmW0LjHD3xcATQ2xbl2aViEjbSJo3bHoncIOZ/QoYDcwAngVmmtmY7OPb3b1qZnOBR0kL4dOzrz9zKG0bda6SJEmjNoM2cdNdyrt4P2sDPp/g8R23Dxbr+MXh5qS+b1S4MVqNcQ/Nk0teDBYr5Bj33NceqAz1GqsevDZ3zhl70MlDjjdcVHGLSPvQAhwRkchoybuISGRUcYuIREYVt4hIZPq0kYKISFxUcYuIREZj3GBmDwJdA05XgMTd9y+tVyIig6GKG4CzSdfof5L0+bEiIq1LFTe4+2/M7EfAbu6uDYNFpLWp4k65+7dCdEREZMg0q0REJDIlPnuplShxi0j70Bi3iEhklLhFRCKjm5MiIpGphtsTdjiVmrg3GtVd5uX7xQF/7o4gsbbd6YggcQDO69krWKwZbzwYLNa0rT4YJM5pvVsEiQNwxuqngsW6dbMDgsV6ZHRkFayGSuIRKmmLSItT4hYRiYzGuEVE4pLUNI9bRCQuGioREYmMZpWIiERGFbeISGSUuEVEIqOHTImIRKZDKu4RRb/AzAZuZSYi0hpqSf4jYu9acZvZkcCVQC9wnrvfln3qHuAjAfomIlJMh8wqqVdxnwfsAUwCTjWz47PzldJ7JSIyCEmtlvuIWb0x7rXu/lcAM/s4MMfMXgXifo8hIu0r8iGQvOpV3C+b2eVmtrG7LweOBq4C3h+mayIiBSW1/EfE6lXcJwLHklXY7v4nMzsIOCdEx0RECmtSxW1mo4Hrge2BLuAid787+9yngS+6+37Z61OAU4G+rN1sM5sAzALGAq8DJ7j7yiJt6/XvXStud+9z9xv7X8Dd/+zuMwbzFyEiUrq+av6jvmOBJe4+FTiEdKIGZrYncBLZvT4zmwicAUwGpgOXZjPvvgrMyr7+CdL7hLnbNupc4emAIiItq3lDJT8BLsg+rgB9ZrYFcAnQv3jdB5jn7mvcfSnwArAbMAW4N2tzD3BwwbZ1aQGOiLSPJg2VuPsKADPrAW4nTeLXAV8GVvVrOg5Y2u/1cmD8gPMbOteobV1K3CLSNpo5zc/MtgHuAq4Gngf+EfgB0A3sambfBeYAPf2+rAd4G1iWfbxqA+fytK1LiVtE2kfzbk5uBdwHfMHdH8hOfyD73PbAre4+Ixu3vtjMuklvYu4CPAPMAw4DbgQOBeYC8wu0rUtj3CLSPpq35P1cYDPgAjN7KDvGDmzk7ouAK0iT7RzSVeargYuAY8xsHrAfcGWRto06V0lKfJpWV/c2QWbD7zBuYogwALy0bFGwWF0jRweL9R8T9gsW62lWBImzKukLEgfgqRWvBot1SffuwWL9ZFTDd+1N87NXZw95VfaKLx+VO+dscvnd0a4C11CJiLQN7TkpIhIbJW4RkchE/vCovJS4RaR9qOIWEYmMEreISFySqoZKRETioopbRCQumg64AdnKoZq7rympPyIig6fEDWa2K+ljDP8K3AJcC1TN7N/dfXaA/omI5NcZQ9wNK+5rSB9nuD3pow13BlaTPjNWiVtEWkrS1xmZu1HiHuHuDwMPm9lB7v4mgJmFewiEiEhenZG3GyZuN7Nrgc+5+2cBzOxsINyTlkREctLNydQpwJHu3v/n2P+SPppQRKS1qOKGLGH/bMC5H5faIxGRQVLFLSISG1XcIiJxCbh3xrBS4haRtpGo4hYRiYwSt4hIXFRxi4hERom7CXYc994yL7/e6yuXBIkDMGfTScFifWr1H4PF+uYbDwWLNXZ0V5A4q/vWBokD8LGtwu28fkPtzWCxHlu8IFisZkiq0W7cXogqbhFpG6q4RUQik9RUcYuIREUVt4hIZJJEFbeISFRUcYuIRKamWSUiInHRzUkRkch0SuIekbehmW1ZZkdERIYqSfIfMXvXitvMdh5w6mYzOw7A3eNaTiUiHaHZFbeZTQIuc/dpZrYH6QbqfcAC4GR3r5nZKcCp2fmL3H22mU0AZgFjgdeBE9x9ZZG29fpVr+K+H7g76+gPAct+v2ZwfwUiIuVKkkruoxEzOwu4FujOTl0IfMPdpwBdwOFmNhE4A5gMTAcuNbMu4KvALHefCjwBnFqkbaO+1UvcewN/BC5194OAJ939IHf/SMM/sYjIMKhWK7mPHF4Eju73+glgczOrAD1AL7APMM/d17j7UuAFYDdgCnBv9nX3AAcXbFvXuyZud38T+FfSnyrn5vlTiogMp2ZW3O5+B2lyXud50o3SnwW2Ah4CxgFL+7VZDowfcH5D5xq1ravuzUl373P3GaTjLrlvZIqIDIekVsl9DML3gKnu/n7gZuDbwDLS6nudHuDtAec3dK5R27pyJWN3v9HdD8zTVkRkuJQ8q+QvpEkW0mJ2M2A+MNXMus1sPLAL8AwwDzgsa3soMLdg27pURYtI2yi54j4ZuNXMHgZOA85190WkwydzgTnAee6+GrgIOMbM5gH7AVcWaduoI5WkxAmNu2y5T5DZkiE3UvjFJnsGixVyI4XFK5c2btQk2khhaFbVehs3apLH/hJu5u87K18e8ly+p3c4MnfO+dBLP492tY5WTopI24h9YU1eStwi0jZqeqyriEhc9DxuEZHIaKikCV5c+nqZl19vkzFjg8QBOGLFk8FiVQM+Ff77Wx0ULNalK58KEufYng8GiQPw7UWPBIu19SZbBIt1wpaTgsVqBg2ViIhEplrrjBnOStwi0jY6ZKREiVtE2oeGSkREIqNZJSIikemQTd6VuEWkfSSo4hYRiUqfhkpEROKiinsAMxsBvBd4w907ZShJRCLSKYmp7mx1M7su+30S6a7GdwLPmNm+AfomIlJIQiX3EbNGy4x2yH6/GDjU3SeRbmR5Wam9EhEZhFqBI2Z514dW3f15AHfX/pMi0pKqVHIfMWs0xj3ezB4HNjazk4BbSDfIfKX0nomIFDS4HcniUzdxu/uHzawL2B1YSfoO42ngugB9ExEppBZ5JZ1Xw1kl7r6GdHfida4przsiIoOnh0yJiEQm9puOeSlxi0jbqFU0VCIiEpXqcHcgECVuEWkbmlUiIhIZzSppgvHdG5d5+f+nWgtzW+Kd3tVB4gDUAm5Z/d01HizWbV07BYnzA9YEiQNhN/AdO3JMsFh/6F0cLFYzaFZJREIlbRFpbRoqERGJTKeUcErcItI2qqq4RUTioopbRCQyStwiIpHpkC0nlbhFpH00u+LOdv+6zN2nmdkewPdJF2iuAY5z9z+b2SnAqUAfcJG7zzazCcAsYCzwOnCCu68s0rZev7Qhgoi0jWqBoxEzOwu4FujOTn0P+KK7TyPdxvErZjYROAOYDEwHLs0ehf1VYJa7TwWeAE4t0rZR35S4RaRt1Cr5jxxeBI7u9/oYd38y+3gUsBrYB5jn7mvcfSnwArAbMAW4N2t7D+mWj0Xa1lUocZvZBDPrkFEkEYlNM/ecdPc7gN5+r98AMLP9gS8A3wHGAUv7fdlyYPyA8xs616htXXXHuM3sBGAbYDbpGMxqYCMzO83d7290cRGRkMqeVWJm/wacBxzu7m+Z2TKgp1+THuBtYN35VRs4l6dtXY0q7tNI95j8FnCUu+8BTAMubXRhEZHQkgJHUWZ2LGmlPc3dF2an5wNTzazbzMYDuwDPAPOAw7I2hwJzC7atq1Hi7nX3d0jL94Wwfpf3TnmWi4hEpMlj3OuZ2UjgCtKK+E4ze8jMvu7ui7Lzc4E5wHnuvhq4CDjGzOYB+wFXFmnbqD+NpgPebWY/I/2pMNvMfgkckgUVEWkpzd5Iwd1fBvbNXm7+Lm1mAjMHnPszaa4cdNt66lbc7v5fwOVABXgV2BK4wt3PLhJERCSEGknuI2Z5dnl/GHg4QF9ERIZES95FRCITdx2dnxK3iLQNVdwiIpHpq3RGza3ELSJtozPSthK3iLQRDZU0wac327PMy6+3KAm387pVwu1c/9bfHpNQut6A3/K3BaoXJhLusTq/3rXh4yWaZufHXw0Wq2vc6GCxmiH2aX55qeIWkbbRGWlbiVtE2oiGSkREIlPtkJpbiVtE2oYqbhGRyCSquEVE4qKKW0QkMpoOKCISmc5I2w2ex21m40J1RERkqPpIch8xa7R12SIzOylIT0REhigp8CtmjRL3H4A9zWyOmR0YokMiIoNVK3DErNEY9yp3/4KZ7Q2cY2ZXAg8AC939ivK7JyKSX+yVdF6NEncFwN1/B/xztqX8AYCV3TERkaJir6TzapS4b+z/wt2XAj/PDhGRllJNVHHj7jeF6oiIyFBpHreISGQ0xi0iEhmNcYuIREZDJSIikdFQiYhIZDSrREQkMhoqaYLf975V5uXXe3b5n4LEAXh0VFewWEeN2zVYrLeStcFiTap2B4kT8kbV+OP2ChbrY6+Ee/bbY8sXBovVDLo5KSISGY1xi4hEpplDJWZ2DnAUMAa4GniYdDV5AjwDnO7uNTO7EDgc6ANmuPt8M9spb9vB9K3R0wFFRKKRJEnuox4zmwbsD0wGDgS2AS4Hznf3qaTPcfq4me2VfX4ScAxwVXaJIm0LU+IWkbZRJcl9NDAdeBq4i/TZTLOBD5NW3QD3AAcDU4D73D1x91eBUWb2noJtC9NQiYi0jSYOlUwAtgOOAHYA7gZGuPu6AMuB8cA4YEm/r1t3vlKgbeFZHErcItI2Gg2BFLAEeM7d1wJuZqtJh0vW6QHeBpZlHw88XyvQtjANlYhI26iR5D4aeAQ4xMwqZrY1sDHwQDb2DXAoMBeYB0w3sxFmti1pVb4YeKJA28JUcYtI22jWdEB3n21mBwDzSQvc04GXgJlmNgZ4Frjd3atmNhd4tF87gDMLtC2sUOLOOjHS3VcNNqCISFmaueTd3c/awOm/23vX3b8GfG3AuQV52w5G3cRtZjsDlwBrgSuAm0nvhJ7j7rcNNbiISDNpyXtqJvBN0jufs4HdSQfT7weUuEWkpXRK4m50c3KUu98P3AkscffX3P0doLf8romIFNOsBTitrlHF/bKZ3Zq1W2FmFwNLgTdK75mISEGdUnE3StzHA4cBC4AVwJeAlcCJJfdLRKQwPWQKcPc+0hVD65xZbndERAavmnTGg101j1tE2kbsY9d5KXGLSNvQGLeISGQ0xi0iEpmahkpEROKiiltEJDKaVdIEb/UuK/Py663qC7dD+Sc23y1YrDPGLg0Wa69XPFis+d2bBImz3dgtg8QB+MU3BvVY5UG5762ng8UaPWJksFjNoKESEZHIaKhERCQyqrhFRCKjiltEJDLVpDrcXQhCiVtE2oaWvIuIREZL3kVEIqOKW0QkMp0yq6TR1mXrmVmlzI6IiAxVUuBXzBrt8v4PwFXALsDWZvY4sBD4srsvCtA/EZHcOmXJe6OK+yrgDHffDpgKPAh8G7iu7I6JiBTVKZsFN0rc4919AYC7PwZMdvfHgc1K75mISEG1JMl9xKzRzcmFZnYNcA9wBPA7MzsceKf0nomIFBR7JZ1Xo4r7BOBp4GPAfOA/gSXAMSX3S0SksBpJ7iNmjXZ5X0s6zt3fY+V1R0Rk8Dql4tY8bhFpG50yq0SJW0TaRuw3HfNS4haRtqGhEhGRyDRrRaSZjQCuBnYH1gAnu/sLTbl4E+Re8i4i0uqauADnE0C3u+8HnE268LBlKHGLSNto4gKcKcC9sH7x4d5l972IUodKnnvzt3owVSRWDHcHRJqgb+1rzco544Cl/V5XzWyUu/c16fpDoopbROTvLQN6+r0e0SpJG5S4RUQ2ZB5wGICZ7Uu6grxlaFaJiMjfuwv4qJn9GqiQPv6jZVQ6Zd6jiEi70FCJiEhklLhFRCKjxC0iEpmWuTk5HEtMzWwScJm7TysxxmjgemB7oAu4yN3vLinWSGAmYEACfN7dnykjVhZvS+Bx4KPu/lyJcX5POj0L4CV3L+1GkZmdAxwFjAGudvdStukzs88Cn81edgN7ABPd/e0mxxkN3ET6/VcFTinr38rMuoAbgB1J/71Od/fny4jV6Vqp4g66xNTMzgKuJf1PU6ZjgSXuPhU4BLiyxFhHArj7ZOB84OKyAmUJ4YfAqrJiZHG6gYq7T8uOMpP2NGB/YDJwILBNWbHc/cZ1fybSH35nNDtpZw4DRrn7/sA3KPF7AjgFWOHu+wJfpNzv9Y7WSok79BLTF4GjS44B8BPgguzjClDaJH53/ynwuezldkAZiWCd/wauAV4vMQak78A2MrP7zGxONqe2LNNJ5+veBfwcmF1iLADMbG/gA+7+PyWFWACMyt7RjgN6S4oDsCvpNoe4uwO7lBiro7VS4t7gEtOygrn7HZT7Tbwuzgp3X25mPcDtpJVwmfH6zOwm4PvALWXEyN7mv+Xuvyzj+gOsJP0hMR34PHBLid8XE0gLhn/pF6vsxzacC3y9xOuvIB0meY50GO2KEmM9CRxhZpXsB+z7suE7abJWStwtvcR0KMxsG+BB4EfuPqvseO5+PLAzMNPMNi4hxImkixMeIh2bvdnMJpYQB9KK8cfunrj7AtI9T99bUqwlwC/dfW1WMa4G3lNSLMxsU8Dc/cGyYgBfIv0z7Uz67uWmbPipDNeT/j+eC3wSeNzdqyXF6mitlLhbeonpYJnZVsB9wFfc/fqSY30mu7kGaaVay46mcvcD3P3AbHz2SeA4d1/U7DiZE8nud5jZ1qTvzN4oKdYjwCFZxbg1sDFpMi/LAcADJV4f4K/87Z3sX4DRQFlV8D8BD7j7FNIhwoUlxel4LTOrhBZfYjoE5wKbAReY2bqx7kPdvYybencCN5jZr0j/g84oKU5I1wE3mtkjpDNlTizrnZi7zzazA4D5pEXN6SVXjEb5ye07wPVmNpd0psy57v5OSbGeB75pZueR3l85qaQ4HU9L3kVEItNKQyUiIpKDEreISGSUuEVEIqPELSISGSVuEZHIKHGLiERGiVtEJDL/B8IwouVl0TLzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10cd820b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(pb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10d5f6c18>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD3CAYAAADv7LToAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFoFJREFUeJzt3X20XXV54PHvzQtJgQR1iWIUaZHyVG1LrbBAXkK0IO/y4jh1bKsGoaVmpgPD2AGBBdTqtOPILCkgNoLRli5pqVTMFKUopDGAVESFteQJiJZB3ikQXszLPefMH/tEb6/ce8+5OXvn3N/9flh7kbPzu+f53eTmOc959u/s30in00GSNLPN2d4TkCRtO5O5JBXAZC5JBTCZS1IBTOaSVIB5dT75U+9c1shSmfbG5lbk7PTRDzUWq33jPzQWa9Ot9zUWa/6eL20kTmfTlkbiALQefb6xWB/75m6NxVo+/5nGYr3+3n8c2dbn2PLE/T0ng/kv33Ob4w0TK3NJKkCtlbkkNard2t4z2G5M5pLK0Rrd3jPYbkzmkorR6bS39xS2G5O5pHK0TeaSNPNZmUtSAbwAKkkFmMWVec/rzCPCNemShlqnNdrzUZpJK/OI2BO4CNgXGO0m9LuAMzJzfQPzk6TeeQF0Qp8Bzs7Mb249EREHAJ8FDqpzYpLUN9ssE1o4NpEDZOZtNc5Hkqav3er9KMxUlfl3I+JK4CvAM8Ai4Gjge3VPTJL6Nosr86mS+QeBE4CDgcXABmA1cG3N85Kk/hV4YbNXkybzzOxQJW6Tt6Th5wVQSZr5Op3yeuG9MplLKoc9c0kqgG0WSSqAlbkkFaDV3L6vw8ZkLqkctlnqMWdxM68VI/Oau4K96rgvNBbrP73t4cZi/fD2XRqL9brFzzYSZ8tDzVVpC355p8ZiXXjhHo3FeuyT32os1kDYZpGkAliZS1IBTOaSNPN1vAAqSQWwZy5JBbDNIkkFsDKXpAJYmUtSAazMJakAo25OIUkz3wAr84g4G3gHsANwGbAGWAV0gLuBFZnZjojzgWOAUeD0zLw9Ivbqdeyg5jvVhs6SNHO0270fk4iIZcCBwEHAocDuwEXAuZl5CDACHB8Rv9n9/f2BdwOXdp+in7EDMWllHhE3AQvGnR4BOpl54CAnIknbbHCV+RHAXVRbZi4GPgScSlWdA1wPvB1I4IbuFpsPRMS8iNgVeHOvYzPz8UFMeKo2y1nASuBEqrcFkjS8Brea5eXAHsCxwC8B1wFzuokY4FlgF6pE/+SYr9t6fqSPsfUn88z8ZkT8FfDrmemmzpKG2+Aq8yeBezJzM5ARsZGq1bLVIuBpYEP31+PPt/sYOxBT9swz8+Mmckkzwuho78fkvgEcGREjEbEE2An4WreXDnAUsBZYBxwREXMi4rVU1fsTwJ19jB0IV7NIKkenM/WYHmTm6ohYCtxOVfSuAH4IrIyIHYDvA9dkZisi1gK3jhkHcGYfYwfCZC6pHAP8BGhm/vGLnD70RcZdAFww7tz6XscOislcUjn8OL8kFcCP80tSAVrN7Qc8bGpN5ht/1NyuH69ed18jcR4+tLlX/o+sWdJYrAtPm99YrB+sbObP8MudVzYSB+CUzQ82FuvvvtRcrPeduLCxWANhm2VmayqRSxpyJnNJKoA9c0ma+Trtwawzn4lM5pLKYZtFkgrgahZJKoCVuSQVwGQuSQUY0I22ZiKTuaRyzOLKvO89QCNi/DZykjQc2p3ej8JMWJlHxHHAJcAW4JzMvLr7W9cDb2tgbpLUn1m8mmWyyvwc4DeodpL+g4h4X/f8SO2zkqRp6LTbPR+lmaxnvjkznwKIiOOBr0fEA0B5708klaHA9kmvJqvMfxQRF0XETpn5LHAScCnwK81MTZL61Gn3fhRmsmR+MvA9upV4Zv4/4K3A3zYwL0nqnxdAf15mjgKrxp17FDi95jlJ0vSMzt4LoK4zl1SOAtsnvTKZSypHge2TXpnMJRWjxCWHvTKZSyqHlbkkFcBkXo8FS+bW+fQ/9chv7dVIHIB5L2vmewI4b8ljjcVac/GujcXabYe+bwk0Lb+766ONxAFYsFsz3xPA0U8093MBOzQYawBm8cf5rcwlFcM9QCWpBCZzSSqAq1kkqQBW5pJUAJO5JM18ndZg2iwRMQe4DNgH2ASckpn3DeTJa9LceipJqtvg7pp4ArAwM98CnAV8ova5byOTuaRidNqdno8pHAx8BSAzbwP2rXvu26qvZB4Rv+CGzpKG1uAq88XAM2MetyJiqNvSk04uIt4AfAx4CrgK+AzVN/VfM3N1A/OTpN4NbmXiBmDRmMdzuns8DK2pXmkuB84DfhG4Btgb2AhcD5jMJQ2VzujAsvk64DjgbyPiAOCuQT1xXaZK5nMycw2wJiLempmPAUTEUL9CSZqlBleZXwscHhG3ACPA8oE9c02mSuYZEZ8Bfj8z3w8QEWcBj9Q9MUnq16DuzZKZbeC0gTxZQ6ZK5qcCx3W/sa0eBC6ub0qSNE2z99P8kyfzbhL/0rhzf13rjCRpmrxroiSVwMpckma+zixemmEyl1SMjpW5JBXAZC5JM5+VuSQVwGRek9aGZnbKnrvzSCNxAP5k3Ssai3Xefs19Nuugtz/eWKy5r9y5kTitx5tbpjb6b81lkd2W79FYrC3fvb+xWIPQaTWXC4aNlbmkYliZS1IBOm0rc0ma8azMJakAnY6VuSTNeFbmklSAtqtZJGnm8wKoJBVgNifzOb0OjIjmPi0jSdPQ6fR+lGbCyjwi9h536vMR8V6AzFxf66wkaRpmc2U+WZvlRuAF4CGqDU0D+DTQAd5W/9QkqT8uTXxx+wKXA5/KzH+KiJsy860NzUuS+taaxatZJuyZZ+ZjwH8EjomIDzc3JUmank5npOejNJNeAM3M0cw8narV0vPFUknaHjrtkZ6P0vS0NDEzVwGrap2JJG2jElep9Mp15pKKUWLF3SuTuaRitNqztxtsMpdUDNssklSAdoGrVHplMpdUjBKXHPbKZC6pGLZZajJ352YuRsx//ZJG4gCcs/HBxmLNfcmCxmJ1NrYaizV36SGNxFlw4EmNxAF44YxTG4u1ad29jcV68NuLGov1awN4DtssklSAOlezRMRc4CKqW50sAC7IzNURcQDwSWAUuCEzL4yIOcBlwD7AJuCUzLyvn7H9zm/2ruORVJxOH8c0/B4wPzMPAo4H9uqevxx4D3AwsH9EvAk4AViYmW8BzgI+MY2xfTGZSypGuzPS8zENRwA/joj/C6wEvhwRi4EFmfmDzOwAXwUOo0rWXwHIzNuAffsZO53J2WaRVIxBrWaJiA8AZ4w7/TiwETgWWAp8lqrK3jBmzLPAnsBi4Jkx51vdcz2NjYh5mTnaz5xN5pKK0R7Q82TmFcAVY89FxBeA1d2qek13A58NwNirxIuAp4Edx52f08/YfhP51gCSVIQOIz0f0/AN4GiAiNgHeCAzNwCbI+J1ETFC1YpZC6wbM/YA4K5+xk5nclbmkooxWu/SxJXApyLiNqrd107rnj8NuAqYS7VC5ZsR8S/A4RFxS3fs8mmM7YvJXFIxpllx9yQzNwEnv8j524ADxp1r87NkP62x/eo5mXfXQr4KeLgbXJKGymxOTJP2zCPiiu7/9wfWA18E7u72dSRpqNTcMx9qU10A/aXu/z8KHJWZ+1Oti/zzWmclSdPQ7uMoTa+rWVqZeS9AZrofqKSh1GKk56M0U/XMd4mIO4Cduovor6L6qOm/1j4zSerTLN41bvJknplvjogFVDeAeYHq3cldjFtML0nDoF1gxd2rKVezdJfj3D7m1OX1TUeSpm8W387cdeaSylHihc1emcwlFaM9YptFkma85vbLGj4mc0nFcDWLJBXA1Sw1aW9s5nLEpjsfZO7iZl6XmtqkGqCzubk3jRsf7Pv2ydP2zhX/1Eicv1j49UbiALxq3+YuvXXaza3ZWPKGDVMPGiKuZpnhmkrkkoabbRZJKoBLEyWpAC0rc0ma+azMJakAJnNJKkC9W4AON5O5pGJYmUtSAfw4vyQVwHXmPYqIlwNPZuZs/qCVpCFlm2UCEbEc2B1YDfwNsBHYMSI+mJk3NjA/SeqZyXxiHwSWAdcB78jM9RGxBPgSYDKXNFRmc8tgqrtGbcnM54FngfsBMvMhZvefmaQh1R7p/SjNVJX5dRHxJeBuYHVEfBU4EmjudnSS1KPZvJpl0so8M/8MuAgYAR4AXgFcnJlnNTA3SepLm07PR2mmXM2SmWuANQ3MRZK2iRdAJakA5dXbvTOZSyqGlbkkFWB0pP7aPCJOBN6Vme/pPv4t4E+BLcBjwHsz84WIOB84BhgFTs/M2yNiL2AV1ZuIu4EVmdnuZ+xE82puQ0tJqlmnj2M6IuKTwP/k3+fOy4ATMnMpcC9wSkT8JnAosD/wbuDS7tiLgHMz8xCqhSXH9zN2srmZzCUVo93HMU23AH847tyyzHy0++t5VJ+UPxi4ITM7mfkAMC8idgXezM8WlFwPHNbn2AnV2ma5+dZX1/n0PzW3kSiVpYc81Fis1obRxmL98z2vaSzWJYueaiTOHpcc20gcgL1++9KpBw3I3Yft2lis0Q0z65LioJYcRsQHgDPGnV6emVdHxLKxJzPz4e7XnAS8FTgP+O/Ak2OGPQvsAoyMubfV1nOL+xg7IXvmkooxqJeezLwCuKLX8RFxBvAfgCMzc2NEbAAWjRmyCHiaf/+mYOu5fsZOyDaLpGI00Gb5ORFxDnAIcFhmPtE9vQ44IiLmRMRrgTnd37tzTGV/FLC2z7ETsjKXVIxWwyvNI+KVwPnAt4HrIwLg6sz8VESsBW6lKppXdL/kTGBlROwAfB+4JjNbvY6dbC4mc0nFaGKdeWbeDNzc/fWjwA4TjLsAuGDcufVUK1emPXYiJnNJxejM4s+AmswlFcNPgEpSAUq8G2KvTOaSijF7U/kUSxMjYnFTE5GkbTVKp+ejNFOtM3+k+0koSRp6nT7+K81Uyfy7wJsi4usR0fMSGUnaHrbHh4aGxVQ9859k5n+OiH2BsyPiEuBrwP2ZeXH905Ok3pVYcfdqqmQ+ApCZ3wLeGRG7AEuBqHtiktSvEivuXk2VzFeNfZCZzwBf7h6SNFRaHSvzF5WZn2tqIpK0rVxnLkkFsGcuSQWwZy5JBbDNIkkFsM0iSQVwNYskFcA2S0323nFDnU//U6/Z//lG4gDM2XF+Y7Eeu6u519o37jjpXrED9enNzdy/7VdPvqWROAC7/8KujcV69Ls7NhbrNccvaCzWIHgBVJIKYM9ckgpgm0WSCtDxAqgkzXwtK3NJmvlss0hSAWyzSFIBrMwlqQAuTexRROwAzM3Mn9Q0H0maNj/OP4GI2Bv4GLAZuBj4PDAvIs7OzKsbmJ8k9cw2y8RWAh8BdgFWA/sATwM3AiZzSUNlNifzOVP8/rzMvBH4IvBkZv44M58HttQ/NUnqT6fT6fkozVSV+Y8i4gvdcc9FxEeBZ4CHa5+ZJPVpNlfmUyXz9wFHA+uB54AzgBeAk2uelyT1rYnVLBFxIvCuzHxP9/FhwJ8Bo8CNmXlu9/z5wDHd86dn5u0RsRewCugAdwMrMrPdz9iJ5jVpMs/MUeC6MafO7Pcbl6SmtDr13gQ3Ij4JHAF8Z8zpjwO/A3wfWBsRvwbMBw4F9gd2B/4e2A+4CDg3M2+OiMuB4yPiX3sdC1w70dym6plL0ozRQM/8FuAPx527E3gZVQJfCLSAg4EbMrOTmQ9QrQLcFXgzsKb7ddcDh/U5dkJ+aEhSMQbVM4+ID1C1lcdanplXR8Sycefvolrt9yTwPeAe4KTu462epVoVOJKZnXHnFvcxdkImc0nFGFTPPDOvAK6YalxEvAQ4G3hjZv44Iv4XVTt6A7BozNBFVMu62y9yrp+xE7LNIqkY7U6n52NAfkK1OOS57uOHgZcC64AjImJORLwWmJOZTwB3jqnsjwLW9jl2QlbmkorR9L1ZMnNTRJwJ3BARG6mq5/dn5lMRsRa4lapoXtH9kjOBld1bo3wfuCYzW72OnWwuI3Uunr97z2Mb+ZNtdkPn5l7/Hl7XXKzWlubepK1sLZp60AD86mhzf34r2w82FmvVzmVu6Lzzx68d2dbn+JVX7NdzzrnnsX/Z5njDpNaf9t0PauZ+XJ12c38nz+doY7EeeOpljcXab9ljjcU6f14zPxebH2k1Egfgt9/0ksZitR5/bupBA7LlRy80FmsQBtg+mXFss0gqhrfAlaQCWJlLUgGszCWpAK1Oc9dJho3JXFIxSry1ba9M5pKK4S1wJakAVuaSVIDZvJql54/9RURRn5aSVJ5OH/+VZtLKPCJeB1wKvB5YEhF3APcD/y0zH2lgfpLUs7o3pxhmU1XmlwJ/lJl7AIcANwGfoIdbQ0pS02bzhs5TJfNdMnM9QGbeBhyUmXdQ3eJRkobKdrgF7tCY6gLo/d29564HjgW+FRHHAM3dplCSelRixd2rqSrz5VRbIr0duB34ENX2Ru+ueV6S1Lc2nZ6P0kxamWfmZqq++Vi31TcdSZq+2VyZu85cUjFm82oWk7mkYpR4YbNXJnNJxbDNIkkFKPGTnb0ymUsqhpW5JBVgNvfMR2bzK5kklaLnuyZKkoaXyVySCmAyl6QCmMwlqQAmc0kqgMlckgpgMpekAgzNh4YiYg5wGbAPsAk4JTPvqznm/sCfZ+ayGmPMB64EfhFYAPxpZl5XU6y5wEoggA5wWmbeXUesbrxXAHcAh2fmPTXG+Tawofvwh5m5vMZYZwPvAHYALsvMWrZIjIj3A+/vPlwI/AawW2Y+PeA484HPUf38tYBT6/q7iogFwGeBPan+vlZk5r11xNLPG6bK/ARgYWa+BTiLaq/R2kTEHwOfofqHVKffBZ7MzEOAI4FLaox1HEBmHgScC3y0rkDdJPFp4Cd1xejGWQiMZOay7lFnIl8GHAgcBBwK7F5XrMxctfV7onpB/KNBJ/Kuo4F5mXkg8CfU+DMBnAo8l5kHAP+Fen/WNc4wJfODga/AT/cb3bfmeD8ATqo5BsDfAed1fz0CjNYVKDP/Afj97sM9gDqSw1b/G7gceKjGGFC9U9sxIm6IiK9HxAE1xjqCameta4EvA6trjAVAROwLvDEz/7KmEOuBed13vouBLTXFAXgD1RaTZGYCr68xlsYZpmS+GHhmzONWRNTWBsrMv6feH+ytcZ7LzGcjYhFwDVXFXGe80Yj4HPAXwFV1xOi2CB7PzK/W8fzjvED1wnEEcBpwVY0/Fy+nKiLeNSbWSE2xtvowcGGNz/8cVYvlHqoW3MU1xvoOcGxEjHRfdF/dbf2pAcOUzDcAi8Y8npOZtVWxTYqI3YGbgL/KzL+pO15mvg/YG1gZETvVEOJk4PCIuJmq1/v5iNithjhQVZZ/nZmdzFxPtQftq2qK9STw1czc3K0sNwK71hSLiHgJEJl5U10xgDOovqe9qd7lfK7buqrDlVT/jtcCJwJ3ZGarplgaZ5iS+Tqq/h7dV/W7tu90BiMiXgncAPyPzLyy5li/172AB1VF2+4eA5WZSzPz0G6/9zvAezPzkUHH6TqZ7vWTiFhC9Q7u4ZpifQM4sltZLgF2okrwdVkKfK3G5wd4ip+94/03YD5QV7W8H/C1zDyYqr14f01x9CKGZjULVZ/y8Ii4haq3XNuFroZ9GHgpcF5EbO2dH5WZdVw4/CLw2Yj4Z6p/tKfXFKdJVwCrIuIbVCt0Tq7rHVtmro6IpcDtVIXOipory6D+hPd/gCsjYi3VCp0PZ+bzNcW6F/hIRJxDdb3mAzXF0YvwFriSVIBharNIkqbJZC5JBTCZS1IBTOaSVACTuSQVwGQuSQUwmUtSAf4/iph0cSBeXwgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d474978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(pa - pb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-23424.349999999999"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(pa - pb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "94580.070000000007"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(pa + pb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9458007, 945800.69999999995)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_population = np.sum(pa + pb)\n",
    "avg_district_size = total_population / 20\n",
    "(total_population, avg_district_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(286356.4212655271, 3285084.9639999997)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# maximum standard deviation of district populations\n",
    "max_std = (164e10 / 20) ** 0.5\n",
    "# minimum value of <total population of b's won districts> - <pop. of a's won districts>\n",
    "min_victory_size_gap = ((-0.074 * total_population) + np.sum(pb) - np.sum(pa)) * 2\n",
    "(max_std, min_victory_size_gap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 1, 2, 3, 4],\n",
       " [5, 6, 7, 8, 9],\n",
       " [10, 11, 12, 13, 14],\n",
       " [15, 16, 17, 18, 19],\n",
       " [20, 21, 22, 23, 24],\n",
       " [25, 26, 27, 28, 29],\n",
       " [30, 31, 32, 33, 34],\n",
       " [35, 36, 37, 38, 39],\n",
       " [40, 41, 42, 43, 44],\n",
       " [45, 46, 47, 48, 49],\n",
       " [50, 51, 52, 53, 54],\n",
       " [55, 56, 57, 58, 59],\n",
       " [60, 61, 62, 63, 64],\n",
       " [65, 66, 67, 68, 69],\n",
       " [70, 71, 72, 73, 74],\n",
       " [75, 76, 77, 78, 79],\n",
       " [80, 81, 82, 83, 84],\n",
       " [85, 86, 87, 88, 89],\n",
       " [90, 91, 92, 93, 94],\n",
       " [95, 96, 97, 98, 99]]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We need to give big districts to B and small districts to A; the efficiency gap metric doesn't really hurt us?\n",
    "# Careful about the imbalance though\n",
    "\n",
    "# Small naive test\n",
    "naive = [list(range(x, x + 5)) for x in range(0, 100, 5)]\n",
    "naive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's write the grader, at least.\n",
    "from scipy.stats import norm\n",
    "\n",
    "population_a = pa.ravel()\n",
    "population_b = pb.ravel()\n",
    "population = population_a + population_b\n",
    "\n",
    "def win_probability(na, nb):\n",
    "    mean = .6 * na + .4 * nb\n",
    "    variance = (na + nb) * .4 * .6\n",
    "    z = ((na + nb) / 2.0 - mean) / (variance ** 0.5)\n",
    "    return 1 - norm.cdf(z)\n",
    "\n",
    "def expected_districts_won(soln):\n",
    "    ret = 0\n",
    "    for district in soln:\n",
    "        na = np.sum(population_a[district])\n",
    "        nb = np.sum(population_b[district])\n",
    "        ret += win_probability(na, nb)\n",
    "    return ret\n",
    "\n",
    "win_probability(50e3, 50e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.73965876240101747"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "win_probability(51e3, 50e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.89942375697270061"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "win_probability(52e3, 50e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99953573328876688"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "win_probability(50e3, 45e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4191282431929132"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "win_probability(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0000518963627822"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expected_districts_won(naive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "655328510334.54993"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def district_population_imbalance(soln):\n",
    "    sizes = [np.sum(population[district]) for district in soln]\n",
    "    return np.var(sizes) * len(sizes)\n",
    "\n",
    "district_population_imbalance(naive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.20981211954077811"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def expected_efficiency_gap(soln):\n",
    "    # This is approximate only\n",
    "    # Who tf wants to integrate this actual mess :(\n",
    "    total = 0.2 * sum(population_a) - 0.2 * sum(population_b)\n",
    "    for district in soln:\n",
    "        na = np.sum(population_a[district])\n",
    "        nb = np.sum(population_b[district])\n",
    "        p = win_probability(na, nb)\n",
    "        ea = 0.6 * na + 0.4 * nb\n",
    "        eb = 0.6 * nb + 0.4 * na\n",
    "        used = p * eb - (1 - p) * ea # used for a - used for b\n",
    "        total -= used\n",
    "    return total / sum(population)\n",
    "\n",
    "expected_efficiency_gap(naive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10dab0940>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD3CAYAAAC+eIeLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0XGed//H3FPVR79W2JOu6yj3uJQkxCUlID2SBFBZY2Ozuj8AuZ2GpC7vZ5bDwY38QOgQInfSQhDTHcW9yL9eWi7pk9d6m/P4YaayRZFuOR1eW/Xmd43M8M/fOfHVn5jPPfe5zn2vz+XyIiIg17BNdgIjItUShKyJiIYWuiIiFFLoiIhZS6IqIWMh5oQfr69uvmKENiYnRNDd3TXQZl2Qy1gyTs+7JWDNMzronY81gbd2pqbG28z02aVq6Tqdjoku4ZJOxZpicdU/GmmFy1j0Za4Yrp+5JE7oiIlcDha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIhha6IiIUUuiIiFlLoiohYSKErImIh50QXIHKt8Xp9PLvpFFsP1ZKcEMVdq6Yxc0riRJclFlFLV8Rib5VU8pdtZTS391Ja0cL/Pn2Arh73RJclFlHoiljsaFlz0O3ePg+nalonqBqxmkJXxGJT0mODbjvsNnJSXRNUjVhNfboiFuh3e3lu0ykOnGwkLSmauflJHDrdRGx0OPetKyDBFTHRJYpFFLoiFnh20yle3VEOQFVDJzmpLp54bC2ZGXE0NXVOcHViJXUviFhgf2lD0O3K+g46uvtxOPQVvNaopStigayUGGoauwK3XVFhxMWE89buct7cWU5yXCS3Lp9CUlzkBFYpVlDoiljgvusLqW3soqqhk5hIJw/dPIPtR2r5xcvHAsscLWvmGx9fit1mm8BKZbxp30bEAnYbFBcks2JOBo/dP59FRio7j9QFLVPb1EVZbfsEVShWUUtXZJx19vTzjV/upq2rH4Bdx87ylYeXkDisK8Fus2kUwzVALV2Rcbb3eEMgcME/fGzb4VpuWzGVtKRoAGw2eP/KqSTGKnSvdmrpioyzmMiRX7PIcAf7jtczIy+RRdOTWbcwh9T4qAmoTqym0BUZZ3MLkinKTeB4RQsA6UnRVJztYOfRs4Fl7HY796wtmKgSxULqXhAZZ06HnYduMVi/JJcP3FjIlx9azB6zPmiZzQdqJqg6sZpauiLjbPexs/zw+cN4fT4A3G4v0ZFO2of087qiwyaqPLGYWroi4+yFLacDgQvw0rYy7l6Tz+BwXKfDpq6Fa4hauiLjrK/fG3Tb7fZSkB3PbcumEB4ZxopZ6Rq1cA1R6IqMsxsWZvP7t0oDt+fkJ/HvT+7G7fGHcXlNG5+6c85ElScWU+iKjLP11+UR7wrn4KkmZuQlUnK8PhC44D9Z4p7mLtISoyewSrGKQldknL22s5yn3zlFv9tLbVMXkeGOEct4vL5R1pSrkQ6kiYyjprYe/rjhJP1uf8v2VHUbUeFOHPZzk9oUFySTmRwzUSWKxdTSFRlH9S3dQSMXAPo9Xr788BJKjteTn5vAzJz4CapOJoJCV2Qc5WfFER8TTmtnX+C+BdNTyE1zkZvmIjU1lvp6zSx2LVHoioyjMKeDz35wPs9tOk1zew9LZ2Wwdn72RJclE0ihKzKOjpU18+vXTOqaullQlMKquZnUt3Tz9MaT1DZ2sXJ+NjfMz8Rh1+GVa4VCV2Sc9Lu9/OD5Q4HTffeY9cRHh3O0vDlw6Z7y10y6uvq4Y9W0iSxVLKSfV5Fxcra5K2h+BYBjFS1B10oD2HsiePIbubopdEXGSXpSNHEx4UH3FeXGjxinq+Fi1xaFrsg4cTrsPHrXHHJSXYQ77Sydlc596wr5yHuNQPDmZcRyz5r8Ca5UrKQ+XZFxVF7XQVSEg4VFqdyxehpREU5mTU3ipsU5tHb18cH1M4lU0+eaotCVq55Z3szxihYKsuOZNTXJstfdsLeK37x+HIATtFJa1cpXH1nCN365m8a2HgB2HDnLFx9cTHaKuhiuFQpduar9dWc5fxgyw9c9a/O5dfnUcXmtmsZOXtleTmdPP6vnZVFing16vKG1hzdLKgOBC9Db52HLgRruv6FwXGqSK49CV65qr+woD769vTwkoVvV0Mlbeyrx+XysW5BNakIU//WbksBohX0nGiguTA5ax2G3kTzssusA4WHqX7iWKHRFLlFzey//+evddPd6ANh6uJb71hUGDQ/zARFhDrJTY6iq78TpsHHv2gKWzkrn7X3VlFa2ApASH6kz1K4xCl25qr1vaV7QBOK3LMu77OfcbZ4NBC74rwxR1dAxYrmMpGj+7v2zqWroJMEVgSsqjK6efpbPSidvYO6F960uoLO9Z8S6cvVS6MpVq7yundauPm5YmENMpJMZeQnMDMGBtNhRLiKZlx7LmnlZvLO/GoCcVBfvWZyLzWYjJ9UFQEd3P1/7xU4a23oBSEuI4pZVujbatUahK1elE5UtfPO3ewOTg0/JiOWO1aE51XZRURrTc6o4MdBFMCUjlvkFKfT1eVhTnMn03ASWz8mgp9fNC5tPc7alm0VGKg2tPYHABTjb0s2mfVUsGtb3K1c3ha5clTbsrQq6GkNZbTulla0U5SZc9nM3tHZTlJtATqqL2flJFBck841f7qa8zt/FsOVQLWmJUfxxQyknq9oA2HqolqWz0kc816VcMaLkeD1natsw8hKZbeHQNwktha5clcKdI0cEhGKUQOXZDr7xq930DVwJ4mhZM1HhjkDggj9IX95eFgjcQY2t3UFz6ybGRrBmQTY9nb1czJ82lAZGYry0tYwHbpzOTUtyL/vvEespdOWqU3G2g/mFKew6Vk93rxuAeQXJTM2Iu+zn3ri/OhC4ALVNXVTUjTyIFhXhxGaDoReNiHdF8OjdxWw7VIvNBsvnZBAbHX7R0PV4vby5pzLovtd2lSt0JymFrlw1+vo9fPfPBzha1gxAcX4yC4pSSHBFMDc/NP2mYaO0oLNTXcwvTGFfaQMA0RFObls+lURXRKB1GhPp5IaF2Ty/+TRlte3MmppI1CgXqByNDRv2IddUAzT/7iSm0J2kqho66e51k58Vh91mu/gKVxmfz8fWQ7WY5S3kZ8Vx941FbDtcGwhcgAOnGlm7IIt5hSkhe90bFmSz5WBNYEzu9Jx4uvvcJLjCuXFhNtlpLhYWpdLT58HhsLF2XhZTM2NZMiOdJ547yJEz/vpO17TR3evmsQ8tPu9rdXT3887+ajq6+1kxJ5O3Ss61dm9dPuVd1d/X72HjvmpqmrqYX5hCcYEO4llNoTsJ/eTFI2w7XAtAbpqLf3lgAa6okcOYLpd/t7aKExUt5GfHcdPiXJyOsbWw3B4vv3vzBDuP1JEYG8kDNxYGDdeqb+kmMtxBbHT4BZ7l/J555xR/2VYGwOaDNdS39476YW5qu3h/6aVo7epj3fwsevs8TM9NoL6lmyeePRR4fMWcDIrzk/n3X+yia6Bro+REGDOnJgYCd1DJ8fPPo+v2eHn8qT2BuXedDhsP3WzQ0+dhRl4iUzJiL1prv9vL9sO1NLT2sMhIJS89lieeO8SBk40AvL23ir+9dSYr52Ze8naQd0+hO8kcr2gJBC74+y837K3i9hVTQ/5af3izlDcG+hL3HK+ntrGLR943c0zrvrKjnA0lVQB09nTwvWcP8j+PrsTng+89c5CjZc047DZuXprHPWsvfazq4HjYQW/sLOcLH17Ey9vLAiMCwsPsGLkJvLmnkp4+N8tnZ5A0ymm4mw5U89ym0/T0ubl+QQ73rM3HNsrewzv7q3nylWOB22FhDnabwcG5/XAdqQlRgcAFaO/q5+DJRhJc4bR0nLtAZVpi9Hn/vqNlzUGTnbs9Ps7UtvPQzTPOu85w//v0AQ6fbgLg5e1lfPTWmYHAHfT23iqFrsXUMTTJtHSMbLm1tIe2NTdo66HaoNvbDtfi841tiNPx8uBWXXevh/K6Dt4qqQx0AXi8Pv6yrYyKsyMPRF1MdERweyE60kluWgwPv28G8wqTuW5mGo/dN58nnjvEb14/ztMbT/Hln+3kbEt30Ho1jZ08+fIxmtt76e718PL2MnYcqRv1NV/ZXhZ0+/VdFUQMGxER5rSP2lcbGe7kwffOIGqg7sTYCB64cfp5/77RRl9EhI2tDxj83U+DgQv+bb3tUC2OYX3DkRFqd1lNoTvJzJmWHNSVYINRx3+GwvCrHsRGh4/aAhzNtKzgkQLhYXayU2Ooa+oesWxdU9eI+y7mnrUFgQCxAXeuLeDLP9/Fz146yrHyFubmJ9PW1UftkOfu6nWzaVgLubSqleE/IyeqWsdcx+0rpgYF2W0rprCqOJP0pHOt2MFJzN/cU0FhVhx/e+tMvvmp5RfsIijKTWD21MTA7djoMG5YlENzey9lte0X/fEbHq4AEeEO1l93bsRDuNM+LntIcmHj8jN3uqaNstp2ZkxJJCPp/LtQY9Xv9rLlQDWNTZ3ML0wJtBZCqaqhkxMVLUzLjBv1y9Dv9rKvtAGPx8v86SlEhlvfQvD6fJyqbuWmxTlUN3bh8XhZMz9rTAP+e/rc7DvRgMNhZ35hMmHO87eaaho62byvitXFmTy76TRujxe7zd9f+c7+auYVJBPvihh1XZ/Px5EzzcTHRLBgegr7SxuJd4Vzw8Jsdh07S35WHJsP1gSWjwiz4/X52LS/mnnTU4gbYx/v3PxkPnBDIWdbulk7L4s3SqqpbugE/NMlPvX6cR5cb4xYb/Cg45naNs7UtJMUF4ENgoK3MDs+aJ1T1W2U17WzYk4Gz246Hbj/xsU5FOUmcPfafFo6elk+K4OpmXG4PV7uWj0tcJAvNSGK//5tSWD42LGKFoy8BFLiozhe3sz+Y3WjflfWX5dHRnIMaYlRrJiTwas7ynl5exk+nz/I//mD80f8MHq9Pg6caqSjq58F01PYe8I/oiLMaee91+URGx2GzWbDabdx/cIc4mPGtr3PtnRz5EwTOSkuUlMv3p88qK6pi6PlzeSmuSjI8m9Xt8fLgZONdPe6WTA9lejIS/su9fV72FfagM8H86enXNIewFi0dfaxv7SBhNgIZk9LCvmB6pAnx0tbz/DMO6cA/wf8k3fMZvGMtHf9fP1uL//51B7KatsBSI6L5EsPLR7xYbscWw/V8LOXjga+eMMHnvf1e/iPX+8J7AanJkTypYeWjMvBqwv5wXOH2DPQhxgV4eDzH1pETprrout1dPfz9V/uor7FP7FKTqqLf3tw0agf1n2lDXz/mYOBftGbluQye2oiG/dVBw5cRYY7+NcPLSQvfeSX78cvHgnsnkeEO/i3jyxiz/F6nt7o/0w47DbWL8nlVHUb0ZFOunvd/PD5wwBEv+Xk8x9ZdNEJvTt7+vn6k7sDXQXHypqJGPYj2NvnIScthtw0V+B9i40OY+38LF7dUc4fN/gnwbHZYM38LPaXNtDT5+H6hdksG7Ln8OLWMzw75PN895p8vD4fU9JjSUuM4gs/3k5nj7//tq6pm3+6p5hv/m5vYBaxkhP1LJqeGjRet999LnQGt8vw78pPXzoa6LuPCHOQ6IoIbH+AyvoOXt1RHjQPr8/n4zt/2h/oVoiJdPLgew3cHi/zC1OoqO/gv39TEnQW3J2rL36poAMnG/h/T5/7TNxf287Ni3Muut4es54fPHcI78Aff+fqady6fArf/O1eSgf2JhJcJ/nSQ0tIjB39R3y47l43//HrPYEf2IykaL744OJLDu7zqazv4PGnSgLjuxcVpfLo3XND8tyDQtq90O/2Bn0wvD4fL2w5c1nPua+0IRC4AI1tPWw6UH2BNS7dC5vPBLV0XthyOvBBAf9BpKH9jvUtPWwZ0lqzQuXZjkDggr+P9K+7yi+wxjlbD9YEAhf8H6wSc/Qj5y9uORP0pdxQUkVsdHigxQTQ0+fh1Z0jX7umsTOoP7S3z8MrO8p4bUidHq+PUzVtfOEji7hrdX5g/gLw7/6/vqtiDH9PbVDfbGV9Z9DuPEBGUhRmeQvTsmJZvySXD68v4usfW0q8K5wXt55rrfp8/oOT3/6HVTzxmbXct64w0IXS7/by8rDP865jZ3n/ymnMK0zhjd2VgcAFOHCykTf2VAYCF6C1o4+GtpFdKkmxEby0dfh3xV9XXXNX0MHS3n5P4IDmUA2twc97orI1qB+3s8dNdUMn71mcS0pC1Ij39pUd5fT2ebiY4es9+3ZpIJQuvF7w9+jl7WXsL20IBC5AS0cfb++tuuhzDdp17GwgcMF/csqOI7UXWOPSvLazIuhvG/7dD4WQtnS9Ph9ujzfovj73xd/UC+nrH7l+X793lCXfvd5hNfa7vf79TduQ22OoazwNr9Ffw9i2Q98o9Y/2fAD9w+73eL1jfg9Gu6+334PH4xu2nCfw2MhaL75dR1umeHoKaQmR7D3RQEZSNO1dffz2jROBxx+5ZQZx0eG4PV763SPr6ejup6axk7z02MAegNc7yue53zPq/wf1jBJGaQnRFBfYOHCyERuwcm4mM6YkjvLc/tv9o2zHCKedeFc4rUNGPyyZGdyXP+rndMi2Gv64x+Mb09wPwz8/Ho/3Xa7no2eUkL+UjBhtm/eGMA9GqyXU3/WQtnQjwhysKg4efnLjwovvhlzIgumpQbsekeEOVs7NuKznHG54jesWZAedAbSwKJV417nujKgIJ8vnhLaGi8nPjGNa5rmDU3abjQVFqRwrax71yzbU8tkZQUf7413hLDZG7/K5YVHwtlg2KwMjLzGon9Nus3H9gpETb0/JiKUw59xyNpu/e+K6YQf6Brd3YU48U4Z0UTjsNtaNYULv5bMziBmyOxkXHUbl2Q52H6snMzma9y7JHTEmdsNAa8rpsLN2XlbQY9NzEvjs97fw+FMl/PP3t3CisgXwd48MH041tyCZH794mJ//5SgzpyYGHbDKSXWxfkkuqQnnhqWFO+2smZ/Fp++bxxc+spAVczJo6exlx5E6Vg/7rgxu+5w0V1A/vc0GNy7O5XMPLGD57AxmT03k47fNYsmwbrsZUxKCumacDlvQBOnD39vls9PHtFt+47D11izMGVPX2vD1VszJYJGRRkr8kO0TZmdVcdbwVc9rycz0oKk1XVFhIT2QvG5+dlAf7rTMWPKzLv/08aFsFzoKWl/fPvYpkAZ4vT62Ha6lvK6D2dMSKS64/LOBWjt62XOykZbWblbOzST9AuMb362S4/X+3dHMWJbOSh9xlL65vZdNB6rxeHysKs4kNSHqos+ZmhpLfX37RZcbq+5eN5sO1NDS3kuv28Pbe6vw+fwh+i8fXEDWBfpC61u62XygBofDxurirAv2oVU0drNlXyXZqTGsmJOBw26np8//2s1tvSyZmRb0AzBUT5+bzQdqaGrrZfGMNPKz/AeWth6qpbqhk7kFyUEzZHUPjCho6exj6cz0MQ36B2ho6WbTwN/T3NbLxiGjEnLTXFQ3dAa1xmZOSeRfHlgA+PfIdhyp40xNO0V58fz6VZO2IVd9KMiK498e9J8pNvTznJ4Uxe/fPIF7oOUeGe7g0bvmcuRME339Xux2G9mpMcydlsTWw7X09HlYMSeDzGT/+/KVn+8M2lV94D3TyUh1cehEAx6vl+qGTpLiInn/yqnEx0Sw+WANFWfbyUqJYXVx1pgOIHf29LNpfw3t3X0sn5Uxos//0KlGDp5qIjfNxfI56WM+nfjQ6UYOnWoiOyWG918/neamzouvhL/L5fDpJvLSXSyfnYHdbqOtq49N+6tHbJ+xamrrYdOBGnw+H6uLs0iOHzn2ejRj/T6eqW1j59GzJLgiWF2c+a4O3Kemxp736FvIQ3e8hDrArDBeNbd19fHZ720JCpXrZqbxyTvmhOT5J9u2/vyPtlHXHNy/eeOinMAkMWFOO5++bx4zpySOWLe3z8Pff3tjUJ9+giucb//DqhHLPvPOKV7aeibovo+sLyLBFcH3njkYeI55Bcn8n/vmsf1ILQdPNpGTGsPMqYn8+5O7g9YtzI7nO59ZxzNvmPxiyEkXKfGRPP53y3irpIrfv3kCn88/LvkzH5gf8lbXuzHZPh+DrKz7QqGrkdGTUHtX/4g+tfE6QWIyyEiKDgrd2Ogw7ltXwNKZadQ1dzNratJ5W/YR4Q7mDZmsBuC6maPvria4Ro6YiXdF8MaeyqDQ3n+ycURAz5mWhNNhC7SSAZLi/DXtGuWqwWZ5C09vPBkY9dDV6+bZd07y2Q8uGLU2mTx0csQklJ0SQ1568G7jstnW9jFfSe6/oTDQtRIT6WReQQqPfW8L3/zdPs7UtAf1x4/m47fP4pZlecyamsi96wq4d93opyWvnJNJwZCWZnFBMvMKk0fMPGYDdh8LDtJDp5t43/Ipgf7CxNgI7ljlv5JFanxwV5XdZiMmMmzEgcnmjl5e3HqGz/94O//11B6OV7Rc8O+SK5Pjq1/96nkf7OrqO/+DFouJiaCrq+/iC15BxrPmhUWpeH0+YiLDuGXZlJBeUXaybevY6HDuXz+D4mmJLCxK5TevH8ft8eL1+Thd00ZaYhS5aefvKw5z2pk9NYkVczKZnpMwYhrFQU6HnaUz03DYbRTmxHP/9YVEhjtJjI1gx5GzgeFRq4szcXt8NLadG6YXHmbnH+8uZmFRCjY7zJ2WTGFOPIkJ0SS7wjl4qpGO7n4cdht3rZnGkpnpnKhsCRrqV5SbwFslVXR299PY1stu8yzXL8gZdbrJ8TTZPh+DrKw7Jibia+d7TN0Lk5QrKox+t5fDZ5owK1to7+rj1uVTJ7qsCWOz2UhPjMYsHzmG+3RNOyvmXP6kLv4TdUoCB8M2H6jhyw8vwchL5GsfXcKbeyrJTXOxZl4Wxyta+M6f9gdaq3esmkZXr5tv/3E/Hd3+g3ZvlVTy/c/dSGJsBF//2FIqz3YQ74oInCX26F1zeXVHOZVnO5g3PYWDp4Inq+np83CyujVkcwWLNRS6k9T2I7WBYVCePg9PbzzFrKlJ5x1VcK0oyk0YccWG0Q6gvRv7SxuCRh+0dvbxzv5qls9O51u/30fzQL96dUMXD7xnOt/85HLe2O0/mLeoKJWth2oCgQvQ2NbLlgPVLMj3n2o6/Ay/uuYu9hyvp7qhk/bufqZlBj9ut9nIusQj/zLx1Kc7SZXVjjxLZuiZe9eqjKRoPn77LDKSokmMjeDGRdm4PV7aQrBbOdoJAR6vl5e3lwcCF+CN3RWcbenmzxtP8dK2Ml7aVsYXf7qDhtaeEes7Hec/r/8nLx4JnH1VWtVKdWMX8wYmHY+KcPChm6aPebiUXDkUupPUzKnBrTebDYy8y7/S7dVg2awM/vMTy1hdnMmbe6r44fOH+dwPtl72gaf5hSmkDRmfHRXhZFVxFm2dwYHuA6rqO9hy4Nyp4m6Pj7qmrsCIBYDM5GhWzB39xIDuXnfQfLoAFXXt/OO9xdx3fQFGbgKtnX309F38dFy5sqh7YZKaX5jCB28o5I09lUSEObh95dRLHmR+Nevs6eflIfPf9vV7eXHL6csachUR7uCLDy1m66Fa6lu6iY5w0tDSzYo5GUFDzjKTo8lJdY2YMtJms/G1j17HWyWVdPe4uXlpHpERTkbbP4mKcDIlPZayunOPGnmJvLD5dGA+k32ljZTXdfBP9xa/679JrKfQncTWX5fH+uvyJrqMK1JfvzdoTCwQNDnNu+WKCiMpNoI/vOU/aeHFrXDz0jz+8e657DhaR3iYnTXFWaQmRLGwKDVwSR6bDd6zKIc3d1fy/ObT+ICth+t4/O9XEnWeLoZP3jmbX71qUl7XzswpiXx4fRHf+FXwCRb7Sxvo6uknOtLaGe/k3VPoylUpMTaC4oLkoMvTrJ0/9nP8L+QvA3PaDnpjdwXvWzaFrYdq2Xy0ls0HallYlMrHbpvJ7mMp1DV3sWB6KsnxkTzx3KFAC7its4/fv36cR24eOe8vQHpidODU5UHxroigYWRREU7CQzyfrIwvha5ctT515xw2lFRR0+if/H5BUWpInnf4rGleL+w70cCeIReaLDlez/7SNBYZqZgVLcREOunoHnkmYVPbyINrF3Lv2gK+++f9dPd6cNht3H9D4ZgvFipXBoWuXLUiwhzcvDT03S/rl+Ty85ePBm6vKs4cdXTE8YpmfvVXMzA/691r8pmSERs0yuSGxbkj1ruQotwEvvX3KzlZ3Up2imvMk3/LlUOhK3KJVhVnkpYYxaHTjeSkulg8I426pi6efedUoCXrsNuoaugKmhD7hS1n+MbHrmPjvmpqm7ow8hJYv3TKJU/CEhXhZM40nRAxWSl0Rd6FotyEoDlvM5Nj+MwH5vP6rgp8Ph/rl+Ty540ng9Zxe/zTP+akunh7XxV7TzSw4+hZHr1zzqiXhperkzqDREKks7uf6sZOSqtaOXS6acSpxzPyEoiJDBvocvBfjeB0dVvgmoJybVBLVyQEmtp6+NELhwPdC6/sKOfhW2bwd++fzR7zLCkJUdy+YipN7b0jLlM0/CQIubqppSsSAqVVrSNGJhyvaMHj9XK0rJlXd5TzvWcOEhsdRlpi8FSO8wvVP3stUUtXJASmZcaNmGgnKyWGJ18xAxegPFrWzMvbynjsvnk8vfEktU3drF6QzQ3zL38GNJk81NIVCYHUhCgeuWUmcTHhOB12rl+QTUF23Igr/lbWd+AYuIJEZ08/tY2dY7oMulw91NIVCZFVxZmsKs7E6/Nht9nod3uIjQ6jfchFL+dMS+ZHzx/mZHUbABv2VNLd088nbp89UWWLxdTSFQmxwUvyhDnab7OQAAAJj0lEQVQdPHa//4KYaYlR3LZiCmvmZQYCd9CR000TUaZMEIWuyDiKjnASHeHEYbfh8fhwOuwjDqTlpo/tsvNydVD3gsg48fl8fPfPBwJDwmoay7HZbHzstln8+IXDNLT2kJ8Vz4dvKprgSsVKaumKjJP61p4RY3APnGygICuO21dOZbGRyrpFOUETm8vVTy1dkXGS6AonJtIZNI9vVkoML249w3ObTgOw26znUGkan7xjzkSVKRZTS1dknIQ5HTx8ywxiIv1tm9w0F/euK2DzkMv4AOw6dlaX3bmGqKUrMo4WGWkU5SbwwubT1LV0s/tYPTFRYUEXqYwKd2pO3GuIQldknP3i5WOBa6gdOtXEYiOV6oZO+t1ebDa4e22+QvcaotAVGUc9fW72D7loJcDJ6ja+9OBiXtx6hoS4SIpydBXna4lCV2QchTntuIadlZYQE853/3yAxoFL9bxdUskXH1xMbpprosoUC2mfRmQcOex2HrhxeqD7IDrCybzpKYHABeh3e9l0oHqiShSLqaUrMs6Wzc5gSmYsb+2uJDY6nNSEqBHLRIbrq3it0DstMs56+tz83z/uD1w6PTrCQUFWXGAOhsTYCK5fkD2RJYqFFLoi42zv8YZA4AJ09XqYMy2J21dOo6Khk/qmTvaeqGd1cRZhTvX4Xe0UuiLjzDZKjjocdprbe3h6Q+nAPTUcLWvm0bvmWlqbWE8/qyLjbMH0VLJSYgK3413hrCrOZMPeqqDlSsx6Wjv7rC5PLKaWrsg4iwhz8KUHF/Pargq6evu5ZdkU4qLDRxw8czhshOkkiaue3mERC/zqrybPbjrFX3dW8K3f7aOju5/bV04N6sO9eWke0ZFqB13t9A6LjLMTlS1sO1wbuF1Z38GGkkpuXzmNn3zhPWzaU0FWSgzTMuMmsEqxikJXZJw1t/eOvK/D33ebHB/Fyrm6GvC1RN0LIuNszrRkXFFhgds2IDMpit+9cYJXtp6m362rAV9L1NIVGWfRkU4+/+GF/HVnOV09btITo/jdm/6hYq/vrmB+YQr/dG/xBFcpVlHoilggMzmGh2+ZCcDjT+0JemxfaQONrT0kx0dORGliMXUviFgsPMwRdNtus+lMtGuI3mkRi922fEpQyF6/MJu4mPAJrEispO4FEYsZeYk8/ollHDrdxIz8FNJiFbjXErV0RSZAUlwka+ZlMTs/eaJLEYspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkM3n8010DSIi1wy1dEVELKTQFRGxkEJXRMRCCl0REQspdEVELKTQFRGxkEJXRMRCTqtf0DCMeOApIA4IBz5jmuY2wzCWAd8F3MBrpml+zTAMO/AEMA/oBT5mmmbp5S57GbXfBdxnmubfDLn9LaBiYJGvAJuupJrPU/cVv60H6rQBlcCJgbu2mab5ecMwbge+PPA6PzdN8yeGYUTh/1ylAe3AQ6Zp1l/KspdT6yX+XaNuP6tef5R6SoC2gZungR8xDp+PENW6FPhv0zTXGYZRCDwJ+IBDwKOmaXoNw/gKcOtATZ82TXNnKJYNRf0wMS3dzwBvmqa5FngY+P7A/T8E/gZYBSw1DGMBcCcQaZrmcuBfgf8J0bKXzDCM7wKPE7zNFgGfM01z3cC/jVdSzReo+4re1kMUACVDtu/nDcMIA74DrAfWAp8wDCMd+BRw0DTN1cCvgC9eyrKXWeelOt/2s5xhGJGAbcg2foTx+3xcbq2fA34KRA7c9W3giwPvow24wzCMhfjf66XABzmXL5e1bCjqHzQRofsd/L+k4G9p9xiGEQdEmKZ50jRNH/BX4D3438hXAUzT3A4sDtGy78ZW/F/WoRYBHzUMY5NhGP9jGIbzCqt5RN2TZFsPWgRkG4axwTCMlw3DMICZQKlpms2mafYBm4E1Q2sCXhl47UtZ1kojtp/Frz/UPCDaMIzXDMN4yzCMNYzD5yNEtZ4E7h5yexGwceD/g+/jKvwtbp9pmuWA0zCM1BAsGzLj2r1gGMbfAo8Nu/sR0zR3GYaRgX8X79P4uxrahizTDuQP3N865H5PiJZ9NzX/wTCMdcPufx14Dv8u2Q+BT05EzZdY9xWzrcdQ/6PA46Zp/skwjFX4Py+PDXvtdiB+WE2j3XexZa00YvsZhuE0TdNtcR0AXfi7yH4KTMcfMi1DHg/J5yMUf59pmk8bhjF1yF22gbAffO3B97ZxWE3xIVg2ZMY1dE3T/Bnws+H3G4YxF/g98M+maW4c+MWMHbJILP43PnrY/Xb8b/LlLnvJNZ/Hz03TbBn4m54H7sH/YbO05kusOxS1hKzuQaPVbxhGNP6+NkzT3GwYRhb+L8ForzP09Ue772LLWml4XfYJClyA4/j3BnzAccMwWoGkIY+H5PMxTn/f0H7Wi73nl7tsyFjevWAYxizgT8DfmKb5CoBpmm1An2EYBQMHT96L/4DUFuB9A+stw98PF4plQ/F32IADhmHkDNx1I7DnSq4ZJt22/gr+PSEMw5iH/4DlEWC6YRhJhmGE4+8u2Da0JuCWgdc+egnLWmnE9rP49Yf6KAN9rgM/atFAZ6g/H+NU+94he3GD7+MW4L2GYdgNw8jDH/gNIVg2ZCwfvYD/oE4k8F1/Fx2tpmnegX/X/DeAA38/yw7DMHYBNxmGsRV/h/YjA89xWcuG4o8wTdNnGMbHgGcMw+jGHwY/wb/bdUXWPMRk2db/BTxlGMbg0eWHTdPsNwzjM/j7D+349zaqDMP4AfBLwzA2A334f9THvOxl1nmpnmX07TcRfgY8ObAtfPhD2Mv4fD5C7bPATwZ+UI8CfzZN02MYxib8P652/F1Ul71sKIvW1I4iIhbSyREiIhZS6IqIWEihKyJiIYWuiIiFFLoiIhZS6IqIWEihKyJiof8PSwsHN4i5cF0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d7b1a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.swarmplot(population_a - population_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.JointGrid at 0x10d9923c8>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAGoCAYAAAA991BSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X143GWB7//3ZNIkTZO0STtdoQ/0Rx9udhH6QKUgYkstsAhddNEt9sDxyFkELAo/UXShLOjCXsiR+lsVOVpkQZRzqrAo6yVSuzxsKVQEWh7Okbu0lVYq2CTNY9NOksn8/shMOpPMTOZ5vjP353VdXFfmO/dMvt+bdD5z39/7wRcOhxEREal0VaU+ARERkWJQ4ImIiBMUeCIi4gQFnoiIOEGBJyIiTqgu9QnkorW1J+shps3N9XR09OXzdMqa6uMY1UU81Uc8L9ZHINDoK/U5lANnW3jV1f5Sn4KnqD6OUV3EU33EU32UL2cDT0RE3FLWXZq5+PULb9PTe7Qg771i0YyCvK+IiGRPLTwREXGCAk9ERJygwBMRESco8ERExAkKPBERcYICT0REnKDAExERJyjwRETECQo8ERFxggJPREScoMATEREnKPBERMQJCjwREXGCAk9ERJygwBMRESco8ERExAkKPBERcYICT0REnKDAExERJyjwRETECQo8ERFxggJPREScoMATEREnKPBERMQJCjwREXGCAk9ERJygwBMRESco8ERExAkKPBERcUJ1qU+gEj2z80BB3nfFohkFeV8REReohSciIk5Q4ImIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBM3DKyOFmt8H8MlzTyrYe4uIeIFaeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBI3SFAB+/cLb9PQezfv7aocHEfEKBZ6UrUJN09AUDZHKpC5NERFxggJPREScoMATEREnKPBERMQJGrQiBVXI9T8LpVAjVkGjVkVKSS08ERFxglp4IkVUji3eQk7TKFR9FLIlrTmr5UstPBERcYIvHA6X+hxEREQKTi08ERFxggJPREScoMATEREnKPBERMQJCjwREXGCAk9ERJygwBMRESco8ERExAkKPBERcYICT0REnFDWi0e3tvZkvS5ac3M9HR19+Tydsqb6OEZ1EU/1Ec+L9REINPrSLZvL52Y5SFUXzrbwqqv9pT4FT1F9HKO6iKf6iKf6KF/OBp6IiLhFgSciIk5Q4ImIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBAWeiIg4QYEnkifBgRAHO/oIDoRKfSoikkBZr6Up4gWhoSE2PbWbHbtaOdQdpKWplsULAqxZOQ9/lb5TiniFAk8kR5ue2s2Wl94ZedzeHRx5vHbVglKdloiMoq+fIjkIDoTYsas14XM7drWpe1PEQxR4Ijno6g1yqDuY8LmOnqN09SZ+TqRUntl5oNSnUDIKPJEcTG6opaWpNuFzzY11TG5I/JyIFJ8CTyQHtRP8LF4QSPjc4gXTqJ2gvdNEvEKDVkRytGblPGD4nl1Hz1GaG+tYvGDayHER8QYFnkiO/FVVrF21gEuWz6WrN8jkhlq17EQ8SIEnkie1E/xMb64v9WmISBK6hyciIk5Q4ImIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBAWeiIg4QYEnIiJOUOCJiIgTFHgiIuIEBZ6IiDhBgSciIk5Q4ImIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBAWeiIg4oXq8AsYYP7ARMEAYuBo4CjwQefwGsM5aO2SMuRW4EBgErrfWvmiMmZdr2fxdroiIuCqdFt5qAGvtWcB64A5gA7DeWns24AMuNsYsAZYDy4BLgXsir8+pbM5XKCIiI1YsmlHqUyiZcVt41tqfG2N+GXl4AtAJrAKejRx7AjgPsMBma20Y2G+MqTbGBIDTciz7WLJza26up7ran/bFjhYINGb92kqk+jhGdRFP9RGvnOsj18/NcjZu4AFYaweNMQ8CHwc+AZwbCSuAHmAy0AS0x7wsetyXY9mkOjr60jn9hAKBRlpbe7J+faVRfRyjuoin+ojnxfrIJIBz+dwsB6nqIu1BK9baTwMLGL6fNzHmqUaGW33dkZ9HHx/KsayIiEjOxg08Y8zlxph/iDzsYziUXjLGrIgcuwDYCmwDzjfGVBljZgNV1to2YEeOZUVERHKWTpfmvwH/aoz5T2ACcD3we2CjMaYm8vMj1tqQMWYr8ALDQbou8vobcimbj4sUERHxhcPh8Ut5VGtrT9Yn78V++FJSfRyjuoin+ojnxfoIBBp96ZbN5XOzHKSqC008FxERJyjwRETECQo8ERFxggJPREScoMATEREnKPBERMQJCjwREXGCAk9ERJygwBMRESco8ERExAkKPBERcYICT0REnKDAExERJyjwRETECQo8ERFxggJPREScoMATEREnKPBERMQJCjwREXGCAk9ERJygwBMRccgzOw+U+hRKRoEnIiJOUOCJiIgTFHjijOBAiIMdfQQHQqU+FREpgepSn4BIoYWGhtj01G527GrlUHeQlqZaFi8IsGblPPxV+s4n4goFnlS8TU/tZstL74w8bu8Ojjxeu2pBqU5LRIpMX2+logUHQuzY1ZrwuR272tS9KeIQBZ5UtK7eIIe6gwmf6+g5Sldv4udEpPIo8KSiTW6opaWpNuFzzY11TG5I/JyIVB4FnlS02gl+Fi8IJHxu8YJp1E7wF/mMRKRUNGhFKt6alfOA4Xt2HT1HaW6sY/GCaSPHRcQNCjypeP6qKtauWsAly+fS1RtkckOtWnYiDlLgiTNqJ/iZ3lxf6tMQkRLRPTwREXGCAk9ERJygwBMRESekvIdnjJkA3A/MAWqB24E/Ar8E3ooUu9dau8kYcytwITAIXG+tfdEYMw94AAgDbwDrrLVDmZTN47WKiIjDxmvhXQa0W2vPBv4a+C5wGrDBWrsi8t8mY8wSYDmwDLgUuCfy+g3A+sjrfcDFmZTN10WKiIiMN0rzZ8AjkZ99DLfITgOMMeZihlt51wMfAjZba8PAfmNMtTEmECn7bOT1TwDnATaDso/l4RpFRERSB561thfAGNPIcPCtZ7hr8z5r7cvGmJuBW4FOoD3mpT3AZMAXCbbYY00ZlE2pubme6urs51MFAo1Zv7YSqT6OUV3EU33EK+f6mFRfU9bnn4tx5+EZY2Yx3NL6nrX2YWPMFGttZ+Tpx4DvAL8AYmuwkeEQHEpwrDuDsil1dPSNVySpQKCR1taerF9faVQfx6gu4qk+4nmxPjIJsMN9/Z47/3xKVRcp7+EZY/4C2Ax8xVp7f+Twk8aY0yM/fwR4GdgGnG+MqTLGzAaqrLVtwA5jzIpI2QuArRmWFRERyYvxWng3Ac3ALcaYWyLHvgh8yxgzALwHfNZa222M2Qq8wHCIrouUvQHYaIypAX4PPGKtDaVbNi9XKCIiAvjC4fD4pTyqtbUn65P3YrdEKak+jlFdxFN9xPNifQQCjb50y/7sN2+GVyyaUcjTKalUdaGJ5yIi4gQFnoiIOEGBJyIiTlDgiYg4pJLv341HgSciIk5Q4ImIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBAWeiIg4QYEnIiJOUOCJiIgTFHgiIuIEBZ6IiDhBgSciIk5Q4ImIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBAWeiIg4QYEnIiJOUOCJiIgTFHgiIuIEBZ6IiDhBgSciIk5Q4ImIiBOqS30CIiJSPM/sPBD3eMWiGSU6k+JTC09ERJygwBMRESco8ERExAkKPBERcYICT0REnKDAExERJ6SclmCMmQDcD8wBaoHbgf8LPACEgTeAddbaIWPMrcCFwCBwvbX2RWPMvFzL5vVqI4IDId5tO0xoIETtBH8hfoWIiHjMeC28y4B2a+3ZwF8D3wU2AOsjx3zAxcaYJcByYBlwKXBP5PU5lc3PJR4TGhri4S27WL9xO1fduYX1G7fz8JZdhIYKkqsiIuIh4008/xnwSORnH8MtstOAZyPHngDOAyyw2VobBvYbY6qNMYE8lH0sx+uLs+mp3Wx56Z2Rx+3dwZHHa1ctyOevEhERj0kZeNbaXgBjTCPDwbce+GYkrAB6gMlAE9Ae89LocV+OZVNqbq6nujq9Lsmj/YO8tqc94XOv7WnnqksmUlfj9sIzgUBjqU/BM1QX8VQf8cq5PibV11BVdaxzr5yvJVPjfsIbY2Yx3NL6nrX2YWPMXTFPNwKdQHfk59HHh3Ism1JHR994RUYc7OijteNIwufaOo+w5+12pjfXp/1+lSYQaKS1tafUp+EJqot4qo94XqyPTELrcF9/3GOvXUuuUtVFynt4xpi/ADYDX7HW3h85vMMYsyLy8wXAVmAbcL4xpsoYMxuosta25aFs3kxuqKWlqTbhc82NdUxuSPyciIhUhvFaeDcBzcAtxphbIseuA75tjKkBfg88Yq0NGWO2Ai8wHKLrImVvADZmWzYvVxhRO8HP4gWBuHt4UYsXTNNoTRGRCucLh8Pjl/Ko1taejE4+NDTEpqd2s2NXGx09R2lurGPxgmmsWTkPf5XbUxK92E1TKqqLeKqPeF6sj0Cg0Zdu2Z/95s24z81K2y0hVV04NUrDX1XF2lULuGT5XPw1Ewj1D6hlJyLiCCebNbUT/Bw3bZLCTkTEIU4GnoiIuEeBJyIiTlDgScULDoQ42NFHcCBU6lMRkRJyatCKuOXYqNxWDnUHaWmqZfGCgEblijhKgScVS2unikgsfc2VihQcCLFjV2vC53bsalP3poiDFHhSkbp6gxzqDiZ8rqPnKF29iZ/LJ907FPEWdWlKRYqundqeIPQKvXaq7h2KeJP+9UlFiq6dmkih106N3jts7w4S5ti9w01P7S7Y7xSR8SnwpGKtWTmPVUtnMrWpjiofTG2qY9XSmaxZOa9gv1P3DkW8S12aUrFi107t6g0yuaG24MvJpXPv0OV9F8VbKm3h6PGohScVr3aCn+nN9UVZO1X7Lop4lwJPJI9Kee9QRFJTl6ZInkXvESbad1FESkeBJ5Jnpbh3KCLjU+CJFEj03qGIeIPu4YmIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLiBAWeiIg4wcnACw6EeLftsFauFxFxiFMTz+M25uwJ0tJY3htzBgdCWslDRCRNTgVedGPOqOjGnABrVy0o1WllTDtqi4hkzplPx0ramFM7aouIZM6ZwEtnY85yUEnBLSJSTM4EXqVszFkpwS0iUmzOBF6lbMxZKcEtIlJszgQeDG/MuWrpTKY21VHlg6lNdaxaOrOsNuaslOAWESk2p0Zpxm7M6a+ZQKh/oCwDQjtqi4hkzqnAi6qd4CcwbRKtrT2lPpWsaEdtEZHMORl4lUI7aouIpM+pe3giIuKutFp4xphlwDestSuMMYuBXwJvRZ6+11q7yRhzK3AhMAhcb6190RgzD3gACANvAOustUOZlM3XhYqIiNvGbeEZY24E7gPqIodOAzZYa1dE/ttkjFkCLAeWAZcC90TKbgDWW2vPBnzAxZmUzccFioiIQHotvD3A3wIPRR6fBhhjzMUMt/KuBz4EbLbWhoH9xphqY0wgUvbZyOueAM4DbAZlH8v1AkVEJLFndh5I+tyKRTOKeCbFMW7gWWsfNcbMiTn0InCftfZlY8zNwK1AJ9AeU6YHmAz4IsEWe6wpg7IpNTfXU12d/ejEQKAx69dWItXHMaqLeKqPeOVcH5Pqa6hKY5H5cr7GZLIZpfmYtbYz+jPwHeAXQGztNDIcgkMJjnVnUDaljo6+TM99RCDQWLbTEgpB9XGM6iKe6iOeF+sjk3A63NefVjmvXWO6UtVFNqM0nzTGnB75+SPAy8A24HxjTJUxZjZQZa1tA3YYY1ZEyl4AbM2wrIiISF5k08K7BviOMWYAeA/4rLW22xizFXiB4RBdFyl7A7DRGFMD/B54xFobSrdsthdVSXLd5FWbxIqIDPOFw+HxS3lUa2tP1ifvxW6JWLlu8prp671eH8Wkuoin+ojnxfoIBBp96Zb92W/eTOtzs1wHraSqC6204lG57s5eKbu7i4jki1Za8aBcN3nVJrEiImMp8Dwo101etUmsiMhYCjwPynWTV20SKyIylgLPg3Ld5FWbxIqIjKVBKx6V6yav2iQ2M7HTN0SkMinwPCrXTV61SWx6Ek3fOGvhDFafOTut6R8iUj4UeB6X6yav2iQ2tUTTNx7fupe+I/2aviFSYfQVVpyl6RsiblHgibM0fUPELQo8cZamb4i4xenACw6EONjRp64rR2n6hohbnBy0EhoaYuPPX2fbqweyWpi5VLTzQf4lmr5x1sLjWX3m7BKfmYjkm5OBV24LK+e6c4Ikl2j6xszjp5R8NXx9uRHJP+cCb7yReZcsn+u5D5hyC+hy5JXpG/pyI1I4zv0LKreReRo675bol5v27iBhjn252fTU7lKfmkjZcy7wJjfUMqWxfEbmlVtAS/b05UaksJwKvP7BQf75oZfp6EkcEl4cmaeh8+7QlxuRwnIq8O740Sv88WDvmOP+Kh+rls705MLKGjrvDn25ESksZwKvp6+fA61jww5gaCjM6g/O8eyggDUr57Fq6UymNtVR5YOpTXWeDWjJnr7ciBSWM6M03znYy1A48XPhyPN/OaelqOeULu184A5t6yRSOM4E3szpDVT5SBh6Vb7h573OK0PnpXD05Ua84pmdB1ixaEapTyOvvNmHVwCN9TXMCCQOtRmBBhrra4p8RiLJRb/cKOxE8seZwAO4+b8uYVakpQfDLbtZ0xu4+b8uKe2JiYhIwTnTpQlQU13N1644nZ6+fnr6h2isqVLLTkTEEU4FXlRjfQ0nntBY8vUSRUSkeJzq0hQREXcp8ERExAkKvDKmDWxFRNLn5D28cqctZEREMqfAK0PaH09EJHNqDpQZbSFTWMGBEO+2HVY9ilQgtfDKTDpbyGj5sczFdRP3BGlpVDexSKXRv+Qyoy1kCiNup/GwdhqX4tIAtOJQC6/MRLeQib2HF6UtZLIzXjfxJcvnql6lIDQArbgUeGVIW8jkl7qJpVQ0AK24FHhlSFvI5Fe0m7g9Qeipm7h4ggMhp/6e1bNQfGkFnjFmGfANa+0KY8w84AGG9019A1hnrR0yxtwKXAgMAtdba1/MR9n8XWrl0f54+aFu4tJytVtPPQvFN+5fkzHmRuA+oC5yaAOw3lp7NuADLjbGLAGWA8uAS4F78lE298sTSc+alfNYtXQmU5vqqPLB1KY6Vi2dmbKbWAMN8iNuwBDuDBjSALTiS6eFtwf4W+ChyOPTgGcjPz8BnAdYYLO1NgzsN8ZUG2MCeSj7WKoTa26up7o6+2/fgUBj1q+tRLH1cbR/kI7uIM1NtdTVuNHzfd2nTkvrukOhIe7/9//D9jfepbXzCIEpEznj/cdxxeqT8fsrs0VSqH8rR/sHeW1Pe8LnXtvTzlWXTPTk31++6uOshTN4fOveBMePZ+bxU/LyO0abVF9DVZot50r7jBz3L8la+6gxZk7MIV8krAB6gMlAExD7Vxs9nmvZlDo6+sYrklQgoO2BYkXrw9XupVjHReoi2V/Hw1t2xXV/Huw4wuNb99J3pL8iBxoU8t/KwY4+WjuOJHyurfMIe95u91y3Xj7rY/WZs+k70j9mANrqM2dn9DsyCabDff1ply3Hz8hUdZHNV6fY+2qNQCfQHfl59PFcy5aVSrjprlFjqWmgQX65PmBIA9CKK5uv7DuMMSsiP18AbAW2AecbY6qMMbOBKmttWx7KloXQ0BAPb9nF+o3b+Yfvb2f9xu08vGUXoaHyGnOjZcvGl85AA0lfdMBQIi4NGIoOQHPlekslmxbeDcBGY0wN8HvgEWttyBizFXiB4RBdl4+y2V5UsXm9VZROy/No/yB7D3Ql/KYNGjUW5XqLpBA0r1SKxRcOh8cv5VGtrT1Zn3y++uGDAyHWb9ye8ANwalMdt1+5rGTf2tK5Hxct89qedlo7juDzwVCCWi31tRTTeH8bo+/hRa1aOtMTX3DyrVj3u8vlloAX7/8HAo2+dMv+7Ddvpv25uWLRjOxOqIRS1YX3hj+VGS/PpUmn5Tm6TLLvPy51L41HLZLC0LxSKTQFXo682sWVzuCK4Z8Tl6nyDYdfS1PxPszL5Ru+BhqIlCcFXo68ukpHuoMrkpUJA1+6dBEnzphc8Gso16kQapGIlBcFXh54sYsr3ZZnsjItjXVFCTvw/qAfEVc9s/NA0ufK8f6eAi9Gtl1qXuziSrflmWvrNNduSM1rE5FiUeCRvy41r3VxpdPyjP782p522jqPpN06zVedeXnQj4hUFgUeldullk7LM1rmqksmsuft9rRbavmqM68O+hGRyuPdEQFFUu6ri6SzYn86qzjU1VSnvdJDPutMK22ISLE438Ir1y61Uo5szHedeXHQj4hUHucDr1y71ErZDZvvOvPioB8RqTzOd2mWY5daqbthC1VnWkBXRArJ+RYelFeXWnAgxN4DXSXvhi2nOhMRAQUeUPwutWzmrsXes2vvDo4s/TVasbph1Q0pIuXGycALDoR4t+0woYFQ3Id0oefR5TLQZPQ9u0Q7GkDxu2G9NvdQRCQZpwIvLnB6grQ0Fm9kY3AgxENPWp5/472RY+kONEl1z67KN7zuZYu6FEVEUnIq8EoxsnF0V2QisbsXJOoeTDUNIBwu3iLPIiLlzJnAK9WajaNDNpFD3Uf58ZOWN/d3jHR1njp3KquWzqKlqS7lNICWpuIt8iwiUs6cCbxSTDBPFbKxamv8bBvV1fn0jj/x9I4/MTVyn2/h/Gk89fLYlcu9OnVCRMRrnJmHF20lJVKokY2pQjZd0W5XH7Bq6UymNtVR5YOpTXWsWjrTiXt26SyfJiIyHmdaeKXYqDVVVyRAS2MtJ53QzAsxrbtkdr7Vzu1XLnNqGkCqUa2DobAz9SAi+eFM4EHxJ0unCtmz3v8+LjvfAGD3dyQNxajYbldXpgEkG2Rk93fSd3SgrHZHF5HScyrwYidL+2smEOofKHjrIFXIRj+gk4VirGi3a64brmaimL8r0e9Odv/zjwd7R36ulK2cRKTwnAq8qNoJfgLTJtHa2lPw35XOiiSxodjefTTh+yyaP5VHn91TlN0RSrkTQ1Sm9z+1O7qIjMfJwCuFVCuSxIbioe6jbHn5HV7b3R7XIhwKh/mPIs0h9MKGuOPd/xzNy1s5iYg3KPA8pHaCn+OmTuLy8wzBc451JwKs37g94Wvy3bIp1XzF0VLd/0zEy1s5iYg36C6/R8VulZPOHMJ8KebvGs+alfPGTMWYNb0hYVnNRxSR8aiFVwaKuUmtlzbETXT/s9rvi9xf1LZEhVLKwUoihaTAy1ApPgyKOYewFPMV0zmn2Htz2paoMJINVrr27xaX+tRE8kKBl6ZSj1ws5hzCctjcNZNtidRiSU+ywUr1E2v42FlzSndiInmiwEtTqUcuFnPD1UrZ3LXUX1LKSarBStvfeJcLTp9Vln8DIrH0rz4N441cLOYaj7GDWSrpdxVC9EtKe3eQMMe+pGx6anepT81zUg1Waus8UtTBSlIentl5YOS/cqHAS4OXRi5Kerz0JaUcpFpcfdqUiZryIRVBgZeGUuy04AXlvEuBvqRkJjpYKZEz3n9c2bbyRWLpHl4Sowc65DpyMZ8DJxK9Vz7fvxLufXlpekW5SDZY6YrVJ3Po0OESn51I7hR4oyT7sP/EihOBzEcu5jM8Er3XovnTCAOvvtWWt3Aq9QCdfPDi9AqvSzZYye8vjy85IuNR4I0y3od9piMX8xkeid7rP0btgp5rOHllabF8KIfpFV6UyZQPkXKSdeAZY14BuiMP/wB8H/gXYBDYbK39mjGmCvgesBAIAn9vrd1tjDkj3bLZnl820v2wz2T+V77CI9V75eP9o9K591UuH4aVMr1CRPIjq74KY0wd4LPWroj89xngfwJrgQ8By4wxi4GPAXXW2jOBrwJ3R94ik7JFk++BDvl8v0y3y8l2YIaXB+hkO4im3KdXiEh+ZNvCWwjUG2M2R97jNqDWWrsHwBjzJLAKOA74NYC1drsxZqkxpindsllfVZbyPdAhn++X6XY52YZTqntfZvaUjN8vHyphEI2IlF62gdcHfBO4D5gPPAF0xjzfA5wINAFdMcdDkWPd6ZQ1xlRbaweTnURzcz3V1dl/aw8EGsccO2vhDB7fujfB8eOZeXzmH/j5fL9k75W4bObvH62Pa/9uMfUTa9j+xru0dhyhrtYP+Hjh/7zH7gNdnPH+47hi9clFG8yw8eevJ13y6sqPnVKQ35nob8Nlqo945Vwfk+prqMrzF8VyqY9sA28XsNtaGwZ2GWO6gJaY5xsZDsD6yM9RVQyHXWM6ZVOFHUBHR1+Wpz/8PyjRjuerz5xN35H+MQMdVp85O6sd0vP5fonea9H8qZFRmu05vf/o+vjYWXO44PRZ/PhJy7Y33hs5frDjCI9v3Uvfkf6ijNgMDoTY9mrilRy2vfqngix5lexvw1Wqj3herI9MAudwX3/ef7+X6iNVXWQbeFcApwCfM8Ycz3BYHTbGzAX2AucDXwNmAquBn0YGqrxure02xvSnUzbLc8tJvgc6JHu/4ECI9q6+jN4/1bl9ckVhFkh+c39HwuPFGrFZSYNoRKS0sg28HwIPGGOeA8IMB+AQ8BPAz/DIy98aY34HnGuMeR7wAZ+JvP7qDMqWRL6HZkffLzQ0xMNbduV0PyrRuRViKLkXwkYTyEUkX7IKPGttP8OjLEc7Y1S5IYbDbfTrt6dbttKU06RuL4SNJpCLSL5oiFsRZbugcanWtEy1vmIxw2bNynmsWjqTqU11VPlgalMdq5bO1ARyEcmIVlopoky7CL0wHN8Lq5VoArmI5IMCr4gy7SIsZfdn7GLUXgkbLXklIrlQ4BVROvejokEzsba6JGtapmpVKmxEpJwp8Ipszcp5hIbC7NzVRufhIC2RLsJPrDgxbvTmlIZaOpIsDVbIEZLlNKhGRCQTGrRSYLEDTqKtp9d2t9HmIhEkAAASaklEQVTRG2TypBpOndvCmpXzeOSZvWx56R3au4OEIWnYATRNqmFibf6/q2iXcBGpZGrh5UGizVcTdQ3W103gjwd7R17X2dvP0zv+BMBre9rT/n2dvf18/YHf5X0Aixfm3YmIFIoCLwep7ncl6hpMtvDzjrfa6OpNvtzPlIYaOkc9X4iuRi/MuxMRKRR1aeYgGmrRbshoCD285a2M9q7r6u1nSpIwmdpUx82Xn0Zzkufz2dXolXl3IiKFoMDLQOz9uFT3u3buakt7Gx+AlqY6Fi2YlvC5xQumERoK0znOAJZ80SRvEalU6tJMQ6KuSzO7Oen9rs7DwYTdkMmcOreFtavm46/yJZzgPRgKF62rUZO8RaRSKfDSkOh+3PNvvIe/ykdoKDymfO0EP4vnTxsZkDKaj+EVt6t8MBQeHrDi9+9mzcp5CYPGX0XR15PUJG8RqTQKvHGk6rpMFHZRl6yYh99fxY5dbbR3H417Lvqq6MtHD0BJFDReWOKr3CQaPSsi7lLgjSPVUP1kgv0hevuGN0hd/cE53Hb/71LOq4t66c2DrP7gHBrra8Y8p67G9HlhDVIRlzyzM/EmzcWyYtGMtMrpX/84okP1E6nyJX5NS9Ox+2pHgoNJB5yM1tnbz233/46Ht+wiNDSUsEy0qzFR2JVqVwWvSTZ6dtNTu0t9aiJSQmrhRSTr/kq1/uWMQEPcRPKo2Ptqkxtqqa3xc7Q/vRDq6M18fp1aNMeMt1pMMXZpFxFvcj7w0gmLVOtfPvLM3nHvq4XDiVtrqWTy4az1L4/RajEikozzgTdeWIxe/3JKw7H1L2Pvq7V29IHPR2DKxLhWVVdvkOBA8sEtyaT74awWTTytFiMiyTgdeKnDopVLls/l0Wf3xAVidP1Lv79qJBAffXZP0hbi5IZaWhprONSTeE6ezwfhBHmY7oezWjTx0tmCSUTc5NYNnlG6epOvb9neHaS188i4uweMN0CidoKfJWZ60nOYEZiU8Hi6H86pBtW42qLRajEikojTLbyJtdUjk79Hq/JB/8BgytZTa0dfWt2Ja1bO4819HbzTenhMuXkzJ3PS7Oas59epRTOWpnCISCJOB96R4GDCsIPhEPT7/dTWVHG0f+ygkykNtRzqSd5CjO1OHAyFORIcTFju9d2HuP3KZTl9OGtSemJaLUZEYjkdeKnur7U01vKfr/4pYdgBHD4ywL/87LWkLcTY7sR077Nl++EcbdGs/uAc3jnYy8zpDQknr4uIuMzpwIveX0vUHbhw3lRe292W9LXBweEgTDTgBMbOxSvkyEHNwxMRGZ/Tn4bBgRDnLJ7BOYuPHzPAYdXSWRkvKQZQV+Nn5Wkz4roTC73PnFYWEREZn5MtvNDQEBt//jrbXj0w0iI6dd40Vp02k5amOgBaO/qStspSOdofosrnG9OyKtR9Ns3DExFJj5OBl2iy+dOvHGAwFOJw3yB/eLeLzt4BamuyawAnCppCjRzUPDwRkfQ416WZqkW09dX3eOWtNjp6BwhD0gEr42nvPsqhUVsCRcWOHMzHQs+ahycikh7nWnjZbPeTjS0v/ZHLzz9pzPF8DzDRPDwRkfQ4F3ipRkzm02t7DhEcCI0JnEIs9Kx5eCIi43Mu8Gon+Dl13jSefqWwGxYmun9WqAEmWllERGR8Tt3DCw0N8fCWXbz61nDoJNvANR8S3T9LZ4BJLlJtDisi4jqnAi/anRhdWSXZsmL5UF9XTbU/PlE1wEREpHScCbzgQIhX7MG8vNfkSRNGJqmvPG0GM6eP3fHgjwd7x0z8LvQE9FwFB0J5GTkqIuJFztzD6+oNJt2TLlNfuORUJk2cMNIiW79xe8Jyie7LeXGASaKJ+FqaTEQqjTOBN7E2P5daV+Pn+EDDSIgd7OjLaOK3FweYFGLkqIiI1zjz9b3rcH5ad2ed8r64gMr2vlyiASal6FIcb+SoujdFpFJ4qoVnjKkCvgcsBILA31tr87MCcrJtDVKYNb2BvqMDHOoJ0tJ4rJsvVj4mfpdytwMtTSYirvBU4AEfA+qstWcaY84A7gYuzscbB5rr8VdBKM3VwmonVPGV/7IEf5Vv3K7HXO/LlbJLsdBbF4kkcv/9P+CFF57D76/mC1/4In/1V+9PWO4b37iDpqYmrrnm8/zqV//Or3717wD09/eze/cufvGLJ2lsbCzIOT700AP89rfPA9Db28uhQ+08/viTY8odPXqUq6++gquvvpYzzvhgWu/9zjt/5I47bsPn83HiiXP54he/wosvbufHP35gpMxrr+3kRz/axJw5/09erke8F3gfAn4NYK3dboxZms83r66qIjSUXuINDA7R29ef1sasudyXK/VuB1qaTIrN2jfZufMVfvCDB/nzn//M+vU3ct99PxpT7uc/f5S9e3ezaNESAD760dV89KOrAbj77m9w4YV/U7CwA7j88v/G5Zf/NwBuvPF6Pve5LyQst2HDN/BlOKf3O9/ZwJVXXsOSJUv5H//jn9m69VmWLz9nJDAffvhHnHLKQoVdnnkt8JqArpjHIWNMtbV2MFHh5uZ6qqvT+0B+t+0w/ek274BpUyYyd85U6moyq6KZGZUePq9DPcm7FP01EwhMGzvtIZ+u/bvF1E+sYfsb79LWeYRpUyZyxvuP44rVJ+P3O3ObN04gULgP0mz827/9G1u2bOHw4cN0dHSwbt06zj//fF588UW+9a1v4ff7mTVrFl//+tcJBoPcfPPN9PT0cPDgQdauXcvatWu5/PLLaWlpoauri3/8x3/kpptuorq6mqGhIe6++26OO+447rzzTl5++WUALrroIj796U/z1a9+lZqaGg4cOMDBgwe58847OfnkkznnnHM48cQTmTt3LjfddNPIuV511VX09fWNPJ47dy633XbbyOMnnvg955yznOnTm5g+vQmfL4zfP0BLS8tImVdeeYW9ey2XXbaWvXv3xv3/eP311zlwYB933nk7AD/4wQ846aST+PCHPzxufUXt27eP9evXx9XxRRddxJo1a8bU/ebNm5k2rYULLzx35Fj0fH74wx9y5pmnU1tbzeTJEwkEGnn33Xe55ZZbCAaD1NbW8k//9E8cd9xxce/51luW885bgc/n47zzPsK2bdv4xCf+BoD33nuPLVt+zaOPPkpNTU2yP4msffyc+Wl/blYarwVeNxD7SVOVLOwAOjr6kj01RmggREtj+mtonjp3Kj1dR+hJ+zdkJ9V5NTfWEeofoLW10GcBV37sFC44fVZcC/XQocMF/71eFAg0FqXOM9HTc5Surh42bPgunZ0dXHnlpzn11NP5h3+4mXvvvY/m5hY2bryXH/3of2HMX3L22StZvnwlbW2tXHvtZzn33NX09w9y9tkfYfnyc3j00Z8yf/5JfO5z1/HqqzvYv/89fvvbHezZ8zb33PNDQqEQ11zz3zHmVAAmT57K5z//ZR5//DEefPDHfPnLN/Huu++yceOPmDx5Slx93X77N8ecf+zzf/5ze9xramrq2LfvPUKhCQC0tbXxrW/9C//8z9/kqad+Q19ff9zrv/3te7jssitGjn38458a8zuS1Vd19fBHXn19Cxs2fC/leUbdc8+93HbbHSPPRf8+XnrpRd588y1uvPFmnn/+t3R1HaG1tYevf/0O/uZvPsGZZ57FSy+9yB133Mmtt94e956h0BBtbb0ADA5W0dbWMfL+3/ve97nkkkvp6goyPJRhfJl8Qcvkc7McpaoLrwXeNmA18NPIPbzX8/XGqbruGiZWU1Ptp7M3WPR5cV7qUozduki8Z9GiJVRVVdHSMpXGxiba2lppb2/jllu+CkAwGOQDH1jGmWeexU9/+jDPPvs09fWTGBw89p1x9uwTALjooov5yU8e5IYbPs+kSQ1cddU69u37AwsXLsLn81FdXc3JJ5/C22/vBWD+fAPA9Ol/weuvvwrA5MlTmDx5ypjzvPHG6+NaeHPmnMiXvvTVkceTJjXQ13fsy1Rf32EaGo59SD399BY6Ozv50pe+wKFD7Rw9epQTTpjDRz+6mp6eHvbv38eSJePf7RhdX52dnUybNg0Yvod2553/FFf+3HP/mosv/tu4Y3/4w14aGhqYOXPWmPf/5S9/wZ///C7XXvtZ9u/fx65dbzJ16lT27t3NQw/9Kz/5yYMA+P3VPP30Fh599KcAXHvt/0tVzGC04etvAGBoaIht257js5/93LjXJ5nzWuA9BpxrjHke8AGfyeebjx5cMnlSLYsWTGPtqvkMhsIlmxfnxcno4j3WvgnAoUPtHD58mEBgOtOnT+fOOzfQ0NDAc889y8SJ9fzv//1j3v/+U/n4xz/BK6+8xAsvPDfyHtEP2ueee5aFCxdzxRWf5Te/+TU/+cmDLF++kl/96nHWrPkvDA4O8sYbr3HBBRexc+fv8CW4SVWVZATxXXf9fymv45RTFnLvvd/mU5+6nIMHDzI0FGbKlGPB+clPXsonP3kpAL/61b+zb9/bI/fuXn31FZYu/UBW9dXc3Dzy3MyZs/jud38w7nu89NKLSQei3HbbHSM/33HHbXzkI+cxf75h9uw5fOpTl3HKKQvZt+9tdux4mXPOWcU556waKT9/vuGVV15iyZKlbN/+/EiA7927hxNOOIHa2rq0rlEy46nAs9YOAVcX6v1jB5f4ayYQ6h8YCTd/FSVr3XhxMrp4z6FD7Vx33TX09vZyww1fwe/3c911X+LLX76OcDhMff0kbrnla/h8Pr71rbv4j//YTENDA36/n/7++HmoJ530V9x++608+OAPGRoa4vOf/yLGnMSOHS9z1VWfYWBggJUrV2HMSfzyl/m9jpNO+ktOPXURV131GcLhMF/84lcA2Lz51xw50jemlRVr//59HH/8jLhjDz30APPnLxgTTInqK1P79+/jAx9YFnfsrrvuYtmys5OOLF237jruvvtO+vv7CQaPct11XxpT5tprr+euu+7g+9+/hxNOmMOKFR+Jub5MRwJIunzhLOaneUVra0/WJ+/F+zSlpPo4xot1EW3pXHPN54v+u71YH7GiLdvTTjvW8itkfXmxPgKBxrTHiebyuVkOUtWFp1p4IiKZmjfP8L73va/UpyFlQIEnUgai97BkrERhp/qSRNycZCUiIs5R4ImIiBMUeCIi4gQFnoiIOEGBJyIiTlDgiYiIExR4IiLihLJeaUVERCRdauGJiIgTFHgiIuIEBZ6IiDhBgSciIk5Q4ImIiBMUeCIi4gQFnoiIOMG5/fCMMVXA94CFQBD4e2vt7tKeVX4ZY14BuiMP/wB8H/gXYBDYbK39WrJ6MMackW7Zol5UFowxy4BvWGtXGGPmAQ8AYeANYJ21dsgYcytwIcPXe7219sV8lC3mdaZrVH0sBn4JvBV5+l5r7SYX6sMYMwG4H5gD1AK3A/8Xx/8+XOBiC+9jQJ219kzgq8DdJT6fvDLG1AE+a+2KyH+fAf4nsBb4ELAs8mGXrB4yKetZxpgbgfuAusihDcB6a+3ZgA+42BizBFgOLAMuBe7JR9lCX1s2EtTHacCGmL+TTQ7Vx2VAe+Qc/xr4Lo7/fbjCxcD7EPBrAGvtdmBpaU8n7xYC9caYzcaYp4wxHwZqrbV7rLVh4ElgFQnqwRjTlG7Zol9V5vYAfxvz+DTg2cjPT3DsujZba8PW2v1AtTEmkIeyXpSoPi40xvynMeaHxphG3KmPnwG3RH72Mdwic/3vwwkuBl4T0BXzOGSMqaSu3T7gm8D5wNXAv0aORfUAk0lQD5Fj3emU9XqdWWsfBQZiDvkiIQ7Jryt6PNeynpOgPl4Evmyt/TCwF7gVR+rDWttrre2JhPwjwHoc//twhYuB1w00xjyustYOlupkCmAX8OPIN81dDP8jbIl5vhHoJEE9JDiWtGwZ1lnsfZNk1xU9nmvZcvCYtfbl6M/AYhyqD2PMLOBp4CFr7cPo78MJLgbeNuCjAJEBGq+X9nTy7goi99iMMccD9cBhY8xcY4yP4ZbfVhLUg7W2G+hPp2xxLykvdhhjVkR+voBj13W+MabKGDOb4SBvy0PZcvCkMeb0yM8fAV7GkfowxvwFsBn4irX2/shh/X04wNPdUgXyGHCuMeZ5hvvvP1Pi88m3HwIPGGOeY3hk2BUMf8v8CeBn+D7Db40xvyNxPVydQdlycgOw0RhTA/weeMRaGzLGbAVeYPjL37p8lC3aFeXmGuA7xpgB4D3gs9babkfq4yagGbjFGBO9l3cd8G39fVQ2bQ8kIiJOcLFLU0REHKTAExERJyjwRETECQo8ERFxggJPREScoMATEREnKPBERMQJ/z/9GC3Qgt5REQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d992198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.jointplot(population_a, population_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.JointGrid at 0x10ddb7ba8>"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAGoCAYAAAA991BSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4nFWh7/FfMpNMmiZt0zYVaAuVli5uQktbCoK0liIbFQFBCmxvgIAICogUkWJhC4KIgm4FpNgjHkEr4AX3EQUscEqBzQZaLke6YkEFCtIkTZtbM0kmc/6YTJqkk7lk3pl3Ztb38zx9nuSdNTNr3kznN2u961IWjUYFAECpK/e7AgAA5AOBBwBwAoEHAHACgQcAcAKBBwBwQtDvCmSjsbHN9yGmdXXVamnp9LsaaaO+uUV9c4v6JlZfX1uW8ycpAbTwshQMBvyuQkaob25R39yivsgGgQcAcEJRd2kCGNkTG7ekLFNbU6W29q6Ety2eM9XrKgG+ooUHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHBC0O8KABjZExu3+F0FoGTQwgMAOIHAAwA4gcADADiBwAMAOIHAAwA4gcADADiBwAMAOIHAAwA4gYnnQA4V88TxbOq+eM5UD2sCeIMWHgDACQQeAMAJBB4AwAkEHgDACQxaAVLIZPBGbU2V2tq7clgbAKNFCw8A4ARaeEhbtkPssxmqXszD+wEUBgIPgOf8/HIEjIQuTQCAE2jhASgptC4xEmcDz6trQsU2Ku9Tx+3v23OP5pwX2/mFN0Z6r+Tj/eDl9eJM6kvQ5h5dmgAAJ5RFo1G/6wAAQM7RwgMAOIHAAwA4gcADADiBwAMAOIHAAwA4gcADADiBwAMAOIHAAwA4gcADADiBwAMAOKGoF49ubGzzfV20urpqtbR0+l2NtFHf3KK+uUV9E6uvry1Lt2whfG7mUrJzQQsvS8FgwO8qZIT65hb1zS3qi2wQeAAAJxB4AAAnEHgAACcQeAAAJxB4AAAnEHgAACcQeAAAJxB4AAAnFPVKKwCK06uvvqIf/OAWBYMBLVhwhM455/yE5TZseEHf+tY39Zvf/J+BY11dXbrssi/p61//pvbZZ0ZO6/noo3/Sr3/9SwWDAe277yxdfvnXVV6+q53Q2rpDZ575Sb3//TMlSccc82GdfvqZaT326tV36ZlnnlIgENRXvvJVHXjgwWpo2KTlyy/TtGnTJUmnnHKajj32I96/MEc5GXjhnoh2tIc1viakUAUrIQD5dsstN+qGG27WXntN1RVXXKKGhk2aPXv/IWXee+9fWrPmXvX29g4ce+WVV3T11deosXFrzusYDndp1ao79POfr1FVVZVWrvyGnn56nY4+etFAGWs3aenS43XZZcszemxrN2njxhd111336L333tOKFct1990/l7WbtGzZv+vMMz/t9cuBHAu8SF+f1qzdrA0NjdrWGtbEcSHNnV2vZUtmKVBO7y4Kyx//+AetW/eEOjs7tX37dp199he0ePGx2rDhBd111+0KBALaa6+pWr78aoXDXbrppuvV3t6mpqZGffKTp+uUU07TxRefr7q6iWptbdXlly/XjTf+hwKBoPr6+rRy5fWqr6/Vf/7nrXr55Y2SpOOO+zedfvqZuuGGa1VRUaF//etdNTc36RvfuFbG7K9TT/249tlnhmbMeL++8pXLB+q6fPml6uzctWbkjBn76mtf+3rC19XR0a6enm5NnTpNknT44Ufq+eefGxJ44XBYt9xyo5Yvv1rnnvuZgePd3d369re/q29965sDx5qbm/TDH35P111345Dn+fSnP6VDDpmjv//9DY0bN07XXvttjRkzZuD2u+66feB1x916649VUVEhSaqoqNSdd65WVVWVJCkSiaiyMjSkvLWvydpNuvji8zVhQp0uvfQKTZ48WQ888Cs9+uifVVkZ1DHHHKtPfeqMIfd7+eWNWrDgCJWVlWmPPfZQJNKrlpYWWfua3nzzn3rqqSc1bdp0XXLJ5aquHpvwPCJzTgXemrWb9djzbw/83twaHvj9rKWz/aoWMKKdO3fq1lt/rO3bW3TeeZ/T0Ucv0ne+c4PuuONu1dVN1KpVd+iPf/yDjDlAS5d+RIsWLVFTU6Muvvh8nXLKaZKkpUuP16JFH9aDD/5aBxxwkL70pUv00ksb1NHRrscff1zvvvuO7rrrZ4pEIrrwwnM1b94CSdIee+yp5cuv1kMP/VYPPfQbXXHFN7R163tavfoXGj9+wpB63nzzbWm/po6OjiEf4tXV1XrnnS1Dytx6680688zPqL5+ypDj8+bNU2Nj25BjkyZN3i3spFjX50c+coLmzDlMt9/+A/3+9w/qjDN2tZzOP/9LSetZXl6uiRMnSZIeeOBX2rlzpxYsWDikzD77zJAxB2jBgoV65JGHddttN+vcc7+ov/zlUd1++92qr6/Vpz/9WS1ceIT23nvGoHPQPuQcVlePVUdHuw444CB9/OMna//9D9A99/xUq1ev0sUXX5q0nkifM4EX7oloQ0Njwts2NDTp1EUz6d5EwZkz57CBD97a2nFqampUc3OTrrkm1noKh8NasGChjjzyKP361/fpyScfV3X12CHdgHvvvY8k6eMfP0n33nuPLr/8yxo7tkYXXHCRXn/9dR166ByVlZUpGAzqoIM+oH/84w1J0n77GUnSlCnv0yuvvCRJGj9+wm5hJ6Vu4T344Bo9/vhfJEkrVlynnTt3le3s7FRNTe3A701NjXrppQ16++23tHr1XWpt3aGVK69KGGrJBINBzZlzmCTp4IMP1bPPrh9ye6oWniT19fXp9tt/qLfe+qduuOFmlZUNXYh/3rwFCoViLcBjjvmw7r77Tr3xxut6771/6ZJLLlRlZVA7duzQW2+9pR/96DZ1dnZq5sxZmjZtb3V2dgw6Bx2qqanVMcd8WLW1tQOPd9tt383oNSM5ZwJvR3tY21rDCW9raevSjvawptRV57lWQHLWbpIkbdvWrI6ODtXXT9GUKVN0003fV01NjZ566kmNGVOtX/3qFzr44EN0yimn6cUXn9czzzw18BjxQRZPPfWkDj10rs4553w9+uifdO+99+gTn/iYfvnLNVq27N/V29urV199WSec8HFJT+/24T74sYZL1cI79dRlOvXUZQO/B4MV2rLlbe2111Q999wzOvvsXYNWJk+u1y9/+ZuB3z/xieMzDjtJ6u3t1d/+1qD99putV155aWBgSVyqFp4kffe731ZFRYVuvPF7CV/7TTddr0WLlujYY4/T888/J2MO0N5776MZM/bV9773Q02ZMk4/+tGdmjlzvyHnaNOm13THHT/UmWd+Rlu3blVfX1QTJkzQeed9TpdddoUOPPBgvfDCczJm/92eM1tPbIy1phfPmer5Yxc6ZwJvfE1IE8eF1Jwg9OpqqzS+JpTgXoC/tm1r1iWXXKj29nZdfvmVCgQCuuSSr+mKKy5RNBpVdfVYXXPNdSorK9Ott96sv/zlEdXU1CgQCKi7u3vIY+2//4G6/vqVuueen6qvr09f/vJXdfTRC/TEE0/pggvOVk9Pj5YsWZqTD9nhvva1q3TddSvU19enBQsW6qCDDlZr6w7ddNP1+va3M2vVjHQNT5Luvfcevffev/S+9+2h8867MKPHtXaT/uu/fq9DD52rr3zli5KkT33qTM2de9hAPb/4xYt1443/od/+9n6NGTNGV155jSZPnqz58xfoS186V319Ee233wGqr68f8tj773+ADjlkji644GxFo1F99atXDpyX2267WYFAUJMmTdLy5VdnVGckVxaNFu9egJluZHjfYw1DruHFLZ0/bdTX8Orra3e7plDIqG9ueVnfP/7xD/rnP/+hCy/8siePl0gpnN/e3l7dccd/6stfvmzI8dNOO1H33vuAQiH/vszm6/xmsgHs/Y9uikql28JLdi6caeFJ0rIlsyTFrtm1tHWprrZKc2dPHjgOoDidddZnUheC85wKvEB5uc5aOlunLprJPDwUvI9+9ES/q1AUgsGgJk2avNvxBx74gw+1QSFzKvDiQhUBBqgAgGOYbQ0AcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwAoEHAHACgQcAcAKBBwBwQjCdQsaYKZJekHScpF5JP5MUlfSqpIustX3GmJWSPtZ/+6XW2ueMMbOyLevVCwUAuC1lC88YUyHpJ5J29h/6vqQV1toPSSqTdJIx5jBJiyQtlHSGpB97UTb7lwcAQEw6XZq3SLpT0jv9v8+T9GT/zw9LWirpaEmPWGuj1to3JQWNMfUelAUAwBNJuzSNMZ+X1Git/bMx5qr+w2XW2mj/z22SxksaJ6l50F3jx7Mtm1RdXbWCwUCqYjlXX1/rdxUyQn1zi/rmFvXNztjqSpWXlxdcvfIh1TW8cyRFjTFLJc2R9HNJUwbdXitpu6TW/p+HH+/LsmxSLS2dqYrkXH19rRob2/yuRtqob25R39yiviM/T7o6OrslqajOYyaSnYukXZrW2mOstYustYslbZT0WUkPG2MW9xc5QdI6SeslHW+MKTfG7C2p3FrbJGlDlmUBAPBEWqM0h7lc0ipjTKWk1yQ9YK2NGGPWSXpGsRC9yIuyo31RAAAMVxaNRlOXKlCNjW2+V54ultyivrlFfXMrj12aZemWvf/RTVFJWjxnau4q5KNk54KJ5wAAJxB4AAAnEHgAACcQeAAAJxB4AAAnEHgAACcQeAAAJxB4AAAnEHgAACcQeAAAJxB4AAAnEHgAACcQeAAAJxB4AAAnEHgAACcQeAAAJxB4AAAnEHgAACcQeAAAJxB4AAAnEHgFLtwT0daWToV7In5XBQCKWtDvCiCxSF+f1qzdrA0NjdrWGtbEcSHNnV2vZUtmKVDO9xQAyBSBV6DWrN2sx55/e+D35tbwwO9nLZ3tV7UAoGjRVChA4Z6INjQ0JrxtQ0MT3ZsAMAoEXgHa0R7WttZwwtta2rq0oz3xbQCAkRF4BWh8TUgTx4US3lZXW6XxNYlvAwCMjMArQKGKgObOrk9429zZkxWqCOS5RgBQ/Bi0UqCWLZklKXbNrqWtS3W1VZo7e/LAcQBAZgi8AhUoL9dZS2fr1EUztaM9rPE1IVp2AJAFAq/AhSoCmlJX7Xc1AKDocQ0PAOAEAg8A4AQCDwDgBAIPAOAEAg8A4AQCDwDgBAIPAOAEAg8A4AQCDwDgBAIPAOAEAg8A4AQCDwDgBAIPAOAEAg8A4AQCDwB8Fu6JaGtLp8I9Eb+rUtLYDw8AfBLp69OatZu1oaFR21rDmjgupLmz67VsySwFymmPeI3AAwCfrFm7WY89//bA782t4YHfz1o6269qlSy+QgCAD8I9EW1oaEx424aGJro3c4DAAwAf7GgPa1trOOFtLW1d2tGe+DaMHoEHAD4YXxPSxHGhhLfV1VZpfE3i2zB6BB4A+CBUEdDc2fUJb5s7e7JCFYE816j0MWgFAHyybMksSbFrdi1tXaqrrdLc2ZMHjsNbBB4A+CRQXq6zls7WqYtmakd7WONrQrTscojAAwCfhSoCmlJX7Xc1Sh7X8AAATqCFBwAOemLjloGfF8+Z6mNN8ocWHgDACQQeAMAJBB4AwAkpr+EZYwKSVkkykqKSviipS9LP+n9/VdJF1to+Y8xKSR+T1CvpUmvtc8aYWdmW9e7lAgBclU4L70RJstYeJWmFpBskfV/SCmvthySVSTrJGHOYpEWSFko6Q9KP+++fVdmsXyEAAEoj8Ky1v5N0fv+v+0jaLmmepCf7jz0saamkoyU9Yq2NWmvflBQ0xtR7UBbIGzbiBEpXWtMSrLW9xph7JJ0i6TRJx1lro/03t0kaL2mcpOZBd4sfL8uy7Ijq6qoVDPq/KkF9fa3fVcgI9d1dJNKn1X/4f3r21XfVuH2n6ieM0REH76lzTjxIgUBml7o5v7lFfbMztrpS5cM2ly20OuZK2vPwrLWfM8ZcKem/JY0ZdFOtYq2+1v6fhx/vy7LsiFpaOtOtfs7U19eqsbHN72qkjfomdt9jDUM24tzaslMPrXtDnTu7M9qIk/ObW9R35OdJV0dn927HiumcppLsXKT86mqM+Ywx5qr+XzsVC6XnjTGL+4+dIGmdpPWSjjfGlBtj9pZUbq1tkrQhy7JATrERJ+CGdFp4v5H0v4wx/1dShaRLJb0maZUxprL/5westRFjzDpJzygWpBf13//ybMp68SKBZNLZiJN1DoHilzLwrLUdkk5PcNOiBGWvlXTtsGMN2ZYFcim+EWdzgtBjI06gdDDxHM5jI07ADSweDYiNOAEXEHiA2IgTcAGBBwzCRpxA6eIaHgDACQQeAMAJBB6KCmtdAhgtruGhKET6+rRm7WZtaGjUttawJo4Lae7sei1bMkuBcr63AUiNwENRWLN285C1LptbwwO/Z7LWJQB38dUYBY+1LgF4gcBDwUtnrUsASIXAQ8GLr3WZCGtdIhMMenIb1/BQ8OJrXQ6+hhfHWpdIB4OeIBF4KBKsdYlsMOgJEoGHIsFalxitVIOeTl00k/eSI2jLo6jE17rkAwrpYtAT4gg8ACWNQU+II/AAlDQ2+EUc1/AAlDwGPSW2eM5Uv6uQVwQegJLHoCdIBB4Ah7DBr9u4hgcAcAKBBwBwAoEHAHACgQcAcAKBBzgq3BPRu00d7BwAZzBKE3DMkJ0D2sKaWMvOAXADgQc4hp0D4Cq+zgEOSbVzAN2bKGUEHuAQdg6Aywg8wCHsHACXEXiAQ9g5AC5j0ArgGHYOgKsIPMAxg3cOCFRWKNLdQ8sOTqBLE3BUqCKgPSePJezgDAIPAOAEAg/wUbgnoq0tncx/A/KAa3iAD4Ys79Ua1sRxLO8F5BqBB/iA5b2A/OOrJJBnLO8F+IPAA/KM5b0AfxB4QJ6xvBfgDwIPyDOW9wL8waAVwAcs7wXkH4EH+GDw8l472sMaXxOiZQfkGIEH+ChUEdCUumq/qwE4gWt4AAAn0MIDAEc9sXHLiLctnjM1jzXJD1p4AAAnEHgAACcQeAAAJxB4AAAnEHgAACcQeACQBjbrLX5MSwCAJNist3Tw1wIcMNrWCa2aXZv1NreGFdWuzXrXrN3sd9WQIVp4QAkbbeuEVk1Mqs16T100kzVQiwiBB5SweOskLt46kaSzls72/H6lJp3NelkLtXi481UNcEyq1slI3ZSjvV8pYrPe0pK0hWeMqZC0WtIMSSFJ10v6q6SfSYpKelXSRdbaPmPMSkkfk9Qr6VJr7XPGmFnZlvX01eZIuCfCFi8oOOm0TqaN8n6utGrim/UObu3GsVlv8UnVwvu0pGZr7Yck/ZukH0n6vqQV/cfKJJ1kjDlM0iJJCyWdIenH/ffPqqw3LzF3In19WvW7V7Ri1bO66ifPasWqZ3XfYw2K9BVFTnuOAQ6FZbStE1o1Qy1bMktL50/TpHFVKi+TJo2r0tL509istwiluoZ3v6QH+n8uU6xFNk/Sk/3HHpb0EUlW0iPW2qikN40xQWNMvQdlf5vl68sprnPEMMChMI22dUKrZig26y0dSQPPWtsuScaYWsWCb4WkW/rDSpLaJI2XNE5S86C7xo+XZVk2qbq6agWD/rzxurp79fLrzQlve/n1Zl1w6hhVVRbmmKD6+lpPH2/V715JGPzVYyp13skfyPrxva5vrhVSfS8+fa6qx1Tq2VffVdP2nZo8YYyOOHhPnXPiQQoEYl9GEtU3nfv5xc/zm6gLOJVCej9I0tjqSpWn8UW00OrthZSfyMaY6Yq1tG631t5njLl50M21krZLau3/efjxvizLJtXS0pmqSM5sbelUY8vOhLc1bd+p1//RXJDXOerra9XY2ObZ44V7Ilr/UuI9tda/9I5OOHx6Vt+Gva5vrhVifU8+aoZOOHz6kNbJtm0dkpLXN9n9/FKI5zeZfNU3k3Dq6OxOq1wxnefBkp2LpDFvjHmfpEckXWmtXd1/eIMxZnH/zydIWidpvaTjjTHlxpi9JZVba5s8KFuwuM4Rk84AB/gvVBHQlLrqjL98jPZ+QCFK1cL7hqQ6SdcYY67pP3aJpB8aYyolvSbpAWttxBizTtIzioXoRf1lL5e0arRlPXmFOcJ1jph48DcnCD2Xgh9A4SuLRqOpSxWoxsY2Xysf6evTH555U+tfekctbV2qq63S3NmTC3qwRi66WO57rCFh8C+dPy3rwTt0YeUW9c2tPHZplqVb9v5HN6X1ubl4ztTRV8hHyc5FYY6qKBKB8nKdd/IHdrvO4Zr48OwNDU27BT8AFAoCzwPx6xyuyvWwbSb2A/ACgQfPeB388Yn961/awvw+AFkj8FCwmNgPwEt8TUZBYgFjAF4j8FCQmN8HwGsEHgoSE/sBeI3AQ0GKT+xPxKWJ/QC8w6AVFKxlS2apekxlwon9AJApAg++STW/jon9ALxE4CHvMt0/z/WJ/YAfntiYeBeUwYpt+TECD3nH/DoAfmDQCvKK+XUA/ELgIa+YXwfALwQe8or5dQD8QuAhr5hfl1vhnoi2tnTSNQwkwKCVElBs2+ewf573Mh35CriIwCtixfohl+v981zEyFcgtcL9VERK8Q+55tawotr1Ibdm7Wa/q5aW+Py60YYd3XcxjHwF0kMLr0il+pA7ddHMkm01FWvLNlfSGfnKxH2AFl7Rcnl4f7G3bL3GyFcgPQRekXL1Q47uu90x8hVID4FXpFz9kHO5ZZvMsiWztHT+NE0aV6XyMmnSuCotnT+Nka/AIFzDK2IuDu+Pt2ybE4ReKbdsU8nXyNdimwIDDEbgFTEXh/fHW7aDh+DHlXLLNl252lnCy4FChCb8QuCVANe2zynGlm2xf8h7Mc+P0bXwG4GHolNMLdti/5AP90TU2NLpyRQYJsfDbwQeilYxtGy9/JDPZytxcFAnul4al+48P5fnjaJwEHhwXq6CxKsPeT9aicODeiTpDhRicjwKAYEHZ+U6SLz6kM93V2CyoB4u3YFCjK5FISj8iwhAjuR6xZZkiwNUVgRUU12Z8jH8mGifLKglqWwU8/xcnTeKwkLgwUn5CJJkH/Jd3RH9bt0bKR/Dj4n2yYJ6Ym1I151zuK4/b6HOWjo7o5Ywk+PhN7o04aR8XVM6+UPv11Mvv6uu7t0DNJ3reH50BSab63iYqde0+ppRPW4xja5FaaKFByflay3S9s4ehROEnZReC82vrsBctsay3RYKGC1aeHBSrldsiY/8HBMKZt1C82OiPa0xlCICD87KRZBEIn2677GGISM/q6sqEgZeusHqZ/gUw1xHIF0EHpyViyBZ/Yf/t9sUgubWsKZPqVFnV29WwUr4ANkh8OA8r4Ik3BPRs6++m/C2zq5effPz87Uz3Ev3IOATBq0ga+GeiLa2dDq5+epgO9rDaty+M+FtLW1d2hnudX6wBu8V+IkWHkat2BdG9tr4mpDqJ4zR1pbdQ8/11UR4r6AQEHgYNVa/HypUEdARB++phxJMKHd9NRHeK6XpiY1b0i67eM7UHNYkPXy1wqj4seRVMTjnxINYTWQY3isoFLTwMCqsfp9YIMD8teF4r6BQ0MLDqCRbqaSmukKB8rI816iwsJrILvla1QZIhcDDqCRb8qq1o0dX3vmMVq5+Tt29vXmuGQoNOyWgUBB4GLXB6y0O1xeV3trarht+/qIPNSss+R6KX4hD/9kpAYWAa3gYtfhKJccvmK7ldzyjaIIyWxrb1dbZrdo09n4rNfkeil/IQ/9ZmxOFgBYesra1ZWfCsJNiLb23t7bntT6FItcbzPr9fKORz2ubhdjShb9o4SFr06bUqLwsFm7DlZfFbndNqqH4qfbBK/TnK2SF3NKFv/jrI2u11ZWaOsKmoFPra5zszsz3TuV+7IxeqIqhpQt/EHjwxNWfPUzT+1t6UqxlN31Kja7+7GH+Vswn+R6Kz9D/GCa5Ixm6NOGJymBQ151zuNo6u/X21nZNm1J4Lbv4pqz5GDCRbIPZ6qqgggFv5ynmekPbYsEkdyRD4MFTtdWVOmDGRL+rMYRf13SWLZkl++Z2vTVs0M5bW9u1Zu1mz9eQzHRD23BPRO82dSjSEymZQIy3dLPZYR6li8BDyfNr4eLeSFSdXT0Jb8vFQJJ0h/4P+QLQFtbE2tIZ1EFLF8kU97sbSMHPazp+DSRJNfR/yKCOaOkN6mCSO0ZCCw8lzc9rOoXYvZbJ9IV8XvP0EpPcMRICDyXNz9ApxO61dL4ATBpfVRLz2OItXSCueN69wCj4vXDx4O61MkkTair14bl7+da9ls70BeaxoVQReCh5ub6mk2wJq95IVB+eO1UH7ztRE2pC2tHerZdfb9aatZsV6evL6LG8kOoLgCTmsaFkpdWlaYxZKOk71trFxphZkn4mKSrpVUkXWWv7jDErJX1MUq+kS621z3lR1ruXClfl6ppOoukORx06VSceubckDdw2vDs10SjRfE6dSDZ9oXlHF/PYULJSBp4xZrmkz0jq6D/0fUkrrLVPGGPulHSSMeafkhZJWihpuqQHJS3Itqyk33r3UuE6r6/pJJru8NC6N9S5s1uSEl67G2zwIJF8Tp0Y/AUgUFmhSHfPwBeAQhxoA3glna+Or0v65KDf50l6sv/nhyUtlXS0pEestVFr7ZuSgsaYeg/KAgUp2WjHF23jiLcNFm8x+TV1IlQR0J6Txw5p7fp9zRPIpZQtPGvtg8aYGYMOlVlr4+vit0kaL2mcpOZBZeLHsy2bVF1dtYJB//8D1tfX+l2FjFDf7L3b1KFtbYm7/ra1hZXOwmGTJ4zRzBmT1NIaHvGxWtq6FKisUP3ksVnUNrnh5/fi0+eqekylnn31XTVt36nJE8boiIP31DknHqRAwP/L/oX4fkim0Oo7trpS5T6Mti2E8zCfldZhAAAQt0lEQVSaaQmDr6vVStouqbX/5+HHsy2bVEtLZyb1zon6+lo1Nrb5XY20UV9vRHoimlibuOtPkioryhXuSX4J+pCZk9S2Y2fSx6qrrVKkuyetczCaeXMjnd+Tj5qhEw6fPuTxtm3rSPAI+VWo74eR5Ku+mYRJR2d3Dmsysnz93ZKdi9HE/AZjzOL+n0+QtE7SeknHG2PKjTF7Syq31jZ5UBYoSKGKgA6ZOWnE28vKRm7jDR8lmm03YqSvT/c91qAVq57VVT95VitWPav7HmtIOAo0EyOt2MLGqihWo2nhXS5plTGmUtJrkh6w1kaMMeskPaNYiF7kRdnRviggLperhSydP12Pb3gn8fN2R3TUwXto05vbB0ZCHjJzopbOn66J46p2q0umCz8Plq8BL2ysimJXFo0m2Ka6SDQ2tvleebpYcmu09c3Hh3O4J6IVq55N2BU5aVyVrj9voSRlFLiZBnRbZ7dWrn5O29t376aK1yHZ42Ryfu97rCHhyNOl86fldBHuwVx5/47iedLeb+r+Rzf58rm5eM7UvDxPsnPB1zKUpHysFpJOV2SqhZwTPWY65ePdmNeu/p+EYSd5u0A1G6uiFBB4KDn5/HAevorLlLoxeVmZPx7oLUkCzct5c37t/AB4icWjUXLyuUPC8FVcZs6IjbzMpWSBPpiX8+aYkI5SQAsPJSedBZK9Fu+KrKrM/XfIZIEuxRao9rqVyYR0ZOuJjVsG/vmFFh5KTiFuy+OlpK2tmpCuPWeBaqsrPX/ebEaSAoWAwENJKuUP52SBPm//+pyEncTGqih+BB5KUql/OPsZ6GysimJF4KGkleqHc6kHOpALDFqBc0ppaaxM5/kBLqOFB2cU+tJYuVgGLZdLqwHFhsCDM/K15mSmIZOLIC70cAf8QODBCalWX4nvPJ6NSF+fVv3uFa1/aUtGIZOLIM7nDupAseCrHpyQj6Wx1qzdrIfWvZHR+p25WAYt3ccM90T0blNHSVzLBNJBCw9OyPXSWKNtQWazDNpIXaepHnNba5ce37Al1t3ZFtbEWro74QYCD07I9eorow2u0QRxqutzqR7zsRfe1uMv7lreie5OuIKvc3DG8J0Nhu88no3Rrt85mjUqU219FKoI6ND9Jid8zINn1unlzU0Jb0u3C7WUpnXALbTw4IxcTtbOpgWZyaop6XadjrQDZndP36i7UBn5iWJH4ME5uVp9ZdmSWaoeU6n1L72T0XJfmQRxOl2n42tC2vi3xK04+8/to76WychPFDsCD/BIoLxc5538AZ1w+HTtaA9rTCioneFe9UaiCqTRAEoniNO55pcsFLe3h3XkQXto/av/2u22ZC3RfEzrAHKNwAM8FgyU6bEX3s5J1186XadjQkFNqAkl3A29rrZKZx43W2OqghktPJ3PTXWBXCHwAI/luutvpGt+py3eV/c91qANDY0Jw06KhWJ1KDjQhRqorFCkuydl64wdz1EKCDzAQ13dvTnv+hvpmt99jzUkbPlJsRGpw1txoYqA6iePVWNjW8rnzGRQjtfrdw5+PCAbBB7goZbW/HX9Db7m19bZrRc2JQ7aCTWV+ubn56uyIqDmHV2jDqJUo0m9HsWZ6PGOOnSqTjxyb0aFYlQIPMBDdePy2/UXD4XnN23V9vbuhGV2tHfr12s3a9ObLbsFkZR+iyzVaFKvu3ITPd5D695Q585uRoViVAg8wENVlcGcrugy3PBQSKSyonzIqMx4EPVFo6qpDmW82HWi0aRej+JkVChygcADPJbJRPJsJAuFoeX6Eh5/+pV/qat712op2bTIvB7FyahQ5AKBB3gslyu6xIV7Inpjy46EXadxoWC5wr2Jw07SkLAbbDQtKK9HcTIqFLlA4AE5kosVXYYP5Cgvk/qiu5cbP7ZC5eXlCrdlvu3RaFpQXi/OnevFvuEmAg8oIsOv2UUThJ0kHfT+SXomwWoq6RhtC8rrrtxEj3fUoXvpxCP3HtXjAQQeUCTSuWZXXiZNra/RGcfOlH2zJWmX50hG24Lyuis30eNN22tCWvMGgUSYzAIUiWQDOeL6otJbW9v10Pp/jrjt0EjKy6TpU2p02uJ9s6nmQFeulztRePl4cBeBBxSJZHvuDbehoUknf2jf/v3/0rtPPCwfeOKNbKoJFCy6NIEiEaoIaM5+k/WXF7akLNvS1qX2zu6BLsH//Werp9O8psc8N+TaExtTv4czsXjO1LTK0cIDisgIY1R2M3jgSagioLOO209HHvQ+hSpS/5ePj9IESg0tPKBIhHsiemmEjV2Hiw88Sbi+5X4TtGTeVN35u1e1rW335chGGqXp9aLQQL4ReECRSDVopUzSxGG7IiRaj7L5r++pprpCh5kpac1z83pRaMAvBB5QJJKtPjJpXEiXnHaI6geNZky1HuV15y5Q9ZhKrX/pnaTz5nK9vx+QLwQeUCSSrz5Sr2lTahXuiWhrS6fG14RSrkfZ3tmj807+gE44fPqIXZUs4oxSQuABRSSd3c7j3Y6HzJw0YotwQk1IY0JBvdvUIUkjLiPGIs4oJQQeUETS3e28uTWsxze8o+lTahIG3o6Obl1559Pq6u7TpCTX5LJdxJmBLigkBB5Q5JJ1O7Z2dOuYOXvqub9uHbI7QqQvqkh3bJJD/JpcpC+q4xdMHxJOo13EmYEuKEQEHlBEEgXJPnvUjrhm5o6Obr30tyb19CbeCmiwJzds0eMvbtmtxTeaRaEZ6IJCROABRSThNIMU62vu6OhJ67Hj2wwND6dMF4VmoAsKFX0LQJFId4dzr2xoaFK4Z1fLMN1FnNMZ6AL4gcADikQ6uyWUefh8LW1daty+U1tbOocEXyrJFrlmt3L4iS5NoEgkGzEZF5U0rrpCrZ3pdWMmU1kR0G2/3qiWtu6MBp2wWzkKFS08oEjEgySZSeOqNG//KZ48X1d3RNvauhXVrut6a9ZuTuu+y5bM6t+aqErlZbF6LZ0/bdS7nwNeoIUHFJF4YDz18rtDphnExUdPBsrLtKGhSdvauhRNd4sFSWVl0sTakDq6etTV3bfb7ekOOvF693PAC7TwgCISD5JbLvqgPnjwHppYG9qtBdUbiWrpvGn65ufn68iD9kj7setqKnXdOYfrktMOUThB2EmZDzpht3IUElp4QBGqDlXoCx8/cMhKJsFA2ZA5ehNqKtTSnv61vJ5In/acVK3eSDSr1VWAQkULDyhig1tQ8Tl6za1hRaWMwk6S2nf26r7H/pb0WiGDTlDMCDygBHg1R29j/9w7Bp2gFNGlCZSAdObopWN7R3hgBwQGnaDU0MIDSkCyyd6ZmDjsGh2DTlBKCDygBKQzRy8dXKNDKaNLEygRw3c1CAbK1N2b/iS8hQdO4RodShotPKBExOfoXX/eQl17zuGqGVOR0f03/q1Ra9ZuVqQv8Rw8oNgReECJCVUEVBksV0tbd0b3C/dEM1o+DCg2BB5QgrIZxDJ8WyCgVBTUNTxjTLmk2yUdKiks6QvWWr5uAhlKtmNBKvHlw6bUVeegZoB/Cq2Fd7KkKmvtkZK+Lul7PtcHKFrxyePjx2Z2LY/lw1CqCi3wjpb0J0my1j4rab6/1QGKV3wQy3+cu1ChYPr/1ZmagFJVUF2aksZJ2jHo94gxJmit7U1UuK6uWsGg//8x6+tr/a5CRqhvbhVafeslHXfEPvqvp/6esuyx86fr4tPnKBAotO/CuxTa+U2l0Op7yof3K4jPTT8UWuC1Shr87igfKewkqaWlM/c1SqG+vlaNjW1+VyNt1De3CrW+J31wH3V19eiFTe+NuKj0pHEhnbZoX23b1pHn2qWvUM/vSPJV30xCtRA+N3Mp2bkotK9x6yV9VJKMMUdIesXf6gClId69+e0LPqijDk68R97c2fV0ZaKkFVoL77eSjjPGPC2pTNLZPtcHKCmhioA+/9H9NaYqOLAiS11t1cBO6UApK6jAs9b2Sfqi3/UASlm8tXfqopkKVFYo0t1Dyw5OKLQuTQB5EqoIaM/JYwk7OIPAAwA4gcADADiBwAMAOIHAAwA4gcADADiBwAMAOIHAAwA4gcADADiBwAMAOKEsGo36XQcAAHKOFh4AwAkEHgDACQQeAMAJBB4AwAkEHgDACQQeAMAJBB4AwAlBvytQKIwx4yX9QtI4SZWSvmqtfcYYc4SkH0jqlfSItfY6Y0y5pNslHSopLOkL1trN2ZYdZb1PkfQpa+1Zg36/RdJb/UVWSlpXwPUt6PPbX8cySW9L+lv/oWestVcZY06U9M3+51htrV1ljBmj2PtoiqQ2SZ+z1jZmUna09czwNSU8b/l47mH1eFFSa/+vf5f0E+Xg/ZBlHRdK+o61drExZpakn0mKSnpV0kXW2j5jzEpJH+uvy6XW2ue8KJtNvbE7Wni7fFXSX6y1iyR9XtKP+4/fKeksSUdLWmiMmSvpZElV1tojJX1d0vc8KpsRY8wPJN2ooX/HeZKWW2sX9/97ssDrW7Dnd5CZkl4cdE6vMsZUSLpV0kckLZJ0vjHmfZIulPSKtfZDkn4uaUUmZbOoY6ZGOm95Y4ypklQ26Lyerdy9H0Zbx+WS7pZU1X/o+5JW9P/NyiSdZIw5TLG/60JJZ2jXZ0dWZbOpNxIj8Ha5VbFvl1Ks5dtljBknKWStfd1aG5X0Z0lLFfsP9idJstY+K2m+R2Uz9bRiH5qDzZN0jjFmnTHme8aYYKHWtwjOb9w8SVONMY8bY/5ojDGSDpC02VrbYq3tlvSUpGMG10fSw/3Pm0nZfNntvOXxueMOlVRtjHnEGLPWGHOMcvB+yLKOr0v65KDf50l6sv/n+N/saMVamFFr7ZuSgsaYeg/KwmNOdmkaY86VdNmww2dba//HGLOHYt1MlyrWvdk6qEybpH37j+8YdDziUdlM67vGGLN42PFHJf1Ose6hOyV9sYDrWxDnN426XyTpRmvt/caYoxV7f1w27HnbJI0fVp9Ex1KVzZfdzpsxJmit7c1jHToV636/W9J+in3Qbx90uyfvh2xel7X2QWPMjEGHyvoDNv6c8b9j87C6jPegLDzmZOBZa38q6afDjxtjPiDpV5K+Zq19sv9bZO2gIrWK/YesHna8XLH/fNmWzai+I1htrd3e/3p+L+lUxT4ACrG+XtQh6/qmqrsxplqx6y2y1j5ljNlLsQ+lRM8x+LkTHUtVNl+G16k8z2EnSQ2KtXyjkhqMMTskTRx0uyfvB49f1+Draqn+vtmWhcfo0uxnjDlQ0v2SzrLWPixJ1tpWSd3GmJn9AxeOV2wAyHpJH+2/3xGKXYfxomy2r6FM0svGmGn9h46V9EKh1reIzu9KxVr8MsYcqtiAoL9K2s8YM9EYU6lYF+Uzg+sj6YT+530tg7L5stt5y+Nzx52j/mts/V8iqiV1eP1+8LjOGwb1UsT/ZuslHW+MKTfG7K1YyDZ5UBYec7KFN4IbFbsw/YPYJRrtsNaepFiX4L2SAor1vf+3MeZ/JB1njHlasQvMZ/c/RlZls30B1tqoMeYLkn5jjNmp2IfyKsW6gAquvl7UIU/1vUnSL4wx8ZF1n7fW9hhjvqrYtaNyxVrWW4wxd0i6xxjzlKRuxb5ApV02izpm6rdKfN7y6aeSftb/+qOKBWCfcvN+8Mrlklb1f3F5TdID1tqIMWadYl9iyhXrAs+6rMf1htgeCADgCLo0AQBOIPAAAE4g8AAATiDwAABOIPAAAE4g8AAATiDwAABO+P/mzLe5BibIUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ddb79b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.jointplot(population_a - population_b, population_a + population_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[8, 19]\n",
      "[19, 8]\n",
      "[19, 8]\n",
      "[19, 8]\n",
      "[8, 19]\n",
      "[19, 8]\n",
      "[19, 8]\n",
      "[8, 19]\n",
      "[8, 19]\n",
      "[19, 8]\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "def neighbors(idx):\n",
    "    x, y = divmod(idx, 10)\n",
    "    ret = []\n",
    "    if x: ret.append(idx - 10)\n",
    "    if y: ret.append(idx - 1)\n",
    "    if y < 9: ret.append(idx + 1)\n",
    "    if x < 9: ret.append(idx + 10)\n",
    "    random.shuffle(ret)\n",
    "    return ret\n",
    "\n",
    "for i in range(10):\n",
    "    print(neighbors(9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ideas\n",
    "\n",
    "- Genetic Algorithm\n",
    "- Simulated Annealing (needs a fitness function though)\n",
    "- Some kind of incremental/greedy approach, possibly starting from a seed\n",
    "- Stupid approaches such as 10 2x districts and 10 virtually invisible districts\n",
    "- Do it by hand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's try the stupid approach!\n",
    "\n",
    "First, we need to use some quick calculations to see how possible it is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.7272529393917373"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_pop = np.sum(population)\n",
    "sos = ((total_pop / 20) ** 2) * 20\n",
    "sos / 164e10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "472900.34999999998"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "total_pop / 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, nevermind, maybe not so easy..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "0.0\n",
      "0.0\n",
      "0.999999854587\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "1.0\n",
      "0.999999537809\n",
      "1.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "5.25039667882e-05\n"
     ]
    }
   ],
   "source": [
    "for district in naive:\n",
    "    na = sum(population_a[district])\n",
    "    nb = sum(population_b[district])\n",
    "    print(win_probability(na, nb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulated Annealing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.4999481036372178"
      ]
     },
     "execution_count": 126,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def loss(soln):\n",
    "    # This is completely arbitrary\n",
    "    ret = 0\n",
    "    x = expected_districts_won(soln)\n",
    "    y = district_population_imbalance(soln)\n",
    "    z = expected_efficiency_gap(soln)\n",
    "    if x < 10:\n",
    "        ret += 10.5 - x # screw you 0.99999\n",
    "    if y > 164e10 * .99:\n",
    "        ret += (y / 164e10 - 1) * 3\n",
    "    if z < -0.07:\n",
    "        ret += (-0.074 - z) * 2\n",
    "    return ret\n",
    "\n",
    "loss(naive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "def connected(d):\n",
    "    if not len(d):\n",
    "        # We can't have empty districts...\n",
    "        return False\n",
    "    vis = { city: False for city in d }\n",
    "    def dfs(city):\n",
    "        if vis[city]:\n",
    "            return\n",
    "        vis[city] = True\n",
    "        for n in neighbors(city):\n",
    "            if n in d:\n",
    "                dfs(n)\n",
    "    dfs(d[0])\n",
    "    return all(vis.values())\n",
    "\n",
    "def valid(soln):\n",
    "    return all(connected(district) for district in soln)\n",
    "\n",
    "assert connected([0, 1, 2])\n",
    "assert connected([2, 1, 3])\n",
    "assert not connected([3, 4, 6])\n",
    "assert not connected([9, 10])\n",
    "assert connected([3, 4, 5, 15, 25, 35, 36])\n",
    "assert not connected([3, 4, 5, 15, 25, 36])\n",
    "assert valid(naive)\n",
    "assert not valid(naive + [[]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 1, 2, 3, 4],\n",
       " [5, 6, 7, 8, 9],\n",
       " [10, 11, 12, 13, 14],\n",
       " [15, 16, 17, 18, 19],\n",
       " [20, 21, 22, 23, 24],\n",
       " [25, 26, 27, 28, 29],\n",
       " [30, 31, 32, 33, 34],\n",
       " [35, 36, 37, 38, 39],\n",
       " [40, 41, 42, 43, 44],\n",
       " [45, 46, 47, 48, 49],\n",
       " [50, 51, 52, 53, 54],\n",
       " [55, 56, 57, 58, 59],\n",
       " [60, 61, 62, 63, 64],\n",
       " [65, 66, 67, 68, 69],\n",
       " [70, 71, 72, 73, 74],\n",
       " [75, 76, 77, 78, 79],\n",
       " [80, 81, 82, 83],\n",
       " [85, 86, 87, 88, 89],\n",
       " [90, 91, 92, 93, 94, 84],\n",
       " [95, 96, 97, 98, 99]]"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import copy\n",
    "\n",
    "def all_neighbors(soln):\n",
    "    dnum = {}\n",
    "    for i, district in enumerate(soln):\n",
    "        for city in district:\n",
    "            dnum[city] = i\n",
    "    for city in range(100):\n",
    "        for n in neighbors(city):\n",
    "            if dnum[city] != dnum[n]:\n",
    "                # Attempt to switch city's district to n's district\n",
    "                soln2 = copy.deepcopy(soln)\n",
    "                soln2[dnum[city]].remove(city)\n",
    "                soln2[dnum[n]].append(city)\n",
    "                if valid(soln2):\n",
    "                    yield soln2\n",
    "\n",
    "def random_neighbor(soln):\n",
    "    dnum = {}\n",
    "    for i, district in enumerate(soln):\n",
    "        for city in district:\n",
    "            dnum[city] = i\n",
    "    while True:\n",
    "        city = random.randint(0, 99)\n",
    "        for n in neighbors(city):\n",
    "            if dnum[city] != dnum[n]:\n",
    "                # Attempt to switch city's district to n's district\n",
    "                soln2 = copy.deepcopy(soln)\n",
    "                soln2[dnum[city]].remove(city)\n",
    "                soln2[dnum[n]].append(city)\n",
    "                if valid(soln2):\n",
    "                    return soln2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hill Climbing\n",
    "\n",
    "Before we bring out the big guns with simulated annealing, let's see how a simpler approach does with less parameters to adjust: hill climbing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1, loss = 4.999948103637218\n",
      "Iteration 2, loss = 4.054051826191645\n",
      "Iteration 3, loss = 4.004477346080161\n",
      "Iteration 4, loss = 4.000000607604006\n",
      "Iteration 5, loss = 4.000000145413409\n",
      "Iteration 6, loss = 4.0\n",
      "Iteration 7, loss = 4.0\n",
      "Iteration 8, loss = 4.0\n",
      "Iteration 9, loss = 4.0\n",
      "Iteration 10, loss = 4.0\n",
      "Iteration 11, loss = 4.0\n",
      "Iteration 12, loss = 4.0\n",
      "Iteration 13, loss = 4.0\n",
      "Iteration 14, loss = 4.0\n",
      "Iteration 15, loss = 4.0\n",
      "Iteration 16, loss = 4.0\n",
      "Iteration 17, loss = 4.0\n",
      "Iteration 18, loss = 4.0\n",
      "Iteration 19, loss = 4.0\n",
      "Iteration 20, loss = 4.0\n"
     ]
    }
   ],
   "source": [
    "soln = copy.deepcopy(naive)\n",
    "\n",
    "for iteration in range(20):\n",
    "    closs = loss(soln)\n",
    "    print('Iteration {}, loss = {}'.format(iteration + 1, closs))\n",
    "    for neighbor in all_neighbors(soln):\n",
    "        nloss = loss(neighbor)\n",
    "        if nloss <= closs:\n",
    "            soln = neighbor\n",
    "            closs = nloss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alright, that clearly didn't work. Back to SA!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1000, loss = 7.6353325936022385, temp = 1.8002\n",
      "Iteration 2000, loss = 9.0120319761963, temp = 1.6002\n",
      "Iteration 3000, loss = 8.05338906700275, temp = 1.4002\n",
      "Iteration 4000, loss = 10.248659007643331, temp = 1.2002000000000002\n",
      "Iteration 5000, loss = 8.954673145385021, temp = 1.0002\n",
      "Iteration 6000, loss = 4.636728910036025, temp = 0.8002\n",
      "Iteration 7000, loss = 5.446503253801343, temp = 0.6002000000000001\n",
      "Iteration 8000, loss = 0.693513817307104, temp = 0.4001999999999999\n",
      "Iteration 9000, loss = 1.1737058650000527, temp = 0.20019999999999993\n",
      "Iteration 10000, loss = -0.026199497391676818, temp = 0.00019999999999997797\n"
     ]
    }
   ],
   "source": [
    "def P(e, e2, t):\n",
    "    if e2 < e:\n",
    "        return 1\n",
    "    return np.exp(-(e2 - e) / t)\n",
    "\n",
    "def accept(soln, soln2, temp):\n",
    "    return np.random.rand() < P(loss(soln), loss(soln2), temp)\n",
    "\n",
    "def temperature(completed):\n",
    "    # Simple linear temperature function\n",
    "    return 2 * (1 - completed)\n",
    "\n",
    "soln = copy.deepcopy(naive)\n",
    "iterations = 10000\n",
    "for iteration in range(iterations):\n",
    "    temp = temperature(iteration / iterations)\n",
    "    if iteration % 1000 == 999:\n",
    "        print('Iteration {}, loss = {}, temp = {}'.format(iteration + 1, loss(soln), temp))\n",
    "    soln2 = random_neighbor(soln)\n",
    "    if accept(soln, soln2, temp):\n",
    "        soln = soln2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10.35391904994836, 1625677608092.55, 0.078676720921604826)"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def properties(soln):\n",
    "    return expected_districts_won(soln), district_population_imbalance(soln), expected_efficiency_gap(soln)\n",
    "\n",
    "properties(soln)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[56, 66, 46, 57, 45, 68, 67],\n",
       " [9, 8],\n",
       " [53, 63],\n",
       " [26, 38, 28, 36, 25, 16, 35, 37, 27, 6, 15],\n",
       " [3, 2, 14, 24, 34, 33, 13, 43],\n",
       " [18, 17, 19, 7],\n",
       " [76],\n",
       " [12, 11, 10, 20, 0, 1],\n",
       " [5, 4],\n",
       " [54, 44],\n",
       " [93, 92, 91, 81, 90, 71, 80, 83],\n",
       " [23, 22, 21, 31, 30, 32],\n",
       " [48, 47, 49, 39, 58],\n",
       " [89, 99, 98, 97, 96],\n",
       " [65, 74, 64, 75, 55, 84],\n",
       " [69, 59, 79, 78, 77],\n",
       " [87, 86, 88, 85, 95, 94],\n",
       " [72, 70, 60, 62, 52, 42, 40, 73, 41, 61, 82],\n",
       " [51, 50],\n",
       " [29]]"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "soln"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# We did it!\n",
    "\n",
    "Totally did not expect that simulated annealing approach to work! I could cleaned the notebook up, but I guess it gives a better feel for the solving process if it's left like this :)\n",
    "\n",
    "This is also the first time I've implemented simulated annealing, so a lot is probably not optimal. I was really surprised by how well it worked though after just a little bit of tweaking, even for a newbie like me who was blindly guessing when it came to tuning the parameters. Hope this was helpful & let me know if you have any suggestions!\n",
    "\n",
    "-Eric"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}