class09_similarity_python3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CSX46 - Class Session 9 - similarity and hierarchical clustering\n",
    "\n",
    "In this class session we are going to hierachically cluster (based on Sorensen-Dice similarity) vertices in a directed graph from a landmark paper on human gene regulation (Neph et al., Cell, volume 150, pages 1274-1286, 2012; see PDF on Canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Pandas `read_csv`, read in the ifle `shared/neph_gene_network.txt`, which has two columns of text (first column is the regulator gene, second column is the target gene), into a data frame. The file has no header and is tab-delimited.  Assign the column names of the dataframe to be `regulator` and `target`, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load the Python packages that we will need for this exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import igraph\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.cluster.hierarchy\n",
    "import scipy.spatial.distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `pandas.read_csv`, read the file `shared/neph_gene_network.txt`; name the two columns of the resulting data frame, `regulator` and `target`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_list_neph = pandas.read_csv(\"shared/neph_gene_network.txt\",\n",
    "                                 sep=\"\\t\", \n",
    "                                 names=[\"regulator\",\"target\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the edge-list data into an undirected `igraph.Graph` object `neph_graph`, using `igraph.Graph.TupleList`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH UN-- 538 47945 -- \\n+ attr: name (v)'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "neph_graph = igraph.Graph.TupleList(edge_list_neph.values.tolist(), directed=False)\n",
    "neph_graph.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the igraph `Graph.similarity_dice()` method, compute a similarity matrix and assign it to name `S`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = neph_graph.similarity_dice()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the `numpy.matrix` constructor, compute a distance matrix (1-S) and assign to object `D`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = 1 - numpy.matrix(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use scipy.spatial.distance.squareform to make a vector-form distance vector from the square-form distance matrix `D`; call the resulting object `vD`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "vD = scipy.spatial.distance.squareform(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `scipy.cluster.hierarchy.linkage` on `vD` (with `method=\"average\")`, perform hierarchical agglomerative clustering.  Assign the resulting object to name `hc`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "hc = scipy.cluster.hierarchy.linkage(vD, method=\"average\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a dendrogram using `scipy.cluster.hierarchy.dendrogram`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.00000000e+01, 1.10000000e+01, 1.78970917e-02, 2.00000000e+00],\n",
       "       [6.40000000e+01, 6.60000000e+01, 1.82926829e-02, 2.00000000e+00],\n",
       "       [6.50000000e+01, 6.70000000e+01, 2.02774813e-02, 2.00000000e+00],\n",
       "       ...,\n",
       "       [3.23000000e+02, 1.07000000e+03, 8.70781231e-01, 5.35000000e+02],\n",
       "       [1.07100000e+03, 1.07200000e+03, 9.06844511e-01, 5.37000000e+02],\n",
       "       [3.71000000e+02, 1.07300000e+03, 9.25471476e-01, 5.38000000e+02]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAHRCAYAAABU26HZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dfbycZX0n/s+VSSaEnOQECJycQAQEfMAHTgChSltO1Vq0rnRb60+z9lf7YLbbtbXdtdTWakF392VX1La/n7v1tNsHa0+turvCKooVDaIoQpoJSJCYBDBhJnMCJJNMOJw7Z+baP77f69z3mczz3DNz3XN/3q9XXjlnzsw91/34vZ4vY60FERERDd+KYSeAiIiIBIMyERGRJxiUiYiIPMGgTERE5AkGZSIiIk+sHNYXb9y40V500UXD+noiIqKB2rlz51PW2nObvWdoQfmiiy7CAw88MKyvJyIiGihjzBOt3sPqayIiIk8wKBMREXmCQZmIiMgTDMpERESeYFAmIiLyBIMyERGRJxiUiYiIPMGgTERE5AkGZSIiIk8wKBMREXmCQZmIiMgTDMpERESeYFAmIiLyBIMyERGRJxiUiYiIPMGgTERE5AkGZSIiIk8wKBMREXliaEH50UeB6Wn5NzMzrFQQERH5Y2hBeX5e/s/lgNnZYaWCiIjIH0MLymvWADt2AFNTw0oBERGRX9imTERE5AkGZSIiIk8wKBMREXmCQZmIiMgTDMpERESeGGpQnpmRIVG5HMcrExERDTUou/HJU1Mcr0xERDT06uupKY5XJiIiAjwIykRERCQYlImIiDzBoExEROQJBmUiIiJPMCgTERF5gkGZiIjIEwzKREREnmBQJiIi8gSDMhERkScYlImIiDzBoExEROQJBmUiIiJPMCgTERF5YuhBuVCQtZS5rjIREaXdymEnoFiUf27pxlxO/t++fXhpIiIiGoahB2UgXFMZkJIyERFRGg29+pqIiIjE0ILyqVNSVR0Ew0oBERGRX4YWlBcXgVIJyGaHlQIiIiK/eNGm7HpgA2FHr+lpYNs2dvgiIqL08CIoR3tgsxc2ERGllRdBGVjeAxtgL2wiIkof73pfz8xwMhEiIkon74Ly7KwE44kJ+f3uu4GbbpLgzABNRESjzJvq6yjXrlwsAtdfH77OdmYiIhplXgZlh+3MRESUJt5VXxMREaWV1yXl6PhlYPkYZodjmYmIaFR4HZRrV5By/ztsYyYiolHiRfV1EEiArdez2rUru3/bti3/O4dNERHRqPAiKGezMg/27Gzr97ohU0D9YVMMzkRElFReBOVOudLz5KT8fv318lou115gJyIi8lEig3LU1NTyKm1WZxMRUVJ53dGrXa5Km4tZEBFRko1EUAaWTzTCSUaIiCiJEl99TURENCoYlImIiDzBoExEROQJBmUiIiJPeNXRK58H5ubkZzcxCBERUVp4VVJ2ARmQOa+JiIjSxKugDMjQptqFJ4iIiNLAu6Dcq5kZmTyk0QIXREREvmorKBtjbjDGPGqM2WeMeW+dvz/PGPMNY8wuY8yDxpg3xJ/U9szOyqpTpZIsUsHATERESdEyKBtjMgA+AeD1AC4H8DZjzOU1b/sjAJ+11m4F8FYA/y3uhHYimwXGx+VnLlBBRERJ0U7v62sA7LPWHgAAY8xnANwIYE/kPRbAev15HEA+zkS2UihIx7BoqZjt0kRElDTtVF+fD+Bg5PdD+lrUzQDebow5BOAOAL9Vb0PGmO3GmAeMMQ9Uq7aL5NY3Odn+esxERES+iquj19sA/K219gIAbwDw98aY07ZtrZ2x1l5trb16xQoT01cTERGNhnaC8pMAtkR+v0Bfi/o1AJ8FAGvtdwCcAWBjHAlsxPWyLhT6+S1ERESD005Qvh/AZcaYi40xWUhHrttr3vMjAK8BAGPMiyFB+UicCa01OytV1pOT/fwWIiKiwWkZlK21iwDeBeBOAI9Aelk/bIz5oDHmTfq2/wjgncaY3QD+EcA7rLXxNRoTERGlQFtzX1tr74B04Iq+9oHIz3sAXBdv0nrjemQDUtW9fftw00NERNTKyM3o5UTnzmavbCIiSoKRDcpAuNIUp9wkIqIkGImg7Hpg1/bGZmmZiIiSZCSCcrEIZDLhz1EsLRMRUVKMRFAGgLGx+lNrsrRMRERJMTJBuZloaXl6miVmIiLyUyqCsistT01JYGaJmYiIfJSKoAxIQN6xg6tHERGRv1ITlImIiHzHoExEROQJBmUiIiJPtDX3tU/cnNauRzUREdGo8C4ou6AbBEA2e/rfJyeBvXsHny4iIqJ+8676uliUdZLrBWQiIqJR5l1QJiIiSiuvg3IQcM5qIiJKD6+DcjYrVdmcgYuIiNLA66BMRESUJqkJyoWCLEaRy7FKnIiI/JSaoFwsSjB2c1+zSpyIiHyTmqAMLJ9whMs4EhGRb7ybPKSfikX550rLuZz8v3378NJERETkpCooA+ESjoCUlImIiHyRqurrqJmZsNMXq7GJiMgHqQ3KrqPX1JQEZnb8IiKiYUttUAbCqmzXxkxERDRMqQ7KREREPkl1UK6dUIRty0RENEyp630dVSwCBw+Gy0Tefffp7cvbtnHIFBERDUaqgzIQBuR67cocx0xERIOUqqDcaCnI6NjlKI5jJiKiQUpVUG62FOTMzOmvu5JybXBmlTYREfVDqoJyM7Oz4YIVhYK0NzsuOAMS1JuNa2bAJiKibjEoR7hq7Onp5XNkt4tt0ERE1AsG5QYatTM3wzZoIiLqBYNym+q1Oddq1AZdi1XcRERUT6onD+mEa3NuZmqqdZU359kmIqJGWFLuQLMq7XZK0o6bPawdLFUTEaUHS8oxaackDbRXmnZYqiYiShevSsrlcjjDVqfcxCBB0P02etVN57Ba0RJ3EADf+hawYUO4fYclaCKi0eNVUK5UJKDOz3f+2SCQz42Px5+uqEbV1M06eXUSQKPjpbNZ2a9G38WgTEQ0WrwKykkQDZpRjaqkuwmg0fHSwOmlbw69IiIaTQzKXeikmjquABotoUdL5azGJiIaHanq6OWqgnM5mUqzls/rK0c7krnOYuwIRkQ0WlJVUg4CIJORn6NzWzvF4vLpNX1ru60tobMam4hotKQqKAPA2FhYyqwnGvg6CXqddgBjtTMREdVKVfV1PzUap1xvXDKrnYmIqJ7UlZT7qd0OYKx2JiKiehJXUnYdtOqN3yUiIkqyxAVl10FrWLN2ERER9UvigjJJpzJfh24REVH32KacQLWzivk2dIuIiLrDknIDMzNhSTTOUmh0u72UcKemZFiVwxIzEVHyMSg3EB2yFOfwJbetOGbkipaYJyaAu+8GbrqJwZmIKKkYlJvoZO3jTre7Y0c823bbmpyUFbI4/SYRUXIxKI+QOIM9EREN3sgG5SAAyuX0jWeOq82aiIgGb2R7XwcBUKkkezxzoSDjsjsJrLOzYUbk7rvl32/+psz5DSwvRXP+bSIiv4xsUO6EC36tStXufRMTg0nX5CSwd2/n7cPZrNQSOJVK+LsbPlUqNW97ZsAmIho8BmVIoC2VpKNUO+9LQputKxk7tWl2GYx6i2gwYBMRDQeDcgrUWyhjenr52tH1uMAd1ShgM1ATEfWOQXlIXIcs9/MwAlqrVa3aCdwAZxQjIooLg/KQ1E5O4mtAa2c5Si5FSUQUD6+HRAWVBQBA7nAOMzvjGdsT2Cpy5RMIbDWW7fWiX5OTEBFRMnldUg4qp5C5+H7gHW/E7ENT2H5V78XJoFrF/K07Mf6eq4B962JIZXrkZ/IozhZPe72cuxQAsGt632l/m9g2gc3bN/c9bUREo8DroAwAY9kxTG1KaXGykAewWRptJ84DJocb3IqzRZz43gkAQDUIaxpuxU4AQOlbNR+oAKVvlbD/pv1tf8fY1FjrN3WBmQMiSgLvg7LvZvJ5zBaLyJWltDi9S0qL2yYmsH1zj0GgOIfrMwewA2/EdHHH0IMyAJisQaVcAQBkxjJDTk17yjkZpM2gTES+Y1DukQTkMqb+Mqy6zelMHT0HZUAGHE9NAXXGEzszMzPI5a4BAExMnAegv8EnM5bB2NQYtu7YuvRao6ptX5RzZeya3jXsZJyGJXgiimJQjsHU2Bh2bA0D1PSuwT78Z2dnAdwEACgWdyCOoDyTzyNXXqc/n2iZwSjOFlHOlftW/dwLH9MEsARPRKdjUB4RU9qNu94MXd2YLRaBj+/Xn8faKvXXlp4B/0vQw+ZrCd4XrEmgtGFQTpj8TB7lnJRgsxNZZCdX9+27pmrn6uyCzyXoYeMxaY41CZRGDMoJU5wt4uOQEuzvFqf6GpTjUq8ETdQKaxAojRiUE2iphBVTVTUREfmBQTlm0kFKqt2md+2KZ2hUCrDtmWq56muWmClq1PsZeD3NZhLN6rJKU2NjyJXLS7/XUyjkkcvlkMvlMDMTzzSiSeXanomcsakxtrvTMuVceeQz7ywp94EbItVqaFSxOAdgGgAwOzuF7b6uSjEgbHsmombSUGvCoDxkU22uSFFYN4Zi+QQK2QCT2WyfU9U/jefPrl9VOepVVUQk2mnCardJI8nPDVZfJ8Tkh76F0q07MRcEw05KTxpVU9erqkxDVRURiXaasNpp0kj6c4Ml5Tqi01auHl+NhTULmFgbmb4ylwNmvtfxIshujHG/xxf7zlVTt5szbpYrTnKOmIiWi6MJK+lV3Cwp1yHTVk4DmMbCmlei9LatKJ6cAwBcP57DjtJWYHa24ecLQYDpXbuQK5eRK5cxk88DkJzgraWd+MTknv7vRAK0yhm3yhUnPUdMRFSLJeUGlqatPNz5YOBiEKAYBMt6YHNYVH295IyTniMmIqrVVlA2xtwA4M8AZAD8lbX2w3Xe8xYANwOwAHZba7fFkcCgsrAUGGd2zgBIRg/ldntgU31xdvoAWM1NlFSdzmHQ6fh2354NLauvjTEZAJ8A8HoAlwN4mzHm8pr3XAbgDwBcZ619CYDfiSuBQeXU0s+zDzWuMh60mauuqltF3alCoYBcLodCoRBzCpMtrk4fAKu5iZKs0zkMOhnf7uOzoZ2S8jUA9llrDwCAMeYzAG4EEG0YfSeAT1hrjwKAtXYuzkRObWpv2NAgzb7sZbKOco9V1MViEaVSqe2hUWkS17hlVnMTJVu/5jDw8dnQTlA+H8DByO+HAFxb854XAIAx5tuQKu6brbVfqd2QMWY7tP7ZmCu7Sa9XWEXdf3FMvxnHdI2+VXER0WiKq6PXSgCXQbosXwDgm8aYl1lrj0XfZK2dATADACtXXm0rlZi+nUZWHEs/9jpVI5cQJPJbt5n3XjLs/cqotxOUnwSwJfL7Bfpa1CEA91lrTwF4zBizFxKk748llZRqw55+08cqLiIKdZt57zbD3s+MejtB+X4AlxljLoYE47cCqO1Z/QUAbwPwN8aYjZDq7ANxJpSaKxQKKBaLCIIFJHgWTiKirvQr895sauB+TAvcMihbaxeNMe8CcCekvfivrbUPG2M+COABa+3t+rfXGWP2AKgA+D1r7dM9pYw64jqMjY+nd6awTnRS3dVNFRfboIl64+7R6P03jPuqXim8Xgk7rtJzW23K1to7ANxR89oHIj9bAP9B/6VOQScL6XZYFA1eJ9VdnVZxsQ2aqHe19+gw76t2SuFxNXNxRq8YTGaz2Ds/33TtZPJPHNVd7VZtseRM1LnoPZqWvh0MykQ9aKdqiyVnonjVZoZHaelXBuW0mJlZvohG7k/l/2mdfG3bto5XvSLRqsSdlhw+0aDUZob72cY7aAzKaTE7K0tO6sxhO6YiM6HmdNENBmUiSohRzQwzKKfJ1BSwY8fpr09PDzolRERUB4MyUcxq27tOfO8EqkEV92y4B8DpVW1JbPciilN+Jn/a0Ke0YlDuo5l8HrlyGRMjNpuHDSzKuTLyM3kGkzpq27tM1gBB/fcmtd2LKE4uEzs2NdbRilDDUm/URVydzRiU+2i2WESpUsFUzEE5KAQIigGqgY11u+0yWYNKqYLibJHBpIHaoRxBIUB2sv51UG9mIIelaEoLd88koS24nxOKJCIoF04UUDxZxMRaf6o0CmNjKNaso5wPAswFAQLb32AZFANUShWsGDd9/R5qLlrl1qrWIChKRsrduC5j5dQrHVRKlbbWe2XgJhq8fk0okoigXDxZRGmhhKlNU/Bleo7Jchl7N25cNmHIXBCgVKlgPJMZYsr8UDtF3ihWdUeDZTu1BqeVnpsEaQDIjMt11Kw6r93AXQ+DOZF/EhGUUy1YkCFLCwt44/8BNj1aga2sGHaqWopW77igMYoBoLbKqpPSc6MgXS9AN9JO4K6nl2DeCQZ+os4wKMeo3I8FooNTwHwJWL0ar7kLOKMMmHEDzMufC4U8isUcgqC9h/ggJamNKC6dlp6joscrKAbITmTbDs6d6jaYd2JQgb/fmLGgQWJQjsl4JoNSP4JyC8XinK4ONY75+fmBf78PulngvJ+Lm3e7Rmu9bUSruGnw2Ds+GfrZG3rQGJTJSyd3nwQArL1ibcv3drPAuY+Lm9fTrzViqT1pquVJskEvr9iOTpqyohiUyUuLxxY7ev+ggldSHtLd1B7Q6XqpUaHl+l1CrfcMaFSCHsQKbt02ZTEoN1AoFFAsFhGcEX+bXnDV/0T5bgNrF2A4qinxXMes/Iw/62l3U3tAp+Pxi8ewmgFaLVzR73R1c/0wKDdQLBalrXbTOOYRb1tt8LKv4NZ/9Z/xnveMY9++WDdNQ5CdzGJ+77x3JVNWfZMvhlnTUHsf1Fv20ae1zxmUByCnk4y8YtgJISIagHbWOx5W8IuWnn1c+5xBuUetZvFyvbJni8WWQXlm5wxyh6/Rn78H4Jp4E0tENADDrjZupVEtkg99BxiUI6xdQLn8CLLZQtufic7i1Wsl9+xDs8A7btKfp8CgPHrSMNMZEdC8+SSu4Fd7P7kVppJ8TzEoR1x22SO49dateM97xoeWhqlNU0P7buq/OGY6a6dndSe9hpP+EKP0qi2RR0vg3Q5JaqZeJgBofA+5TqDZifYXJfJ/vkaiEeNKEN327I0+FJp9RzvbH4UZtyjd3P1Ue0/VDkmKQ7326Gb3kFs8qNEKcfWwpOyJXC6H83ach83TLLH4yuW8O8n1NtPLUKq4elb70IZGVCvIBwjmZDiqK+V2c7/0Y0ibu/eiNVZxri/PkrIHxsfHUSqVMPfduWEnJXYukLlhB0HBvzm621WcLXac620mO5ldWpeaiEIuIANhKde3+6W26jyudDEoU1+5C9W1ofZrgQUiGi3tNsEMU69NUfUwKFPf9ePCHSZXVUVEoyE/k8eu6V1LtXrDvL/ZpkzUgcx4xqsqtEbanfu6Xi/teus528CiGlSRGcs03FarTBd7edOwRXtkR4dPtRoVERSCpaDtttMvDMo0EO5msEH9SVYoXu3OfV3v70ExQKVcWRaAq0EVaHNl0npBvdHaygzU1Ej0OoorCNY2pwHhBCbN1n8PisHSMqr9HrHAoEwD4TpJZcYbl7QoXt320HYPpehn673W7PPtrAM97FmdyG/RjF07QbB2XHIjzYJvs7HNzT4XJwZlaku0pGuyw1naqlGVbLOJMnwqifUj5++rdjIEHI5FrXTSD6XbccnRQBwdHdLNxD5xYFBuUzkoYzS6KXXHh5JuoyrZRjeubyWxTnP+RNSZbjqTRu/Fdmp4utFJ5zEG5TZVqovDTgKhsypZH0tio9IDnWiU1E7TGadOO4cyKCdYdd1TqBzr8iKamQFmZ+XnXE7+n54O/37TTcDkZE/pi1bXxjULFiVTfiaPgx89uKy2oFKWnmP3bLhn6TX24KakadSju1sMygn2Zx/6RfzWn+7E7t1dfHh2VoLx1JT8i8rlgLm53oNy5AGcpklDupmEvludtrMPcw3b+f2yjprr1d1seFU9vjVHUPeiU2Ym/Xw269HdDQblNJuaAnbsOP31aIm5R/2sFvLV0iT0U/0Pyp20sw87qGXGMj3N2e1jcwSdrp2SY3Yyi/m980PrTBW3OHtmMygTJVyjQFevFO3mII9ilTDFqZeSYyejPEZ1NAOn2SQaoKAQLE3l1+/FOWqXeKw3lzCXbqR+6HZqXTfKo51hl6M6moEl5T4IrEWuXMZElp2bRpVbwrHTFaMG3c7eqro4zirhZgvAt/vZZpqNR3dY6h8t0f4Z9e61URzNwKDcB1ljUKpUgCA9nZvSxA1xCBB0tYzjqLaz12vfbncf25kWtNHf3IO7dipPBujkq9c/I841zX3EoEwUk1F/WLSjtmTeSUm8l2lBg2KA8evHl14bdqc26o9ohnhUMSgTxSANDwuf9ZIZoNHRqrq7n1wntV4z5okKyrnDOQSVBQCrh50U8kCrdshWbZCs3hwdbmk9p9655/n2UzSY9RpIBzkcsZbrpNbrdyem9/X46nGUFkoIKqeGnZSeuZmMqDe1vYtr1ett7LDX8WgJikHTnuY83/5aCmYDLtn6KlEl5VGQGc+gssg1hePS6/KENDqi10KrMdosNSdbpVzBimzrMmVcY5ltYFHOlbHn7XtiK9U3kpiS8iBVKqPVK5ZC+Zk8dk3vWhorvGt610hNPNCu2uMwaseg2RhtlppHQAUDHctssgaVUgVPf/HpvpfqWVKuw9r4qpcDK6XiXLmMQpDFuibvLR9kZqDfaofeJLWXrmtDbbQgeyvR4+CCVNKOQSv1alFc+6XLkLHEPPrqNWHF2Y4dN5aU+yyoVuGm3Z9rMW658myf2poLeVlkIpeT1aFSzj2su5lxyBeuDdWlP1oKaLc2oNtZl5KsdgpIlpjTyed2bAblARjLZDA1NsQHX3Eu/Nkt10iJ1yioRkvBaQw+0dIwMyOUNKy+Tova5RlTrt6UkKNUlRmtuk1bp7ZoBmQUq+Upfq5DmA2G3wmXQZkaKgQBilrlXgiyOHPI6QGWj0ftZSzqqLQtU30sBVMn3PjmzHhna3z3A6uvqaFipA28VXt4I25VpLhyoNHxqL2ORR2FtuVmVbWj3sO6W9FOcmnugT8q3DOm36uudcs9+9qdB54lZWoqbAvvrmf4shzofDxpajQ2OW3VtEDzqtp2elhHF6TvtAd3UgVFqapkLUmyubHD2Yms11PcVoMq0EEBnCVlooRrNnNZq05NtUE9LcamxpYtKckSc/K4scO9LIE6qBJ2ZizTdm0cS8pEKZfUqvtesV+B32xgpZTZpZO7Tzad0rjXRWRO7j4JAFh7xdquPt8IgzLRCPF5UgQfpbmXuu9M1vTU5LV4bDG+xHSw/V77bzAoU6qN2tCoXleqSVsbc7Pe/Em+Dqh7vcyX7UrfvTQFMShTYtUuOtDNEClWYS5XnC1KteF8FQc/enDkj4Pr9BVdA7ecK6NSqtTtzc9APfrimi+7WwzKlFi1AbW2bbTdAJvEKsxoz9NW1dSdln5dtWFaqr/ddRPtkV1P2jNsvmt35ah2DLOfBYMyLTczI3NkLywMOyVtabZ0Y1ICbDeWep4iaBk8OcNVe2qvpVbLPzosPXui0ns7tA8YlLtU0PG7bhWokTE7C5RKwOrVw04J9SDIBwjmpBouO5FNbQ/rXrSqiQFYeo6Tq/lJOwblLhU1KGeNSXrGjEaQC8iAVMumpSo6bo1qYvIzeRz86EFUyhWUvlXqaqpXCvU6PKlWu7Nn1XLDsILC8O4ZBmXPFE4UUDxZxMTaCUzGuuECUCxKCXh6WqqoAfl52zZg+/Y4v408UNt5Lcr1MOXQqe4UZ4uY3y/Z8cxYOF0TS87D5wJ8N1z19zAzsgzKnimeLKK0UMLUpphXdSoWpVp6fFwCsls1ygVnT4NytIOSL9xctj6sKNOt7GQW83vn22qTpvrcLE3RkvQo92MYBf2a8CNODMppNDUF7NghP09P9+1rms2mE+0R3M7YYJ+mgHSzDJmsGXJKRHRcZS9tcq4dOsmZDRpNcV2b3U4oMsiMeOqDchg4hr9k18hpUoPkgqxbLAGoX+XXaQelaI/ZUZkQpJXouMpe5gEO5uJfPGRp24UAlXLFm5V8fFo/dxRFV+ICep+Ipp/XZjuWZcT7/P2pX5DiksoJXFI5MexkpFKrxRK6EZ2dyy3U0OmSjv0W93KWQPNFKXyQncwCld4yDXFyq5f5UtsxatwSq+6a9On+i5tbAjSuDGfqS8qDEFiLXLmMhaoFwIdAv9W28933wvtazkfbzXSb0Qk8OtGP5SxpudqpEke1lsRn7j5MSzt7XBlOBuUBCKpVzANYvYJV5HGJLrzQiuvU1Cy33s10m9EJPJLG5epHtfq2dqrEfgXl6NzZAOfPTqPoM8PdT72MuWZQTquZGZkoxPW+npnxtge249olXVVRLwsv1NPv6TbdGMi4pgJsR7QtN9rL2gWtUau+dZk1G1isu2Zd37/PzZ1db4IRDo/yQ3T4X79Vg+pS96RuS86pb1NOLReQ3dCo2dnhpqcNcXVoGhaTNUAFPa0R2ynf2nL7za2SNcjMhsvM1f7zuY0/TVxz0aCG/rmhct1iUE4zNzRqKuYx0X3U6wUf5Tpc9atHcK/rqtYTLQkS0XB0O2NYOxiUaahckBnGUJl+5qDjWFe1nmGUBLs1zHNL1C+9zBjWDrYpp51bFcr9PGAuyMTZNtyL2l7Yo9Jz1/UUB+rPjnZy98m6k730cjx8O7f9Fsf63tQ//Szdxokl5bSLtiUnoF25G66k1k6prbYX9jDGV7pqdfcvjirwaDt2vX1aPLZYd7IXH45HUkQzL8DpY8d9Gy+fJv0s3cY97wBLypSoNuVGmo0zDorBsh6Rraqrhz2+MihKj+nxHx9fepDHUbrqti1+2McjSdK6vveoidYsteq1vWzegRhaalhSToNgASgno+qmW9ESXb1ZvOLsIDYImbEMe/ASDUm0ZmnQIxdYUk6D7Gpgvn8dEzpRb8xsM52MMez3OGMiGj2VcqXu3AHNlj7tJ5aUaWAy45mOx8x22kM6P5Nfmgi/H0OSBsn1Xm62H81W4iKi5twzaZBzB7TCoEwjZZQ6JkXT3nA/GsRk1ybG4UhEycKgTF6Is1Tbj9WnhqXb1Z+W5uVOyUxevovWeuya3pXoGhzqr7aCsjHmBmPMo8aYfcaY9w+jF40AACAASURBVDZ53y8YY6wx5ur4kkh9USjI+ORCYdgp6dtEGxQfBpXeuGs7O5FF6e4S9t+0n8cwRaK9uVtpGZSNMRkAnwDwegCXA3ibMebyOu9bB+DdAO7rKLU0HJOTQKkEFEczEHI2qXi5oOLj+tRJMTY1huxkFivWrEClVMHe39ybigyOG8c76vvZTDWotj1Oup2S8jUA9llrD1hrAwCfAXBjnfd9CMCfAHiu3YQS9cvSbFIDmoQ+DUapWWCYTNYsjZtPQwYnO5llTVgH2gnK5wM4GPn9kL62xBhzJYAt1tovNduQMWa7MeYBY8wD1Son1KfedZsDbzStJI2Ofi840ovMWAbjPz7ODA6dpudxysaYFQA+BuAdrd5rrZ0BMAMAK1debSt8JlIPemmLXjy22IcUxc+N0wbqz1lNjS0Np/Ng7u2Tu08CANZesXbptei5jc5AR+nWTlB+EsCWyO8X6GvOOgAvBbDDGAMAmwDcbox5k7X2gbgSSpRG0d7TvVb/DXKxd1quXibQnVvXTg+AQZnaqr6+H8BlxpiLjTFZAG8FcLv7o7W2ZK3daK29yFp7EYDvAmBApkTzqbdxt8OiasWxVGU7E5r4xue0sp1+OTYptRGUrbWLAN4F4E4AjwD4rLX2YWPMB40xb+p3Aqlz5aCMcjDac13HKcif3vbI3sb1tTWhiUd8H24XFIJlM9ANOwM4dIzJ7bUpW2vvAHBHzWsfaPDe6d6TRb2onKeVFMemh5oOH7gl1UzWNHxPMFe/7ZGrI9XHUl18gqI0KdTOs8xq7PTighQ0UqLtpktrqGaATDbT+rP5AMGcfJZDqSRA8Fj0HxdSoSgGZRop2cks5vfOI+hiYdNGJeY0chkadxxdRocBmqi/GJSppd0nT6LcZPya6/wDhEM7RkV037rptexK362q0H3m09AiolHHoEwtHVtcxGKTuV5qO0WNkmgHoW4Wd3Cl78x46+rztMjP5E/reBXN1EWlZeyuOya1mds07Dstx1WiqGtuQgRguEM76vWejlNcQ5K65fPMVN2IBh+n3jFOU6/36JKj7PGfbiwpU9d8mRVr1NuCR7H6ONq5qZG0dXoamxrDxLaJpWDshkixxJwuDMpEfcTJEKgT0RIzwCFSacSgTInndeDzOGm+qW1rrtfOnIZSI4dIpRuDMgGFQriu8oTfPafdYuHLejPHHPg4Xnk4akuJ9dqYgdEvNbpZvoD0ZkzSjEGZwoBc+7OHqkEVmEdfezMntY06KAR+1xq0oVlbczelRhvIsAGf57+uFZ3lK60ZkzRjUE6bIAByudNLxFNT8n8uN/Ak1T44+cDpjps4hULVoApo/s233sw2sKgGVQSF4LQamUaZE1Znjz4OiUqbbBYolcIScXn444p9fnBS8mXGMk2HtA1ryFk1qAIVsImElmFQTrsmM3UNUqsHJ1G/xLGkJVFcWH1NlHCu81s304ASkV8YlIkSzmTNssUjiCi5WH1NROS5/Eweu6Z3oZwrJ6oneRpFh212g0GZiDriOkYxOAxOcbaIE987AQColCrYf9N+HntPVYMqKqVK16vCMShTQ4G1yHnQO5vi1UsuHli+WhZ7yw+Oe8iPXz8OgMd+VDEoU0NZY1DypHc2xSMznkGlVJHhOD2ITmxRW7W6a3pXw1Kcq9pjKbs7btEKIFkTotRyw89GafWzuLCjF1ELQUFmWFo2tSctabSIQj3VoIpMNrP0OU4U0z6XoYkGsaQew6AYLM1N0M065aOMQZlGRr+GBblxrP2a2jMaxJqVfHwe9tTJIgocj96daC/7JB7D2hn73NwEzTJxacTqaxoJrlo2zlz3oKrWKuXK0kO2UTthP/aPkitpne3c9ct28NYYlCk2+Zn80oMi6e1EgwyCmbEMtu7YmsjSDw1H0jrbJXFhkGFhUKbYFGeLsIFFpVTBwsGFjj6b9A5A7qHTS69mok64powk3DOc3759DMoUK9cRqtMOUdHewEm8aV36R7kj2CjUgIySpJWWOb99e9jRi7zBG9ZfS9X5nMrTK7xnRg+DMoUKBS+WciTqFBfloFHB6msKTU56s5Rj2vU6f+4wDWN94qXhQk065gX5oO1jGu202GwyFKK4MSgTNVEpV1ApDz6jUjt/bpI6kmUns14O3wrmgrbnJHZttG4cbRLabGk0sPqaqBlPKg6WdSSbH3JiBig/k18KiG6SCTc5ycS2ib7OZuUmRGk1GQpRnFhSppE0jNItxc9N4Qksn297EKXXoBC0Pad3tLqbVd3UC5aUaTQxJveVmw88O5FFdrK/nauiU3g6gyi9BkXZx9o5veuVzqMZhKTOR01+YFAmIAiAXA6YmBh2Sigh3Hzg2anR7u3s+5zeLnMEYNm80pRcrL4mIJsFSiWgyM4s1H+9ztvsqoo5kUnyJhCh1hiUiWigeg0kxdmilNL7XG2eFEmabpNaY1D2VO5wDkGls/mjR1lSx+y2K22lv2inLeoNS8ujhUHZQ+Orx1FaKCGonBp2UrzgpniMzo89atJY+nPV2O30aOYEHs0xkzM6GJSJaCjcJCONSnecwKN9Pq2vHBQCVMqV09LiZqlLS21Qt5IflHM5IGA1L1E/DHu8t+v9zFJgc0ExWFo2df9N+4camBtVp7tZ6tJUG9SNRAfl68dz2FHaCgSs5iXqC473jlU/+0aYrEFmXBYtHnatApdp7F6igzIR9cZVKVJj+Zn8spm9eimFtjv3djeiwZ7t8MnFoEw9SVuv4VHjVlfqRasOW0nnpvp0Jb9hl0IbqQZVVMqVpU5fbIdPJgZl6kkaew3Tcq06bLXD97mjk9K2nRnLYOuOrYlIK9XHoExEQ1c7d3Q7fA3gRL3g3NdE1DYbWFSDKlZk48/P1yvZRVdqArAUhN3YdZ+qZzkPdWdcfwYb2L61sycRS8pE1LZqUAUqGNhDNCgGiWjPBfyZWSvOjmn95IZIMSAvx6BMRF5LSnsu4MfMWknpmEb1MShTLDrpge2qrWo/43rxsic3UW/GpsYwsU2WYvW5tEynY1Cmnrn2vWj1XTPROayjn1lao5c9uYl61k3nORo+dvSioXBVa5y4YrlRXw0r6VxnLrdcou+GVZXuasOyE1lmsjvEkjKRJ9KwGlYzSWi+cGOy260VGrZhLVTBea67x6BMRF7wsfnCp9WXuuFLj3BqH6uvyTus0qZabsYvAAOtOvYxqEXHircznGjYvcGpMwzK5BVXhUsE1F+Dd9BVx74FNZM1wLxUEWeymWEnh2LG6msi8pZbMCMoBl6MAU6iaBU8V47yH4My9YS9hYn85moWuHJUMjAoU084TV66+d5bmkSSZkUbBJ8LE2xTpthUyhVkxtjGlRZLk8ZgNINypcy+DaOqUqogM+7ns4pBmeLDZ1gquRJHt8OGgnyAYC7wr+RS53p2tQJc2Yj6hdXXRNSTalAFtNDRTXtlMBd03AwSzQgMsvrctc8yIFO/MCgTUc8yY5m+t1e6XsRunK7LCPR7iFQSZhprJboudVInQkkLBmUiSgQ345crpQ4iIxD9Xp9mGutUktalTjsGZaIhcxNkeNemOqIq5UoqjzWXc0wGBmWiIXOT96e1nbLfvZyjmZ7MeAaoILWLfnA5R/8xKBPRcPW5136SMj1BIUClXOlrSZYzo/mNQZmIqA2DaGLwcQGMXiS9g9wwcJwyEVELSwul9GG+CTdO261+NagObP3WaHIZ9p9ojiVlIqIhcuO0B9W7e5hrRLtAndY2/XYwKNNQRMecEtHgjFoV+ahh9TUNhRv7mRnPYESnTk4V18MZwFI1bBoFhWAp6OVn8ti8ffOQU1TfKFSP95s7l4MuOLCkTEQ9i1ZH9mOGraQsDpGUUugwq7CTonaymkFhUCaiWMQ11KZuoEhGTAaQjCFHSck8pBGrr4nIG64jULeBIiklah/4nnFIK5aUiWh0MCZTwjEoExEReYLV10TUV+VcOdU9sok6wZIyEfXN0qxOfV7zmDrHHth+YlAmUuwkRGnCHth+YvU1ESJzG5MX2qnyji7JmIQVoHzEHtj+YUmZiLzSbpV3kpZk9BWrsP3DoExElFKswvZPW0HZGHODMeZRY8w+Y8x76/z9Pxhj9hhjHjTG3GWMuTD+pBIRjT5XLT+o0msSZiBLk5ZB2RiTAfAJAK8HcDmAtxljLq952y4AV1trXw7g8wD+a9wJJSJKgl6rhKPziLP0mj7tlJSvAbDPWnvAWhsA+AyAG6NvsNZ+w1r7rP76XQAXxJtMIqJk6LZKOLqUKUuv6dVO7+vzARyM/H4IwLVN3v9rAL5c7w/GmO0AtsvPV7aZRCKiZOk0oC71/s/0Jz2uSpyTuPgv1iFRxpi3A7gawPX1/m6tnQEwAwArV15tKxyBQkTUdyZrpEc7Fy/3XjtB+UkAWyK/X6CvLWOMeS2A9wG43lq7EE/yiIiI0qOdNuX7AVxmjLnYGJMF8FYAt0ffYIzZCuCTAN5krZ2LP5lERESjr2VQttYuAngXgDsBPALgs9bah40xHzTGvEnf9hEAYwA+Z4zJGWNub7A5IiLyCCcQ8UtbbcrW2jsA3FHz2gciP7825nQREdEA1PYW37x98xBTQ5zRi4hi4UpcblgPdW/QpVfXK5ul5eFjUCaiWATFgHNRq16D6qCnv+R0m/5gUCYiilkcQW7QE4iwtOwHBmUioj5I2qxcLC37gespExERAK6v7AOWlImIiDzBkjIRUR8EhQBBMeB809QRlpSJqC8q5Uqqh0e53ujZycEF5ehazEGh83muOZHI8DEoE1F/VMDhUQMWXYs52nGrXezsNXwMykREGJ2hQL32+k5ar/FRw6BMRKnn1jOOo3RoAxtL1X20KrmbqmhKJgZlSgXX1pbmNk5fuHMxqoGmGlRjqbqPViV3UxXd03drhmAUag6ShkGZUqEaVDkFpCdM1qBSqgw80CTRIKqSK+UKKuXKsteyk9nYag6oMxwSRUSUZpXWb6HBYUmZiIjIEwzKRES0jGvvZ7vy4DEoExHRMkExADLyM9uVB4tBmYi84nrIs6f8cGXGMhyvPAQMykTkFTcrFXvKD1d0yk5WYQ8OgzIRpZ4rlTMAhaJTdrIKe3A4JIqIUq8aVNmGWgerrwePJWUiIrANlfzAoExEROQJVl8TEfVROVfu6nNL87WfsqguVJEZy8ScstaCQrA0HWp+Jo/N2zcPPA1pw5IyEVGfuNWnuuHmCK8+Wx3YVJi1C7dwfeXBY1AmIiIAkYxApOd1diILgD3TB4VBmYiIGmJpebDYpkxERE2xV/rgsKRMRETkCZaUiYi65DpGuXZXol6xpExEhDDAumUL2+E6RkXbXX1mA4tKucLFPjzGoExEBJlqs1KqIDs5uqXeaiDDq7jYh78YlClVXGmIiMhHDMqUKq66kYjIRwzKREREnmBQJiIi8gSDMhERNRUUApRzZZRzZeya3sXpNvuIQZmIiJpyQ77GpsZQzpU53WYfMSgTEXXJjfdNw7jfsakxbN2xlVNu9hmDMhFRl9xqSmkY9xsUAuya3sVq7D7jNJtERCMkKAQIigFsYGPNLARF2a4rKbvx/pu3b47tO4hBmYhopATFAJVSBZnxTOzbdlXYALBrelfs2ydWXxMREXmDQZmIiMgTDMpERBGdrhTlK64GlUwMykREKjOeSdRSjI1kxjNAJewdHqf8TJ69sPuIHb2IiKhtxdkiyrkye2H3CYMyERF1hL2w+4fV10RERJ5gUCYioo7Vti2zXTkeDMpERNSx2rZlLlIRDwZlIiLqChepiB+DMhERkScYlImIiDzBoExEROQJBmUiIiJPMCgTERF5gkGZiIjIEwzKREREnuDc10RE1LX8TH5pUQrO6tU7lpSJiKhr0Zm8OKtX71hSJiKinnBGr/iwpExEROQJBmUiIiJPMCgTERF5gkGZiIhasoFFOVdGUAiGnZSRxo5eRETUkskaVEoVBJCgHBTl/+xEFtnJ7DCTNlJYUiYioo64gFz7M/WOQZmIiDo2NjXGoVB9wKBMRETkCQZlIiJqyQZ22f/UHwzKRETUUjWoApAOX9Q/7H1NREQ9CQoBgmKA7AR7YfeKJWUiIurKyd0nUSlXEBQDVEoVDo2KAUvKRETUETeRSKVUWfZ6OVdmm3OPGJSJiKgjbiKRqMx4Rl7LDClRI4LV10RERJ5gUCYiIvIEgzIREZEnGJSJiIg80VZQNsbcYIx51Bizzxjz3jp/X22M+Sf9+33GmIviTigREdGoaxmUjTEZAJ8A8HoAlwN4mzHm8pq3/RqAo9baSwF8HMCfxJ1QIiKiUddOSfkaAPustQestQGAzwC4seY9NwL4O/358wBeY4zhXGxEREQdMNY2H+htjHkzgBustb+uv/8SgGutte+KvOf7+p5D+vt+fc9TNdvaDmC7/vpCAI/GtSNERESeu9Bae26zNwx08hBr7QyAmUF+JxERUVK0U339JIAtkd8v0NfqvscYsxLAOICn40ggERFRWrQTlO8HcJkx5mJjTBbAWwHcXvOe2wH8sv78ZgBft63qxYmIiGiZltXX1tpFY8y7ANwJmdX0r621DxtjPgjgAWvt7QD+B4C/N8bsA/AMJHATERFRB1p29CIiIqLB4IxeREREnmBQJiIi8gSDMhERkScGOk45yhgzZq0tD3PbcaWhm+30c/+7UZuedtLX73MIAM22rxPZGABrAcwD+HzcaR7mNdLu54wxP2OtvbPVe93f65zrt+iP1wH4trX2s72kuRVjzNsgQyv/O4BXWmv/Oabt1u7XzwP4dQBHIZ1SP97GZ2LfZ2PMOwCcp78WrbV/1+TtTdMXY5r+CMB6AOsAXATgU9baf2z3+yPX0q9DZn0EgNustV/qIU1DuTe73X7fzs0gO3oZY84D8D4AJQD/FjKW+V8AXAxgAUAVwFWQm+g5AGP6PwB8E8ALAByGTOn5/wB4HoBNAJ4AUARwNoBrAawB8AMAzwJ4KYAHIWOnjwP4fQDvBPAUZM7uByE1BpMADgG4EsAJAI8AOAvADgCPA/h73e5/1LT+PIA5AN8A8GOa1i0A/g+AjQC+BuBlABY1TQGALwL4GQCvgNwQhwHs1W1cA2AngP8E4LcBnAMJOOt1/zfr+68FkAfwIkhv+Ed0Hz4K4MWa5h9AHnj3AKgAWAXgUwBOAThD9+EN+rPV771O33tK07oWMqXq9wGst9b+iTHm6wAeAnAugA2atjMBPKCfz0B64v+YpvEEgJ/U/VjUdF0AYLd+/zHIGPev6TYWANwC4Pl6LC3k/BchvfrXQ2aCexDA70GG4Z2j27gecn18H8Bj+j4A+AokYF8I4HWQa2yHfucNAF6t56wCeTg9CXmAXqfHbZW+fz9k6J/R/TCQ6/Ag5Dr8qqZ7AcA7INfSt/S8fhPAr0KGF74TwBEAE5BrMKOfWwngp3V/12jaztHv/45+JgBwG2TY4Ysg53sfgClNx/m6z9/V9D8L4IcALgVQ1mNwUI/B32kaFgFkIdfzPsg180ZNz8V6Hu+EjKio6vsfhNwvz0Hu1Rxket136zn7O8h1/D59rQy5bh8H8BbI9fakpnsOck2dAnCJ/q0KuUcegFwjWX3fjG7rFyH3/L0A3q/H9indvzLkPjkI4Kf0uH0NwDbItfFxPX9/rPtwUo/vpZDr67WQ6+2fAfw3TeflkWP5fsgzYgFynbiCjTtX6zS9AeQ5t1L/zUHuKQu5njbpObwXcv//G03DGj2uVvfhxfr/0wC+B7mnnge5PvYAuANy7R+HPFvmIc/UyyDPjAf084t6ft8Mueay+ve7AbxGj/NRyL3zIsh9+nU9HhU9d+74PgXgJfrzswD+EsB7ABQg99urIffKTwG4C8D/BnAAwFYA/163d7Yes3X686Iev88B+M8APq37+0PIc/SfIPf7UT1vr9M0fw3A32rafkaP2xzk2n0G8lysQJ6LPw+5x68F8CPdv0v0/RcC+Bjkefa7+vdLdT8zeu4WIc8BA+DLkHv72wD+Qo/nCwD8kp6LT+t5eVA/dxzAagB7rbUn0MSgg/L/hlwYr4U8iCzkgC1CguZzkJskq6+vgFzY7mJY1NfWQm6KlZAb+GnIQV+pn5+HHLinIA+qdfr6GZAbdh3kBjSQoHiRbmelfsca/b41+p4K5AJYBzlB7kF2DBKMV+n2je7DefoZA3nYnNT9c0HvEORCPAS5KDbo557TbT0LCewnEF4MAeSkHodcOGv09Yr+bRFykZd1Gxk97Ccixzarr83r3w/oPp2rf7sLcjOv1vQd13RY/Y6Vmv6yHu8VejzP1fNRgtzo7pxWIA/F5+mxOlPTdkS/dzHy2n2Qh90p/f2UprUACdCHIQ+yU5AH0Ws0fRfod5c1HUf1e8/Q/Vih2zuur7trxp2rjL5nXeT7Nun/K/XnI7rtp/UYP6XpXq+vbdBjtAh5eF2q239WXz9Tf3cZuBX62nFN5w8hAXKDfmde/74Ccu24gHpCX8/qPj8DCRoXaVqKep7W63as/p+BPGzPgQS91+rn1+t7juj5zOvnn69ph6bvOOQ6hB7fjQjvFZdZegZyr1X0fe76cwGmrJ8L9D028h0r9HsM5Noc19eMHs8JyH0WIDx3KyHX2z6EGVT3PFmF5fdeVtN9EvKQtXqsNujxnNdtj+l79kLuy3H9vucg19Jqfd88wuv2qB6bOcg9U9b33avH9OW6j09Dngv3QYJCRf+t0eNwhr7vSd3Ogv7NPYsymh73/av1Oy/Q81BFmEFcQHjdPQK5dtbpe9x9cEq3vQ4SaM5CeE8/oZ9119oeyH09AeBPIVMlr4IEtZchfC4+AsnAnAEJaudoWlbrdub1ezdAruVTCAPynB6fA/pdY7rfVre/UrflJqXaFDm3VX1f9Bpa1PdnIdene57MQe7lE7q/gNzr2cgxquq2Srq9dZB7fn3kb/Oaxh9pWqBprOjxcxltd83OQTLnR6y1v48mMjfffHOzv8fqlltu+a619nO33HLL+ZAc1W7ITXIUcpCqkPmwvw/gDyAncxYSKFywWgU54IchJ+hcALsgJ6gEKSXkIbmpP0VYGnGBq4Lw5I5BDtz9kJt0FeSCOam/Pwi5QS+GnNSjCIPjAf38ak33dyAn8Bnd3cd1GyXIxVuB3BCf0v2pQB7G74Tkhr8BOcHuQbkWcjO4wLKoab1Nv+dT+vkqJFf2U/qeWyE5uO9DbiqL8EI1um85yIPm3ZAL888gF94BSAl0LeRBl4fcNPv1eC/qvn0BUrqZBPDnmv5DkNzx2QgDyQP6uUchuV43jn2zpme/bmejnsdVkBLnOXosVur5c0EP+rnLNP1fhTwUfoTwprwK8vDaCQk2LhO1Xz/jakxczcyzenz36n4fhdygf67ftV6Pc6DpcNfYOsg14Wpyirp/X4E8kDfrOfga5KHogqJ7UFUhJdAX6XbdNlxQeQByfWX0WD6GMIAcQHidzENy9mv0POzXbZ/U434ewutkJSQ3fypybJ/Vf2U9XvfrcXWl/y9AroN7IQFlC4CH9XhkIdf5GMKM64K+32W8ntLjeq6eh7yegyOQledW6v66mikbSaurNTih52dCz9URhMHjXD0GTyOsqXkB5AH4A0htxW7ItZrVfXXpzCK8l9doOlfpd6/R81HW/Vihn/0+JPOxWtN3TN97HOGzYY2e5xX6/n16ju8B8BFITcSX9LUx3e5Zuj+7IcH1j3F6IcLqOV2t+2z18zsgQaCkx2qDfuaInpczdT9W6meh+3YEEpAW9D2uUHMW5H6eh9yba/T/POReuzxyvKwe8zP1u5/S7zyBMPhtQBic8/qe45Aau9sgpfWnINfvjyAZjSOahmORfbgHcj1X9e/HIM+To5r2JyLn9iE9jpMIg+hxhAWUA/pdgATMOUjmdhHy7Niof3MZuznINbZK03VCz8kxyLk+DHm+Gn3vBk33Ud3vMyAZrqdvvvnmu9DEoEvKn4Hs1EFIFcdFkIPqqqeqkB0dgxy8E5ADMgapPnYPyOsQVg0vQk7GSyAPlA0ISxI/gFRzP6Lv+yLkovgd3eZjkFLNcchDdA/kInuJfvc5kKqgRQDvhVSBueqcu/QzX4FUw70acnNtQlh99Zxu+7i+/hikCmcRchEWIQ/sMchNegDAPwD4MOShOAl5yK7T47QI4FWQB82LIQ/HCyHVQ18F8Ne6zTykKufTAP4VpKplC+Th8Cbdr7MRlpwzCAOUy939UNP1kG7v/4dUE52L8Cb5PqTkdxuAKzS9/6yfvRaS0VgPCf4V3f53IZmGcf35xbrdr0GaGsYhVXkX6v49BbmgX6DbmIOUYv5R0/9fIFWdlwJ4JeRh8UJN+6cgQePfQDIFV+t5eAjyMLoU8oB9EvJQq+r3H4JcA2sgwfFJPSafipzn2/U8fVW3+yUAf6iffRXCkv+Zek6f1GP0bgA36bm4T/fzMOR6+yikOuyLer5+GfKQfiHkAbgX8uB9HuR6vBnyoPh1PW9fhzQXuEzKuZDAsALycDgDElxfo8fgx/S7Pw3g7QA+q5/7bU2vK3GsQ1hjAj0WGUiV4oshmeKdekzXQ+6Zv9Dz+h7I9Q9Nw+chVa8/Brlv90MC/yv1eD0Kua826nbXQR5q/6JpzujxPQLJsN8KqRJ3wfBZTde0HuPnINfZej2n5yN8eF8I4Lf03zqEJcyrdVtPQZ47RyDPqocA/Aok4/co5PnyryE1N9sRZoDP1PSejzAgHoJcezfpdv4IUop+BvKMczUKL9LPuFqzM/V7/hnS5FeBnPPNCKvCnw859z8LuaZPQM55UT//LOQaOqn7fhxyX7xYz9sNkGtjjR6nHZCq18N6Hp7W113HYFf9uglhBnGt7vMl+p33QoL3IU3Tefr9OyAFmFcjzJzdCbmvfx5yT38T0kRxF+QZ85Oazuc0Pb+h58qV+h+F3CfufjkOecZthFxfT+l7XgUpvJwHyXjcqcfINZV95Of7tQAADs1JREFUXvfpck3jyyDPqoch18GHIffZId3vcd32MYT3ypWQ+2NS/36bvv8U5Pn/oLW2iiYGHZTfDTmYb4XcIC6nux4SUF8ICRYuN52D3ISbIQfZVZmthlzMCwir5DZCTr6rngTkQKyFPKQ2Qi68Z/Q7XLWUq2K9FPJQOQa5sI5CTtyzkAe81e25HPom3Y4rraxCWNJZod/3mKbjPN3GCt2fJyEnzVWzHUJYDZ6JfNcC5GLdrPvh2nPcXORrIBfASf37Cv3sYYQZjvHIPqxG+OA4rsftLP3OTfr3+3V/LkF48T6r6XHVivsQPnBOIqyaO6BpczloCwkia/V9Z+vxfxph1Vy0Ogv68wrINbAIKR25EovL6a/Xc+lKjmUsr9J/WtNXRJjjParp/zokAGQ1TXk9Tqt0W7fr+Xmjvj+r6V+DsMbiiP7uSrKLuh8VPT9jkIzTMQA/ods5V4/7ft3HVZHz4TqLrNX3QPf1Wf2eVXp8vgnJ7JwNuXaPQB5+V2ra90CChPtcQff/cT1n1+h3rNd9P65pCSDXWk6PxSZIO+G/1/27CxJEXZPMocg5cOesDLkGXMZuNcLrxgUcV/3uapee0fecD7mPF/V31zzlasZcjdYiwgykgVxbGX3vMUjw+hyklmM15D5YgNyvmzUtz9fzWNbz4u4398xwmegjkMy2RVhNmkHYtOKaXtwzZb+m+7cgbd+ub8PbNV0v1b8X9Py5Z8oRyH3zQj0+C5CS4gl9j6uq/5H+ewxhtekNkAya0XM4p/s1psf0dkgGYgvknngOYRXyY3p8z0H4/LkB8kxYp/s5gfA59ENIJnIBUjOzQY/TCUgG7N8hbAPfjOXNHycQZmzc9e2uCddsZPT3o5DzugVhU8AZCEvaBnJez9Z9ciVx16Tm2o8DhNffmsj2XXWyu3fP0HP1AoTPXXdfuqahBf3cmQibTTOQe91Vjy9CrhP3TDoFuR+OI2w+2QfgXmvtB9DEoHtfjyO8AX8ICUyuWmYrJAd5NsKi/3UI26xctTX0/3+BPGROQR602xDeXHOQg+2qt11b5HOQg+tKAWfov9sgF9SbAfwcJNf7bsgDuwi5OLMIq4NeCjmReyAPxEXd3jmQk3VS030uwpu+BDlJeyE59FMI2+pc+7HLcEDTea5+9hnIg+SEfiYDOckWYa3ACv3sQcgNvhNygWQgJdZX6u/uAbUZ0rHI5X6r+vO5kOvCtXn+FSSnWoUE6o2RY3GWHsvzEVbPuraf70Ae5O4hu1aP3Y8gzQGu3Tej23LX4nP6GZcZc9Xwrtq0Crk+1kA6aL1aP/ug7v9myI1+HHLdPK6fzevnXqznYgJyTV2i6TgGCUoTkOuurO97of5/Ure7U7d7KSQDc4X+fUHTbjTNl+trfwLgv0I6P70Rcn0/rp/fA7kmXDXahQj7DJwNOd8uk3UW5HpfDbm5N+t736jH5ClIQH5Yj83Zui+uFugYwozPSgCfhNQYHdH93qrH4RP6+jV6TPfofnwZEiwugAS+j+hn74DcMy6T7Kp/Xd+DH+ox3qfHchWkunMaUitwHFKb8zCklAiEVd/uXvgWpJTjOnM+AODHIc+LzZBqzbfo8XsJwrbhQwivydWQ69RlRCqQa/RySE2B1W0W9LXzIef1LMh5dcFjN+T6Pw9yb1wBCd4nIB0j11hrL9f3whhzDHJf3wl5Dl2ix/YZ3f89kELKPsh1cwPkHt0BKbH/E6SJawvk/lwJKcHt1uN0MaSGzPVxeVL3eQ+k5Pwo5H66QI/LFyC1aE9DSv3vs9beoj3Vb9M0LkCu8xfqeXwdgPuttZ/SffoGpFR7NoDPWWu/Y4w5F3K9AZJBKSNsGz8b8lw5EDn+5+rPruDkmjGO6d8O67ZWQ66hMyCZm5/U9y3ouVmDMGPsnj2uL8UPINdDVl/7F8j9Pw6pMn+Jpm0T5Png+mLcA3mur4U8T13fGndv/gCSGRpHmOH6G0gMOgWpmb0SEodcn5ajkExN01IyMPiS8u9B2pF+BLnJLoBUUZ0BCdDHITfWUYTV3E9B2hBdx6xx3Vxe/38IctE9D3LA5yEH+hTkAt2s21ujny1CDswu3e4+fe8FkBvla5Cb85WQC+R5kJNaQJhLfwhyQjdBbthndPsW8pB6C+RmeBHkonoYYWAvIex45npXn9T3uCrTFQh767lOHVWEVSSrIRfufshN6UoOruRagVSnvhRh29mFuv+u3eOI/ose11MIc+vuIWkR9tJdC3nQbdFtus5p+zXtD0OCjeuk9CzCjj+uurik3+syUBfo3wP9/wz9LldCcDnjR/U4X47whjCQa8R1kNqgrz2h23I9s8/XY+dKx+5mPoywJsFCHlQ3QYLShZCH/6t1X1yHFNeGVYLkol0HFddutQB5YD4LuR6fhFSZlRBWKbqe664z1En93w3vcp2pXM1MFWGHPldjkdc0fhlyrZ4BuXan9VhNIGxX3AIpNV2J8B5zJeWzIdf5FshDYxXk+nelWtfR60xN47N63C7Q/V2BsG1tjZ6HA3rM5yDXRknPxUcgTQnP0/du0bReq9t2meb1+pkn9PdxhBmNA5DS7g8gQeYshJ39ViLMpLo+Bjv0vedomiZ0X/ZC7qHHIYFyPcJM9z79jksgwWUXpMR5StO1CpKR2QLJmL0eksnZCCBrrX0v3Ak15gOQXr0f0+PqekL/BMJRHicR9v4e03+uU6vroLYWEsQndP8+A7nf3gJpknLtw0a/6zc03ddpUt4DaebII+zlvxvAYWvtX6KGMeYf9fMGwKustTfWvqfm/b+nx/hTkELSXsi5ORtyT10Eafp7CeR8HtM0uPtqUX93HarWIby+ViDsd7ESYW3DSYSd71wNw0LkeLn+Q0f1mM7rd7gRJucgDOrumnbV9McQdkrcBHk2lBC2sWch97IrDLntnYh8dhJhAcAVzH7HWvsPzY7loEvK34ZUt10LSeh9kBvgIoTd818HuZGugxzYXZCb2PWgLei2noEcMHfh3QK5IB6BHPxLIW1J5yGsInc9Gx+AXDAP6/8fg1w434cMJXop5CQ9o+95vab7VZCL1D3cFyA32VpIqeLlCId1XIawB+89kAsSkAtkHYD/D5Ir3oSwynYOcqOfguQOXbvyUU3bFQhLjyVN0yTkweGqQF0V4kX6N1e74Dr2uB7Xj+g+XoUw43MvJNddhnTCsJD26rcjfMgdhlyMV0BKMJfo9+R1uyf0OL8A8jDcjLBKzOrxXoUwCD6q++m4asqibtflSs/T4/sEwj4EX4RcL0/o9i6BPJx+GpI5Oku/+2zdt+t0/8+FZGa+ATm3j+nrV0CGzI1BSo5fgFxjTyLMFK1HmKM+BAm8T+m5ux9SCtoPKbVCv/e4/n03pFR4GGEQOEu3C8h1kIfc+Cd1O+fp+1yHrash15zrTLMH8pDYCQnOd+l2LobcW64J5TZIhsZV31+NsMOXawJxVfP3Qa6HvQg7av0P3c4v6jl5DeTcPoSw89gJPSYZhD2bz4LcT/OQEqlrypjX83ZEt3EUkqGoanoCPQcndF/dwvCuyvlVCKu7n9bPu06ZFpIpPRsSPN8Oub8uRVhzM6Hb+RmENRzrIffxKxGWrjfosXfXrKu9ej7k4TsBuT5fAznnH8Jyu3U8rwsyh/UYHkE4XPJ+TdMX9bxkIM+Hfwu5ft6g6bkS8jx6BaRUdpbu7yOQ+8tAhmz9gabNVRW/H8AHIddqXr/jewhriur5C2vt3QBgjHl5g/dEfRtSdT8LeTZdAXk+fhbAByDn0HWUzGo67oA8714BuVb2QDJOrlbQ9cp/TD/zhO7vFyHX55Qew8cQdtJ9N2RIVVm3cQJy3VvIc9M1bz6un3f9EHZCnplHIPfNekiQLUIKkR+HXCs/q9u5EtKf6BWaZleq3qDbm9T0btTj8DVN5z2tDuRAg7K19l4AMMZcDeAfrLWfNsa48ap3QXbeWGu36/uuj1wY74AEkKgDkAfI+yG50Nsh1QhXAfjv1trbjDErIDdKGRIAd0Auyl0Atlpr32eM+UnoWDq9gX4FEpx+FXJjHNPtrtX3vMsY8+8gF8jTkJv/DyAXYNVa+xZjzPWQTij3QtpgPgk5Wb8LuXC/BeBWa23VGPOHkAwJIOPc7tJ9cGNqXQbhQX3/pwF8WH82CHtk/oQejwxkDPbfIiw1nYOwg0jGWvuL7iDq/l8VnVzBGPMXkIeRqwqaQLiO9jzkwr0dUoX8GYRjOD+j6T2m++gyM99zm4ZcyIDcWP+gr23T1xYQVucWIRe2+8wr9LvfBznfWwB8z1r7Bt2H34BUsz1hjHknJFP2+5Dz/guRdD4JeYB+EvKQ/yuEHT22ArjMWnvCGPPTCK+d6xCWXI9DbrCv6D78FORBYQE8bq39kDHmNgBHrbUPGmOOQB6K10Ju4Mch1+6Cft/Dur8uQwJIUPlZ/fluyE3+Aj22v6zH4M36vvfrtfBha+17jTG/A+A5a+2NxpibAOy01n7eGOOGhdwBedA/D5LxOALtERy5337FWvs3+vPLrbUP6s8HrLV3G2NeBQly+yGdggAJilfrzx/Tc/pDPUYLkMD+a5ASs3MWgA16r/6qHsvP6/G9Wn//K0gQ3Iqw9PavAfwvyMPU2QAJeAf09austZ80xvxAz10Wch/dCMmknAXJlH0ZEtQv0/ROQzJQgFSXHoJU5T6KMMP9CgBf1mPxCwB+qB143hdJD6y1t+mP0SBXAPAd3ec367H/oLV2R2RbD0aO9W16rk7ovv0c5H7/KoDnW2v/xgVO/dyygOpe07/fbYz5bchz5fs4PRPh0n135OcH672n5v33GmNWRb7jVwCss9b+uTHmCkhJ3d1nbwTw6eh2jTE36vH4QwD/LySw/hLkWZaHBLMxyHPsu/oMvk2Pwxsgz++Teny+BMkA3gS5Pu+HZJr2QO7lirX2A/qM/gikme3HIfen6zM0C8lMboE88z9rjDkEaUr4Tciz8QuQ6+jnILVPayGjV35W0+3u01+GBOwfWWvds60hrhJFRETkCc59TURE5AkGZSIiIk8wKBMREXmCQZmIiMgT/xeAweyUFODDcAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(8,8)\n",
    "res = scipy.cluster.hierarchy.dendrogram(hc)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many clusters are there if we cut the dendrogram at distance of 0.65?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "62"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(numpy.unique(scipy.cluster.hierarchy.fcluster(hc, 0.65, 'distance')))"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}