3_WMD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Word Mover's distance\n",
    "\n",
    "In this note notebook we will see an application of Optimal Transport to the problem of computing similarities between sentences and texts. The method under the lens is called 'Word Mover's Distance' in reference to 'Earth Mover's Distance', another name of the Wasserstein $1$ distance, mostly used in computer vision. \n",
    "\n",
    "Traditionnally, portions of texts are compared by Cosine similarity on bag-of-words vectors, i.e. histograms of occurences of words in a text. It captures the exact similarity in terms of words, but two very related sentences can be orthogonal if the words that are used have the same semantic but are different. Such a semantic distance can be obtained by using *word embeddings*, that are embeddings of words in a Euclidean space (of potentially large dimension) where the Euclidean distance have a semantic meaning: two related words will be close in such embeddings. A popular embedding is the *word2vec* embedding, obtained with neural networks. A study of those mechanisms is not in the scope of this notebook, but the interested reader can find more information on [the corresponding Wikipedia page](https://en.wikipedia.org/wiki/Word2vec). Throughout the rest of this tutorial, we will use a subset of the [GloVe](https://nlp.stanford.edu/projects/glove/) embedding.\n",
    "\n",
    "The key observation made by Kusner and colleagues [1] is that when confronted to a sentence/document, the optimal transport distance can be used between histograms of occuring words using a ground metric obtained through word embeddings. In such a way, related words will be matched together, and the resulting distance will somehow express semantic relatedness between the content.\n",
    "\n",
    "[1] Kusner, M., Sun, Y., Kolkin, N., & Weinberger, K. (2015, June). From word embeddings to document distances. In International Conference on Machine Learning (pp. 957-966). http://proceedings.mlr.press/v37/kusnerb15.pdf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A basic example \n",
    "\n",
    "We will start by reproducing the Figure $1$ in the original paper\n",
    "\n",
    "<img src=\"https://remi.flamary.com/cours/otml/fig1.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two sentences are considered: 'Obama speaks to the media in Illinois' and 'The president greets the press in Chicago'. It is clear from this example that the Cosine similarity between the two sentences indicates that the two sentences are totally not related, since there is no word in common.  We will start by some imports and creating a list of the two sentences as words without stopwords that are not relevant for our analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "\n",
    "\n",
    "s1 = ['Obama','speaks','media','Illinois']\n",
    "s2 = ['President','greets','press','Chicago']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use a subset of the GloVe word embedding, expressed as a dictionnary (word,embedding) that you can load this way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "model=dict(np.load('data/model.npz'))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the embedded representation of the sentences can be obtained by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1_embed = np.array([model[w] for w in s1])\n",
    "s2_embed = np.array([model[w] for w in s2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the multidimensional scaling method in Scikitlearn, try to visualize the corresponding embedding of words in 2D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAGpCAYAAADvIG1RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu8HlV97/HPNwkQDREVFUGKHEDkoiIIqCFookY5iOCtjajVqLVejse7xWuJd3ujVD3aYjkgLZQcQSpUEJBbhYjKQXoKiDZoEMNNVCAEAiT5nT9mtjw8eXb25MK+ZH/er9fzmuw1a9asZyv5ZtZaM5OqQpIkjWzKWHdAkqSJwtCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSNrEkOyepJCf2lZ/Ylu88Ut3RMKg/WjdDU5I2Qhs63vA+SUwb6w5I0iS3DNgTuHMMzv0R4AttH9SBoSlJY6iqHgCuG6Nz3wzcPBbnnqgcnpWkMdRl/jPJ25L8Z5KVSW5NclySbYZp71lJTk9yW5L7ktyQ5CtJth9Qd+CcZpLDk1yQ5Oa2jZuSXJLknZvwq09IXmlK0vj2l8BLgLOA84C5wFuB3YAX9FZMchhwOhDgNOAG4FnAO4AjkhxUVUvXdbIkfwr8A3BLe87bgScAzwDeBHxlE32vCcnQlKTx7TnA06vqlwBJpgEXAnOTHFhVP2zLtwZOpPl7fU5VfW+ogSRH0cxdHge8eITzvQ24H9inqm7r3ZHkcZvkG01gDs9K0vj2qaHABKiqVcAJ7Y8H9tQ7AtgWWNQbmK2/AZYC85Ls1OGcq4AH+gur6vb16PdmydCUpPHtigFlN7bbx/SU7dduL+yv3Abtv7c/7jvC+U4GHglck+Rvk7w8yePXo7+bNUNTksa3OwaUrWq3U3vKhhYGDbcadqj80es6WVUdA7wR+CXwbuAM4NYkFyXZv1OPN2OGpiRtHobu83ziMPu376s3rKo6qaqeQzPc+1LgeOB5wLlJnrCxHZ3IDE1J2jz8uN3O6d/RLh6a3f54ZdcGq+qOqjq7qt5Ks8joscDBG9fNic3QlKTNw78CvwWOTPKcvn3vBXYBvtu7qGiQJIe0Idtv6Arzno3u6QTmLSeStBmoqruTvBn4BnBJkm/QzEs+i+Y2k1tobicZyanAyiSX0qy4Dc3V5QHA/wW+u+l7P3EYmpK0maiqbyU5CPgozQMRtqEJy78HPl1VN3Vo5sPtsfsBhwIraR6ScBTw1faxf5NWqnw4vyRJXTinKUlSR4amJEkdOacpSaMoyUxg/jTYbRUsoXns3fKx7pe6cU5TkkZJktnT4ey5MOUgmHEZrLgI1qyEQ6vq0rHun0ZmaErSKEgyczosOxNmzuspPx84HJavhB2q6u6x6p+6cU5TkkbH/LkwZV5f4TxgbvN38fwx6JPWk6EpSaNgGux2EMwYtG8WzJgKu452n7T+DE1JGgWrYMllsGLQvsWwYjVcP9p90vpzTlOSRoFzmpsHQ1OSRknv6tlZMGOxq2cnHENTkkZRkq2B+VNh13ZIdpFXmBOHoSlJUkcuBJIkqSNDU5KkjgxNSZI6MjQlSerI0JQkqSNDU5KkjgxNSZI6MjQlSerI0JQkqSNDU5KkjgxNSZI6MjQlSerI0JQmgSQ7J6kkJ67HMQvaYxY8fD2TJhZDU5rAkuyR5EtJrk5yZ5L7k9yU5NtJ3pJk+lj3UdqcTBvrDkjaMEn+HDia5h+/lwNfB+4GtgPmAP8IvAPYfwNPcUbb7s0b21dpc2FoShNQko8CnwRuBP6wqn4woM5hwAc29BxVdSdw5wZ3UtoMOTwrTTBJdgYWAg8Ahw4KTICq+jfgkEHHJzk1ye1JVia5og3Y/nrDzmkm2THJF5P8V9vGb5P8MMkn+urNTXJckmuT3JXk3nYo+ejhho6TbJ/khCS3tfWvSvLGJHPa/iwccMxTkpyUZFnPEPVJSZ4y6BzShvJKU5p43gRsAZxaVVevq2JV3ddX9GTgh8DPgX8CHgvMB76V5EVVddFIJ0+yP3Bue+y/A98EHgnsRRPmn+6pfhSwB7AY+DYwHTiorTenPefqnraf0NbduW17MfBE4CvAecP05wDgu8BM4Ezg2vacrwOOSPLCqrpipO8ldWFoShPP7HZ7wQYcOwdYWFWfHCpIcgrwHeBDwDpDM8mWwDdoAvN1VXVK3/4/6DvkncAvqqr66n0a+DjwamBRz67P0wTmX1bVUT31j6UJ+/7+BDgJeBTw+qo6uWfffOBU4J+T7FVVa9b13aQuHJ6VJp7t2+2vNuDYG4DP9BZU1bnAL4EDOxz/MppQO7M/MNu2buz7+ef9gdk6tt2+ZKigDeQjaeZR+/v4HzTh2G8WzVXl93sDsz1mEXAp8FQe/IeGtFEMTWniSbsdFEYjuap3OLTHjcBjOhz/nHZ7TpeTJZmR5KNJftTeErMmSQG3t1We1FP9qcAjgP9XVcsHNHfpgLL92u2Fw3RhqHzfLv2VRuLwrDTx3ERzdbXjBhx7xzDlq+j2j+hHt9tlI1VMsgVNaB0IXE0zDPtrmgVM0Nwus1XPIdu021uHaXJQ+dAxw90WM1T+6GH2S+vF0JQmnkuBFwAvBI4f5XMPhe6T1lmrcQRNYH69qhb07kiyPU1o9rqr3W43THuDyoduiXniMMds31dP2igOz0oTzwk0V2uvSrLXuiom2Wpd+zfA5e32v3eou1u7PX3AvucPKLsOuBd4RpKZA/YPmpf8cbudM0wfhsqvHGa/tF4MTWmCqaqlNLdsbAl8u70FZC1JDqHj3ON6OAtYChye5MgB5+y9Al3abuf01dkF+Iv+Y6vqfpoh3G1oVtb2HrMP8IYB/bkM+CkwO8mr+455NfA84GcMng+V1pvDs9IEVFWfSzKNZojzR0kWA1fw4GP0ngc8pS3blOe9P8kf0twzeUqSt9FcfU4H9qQZMh76e+UsYAnw/iRPp7kq3Ak4jOaezZ0GnOLDNEPPf5bk2TT3aW4P/BFwNvBy4Pe3jlRVJXkjcD6wKMm3aK5Yn9rWXQ68wdtNtKl4pSlNUFX1KeBpwJdprs7eRHOv5UuB64E/4WG41aJ9UMAzga/SPCzh/cAf0yy2Obqn3gqaADwF2Bt4N/AMmocfvH6Ytm+luY3kpPaY99GsfH0nMHRLyV19x/wAOKA9z3NpfgezgH8BDhjuiUnShsjgW6gkaXxJ8lngo8Ah7b2l0qgzNCWNK0l2qKqb+sqeTjNUez/wpKpaOSad06TnnKak8eaKJEto7u1cQTM3+1Ka6aS3G5gaS15pSpNAewvHfJi2G6xaAiwa5qk7Yy7J0TSLeHameQj7HTSLjf66qi4eu55Jhqa02UsyG6afDXOnwEEz4LIVcNEaWHloVXkrhrQeDE1pM9ZcYU5fBmfOhHk9e84HDl8OK3eoqrvHqn/SROMtJ9LmbX5zhTmvr3geTTnzx6BP0oRlaEqbtWm7NUOyg8yaAVN3Hd3+SBOboSlt1lYtaeYwB1m8AlZfP7r9kSY25zSlzZhzmtKmZWhKm7mHrp6dNaO5wnT1rLQhDE1pEkiyNTC/mcNcfT3NfZpeYUrrydCUJKkjFwJJktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhKktSRoSlJUkeGpiRJHRmakiR1ZGhOAkl2TlJJTkyyR5J/TfLbJCuSXJrkxX31F7T1FyQ5JMnFSe5MUn319mjbvDHJfUluTXJKkqcO6MN2Sf46yU/b897R/vnEJLv01EuSNyZZnOTXSVa27Z+bZP7D91uSpJFNG+sOaFT9N+D7wNXAPwDbA/OBc5K8tqoW9dV/NXAIcA7w98DOQzuSHAJ8E9gCOAtYAuwIvBJ4aZK5VXVlW/eRwGXArsD5bf0ATwaOAE4Dft42/VngI8AvgP8D3Nn28wDgD4H+PkrSqElVjVxLE1qSnWlCCOCvq+pDPfv2pwnSu4EnV9VdSRYAJwAFHFpV3+lr7zE0IbcaeF5VXduzb2/gB8DPqmq/tuxlwJnAsVX1vr62tgS2qqrl7c+/Ae4Fdq+qe/rqPq6qbt+IX4UkbRSHZyeXO4FP9RZU1RXAycCjgVf01f9Wf2C23tDWP7o3MNv2rgG+BuybZK++4+7tb6iq7h8KzB4P0ARyf10DU9KYcnh2crlyQEABXAy8EdgX+HpP+Q+Haee57XafJAsH7N+93e4JXAtcAiwDPpxkP+BsmuHaq6qqPxxPBv4ncE2Sb7THfr+q7lzH95KkUWFoTi63DlN+S7vdZpjyftu227eOcL6tAdoh3+cAnwQOB17S7r89yVeAz1TVA23Z+4DrgTcDH24/q5KcDXygqpaMcE5Jetg4PDu5bDdM+RPbbf/V3HAT3kP19qmqrOPz+6vWqvpVVb0FeALwNODdwG+AP28/Q/VWV9XfVdU+bX9fBZxBE7bfSbJV968rSZuWoTm57Jdk5oDyOe32xx3bubzdHry+HajGNVX1JWBeW/zyYereVlXfrKo/Ai6kWX37tPU9pyRtKobm5LINPVd18PvVs6+juXo8o2M7JwB3AEcnObB/Z5IpSeb0/Py0dgVvv6Er33vaelsleWGS9LW3BfDY3rqSNBac05xc/h34kyTPplmIM3Sf5hTgbVV1V5dGquo3SV5NE7KXJ7kAuAZYA+xEs1BoW2B6e8iLgGOSLAauA26juafziPaYv2rrPQL4LrA0yQ+AG9o25tEsKjqzqn6y4V9fkjaOoTm5/AJ4O/CFdrsVcCXwqao6d30aqqoLkjwD+CDNwp6DgfuBm2iGUk/vqX4ucCzwPJqgfBRwM82DDo6pqsVtvRXAUcBcYBbNsO1ymoVB7wD+9/p9XUnatHy4wSTQ83CDr1fVgjHtjCRNYM5pSpLUkcOzeoh2de18mLYbrFoCLBrmgQiSNOkYmvq9JLNh+tkwdwocNAMuWwEXHZPk0Kq6dKz7J0ljzTlNAUNXmNOXwZkzH7x9Epq1Oocvh5U7VNXdY9U/SRoPnNPUkPnNFea8vuJ5NOX4LktJk56hqda03Zoh2UFmzYCpu45ufyRp/DE01Vq1pJnDHGTxClh9/ej2R5LGH0NzAyVZmmRpX9mCJNW+xHmiWQQXrWnmMHudT1POojHokySNK66eFQBVtTzJoXB4u3p21ozmCvOiNbDyUBcBSZKrZzfY0FVmVe3cU7YNzfNcb56oL01OsjUwv5nDXH09zX2aBqYkYWhusEGhKUnavI3rOc0kO7dzhCcm2TXJaUl+k2R5kvOSPK2t9/gkxyW5OcnKJD9KMndAe9OSvDPJ5UnuSnJPkh8neVeStX4XabwryTVtu8uSfLm9ohzU34Fzmknmtv27tj3vvUmuTnJ0kumD2pIkjT8TZU5zZ+AHwE+AE9ufXwFcnOS5wHeAu2gWqzwWeA1wTpLdq+qX8Pt3Mp5F80aOnwKnACtp3qjxJeDZwB/3nfdY4N00b+Q4DniA5i0dzwa2pHmrRxdHAXsAi4Fv07zu6iBgITAnyYuqanXHtiRJY6Wqxu2HJhyr/Xysb98n2vLfAn8PTOnZ98ftvr/tKVvYln0JmNpTPhU4vt13RE/5rLZsCfDYnvLpwPfbfUv7+rSgLV/QV74L7VB4X/mn2/rzx/p37cePHz9+Rv6M6+HZHktp3gHZ6+vtdivgQ1W1pmffKcAq4JkA7dDru4BbgPdVz1Vd++cP0ITX63raeFO7/WxV/ban/krgI+vT+ar6eVUNmjw+tt2+ZH3akySNjYkyPHtVrT18eVO7/Vn1vYWjqlYnuRXYsS3aHdgW+C/g40kGneNeYM+en/drt5cMqPs9mlDuJMkM4D00Q8q7AzOB3k48qWtbkqSxM1FCc63bN6pqVRt+w93asQrYov3ztu32KcDR6zjP1j1/Hlrsc+uAc69O8pt1dXhIO5d6IXAgcDXNvOuvaeZHafuzVZe2JElja6KE5sYaCtYzquqV63nMdsDPe3ckmUoTxMs6tHMETWB+vaoW9LWzPesOcUnSODJR5jQ31nXAHcBz2iu/Lq5st88fsO9guv+DY7d2e/qAfYPaliSNU5MiNKtqFc2q2e2BLyZ5RH+dJNsn2aun6MR2+7Ekj+2pNx34/Hqcfmm7ndN3vl2Av1iPdiRJY2yyDM9Cc3vHPsDbgZcluZBmePUJNHOdBwEfA64FqKrLknwJ+J/A1UlO48H7NH9Hc+9mF2fR3Lby/iRPB34M7AQcRnPP5k6b5NtJkh52k+JKE6CqHgBeDryB5uEGh9HcanIIze/hE8DJfYe9hyY07wTeBhwJnAu8iI4PNqiqFcALaG6D2ZvmYQnPoAnx12/Md5IkjS6fPStJUkeT5kpTktRN+wzti9ej/pz2mIUPX6/GB0NTkkZBGyq9n9VJbk9yYZLXjdyCBmlf6FFJdh6N8202C4GSzATmT4PdVjULbxb1PylIksaBT7bbLYCn0qy1mJvkWVX1/rHr1kPsCdwz1p0YjzaL0EwyezqcPRemHAQzLoMVF8ExSQ6tqkvHun+SNKSqFvb+nOSFwPnAe5N8saqWjkW/elXVdWPdh/Fqwg/PJpk5Hc4+E2aeDTM+BpwNM86EmdPh7CRbj9iIJI2RqrqA5gEsAQ6Atd4lvHuSRUluS7ImyZyhY5M8Nsnnk/ykfU/vnUkuSPLi/vMk2TLJu5NcmeR37fuElyb5VpIX9dUdOKeZZLskxye5tT3fVUneuK7vt559/P07idv3EF+c5v3JdyX5dpI9++oXMHT+X/QMfS9dV582xuZwpTl/LkyZ11c4D5gLU86B+TSv/pKk8WroBQ79tzPsSvMu4Z/R3BL3CJp3B5PkycDFNK9Q/B7Ne4Vn0NxO950kb6uqr/W0dSLNbXNXAyfRvKRiB2A2za13311nB5Ntad4JvAtwafvZnubVjOcNc8z69nHIYTT3xJ/Ttr8XcChwQJK9qur2tt4naYa39wH+jubJb/RsN72xfjfZxn6mwRc+A1UDPp+GmgqfG+s++vHjxw/tu4EHlL8IWNN+ntyW7TxUn2H+DqMJozXAa/rKHw1cRROK27Vl27R1r6DnfcI9x2w7oK8X95UdR997itvy/Wke/FLAwg3tY1u+oG1nFfDCvmM+3+77s77yE9vynUfjf8cJPzy7CpZcBisG7VsMK1bD9aPdJ0kaTpKF7eez7ZPGvkNzpXlsVd3QV/1WHlw41NvGPjTPrj69qk7t3VdVd9C8CGI68Kqh4vYc99GEGH3HrPOtTe0zu18HLAcW9h17BWs/GGZD+tjr1GqGrXsd124PXFdfH26bw/DsoovgmPNphmSHnA9c1PyfY9HYdEuSBhp6s1HRDCN+Dzi+qv55QN3/qKr7BpQ/t91uM8y9kY9vt3sCVNVdSc4CXgZcleT09rw/qKouq2T3AB4JfK+qBr2O8WIenFvcoD72uWJA2Y3t9jHr7OnDbMKHZlUtT3Lo4e3q2VkwY3GzenbNSji0qu4e6z5K0pCqysi1fu+WYcqH3hE8j4deL/TrXQg5HzgKeC0PXr2ubK92P1hVa707uMew7xdeRz83pI9D1pqTrAffoTx1HW097CZ8aAJU1aVJdjgH5p8Hu7ZDsosMTEkT3HDPOR262ntPVX2xU0NV99IMrS5M8gfA82jmEF9PM4d68DoO732/8CBP3BR9nAgm/JzmkKq6u6qOX1X10ao63sCUtBm7vN2uK+iGVVU3VtXJwEuA/wJmt6tjh3MdzcMOnplkmwH752zqPq6H1e12VK5AN5vQlKTJol188z3glUnePKhOkqcneUL758cnefaAajOAmTSrVYd9c1M1b4k6ua27sO88+9MsEtqoPm6EoUVMo/Kaxc1ieFaSJqHXAhcCxyd5N839nHcAO9K8fvBpNItxbgOeBFye5CfAlTSLah5Fcz/kE4Ev1siPHf0o8EKaJxftz4P3ac4HzgYO38g+bqgLgA8BX2vnZ+8G7qiqL29Em8MyNCVpAqqqXyV5Fs07f19Fc7U3lWZRzrXAl4D/bKsvpVm1OweYCzwO+C3Nu4U/DDzklpBhznd7koOAz9Gswt2/Pf4dbftrheZ69nGDVNW5ST4AvBV4H7AlcAPwsISm79OUJKkj5zQlSerI0JQkqSPnNCVpHBp6RzBM2w1W+Y7gccI5TUkaZ5LMhulnw9wpcNAMuGwFXLQGVvqO4DFmaErSONJcYU5fBmfOXPuJ2ocvh5U7+PCWseOcpiSNL/ObK8yBbwme0uzXWDE0JWlcmbZbMyQ7yKwZMHXX0e2PehmakjSurFrSzGEOsngFrPYdwWPIOU1JGkec0xzfDE1JGmceunp21ozmCtPVs+OBoSlJ41CSrYH5zRzmat8RPE4YmpIkdeRCIEmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OapJLMSVJJFo51X6SJwtCUJKkjQ1OSpI4MTUmSOjI0pY2U5PAkFyS5Ocl9SW5KckmSd/bUubidP9wqyWeS/KKte32So5NsOUzbeyQ5McmNbf1bk5yS5KkD6u6e5AtJrkjy67b+DUmOS7Ljenyf6UlOa/v7v5JMacu3S/LXSX6aZEWSO9o/n5hklw353UkTzbSx7oA0kSX5U+AfgFuAs4DbgScAzwDeBHyl75D/AxwAnAY8ABwBLAT2T3J4VVVP24cA3wS2aNteAuwIvBJ4aZK5VXVlT9uvBN4OXAQsBu4H9gb+BHhZkv2ratkI3+cxwJnAQcBHquoLbfkjgcuAXYHz2/4EeHL7HU4Dfj7iL0ya4AxNaeO8jSac9qmq23p3JHncgPp7AntX1e/aOh+jCbnDgNcD/9SWPwb4F+Ae4HlVdW1Pu3sDPwD+Edivp+1/Av62qu7r68eLgXOAjwPvGO6LJHlyW2834A1V9c89u19IE5jHVtX7+o7bEthquHalzYnDs9LGW0Vz1fgQVXX7gLqfHgrMts5K4CPtj2/uqfcG4NHA0b2B2R5zDfA1YN8ke/WUL+sPzLb8POAa4CXDfYEkzwS+DzwJ+O99gdnr3gHt319Vy4drW9qceKUpbZyTgb8BrkmyCLgEuKyqfj1M/UsGlH2PJnj37Sl7brvdZ5j7KHdvt3sC1wIkCfA6YAGwD/AYYGrPMfcP06fZwPuB5TRXtf8xTL+XAR9Osh9wNs1w7VVVtXqYdqXNjqEpbYSqOibJ7cA7gXcD7wUqySXAh6rqir5Dbh3Qxuokv6GZCx2ybbt96whd2Lrnz8e0578ZOJcm5IauDBfQzD8Osi8wk2Ye9LpBFarqriTPAT4JHM6DV623J/kK8JmqWutqW9rcGJrSRqqqk4CTkjwamAW8gmao9dwke/bNdW4H/LL3+CRTaULyrp7iO9vtPlX1/0bqQ5In0IT21cCs/uHSJEeu4/AvA4+nme88M8nLq2rQMOyvgLe0V7R7AS8A/gfw5zRTPZ8YqZ/SROecprSJVNUdVXV2Vb0VOBF4LHBwX7XnDzj0YJp/wP64p+zynn1d7ELz3/N5AwJzx3b/Orpe7wSOBV4MfDvJjHVVrqprqupLwLy2+OUd+ylNaIamtBGSHJJk0IjN0FDrPX3ln2hXxg4dPx34fPvjCT31TgDuAI5OcuCA805JMqenaGm7nd1euQ7V25pm0dCIo0rtqtjPA3NprpIf1dPO05LsPOCw7dpt//eUNksOz0ob51RgZZJLaYIrNFeHBwD/F/huX/2f0Cwa6r1Pc1fg27S3mwBU1W+SvBo4A7g8yQU0K2DXADvRLBTaFpje1r8lyanAa4CrkpwHbENzJbgSuAp45khfpqo+mmQlzdzl+UkOaVf7vgg4JsnQvOdtNPeMHtH26a+6/sKkiczQlDbOh2kWxewHHEoTUDcARwFfHbA45o9o5v5eB+xAs1hnIfCF3gcbAFTVBUmeAXywPcfBNCtgbwIuBE7va/stNA8YmE8z1/hrmgcV/PmAusOqqk8luRf4S+CC9j7Pc2mGb59HE5SPollwdD5wTFUt7tq+NJGl779TSQ+DJBcDz6+qjHVfJG045zQlSerI0JQkqSPnNKVRlGQmMH8a7LaqeQD7Ih9BJ00czmlKoyTJ7Olw9lyYchDMuAxWXARrVsKhVXXpWPdP0sgMTWkUJJk5HZadCTPn9ZSfDxwOy1fCDlV191j1T1I3zmlKo2P+XJgyr69wHjC3+e9w/hj0SdJ6MjSlUTANdjsIBj6abhbMmNo84EDSOGdoSqNgFSy5DFYM2rcYVqyG60e7T5LWn3Oa0ihwTlPaPBia0ijpXT07C2YsdvWsNOEYmtIoat86Mn8q7NoOyS7yClOaOAxNSZI6ciGQJEkdGZqSJHVkaEqS1kuSBUkqyYKx7stoMzQlaZJLsn+SE5L8PMm9Se5K8p9J/irJk8aCbIGTAAALCUlEQVS6f+OJoSlJk1QafwH8CHg9cB3wReB44B7gg8DPkrx67Ho5vvhqMEmavD4B/BmwFDisqq7p3ZnkVcA/A6cmmVdVF41+F8cXrzQlaRJKsjNNaD4AHN4fmABVdTrwPmAq8NUka2VGkpcmWZxkRZLfJTktyVMG1Ns9yReSXJHk10nuS3JDkuOS7Dig/px23nRhO3z8nSR3tuc4PckftPV2SXJq2+a9SS5Kss/Gnn84hqYkTU5vohltPKOq/nMd9f4RuAl4KvD8vn2vBP4V+BXwd8D3gVcBlyd56oC6bwduBP4F+BJwLfAnwI/WMXd6APC99s9fA37YtnVBkj3an3cETgK+3fbx/PZBIpvi/A9VVX78+PHjZ5J9gAuAAt7aoe7Jbd2Ptz8vaH8ummHd3rrvacsv6Ct/ErDVgLZfDKwGvtpXPqfnHK/r23d8W/5b4GN9+z7R7nvPxpx/uI9XmpI0OW3fbm/sUHeozg595RdW1b/1lX2Z5hGRL0jy5KHCqlpWVff1N1xV5wHXAC8Z5tyXVtXJfWVfb7d3Al/o23dSu31m33k29PwPYWhK0uSUdtvlWarD1b2kv2JVrQaGXkCw7+8baLw+yXfbOcVV7ZxlAU+nuRIc5IoBZTe126va8/Va1m4fMk+5Eed/CFfPStLkdDOwB7BTh7pDAXRzX/mtw9S/pd1u01N2DPDeto1zacLt3nbfAuDJDHbngLJVw+2rqlVJALbo27Wh538IQ1OSJqdLgbnAi2gW2AyUZCrN/CLAZX27txvmsCe22zvbNp4AvBu4GphVVcv7znHk+nR8fW3K8zs8K0mT04k0C2BekWTvddR7M81c5k9Zezi2fzXtUMjObn/8cbvdhSZvzhsQWDu2+x9Om+z8hqYkTUJV9XPgczTDmGcm2au/TpKX09xKshp4Z1Wt6avygiSH9ZW9C9gVuKiqbmjLlrbb2W2oDrW/Nc1V7sM96rnJzu/wrCRNXguBGcD7gf9Ici7NStItgFnAs2nm/Y6sqgsHHH8WcEaSM4AlwD7AoTS3grxzqFJV3ZLkVOA1wFVJzqOZ75wHrASuom+166a0Kc/vlaYkTVJVtaaqPkATjqcAe9PM/f0psDXwN8DuVfWNYZr4JvAK4A9o7s88qC17blVd11f3LTRXto8A/gfNLR7/RhPOgxb7bGqb5Pxpb+6UJEkj8EpTkqSODE1JkjpyIZAkaURJZgLzp8Fuq5pFP4v6b9+YDJzTlCStU5LZ0+HsuTDlIJhxGay4CNashEOr6tKRW9h8GJqSpGElmTkdlp0JM+f1lJ8PHA7LV8IOVXX3WPVvtDmnKUlal/lzYcq8vsJ5wNwmQ+aPQZ/GjKEpSRrWNNjtoOYBCGuZBTOmNk//mTQMTUnSsFbBkstgxaB9i2HF6ubdmZOGc5qSpGE5p/lQhqYkaZ16V8/OghmLXT0rSdLw2jeCzJ8Ku7ZDsosm0xXmEENTkqSOXAgkSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpSVJHhqYkSR0ZmpIkdWRoSpLUkaEpTSBJ5iSpJAvHui/SZGRoSpLUkaEpSVJHhqYmtTTek+TaJCuTLEvy5STbJFmaZGlP3QXt0OiCJIckuTjJnUmqr809kpyY5MYk9yW5NckpSZ46TB8emeQjSa5KsiLJ3Um+n+TIvnonAhe1Px7d9mXoM6ets2WSdye5MsnvktzTfo9vJXnRJvzVSZPStLHugDTG/hfwDuAm4DjgfuBw4EBgC+CBAce8GjgEOAf4e2DnoR1JDgG+2R57FrAE2BF4JfDSJHOr6sqe+o8GLgT2Ba4E/jfNP2ZfApySZO+q+nhb/V/b7RuBS4CLe/q0tN2eCBwJXA2cBNwL7ADMbvv83U6/FUkDpapGriVthpIcDPw78DPg2VV1R1u+JU24HAzcUFU7t+ULgBOAAg6tqu/0tfcY4OfAauB5VXVtz769gR8AP6uq/XrKT6QJwaOq6i97yqfThOSLgf2q6qq2fA7N1eYnq2ph3/m3AX5HE77PrqrVffu3rarfrOevSVIPh2c1mb2x3X52KDABqup+4CPrOO5b/YHZegPwaODo3sBs27wG+Bqwb5K9oAkx4PXAFb2B2dZfCRwFBHhtx+9Tbf37gDVr7TQwpY3m8Kwms33b7aUD9l0OrBrmuB8OU/7cdrvPMLeE7N5u9wSuBQ4ApgLD3UKyRU/9EVXVXUnOAl4GXJXkdOB7wA+q6p4ubUhaN0NTk9k27fbW/h1VtTrJcFdmtwxTvm27fesI5926r/4B7Wek+l3Mp7lCfS3wybZsZZLTgA9W1VrfVVJ3Ds9qMrur3W7XvyPJVB4MtX7DLQS4s93uU1VZx+frffX/doT6c7t+oaq6t6oWVtXuwE40w7+XttvTurYjaTBDU5PZj9vt7AH7nsP6j8Rc3m4P7lj/hzRzj13rQ7PICJph3XWqqhur6mSalbj/Bcxu51ElbSBDU5PZSe32Y+3KU+D3q2c/twHtnQDcQXMP5YH9O5NMGbqfEqCqbgNOBvZP8okka4V0kl2T/LeeoqEh450G1H18kmcP6NcMYCbNHO396/F9JPVxTlOTVlVdkuQ44E+Ba9qFMw/QLKS5k+bezbVWoa6jvd8keTVwBnB5kguAa9o2dqJZKLQtML3nsHcBTwE+Bfxxkktp5lh3oFkAdADNfZe/aOv/FFgGvCbJ/cAvaYaL/wl4THven9DcdnIj8CjgMOCJwBerannnX5CktRiamuzeAVwHvA14O82V3BnAR4FfAdevT2NVdUGSZwAfpBkWPZjm6u4mmocYnN5X/64kz6cJ7tcCr6IJ1VtphlTfB5zfU391klcAXwD+iOYKMjTzllcBRwNzgLnA44Df0gTth4FT1+e7SFqbDzeQBkjyFJqHHpxaVUeOVF/S5OCcpia1JE9MMqWv7JHAse2PZ4x+rySNVw7ParJ7L3BkkouBm2nm/l5I87zYc4BvjF3XJI03hqYmu/OBfWie8fpYmhWmPwO+CBxbPfMXSWYC82HabrBqCbDIhTXS5OKcptRBktkw/WyYOwUOmgGXrYCL1sDKQ6tq0GP4JG2GDE1pBM0V5vRlcOZMmNez53zg8OWwcoequnus+idp9LgQSBrZ/OYKc15f8TyacuaPQZ8kjQFDUxrRtN2aIdlBZs2AqbuObn8kjRVDUxrRqiXNHOYgi1fA6vV6AIKkics5TWkEzmlKGmJoSh08dPXsrBnNFaarZ6XJxtCUOkqyNTC/mcNcfT3NfZpeYUqTiKEpSVJHLgSSJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOjI0JUnqyNCUJKkjQ1OSpI4MTUmSOvr/qcFGT3Lhe8oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn import manifold\n",
    "\n",
    "C = ot.dist(np.vstack((s1_embed,s2_embed)))\n",
    "\n",
    "nmds = manifold.MDS(\n",
    "        2,\n",
    "        eps=1e-9,\n",
    "        dissimilarity=\"precomputed\",\n",
    "        n_init=1)\n",
    "npos = nmds.fit_transform(C)\n",
    "\n",
    "pl.figure(figsize=(6,6))\n",
    "pl.scatter(npos[:4,0],npos[:4,1],c='r',s=50, edgecolor = 'k')\n",
    "for i, txt in enumerate(s1):\n",
    "    pl.annotate(txt, (npos[i,0]-4,npos[i,1]+2),fontsize=20)\n",
    "pl.scatter(npos[4:,0],npos[4:,1],c='b',s=50, edgecolor = 'k')\n",
    "for i, txt in enumerate(s2):\n",
    "    pl.annotate(txt, (npos[i+4,0]-4,npos[i+4,1]+2),fontsize=20)\n",
    "pl.axis('off')\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now compute the coupling between those two distributions and visualize the corresponding result \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAGoCAYAAAAkfL70AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xm8VVX9//HXG1BRJMcccEycB0xzCjFBpcwx0yRNC6fSvmVZmqaZaOO3X98y7dugX5MsTcqhNHFABBJxyNQcyAynFBVzQBFBBT6/P9Y6cticc+/mTufee97Px+M8Nnfttfde5wD3c9asiMDMzMxa16fRBTAzM+spHDTNzMxKctA0MzMryUHTzMysJAdNMzOzkhw0zczMSnLQtBZJGiPpt13wnI0lhaR+HXS/pyTtU+fccEnPVv38iKThHfHczpQ/n00bXY4KSZMlHZ///ClJtzS4PG9I2qSF83X/TZiV5aDZg0j6uqTxhbR/1Un7ZNeWrueKiG0iYnKjy9GTRcTlEfHhBpdh5Yh4AkDSWEnfbs/9JK0v6XJJL0uaK+keSQfkcxvmIF15Rc5T+XmPjnhP1v04aPYsfwF2l9QXQNI6wHLAjoW0TXPe0pT434MZIGl1YCrwNrANsCbwY+AKSYdFxL9zkF45IlbOl21flXZ7g4puncy/JHuWv5KC5Pvzzx8CJgH/LKQ9HhHPAUgaKumvkl7Lx6GVm+Xmte9IugN4E9hE0vskTZE0R9IE0i+LuiQdIOkBSbMlTZM0pOrcU5JOk/Rg/hZ+iaS1Jd2Y73+rpNUKtzxW0nOSnpf01ap79ZF0hqTH8zf/3+dfbJXzR0t6Op87q1DGFXPN41VJ04GdC+ffbbbLzdG/l3RZLuMjknaqyrujpPvzuT9IGlevRiNpsKTbcpleyrWWVQvPPTV/Pq/le/WvOn9a/hyek3RsK38Pq0u6NOd9VdIfq86dIGmGpFckXSdpUE5fqklcSza5jpZ0h6QLc/kelbR3neePljS16ueQdKJSq8erkv5XkvK5vpL+J38mT0r6QrEcVfc5RtL1VT/PkPT7qp+fkfT+qmduKumzwKeArynV+q6vuuX7633eBacAbwDHRcQLETEvIn4HfAf4n8p7sebjoNmDRMTbwN2kwEg+3k76Rlyd9hd499vyDcAFwBrAj4AbJK1Rddujgc8CA4GngSuAv5GC5beAz9Qrj6QdgV8Bn8v3/yVwnaQVqrIdCowENgcOBG4Ezsz37wOcXLjtCGAz4MPAGVrcB3Uy8DFgT2AQ8Crwv7kcWwM/z+9lUC7L+lX3PAcYnF8faek9ZQcBVwKrAtcBP83PWR64FhgLrA78DjikhfsI+F4u01bABsCYQp7DgX2B9wFDgNH5WfsCp5I+u82A1vrifgOsRKoVrUWqFSFpr1yGw4F1SX/HV7Zyr2q7Ak+Q/r7OAa6p/rLSigNIX1C2z8//SE4/Afgo6YvejqS/13qmAHvkL03rkr407g6g1H+5MvBg9QURcRFwOfCDXOs7sOp0zc+7hpHA1RGxqJD+e2BD0r9na0IOmj3PFBYHyD1IQfP2QtqU/Of9gX9FxG8iYkH+pvwoKXhVjI2IRyJiAemX6s7A2RHxVkT8Baj+ll50AvDLiLg7IhZGxK+Bt4DdqvJcGBGzImJmLufdEXF/RLxFCkA7FO55bkTMjYiHgEuBI3L654CzIuLZfO0Y4LBcOzkM+HNE/CWfOxuo/mV3OPCdiHglIp4hfYloydSIGB8RC0nBaPucvhvQD7ggIt6JiGuAe+rdJCJmRMSE/Fn+h/SlZc9Ctgsi4rmIeIX0WVdaDA4HLo2IhyNiLksH23flYPJR4MSIeDWXrfJv4FPAryLivvzZfB34oKSNW/kMKl4Ezs/3HEdq1di/5LXfj4jZEfFvUotI9Xv7Sf67fBX4fr0b5D7KOfnaPYGbgZmStsw/314jsLWk3uddtCbwfI3056vOWxNy0Ox5/gIMy82a742IfwHTgKE5bVsW92cOItUsqj0NrFf18zNVfx4EvJp/SVfnr2cj4KtKTbOzJc0m1aYGVeWZVfXneTV+XpklVZfn6ap7bQRcW/WcfwALgbVznnevy+V/ufC+ivdtyQtVf34T6J+D8yBgZiy5y8Ez1CFpLUlXSpop6XXgtyz9y7b4rMrnsSxl3gB4JQegoiX+DUTEG6TPZr0aeWspvt/qv5PWlH1vdT/DbAownPTFcAowmRQw92TxF8Sy6pWp6CXSl8iidavOWxNy0Ox57gRWITWp3gEQEa8Dz+W05yLiyZz3OVKwqbYhMLPq5+pfiM8Dq0kaUMhfzzOkGtyqVa+Vco22rTYoPPu5qmd9tPCs/rkG+3z1dZJWIjXRVr+v4n3b4nlgvUJ/1gb1MpOaRQMYEhHvAY4iNdmWfVbZMj8DrF7dX1pliX8D+e92DdK/gcqXo5Wq8q9TuL74fqv/TtrqeZZsPm/pM4TFQbPSijKF1oNme7dvuhU4VEsPjjuc9Hk/1s77Ww/loNnDRMQ84F7gK6TmzoqpOa161Ox4YHNJR0rqJ2kUsDXw5zr3fjrf+1xJy0saxpJNuUUXAydK2lXJAEn7SxrY5jcIZ0taSdI2wDHAuJz+C+A7kjYCkPReSQfnc1cBB0galvsdz2PJf9u/B74uaTVJ6wNfbGPZ7iTVbr+QP8+DgV1ayD+QNJhktqT1gNOW4Vm/B0ZL2jp/CTinXsaIeJ7UV/yz/B6Xk1Rprr8COEbS+3Nf83dJTeRP5SbjmcBReXDOsaR+32prASfne36C1Dc7nvb5PfAlSevlQH96K/mnkPq6V4yIZ0n/7vclBf/761wzC6g7Z7OEHwPvAS6RtI6k/pKOAM4CTivUvq2JOGj2TFNIv8ymVqXdntPeDZoR8TJpMMZXSU1yXwMOiIiWmpaOJA3+eIX0i/qyehkj4l5Sv+ZPSQNzZlB/YEVZU/J9JgI/jIjKhPmfkAbl3CJpDnBXLicR8QjwX6QA8Xwuy7NV9zyX1Kz4JHALqZ9ymeWBWB8HjgNmk2qOfyb149ZyLmmgy2ukAVnXLMOzbgTOB24jfR63tXLJ0cA7pD7rF4Ev5/tMJPXxXk36bAYD1XN4TyAF85dJg4imFe57N2kg0kukkaOH5X9X7XEx6e/hQVLQGw8sIH0hWUpEPEb68nF7/vl10uCkO3K/cy2XAFvn5vw/1slTV36Pw4D+wHTS5/MV4Ojct2tNSv7CZNZ2ku4GfhERlza6LB1N0mjg+IgY1snP+SjpMyx2JZh1O65pmi0DSXvm5rp+kj5DmrZwU6PL1ZMozZvdL3+G65FaNK5tdLk6khbPgR1bSB+b0zduLW9XqFUea5mDptmy2QL4O6nJ9auk5spaUxOsPpGarl8lNc/+A/hmQ0vUDjnouMmuSbh51sysHSoBMyJUlbYxqQ/91xExuip9LGlxjfdFxFM5bTlSX/NrXf0FLM/xXYW0itg7XfnsnqpDdpQwM7O2ycHq0QY9+3lqL+Jgdbh51sysgcr0f0r6nKSHJM2XNEvSRZJWqXO/D0i6WtKLkt5SWpP5Z7lWWcxbs09T0kGSJiqtffyW0prGUyR9vgPfeo/kmqaZWff2A9K6vdeTpuqMIE0V2hTYqzqj0tZlV5P6ja8iTbX6AHAScLCk3SvNwvUoLXj/S9LqSdeTphutRRr0dgzwsw56Xz2Sg6aZWfe2G7BdXsOXvKTjbcAISbtExD05fWXSZgL9gOHV25NJOp20xu9FpM0QWvI50pZo20fEi9UnJDX9mrtunjUz697OqwRMgLy5QmVecPWKVAeTVkkaV2M/z/8BngJGSiqzjOQC0mIZS2hlYZSm4KBpZta93VsjrbLIffV+tDvm41KrR+VAW1ktrLizUNHlpPWIH5H0Y0kfk/TeZShvr+agaWbWvc2ukbYgH/tWpVUGBtUbDVtJr7Ww/7si4kekaTH/Ju1jey0wS9IkVW3I3qwcNM3MeofX8rG4U03FuoV8dUXEZRGxG6m5d3/SWr4fAm6WtFZ7C9qTOWiamfUOlR1fhhdP5MFDlTWE7yt7w7yJ+PiIOIE0yGh10hZtTctB08ysd/gjaXeiIyTtVjj3ZdJWabdWDyqqRdK+OcgWVWqYb7a7pD2Yp5yYmfUCEfFG3hP1D8AUSX8g9Ut+gDTN5AXSdJLWXAnMlzSVNOJWpNrlzsDfSBt0Ny0HTTOzXiIi/iRpd+BM0oIIq5CC5S+Ab0XEcyVuc0a+dkdgP2A+aZGE04GfN/satV6w3czMrCT3aZqZmZXkoGlmZlaS+zTNzLqQpIHAqH6w6QKYQVr2bk6jy2XluE/TzKyLSBrWH8aPgD67w4A7YO4kWDQf9ouIqY0un7XOQdPMrAtIGtgfZl4HA99HWstuTWACcBDMmQ+DIuKNxpbSWuM+TTOzrnHU5rD894DNSBtWAowERqTfxaMaVjIrzUHTzKwTSdpe0oXAjx+EFZ4Cvg2MrsozFAb0hcFdXzpbVh4IZGbWwSStAhwBHE9akedt4IGdYchd0L9YW5kGcxfC411dTlt2rmmamXUAJcMkjSVtw/VzYHngS8AgYJ+H4J2JhesmAJNgETCuK8trbeOBQGZm7SBpbeDTwHHAFsAc4HfA/wH3RtUv2erRs0NhwDSPnu1xHDTNzJaRpL6kRdCPBw4idXXdQQqUf4iIuS1cuzIwqi8Mzk2y4zxqtudw0DQzK0nSxsCxwDHA+sBLwK+BSyLiH40rmXUVB00zsxZIWgH4GKn5dZ+cfDOpVnl9RLzdqLJZ13PQNDOrQdK2pEB5NLAGaW/KXwGXtraRs/VennJiZpZV1oUlBcvdgHeAP5JqlRMjYmEDi2fdgGuaZtbUJAnYlTSo55PAAGA6KVD+NiL+08DiWTfjmqaZNSVJa5KaXo8HtgbeBK4kBcu7wjUKq8E1TTNrGpL6kAbzHAccAiwH3A1cQpr68XoDi2c9gIOmmfV6kjYgTRM5FtgIeAX4DWmqyEONLJv1LA6aZtYrSVoeOIDU/LovIOBWUvPrnyJifgOLZz2Ug6aZ9SqStiQ1v34GeC8wE7gU+FVEPNnIslnP54FAZtbjSRoAfIIULIcBC4DrSbXKmz1VxDqKa5pm1iPlqSIfIDW/HgkMBB4jDeq5LCJeaGDxrJfy1mBmTUDSxpIib1tV9prR+ZrRnVeyZSdpdUlfBO4H/kraYeRa4EPAlhHxAwdM6yxunjXrwXL/3X8BI4ANgBVJi4jfD1wDXN4bBrzkqSJ7kmqVhwIrAH8DTgJ+FxGvNbB41kQcNM16KEnfBM4htRjdRdpt4w1gbWA4qT/vJGCnNj7i2nzf59tb1raSNAgYTeqr3AR4jfS+LomI+xtVLmteDppmPZCkM4FzgWeAT0TE3TXyHAB8ta3PyLW3Lq/BSVoO2I8UKPcnfSmYTPqCcHVEzOvqMplVuE/TrIfJezqOIS0mvl+tgAkQEX8mzU9c6npJV0p6SdJ8SffmAFvMV7dPU9L6ki6Q9K98j1ck3SPp7EK+EZIukjRd0uuS5kl6WNI5kvoX8m4q6XvAs6RF0vcHFgH/AMbm9DcljalRns0kXSZppqS3JT2Xf96s1mdj1lauaZr1PMeQln+7MiIebiljRLxVSNoIuAd4grQizuqkXT3+JGmfiJjU2sMl7UTaT3J14C+kvtOVSOu3jgG+VZX9dGBLYBpwA9Af2D3nG56D9cdIfZXDSUGy0gd7R36tA/wMuKVOeXYmLVowELiOtNj6lsCngIMl7R0R97b2vszKcNA063mG5ePENlw7HBgTEedWEiRdAdwEnAa0GDTzKjt/IAXMT0XEFYXzGxQu+TzwZHHxc0kXAScAs0i7ijwBnAVsR9pp5AcRcXpV/vNJwb5YHgGXAe8BjoqIy6vOjSItwP5bSVtHxKKW3ptZGW6eNet51s3HZ9tw7dPAt6sTIuJm0gbLu5S4/kBgY+C6YsDM93qm8PMTlYApaRVJJ0m6lxQwIQXNvYDNgB8CB5P6UYtl/DspOBYNJdUq76wOmPmaccBUYAsWf9EwaxcHTbOeR/nYlpVJHqizOs4zwGolrt8tH28s8zBJAyT9UtJLwGxSM+sHqrLMiIhJuRa4BWnKzIMRMafG7abWSNsxH2+rU4RK+g5lymvWGjfPmvU8z5FqV+u34drZddIXUO5L9Kr5OLOlTJLWJvW9nkPqx6wM6LmTFKDJ51aoumyVfJxV57a10ivX1JsWU0lftc55s2XioGnW80wlNWnuTVoyritVgu56xROS+gIfIQ3qOZDFv1+mAvtGxNyqvOuSgma1yl6Wa9d5dq30ypSYdepcs24hn1m7uHnWrOe5lDTd5FBJW7eUUdIKLZ1vg7vy8aNVz3ifpG+R+ktvIPUfng/8KGf5QXXAzPasce9HgXnAEEkDa5yv1S9ZWeBgeJ3yVtLvq3PebJk4aJr1MBHxFGnKxvLADXkKyFIk7UvJvsdlcD3wFHBQnqc5gcUjXx8EPgusHxGnkdaFhUJAk7QJ8N/FG0fE28A4UpPrNwrXbE9aY7boDuCfwDBJhxWuOYy0Hu1j1O4PNVtmbp4164Ei4ruS+pGaOP8qaRpwL4uX0fsQaURqR89P3JzUL7kR8EXgLVJAeozUx/rziLg4570emAF8RdJ2pFrhhqSNoW/Ify46g9T0/DVJu5Lmd64LHA6MJ83pfHfqSESEpM8AE4Bxkv5EqrFukfPOAT7t6SbWUVzTNOuhIuI8YFvgp6Ta2TGkuZb7A4+T+hY7YqpFf0nHS7oLeAg4jBT0ricNtNmVFKBWpaqfMjfJ7gVcAWwDnAwMIS1+cFSd9zSLNI3ksnzNKaSRr58HKlNKXi9cczewc37OB0mfwVDgd8DO9VZMMmsL76dpZkvJiwbsRgq8o0gLEEwnLZb+m4h4qQFl+g5wJmlQ0c1d/XwzcNA0syqS3gscTVosfWtgLmlVnf8D7i6u7NNJZRgUEc8V0rYjNdW+DazXG7Y7s57JfZpmTSCPRh0F/TaFBTOAcZUFBPJelfuQapUfI61rezdp1Z5xdRYa6Ez3SpoBPEwK2puxeLeTEx0wrZFc0zTr5SQNg/7jYUQf2H0A3DEXJi2C+aNJa70eQxrY8wqpL/GS1haC7+TynkMK3huTFmGfTZrq8sOImNyocpmBg6ZZr5ZqmP1nwnUDYSSpdfN64PtUDaydQFok4Y81dkUxsypunjXr3UalGua2wKmkiuR/SAv6DH4HHj87IpaaM2lmtXnKiVmv1m/T1CS7gDQzZQ/SbJGngdHLQd9VWrzczJbgmqZZr7ZgRurDPGtAWu+8OkZOmwsLH29Uycx6IvdpmvViS/dpVkwADpoD8wdFxBuNKp9ZT+OgadbLLTl6duiAVMOctAjm7xcRXpPVbBk4aJo1AUkrA6Og7+DcJDvONUyzZeegaWZmVpJHz5qZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTloNgFJG0sKSWMlbSnpj5JekTRX0lRJHy7kH53zj5a0r6TJkl6TFIV8W+Z7PiPpLUmzJF0haYsaZVhb0g8l/TM/d3b+81hJm1Tlk6TPSJom6T+S5uf73yxpVOd9SmZmrevX6AJYl3ofcCfwMPBLYF1gFHCjpCMjYlwh/2HAvsCNwC+AjSsnJO0LXAMsB1wPzADWBz4O7C9pRETcl/OuBNwBDAYm5PwCNgIOBq4Cnsi3/g7wdeBJ4PfAa7mcOwOfAIplNDPrMoqI1nNZjyZpY1IQAvhhRJxWdW4nUiB9A9goIl6XNBq4FAhgv4i4qXC/1UhBbiHwoYiYXnVuG+Bu4LGI2DGnHQhcB5wfEacU7rU8sEJEzMk/vwzMAzaPiDcLedeMiJfa8VGYmbWLm2eby2vAedUJEXEvcDmwKnBIIf+figEz+3TOf051wMz3ewS4GNhB0taF6+YVbxQRb1cCZpV3SAG5mNcB08ways2zzeW+GgEKYDLwGWAH4NdV6ffUuc8H83F7SWNqnN88H7cCpgNTgJnAGZJ2BMaTmmsfiIhicLwc+CLwiKQ/5GvvjIjXWnhfZmZdwkGzucyqk/5CPq5SJ71ojXw8oZXnrQyQm3x3A84FDgI+ks+/JOlnwLcj4p2cdgrwOHAscEZ+LZA0HvhqRMxo5ZlmZp3GzbPNZe066evkY7E2V6/Du5Jv+4hQC693a60R8WxEHAesBWwLnAy8DHwzvyr5FkbETyJi+1zeQ4FrScH2JkkrlH+7ZmYdy0GzuewoaWCN9OH5eH/J+9yVj3ssawEieSQiLgRG5uSP1cn7YkRcExGHA7eRRt9uu6zPNDPrKA6azWUVqmp18O7o2U+Rao/XlrzPpcBs4BxJuxRPSuojaXjVz9vmEbxFlZrvmznfCpL2lqTC/ZYDVq/Oa2bWCO7TbC5/AY6XtCtpIE5lnmYf4HMR8XqZm0TEy5IOIwXZuyRNBB4BFgEbkgYKrQH0z5fsA/xI0jTgUeBF0pzOg/M1/y/nWxG4FXhK0t3A0/keI0mDiq6LiH+0/e2bmbWPg2ZzeRI4Efh+Pq4A3AecFxE3L8uNImKipCHAqaSBPXsAbwPPkZpSr67KfjNwPvAhUqB8D/A8aaGDH0XEtJxvLnA6MAIYSmq2nUMaGHQS8Ktle7tmZh3Lixs0garFDX4dEaMbWhgzsx7MfZpmZmYlOWiamZmV5D5NW0KekjIK+m0KC2YA4+qsImRm1nTcp2nvkjQM+o+HEX1g9wFwx1yYtAjm7xcRUxtdPjOzRnPQNKBSw+w/E64bCLuQZpAMJQ1wPWgOzB8UEW80tpRmZo3lPk1DUh/ga7DBCmkTlDWAD5NmkIwk1TzxBtBm1vTcp9mkJK1Nioz75uOa8C/SFMozSFMv++bcQwfALYMbUlAzs27EQbONJD0FEBEbV6WNJi0xd0xEjG1EuerJmz0PJUXDfYH351MvAjcCb8FeR8LElZa+etpcWPh4FxXVzKzbctDsxSRtwuIguRdpq64FpCX0zgRuAv4eEYtSn+a0UakPc2TVXSaQBgMxrksLb2bWDXkgUBvVqWmuQlrP9flGbJosaQBpx5J9ScFys3zqSdJSdjcBk+qtMbvk6NmhA1IN06NnzcwqHDTbqFbQbEAZRNoqqxIk9wCWJ+0EMonFgXJGlPyLlrQyMAr6Ds5NsuM8atbMLOnWo2clbSwpJI2VNFjSVZJeljRH0i2Sts353ivpIknPS5ov6a+SRtS4Xz9Jn5d0l6TXJb0p6X5JX8gjSIv5lc89ku87U9JPc42yVnlH5/KOLqSPyOWbnp87T9LDks6R1L/WvVr4TNaQNErSr4BngQeBH5A2d76AtKPI6hFxQERcGBH/KhswASLijYi4JGLBmenogGlmVtFT+jQ3Bu4G/gGMzT8fAkyW9EFSbep1Ur/b6sAngRslbR4R/4Z392S8nlQj+ydwBTCftKPGhcCuwNGF554PnEzakeMi4B3SLh27kmp0b5cs/+nAlsA04AbSdle7A2OA4ZL2iYiFtS6U1A/YmcW1yV0AAa+SOhxvAm6JiJkly2JmZm0VEd32RQqOkV9nFc6dndNfAX4B9Kk6d3Q+9+OqtDE57UKgb1V6X+CSfO7gqvShOW0GqeZWSe8P3JnPPVUo0+icPrqQvgm5KbyQ/q2cf1QhfX3gOOD3pOAYwML83HOA3arfg19++eWXX13z6tZ9mlVbWj0FbBpVtTFJG5I2KX4TWCeq1keV1JdUi5waESNy0+uLpJriBhGxoPCcVUnB96qIODynXQwcDxwbEZcW8g8n9Rk+He2YciJpDeAl4DLgtywe6bpNzjKTxf2SEyPildbuaWZmnaenNM8+EEs3Xz6Xj49FYUHxiFgoaRapxgawOWmZm38B30jjZ5YyD9iq6ucd83FKjby3k6ZulJJHtX6J1KS8OTCQ1MRacRTwaVJz719Igfdm4JHozt9qzMyaTE8JmktN34iIBTn41ZvasQBYLv95jXzcjNS8Wc/KVX+uDPaZVePZCyW93FKBK3Jf6m2kvsh/k2qPg6ruD6nP9LPAlIiYW+a+ZmbW9br16NkOVAms10aEWni9r8Y1axdvlpt/1yimF/L0kfQB0sClXUj9khuSar+TgJNIA4ogTQkZ74BpZta9NUvQfBSYDeyWa35l3JePe9Y4twe1a+kD8/GzwAvAvcCROe2qfK81IuKQiPgFaYCQmZn1EE0RNPPAnwtJq/VcIGnFYh5J60rauippbD6eJWn1qnz9ge9V/bynpO9K+htpniSkBQduIY3iPTGnPRMRf4mId/J1mwD/3RHvz8zMukZP6dPsCN8Ctid5Z8b+AAAdbElEQVQFsQMl3UbqX1yL1Ne5O3AWMB0gIu6QdCHwReBhSVeRRt8emu83n9TcOpnUf3oncHU+/6XKiNs8COhU4CuStgPuz9cdQJqzuWGnvmszM+swTVHTBMg1vI+RRqn+kxS0vkqa4tGHNO/z8sJlX8p5FgCfzz9vlF/LAXNJI2LXiIgPAX+uPK7quXNJi6VfQZpKcjIwhBTEj+rgt2lmZp2oW8/TbISq9VwrcyYr67nOI9UqbyJNB3nM00HMzJqLgyaQ+yz3YfFSdYPyqUdYHCRvj4j5jSmhmVnXkRSkKXDDS+YfTpoVcG5EjOm8kjVer+nTTPtBMqofbLogLX03rrjoQVXeviy9nmsf0gjbCaQgeXNEPNslhTezXi8HomqLSMtkPghcEhHF7iErQdJY4DPA+yLiqc5+Xq8ImpKG9YfxI6DP7jDgDpg7CX4k6d19ICWtRwqQHyHtsrwaqe/xHlL/4s3AX4tL7JmZdbBz83E5YAvSWIsRkj4QEV9pXLGWsBVpiVIr6PFBU9LA/jD+Ohg4cnHygBuAQ+BmSf9HGoizbT73PPBHUpC8NSJKrexjZtYRis2XkvYmtXB9WdIFXVFbak1EPNroMnRXvWH07KgR0GckaTHaC4H9gU8A78BKpFGvs4CvkUatrhcRx0bEOAdMM2u0iJhIWoBFpG6j4l7Cm0saJ+lFSYty/yE53+qSvifpH3mf3tckTZT04eJzJC0v6WRJ90l6Ne8n/JSkP0nap5A3JE2ucY+1JV0iaVZ+3gOSPtPS+1vGMr67J3Heh3iy0v7Jr0u6QdJWhfxBapoFeDJfG5KeaqlM7dHja5r9YNPdYQCkjoGTSZMujyN1UF4B5y+MOK2BRTQza01lA4div+dg0l7Cj5GmxK1I2jsYSRuRRvRvTNpE4ibS78IDgJskfS4iLq6611jgCOBh0s5K80iDHoeRxnfc2mIB065M00grmU3Nr3VJWzPeUueaZS1jxQGkvYtvzPffGtgP2FnS1hHxUs53Lql5e3vgJ6Rf+1QdO16j9yZr7ws4/qPwRkDMg3gcIvLro/AGcFyjy+iXX375Rd4buEb6PqRBQYuAjXLaxpX8wHfr3G9yvuaThfRVgQdIQXHtnLZKznsvNfbiJc01L5Z1ciHtIgr7FOf0nUgLvwQwpq1lzOmj830WAHsXrvlePve1QvrYnL5xV/w99obm2XGTYNEE0u7QlcVcJwCT0l/WuIaVzMysQNKY/PpOXmnsJlJN8/yIeLqQfRaLBw5V32N70lrWV0fEldXnImI2aTen/ixewSzyM94i/V6kcE2LXVV5ze5PAXOAMYVr72XphWHaUsZqV0Zqtq52UT7u0lJZO1uPb56NiDmS9jsoj54dCgOmpdGzi+bDfhHxRqPLaGZWpbI9YZCaEW8nTTn5bY28f4+It2qkfzAfV5E0psb59+bjVgAR8bqk64EDgQckXZ2fe3dElBkluyVpjMjtEVFrO8bJLO5bbFMZC+6tkfZMPq7WYkk7WY8PmgARMVXSoBth1C0weCE8Tpqn6YBpZt1KRKj1XO96oU56ZWvCkflVT/UewaOA00k7L1Vqr/NzbffUiFhq7+AqdfcXbqGcbSljxVJ9krF4D+W+Ldyr0/WKoAmQA+QljS6HmVkHqrdkW6W296WIuKBOniVvFDGP1LQ6RtIGwIdIfYhHkfpQ92jh8rr7C2frdEQZe4Le0KdpZtZs7srHlgJdXRHxTKQViD4C/AsYlkfH1vMoabGD90tapcb54R1dxmWwMB+7pAbqoGlm1sPkwTe3Ax+XdGytPJK2k7RW/vN7Je1aI9sAYCBptOrbLTzvHdJgn4EUBgJJ2ok0SKhdZWyHyiCmLtlmsdc0z5qZNZkjgduASySdTJrPORtYn7SQy7akwTgvAusBd0n6B3AfaVDNe0jzIdcBLog6a3VXORPYm7Ry0U4snqc5ChgPHNTOMrbVROA04OLcP/sGMDsiftqOe9bloGlm1gNFxLOSPgB8kTRt41OkJsoXgOmkBdIeytmfIo3aHQ6MANYEXiHtLXwGsMSUkDrPe0nS7sB3SaNwd8rXn5Tvv1TQXMYytklE3Czpq8AJwCmkrRyfBjolaHprMDMzs5Lcp2lmZlaSm2fNzLqhyh7B0G9TWNDiHsHWddw8a2bWzUgaBv3Hw4g+sPsAuGMuTFoE89/dI9gaw0HTzKwbSTXM/jPhuoFLLqQzAThoDswf5NXOGsd9mmZm3cuoVMPcDfgccG1OHklKZ1TDSmYOmmZm3Uu/TVOT7EqktQG+weJFb4YOgL6DG1c2c9A0M+tWFsxIfZh9SYvvTGfxDofT5sLCxxtWNHOfpplZd7Jkn+bewPtJ22D+BDjEfZoN5qBpZtbNLDl6dtUB8Dug33xYMNKjZxvLQdPMrBuStDIwCvoMhkVHAouALfLi6dYgDppmZt2cpP2BPwOfjYiLG12eZuagaWbWzUkScCcwCNgsIt5qcJGalkfPmpl1c5FqN2cDGwDHN7g4Tc01TTOzHiDXNqcAmwKDI2Jeg4vUlFzTNDPrAapqm+sCJza4OE3LNU0zsx5E0q3AdsAmETG30eVpNq5pmpn1LGcDawFfaHRBmpFrmmZmPYyk8cCuwPsi4vVGl6eZuKZpZtbzfBNYHfhSowvSbFzTNDPrgST9ERhOqm2+2uDiNA3XNM3MeqZzgFWAUxpdkGbimqaZWQ8l6Q/AR0i1zZcbXZ5m4JqmmVnPNQZYGTitweVoGq5pmpn1YJIuBz5Gmrc5q9Hl6e1c0zQz69nOBfoDpze6IM3AQdOsSUkaLikkjWl0WaztIuIx4DLgJEmDGl2e3s5B08ys5zsP6Aec2eiC9HYOmmZmPVxEPAn8CjhB0oaNLk9v5qBpZtY7fDsfv9HQUvRyDppm7STpIEkTJT0v6S1Jz0maIunzVXkm5/7DFSR9W9KTOe/jks6RtHyde28paaykZ3L+WZKukLRFjbybS/q+pHsl/Sfnf1rSRZLWX4b301/SVbm8/yupT05fW9IPJf1T0lxJs/Ofx0rapC2fnXWciHgGuAg4xn8fncdTTszaQdJngV8CLwDXAy+RdqAYQvr/tXPONxnYE7gO2Bm4CngHOBgYDPwZOCiq/kNK2he4Blgu33sGsD7wceAtYERE3FeV/wzgDGAS8AzwNrANafL7LGCniJhZlX94zntuRIzJaavlMu4OnBkR38/pKwEP5rJOyH8WsBGwN3B0RPy57Z+kdQRJ6wJPAOMiYnSDi9M7RYRffvnVxhfwN1IAW6vGuTWr/jwZCOAxYLWq9P7Anfnc0VXpqwGvkoLw1oX7bgO8AdxXSF8PWKFGOT4MLAR+Xkgfnp87Jv+8ETCdFGyPKuQ9MOf9cY37Lw8MbPTfhV/v/n38T/773qLRZemNLzfPmrXfAlKtcQkR8VKNvN+KqsW1I2I+8PX847FV+T4NrAqcExHTC/d9BLgY2EHS1lXpMyPirRrluAV4hFTjrEnS+0nBez3goxHx2zpZ59W4/9sRMafeva3L/Tcwn7Q2rXWwfo0ugFkPdznpm/0jksYBU4A7IuI/dfJPqZF2Oynw7lCV9sF83L7OPMrN83ErUu0QSQI+BYwGtifVVvtWXfN2nTINA74CzAE+FBF/r1PumcAZknYExgN3AA9ExMI697UGiIgXJV0IfE3Sd/KXLOsg7tM0aydJnwY+T+qr7ENqxpwCnBYR9+Y8k0l9mv1r1QYlvUBq4q0MupkA7FPi8aMj4tf5mh8DXwaeB24jBblKzXA0sFFEqOqZw0l9mq+Q9ma8GTi4Vvly/vVJq88cBKyZk18CfgZ8OyKWqm1bY0haA3gSuDkiPtHo8vQmDppmHUTSqsBQ4BBSU+tsYKv8zX8yKWhuFBH/LlzXl9ScNjciVs1pVwGHAttHxIMlnr0WKVhOB4YWm0sl/RPYvE7QPA94L3AScAvwsYhYqhm26joBWwN7Af8FbEEKmme3Vk7rOpLOA84GdoiIBxpdnt7CfZpmHSQiZkfE+Ig4ARhLqr3tUci2Z41L9yB1ldxflXZX1bkyNiH9f76lRsBcP59voejxeeB80qChGyQNaClzRDwSERcCI3Pyx0qW07rOj0hf3M5tdEF6EwdNs3aQtK+kWmMD1srHNwvpZ+dpHZXr+wPfyz9eWpXvUtIvvHMk7VLjuX1yTbHiqXwclmuulXwrkwYNtTp+ISJOyWUZAdws6T1V99lW0sY1Lls7H4vv0xosImaT+tsPkrRzo8vTW3ggkFn7XAnMlzSVFLhEqh3uTJqOcmsh/z9Ig4aK8zRvAH5TyRQRL0s6DLgWuEvSRNII2EXAhqSBQmuQpqwQES9IuhL4JPCApFuAVUg1wfnAA8D7W3szEXGmpPmk2skESfvm0b77AD+SNA14FHiRNGf04Fym/1f2A7Mu9RNSP/d5wEcbXJZewTVNs/Y5gzRVY0fSYKBjSIsRnE5afKA4OOZw0hqhBwJfIP0fHAMcGoUBBhExkbRIws+AjYETgeOBbUkDfT5ZuPdxwHeBFUl9jR8hLZowFHit7BuKiPOArwG7ABMlrUkaJHQ+KUgfDHwV+BBpoYM9IuKqsve3rpOb6n8A7CtpaKPL0xt4IJBZF6gMBKoeiGPWFXL/9BPAwxGxd6PL09O5pmlm1otFxFxSX/VehX5wawMHTTOz3u8XwHPAt/KUIWsjB02zLiRpoKTjl0u7kRwvaWCjy2S9X16u8Tuk1Z9GtpLdWuA+TbMuImlYfxg/AvrsDgPugLmTYNF82C8ipja6fNa7SVoB+BdpEYzdigPPrBwHTbMuIGlgf5h5HQzcgrQg7CDSfJSDYM58GBQRbzS2lNbbSTqBtOfmgeGt3NrEzbNmXWPUCOgzkjRhrrIy+0hgRPp/OKphJbNmMpY0kvY89222jYOmWRfoB5vuDgMg7d68HWkVBIChMKBvWuDArFPlecOV721e+rANHDTNusACmHEHzF1IWtZnSNW5aSn98QYVzZrP5aTN0M+T5BiwjPyBmXWNcZNg0W9Ii7RWguYEYFJahm5cw0pmTSUiFpBWodoW8LZhy8gDgcy6iKRhy8Et78CKJwJPe/SsNUiuYT5IGpO2rTcSL89B06wLSfoucEYf+OEi+CcwzqNmrREkHQpcBXw6In7TWn5LHDTNupCka4CtI2LLRpfFmluubf4NGEjaLL24uYDV4D5Ns641hNQsZtZQEbEI+CZp5PanG1ycHsNB06yL5A2hB+Ogad3Hn4F7SJujL1/2IkmjJYWk0Z1Wsm7KQdOs62ybjw81tBRmWV5K75vARsBkSU9ImifpdUkPSfp/ktZrcDG7FQdNs66zXT66pmndQl4VaK/84wdJ8zcvAC4hzY46FXhM0mGNKWH346Bp1nWGAHOApxtdELPsbOBrwAv55xsj4vSIOCUidgUOI8WJKyWNaFQhuxMHTbOuMwR4OA/AMGsoSRuTguY7wIeBycDXJa1UyRMRVwOnkOZz/rzWCkKS9pc0TdJcSa9KukrSZjXyba60Jd69kv4j6S1JT0u6SNL6NfIPz/2mYyTtJOkmSa/lZ1wtaYOcbxNJV+Z7zpM0SdL27X1+PQ6aZl0gN4Nth5tmrfs4BugHXBsRD5EC6NrA5wv5/o+0gfUWwJ6Fcx8H/gg8C/wEuBM4FLhL0hY18p4IPAP8DrgQmA4cD/y1hb7TnYHb858vJg1c+jgwUdKW+ef1gcuAG3IZJ+SBdx3x/CVFhF9++dXJr/yfOoDPN7osfvkVEQAT87/JE6rSbgb+Awws5L085/1G/nl0/jmAAwp5v5TTJxbS1wNWqFGODwMLgZ8X0odXPeNThXOX5PRXgLMK587O577UnufXe7mmadY1KsvNuqZp3cW6+fhMVdrZwJrAFwt5K3kGFdJvi6X35fwpaQOCvSRtVEmMiJkR8VaxEBFxC2kfg4/UKefUiLi8kPbrfHwN+H7h3GX5+P7Cc9r6/CU4aJp1jUrQfLihpTBbrLI73bvLwkXEPaS5m6dKWqWlvNmU4k0jrWNbWUu5snUsSo6SdGvuU1yQ+yyD1HVRr3n03hppz+XjA7H0urkz83GJfsp2PH8J/cpkMrN22w74d0TMbnRBzLLngS2BDQvp3wTuIw0AGpPT1q+6ptqsOveujMatDrw/Ar6c73EzKbjNy+dGk+aK1vJajbQF9c5FxIK8v/ZyhVNtff4SHDTNuoaXz7PuZiowAtiHNMAGgIi4P6+RfIqkC0iBaXg+fUfhHmvXufc6+fgagKS1gJNJLS1DI2JOdWZJR7T9bbSuI5/v5lmzTpaXJ9sSB03rXsaSBsAcImmbwrlzSAu5nwocS+rL/CdLN8cWR9MiqS8wLP94fz5uQoo3t9QIWOvn852pw57voGnW+bYktep4+TzrNiLiCeC7pGbM6yRtXXXuYdLG6KeQppIsJI38Ls4x3kvSAYW0L5DWWJ4UEZWFPJ7Kx2E5qALvrsd8MZ3f6tlhz3fzrFnn88hZ667GAAOArwB/l3QzaSTpcsBWQH9S/+EREXFbjeuvB66VdC0wA9ge2I80FeTd+Z4R8YKkK4FPAg9IuoXU3zkSmA88QGG0a0fqyOe7pmnW+YYAb5PW9TTrNiJiUUR8FdgVuALYhtT391lSpWo6KWhOrXOLa4BDgA1I8zN3z2kfjIhHC3mPI9VsVwT+izTF48/AUGoP9uloHfJ8b0Jt1skk3QSsHRE7tJrZrBuRtCnwKPCziDi50eXpDlzTNOt8Xj7PeqSImEEaMPS5ylqvzc5B06wTSVqTNPLQQdN6qm+RFjc4q9EF6Q48EMisc3kPTevRIuJpSf8HnCDpyX6w2oI06GdccfpGM3BN06xzVUbOerqJ9WS3AP3Wg2+PgdM/Cuf3h5mShrVyXa/joGnWubYj7RpRb7kxs25N0sD+cNkhwAvQbxQwHgZcBwP7w/gaW3D1ag6aZp1rCPBgeJi69VyjRkCfnwHLkzo4IU1wHJFiyKiGlawBHDTNOkleeWRb3DRrPVg/2HR3GLAO8FsWB02AoTCgb1r9p2l4IJBZ59mENJHag4Csx1oAM+6AucCAjxfOTYO5C9PemU3DNU2zzuPl86w3GDcJFk0oJE4AJsEi0hq1TcM1TbPOM4T0S2V6owti1lYRMUfSfgfB+BHQZygMmAZzJ8Gi+bBfRLzR6DJ2JS+jZ9ZJ8iLWW0XElo0ui1l75VGyo/rC4NwkO67ZAiY4aJp1GkkzgPsi4vBGl8XMOob7NM06Qf5WPhj3Z5r1Kg6aZp1j23x00DTrRRw0zTqHl88z64UcNM06x3bAHODpRhfEzDqOg6ZZ5xgCPBQRixpdEDPrOA6aZh1MkshBs9FlMbOO5aBp1vHWA1bFg4DMeh0HTbOO5+XzzHopB02zjueRs2a9lIOmWccbAvw7Il5rdEHMrGM5aJp1vO1w06xZr+SgadaBJK0AbImDplmv5KBp1rG2JG255/5Ms17IQdOsY22Xj65pmvVCDppmHWsI8DbwWKMLYmYdz0HTrGMNAR6JiAWNLoiZdTwHTbOO5eXzzHoxB02zDiJpTWBd3J9p1ms5aJp1HA8CMuvlHDTNOo7XnDXr5Rw0zTrOEOA/ETGr0QUxs87hoGnWcbx8nlkv56Bp1gEk9QW2xUHTrFdz0DTrGIOBFfF0E7NezUHTrGN45KxZE3DQNOsYQ4BFwPRGF8TMOo+DplnHGAI8FhHzGl0QM+s8DppmHcPL55k1AQdNs3aStDKwCe7PNOv1HDTN2m/bfHTQNOvlHDTN2s/L55k1CQdNs/bbDpgDPN3ogphZ53LQNGu/IcBDERGNLoiZdS4HTbN2kCRS0HTTrFkTcNA0a5/1gVXxdBOzpuCgadY+Xj7PrIk4aJq1T2XkrGuaZk3AQdOsfYYAT0fEa40uiJl1PgdNs/bp0uXzJA2XFJLGdNUzzWwxB02zNpK0ArAF7s80axoOmmZttyXQDwdNs6bhoGlNTcmXJE2XNF/STEk/lbSKpKckPVWVd3RuGh0taV/gynzqysI9t5Q0VtIzkt6SNEvSFZK2qFOGlSR9XdIDkuZKekPSnZKOKOQbC0zKP56Ty1J5Dc95lpd0sqT7JL0q6c38Pv4kaZ8O+dDMmli/RhfArMH+FzgJeA64CHgbOAjYBVgOeKfGNYcB+wJPAQuBayoncjC9Jl97PTCDNJfz48D+kkZExH1V+VcFbgN2AO4DfkX6MvsR4ApJ20TEN3L2P+bjZ4ApwOSqMj2Vj2OBI4CHgcuAecAgYFgu862lPhUzq0le+cualaQ9gL8AjwG7RsTsnL48KbjsQRoZu3FOHw1cCgSwH/BlYK2I2DGfXw14ghRIPxQR06uetQ1wN2mj6h2r0seSguDpEfGDqvT+pCD5YWDHiHggpw8n1TbPjYgxhfezCvAqKfjuGhELC+fXiIiX2/RhmRng5llrbp/Jx+9UAiZARLwNfL2F6/4UETex9PJ5nyatDnROdcDM93wEuBjYQdLWkIIYcBRwb3XAzPnnA6cDAo4s+X4i538LWLTUSQdMs3Zz86w1sx3ycWqNc3cBC+pcd4+kNYF1WXK6yQfzcfs6U0I2z8etgOnAzkBfoN4UkuWq8rcqIl6XdD1wIPCApKuB24G7I+LNMvcws5Y5aFozWyUfZxVPRMRCSfVqZi9Qe/m8NfLxhFaeu3Ih/8751Vr+MkaRaqhHAufmtPmSrgJOjYil3quZlefmWWtmr+fj2sUTkvqyOKgVBbU3nq6sCrR9RKiF168L+X/cSv4RZd9QRMyLiDERsTmwIan5d2o+XlX2PmZWm4OmNbP783FYjXO70XJLzBDgxULN7a583KPk8+8h9T2WzQ9pkBGkZt0WRcQzEXE5aSTuv4BhuR/VzNrIQdOa2WX5eFYeeQq8O3r2u61cW2v5vEuB2aQ5lLsUL5DUpzKfEiAiXgQuB3aSdLakpYK0pMGS3leVVGky3rBG3vdK2rVGWQcAA0l9tG+3+K7MrEXu07SmFRFTJF0EfBZ4JA+ceYc0kOY10tzNpUahkkaobgP8onC/lyUdBlwL3CVpIvBIvseGpIFCawD9qy77ArAZcB5wtKSppD7WQaQBQDuT5l0+mfP/E5gJfFLS28C/Sc3FvwFWy8/9B2nayTPAe4ADgHWACyJizrJ/UmZW4aBpze4k4FHgc8CJpJrctcCZwLPA4zWuWQtYkRrL50XERElDgFNJzaJ7kGp3z5EWMbi6kP91SXuSAveRwKGkoDqL1KR6CjChKv9CSYcA3wcOJ9UgReq3fAA4BxgOjADWBF4hBdozKKxcZGbLzosbmNUgaTPSogdXRkRxObvDgD8AO0XE3xpRPjNrDPdpWlOTtI6kPoW0lYDz84/X1rhsO1KT6/Qa58ysF3PzrDW7LwNHSJoMPE/q+9ubtF7sjaQaJQCSBpLmQR4BvIj//5g1Hf+nt2Y3AdietMbr6qQRpo8BFwDnR+6/kDQM+o+HEX3grgGwwgKYPVPSfhFRa0UhM+uF3Kdp1opUw+w/E64bmKZvvgf4FrArcNAcmD8oIt5obCnNrCu4T9OsdaNSDXMkaQYJpGmaI0npjGpYycysSzlomrWq36aw+4D058osk8oqekMHQN/BjSiVmXU9B02zVi2YAXfMTX8+jjTtcaN8btpcWFhrLqeZ9ULu0zRrxZJ9miOrzkzAfZpmzcVB06yEJUfPDh2QapiTFsF8j541ayIOmmYlSVoZGJX6MBc+DoxzDdOsuThompmZleSBQGZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlaSg6aZmVlJDppmZmYlOWiamZmV5KBpZmZWkoOmmZlZSQ6aZmZmJTlompmZleSgaWZmVpKDppmZWUkOmmZmZiU5aJqZmZXkoGlmZlbS/x8Ah1Wqe5XVONEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "C2= ot.dist(s1_embed,s2_embed)\n",
    "G=ot.emd(ot.unif(4),ot.unif(4),C2)\n",
    "\n",
    "pl.figure(figsize=(6,6))\n",
    "pl.scatter(npos[:4,0],npos[:4,1],c='r',s=50, edgecolor = 'k')\n",
    "for i, txt in enumerate(s1):\n",
    "    pl.annotate(txt, (npos[i,0]-4,npos[i,1]+2),fontsize=20)\n",
    "pl.scatter(npos[4:,0],npos[4:,1],c='b',s=50, edgecolor = 'k')\n",
    "for i, txt in enumerate(s2):\n",
    "    pl.annotate(txt, (npos[i+4,0]-4,npos[i+4,1]+2),fontsize=20)\n",
    "for i in range(G.shape[0]):\n",
    "    for j in range(G.shape[1]):\n",
    "        if G[i,j]>1e-5:\n",
    "            pl.plot([npos[i,0],npos[j+4,0]],[npos[i,1],npos[j+4,1]],'k',alpha=G[i,j]/np.max(G))\n",
    "pl.title('Word embedding and coupling with OT')\n",
    "pl.axis('off')\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sentence similarity\n",
    "We will now explore the superiority of this Word mover distance (WMD) in a regression context, where our goal is to estimate the similarity (or relatedness) of two sentences on a scale of 0 to 5 (5 being the most similar). Given a set of pairs of sentences with a human annotated relatedness, our goal is predict the relatedness from a new pair of sentences.\n",
    "\n",
    "We will use the [SICK (Sentences Involving Compositional Knowledge) dataset](http://clic.cimec.unitn.it/composes/sick.html) for this purpose.\n",
    "\n",
    "We first load it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A group of kids is playing in a yard and an old man is standing in the background\n",
      "A group of boys in a yard is playing and a man is standing in the background\n",
      "4.5\n"
     ]
    }
   ],
   "source": [
    "   \n",
    "data=np.load('data/data_text.npz')   \n",
    "setA=data['setA']\n",
    "setB=data['setB']\n",
    "scores=data['scores']\n",
    "\n",
    "print (setA[0])\n",
    "print (setB[0])\n",
    "print(scores[0])\n",
    "\n",
    "np.savez('data/data_text.npz',setA=setA,setB=setB,scores=scores)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will only keep 200 sentences for learning our regression model and the rest for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=200\n",
    "testA=setA[n:]\n",
    "trainA=setA[:n]\n",
    "testB=setB[n:]\n",
    "trainB=setB[:n]\n",
    "\n",
    "scores_train=scores[:n]\n",
    "scores_test=scores[n:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the countVectorizer model from ScikitLearn, compute all the bag-of-words representations of the sentences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "vect = # TO BE FILLED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a big data matrix of all the words present in the dataset embeddings\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_feat = # TO BE FILLED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute a big matrix of all pairwise feature distances using the dist() method of POT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = ot.dist(all_feat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now you can write a code that will compute the Cosine and WMD dissimilarities from all the pairs of the training set "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_cos=[]\n",
    "X_wmd=[]\n",
    "Y=[]\n",
    "\n",
    "\n",
    "\n",
    "for i in range(len(trainA)):\n",
    "    s1 = vect.transform([trainA[i]]).toarray().ravel()\n",
    "    s2 = vect.transform([trainB[i]]).toarray().ravel()\n",
    "    # Cosine similarity between bag of words\n",
    "    d_cos=# TO BE FILLED\n",
    "    X_cos.append(d_cos)\n",
    "    # WMD\n",
    "    d_wmd=# TO BE FILLED\n",
    "    X_wmd.append(d_wmd)\n",
    "    Y.append(scores_train[i])\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the corresponding golden similarities / distance from the learning set. Hence you have a first appreciation of how much WMD better captures this similarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+cFeV96PHPlwWEKBUJq5Hlp4Ruq6UB3QqW3NbQmDXE6t6oLVYbbZPYJLZNruk20PCyiSWVlFuT21fSpmrTaLAm9Ue2JJBSUuQm4SLpIiL+KFFUwBUFRQhGVNz93j9mDs6eM2fPj505zzNzvu/X67z2nGfmzHzPzJznzM73mecRVcUYY0y+jHAdgDHGmORZ5W6MMTlklbsxxuSQVe7GGJNDVrkbY0wOWeVujDE5ZJV7zonIoyJyfoPXOVVEXhGRljrf/4qInBE+/4aILB9GLN8Xkavrfb+PRGS6iKiIjCwz/XMisqrRcRm/WOXuERH5PRHpDSu3fWHF9O7hLFNVz1LVjQmFeJyITBaRe0XkRRE5LCKPiMg14Tr3qOpJqtpfz7LD9z6VRJyq+n5VvT2M+RoR+XE9yxGRJSLyw5jyiSLyhoj8ioiMFpG/FZFnw334jIh8ebifwZh6WOXuCRG5Hvgy8NfAacBU4O+BS1zGNYRvAnuBacDbgd8HXnAaUYQEkjy+VwG/LiIzisoXAztU9RFgKdABnAuMA84HHkwwhlwr95+IqZOq2sPxAzgZeAW4fIh5TiCo/J8LH18GTginTQS+BxwCDgI/AkaE054B3hs+/xzwr8AdwBHgUaAjso5JwL3AAeBp4E+HiOcVYE6ZadMBBUaGrzcCy4H/F77vuwQ/CHcCPwP+C5geeb8C7wyffwNYHj4/JfycB4CXw+eTI+/bCHwB2AQcBd4Zln0E+GXgNaA/jOEQ8GsEP0gtkWV8ENhe5nP9B3BDUdlPgE+Gz78HfKqG/f4+YCdwmOCH/P8CHwmnjQCWAbuB/eE+O7nM9p0RvvcIsB74CrAqsp754bY/BGwHzi/aZn8VbrMj4WecWCbeoY6zKcB94b55CfhKDZ/jw8Ae4IeV4rVHDfWK6wDsoQAXAm8Wvqxl5rkReAA4FWgND/6/CqfdBHwNGBU+/gcg4bRnGFy5vwYsAlrC9z0QThsBbAVuAEYDZwBPAZ1l4vlBWCEsBqYWTSuufDYCTwIzCX7IHgN+CrwXGBl+4f858v5ylfvbgUuBtxGcGd8N9ETetzGsJM4KlzsqLCtUmNcAPy6K9THg/ZHX3wE+XeYzXwk8EXndDrwBtIavl4Xr/wQwu7APyixrIsEP2wfDWD8JHIvE+ofhNjsDOImg4vxmme27GbiZ4ATgNwgq6VXhtDaCynZRuI8vCF+3RrbZLuAXgbHh6xVlYo49zgiOpe3Al4ATgTHAu2v4HHeE7xtbKV571FCvuA7AHscrjecrzLMLWBR53Qk8Ez6/Efg3wgqx6H3PMLhy/0Fk2pnA0fD5PGBP0XuXEql0i6adAqwgOPvvBx4Cfi2cVlz5bAQ+G3nv3wLfj7z+beChyOvYyj0mhjnAy5HXG4Ebi+bZyNCV+2eAO8PnE4BXgdPLrO9tBBXyr4evvwD8W2R6C3AdwY/e6wT/YV1dZlkfAjZHXgvBZa5CrP8JfCIyvZ2g8h8Z3b4El+/eBE6MzPsvvFW5f4awMo1MX1eIK9w+yyLTPgH8e5mYY48z4DyCM/aSk5MqP8cZRfujbLz2qP5h19z98BIwscI1x0kE/9oW7A7LAFYSnB39h4g8JSJLhljO85HnrwJjwvVOAyaJyKHCA/gLguv/JVT1ZVVdoqpnhfM8BPSIiJRZb/R6/NGY1ycNETMAIvI2EflHEdktIj8DfgiML2qVs7fScoqsAn5bRE4Efgf4karui5tRVV8l+G/hQ+HnvJLgrLMwvV9Vv6qqC4DxBJX/10Xkl2MWNykaqwa12LNF04v390hK98ckgh+4nxfNWzANuLxov74bOD0yT/ExUW5flDvOpgC7VfXNmPdU8zmi+6yaeE0VrHL3w2aCM72uIeZ5juDAL5galqGqR1T106p6BnAxcL2I/FaNMewFnlbV8ZHHOFVdVOmNqvoi8L8JvsgTalxvLT5NcOY3T1V/geASBARnvcfDGeL9JdNUtY9g+3+QICn8zQox3E7wI3ABwaWh78auSPWoqn6VIDdwZsws+4DJhRfhj8XkyPS4/f0mpUnrfcAp4Y9TdN6CvQRnwtH9eqKqrhjiM8Ya4jjbC0wtc3JSzeeI7pfE4m12Vrl7QFUPE1zr/qqIdIVnqKNE5P0i8jfhbHcBy0SkVUQmhvOvAhCRi0TknWEFcZjgMslAjWH8BDgiIp8RkbEi0hI27/u1uJlF5Ivh9JEiMg74OPCkqr5U8wao3jiCs/xDIjIB+Msa3/8CMFlERheV3wH8OcF18vsqLONHBIm+W4BvqeobhQki8ikROT/cfiPD9vXjgG0xy1kDzA7390iCyznviEy/C/hfIjJDRE4iaEX17eKzY1XdDfQCnw+bYr6b4DJXQeE/k85wn44JY4z+kFRliOPsJwQ/MitE5MRwHQtq+RxpxNvsrHL3hKr+LXA9QVLuAMEZzB8DPeEsywm+xA8DOwia2BVu7plFkOB8heAs9O9V9f4a198PXERwHftp4EXgNoIEaJy3ESQfDxEkXqcRnM2l6csESbcXCZLL/17j+zcQ5AieF5EXI+XfIYj/O+Gll7LCyyd3hPPfUTT5VYJ8wvNhjNcBl2pMm/3wv53Lgb8huCx3JsH+fT2c5esE/0X8kGB/vAb8SZmwfo8gZ3KQ4AcveqloL0Fz2r/greOqm/q++7HHWXjs/DZB66Q9BJeXfreOz5F0vE2t0KLCmKYmIruAP1LVHzha/wiCSvHKWn+YjYljv4am6YnIpQTXfTc0eL2dIjJeRE4gOFMVgv9IjBk2uyPMNDUR2UhwSeT3VbXWPMVwnUfQbHE0QXv7LlU92uAYTE7ZZRljjMkhuyxjjDE55OyyzMSJE3X69OmuVm+MMZm0devWF1W1tdJ8zir36dOn09vb62r1xhiTSSKyu/JcdlnGGGNyySp3Y4zJIavcjTEmh6xyN8aYHLLK3Rhjcsgqd2OMyaGqmkKKyDMEQ3f1A2+qakfRdAH+D8HQWK8C16hqagMD92zrY+W6nTx36CiTxo+lu7Odrrltaa3OGGMyp5Z27u8JuymN836C7kBnEXQ9+g/h38T1bOtj6X07OHqsH4C+Q0dZet8OAKvgjTEmlNRlmUuAOzTwAMHQZ6kMi7Vy3c7jFXvB0WP9rFy3M43VGWNMJlVbuSvBuIlbReTamOltDB4H8dmwbBARuVZEekWk98CBA7VHCzx3KL7TvHLlxhjTjKqt3N+tqmcTXH65TkR+o9Ib4qjqLaraoaodra0Vu0aINWn82JrKjTGmGVV1zT0cRBhV3S8i3wHOJRg2q6CPYAT0gslhWeK6O9vpvmc7x/rf6qp4VIvQ3dmexupSZ8lhUw87bkwlFc/cwwFvxxWeA+8DHimabTXwIQnMBw6r6r7Eoy0o7oI+o13SF5LDfYeOoryVHO7ZlsrvoskJO25MNaq5LHMa8GMR2U4wyvkaVf13EfmYiHwsnGctwSDJTwK3Ap9IJVqChOqxgcG1+bEBzWRC1ZLDph523JhqVLwsE47c/q6Y8q9FnivBSO+py1NCNU+fxTSOHTemGpm7QzVPCdU8fRbTOHbcmGpkrnLv7mxn7KiWQWVjR7VkMqGap89iGseOG1MNZyMx1atrbhu9uw9y15a99KvSIsKl57QN2VLA15YFhRh8jM34y44bfy3r2TGobrpi3hSWd812EosEl8sbr6OjQ+sZZq+4+wEIzlpu+uDs2IO71vmNMaYey3p2sOqBPSXlV82fmmgFLyJbi/v3ipO5yzK1thSwlgXGmEa4a8vemsrTlrnKvdaWAtaywBjTCP1lroKUK09b5ir3WlsKWMsCY0wjtIjUVJ62zFXu3Z3tjGoZvLGG6n7A95YFy3p2MHPpWqYvWcPMpWtZ1rPDaTw92/pYsGIDM5asYcGKDXbXo6lZsx5DV8ybUlN52jLXWgaoqfsBn1sWFCdg+lWPv3aRYbe+8s1wNfMxVPjOWmuZOlvLLFixgb6Y6+Vt48eyacnCJEJrmJlL18Zej2sRYddNixoeT562rXHDjqH05ba1TJ4SpL4lYPK0bY0bdgz5I3OVe54SpL4lYPK0bY0bdgz5I3OVu+8J0lr4loCpNVlt3PE1aZmn72fWZS6h6nOCNBdy0ld+nvmctLTvpz8yl1DNE0uomnrYfmpuuU2o5oklVE09bD+Zaljl7pAlVE09bD+ZalRduYtIi4hsE5HvxUy7RkQOiMhD4eMjyYaZT5ZQNfXo7mynZcTg/dQywvaTD3xKdNdy5v5J4PEhpn9bVeeEj9uGGVdTWN41m6vmTz1+pt4iknj3oDWzhKr3encfpL9oHOH+AaV390FHERnwb+Dyqip3EZkMfACwSjthy7tms+umRTyz4gPsummR04o9T4OP55lvXcuagG/di1d75v5l4M+BgSHmuVREHhaRe0Qk9rqCiFwrIr0i0nvgwIFaYzUps0RdNviWiDcB374/FSt3EbkI2K+qW4eY7bvAdFX9VWA9cHvcTKp6i6p2qGpHa2trXQGb9FiiLht8S8SbgG/fn2rO3BcAF4vIM8C3gIUisio6g6q+pKqvhy9vA85JNErTEHZ3YTb4log3Ad++PxXvUFXVpcBSABE5H/gzVb0qOo+InK6q+8KXFzN04nXYLrh5I0/s//nx17NOPZH115+f5ipT49Pg3fUMPm4az7euZU3At7tza7pDNVK5XyQiNwK9qrpaRG4iqNTfBA4CH1fV/x5qWfXeoVpcsRdksYL3bfBu3+IxxpRK5Q5VVd2oqheFz29Q1dXh86WqepaqvktV31OpYh+OuIp9qHKf+ZZd9y0eY0z97A5Vh3zLrvsWjzGmfla5O+Rbdt23eIwx9ctc5T7r1BNrKvdZd2c7o4puIx/l8DbyoQYZN37x6Tb3Yq5jc71+X2Sucj913Ak1lXuvuGmyw6bKd/fuqancuOHbbe4+xeZ6/T7JXOW+aVd8/xnlyn22ct1OjvUX3e7f7+52/zxt2zzzOfHtOjbX6/dJ5ir3PLEEpqmHz8eN69hcr98nVrk7ZAlMUw+fjxvXsblev08yV7kvmDmhpnKf+Xa7cp62bZ753O++69hcr98nmavcZ7SeVFO5z7rmtnHTB2fTNn4sQjAGpsu7Qe/86HklFfmCmRO486PnOYnHDMHnfvddx+Z6/Z7I3ADZvg0qbUyj+TxAtuvYXK+/EXI7QLb1ZW2anc9JQ9exuV6/TzJXuVtf1qbZ+Zw0dB2b6/X7JHOVu/VlbZqdb4n4KNexuV6/Tyr25+6b5V2z2fLUSyX9uSfVl3Wj+1e/8tbNg24SsgSmqcS3fsOjXMfmev0+yVxCdVnPDlY9UHo7/FXzpw67gm90f+bFFXuBVfDGmHJym1BNc+T3Rt+6bLf7G2PSkrnKPc3WMpZpN8bkReYq9zRby1im3RiTF1UnVEWkBegF+gpD7UWmnQDcAZwDvAT8rqo+k2Ccx10xb0rsNfckWst0d7bHXnNPK9O+YOaEstfcXVnWs8MGXjbDYo0ESrn4XtVy5v5J4PEy0z4MvKyq7wS+BHxxuIGV0zFtQknQI8Ly4Wp0dwC+3e5fSFYXLnH1q7LqgT0s69nhJB6TPXGNBDbtOsiVt252FJF7rr5XVbWWEZHJwO3AF4DrY87c1wGfU9XNIjISeB5o1SEWXm9rmWa4vdgV69rBDNf0JWvKTntmxQcaGIk/kv5eJd1a5svAnwMDZaa3AXsBVPVN4DDw9pigrhWRXhHpPXDgQJWrHsySnumxrh2MSZ6r71XFyl1ELgL2q+rW4a5MVW9R1Q5V7Whtba1rGZb0TI917WBM8lx9r6o5c18AXCwizwDfAhaKyKqiefqAKQDhZZmTCRKrievubKelaFDpFoeDSueJde2QHVfeupnpS9Ycf/hyTdvGBCjl6ntVsXJX1aWqOllVpwOLgQ2qelXRbKuBq8Pnl4XzpPI/R+/ug/QPDF50/4DSu9tu/Bmu5V2zuWr+1ONnFC0iidz5a5Llc9LSt0YCPkizEchQaup+QETOB/5MVS8SkRuBXlVdLSJjgG8Cc4GDwGJVfWqoZVl/7sbUx5KW2ZJ0I5BqE6o1dRymqhuBjeHzGyLlrwGX1xZifSzpZ4zJEleNQOwOVWOMSZGrRiCZq9wt6WeanSUts8VVH/OZ7M99/aPP88KRN46XnTZutCX9EnLBzRtL+spff/357gIyJe786Hmx+8mXpKXrLiwaPSZDtetudEyZq9yvvHXzoIod4IUjb3DlrZu9ObizqrjCAHhi/8+54OaNVsF7pGdbH8++/Nqgsmdffo2ebX3OB6UoHm+hcKs90JAKvnhMhr5DR1l6X3Cbf9rbptK6G71vMndZxvpAT09xxV6p3LjR6HEHapHmeAvVcLltfNsvmavcjWl2PnfB4bo1m8tt49t+scrdmIzxuQsO163ZXG4b3/ZL5ir3vLUUWNazg5lL1zJ9yRpmLl3rtHvdWaeeWFO5caO7s52iHjgYIXjRBYfr1mzdne2MKto4oxrUPclQrWJcdBeRucr98o6pNZX7zLf+0697z6yayo0bd/fuoagHDgY0KHfN1a32gxT/k9CgW2DKjQdxd+8eJ91FZK5yL5ec8CGZVCvXyadiedq2eeZzo4KV63aW9As+QOOOoZXrdnKsf/Av37F+bdj6u+a2sWnJQp5e8QE2LVkYvHa0vzJXufuWtBgO18mnYnnatsYN18eQ6/X7JHOVu29Ji+Eovm5aqTxtedq2xg3Xx5Dr9fskc5W7y4RJ0k4YGb/5y5WnrbuzPfZ6aRa3bZ753KjA9THk6lb/objaX5mr3AFnCZOkHT0WP2phufK09e4+GHu91PrK98uM1pNqKm8k18dQowe5r4arPu5r6s89STZAtn990/sWj4nn837yOba8SHqAbG/kKWHiW0LVt3hMPJ/3k8+xNZvMVe55Spi0lYm5XHnaXN9daKrj837yObZmU7FyF5ExIvITEdkuIo+KyOdj5rlGRA6IyEPh4yPphOtnwqRevn0W13cXmur4vJ98jq3ZVHPm/jqwUFXfBcwBLhSR+THzfVtV54SP2xKNMqJrbhuTTxkzqGzyKWOcd3VaD9+SP8u7ZnPauNGDyqyvfP/4PJD58q7ZsclDH2JzqWdbHwtWbGDGkjUsWLGBnm19qa+z1gGy3wb8GPi4qm6JlF8DdKjqH1e7rHoTqnEjv4ONsJ4E27ZmuIr7NIfgv1HXLVZcSnqbJJpQFZEWEXkI2A+sj1bsEZeKyMMico+IpPY/mM+3XmedbVszXL71ae4DV9ukqspdVftVdQ4wGThXRH6laJbvAtNV9VeB9cDtccsRkWtFpFdEeg8cODCcuI0xHspTa7akuNomNbWWUdVDwP3AhUXlL6nq6+HL24Bzyrz/FlXtUNWO1tbWeuI1xngsT63ZkuJqm1TTWqZVRMaHz8cCFwD/XTTP6ZGXFwOPJxlklM+3XtfDp/7c87Zt88yn4ybqPb8Uf9JWrrwZuGoVV82Z++nA/SLyMPBfBNfcvyciN4rIxeE8fxo2k9wO/ClwTTrhWn/uafL5tnbzFt+Om6g1D++rqbwZuGoVZ90POOTbrdq+xWPi+byfpi9ZU3baMys+0MBI8su6H8gA327V9i0eE8/2k6lG5ir3PCVsfLtV27d4TDyf99P4saNqKjfpyVzl7tst+8Ph263avsVj4vm8nz538Vk1lafBxd2gldbtIqaRqa8hYYUkxMp1O3nu0FEmjR9Ld2d7Ju9+K9ySfdeWvfSr0iLCFfOmNP2t2mZoPh835Qbpvrt3T0O+o8V3g/YdOsrS+4JEc9rrL7fu3t0HuXdrX8NjylxC1aTH50SdyQbXCVWXDS7KrbtFJPZ7VW9MuU2omvRYos5kncsGF+XWUe7749UdqibffE7UGVMNlw0uyq2j3Pcn7Zgyd80dYN4X1vPCkTeOvz5t3Gi2fPYChxHVr2dbnzf5gyvmTWHVA6XXTH1I1JnBlvXs8PKa+4KZE8r2LNoI3Z3tsT0wNqLBRbl1X3pOG/+yZQ8DkRP4EZL+oOGZO3MvrtgBXjjyBvO+sN5RRPUrJGD6Dh1FeSvR0sjsvsken+9QdX0HucsxEsqt++kDrwyq2AEGtHzyOSmZS6i6Ttgkybe7bS2hmg0+7yffjmkfJF1nWUI1A3y729YSqtng837y7ZhuZla5O+Tb3baWUM0Gn/eTb8d0M8tc5V48xmelcp/5dretz3c+mrf4vJ98O6Z94Kor7cxV7ls+ewFjWgafoYxpkUy2lvFxgGwb3Nh/Pg+Q3TW3jUvPaRsU26XntGXyDvJ6xHUzcOdHz4v9XqU9LnHmEqqFlgLFfDm4s8wGNzbD1czHUKM+e24Tqndt2VtTuameDW5shquZjyHfPnvmKnefWwpknbV0MMPVzMeQb589c5W7zy0Fss5aOpjhauZjyLfPXs0A2WNE5Ccisj0cJ/XzMfOcICLfFpEnRWSLiExPI1jwu6VAPVz2PV2su7OdUSMG/0iOGiFN3dLB1KaZB8geqqXQBTdvZPqSNccfF9y8MfV4qjlzfx1YqKrvAuYAF4rI/KJ5Pgy8rKrvBL4EfDHZMN/SMW1CSdAjwvKs8bL7geJ/gOwfIlOD+7Y+W1N5npRr/fbV+5/gif0/HzTvE/t/nnoFX7Fy18Ar4ctR4aP4AvclwO3h83uA3xJJ5zrJynU7GSgqGwjLs8a3BMzKdTs51j941x7r10xuW+PGq8eKv51Dl+dN19w2Ni1ZyNMrPsCmJQvpmttWUrEXlCtPSlXX3EWkRUQeAvYD61V1S9EsbcBeAFV9EzgMvD1mOdeKSK+I9B44cKCugH1LWgyHb5/Ft3iMMfWrqnJX1X5VnQNMBs4VkV+pZ2WqeouqdqhqR2trfdfgfEtaDIdvn8W3eEz2lPt/3do7NF5NrWVU9RBwP3Bh0aQ+YAqAiIwETgZeSiLAYnm6vdm3z9Ld2c6oort/R7VYQtVHy3p2MHPpWqYvWcPMpWu96O4X4Mp58V37livPm7gGErNOPTF23nLlSammtUyriIwPn48FLgD+u2i21cDV4fPLgA2a0q2vvt2yPxxefpbivWa3D3jH5/7cfe4aIW3lGkicOu6E2PnnnVFy5TpRFbsfEJFfJUiWthD8GPyrqt4oIjcCvaq6WkTGAN8E5gIHgcWq+tRQy7UBsv1jfXFng8/9uTezct+fcurdX9V2P1BxmD1VfZig0i4uvyHy/DXg8lqDNH6xhGo22F3afqr1e5L2/srcHaomPZZQzQa7S9tPtX5P0t5fVrmb43xL8Jp4ebtLOy/KfX/K9due9v6yyt0c1+x9cWdFMyctfVaugcSdHz2vpGXMrFNPTH1/Za4/d5OeZu6L25i0JD0GRW77czfp8a07BGPywNUYFFa5m+OstYwxyXPVuskqd3OctZYxJnmuWjdVbOdu0tWzrY+V63by3KGjTBo/lu7OdmfXt7s72+m+ezvHBt46o7D+3E2trrx1M5t2HTz+uhGDQfvsinlTYq+5W2uZHLP+3E3eFFfsAJt2HeTKWzc7isi9jmkTKBoDhxGS/hgUVrk75FsC0/pzN8NVXLFXKm8GK9ftZKDo8vqApj8GhVXuDvmWwPQtHmPywNX3yip3h3xLYPoWjzF54Op7ZZW7Q77d7t/d2R57bdASqv7xaWD1qHK32pcrbwauxkmwyt0h3/pzv7t3T+y1wbt7SzP9xh0vE/GhGa0n1VTeNByMk2BNIR3rmutP3y2WDMuGoRLxro+loe7GbNa+b1au2zmoeTHAsQFNfX/ZmbsxGeNz4tv6mi9lCVVjTFV8TnxbX/OlvE2oisgUEblfRB4TkUdF5JMx85wvIodF5KHwcUPcsozfLBmWDb4l4qOsr/lSrvZXNWfubwKfVtUzgfnAdSJyZsx8P1LVOeHjxkSjNA1x50fP47RxoweVnTZudFPfOu4jn/vdX941u+RkYMHMCU17vR2C/XX21JMHlZ099eTU91fFyl1V96nqg+HzI8DjgPujyCRuWc8OXjjyxqCyF468wbKeHY4iMnF6tvVx79a+49ex+1W5d2ufF61lerb18eCew4PKHtxz2IvYXFnWsyO2S4a0v1c1XXMXkekEg2VviZl8nohsF5Hvi8hZCcRmGsxVv9OmNr51WxHlc2yuuPpeVd0UUkROAu4FPqWqPyua/CAwTVVfEZFFQA8wK2YZ1wLXAkydOrXuoE06rKVDNvjcWsbn2Fzxuj93ERlFULHfqar3FU9X1Z+p6ivh87XAKBGZGDPfLaraoaodra2twwzdJM1aOmSDz61lfI7NFVffq2paywjwT8DjqnpzmXneEc6HiJwbLvelJAPNK59uI7eWDtngc2sZ68KilKvvVTVn7guA3wcWRpo6LhKRj4nIx8J5LgMeEZHtwN8Bi9XVyNsZ4ttt5Mu7ZnPV/KmDWmHUO4ivSY9v3VZE9e4+GNuFRe/u5r3LuWPaBFqKfvFaRkjq/bmLqzq4o6NDe3t7nazbFwtWbKAv5lpk2/ixbFqy0EFExgzPzKVrY68lt4iw66ZFDiJyL+nvuYhsVdWOSvPZHaoOWfLJ5I0l5UtZ9wNNyJJPJm8sKV/K2+4HTHp8TowZUw9Lypfq7myPvebuQ/cDJiU+J8aMqYer5KHPencfpL8oy9w/oKknmS2haoxJjDUSKJV0ktkSqsaYhrNGAqW8vkPVGGOqYY0ESnl7h6oxxlTL1WDQPnOVZLYxVI0xyXIwGLTPCnd437VlL/2qtIhwxbwpqd/5bQlVY0xiLKGaPkuoGmMazhKq/rDK3RiTGEuo+sMqd2NMYuyua39YQtUYk5iuuW307j44KHnoy+DdLvVs62Plup08d+gok8aPpbuz3f0A2cYYUy2fB+92xdW4DVa5G2MSYwNkl3K1TaxyN8YkxlrLlLL+3I0xmWetZUp525+7iEwRkftF5DEReVREPhkzj4jI34kqp1rtAAAKZUlEQVTIkyLysIicnU64xhifdXe2M6qoy99RDei73GeuumSoprXMm8CnVfVBERkHbBWR9ar6WGSe9wOzwsc84B/Cv8aYZlPcH1bzDsL0FgddMlQ8c1fVfar6YPj8CPA4UNyG5xLgDg08AIwXkdMTj9YY47WV63ZyrH9wzXWsX5s+oXqsaLCOYwPpb5OarrmLyHRgLrClaFIbsDfy+llKfwAQkWtFpFdEeg8cOFBbpMYY71lCtZT3CVUROQm4F/iUqv6snpWp6i2q2qGqHa2trfUswhjjMUuolvI2oQogIqMIKvY7VfW+mFn6gGjnxJPDMmNME7HuB0p1d7aXVLQjwvI0VdNaRoB/Ah5X1ZvLzLYa+FDYamY+cFhV9yUYpzEmA2zQ91K9uw8yUFQ2EJanqWJ/7iLybuBHwI4wJoC/AKYCqOrXwh+ArwAXAq8Cf6CqQ3bWbv25G2OagasBsis2hVTVH1OhMZMGvxDXVR+eMcY0Bxsg2xhjcsgGyDbGmByyAbKNMSaHbIBsY4wxFdkA2cYY08SscjfGmByyyt0YY3LIEqrGmES5GAzalLLK3RiTmMJg0IUxQwuDQQNWwTeYXZYxxiTGBsj2h1XuxpjEWH/u/rDK3RiTGOvP3R9WuRtjEmMDZMfr2dbHghUbmLFkDQtWbKBnW/rDXVhC1RiTLBsgexBXSWY7czfGJMYGyC7lKslslbsxJjGWUC3l/QDZxhhTiSVUS3k7QLaIfF1E9ovII2Wmny8ih0XkofBxQ/JhGmOywAbILuVqm1Rz5v4NgrFRh/IjVZ0TPm4cfljGmCzqmtvG2VNPHlR29tSTm/ru1K65bVx6TtvxkZdaRLj0nLbUt0nFyl1VfwikO0y3MSYXlvXsYNOuwdXFpl0HWdazw1FE7vVs6+PerX3Hx0ztV+XerX2pN4dM6pr7eSKyXUS+LyJnJbRMY0zG3LVlb03lzcBVa5kk2rk/CExT1VdEZBHQA8yKm1FErgWuBZg6dWoCqzbG+KS/zMhu5cqbQWZby6jqz1T1lfD5WmCUiEwsM+8tqtqhqh2tra3DXbUxxjOF68rVljcDb1vLVCIi7xAJ9pyInBsu86XhLtcYkz1XzJtSU3kz6O5sp6WoS4aWBnTJUPGyjIjcBZwPTBSRZ4G/BEYBqOrXgMuAj4vIm8BRYLG6GnXbGOPU8q7ZQHCNvV+VFhGumDfleHkz6t19kP6BwVVi/4DSu/tgqi1mxFU93NHRob29vU7WbYwxjTJz6drYnEOLCLtuWlTz8kRkq6p2VJrP7lA1xpgUuUoyW+VujDEpcpVktsrdGJMoF32X+8xVktn6czfGJMYGyC7lKslsCVVjTGIWrNhAX8zNOW3jx7JpyUIHEeWPJVSNMQ1n/bn7wyp3Y0xirD93f1jlboxJjPXn7g9LqBpjElNImq5ct5PnDh1l0vixdHe2N20y1SWr3I0xieqam/5AFKYyuyxjjDE5ZJW7McbkkFXuxhiTQ1a5G2NMDlnlbowxOWSVuzHG5JCzvmVE5ACwu863TwReTDCcpFhctbG4amNx1cbHuJKIaZqqVhyE2lnlPhwi0ltNxzmNZnHVxuKqjcVVGx/jamRMdlnGGGNyyCp3Y4zJoaxW7re4DqAMi6s2FldtLK7a+BhXw2LK5DV3Y4wxQ8vqmbsxxpghWOVujDE55G3lLiKXi8ijIjIgImWbDonIhSKyU0SeFJElkfIZIrIlLP+2iIxOKK4JIrJeRJ4I/54SM897ROShyOM1EekKp31DRJ6OTJvTqLjC+foj614dKXe5veaIyOZwfz8sIr8bmZbY9ip3rESmnxB+9ifDbTE9Mm1pWL5TRDrrjaHOuK4XkcfCbfOfIjItMi12fzYormtE5EBk/R+JTLs63OdPiMjVDY7rS5GYfioihyLTUtleIvJ1EdkvIo+UmS4i8ndhzA+LyNmRaelsK1X18gH8MtAObAQ6yszTAuwCzgBGA9uBM8Np/wosDp9/Dfh4QnH9DbAkfL4E+GKF+ScAB4G3ha+/AVyWwvaqKi7glTLlzrYX8IvArPD5JGAfMD7J7TXUsRKZ5xPA18Lni4Fvh8/PDOc/AZgRLqcloe1TTVzviRw/Hy/ENdT+bFBc1wBfiXnvBOCp8O8p4fNTGhVX0fx/Any9AdvrN4CzgUfKTF8EfB8QYD6wJe1t5e2Zu6o+rqo7K8x2LvCkqj6lqm8A3wIuEREBFgL3hPPdDnQlFNol4fKqXe5lwPdV9dWE1l9OrXEd53p7qepPVfWJ8PlzwH6g4h14NYo9VoaI9R7gt8JtcwnwLVV9XVWfBp4Ml9eQuFT1/sjx8wAwOaF1DyuuIXQC61X1oKq+DKwHLnQU1xXAXQmtuyxV/SHBSVw5lwB3aOABYLyInE6K28rbyr1KbcDeyOtnw7K3A4dU9c2i8iScpqr7wufPA6dVmH8xpQfXF8J/zb4kIic0OK4xItIrIg8ULhXh0fYSkXMJzsh2RYqT2F7ljpXYecJtcZhg21Tz3nrVuuwPE5wBFsTtz0bGdWm4b+4RkSk1vjfNuAgvX80ANkSK09pelZSLO7Vt5XSYPRH5AfCOmEmfVdV/a3Q8BUPFFX2hqioiZduShr/Ms4F1keKlBJXcaII2r58BbmxgXNNUtU9EzgA2iMgOgkqsbglvr28CV6vqQFhc9/bKGxG5CugAfjNSXLI/VXVX/BIS913gLlV9XUT+iOC/noUNWnc1FgP3qGp/pMzl9moop5W7qr53mIvoA6ZEXk8Oy14i+LdnZHgGVigfdlwi8oKInK6q+8LKaP8Qi/od4Duqeiyy7MJZ7Osi8s/AnzUyLlXtC/8+JSIbgbnAvTjeXiLyC8Aagh/2ByLLrnt7FSl3rMTN86yIjAROJjiWqnlvvapatoi8l+DH8jdV9fVCeZn9mURlVTEuVX0p8vI2gvxK4b3nF713YwIxVRVXxGLgumhBiturknJxp7atsn5Z5r+AWRK09BhNsDNXa5CpuJ/gejfA1UBS/wmsDpdXzXJLrveFFVzhOncXEJtdTyMuETmlcFlDRCYCC4DHXG+vcN99h+Ca5D1F05LaXrHHyhCxXgZsCLfNamCxBK1pZgCzgJ/UGUfNcYnIXOAfgYtVdX+kPHZ/NjCu0yMvLwYeD5+vA94XxncK8D4G//eaalxhbL9EkKDcHClLc3tVshr4UNhqZj5wODxxSW9bJZUtTvoB/E+C60+vAy8A68LyScDayHyLgJ8S/Pp+NlJ+BsEX8EngbuCEhOJ6O/CfwBPAD4AJYXkHcFtkvukEv8ojit6/AdhBUEmtAk5qVFzAr4fr3h7+/bAP2wu4CjgGPBR5zEl6e8UdKwSXeC4On48JP/uT4bY4I/Lez4bv2wm8P+FjvVJcPwi/A4Vts7rS/mxQXDcBj4brvx/4pch7/zDcjk8Cf9DIuMLXnwNWFL0vte1FcBK3LzyOnyXIjXwM+Fg4XYCvhjHvINICMK1tZd0PGGNMDmX9sowxxpgYVrkbY0wOWeVujDE5ZJW7McbkkFXuxhiTQ1a5G2NMDlnlbowxOfT/ASd8eU+ELdi6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXu0XXV17z/fHI54IsgBOZfiSUIQGUnRKJEjj4ZaQocGEWkG6AWKDzq0ufXqrVaMI+m1GKxe8HJrba+23rRVQShSCaQ8i9iAPAroiQmPgFSQRziAHB4HCAkhj3n/WGsn++yzH2uvvdZejz0/Y5xx9l7r95i/x5577d+cv/mTmeE4juOUi2lZC+A4juMkjyt3x3GcEuLK3XEcp4S4cnccxykhrtwdx3FKiCt3x3GcEuLK3ekISbMkbZLUFzP/JklvCV9/X9JXO5Dlekkfj5s/j0iaLckk7dHg/gpJF3dbLif/uHLPIZKWS7q+5tqvGlw7PXxtkp6pVgKS+sNrVnXtZkmvSnpZ0kuS1kpaJmnPJvLMkLRK0rOSXpR0n6SzAMzscTPby8x2xGlrmPfXcfLWKev9ZnZhKPNZkm6LU07YH7fUub6/pNckvV3S6yT9laQnwi+oRyV9s9M2OE5SuHLPJ7cAv1N5GpZ0INAPzK+59tYwbYUXgPdXvX9/eK2Wz5jZ3sCBwNnA6cB1ktRAnh8AG4GDgDcBHwV+E69pyaOAJOfyxQT9f3DN9dOBe83sPmA5MAIcCewNHAf8IkEZSk2jXyJOcrhyzyc/J1Dmh4fvfxe4CXiw5trDZvZkVb4fAB+rev8x4KJGlZjZK2Z2M3AycAzwgQZJ3w18P0y/3czWmdn1MHXZIPxl8FVJ/xE+0V4t6U2SLgl/Kfxc0uxKwWHet9ZWKGlfSddIGpf0Qvh6RtX9myV9TdLtwGbgLeG1T0r6beA7wDGhDBOS3i3pN9XLR5JOkXR3nX55AlhD8CVWTXV/vhu40syetIBHzaxhX0t6n6QHw18+fyfpp5I+Gd6bJulLkh4Lf2ldJGmfBuUcHOZ9WdKNwP41948O+35C0t2Sjqvps7+UdHuY/8eS9q+tI0y7f9jnE5Kel3Rr5QtU0kxJV4Rj85ykb7VqR9U8+YSkx8P+bSqv0xmu3HOImb0G3AW8J7z0HuBW4Laaa7VLB6uB90galLQvwRfAv0ao73FgNExfjzuBb0s6XdKsCE04nUAxDgOHAHcA3wP2Ax4AvhyhjGlhnoOAWcAW4Fs1aT4KLCF4cn6sqj0PAH8C3BEu+wya2c+B54D31eRvpJAvpEq5S5pD8MX6z+GlO4HPS/rvkuY1+dVDqEAvJ3jafxPBl/TvVCU5K/xbCLwF2KtOWyv8M7CWQKn/JbDLxiBpGLgW+CpBX38BWCVpqCr/HwJ/BPwX4HVhmnqcDTwBDAEHAH8OWPjleA1Bf88mGOMfttGO3wN+G1gUUV4nJq7c88tP2a3If5dAud9ac+2nNXleBa4GTgv/rgqvReFJgg9YPT4c1v0XwCOS1kt6d5OyvmdmD5vZi8D1BL8wfmJm24EfAfNbCWNmz5nZKjPbbGYvA18jUAzVfN/MNoS/Jra1KpNAYX8EQNJ+wCJ2K+targQOkFRRwh8Drjez8fD9ecDXgTMJvhjH1NiYeyKwwcyuCPvgb4Gnq+6fCXzDzH5tZpsIvgROr126CL9Y3w38hZltNbNbCMa7wkeA68zsOjPbaWY3hrKdWJXme2b2n2a2BfgXdv8SrGUbwbLdQWa2zcxutSAQ1ZHAm4Gl4S+5V82sYtuI0o4VYb4tEeV1YuLKPb/cAhwbKqEhM/sV8B8Ea8H7AW9n6pM7BE+iH6PFkkwdhoHn690wsxfMbJmZvY3gKW49sLrJ02r1evyWOu/3aiWMpOmS/l/4E/8lgrYOarJXzsZW5dRwMfBBSW8A/itwq5k9VS+hmW0m+CL6WNjOM6nqTzPbYWbfNrMFwCDBl893wyWhWt5cLWuoJJ+ouf9Y1fvHgD0I+rq2nBfM7JWatBUOAj4cLnFMSJoAjiVQ0hWqv1Q203gsLgAeAn4s6deSloXXZwKPhV9StURpR/WYRZHXiYkr9/xyB7AP8MfA7QBm9hLBE/YfA0+a2SN18t1K8OE4gGAZpyWSZgJHhHmbYmbPAv+H4IPc6Ek/Cc4G5gBHmdkb2f2LpfoLpVlI0yn3zGyMoF9PIVhy+UELGS4k+BJ4L8HSz9X1EpnZFjP7NoHx+rA6SZ4Cqu0Fqn5PMKYHVb2fBWxnqtH6KWDf8MupOm2FjcAPwmWoyt8bzOz8Jm2si5m9bGZnm9lbCGwyn5f0+2Eds2p/VbTRjupxSUxeZyqu3HNK+LN1FPg8k5XubeG1ek/tlafCDwInW4t4zuHT8e8RrMv/DLiuQbqvK3D/20PS3sCngIfM7Lk2m9UOexM85U+Ev1SirNNX8xtghqTX1Vy/CPgiMA+4okUZtwITwErgh6EtBABJn5N0nKSBsF8+Hsq8rk451wLzJC0OleKngd+qun8p8GehsXQv4H8Bl9U+HZvZYwRz4lwFrpjHEox1hcovk0WS+iS9PpSx+oskEpJOkvTW8IvoRWAHsJNgnjwFnC/pDWEdC9ppRxryOlNx5Z5vfkpg+Kp+Ar81vFZXuQOE69AbmpT7LUkvEyjAbwKrgBPMbGeD9NMJ1qAngF8TPJ2dHLURMfkmMAA8S2C8/Lc2868BNgBPS3q26vqVBPJfGS69NCT8crwoTF+7xLUZ+CuCZY5nCRT2qfV89sNfOx8G/jeBUfcwAiW9NUzyXYJfEbcAjxDYSf5HA7H+EDiKYAnty0xeKtoI/AGB8XOc4Ml4KfE+54cCPwE2Efza+Tszuyncz/BBAjfcxwmWl06L0Y6k5XVqkB/W4fQakh4G/puZ/SSj+qcRKMUzzeymLGRwyo9/Qzo9haRTCdZ913S53kWhi+qeBE+qIvhF4jip4LvEnJ5B0s0ESyIfbbIElRbHELhdvg64H1gc2lUcJxV8WcZxHKeE+LKM4zhOCclsWWb//fe32bNnZ1W94zhOIVm7du2zZtYyRENmyn327NmMjo5mVb3jOE4hkfRY61S+LOM4jlNKXLk7juOUEFfujuM4JcSVu+M4Tglx5e44jlNCXLk7juOUkEiukJIeBV4mCPu53cxGau4L+BuCE1Q2A2eZmR8WnCCr141xwQ0P8uTEFt48OMDSRXNYPH84a7Ecx8kp7fi5LwxDl9bj/QQhQg8lCEf69+F/JwFWrxtj+RX3smXbDgDGJraw/Ip7AVzBO45Tl6SWZf4AuCg8Bf5OguPQ/KishLjghgd3KfYKW7bt4IIbHsxIIsdx8k5U5W4EZymulbSkzv1hJp+N+ER4bRKSlkgalTQ6Pj5ee9tpwJMT9YMHNrruOI4TVbkfa2bvIlh++bSk97TKUA8zW2lmI2Y2MjTUMjSCE/LmwYG2rjuO40Racw8PFsbMnpF0JXAkk495GyM4Fb3CjPBaLunEOJmFYXPpojmT1twBBvr7WLpoTqr1RsWNvY6TP1o+uYeH4O5deQ28D7ivJtlVwMcUcDTwopk9lbi0CVAxTo5NbMHYbZxcva71d1EneTth8fxhzjtlHsODAwgYHhzgvFPm5UKBZtUnjuM0J8qT+wHAlYG3I3sA/2xm/ybpTwDM7DvAdQRukA8RuEL+UTridk4z42QrZdlJ3k5ZPH84F8q8liz7xHGcxrRU7uFp7u+sc/07Va+N4PT33NOJcdINm1PxPnGcfNJzO1Q7MU66YXMq3ieOk096TrkvXTSHgf6+SdeiGic7yVtWvE8cJ59kdkD2yMiItXsSU1JeGdXl7DPQjwQTm7dNet2o/Lx5huRBnjzI4Di9gqS1tSFg6qYrinKv3YIPwRNiJ14j9cqsptPy0yaNPnEcJ99EVe6FWZZJYwt+vTKTLD9tPCyB4ziNKIxyT8MroxMPmTzgniqO4zSiMMo9Da+MTjxk8oB7qjiO04jCKPeli+bQP02TrvVPU0deGfU8PaoZ6O9j4dwhFpy/hoOXXcuC89c03Hm5et1YpHRJ0sxTJQt5nOTw8XM6pZ147tmjFu/bpGJ0rOc58+bBARbOHWLV2rGWcdSzirdeK3/FUwXw+O8FxuP3O0lQGG+ZBeevYazOWvLw4AC3Lzs+SdHarjML2ZqRN3mc9vDxc5pROm+ZLIyHUevMm2Ezb/I47eHj5yRBYZR7FsbDqHXmzbCZN3mc9vDxc5KgMMo9i23uUetMQrZqA9r8r/yYw8/9cWxjWj15+qeJza9tj1WmG/e6i4d0cJKgMAbVRsbDLAyWtXV2KlutAe2Fzdt23YtjTKtnKH7lte27ym2nTDfudZ8s5rpTPgpjUC0zjQxo1XRiTOvEQOfGPcfJF6UzqJaZtHfKegx7x+k9XLnngLR3ynoMe8fpPSIrd0l9ktZJuqbOvbMkjUtaH/59Mlkxi0UzA2S9e1F2yjYzprUyeJY5hr0bex2nPpHX3CV9HhgB3mhmJ9XcOwsYMbPPRK24rGvuzcLwAk3vVQxog9P7MYMXtzSOKx+lvtpdtJ0YfPNo3POQx04vkmg8d0kzgAuBrwGfd+XemGYGSCBx42QvGzx7ue1O75K0QfWbwBeBnU3SnCrpHkmXS5rZQKglkkYljY6Pj0esulg0M0B2M2xxLxg8e7ntjtOKlspd0knAM2a2tkmyq4HZZvYO4EaCp/wpmNlKMxsxs5GhoaFYAuedZgbIboYt7gWDZy+33XFaEeXJfQFwsqRHgR8Cx0u6uDqBmT1nZlvDt/8IHJGolAWimQEyDeNk3g2eadLLbXecVrTcoWpmy4HlAJKOA75gZh+pTiPpQDN7Knx7MvBAwnLmnmqj4+D0fvbcY1pDg2i1cXLh3CFWXLWBz122HoB9p/fz5Q++ra3dqKOPPc+ld21khxl9EqceMRwpf9KG0m4bXn0np+M0Jnb4AUlfAUbN7CrgTyWdDGwHngfOSka8YlAvfMBAfx9/fdrhdUMVVK6tXjfG0h/dzbadu43aL2zextLL796VNkrdq9aOsSM0jO8wY9XaMUYO2q9p/qTDCmQZ096VueNMpa1NTGZ2c8VTxszOCRU7ZrbczN5mZu80s4Vm9ss0hM0rcQ+qvuCGBycp9grbdljkQ647qTvJw7X9sG7HyRe+QzUB4nptNLsf1eMj6brjepq454rj5AtX7gkQ12uj2f2oHh9J1x3X08Q9VxwnXxQm5C/kd6fk0kVz6u6UbOW1sXTRnClr7gD9fbsP/m7V5lZ1N8rfLF/Ufq5ON/11U8MnxPVcyYOhtxtzLa/z2SkHhVHueY4rHtdro3J/xVUbmNgSxFqv9paJ0uZmdUfJH/dw7dqyX3lt8nq7ILLXTjV5MPR2Y67leT475aAw8dx7cat5p22Ok7/TQ8HjyNmpzEmX14251ovz2UmG0sVz70WDXadtjpO/00PB200Tt/40y+vGXOvF+ex0l8Io91402HXa5jj5Oz0UvN00cetPs7xuzLVenM9OdymMcs/zVvO0Yoo3ivP+ytbtkepo1GcL5w41lLdZP1e3c/Nr2+mfpoZ1xx2bpMc5annVbXtl63b6+9QyTyfkeT475aAwBtW8bjVP0zBWyX/u1RsmHZo9sWVbpDrq9dnCuUOsWjvWUN6ohtYXNm+jv08MDvTz4pZtbcWgj9LmpMY5Snm1YzixZRv908S+0/uZ2NxZezqRy3E6oTAG1bxSNONb3LLKbAAsc9uc8lE6g2peKZrxLS87WvNEmdvm9C6u3DukaMa3vOxozRNlbpvTu7hy75BuGMaSrCNuWWU2AJa5bU7vUhiDKuRzu3Y9o+eee7T+zqy0ZWxiC30SO8wYbtCmJI1vne6mzTosQBrluHHTKSOFMajm+aT7dmWrlz5KvjKR1HjmeV44ThqUzqCa53jh7cpWL32UfGUiqfHM87xwnCwpjHLPs0dDu7J1Eue9LCQ1nnmeF46TJYVR7nn2aGhXtk7ivJeFpMYzz/PCcbIksnKX1CdpnaRr6tzbU9Jlkh6SdJek2UkKCfn2aGhXtkZhBerlSyu0QdYsnDtEbfCCOOPZTt+36ssk+jqr8WpWb1nnkNOcdrxlPgs8ALyxzr1PAC+Y2VslnQ58HTgtAfl2kWePhnZlWzx/mNHHnueSOx+n2pxdGwO9rDG/K4d6N2t7VKL2fau+TKKvsxqvZvVCtPj8TvmI5C0jaQZwIfA14POVQ7Kr7t8ArDCzOyTtATwNDFmTwssSfiAuUba8l3VbfBbtalVnEjJlNV7N6gVKOYd6maS9Zb4JfBHY2eD+MLARwMy2Ay8Cb6oj1BJJo5JGx8fHI1ZdTqIYAstqLMyiXa3qTEKmrMarWb1lnUNOa1oqd0knAc+Y2dpOKzOzlWY2YmYjQ0NDnRZXaKIYAstqLMyiXa3qTEKmrMarWb1lnUNOa6I8uS8ATpb0KPBD4HhJF9ekGQNmAoTLMvsAzyUoJ5APw1BSMkQxBEaNrZ61Aa3d+rIwjreqMwmZkjTutkOzevPsiOCkS0uDqpktB5YDSDoO+IKZfaQm2VXAx4E7gA8Ba5qtt8chD8bFJGWIYgiMc4h1s3tp9FOcPsnCON6qziRkSsq4m3TbOm2XU0zaCj9QpdxPkvQVYNTMrpL0euAHwHzgeeB0M/t1s7KKeEB2HmRoJQd014CWlz4pEt5nTidENai2FTjMzG4Gbg5fn1N1/VXgw+2J2B55MAzlQYa4cqQlY176pEh4nzndwHeoFkyGVnJ0W8a89EmR8D5zukFhlHseDEN5kKGVHN2WMS99UiS8z5xuUJh47ovnD/Oj0ce5/eHnd11716x9umoYinrYclTjVW3ahXOHuOmX41Py1ivzvFPmZWpAq5ZpcHo/e+4xLdLh2HmMyd8OUeRvlKY6hn916IV9p/fz5Q++rVD94OSfwsRz/9Lqe7n4zsenXP/I0bP46uJ5SYoWm3ZiizeL6V6d99Qjhlm1dixX8crjxlAveuz1KPI3SlNvHBuV4TjNKF0890vv2tjW9SxoJ7Z4s5ju1XkvvWtj7uKVx42hXvTY61Hkb5Sm3jg2KsNxkqAwyn1Hg18Yja5nQTteEFE9Ixq1L0vPiqTj1BfFS6STkBGt5mlR+sApDoVR7n2qDRDb/HoWtOMFEdUzolH7svSsSDpOfVG8RDoJGdFqnhalD5ziUBjlfsZRM9u6XqGbW/EbeUEsnDs0RYZmMd2r855x1Ez6p01WDP3TlKlnRVxvj6WL5tDfV9OWvmzbUk2rudJJyIgzjpoZOYZ/VHmSaFM7fGn1vRyy/DpmL7uWQ5Zfx5dW39s6k5MZhVHuj4xvaus67DZujU1swdi9zTstBb94/jDnnTKP4cEBRLDjsGJIq5UBmJL2I0fPmvT+vFPmMXLQfkw51SLjHyv12hnZIFi7OpGTVbUocyVKuxul+eriebuuw+4n+UZ9l8TcTXL+VxwaKstLO8y4+M7HXcHnmMJ4y8xedm3De4+e/4G61/OwzbtTGfLQhqTIc1vyJlve4ssfsvy6unaDPomHzzuxrbKczkgl/EDRyIMBr1MZ4ubPoz95HsajEXmTLW/x5Yvg0OBMpjDLMnHIgwGvUxni5O/2clRU8jAejcibbHmLL18EhwZnMoVR7gsO2a+t69DZNu9uxm1POn9e/cnbbUvtGHxp9b2pGce7ERKgHYNkK3mizM8k2xTXocHJjsIsyxw8tNek0APV1xsRN0Z3t+O2J50/b0sMFdppS70xqN6hnHSc+rRjzNfusK4YJIG6O6ybyRN1fibZpoqMl961kR1m9EmccdTM3OwOd6ZSGINqNw06eTOutUvR5YfGbailKG1Kcv6WYXyd+JQu/EA3DTp5ffKNShmiDnZqcM4bSc7fos9PpzsURrl306CTN+Nau3Tkh54TOjU4540k52/R56fTHVoqd0mvl/QzSXdL2iDp3DppzpI0Lml9+PfJpAXtpkGnDE++i+cPc/uy43nk/A9w+7LjC6XYof4Y1FKkMUly/pZhfjrpE8WguhU43sw2SeoHbpN0vZndWZPuMjP7TPIiBnTToNOu4e/cqzfwwuZtAAwO9LPi5Naxubvlh55Hf/co1BuDRvHu2yFqDP0oedupP8n522p+FnXMnWRp94Ds6cBtwKfM7K6q62cBI+0o93YNqnlk9boxll5+N9t2TO7D/mnigg+/s6mS6EZc86LHT0+aqDH0o8bfz2NfFkVOJz6JGlQl9UlaDzwD3Fit2Ks4VdI9ki6X1BPOrxfc8OAUxQ6wbac19Snvlh96Xv3dsyJqDP2o8ffz2JdFkdNJn0jK3cx2mNnhwAzgSElvr0lyNTDbzN4B3AhcWK8cSUskjUoaHR8f70TuXNDMOyHOvaS9HdyrYjKdeOAUpS+LIqeTPm15y5jZBHATcELN9efMbGv49h+BIxrkX2lmI2Y2MjQ0FEfeXNHMOyHOvaS9HdyrYjKdeOBk1Zft7pT2MXcqRPGWGZI0GL4eAN4L/LImzYFVb08GHkhSyApn/sMdzF527a6/M//hjkj50orpXi8+ObSOt94tb4eye1W0O65RPHD6+8QrW7dPKTOLvowSI6i2DxbOHYosZzfPOnC6T5Qn9wOBmyTdA/ycYM39GklfkXRymOZPQzfJu4E/Bc5KWtAz/+GOKeEHbn/4+ZYKPs0gWovnD3PBh97JvtP7d10bHOhvakyt5OuGH3oZ/N0bEWdc6/VHdQz9faf3g8HElm1TysyiL1utn9frg1Vrxzj1iOGWcuY1uJyTHIUJPxAnnjv4Vu2yksa45m2uHLzs2rpnmQh45PwPdCRv3trqRKd04Qfi4gamcpLGuOZtrrRaP+9E3ry11Ume0it3NzCVkzTGNW9zpdU6fyfy5q2tTvIUJuTvgkP2qxvyt1k8dwg+IPU2dcQ1hCW9+893E8Yj6XGNU+aXVt+7a8dpheEIY9hqzKvvD07vZ889pvHilm1T0nbSB2n0X6f4ZyFZCrPmDlONqgsO2Y9L/viYlvmSmjRJ7/7z3YSdkYYyiFpmbXz2apqNYasxb3dOdNIHeVKm/lmITtQ190Ip96xJ2gjlRq3i0ig+e4VGY9hqzHt1TvRqu+PgBtUUSNoI5Uat4tIqDnu7Y1u53qtzolfbnSau3NsgaSOUG7WKS6s47O2ObeV6r86JXm13mhTGoAqTDVhZnOEYxQjVzjpmVKNWnLXRrNdi2y2jm+u/nci2z0A/UvMn92aGyVZjnpWhs14bJzZPNuKmOUZ5NPAWnb4VK1ZkUvHKlStXLFmyJHL6igGr8pEy4J4nXuTZTVs5fu4BqchYy9wD38iMfQe4d+xFNr26neHBAc754GGT4mgvv+Jent/8GgAvv7qdn/7nODP2HWDugW9su7w4ZcbNk0TeuGUkUWdUOpXt1e07eXXbzobl1xvDalqNeZQ5kTTN2ljpn6df2sL/XfNQamOURbuLyrnnnvvUihUrVrZKVxiDajcPyI5LXnZNZr1zsd0yumlMS0q2KHmLQpQ29kl1P39Fb3sRKZ1BtZsHZMclL7sms965GNeY2EmdUUlDhqIb/aLI3+hzVvS2l5nCKPduHpAdl7zsmsx652JcY2IndUYlDRmKbvSLIn+jz1nR215mCqPcu3lAdlzSCAsbp8xO5EiiDe2W0c1wuknIFjVvUYjSxjOOmlnq8NFlpDDeMl9dPI9HxjdN2aHaTW+ZVrRzsHaaZcbJU+st8fr+aVO8JdKSOY1+a+TZ0alstZ4kC+cOccEND/Jnl63v6i7PJD1XWrWxUvbIQfvlZker05rCGFR9e3J6lK1vy34AednGy2mP0hlU/eDf9Chb35b9APKyjZeTDoVR7r49OT3K1rdlP4C8bOPlpENhlLtvT06PsvVtXg4gT+uM0rKNl5MOLQ2qkl4P3ALsGaa/3My+XJNmT+Ai4AjgOeA0M3s0SUGXLprD0svvZtuO3TaC/r7mB1EXhaxDr5Zt63ej9iycO8SC89fU7ec4Y9BsTtaui1fOKK3QyXinNV7txJl3g2o8utmHUbxltgLHm9kmSf3AbZKuN7M7q9J8AnjBzN4q6XTg68BpiUtba/vNz/6l2DRTAt364KThrZIl9dqzcO4Qq9aONVS2scegwZxstC6+4qoNbN2+s6PxTsu7qFkf5GGeFp1u92Fb3jKSpgO3AZ8ys7uqrt8ArDCzOyTtATwNDFmTwtv1lilrvOeytitvNOtnINYYNCvzyYktbT17ZD3eHmc+fZLqw0S9ZST1SVoPPAPcWK3YK/IBGwHMbDvwIvCmOuUskTQqaXR8fDxK1bsoqxGprO3KG836Oe4YNMvX7vp31uPtcebTp9t9GEm5m9kOMzscmAEcKentcSozs5VmNmJmI0NDQ23lLasRqaztyhvN+jnuGDTL12gn7L7T+2PVlTYeZz59ut2HbXnLmNkEcBNwQs2tMWAmQLgssw+BYTUxurlFPa6XQ5x89dolYOHc9r78nOY0mz9x51azfIvnD3PeKfMYHhxAwGC46/eFzduojdKSB+N1qz7oZJ6m5TVUNLqpwyCCcpc0JGkwfD0AvBf4ZU2yq4CPh68/BKxptt4eh9oPy/DgQCo78ipGj7FwzbRi9Gg1IePmWzx/mFOPGJ70gTdg1dqxnv0QpEGz+RN3brXKt3j+MLcvO56/Pu1wtm7fyQubtwHB+FbGO6153C5R2hJnnsb9XJSRbumwCi0NqpLeAVwI9BF8GfyLmX1F0leAUTO7KnSX/AEwH3geON3Mft2s3LwekB3X6JF1DHUnv5RlfLt9toBTn6gG1ZaukGZ2D4HSrr1+TtXrV4EPtytkHknDuJZWnU4xKMv4dvtsAaczCrNDtVukYVxLq06nGJRlfLt9toDTGa7ca0jDuJZWnU4xKMv4dvtsAaczChPPvVvE3f3Xya7Bsu0QzYI8b41PYnzz0L40zxbIQ/vKRmHiuTtOI8oe39zb51RTunjujtOIssc39/Y5cXDl7hSesntkePucOLhydwpP2T0yvH1OHNygWoUbdYpJmvHom82JevegvvGwk7mVx3j7SX5W8ti+MuAG1RA36hSbNL6Ym80JYMq9/j6BwbadNin9qUcMT4olX11OVBnz9OCRxmd9ai2JAAANI0lEQVQlT+3LO1ENqq7cQ3ybtFNLnBjw9eiT2FHnc1bUueWflWxJLPxAr+BGHaeWpOZEPcUep5y84J+VYuAG1RA36ji1xIkBX48+1Qb5bV5+3vHPSjFw5R6Sl3jxeYt9nTd5usnCuUMNY6/Xmy/9faJ/mqakP+OomcF6fA2vbN2e6/5sNPYeUqAY+LJMSLdCADQ7JBc6OKg5BXr5UOTV68ZYtXZs0jmoAk49YnhS26N4ywBc9rONU+qY2LItt/0ZZezdAJpv3KDaZdI4qDktetlwlmTbG5XVSZlp08tjn3fcoJpTkoyJnTa9bDhLsu1xzwLIkl4e+7Lga+5dJo2DmtMib/J0kyTbHvcsgCzp5bEvC1HOUJ0p6SZJ90vaIOmzddIcJ+lFSevDv3PqleWkc1BzWuRNnm6SZNvrldVpmWnTy2NfFqIsy2wHzjazX0jaG1gr6UYzu78m3a1mdlLyIpaLKMaovBiq4hrO0tot2s1+qW374PR+zODPLlvPBTc82Fb91WWNTWzZtampT5oU/TBPBskiGk1bzZFe2wXbtkFV0r8C3zKzG6uuHQd8oR3l3qsG1bKT1tb0LENDJF1/1u0pI636tEx9nko8d0mzCQ7LvqvO7WMk3S3peklva6dcpzykEZs763jfSdefdXvKSKs+7cU+j+wtI2kvYBXwOTN7qeb2L4CDzGyTpBOB1cChdcpYAiwBmDVrVmyhnfyShpdF1p4bSdefdXvKSKs+7cU+j/TkLqmfQLFfYmZX1N43s5fMbFP4+jqgX9L+ddKtNLMRMxsZGhrqUHQnj6ThZZG150bS9WfdnjLSqk97sc+jeMsI+CfgATP7RoM0vxWmQ9KRYbnPJSlokemlLfxpeFlk7bnRSf31xj5qeb00bzpl6aI5U0I/9E/Trj7Neg5lQZRlmQXAR4F7Ja0Pr/05MAvAzL4DfAj4lKTtwBbgdMtq62vO6LUt/Gl4WWTtudGJ11C9sT/vlHmcd8q8lp4dvTRvEqE2fE/V+6znUBZ4+IGU8W3cvUsnY+/zpj16qb9S8ZZx2qcXDTlOQCdj7/OmPby/puLKPWV60ZDjBHQy9j5v2sP7ayqu3FOmFw05WZMXQ2QnY99u3kqbZy+7lkOWX8fsHjPC+udsKh4VMmV60ZCTJXkyRHYy9u3krW1z5Vi/XjLC+udsKm5QdUpFLxnWKhQxXrwTHzeoOj1JLxrWihgv3kkfV+5OqehFw1oR48U76ePK3SkVRTSsdWoALmK8eCd93KDqlIqiGdaSMAA3ihc/nPO2O+niBlXHyZBeNAA7neEGVccpAL1oAHa6gyt3x8mQXjQAO93BlbvjZEgRDcBOMXCDquOkTLODmYtmAK6m1w6cLhqu3B0nRaJ4wyyeP1w4pZinMA9OfXxZxnFSpKwHM5e1XWXClbvjpEhZvWHK2q4y4crdcVKkrN4wZW1XmYhyQPZMSTdJul/SBkmfrZNGkv5W0kOS7pH0rnTEdZzsaSdcQD1vGAEL5w6lLGW6uJdP/ony5L4dONvMDgOOBj4t6bCaNO8HDg3/lgB/n6iUjpMTKobEsYktGLsNiY0U/OL5w5x6xPCks5sNWLV2rNAHaSyeP8x5p8xjeHAAEeyoPe+UeW5MzREtvWXM7CngqfD1y5IeAIaB+6uS/QFwkQWxDO6UNCjpwDCv45SGZobERortpl+OUxvko1WeIlBEL59eoq01d0mzgfnAXTW3hoGNVe+fCK/V5l8iaVTS6Pj4eHuSOk4OiGNIdOOjkwWRlbukvYBVwOfM7KU4lZnZSjMbMbORoaFirzk6vUkcQ6IbH50siKTcJfUTKPZLzOyKOknGgJlV72eE1xynVMQxJHbL+JiXg8GLTJn6sOWauyQB/wQ8YGbfaJDsKuAzkn4IHAW86OvtThmJEy6gGyEGfMdo55StD1vGc5d0LHArcC+wM7z858AsADP7TvgF8C3gBGAz8Edm1jRYu8dzd5zk8LjwnVOUPowazz2Kt8xtMMmTq14aAz4dXTzHcZLEjbadU7Y+9B2qjlMC3GjbOWXrQ1fujlMCfMdo55StDz3kr+OUgCLHhc8LZetDPyDbcRynQPgB2Y7jOD2MK3fHcZwS4srdcRynhLhB1XFyih9A7XSCK3fHySFl2wrvdB9flnGcHOIHUDud4srdcXJI2bbCO93Hlbvj5JCybYV3uo+vuTtOisQ1ii5dNGfSmjvU3wrvRlenEa7cHSclOjGKRtkK70ZXpxmu3B0nJeIcpl1NqwOoOy3fKTe+5u44KZG2UdSNrk4zXLk7TkqkbRR1o6vTjJbKXdJ3JT0j6b4G94+T9KKk9eHfOcmL6TjFI+344GWLP+4kS5Q19+8TnI96UZM0t5rZSYlI5DglodP44K08YTopPy0vG/feyQ9RzlC9RdLs9EVxnPLRyijaiKieMHHKT8vLxr138kVSa+7HSLpb0vWS3pZQmY7Ts6QZfiCtsj1kQr5IwhXyF8BBZrZJ0onAauDQegklLQGWAMyaNSuBqh2nnKTpCZNW2e69ky86fnI3s5fMbFP4+jqgX9L+DdKuNLMRMxsZGhrqtGrHKS1pesKkVbZ77+SLjpW7pN+SpPD1kWGZz3VaruP0Mml6wqRVdr1yBSycm/yD3Op1Yyw4fw0HL7uWBeevYfW6scTrKDotl2UkXQocB+wv6Qngy0A/gJl9B/gQ8ClJ24EtwOmW1anbjlMSOvW0yaLsxfOHGX3seS6583EqCsCAVWvHGDlov8SMqm64jYay0sMjIyM2OjqaSd2O46TDgvPXMFZnjX14cIDblx1fmDryjKS1ZjbSKp3vUHUcJzG6YVR1w200XLk7jpMY3TCquuE2Gq7cHcdJjKSMtc0Mph52IRoe8tdxnMRIwljbymCaprG5TLhB1XGcXNHrBtNWuEHVcZxC4gbTZHDl7jhOrnCDaTK4cnccJ1e4wTQZ3KDqOE6ucINpMrhydxwnd8SNg+/sxpdlHMdxSogrd8dxnBLiyt1xHKeEuHJ3HMcpIa7cHcdxSogrd8dxnBKSWWwZSePAYzGz7w88m6A43aKIcrvM3aOIcrvM3aFa5oPMrOXZhZkp906QNBolcE7eKKLcLnP3KKLcLnN3iCOzL8s4juOUEFfujuM4JaSoyn1l1gLEpIhyu8zdo4hyu8zdoW2ZC7nm7jiO4zSnqE/ujuM4ThNcuTuO45SQwil3SSdIelDSQ5KWZS1PIyR9V9Izku6rurafpBsl/Sr8v2+WMlYjaaakmyTdL2mDpM+G13MrM4Ck10v6maS7Q7nPDa8fLOmucJ5cJul1Wctai6Q+SeskXRO+z7XMkh6VdK+k9ZJGw2u5nh8AkgYlXS7pl5IekHRMnuWWNCfs48rfS5I+167MhVLukvqAbwPvBw4DzpB0WLZSNeT7wAk115YB/25mhwL/Hr7PC9uBs83sMOBo4NNh3+ZZZoCtwPFm9k7gcOAESUcDXwf+2szeCrwAfCJDGRvxWeCBqvdFkHmhmR1e5XOd9/kB8DfAv5nZXOCdBH2eW7nN7MGwjw8HjgA2A1fSrsxmVpg/4Bjghqr3y4HlWcvVRN7ZwH1V7x8EDgxfHwg8mLWMTWT/V+C9BZN5OvAL4CiC3Xx71Js3efgDZoQf0OOBawAVQOZHgf1rruV6fgD7AI8QOo8URe4qOd8H3B5H5kI9uQPDwMaq90+E14rCAWb2VPj6aeCALIVphKTZwHzgLgogc7i8sR54BrgReBiYMLPtYZI8zpNvAl8Edobv30T+ZTbgx5LWSloSXsv7/DgYGAe+Fy6B/aOkN5B/uSucDlwavm5L5qIp99Jgwddv7vxQJe0FrAI+Z2YvVd/Lq8xmtsOCn7AzgCOBuRmL1BRJJwHPmNnarGVpk2PN7F0Ey6KflvSe6ps5nR97AO8C/t7M5gOvULOckVO5CW0uJwM/qr0XReaiKfcxYGbV+xnhtaLwG0kHAoT/n8lYnklI6idQ7JeY2RXh5VzLXI2ZTQA3ESxpDEqqnBGct3myADhZ0qPADwmWZv6GfMuMmY2F/58hWAM+kvzPjyeAJ8zsrvD95QTKPu9yQ/Al+gsz+034vi2Zi6bcfw4cGnoVvI7gJ8tVGcvUDlcBHw9ff5xgXTsXSBLwT8ADZvaNqlu5lRlA0pCkwfD1AIGd4AECJf+hMFmu5Daz5WY2w8xmE8zhNWZ2JjmWWdIbJO1deU2wFnwfOZ8fZvY0sFHSnPDS7wP3k3O5Q85g95IMtCtz1gaDGAaGE4H/JFhX/Z9Zy9NEzkuBp4BtBE8PnyBYV/134FfAT4D9spazSt5jCX7m3QOsD/9OzLPModzvANaFct8HnBNefwvwM+Ahgp+1e2YtawP5jwOuybvMoWx3h38bKp+9vM+PUMbDgdFwjqwG9s273MAbgOeAfaqutSWzhx9wHMcpIUVblnEcx3Ei4MrdcRynhLhydxzHKSGu3B3HcUqIK3fHcZwS4srdcRynhLhydxzHKSH/Hx0ARVv23WaMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.figure()\n",
    "pl.scatter(X_cos,Y)\n",
    "pl.title('Cosine Similarity VS golden score')\n",
    "pl.show()\n",
    "pl.figure()\n",
    "pl.scatter(X_wmd,Y)\n",
    "pl.title('WMD Similarity VS golden score')\n",
    "pl.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can learn a simple regression model between those 2 quantities. Use a polynomial of degree 2 to learn the regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.polynomial.polynomial as poly\n",
    "k_cos = # TO BE FILLED\n",
    "k_wmd = # TO BE FILLED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compute from your regression model the estimated relatedness for all the pairs in the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_cos=[]\n",
    "X_wmd=[]\n",
    "Y_test=[]\n",
    "for i in range(len(testA)):\n",
    "    s1 = vect.transform([testA[i]]).toarray().ravel()\n",
    "    s2 = vect.transform([testB[i]]).toarray().ravel()\n",
    "    # cosine similarity between bag of words\n",
    "    d_cos=# TO BE FILLED\n",
    "    X_cos.append(d_cos)\n",
    "    # WMD\n",
    "    d_wmd=# TO BE FILLED\n",
    "    X_wmd.append(d_wmd)\n",
    "    Y_test.append(scores_test[i])\n",
    "\n",
    "# Final regression scores\n",
    "Y_cos = # TO BE FILLED\n",
    "Y_wmd = # TO BE FILLED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use MSE, Spearman's rho and Pearson coefficients to measure the quality of our regression model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error as mse\n",
    "from scipy.stats import pearsonr\n",
    "from scipy.stats import spearmanr\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimate the quality of your regression model for both Cosine and WMD dissimilarities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------- Cosine\n",
      "Test Pearson (test): 0.4391501309646927\n",
      "Test Spearman (test): 0.36660337782004127\n",
      "Test MSE (test): 0.8784692440428779\n",
      "-------- WMD\n",
      "Test Pearson (test): 0.7062364858939475\n",
      "Test Spearman (test): 0.5872407398045967\n",
      "Test MSE (test): 0.5858116558928291\n"
     ]
    }
   ],
   "source": [
    "print('-------- Cosine')\n",
    "\n",
    "pr = pearsonr(Y_cos, Y_test)[0]\n",
    "sr = spearmanr(Y_cos,Y_test)[0]\n",
    "se = mse(Y_cos,Y_test)\n",
    "\n",
    "print('Test Pearson (test): ' + str(pr))\n",
    "print('Test Spearman (test): ' + str(sr))\n",
    "print('Test MSE (test): ' + str(se))\n",
    "\n",
    "print('-------- WMD')\n",
    "\n",
    "pr = pearsonr(Y_wmd, Y_test)[0]\n",
    "sr = spearmanr(Y_wmd,Y_test)[0]\n",
    "se = mse(Y_wmd,Y_test)\n",
    "\n",
    "print('Test Pearson (test): ' + str(pr))\n",
    "print('Test Spearman (test): ' + str(sr))\n",
    "print('Test MSE (test): ' + str(se))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not bad isn't it ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}