path: root/notebooks/plot_partial_wass_and_gromov.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Partial Wasserstein and Gromov-Wasserstein example\n",
    "\n",
    "\n",
    "This example is designed to show how to use the Partial (Gromov-)Wassertsein\n",
    "distance computation in POT.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Laetitia Chapel <laetitia.chapel@irisa.fr>\n",
    "# License: MIT License\n",
    "\n",
    "# necessary for 3d plot even if not used\n",
    "from mpl_toolkits.mplot3d import Axes3D  # noqa\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample two 2D Gaussian distributions and plot them\n",
    "--------------------------------------------------\n",
    "\n",
    "For demonstration purpose, we sample two Gaussian distributions in 2-d\n",
    "spaces and add some random noise.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2de5BjV33nvz9dvVr9nu72vKd7/ARjiMG9GAKb2oDJOg6L2a3UFqShTJKqWcKSdTapJSbzx/T8MbXksSSpIo+aJSROeTaEIrAmKWfBJk6yqQrGPV4He2xsj830eN49r55+6/XbP845OkfSVbfUuuqrK/0+VbcknXvvOUetmZ9++j2JmSEIgiBEl1jYGxAEQRCaQwS5IAhCxBFBLgiCEHFEkAuCIEQcEeSCIAgRJx7GoqOjozwxMRHG0oLD8ePHLzPzWFDzyefaPgT92QrtTSiCfGJiAjMzM2EsLTgQ0Wwd13gAZgCcZeYPr3etfK7tQz2frdA5iGlF2IiHAbwc9iYEQaiNCHKhJkS0B8DPAPhy2HsRBKE2IsiF9fg9AJ8DUAx7I4Ig1CYQQU5E/5WIThDRi0T0F0SUDmJeITyI6MMALjHz8Q2uO0BEM0Q0Mzc3t0W7EwTBpWlBTkS7AfwXAJPMfBcAD8DHmp1XCJ33AfgIEZ0C8FUAHyCixyovYuajzDzJzJNjYxIkIQhhEJRpJQ6gh4jiADIAzgU0b1cwPR32Dqph5s8z8x5mnoD6Yv47Zv5EyNsSBMGHpgU5M58F8DsATgM4D2Cemb9TeZ38BK/N4cNh70DYMo4dAyYmgFhMPR47FvaOhA4gCNPKMIAHAewHsAtALxFVaW7d+hO8HbXtRmHmv98ohlyog2PHgAMHgNlZgFk9HjggwlxomiBMK/cB+BEzzzFzDsA3APx4APN2BLW07elpgEgdgH0ehODvhC+PjuTgQWB5uXxseVmNC0ITBCHITwN4DxFliIgAfBCSQLIh09NKKTN9PczzIISwmGralNOnGxsXhDoJwkb+DICvA3gOwAt6zqPNzhtlWqltb7Su0Mbs29fYuCDUSSBRK8x8iJnfwsx3MfMnmXktiHmjSqPa9qFD1fdvZs3Dh7f+y0NogCNHgEymfCyTUeOC0ASS2bnF+AnWyrHNmEbMHK0w1QgBMTUFHD0KjI+rb9nxcfV6airsnQkRRwR5i6nUttdzfm4GPzNOM/MJLWZqCjh1CigW1aMIcSEARJC3mHoEqmsWadQ04mfGOXRIBLkgdBMiyLeAjZyfRksPKopFhLggdBciyLeAWs5PwAp389x93SiVZhxBELqDUDoECdacUkkzwlg0cUHoTkSQbzFGUE9PW8FrtHCjpQuCIDSCmFa2mHpjyQVBEOpFBHkbIFEmgiA0gwjyNkCEuCAIzSCCXBAEIeKIIBcEQYg4Isg7BDHPCEL3IoK8Q5Aa5ILQvYggFwRBiDgiyCNMWA0sBEFoL0SQR5jKGi6HDkkNckHoRkSQdxBiJxeE7kQEeQcwPS0p/oLQzUjRrIhTWUXR2Msl7V8QugcR5BHHCHJmJcSlgqIgdB9iWoko0qtTEASDCPKIIr06BUEwiCDvIESIC0J3IoK8A5CIFUHobkSQdwD1auKisQtCZyKCvIuQhCFB6ExEkAtCUBw7BkxMALGYejx2LOwdCV1CIIKciIaI6OtE9EMiepmI3hvEvELzbLawFhHtJaKnieglIjpBRA+3eq+R5tgx4MABYHZWhRDNzqrXIsyFLYA4gAwSInoUwP9l5i8TURJAhpmv17p+cnKSZ2Zmml5XaIzKhCEiOs7Mk/7X0k4AO5n5OSLqB3AcwEeZ+aVa83f15zoxoYR3JePjwKlTW72bdT9bofNoWiMnokEAPwHgTwCAmbPrCXGhMcJyUDLzeWZ+Tj9fAPAygN3h7CYCnD7d2LggBEgQppX9AOYA/CkR/T8i+jIR9VZeREQHiGiGiGbm5uYCWLY7CNJBudkwRSKaAPBOAM/4nJPPFQD27WtsXBACJAhBHgfwLgB/xMzvBLAE4JHKi5j5KDNPMvPk2NhYAMsKjbIZ7Z6I+gD8FYBfYeYbleflc9UcOQJkMuVjRMrcIo5PocUEIcjPADjDzEZb+zqUYBc2Sbt0/iGiBJQQP8bM39ja1SPG1BRw9KiyiQPlDglxfAotpmlBzswXALxJRHfooQ8CqOkQEzbGr47KVnf+ISKC8nu8zMxf3LqVI8zUlHJsjo9Xl6FcXgYOHgxlW0LnE1QZ218GcExHrLwB4OcDmlcIj/cB+CSAF4joeT32G8z8RIh7igbi+BS2mEAEOTM/D0BCnVpAs3VUpqc3p8kz8z8BoOZW71L27fMPRdy2bev3InQFktnZ5jRrTpG0/BA4cgRIJqvHb9wQO7nQEkSQC0LQTE0B/f3V47mc2MmFliCCvANpl6iXrubqVf9xsZMLLUAEeQfSDlEvXU+tRCCxkwstQAR5hBBBHCGOHAESierxhQWxkwuBI4I8BDYrkDfjuJTuQSExNQUMDFSPZ7PB2cmlbK6gEUEeAlsZSSJafIi00k4uZXMFBxHkbY44LiNMKwtpHTyoskVdGsgeJaL7iegVIjpJRFW1kYRoIYJ8i9isQBbHZYTxK6SVyajxZmkie5SIPAB/AOCnAdwJ4ONEdGfzmxLCIqgUfWED3AzLygYPQocyNaUeDx5UAnbfPiXEzXgz1MoerU/bfzeAk8z8BgAQ0VcBPIgaNZKSXg/3xAcAz7ODuZx6jDsiJJfXT5x/3CYxKl8oDXFCzUNrOXud/g9R7E+Xhmhe/eIg12nsKd2TV9fsWKZHXVewa6D0XGlOueGUXcp5G6Uxo2A5Wy+at+aouzG9bNFOB8qjGj0PFZ0xKl/LnZudNUr3uPdq1s6fuczMVWVGRZBHCHFcRpCpqY0F97FjVtib8MSrV9cX/EeOKJu4a16pX9vfDeBN5/UZAPe6FxDRAQAHACDt9ePHb/oYeHXVnk8rgVu4eKk05u3arp44WkrxivIT0O32C4bfUL8a6K3jdsGij9TaoyResccK8pVdqtVB7ytX7HxnzqsnrjTTc5t7s0NW8q4NKkleSFmJaoSrl7N7Xx5T0rVgv1sw9Lr6gpi/2X4bJOfVPa7Q9rTATyzZQbOeK7RzGf1F02f3kljksjnce577yq/5fHuLaSUUNiuQxZzSgVQ6La9cUcdGDky3bC6Rejx6NBhtH+V15pOxnkDmFFqHCPIQEIEslPBzWrqs58A0ZXOLRfVYvxA/C2Cv83qPHhMiiphWWkijlQc3W6lQiDD1hCIGn9b/LIDbiGg/lAD/GICfq3l1jMCZNIo7bVZq7FW1p9jtN5fG+Ipu1btmbQJ0szapnLdtAPltt6h7LzqtffPK0LxwrzW39Pz1cQCAd4szph+Ls2fsve+4Xc13fcmue31RPf5oHgAwP/W20rnsoDbZONLPmC5ieWviWB1VZpFi2ppHqKBuWhq39vjsdWPotvPFsmqexJI1wRRS5WsBQD6jbsr32jXiS+peb9Uxpm+gcotG3kIajReXSoVdSD3OyYD7fjJzHsBnAXwbqqn215j5RKCLCFuKaOSCECZ+TkuXoMIVK9ANQuprEpLLg89eAE5bLZR27VBPrtk2rkYhpV6n9/qCel+Fm3eVhjyjiced8BEdjdL36jU7dofW9gtWW80NKc9jat9uu5dzyvFpHKsAgD07AQD5u/YDAAZ/ZENLVrepdV3NmGPaEekowfEV7exM2QsHf6T+BuzsPX3ZhKjYe6mgxlLz1c7OQsJeaByp2X67Rlz/U0gs1R/aJhp5wDQaLy4JP11OpdNyZEQdLXBgCp2LCPKAaTSBRxJ+Opx66qG4TsvLl9XRuANT6GLEtCIIrcKEFhqziQknBCIloJkZnMvDGxm2Y4s6WSdt47OL5y+osZQdi42NqMflrL03rZOErjjOTlYmiMJOu0Z8Vseoxxyzw5I221xfsLcO9qknbsLSvHJ2enp/S7ttMHi2rzqe25hWXPNITi9VSFkTR7Yvps/ZMW+tuiOiSRKigl3EODuLjmklr723eccaVVKvqXreWohG3kIajReXhJ8Oo8l6KIJQL6KRt5B6zCNuyKGYUyKGm5Hpl4XZRD2UtiKTRvHutyJ2zX4p5YdVDZnEOeucLLz37eqJk9mZPKMckGu7bEnf9AkVOrj03lvtGvqW3pculIaW7lah7oUeq28ubVfPBwedkMTn9d/zFhsavzSu1isklVbb9+aKvf6yEnvZQTf+UO933jpFl3aqXw4m+xIA+s+oXxbs2Z6siWX1a8I4OAEgeUM5Rb01x9mZ1qUJio6G36/G8j12jeSCuieWt9fl0+vr3KKRh4yEHEaUesrItrL6oSA4iCAXhM1Qj9mkVvXDBx6QhhBCoIhpJQSmp8s1cePTOHRo/egWMb20EfWYTfyqHz7wAPDoo9FygC6vIvbsyyDH2Rm/rB2VKWti8L6ncooo6VQrHFXZoKkz83ZMVyvMPPO6HdPmhvwd1jzS+wNdNcBx+vXsHVXrnzxXGqN+5Snkk/Zv33tROUB5TK2/8JbB0rnVYaW/Fh3faMnZuduKxDX9dgtpa+KI5dX7XZiw96auVuvDtENNnrxh7/Vzdpo48pzj7DTx6/F1KjdUIhp5CGwm5FBMMG1GvWaTynooTzwhDlAhcEQjF4TNsNkyshF0gBIRKJkAu0X0dW0UN/eQPFNzxBld02GHrpau65ZT0mrzrOuVuw7DEgkrpjhuiqK4cYJ6LznrqIzpuUuzudNq/yO5aixXl6KlIpU9qudcNQafLVPBXO/OZ86xM0a1r3OcorxBKKJo5CGzXsihZH22MZstIysOUKEFiCAPmfWEsmR9tjmbKSPbyvZvQtcSmGlF9wGcAXCWmT8c1LyC0FG0sv1bq0glgVv2opiw3sHYKdWVhxxnJ4ZUYSxatVmcpZZwFy/bMe0ALYxaB6Rp05Z1WrL1zCnxlNtlnazZQWWiSfbaL0NTwCo2bgtpFYaVs9OUtl0bcDIsdcx2vuL7FAA8G25eKjHrdggyJXDNObV308PNXmcclVknVdRtD2fIZcrXcudjx3xU3EBSB6mRPwxVElNoAZL12UH4afL11GQxNHKt0BUEopET0R4APwPgCIBfDWJOoRwxp3QwjdRkCat+SyyG2ILt2Wm0ar7oNIwwGrbj2OTTOkzw5j12TIcJFnZbTds49lLnFktj+R1DAIBcnxVTq8NK+06P9pfGvB+qNpY0ahtf5PrVLwXqU3vpvWAdofmM0l9N3RSX5KLbREJnXTqae+8F9cuhkLK/TkyPTddhacbia052pi5j6zZ4zvYaDd/N7NTzOb2ki46v2I+gNPLfA/A5+PZ9VhDRASKaIaKZubm5WpcJQvfRSE0Wqd8i+NC0ICeiDwO4xMzH17vObeY6NjbW7LLCFkBE9xPRK0R0kogeCXs/HUsjIYkNXEtEXyGiS0T0ojO2jYieJKLX9ONw1Y1C5AjCtPI+AB8hogcApAEMENFjzPyJAOYWQkI7r/8AwIcAnAHwLBF9i5lfCndnHci+fcpE4jfezLXAnwH4EoA/d8YeAfBdZv6C/nJ+BMCvr7u/XA705gUUnA483ltvAwBQX19pzPTJdLvnxLS5I9dnPX2JYWUySVy2ZhST2UnZXGmItXM1fdZex54yqcRW7XWkHZ981RbwSvWo9da2q/3FV6ydwphWkkvWgFD0dPeeZHW8dswxccRXdIEsJy3ULwPTlMh1i2axDlx31yAfG4ZxbCZXnVK5uerryva4/umNYebPM/MeZp6AauL6dyLEO4J3AzjJzG8wcxbAVwE8GPKeOpNGQhIbuJaZ/xHA1YrhBwE8qp8/CuCjm9u00E5IZqdQi90A3nRenwFwb+VFRHQAwAEA2CdJLZujkZDE5sMXtzPzef38AoDtfhe5n2s6PgAMDyLufoHoXpz5s7bmiXfLhLrXybDMz6p/QnHHAVowDtKx2+x6RaWaspPFSTntWByw2nxel7SlvKNNz6u+oW6v0GJGOTvjy0qVXRtxGmDoWie5nmrt23VOGgdjwQkbzPYpTbzgRF3mzI8Sx4mZWNZrZNzGEqYvqFNrpVR/xd7La+oxnw4n/BDM/PcSQ95diO/DoZmwwEaSizaTiOQDq5x73w6/7uea9Ho2Nb+wdYhGLtTiLIC9zus9ekzwIzpt3S4S0U5mPk9EOwFcCntDQvOIIBdq8SyA24hoP5QA/xiAnwt3S23MemGB7SXIvwXgIQBf0I+Pb3hHoQgsLKG4YJ2OpV6crrPTZHQ6pgNvWAXFcI+1RcR6dKrkmuPBK2hTScF6Fln32/SW7XXemrZFuIW5dK9OXrB9PGlF7a84qNaKOU5H0uaO+JqdotS/05nWODmL1lKE+JrpBmSNGTHztp17TVy8l3MKX+lMTfbcsrg6ttxxqJp5Yjmf7NEatH2tFUmECQdmzgP4LIBvQ2Xsfo2ZT4S7qzamDasaEtFfAPhnAHcQ0Rki+kUoAf4hInoNwH36tRBx2l4jP3xYhHlYMPMTAJ4Iex+RoLGwwC2BmT9e49QHG5qnUEDx+rzVpAFAhwlSxtrP+YbWiF2HpTmfsyonDeiszHk3/FBpusWbbFh77Ib+heNo+IlFtQdatEVRaFD158yfs/0+4wuqxgrrDMzCNvuLwDgg2VFj2UQTupV6tW/XdWya3plu/RXjxHQxzslY3i5inJeu4zLvW/elunYL1lfI218jF4RIIFUNhRBpS0EudbiFtqKeaJTN1icXhABoS9OK25+SqNyvIQhbSiPRKFNTHSm4Ke7BGx5C4dr10pjJ2MRlm29E24aq7uVFZeIo9tiCVrE3lQmG9+2y92rTSmx+qTRWHNJ9Nz1rV8gOKpGVGLROVsyqYCrX9GPuLfQpu0hi0XosTRx5wemdacwsbhanMZnEnDjyhC6qlViyOnBi2TQMsNd5WS67HrCOzWLcyewsreeYj5a4fF6Um4H8aEuNXBDaBilSJUSAttTIXaQOtxAqbRiNsuUwg/P5siHSPTbLfiybjM6Yox+axpiF6nBBNwPU/Ox267RAZ28SW201ljclXh3VWfcKLWazznUmnNHMa/dkyqSUOTv187JWnF759WoenzEfddho3excZ+q5lN3roeo633VjEn4oCJtHemwKEaDtBbkghIpEowgRoO1NK4IQKlHssRk0nqdiv8d3lIZY9+zkPTvtdRd0MayCU7p1t6rJ5Z23PTsLt6puQbFTNu4bRWUqWXz/raWhnse/DwCIj9tKET3aVFJ45Q177+Sdao1r1lFKV5RjNnFOVSA4//G3ls5ldSMj9jGPuGVl10Z0bHvarTWrRObizdYslB12JtLEssoUEl+yIrao49HdOPJCj3o/+V5rKoovq814K/Xr2SLIBWEjOjQaRegcRJALgrA+sRiKAxnkBm14X+yucQBAtt/WX80YZ17KipWVnSqz09tje2yaxg4Dl63Jyjg5s/1WC83c8zZ17oxtDUlLKoIo9o47SmOLe9Q8aWfdpK77krt1p96TfTv5XqUFx3JOmdiEGis6WZzFHq2JJ52SuebtppxfHabWi+OPLCSr1yik9a8JZy9mXe6x85WSQYtOPZdyX3MVYiMXBEGIOCLIBUEQIo6YVgRBWJ9cDnTmItKXncxJneWZSFpbhIk1j+Vs2dm+15XZg/ps9x7oGHQ33tz06hx81TosvXNX1LxDA6Wx7E71PPnymdLYwOV5tSenZye2qyYniStqvr4zTs/QBR13nqguMWvi1AFgdUSJR2MSAYDec9o80mNNSsnruhStWwJXh7QnFu2gMam4TlZjZso7nYRMD1Bv1V4nRbMEIUya6RokCHUiGrkgtIrodA1an3gcuGkEWHJKx+68CQDA5y7asf06TDDvNIeYVZoz77GtQYsvvKrG7r3L3qtDFr15u0Zh96iart9q/Us7dC/OBaet4Isn1b2jI6Wh7G7d0EJndKbm7Z5MPZVsr6vHmvomjmNTOzELTm/P1IK6OXndqtXJBa11O1GKSV0vJb7qODFT1XpzTje3MOGK7nxuY4liIuKZnYIQWcKt05IgoqeJ6CUiOkFEDwMAEW0joieJ6DX9OLzRREL7Ixq5ILSKWvVYZmeVtt56rfzXmPk5IuoHcJyIngTwKQDfZeYvENEjAB4B8OvrTbK63cPLvzqA9KC1M6+dV7bv2Ohgaax41bSEt/fGRlTYIZ219vXEzSo5Jztb3dT5LffYxKGT/6y07mLSaqapy0r3vHy3tZunDqokouUFu7/YnNLck9fV9Z/4j98tnbuzR1VLTJO15Zvnq2xt3/emlM192LNhkn92Q/0S+dSAbXV6fE0ZxD2n8swNVnt5I3tTaWwsfqNq3V2eqgQ5HreieFb7Gi4WnAqPmg/+adUQANHIBaF1rFeP5cCBVtvLc8z8HAAw8wJUu77dAB4E8Ki+5lEAH23lJoStQQS5ILQKvzothi0shUtEEwDeCeAZANuZ+bw+dQHA9hr3HCCiGSKaKSwu+V0itBFiWhGEVmFMJ5/4hP/5LSiFS0R9AP4KwK8w8w1y+l8yMxMR+93HzEcBHAWAnh17uf+lJHID1um47ay6bWXUflGlrukSs45pZXVEmU/637SD80vKLDJ0pnrpExlbV2XHC+p8LmP33HNVORtdR+X8mpqvd8Fel76i7s3MKTPFX9xzj5134Ha1X8+mS8b1plcLViS+fegcAGA0YXuL/vXZtwMAzu2yfchfWnTqzWgWc8q0Mrdiwy77k8qzmfSs43Ukpb4kd6dt046zq6pBx7U1+7eNlf6or1atBYhGLgitZWpKtX3zo8WlcIkoASXEjzHzN/TwRSLaqc/vBHCp1v1CdBCNXBBazZEj5WGIwFaVwv0TAC8z8xedsW8BeAjAF/Tj4xtNQgWV2OK2XIvrNmTxFbdFmb6+aDXtvNamvTU7ZhJdUgvVGrnntlBb0u3fnE70qWvGUWidkt6Kp/dk50kuFPWY0n5Xlq0j9GpcaeIJz/5K8GK60qFTMOX8qnLkrjnlCq8vq18YF7LW2XpxRTl0Y46zcyWv9je/bB262byaJxm3vwTyup5K3PkZc2VNafFXV2qY5XwQjVwQgsBN/BkdVYdJAgLCaMzcB+CTAD5ARM/r4wEoAf4hInoNwH36tRBxRCMXhGapTPy5csWeM0lAR48Cp05t5a4WmblWFskHt3IjQutpWpAT0V4Afw7l/WYAR5n595udVxAig1/ij4uJUIlSNqdDMQks7gNyQ9ZJl+9R5gzTfAEAcn3mB779/lgdK5ZdDwDL48q0QG6HBX1LYnyhNHT9XH/ldCgklMN1aac1JixPKHNLdt7OV0ir9VZG1XV37vlR6dxtfdVugR5PzeE5Jo539LwJABjxrLNzIa/i4f/D8PHS2DOpWwC4Dkkgp9/b68s2A3Uoof6NDMRtEZXRhHq/E0lbqvdsbpt6zNpcrZxu4PnPVTtXBKGR5+GTeMDMLwUwtyC0P/VEn3RTs2Zhy2lakOuY1PP6+QIRmcQDEeRCd7BvnzKhbHRNRKEikFgguC61nkumm70dS97Q1xfce9X5PifUsKC1c7cyoPETLl61GaADl9Vgtr/aQmTCCwEgO6TEWHyp2hnbe0FpyWfmbQZqtqDW70us2eu1s9OEDQJAQfd/M1ozAJxeVFry8/3283xjRdWEKbITEplT72NuxWZn9ifV87jzBxpLK23/WtqGKc6uqJoxC3knUxW+UaLO+QCpSDyoPFdKMJibm6s8LQjRZb3EH0CaNQstJzBBXpl4UHmemY8y8yQzT46NjVVPIAhRZWqqPCplZEQdWxehInQ5gUSt1Eg8EITuoYMbNHMMKPQA+Yz9eZ/TvTXz/XYsltcNFpzMztyANo8MWrNHrk9dUEhX65Fen42xzvUq8eT20TQNIPK9dr58n9vlXmHi200P0OGMLY+7O6MaUbiZnUndFDObtCJxf0pZDkyxKwAY61GmkImkLe51Le80zdAMxFNVY8NJ5exMOQ04d2h71HjKzmfMKJdyts+ph+r36NK0Rk4q59cv8UAQBEHYAoLQyN8HlXjwAhE9r8d+g5mfCGBuIQSI6LcB/DsAWQCvA/h5Zr6+/l1Cp8JxxupNecQGbPnV1axy5uVHrHZZjCsnIhUdbXlInV9dsZmYPKrKvi5nrdZq7tg1av+Zzd2kapjkexxHadJo5I7zb1Q5LbNpZw29F5ON+v5h2xruth7VDKPf6aWWILXP6wWrXb+nR4UsbncyQF8ZVE7tH0/bhhoZsk5Tw1VdgnZHar40ts1Tqa/pmP077o6rUrl7Ha1/R1z9Da44ZWxNidzfq1pJ0bRGzsz/xMzEzO9g5rv1IUI82jwJ4C5mfgdUlZ7Ph7yf1iBt2IQOQVL0hSqY+TvMbFSt7wHYE+Z+WoLJxpydBZhtBqYIcyGCSIq+sBG/AOAvw95E4KzXhq1DnZabhgmxtRiKqzY7M7GmTBaFNasLmr6TZXHk+rzbEZ71PN5adXz44pr1bHraYuEmgHq6O71bwGtNz0fuXvTcZt2rOWsyuZxQTkS3G1BCb3o+b8NI55LmHluP/ZqeZ86JnzdmlIKjF1/Nq7HLjsPSxJknnD+QMem4XYPm8qogV7lppcXOTiGa3HfffQDwNiJ6seJ40FxDRAehMndrqqmRzQ+olWkpGZhCBBGNvEt56qmnQEQnmHnS7zwRfQrAhwF8kJlrppW5DQgmJyfXTz9rJ2plY0Y4A7NVxNaAvtkYsk5jid7zurGEo6WbbEu3VcXKsjrfe85xWGqnZP9s9T+Xy6NDpefbZ3WYYspq35k53VtzyK5bSKl9OQmY6JlTc5tGFN8/az/X04MqO9MtHZvQzR5MWVkAuKK175GE1cifvnCb2of5aQDgxYVdVe9jVZexvbhsteq+pLrHDXvcphtL7HEaS5xfU1moppwt4JbX/YeqtYAu0Minp8PeQfQgovsBfA7AR5h5nWpQEcYvG1MyMIWI0vGC/PBh9SgCvSG+BKAfwJO6jvUfh72hwKnMxpQMTCHCdI1p5fBhEeb1wsy3hr2HLaGDszGDhOPA6igjN+g63HQZ21HruGNPmzscH+bqmD5P1hSyulOZFuLL1eInOfmdHWAAABM4SURBVGZ/AC7t0o5CxwLjZdU8bhnbtR3K3FLosfORNkXk+tTY/pGrpXN39KsYcNc5mYkps0ciZt+PXxnb14ZvAgC8t/e1qr27TkxTcOu11E2lsW3aROOWsd0WV2M3p2xc+tm0KmN7em2kNJZj9b7/tmpVRUdq5NPTSskyfWbNowhyoZsgojQRfZ+I/oWIThDRYT2+n4ieIaKTRPSXRJTcaC6hvelYQc4MHDpUPn74sBLqItCFLmENwAeY+ccA3A3gfiJ6D4DfBPC7+pfXNQC/uN4kTEAxAXCC7eEB7EEp5vooJvQRt4c5x2QPeAx4XHaduTeRKJQOs0YhZY98mpBPEziO0oE4A3FG0TkKSaCQhPp1QKpkbemIq2M4vlw6BuMrGIyvIBPLlo7+2Ar6YysYco7++Cr646tlY8PxJQzHlzDoLZeOjLeGjLdWur4/vlpa15zLeGsY9JYw6C2hP7ZaOgZiKxiIrWAwvlw6tsWXStq7Hx0pyA1GYJuYC2Z1iCAXugFWGLtAQh8M4AMAvq7HHwXw0RC2JwRIRwtyoForF4Rugog8XQPpElTphdcBXHcyd89ANYKpvK+UH1Bcqq0JCu1Bxzs7jfbtCvTpadHKhe6AmQsA7iaiIQDfBPCWOu8r5Qek9+xlYpQ5HUux4o7/k/ySD4sV11fcYxc061pPqZmPHOdpaY2y+ahqjcrr3Phw0/8y5tzgcfWmilrPzTr6bl47HY3zEXA6Azn7NGN55zrjAHU7CeU4XjWWNWs4Ka1FVGfBunS8Rm5wBbcJSRSEbkFXr3wawHsBDBGRkRJ7AJwNbWNCIHS8Ri4I3QoRjQHIMfN1IuoB8CEoR+fTAH4WwFcBPATg8XUnKgLeMqEYt3pfQlddLSbsWFxHDrqauTmfumq137Vh3WPT9noosbhkA2iGFtQ9+Uy1lu5mcXo3lAYbX7R7MTVWUtfVDZed3plpnVnZG7flZxN64pWCrb8yrJ2LN+K2j+ilVTXPqdxoaex81vYDNdzI9wAALqzYWitrBfW+52J2L8up6oCh01kVdngla6+L+f7csXSVRl4ZkigRLEKHsxPA00T0AwDPAniSmf8GwK8D+FUiOglgBKoxjBBhukYjd+3iRDaSRRA6FWb+AVQz9MrxNwC8e+t3JLSKrhHkgiBsDvaA3EARhUGbubi2qntijhTKL6wgqzM/Y1kragojKosyd6ParDA4aDM7V7cp84TjB0RCB1NmbW0tFEdUZmc2YdeIZc1e1M2399vMzrf0qixKzzFXmMzOguNUfHtaZXaOxeyenu29os6lztl96nK4MceLu8rqvfV420tjo3rzfU5NX9MP9JaErRw6ElfXnU0Ml8ZyPn9bl64xrbhISKIgCJ1E5DTyIEIHxS4uCA2gMzsRtxon+6iAbKSJj9myLPxQvygkuep83OmPuarXKKTtdTnt+CxTUM3NzlgxpceW1PV9cVt21vTMzMRcZ6f65eA2m+gldU9vzO6pT3e7yDh1VfpjSsN2NfyibjxhrlfrZcvWUnOr82mfZhO9sepeoLWInEYuxa8EQRDKiZwgByQOXBAEwSUSgrwydBCQ0MGu5tgx1fU+FlOP0jC59RCDYvYwxagQ49JRKopVdjBAXFY0q2oOUqYajgFerFg6UDVX9fUcs/MxuUf5XmLg0pGgAhJUQIy4dCQojwTl4aHoHAwPjBhQOjwqwqMiEgTnUPfGUKw6zPUeFRHTR5LypcOu4ewPBSRQ8J2vFpGwkRuBXamJS9OILuTYMdXt3jROnp1VrwGpLS50LZEQ5IB1chqtXOLAu5SDB60QNywvq3ER5K0hxuBUEYmU7TVZSKvwuliPM7aiu9m79Vd6lBOvkLGiJpFW9xSTqaqlBpLWwXdVOzmL1v9Yeu46QJN6vtWcU/8kqZ7ndVLmSNI2hzAhf8bpCdgu9qtsQwNHPTU25tl9jsYXqsZM4wnP8fKa+ZaL9rptOqzQdWKOeWq+HY6jdk0XrCw2YDCJhGnFRUIHu5xaXe5rjQtCFxA5QT49rYS5mFO6lFpd7muNC0IXEDnTinkuzs4u5ciRchs5AGQyalxoCd4KYfDFBHL91sYxdFaZEVbGekpjaV0YyzWtrI4o20b/m3bw+g1VDKrvjGMf1U9fS+0pDe16Xt2T7bX6Zs9VZUbJZezY9SVVmGrAWk+QvqImzFxSsdtfvfue0rnRwTsAAEnPxm7Hdaz4Wt6KxLcOXwAAjCRsPfZvn1FVgE/vsf00f7hgszcNCzn1vq8sZUpjfSkdR+6sO5JWc+9I3yiNnVtRRbjms/ZvGyuZbU5UraXORwQJORQAKDv40aOq6z2Rejx6VOzjQlcTiEZORPcD+H2o3KovM/MXgpi3kunpcoFuHJ9iaukwjh1TzsvTp5XJ5MiRckE9NSWCewspJIHF8fJaK/leJTqyI9bZmb3sV2tFa9B9VtSsTihnH8ccZ6f+vzy039ZEuXZhpOwcABQTao3l7XZw9VY13+oNu0auV123OqKcsvfue7V0zq/WSko7Pl2H5Y/1zAIor7Viap5MDT9TGvteen/VfFmd5vraitXWh3Wd38G4nc84XicSl0tj5/KqxsqbWav1mxowT8GfpjVyIvIA/AGAnwZwJ4CPE9Gdzc7rRqkYgW2EuHF4uj04RZBHiPXiwE144eys+nBNeKHEigtCTYIwrbwbwElmfoOZs1DF6h9sdlKTim+ENbB+82QxvUSEjQT1euGFgiD4EoRpZTeAN53XZwDcW3kRER0AcAAA9gUQYSBhiBHDmEtmZ6vPuXHgEl7YfpAuUhVznJM+vTNLrSvdBET93AnZLvXYZB/pQ1SdIFJ0qt3me6p7VxpFz838NnObWlQJZ1PGBJLyiSN34749vXnP2VMqlq+az6+41aoOeHcLZKXL/ghmDS57dDHFswAg5fN3cdkyZyczH2XmSWaeHBsb871mvS4+lYLbz/QiXX/aFFcLr4UR1BJeKAgNQ9xkiiQRvRfANDP/W/368wDAzP+91j2Tk5M8MzOzwbyNZW9K15/GIaLjzDwZ1Hw1P9eJifWFOKCiT06dqk7BVxtVH+74eLXjU/AlyM+WiOYAzAIYBXB5g8vbnai/h3FmrtKEgzCtPAvgNiLaD9WN+2MAfq6RCcRZ2eFsZBZx48CNkDZmGPcbWuqqhIIRHEQ0E+QXfxh0wnvwo2nTCjPnAXwWwLcBvAzga8zsH7VeAz9HZaM2cLGZtzHrmUX84sCnppR2Pj5e/TNLHJ+CUEUgNnJmfoKZb2fmW5g5kBS7RjV00ejbEBNmaDRrl0wGeOwxJbBradfi+BSEuggts1MclR1OpYOT2X7Y9WZjiuOz3Tga9gYCoBPeQxWhCvJ6Y8SFCOIXD24clutp4S5HjijN3UXqqoQGM0deCHbCe/CjrWutiFAPHyL6NSJiIhpt6MYgzCJSV0UQ6qItBHktR6Vka4YLEe0F8FMAGjdK12MWqadlm3F8Fov1a/KC0GW0hSAPS/MWjX9DfhfA5wCftLON2MgsUitV/zOfkX6cbQYR3U9ErxDRSSJ6JOz91AMR7SWip4noJSI6QUQP6/FtRPQkEb2mH4fD3msgMPOWH/fccw/X4tAhYy0vPw4dqnnLpgGCnzNKAJjhGp8RVL2c39fPTwEYrXHdAQAzAGb27dtXvsBjjzGPjzMTqcfHHrPnxsf9P2ii8teZTPl9Ql2s99k2ckBVNH0dwM0AkgD+BcCdQczdygPATgDv0s/7AbwKVdTvtwA8oscfAfCbYe81iKPpzM7NUE9mJ2BzQVqVtdnt2aBEtAB/s8lBAL8B4KeYeZ6ITgGYZOZ1M+Lq/VwBKI273j++cZAKdRNUZudmMrfbESJ6HMCX9PFvmPk8Ee0E8PfMfEe4u2uetjCtNEozJhEJeyzjVWa+q/IA8AaA/QD+RQvxPQCeI6Idga3cSAihxI2HiV9RvN0h7WVTENEEgHcCeAbAdmY+r09dAFDd3ieCtK0gN4LVT+A24wSVsMeNYeYXmPkmZp5g5gmo/7zvYuYLgS3iZ0OvTBoySNy4sEmIqA/AXwH4FWa+4Z5jZY7oiN/kbSfIjcZcKawPHRKBG3ncKJWDB4GHHioPLfz0pyVuvP04C2Cv83qPHmt7iCgBJcSPMfM39PBFbVKBfrwU1v6CpC0FeaXGbAjaJCL1WepDa+bNVYzzi1J59FElpItF9fjEEyqJyNMtwyRuvB0oFcUjoiRUUbxvhbynDSEiAvAnAF5m5i86p74F4CH9/CEAj2/13lpCGB7W9aJWXExUSWXESrdHmwQFAops4Ho+11pRKiaaJZMpH5dolaYI8rMF8ABU1MfrAA4GNW8rDwDvhzKb/ADA8/p4AMAIgO8CeA2qBea2sPcaxNHWUSu1ytt2e7RJUGxZPXKgdpQKkbKB+9Url2iVTRP0Zyu0N21nWnGpZTYRk0gEWS/TU6ocCkJTtLUgr4U4PCPIepmeUuVQEJoikoK8EhHsEWC9AliNVDmspz6LIHQZHSHIpbhWRKhVAKveKoe16rOIMBe6nI4Q5EIHUE+VQ78a59L6TRCiK8gl1b4LEaeoIPgSaUFugo4BSbXvCsQpKgi+RFaQC12ItH4TBF86QpBLXHmXIK3fBMGXeNgbCAIxp3QRU1MiuAWhgkhr5CLAI47EhAtCIERakEv8eITZTEy4CH5B8CXSglyIMI3GhEsykCDUJHKCXOLHO4RGY8IlGUgQahJJQS7x4x1AozHhkgwkCDVpSpAT0W8T0Q+J6AdE9E0iGgpqY40igjxiNBoTLslAglCTZjXyJwHcxczvgOog8vnmt1Q/bvy4OD4jRqMx4Q880Ni4IHQRTcWRM/N3nJffA/CzzW2nMUQLjziNxIQ/8URj44LQRQRpI/8FAH9b6yQRHSCiGSKamZubC2RBcXx2ERvZyCU0UehiNhTkRPQUEb3oczzoXHMQQB5Azf89zHyUmSeZeXJsbCyQzYvjs4tYz0YuoYlCl7OhIGfm+5j5Lp/jcQAgok8B+DCAKQ6jk7PQHaznHJXQRKHLaTZq5X4AnwPwEWZe3uj6ViKFszqc9ZyjEpoodDnNFs36EoAUgCdJGaq/x8yfbnpXm0DMKV1ALefovn3KnOI3LghdQLNRK7cGtRFB2DRHjiibuGtekTrlQhcRucxOQahC6pQLXU5H1CMXBKlTLnQzopEL4SBx34IQGKKRC1uPifs2Nm0T9w2IVi0Im0A0cqE1rKdxb6YWuWjvglATEeRC8PhlWn7yk8BnPqPONxL3LVmbgrAhbS/IJT48HIjol3WJ4hNE9FsN3eyncTMDf/zHSgBv2+Z/n1/ct2RtCsKGtL0gl/K0Ww8R/SSABwH8GDO/DcDvNDRBLY2bGXj4YWBhofpcIuEf9y1Zm4KwIW0vyIVQ+CUAX2DmNQBg5ksN3b1eRuWVK0A2Wz0+MFA7a7PRNQShy2hLQS7laUPndgD/moieIaJ/IKJ/VetC3/LER47YD69erl71H2+0k5AgdCFtK8ilPG1rue+++wDgbTXKE8cBbAPwHgD/DcDXiPwls2954qkp4NOfrhbmmQwwMuK/oVoatmRtCsKGSBx5l/LUU0+BiE4w82TlOSL6JQDf0GWJv09ERQCjAOrvCPKHfwi8733KKXn6tBLURotutC6KZG0Kwrq0vSCX8rSh8L8B/CSAp4nodgBJAJcbnmU9AVwp4EVQC8KmaXtBLuaUUPgKgK8Q0YsAsgAeCrRpiGjYghAobS/Iha2HmbMAPhH2PgRBqI+2dHYKgiAI9SOCXBAEIeKIIBcEQYg4IsgFQRAiDgUZjFD3okRzAHy65TbEKDYTEtdehP0expl5LKjJmvxcw/5btIqw3legn63Q3oQiyIOAiGb8klmiRCe8h6Do1L9Fp74vob0Q04ogCELEEUEuCIIQcaIsyI+GvYEA6IT3EBSd+rfo1PcltBGRtZELgiAIiihr5IIgCAJEkAuCIESeyAlyIrqfiF4hopNE9EjY+9ksRHSKiF4goueJaCbs/bQDRDRNRGf13+R5Inog7D1tlk75dypEg0jZyInIA/AqgA8BOAPgWQAfZ+aXQt3YJiCiUwAmmbkTk2A2BRFNA1hk5saaPbcZnfTvVIgGUdPI3w3gJDO/oUutfhWq27sgtBPy71TYUqImyHcDeNN5fUaPRREG8B0iOk5EB8LeTBvxWSL6ARF9hYiGw97MJumkf6dCBIiaIO8k3s/M7wLw0wD+MxH9RNgb2gqI6Cmfhs+m6fMfAbgFwN0AzgP4H6FuVhAiQtQ6BJ0FsNd5vUePRQ5mPqsfLxHRN6F+jv9juLtqPcx8Xz3XEdH/BPA3Ld5Oq+iYf6dCNIiaRv4sgNuIaD8RJQF8DMC3Qt5TwxBRLxH1m+cAfgrAi+HuKnyIaKfz8t8jun+Tjvh3KkSHSGnkzJwnos8C+DYAD8BXmPlEyNvaDNsBfJOIAPUZ/C9m/j/hbqkt+C0iuhvKf3AKwH8Kdzubo4P+nQoRIVLhh4IgCEI1UTOtCIIgCBWIIBcEQYg4IsgFQRAijghyQRCEiCOCXBAEIeKIIBcEQYg4IsgFQRAizv8H+Yk5TDl5Y3EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_samples = 20  # nb samples (gaussian)\n",
    "n_noise = 20  # nb of samples (noise)\n",
    "\n",
    "mu = np.array([0, 0])\n",
    "cov = np.array([[1, 0], [0, 2]])\n",
    "\n",
    "xs = ot.datasets.make_2D_samples_gauss(n_samples, mu, cov)\n",
    "xs = np.append(xs, (np.random.rand(n_noise, 2) + 1) * 4).reshape((-1, 2))\n",
    "xt = ot.datasets.make_2D_samples_gauss(n_samples, mu, cov)\n",
    "xt = np.append(xt, (np.random.rand(n_noise, 2) + 1) * -3).reshape((-1, 2))\n",
    "\n",
    "M = sp.spatial.distance.cdist(xs, xt)\n",
    "\n",
    "fig = pl.figure()\n",
    "ax1 = fig.add_subplot(131)\n",
    "ax1.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "ax2 = fig.add_subplot(132)\n",
    "ax2.scatter(xt[:, 0], xt[:, 1], color='r')\n",
    "ax3 = fig.add_subplot(133)\n",
    "ax3.imshow(M)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute partial Wasserstein plans and distance\n",
    "----------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Partial Wasserstein distance (m = 0.5): 0.45151590745848863\n",
      "Entropic partial Wasserstein distance (m = 0.5): 0.46654948274375097\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEtCAYAAADHtl7HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfgklEQVR4nO3dfdxcZX3n8e9XEh4CgRDBbMJTFLHIoga4DT6gUKgWEF/gSxahysKubcTKKrvaBaG7JF1Q3EXpWlpqWDCpWB7kuRRZqIKUtgYSGkIgWAIESQgJzwEDLIHf/nGum07uzJl7rnk4M7nn83695nXPXOfhd53JnF9+c+Zc5zgiBAAAgOa9rdcdAAAA2NxQQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCqgBYPtM2/+nyXnn2T6n233a3Nj+vO1be90PYHNne3fbL9veotd9GYlc2b5BypUUUH3A9grbr6SksibtmNu1uK5DbK+sbYuIb0XE77fZx3GpfwfWtH3edtRpe6idWFWwPdv2Zc3OHxE/johPdrNPQKeMyCnDjwubXPYO223li0Yi4tcRsV1EvNGtGM0gVzaHXFmOAqp/fDoitpO0v6QhSX+cuwLb4zreqyQiNkj6J0kfr2n+uKSH6rTd2a1+NKub7wWwmfh0KlSGH6d2YqVjYd8iV/6rsfDv2SsUUH0mIlZJ+qmkfSXJ9n+wvcz2S7Yftf2l4XmHv0HZPt32U5IuT8tOq/nWOW3kNwjbP7H9lO0Xbd9p+9822b07tXEC+Jik79RpuzPFmWn7n2y/YHu17Qttb5mm2fYFttfaXmf7ftvD23yk7QfTNq+y/Y2avh9le3Fa5z/afn/NtBXpvVgi6Tfpm+DpaR0v2f6V7cNsHy7pTEmfS+/RfWn5HWxfkvq6yvY5wz8z2D7Z9l01scL2KbYfTn35c9tu8n0Eemb4s2z7fNvP237M9hFp2rkq9uELa49apc/7V2w/LOnh1PYR2/ekPHKP7Y/UxLjD9rdt35327xtsT07Tpqf1jUuvJ9v+oe0nU3+ub9Dvf0h55EXbD9k+rGY6uZJcWa2I4NHjh6QVkn4nPd9N0gOS/kd6/SlJe0qypIMlrZe0f5p2iKQNKnbMrSRtk9pWjlj/bEmX1bz+j5ImpmX+VNLimmnzJJ1T0s+DJT2novDeSdLjkiZIWlPTFpJ2T/MfIOlDksZJmi5pmaTT0rTflbRI0qS0be+VNDVNWy3pY+n5jjXbu5+ktZIOlLSFpJPSe7dVzfu4OL2H20j6LUlPSJqWpk+XtGe99yS1XSfpB5K2lfQOSXdL+lKadrKku2rmDUk3pf7vLulpSYf3+rPEg0fExjmlzrSTJb0u6Q/SfvRlSU9Kcpp+h6TfH7FMSLpN0uS0b02W9LykE9P+fUJ6/faadaxS8UVwW0nXDO9vaT8MSePS67+VdGXa18dLOrhBvzdI+s9pvs9JelHS5DSdXEmurPTBEaj+cb3tFyTdJekXkr4lSRHxtxHxSBR+IelWFd9chr0p6eyIeC0iXmkmUERcGhEvRcRrKnaOD9jeoYlFF6hIAu9LfbgrItZLeqymbUVE/DrFWRQRv4yIDRGxQsUOd3Ba1+sqEtPeKhL3sohYXTNtH9vbR8TzEXFvap8l6QcRsSAi3oiI+ZJeU5F4hn0/Ip5I78UbKhLfPrbHR8SKiHik3obZniLpSBVJ6zcRsVbSBZKOb/B+nBcRL6TtvV3SjNHfQqAy16dv/MOPP6iZ9nhEXBzFeUjzJU2VNGWU9X07Ip5L+9anJD0cET9K+/flKn6i+nTN/D+KiKUR8RtJ/03ScR5x4rjtqZKOkHRK2tdfT3muzFpJf5rmu1LSr1JfyJXkyspRQPWPYyJiUkTsERF/OLyD2z7C9i9tP5cKrCNVfHsZ9nREvNpsENtb2D7P9iO216n4JqIR66wrxblbxWHoj0v6+zTprpq2t37Tt/0e2zelQ+DrVBSFO6V1/VzShZL+XNJa23Ntb58W/Wzazsdt/8L2h1P7HpK+XvufgopvUNNquvlETX+XSzpNReJba/sK27Xz1tpDxbfa1TXr/oGKb1dlnqp5vl5SSyf+A10ynFOGHxfXTHvrs5v+Y5dG//w+UfN8moqjKrUel7RLyfyPq9i/RuaZ3SQ9FxHPjxJ72KqI4rBGzXqnSeRKcmX1KKD6mO2tVBz6Pl/SlIiYJOlmFYdxh8WIxUa+Hun3JB0t6Xck7aDiUK1GrLOR4d/2P6Z/TQp/X9NWe1LkRSq+le4VEdur+C39rTgR8f2IOEDSPpLeI+mPUvs9EXG0ih3yeklXpUWekHTuiP8UJqRvv2+ttrazEfHXEXGQip0+VBzC32S+tO7XJO1Us+7tI6LZcx6AsaIsh9S2P6lin6q1u4qf7YbtNmLa65KeGbHME5Im257UZN92GXH+zO6SniRXkit7gQKqv22p4rDq05I2uDjRc7ThoWskvb3BYeaJKj78z6o4xPytzD7dKem3VSTHB1PbP6g4n2CGNk4KEyWtk/Sy7b1VnGshSbL9QdsH2h4v6TeSXpX0pu0tXQzv3SEiXk/Lv5kWu1jSKWk5297W9qdsT6zXUdu/ZfvQlFxflfRKzbrWSJpu+22SlA6J3yrpu7a3t/0223vaPrjeuoExbI2kd40yz82S3mP799IJyJ9T8Z/7TTXzfMH2PrYnSPoTSVfHiEsXpP3up5L+wvaOtsfbrj3ReqR3SPpqmu/fqTgf6GaRK8mVPUAB1cci4iVJX1XxreJ5Fd+IbhxlmYdUjDB5NB1eHXkY9q9UHPZepWKn/mVmt/5RxbexBcOH0iPiGRWJa21EPFwz7zdSn19SsUNfWTNt+9T2fOrPs5L+V5p2oqQV6VD2KZI+n+IsVHHi64VpueUqTlgss5Wk81R8631KRfL9Zpr2k/T3WdvD5w38exWJ+MG0/qtVnBsCbI7+xhtfB+q6Jpf735KOdTEi7vv1ZoiIZyUdJenrKvbd/yrpqJQLhv1IxYnWT0naWkUuq+dEFUenHlJxjtNpDfq2QNJeKvbpcyUdGxHPkivJlb0wPOoCAICOsH2HipFbTV3Vu8l1nqxidOBBnVon0A6OQAEAAGSigAIAAMjET3gAAACZOAIFAACQqa0CyvbhLu6Zs9z2GZ3qFABUgRwGoFUt/4Tn4pL8/yLpE5JWSrpH0gkR8WD5MhOiuB0OqjZVq0unrR680aeo1OpnImLnXvdipNwcNhj5q+wakY2uHflm/ebtyy5kLWndk812SJJ0wLTy/LXoydz81ei4wRYl7a9nxsDYUZ6/xrWx1pmSlkfEo5Jk+woVV20tLaCK5DOrjZBo1SzNKZ02h38TdNWckbf86BeZOWwQ8tf4kvZG/1WU3FbuI7PLF7mlwbQ6Fn6lPH/5rNx/k20aTNu+pH1NZgyMHeX5q52f8HbRxvc6WqmN74MEAP2MHAagZe0cgWqK7Vl662tbMzexBoD+QP4CUKadI1CrtPHNInfVxjeSlCRFxNyIGIqIoeJ2QgDQF0bNYeQvAGXaKaDukbSX7Xfa3lLS8Rrl3kMA0EfIYQBa1vJPeBGxwfapkv6viqELl0bEAx3r2SjOLjkpeo7OrqoLmxXeF2BjleSwM2bXbf6Lb59ct/0P/f4GKyu72f2ykvblDdY1paT9uZL2DQ3WVSLzRPFGfFaj/PXRkvYVJe2b/FBSo+SEeKCOts6BioibJd3cob4AQKXIYQBaxZXIAQAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIFPLNxNuKZinxdi/l9Rg4/IS2NScRcWFKDdv5K9M+84un7a0wbS6pjeYtiJzXeXiC/Xzly8jfw2u8vzFESgAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADI1Na98HJN1WrNqjNKixFaYwf/lhirJh2wtQ5b+O5N2q9xo5v25iq7ya8krelgnDL/paT93gbL3FG/udFIu1NLpl1YtsyKBvHLjK/fvPVZpUv4shbCYGBxBAoAACATBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyVToKb7Wmas4Yv5cU94IDxqYXFr3a4RF3dez65fJpK79TMuGVzsX/xvb128cdUr7MeXfkx5mXv0i+1+s3H99gkXnd6AfGKo5AAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQKa2RuHZXiHpJUlvSNoQEUOd6BQAVIEcBqBVjojWFy6Sz1BEPNPc/NNCY/wyBq3g0gcY2+Ys6tfCJCeHkb/qu6Ekfx1dWf7apqS9g5d30P4l7Y1usoyxoTx/8RMeAABApnYLqJB0q+1FtvlqBmBzQw4D0JJ2r0R+UESssv0OSbfZfigi7qydISWllJh2aDMcAHRUwxxG/gJQpq0jUBGxKv1dK+k6STPrzDM3IoaK3xAntBMOADpqtBxG/gJQpuUCyva2ticOP5f0SUlLO9UxAOgmchiAdrTzE94USdfZHl7PX0fELR3p1RhUNtJOYrQd0COZOcySxtdpL7lpbSuOmV0+7frvl0x4rnPxTy6Jv1P5Ikef30qgkjil7Y1kjrY7vEGMW8qmMdoOm2q5gIqIRyV9oIN9AYDKkMMAtIPLGAAAAGSigAIAAMhEAQUAAJCJAgoAACBTuxfS7Dv9OtqNkXbA5i6UP+Juekn7ivrNJzdY1fFfLWmf3WxnRjdvff32f9Ppa2Bdkzl/2f3upOxReLc8lhkbqI8jUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACCTI6K6YJ4W0qzK4rWrXy+JAGxe5iyKiKFe96Jd5C9gEJXnL45AAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQKYxdzPhThqUkSplo3UGZfuB7im7Ce4RJe3X5of449l1m+ecM750kT3is3XbH/cV+fF77bLZ9du/UNIOdAhHoAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACDTqPfCs32ppKMkrY2IfVPbZElXSpouaYWk4yLi+VGDldxLins2AWNZb++F16kcZg+FdM8m7XFV+fdQH7e55a8pJe3rGizzSgtx/lNJ+5+1sK5c720wbVkF8bF5ae9eePMkHT6i7QxJP4uIvST9LL0GgH40T+QwAB02agEVEXdKem5E89GS5qfn8yUd0+F+AUBHkMMAdEOr50BNiYjV6flTKj/uCwD9iBwGoC1tn0QexUlUpSdS2Z5le6HthdL6dsMBQEc1ymEb56+nK+4ZgH7WagG1xvZUSUp/15bNGBFzI2KoOAlrQovhAKCjmsphG+evnSvtIID+1moBdaOkk9LzkyTd0JnuAEAlyGEA2tLMZQwul3SIpJ0krZF0tqTrJV0laXdJj6sYAjzyJM0666p/GYNe42a6QDf1/DIGHclhreWvyfWbb/pq/fajZmeuX4pf189f3r1R/tq+pL3R5Qr6VdlNk1+vtBcYq8rz17jRFo2IE0omHdZWnwCgAuQwAN3AlcgBAAAyUUABAABkooACAADIRAEFAACQadSTyMcSRtsBaN1EFYP5RlrSYJnP1W8+6sqS+Y9osK5312317g0WKVU2Qu19Lazr/haWKbuhbys38/103dZ947W67Uv96wbramVbMKg4AgUAAJCJAgoAACATBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyDdRlDLhcQWeUXQ5C4j3GGLb1RGn6IZu2P3RHg4UuyosxruSyB5K0oezSB614pX7zjM/Wb3+mwapWtjD0/90l27l8dv66dG3d1qUuucnw0Fnlq1rIZQzQPI5AAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQKaBGoWHPNx8Gajx6nrpoUWbth8zu3yZ6y8umbChpPmy3F61aHL95sXzSuYvu/mwJB1Y0r6gfJHlHRxReOrskhgl899S1XuMsY4jUAAAAJkooAAAADJRQAEAAGSigAIAAMg0agFl+1Lba20vrWmbbXuV7cXpcWR3uwkArSGHAeiGZkbhzZN0oaS/GtF+QUSc3/EeoSlVjJBjtB3GiHnqRA4bN0GadMCm7dff1UKX1pS0799gmXUl7WXDzRp5rqT9oyXtJffOkyT9XQvxy7ZzWf6qLjy3bvMOr36xbvuLW7+/wcpaeS8xqEY9AhURd6p8bwOAvkYOA9AN7ZwDdartJenw+I4d6xEAVIMcBqBlrRZQF0naU9IMSaslfbdsRtuzbC+0vVBa32I4AOiopnLYRvnrzaer7B+APtdSARURayLijYh4U9LFkmY2mHduRAxFxJA0odV+AkDHNJvDNspfb9u52k4C6GstFVC2p9a8/IykpWXzAkC/IYcBaNeoo/BsXy7pEEk72V4p6WxJh9ieISkkrZD0pS72EQBaRg4D0A2jFlARcUKd5ku60BeMUHapAolLDADN6lgO2/Ck9MzstvvT2E+7vP7R/E3H1vTihm+VTtsh9zb2B80un3ZX/Wkvbv2XmUEkfagkzi8bxMfA4krkAAAAmSigAAAAMlFAAQAAZKKAAgAAyEQBBQAAkCl3LAQq1MpIuypuMgwAo9lh3JkNph5St/Vm/Xbd9iNbuV9zKxhthwwcgQIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMfT8Kj1FleXhfgP5xVSyq236cDyhZYkrput4bB9ZtX+Ybc7vVwPiS9kb/VbzSQpw76rYeWZq/tulwfKB9HIECAADIRAEFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmfr+Mgb9OiyfyysAGE355QreV9J+f+m6yi9XcEhJ+x2l67om7q7b/lnPLFni9dJ1VaPBpQpWnlG/fdfzutMVIOEIFAAAQCYKKAAAgEwUUAAAAJkooAAAADKNWkDZ3s327bYftP2A7a+l9sm2b7P9cPq7Y/e7CwDNI38B6BZHROMZ7KmSpkbEvbYnSlok6RhJJ0t6LiLOs32GpB0j4vTG65oW0qzO9BzAZmLOoogY6kVk8heA9pTnr1GPQEXE6oi4Nz1/SdIySbtIOlrS/DTbfBVJCQD6BvkLQLdknQNle7qk/SQtkDQlIlanSU9JmtLRngFAB5G/AHRS0wWU7e0kXSPptIhYVzstit8B6/4WaHuW7YW2F0rr2+osALSC/AWg05oqoGyPV5F8fhwR16bmNen8guHzDNbWWzYi5kbEUPEb4oRO9BkAmkb+AtANzYzCs6RLJC2LiO/VTLpR0knp+UmSbuh89wCgdeQvAN3SzL3wPirpREn3216c2s6UdJ6kq2x/UdLjko7rThcBoGXkLwBdMWoBFRF3SXLJ5MM62x0A6BzyF4Bu4UrkAAAAmSigAAAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADIRAEFAACQiQIKAAAgEwUUAABAplELKNu72b7d9oO2H7D9tdQ+2/Yq24vT48judxcAmkf+AtAt45qYZ4Okr0fEvbYnSlpk+7Y07YKIOL973QOAtpC/AHTFqAVURKyWtDo9f8n2Mkm7dLtjANAu8heAbsk6B8r2dEn7SVqQmk61vcT2pbZ37HDfAKBjyF8AOqnpAsr2dpKukXRaRKyTdJGkPSXNUPEN77sly82yvdD2Qml9B7oMAHnIXwA6rakCyvZ4FcnnxxFxrSRFxJqIeCMi3pR0saSZ9ZaNiLkRMRQRQ9KETvUbAJpC/gLQDc2MwrOkSyQti4jv1bRPrZntM5KWdr57ANA68heAbmlmFN5HJZ0o6X7bi1PbmZJOsD1DUkhaIelLXekhALSO/AWgK5oZhXeXJNeZdHPnuwMAnUP+AtAtXIkcAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADINGoBZXtr23fbvs/2A7bnpPZ32l5ge7ntK21v2f3uAkAechiAbmjmCNRrkg6NiA9ImiHpcNsfkvQdSRdExLslPS/pi93rJgC0jBwGoONGLaCi8HJ6OT49QtKhkq5O7fMlHdOVHgJAG8hhALqhqXOgbG9he7GktZJuk/SIpBciYkOaZaWkXbrTRQBoDzkMQKc1VUBFxBsRMUPSrpJmStq72QC2Z9leaHuhtL7FbgJA61rNYeQvAGWyRuFFxAuSbpf0YUmTbI9Lk3aVtKpkmbkRMRQRQ9KEtjoLAO3IzWHkLwBlmhmFt7PtSen5NpI+IWmZiiR0bJrtJEk3dKuTANAqchiAbhg3+iyaKmm+7S1UFFxXRcRNth+UdIXtcyT9s6RLuthPAGgVOQxAx41aQEXEEkn71Wl/VMW5BADQt8hhALqBK5EDAABkooACAADIRAEFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMoxZQtre2fbft+2w/YHtOap9n+zHbi9NjRve7CwDNI38B6JZxTczzmqRDI+Jl2+Ml3WX7p2naH0XE1d3rHgC0hfwFoCtGLaAiIiS9nF6OT4/oZqcAoBPIXwC6palzoGxvYXuxpLWSbouIBWnSubaX2L7A9lZd6yUAtIj8BaAbmiqgIuKNiJghaVdJM23vK+mbkvaW9EFJkyWdXm9Z27NsL7S9UFrfoW4DQHPIXwC6IWsUXkS8IOl2SYdHxOoovCbph5JmliwzNyKGImJImtB+jwGgBeQvAJ3UzCi8nW1PSs+3kfQJSQ/ZnpraLOkYSUu72VEAyEX+AtAtzYzCmyppvu0tVBRcV0XETbZ/bntnSZa0WNIpXewnALSC/AWgK5oZhbdE0n512g/tSo8AoEPIXwC6hSuRAwAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADIRAEFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMjkiqgtmPy3p8fRyJ0nPVBZ8U4Mcf5C3vdfxB3Hb94iInSuO2XHkL+L3QexBj99X+avSAmqjwPbCiBjqSfABjz/I297r+IO87WNJr99H4rMPD2L8Xm/7SPyEBwAAkIkCCgAAIFMvC6i5PYw96PEHedt7HX+Qt30s6fX7SPzBjD3o8Xu97Rvp2TlQAAAAmyt+wgMAAMjUkwLK9uG2f2V7ue0zKo69wvb9thfbXlhBvEttr7W9tKZtsu3bbD+c/u5YcfzZtlel92Cx7SO7FHs327fbftD2A7a/ltor2f4G8ava/q1t3237vhR/Tmp/p+0F6fN/pe0tK4w9z/ZjNds+o9Oxx7pe5q8Uf2ByWC/zV4rVsxw2yPlrlPj9k8MiotKHpC0kPSLpXZK2lHSfpH0qjL9C0k4Vxvu4pP0lLa1p+5+SzkjPz5D0nYrjz5b0jQq2faqk/dPziZL+RdI+VW1/g/hVbb8lbZeej5e0QNKHJF0l6fjU/peSvlxh7HmSju32to/VR6/zV+rDwOSwXuavFKtnOWyQ89co8fsmh/XiCNRMScsj4tGI+H+SrpB0dA/6UYmIuFPScyOaj5Y0Pz2fL+mYiuNXIiJWR8S96flLkpZJ2kUVbX+D+JWIwsvp5fj0CEmHSro6tXdl+xvERnsGKn9Jvc1hvcxfKX7Pctgg569R4veNXhRQu0h6oub1SlX4oVDxD3Cr7UW2Z1UYt9aUiFidnj8laUoP+nCq7SXpEHnXfkIcZnu6pP1UfIuofPtHxJcq2n7bW9heLGmtpNtUHL14ISI2pFm69vkfGTsihrf93LTtF9jeqhuxx7Be5y+JHCZVnL+k3uawQcxf9eL3Ww4bxJPID4qI/SUdIekrtj/ey85EcXyy6qr6Ikl7SpohabWk73YzmO3tJF0j6bSIWFc7rYrtrxO/su2PiDciYoakXVUcvdi7W7FGi217X0nfTH34oKTJkk6vqj/omEHPYZXmL6m3OWxQ81e9+P2Ww3pRQK2StFvN611TWyUiYlX6u1bSdSo+FFVbY3uqJKW/a6sMHhFr0gfzTUkXq4vvge3xKnb+H0fEtam5su2vF7/K7R8WES9Iul3ShyVNsj0uTer6578m9uHpZ4GIiNck/VC9+fxvznqavyRyWNX7by9zGPlrk/h9lcN6UUDdI2mvdCb/lpKOl3RjFYFtb2t74vBzSZ+UtLTxUl1xo6ST0vOTJN1QZfDhHT/5jLr0Hti2pEskLYuI79VMqmT7y+JXuP07256Unm8j6RMqzmO4XdKxabaubH9J7Idqkr5VnLvQi8//5qxn+Usih0nV7b8pVs9y2CDnrwbx+yuHdeps9JyHpCNVjCh4RNJZFcZ9l4pRM/dJeqCK2JIuV3GY9XUVvxd/UdLbJf1M0sOS/k7S5Irj/0jS/ZKWqEgEU7sU+yAVh7aXSFqcHkdWtf0N4le1/e+X9M8pzlJJ/73mc3i3pOWSfiJpqwpj/zxt+1JJlymNcuGR9d72JH/VfHYGJof1Mn+l+D3LYYOcv0aJ3zc5jCuRAwAAZBrEk8gBAADaQgEFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZPr/qMEntv/9IFkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = ot.unif(n_samples + n_noise)\n",
    "q = ot.unif(n_samples + n_noise)\n",
    "\n",
    "w0, log0 = ot.partial.partial_wasserstein(p, q, M, m=0.5, log=True)\n",
    "w, log = ot.partial.entropic_partial_wasserstein(p, q, M, reg=0.1, m=0.5,\n",
    "                                                 log=True)\n",
    "\n",
    "print('Partial Wasserstein distance (m = 0.5): ' + str(log0['partial_w_dist']))\n",
    "print('Entropic partial Wasserstein distance (m = 0.5): ' +\n",
    "      str(log['partial_w_dist']))\n",
    "\n",
    "pl.figure(1, (10, 5))\n",
    "pl.subplot(1, 2, 1)\n",
    "pl.imshow(w0, cmap='jet')\n",
    "pl.title('Partial Wasserstein')\n",
    "pl.subplot(1, 2, 2)\n",
    "pl.imshow(w, cmap='jet')\n",
    "pl.title('Entropic partial Wasserstein')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample one 2D and 3D Gaussian distributions and plot them\n",
    "---------------------------------------------------------\n",
    "\n",
    "The Gromov-Wasserstein distance allows to compute distances with samples that\n",
    "do not belong to the same metric space. For demonstration purpose, we sample\n",
    "two Gaussian distributions in 2- and 3-dimensional spaces.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eXRkZ3nn/3lr076vpaW09Cr1ot3NCWAclkAgx4Tl2CRmsDEzTkI84WeTCeacIW6f82OOQ0jGWZgJTiDDSRjnhPmROCcDhAB2gBm7u7V1t6TWLrXW0q6SSqr9vr8/1Ldckqqk2qSWut/POX2Qqu597y358L1PPe/zfB8hpUShUCgUxxfD3b4BhUKhUCSGEnKFQqE45ighVygUimOOEnKFQqE45ighVygUimOO6W5ctLCwUFZXV9+NSyvuAzo6OhallEV36fKqDExxkIhwL94VIa+urqa9vf1uXFpxHyCEuH2370GhOExUakWhUCiOOUrIFQqF4piTFCEXQjwjhOgVQvQIIV4RQqQmY12FQqFQ7E/CQi6EKAd+B2iVUp4HjMAnEl1XoVAoFNGRrNSKCUgTQpiAdGAmSesqFAqFYh8SFnIp5TTwVWACmAUcUsofJrqu4t7n8uW7fQcKxb1BMlIrecCHgRqgDMgQQnwyzHFPCSHahRDtCwsLiV5WcQ/wwgt3+w4UinuDZKRW3guMSSkXpJQ+4LvAL+w8SEr5spSyVUrZWlR0t3o1FAqF4t4jGUI+AbxNCJEuhBDAe4BbSVhXcQ9y+TIIsfUP3vpZpVkUivgRyRgsIYR4AXgU8ANdwL+XUnoiHd/a2ipVZ6dCCDiIuSZCiA4pZWvyV44K1aKvOEjCtugnpWpFSvm8lPKslPK8lPLf7SXiiuOPip4ViqOF6uxUxEyyNimffz456ygU9zt3xTRLoQAV2d8raJqGy+VCCIHFYsFoNCJE2AyA4oBQEbkiKtQmpSIcUkq8Xi+BQACv14vT6WRtbQ23200gEEANdz8ckrLZGStqs/N4c1CblMlCbXYeDrqIa5qG3+9HSokQAiklgUAAAJPJhMViwWw2YzCouDEJHB0/coVCcbyRUjIxMYHFYiE/P3/be0IITCYTUspg2sXlcmEymUhJScFkMqnUS5JRQq6IGbVJeX+jR+JOp5PMzMyIxwkhgoKtR+kbGxt4vV4MBgO5ubkqn54klJArYkblxe9fQtMpehrF4/HQ3d1NWloaVquV7OzsXeIcKup6Hl2PzC0WCxaLBYPBoEQ9TpSQKxSKqNgp4gaDAb/fT2dnJ/oM3qmpKTY2NigqKqK0tJS0tLRd6+iibjQagw8Cj8eDwWDYJuqK6FFCrlAo9mWniAsh0DSNsbExzp49S25uLlJKCgoK8Pv9LCwsMDAwgKZplJSUUFxcjNls3rVuqKBLKXG73bjdboxGo9okjQEl5AqFYk+klPh8vm0iHggEmJ6epqCggOLiYrxeb/B4k8mE1WrFarXi8XiYm5vj+vXrpKSkUFpaGrYkcWc+XdM0Njc3EUJgNpuxWCxqk3QPlJArFIqI6CIeCAS2ReLd3d1kZ2eTm5sLEMyX7yQlJQWbzYbNZsPpdGK325mfn8dsNlNQULBnPt1gMASvr2+Q6qKuNkm3o4RcoVCERUqJ3+/fJeLXr1+nsLAweEy0ZGZmcvLkSXJycpifn2d6epqBgQEKCwspLS0lPT191zl66kW/ltfrxePxsLKyQllZGWazOfj+/YwScoVCsQtdxP1+f1DEpZT09PSQk5NDVVUVExMTca0thCAtLY3a2loCgQALCwsMDg4SCAQoKSmhpKRkz3w6wNjYGHl5ecH69Ps9n66EXKFQbCOSiPf19ZGamkptbe22YxPBaDRSWlpKaWnptny6xWKhtLSUwsLCsOKsmo62o4RcoVAEiSTiAwMDGAwGTp06FTw2Ul48XsLl08fGxsjJyaG0tJScnJw969NDm46AYCnj/ZBPV0KuUCiC7BRxgJGREfx+P+fOnTs0QdTz6SdOnGBlZYWZmRkGBweD+fRw7BR1r9eL1+u9L5qOlJArFAogvIiPjY3hdDppaGgIGw1rmhbzdWIRUiEE+fn55OfnB/PpQ0NDbG5uMjU1FVU+/X5oOrp3PonivkDZAxwMfr8fn8+3TcQnJiZYWVnh4sWLYcU3kdRKPOfp+fSGhgZSU1ODFTQ3btxgfn4+4kNFF3W9nNHtduNwOJiamsLj8cT1MDpqKCFXHCuSNZ1I8RbhRHx6epq5uTkaGhqOZORqMBiw2Wy0trZSW1vL+vo67e3t9Pf3s7q6GrHpyGAwBDtJR0ZG2NzcZG1tjY2NDXw+37H1T1epFYXiPiaciNvtdqampmhpadmzRjvZm53xkpmZSWZmJrW1tayuru7Kp4erT4ftlS/Hveno6D1qFYodqOlEB0MgENgl4gsLC4yPj9Pc3IzJdLziPCEEeXl51NfX09LSQnp6OkNDQ3R0dDA1NbXNRkAfgqGfZzQagyWLXq+X9fV11tfXg5OOjjrH67+U4r7k8uW3RPuoTyc6Lvh8PmZnZykuLg4K2vLyMsPDw7S0tITdQNzJUYnIwxFan+71epmbm+PGjRvB+vSsrKyIeX/9W8jO+vSj3HSkhFyhuM/QBer27duUlJQAsLq6Sn9/P83NzVgslqjWOS5pB4vFQmVlJZWVlWxsbGC32xkdHSUQCLCyskJubm7Yz2IwGIIbpEe96UgJueJYoaYTJYamacEyPD2aXltbo7e3l+bmZlJTU2Na76hG5JHIyMjgxIkTWK1Wbt26hd1uZ2hoiIKCAkpLS8nIyNh1zl5NR5qmkZ2dfdfz6UfvO4JCsQd6ikXlx2NH07Rgg4zRaETTNJxOJzdv3qSxsTHsEIi9OMqplWhISUmhrq6OlpYWMjMzGR4epqOjg8nJyW359FBCK18MBgNdXV3BiUd6Pv1u/E2UkCuOJaoMMTZ0EYe3IsxAIMD169e5ePFi2Eh0P+KNQI9COiJ0s9NoNFJSUkJDQwMXLlwA4MaNG1y/fp25ubmIm52hk46EELjd7m2bpIdZn65SKwrFPY4u4lLK4Eadx+Nhc3OTS5cukZWVFffah9kQlExChTyUcPn027dvk5WVRWlpacR8+n6TjiwWy4E+wFRErjg2qDLE2NE9R0JF3Ov10tnZSWpqKjk5OXGvHSpMc3NzLCwsHJsuyUhCHoqeT29ra6O0tBS73c61a9cYGRkJ5sh3sjP1ok86cjgcbGxssLi4GDFtkwhKyBXHhsuXt0oP9WBO/1kJeXjCibjP56Ojo4PTp08npU5cSsnk5CR2ux2Hw8G1a9cYGhpifX094bUPkmiEXEevT9+ZT29vbw8ac0U6z2AwYDKZMBgM+Hw+Pv/5z9PV1ZXMjwKo1IpCcU8SOixZF3F94v2JEycoLCxkcHAwoWsIIXA6nTgcDi5evBi87vLyMuPj43g8nuCgiGhLGg+TeFIdej69pKQEl8tFV1cXN2/exGQyBf3Tw3XD6qkXr9dLSkpKMm5/G0kRciFELvBXwHlAAk9KKd9IxtoKRThUGWJkwol4IBCgq6sLm81GcXFxUq6ztraGw+Hgne98Z/AaBoOBwsJCCgsLtzXipKSkYLVayc/PT8q1EyWWiDwSZrOZ1NRUmpub2djYYG5ubt98usfjibnEMxqSFZH/CfADKeXHhRAWILy5gUKRJFQ6JTyhIq6LiD4sWZ9snwx098Di4mJMJhN+v3/XMaEbh06nk9nZWUZGRkhPT7/rczaTIeShD8qMjAxqa2upqanB4XBErE8/skIuhMgBHgSeAJBSeoHkZ/MVCsWe7BTx0GHJBQUFVFRUJOU6TqeTnp4eTp06xerqalTnZGZmcurUKTRN4/bt29jtdtrb2yktLY3oKX6QJKNqJtzDQAhBbm4uubm5BAIBlpaWGBkZwev1YrFY8Pl8+6ZWJicn+dSnPsXc3BxCCJ566ik+97nPsby8TEFBwb8C1cA48IiUcgWSs9lZAywAfy2E6BJC/JUQYldRqhDiKSFEuxCifWFhIQmXVSgUOuFEXB+WnJ2dTXV1dVKu43a7g7Xn8eR6DQYDOTk5FBUVcfHiRaSUXL9+nZs3b7K4uHhoVS/JjsjDYTQaKS4u5uLFi1y8eJGpqSl6e3t5/PHH6evri3ieyWTij/7oj+jr6+PNN9/ka1/7Gn19fbz44osAP5ZSngJ+DDynn5MMITcBzcB/l1I2ARuhF9CRUr4spWyVUrYWFRUl4bKKg0ClLI4fug3rThG/desWqampnDhxIinX0csW6+vrg6ZT8US2uoDqqZfW1laqq6tZXl7m2rVrDA8PRyzvSxbJEPJY1rBYLHzwgx+koqKCr371q3vuU1itVpqbmwHIysqirq6O6elpXn31VYBv3TnsW8Cv6uckI0c+BUxJKa/c+f1/EUbIFceDF15QYn6c0EU8EAhsE3G9IiV0WHIi+P1+urq6OHXqFHl5eRHvJVph2/kAyMrKIisrC03TgukIn88XrBBJduolWRF5rGt4vV7OnTsX9R7B+Pg4XV1dXLp0ibm5OaSUs3fesgMl+nEJC7mU0i6EmBRCnJFSDgDvASJ/b1AoFElBn3gfKuIAo6OjeL1ezp8/n5RuQn2ztLKyktBv0wfhtWIwGCgqKqKoqAiv14vdbqe7u5u0tLRg1UsyPlOyIvJ4LG2jPcfpdPKxj32Ml156iezs7J3XlkKI4B8/WQ1B/xH4thDiBtAI/Jckras4BFTH5PFDF/Gdw5LHx8dZX19PmohLKblx4waFhYWUlZVte++gTbMsFktwnJvNZmNpaYmrV68yMjKS8LCHuynk0eDz+fjYxz7GY489xkc/+lEASkpKEEJYAf1/5/Xjk3IXUsruO/nvi1LKX9V3UhXHg1g7JpXA310iifjk5CTLy8sRhyVHWmuv927dukVGRkbSNkvjQQhBdnY2p0+fpq2tjczMTDweDx0dHUxPT+Pz+WJe826lVqJ58Ekp+cxnPkNdXR3PPvts8PWHH34Y4PE7vz4OvKq/p1r0FTGjnAfvHroh084RbTMzM9jt9piGJe8XUY+MjCCl5OTJk3GdfxAYDAaKi4tJT0/n/Pnz+P1+uru76e3tZWlpKer7uRtCHu29/Z//83/4m7/5G37yk5/Q2NhIY2Mj3/ve93juuecA3ieEGALeC7yon6Na9BXbUB2TRxu/38+NGzc4efIkmZmZwJZh1eTk5L7DkneylxDfvn0bp9NJQ0NDRLG6W3a0uginpKRQVVWFzWZjfX2d2dlZZl57jYqxMXJKSjB8+MMQwRTsbqZW9rvuO97xjoj/XaSU7wn3uorIFcD+AxsOO4+u0je70dMpodN9FhcXGRsbi2tYciQhn5mZYWFhIaoUzd22o4W3Ui91DgeX/st/wfrNb5Ly5S8T+OAHsd+6FbbrdC8hN77xBqlPPEHqJz+J8cc/jnjdeB4GB/XwU0KuAPZPlxy286BK32zH7/cH0ym6Pery8jKDg4M0NzfHVZ6nrxPKwsICExMTNDY27htt3s2IPByWr3wFhMBQXIy5tJSs1VXSv/c9urq66O3tZXl5OXhuJBE2XLlC6hNPYPzZzzC++Sapv/mbGH/0o7DX268hKNr7TgYqtaJQHHFCRVwX8rW1NSYnJ2MalryTnRH56uoqQ0NDtLa2Rh3dJ9IQlAhhhzusrUHIA00IQYHJRGtrK2tra8zOzjI8PExhYSGWiQlyb9zAlJ2NdukSWl0dAOa/+Rvw+eBO2orNTczf+AaB97531/Vijch9Pt+BWRGoiPw+Jt50yUHl0VUZ5G52irj+2ujoKE1NTQkZMIUKudPppLe3l6ampqgfDIlsdh5EdOp///thY2NLiF0upMlE4B3vQAhBTk4OZ8+epaWlhezlZVK//W3WR0ZYHx/H8O1vI3RL33DCHEGsY43I3W73gVjYghLy+5p40yUHmRdXgyPeIhAI7BJxfcpMbW1tzMOSd6ILscvl4vr16zQ0NMS0Zmg0ephplkiRsO93fgffJz8JRiMyLw/PH/wBWkvLtmOMRiMldjvpRUVkVlbiy8hgZn2d2f/9v1lZWcH7+ONbUb3TufVQMBjw/cZvxHQfkfB4PAcm5Cq1olAcQUIn3utisbm5SXd3N8XFxUkRBCEEHo+Hvr4+zp07F6yCiYWjsNkZxGzG94Uv4PvCF/Y+zmJB+HwYTSZysrPB78dZUsLY3BxDmkblH/4hlf/8zxgB/6c+ReAd7wi7TKwRuRJyxYFz1MoOj9r9HCaapuHxeLaJuNvtpru7m/PnzyfNJVBKSW9vL2fOnCE3Nzfm8+/WZmeiaC0tyNdfxzwzs5VXNxpJffe7OVtaSiAQYGFhgWuVlUgpsVqtFPn9YfcMYo3ID2o6ECghV9zhqKUvjtr9HBbhRFx3HTx79iw5OTksLy8nLOSBQIDV1VVqamooLCyMex09Ip+dnWVzcxOr1XpgYhV6zUQeIrKggKVHHyV7YoLU3Fy0ujrknb+B0WiktLSU0tJSXC4Xdrudzs5OMjMzsVqt26b+xNoQdJARucqRK3Zxv4roUWCnAZbP56Ozs5PTp08Hx6SFKxuMBd0/JT09nYKCgrjX0XPsdrudmZkZTCYTPT093Lhx41C9xePBn52N921vI/DOdwZFfCdpaWnU1NTQ1tZGWVkZdruda9euMTY2hsvlirkh6CA3O1VErtiFsrK9e4SKtD4seWfUrE/+iQc9nZKdnR0Uo0Tw+XyMjY3R1NQEQEVFBU6nk5mZGUZGRigsLMRqtZKe/tb0x8P0AU/GGjun/szPz9Pf34/L5aKoqIiCgoKoOmpVakWhuM/QhyVXVlZSUlKy7T2DwRC3+9/Q0BBGo5Ha2lr6+voSEvLNzU3W19d55zvfidlsDppXZWZmcvr0aQKBAIuLiwwMDABbAxOOylCZeB8GRqMxOPt0YGAAr9dLR0cHWVlZWK1WcnJyIq6rUiuKA0fVcB8ddP/v0tLSXdaxEH9qZXx8HJfLxdmzZ7cNoYgHt9tNX18fmZmZEWvZjUYjJSUlNDU1cfbsWTY3N2lvb2diYgKvN7GxvkchqtdFva2tjdLSUmZmZrh27Rrj4+O43e5dx6uqFcWBc/nyW6ItxFu13IrD58aNGxQUFFBZWRn2/XiEfHp6mqWlJZqamoICFu8DIXRa0MTERFTnpKWlBafMT05OMjU1FRy+XFpaGpNPzEENTo4VvfxQCEFeXh55eXn4/X4WFhbo6+tDCBH8FmI0GlVErlDcL8zPz5OVlbWn/3esAjw/P8/U1NQui9t4InJN0+jq6qK6upr8/PyYz9e7LAsKCrh48WJwvb6+PlZWVg6tLv2g8uwmkyk4c/PMmTNsbm7S0dFBf38/t2/f3rdr9sknn6S4uJjz588HX7t8+TLl5eU0NjYihOgWQnxw53kqIlfs4n6u4b7blJSUBKtTIhHLZufKygojIyO0tLTsinpjFXIpJT09PRQVFWG1WsO6CsaCPgGosrISh8PB7OwsQ0NDlJSUUFpaGjF6jUuENzcx/cM/YBgfR6utRdTXJy0ij0R6enrwW8ji4iLf//73GR0dpbKykk9/+tNhz3niiSd4+umn+dSnPrXt9WeeeYbf/d3fha0JbLtQQq7YhcqLH21CbWz3Yn19nb6+PlpaWsJGgrEK+eDgIBaLZdu3hWRE0KFVIX6/n7m5OXp6ejCbzZSVlZGfn5/YSLVAAMsf/zGGvj5kdjam3l5Krl5NOGKJ9oEihKCoqIiPf/zjuN1uWltbIx774IMPMj4+HvO9qNSKAlDifZyIJrWyubnJjRs3aGxsjLgZGYuQT0xM4HK5OHPmzLbzk43JZKK8vJyWlhZqa2tZXl7m2rVrjIyMsLm5CcT+8BDz8xj6+5GVlZCbi6ysJG1kBOPiYkL3Gk+Lfk5ODhcuXIj5Wn/+53+u+8N/UwiRt/N9JeQKQPl/HxWiEcf9hNzj8dDd3c2FCxfIyMjY81rRiOLc3Bxzc3NhB00cZE5bL2NsbW0lMzOTgYEBurq6WFhYiG2hCLv3IsHByfG06MfjVvlbv/VbjIyM0N3dDTAL/NHOY1RqRaE4Zuwl5Hon6JkzZ8jOzo57HZ2VlRVGR0dpbW3dFX0e1sxOvYyxpKQEl8vFxMQEKysrDA4OYrVaycrK2vN8WVJC4IEHML7xBmRkwMYG6+fPY0mwpv2wbGx39BH8JfDPO49REfl9jKodP55E2uzUm4hqa2ujar3fT4idTid9fX00NTWFHYgQb2olkZRMWloalZWV5OfnU1BQwPj4OO3t7UxPT0fefBUC32c/i++JJwi0teH79KexP/rosYnIZ2dnQ3/9CNCz8xgVkd/HqNrx40m4zU5N07h+/TplZWW7OkEjsZeQu91url+/zsWLF/cUn7tlYyuEoKCggIKCArxe7zZzq7Kyst0dlmYzgQ98AL0fVuvpOfCqlZ1EU0f+a7/2a7z++ussLi5SUVHBCy+8wOuvv053d7d+v78I7DJIV0KuUBwzdqZEdP+U3NxcKioqol4nkpDrDT9nz57dM22xc7BErKWM8bLz3J1ljDMzMwwODu5ZxnjYfi0QndfKK6+8suu1z3zmM6G/PhzuPCXkCkDVjh8ndgr5wMAAZrOZmpqamNYJJ766PUB1dXXMzoiH6U8edmZnDGWMyezsjBbVoq84cFRe/GgQa9XK6OgoPp+P8+fPxyxMQoht5lt6w09BQQFWqzW2Gz9i6GWM5eXlYd0YkyHkENvDS7XoKxSKILqQT01Nsbq6yrlz5+ISpZ259qGhobgi+8MmVhHOlJK6xUUuBQJkmc0MDAzgcDhYWFiI20UyHjweT0LDsvdCReQKxRFjv3yzEAKv18vMzAwtLS1xdz2GXmdiYoLNzU0aGhriWuuoIubmSH36aVhcRAC26mqK//RPaR8YwOVy0d7eTl5eXlRljIkSb9VKNCghVyiOGcvLy3i9Xn7hF34hqoEGkdCFfG5uDrvdTktLy7GYwxnLRqn55ZdhcRFZUoIEDKOjmP7X/8LQ0EBNTQ0nT55kaWmJ8fFxPB4PVquVkpKSmNwYo0XlyBUKBQAOh4OBgQHS0tLC1nbHghCCzc1N5ufnaW1tTeihoGkafr9/X3e/ZBH1dJ/ZWWTIdCJpsWCYnUXe6VIVQlBYWEhhYWGwjLGrq4uMjIzwZYwJcCxy5EIIoxCiSwixq+tIoVAkzsbGBj09PTQ2NiZmInUHj8fDwsICjY2NCT0U9EoX3Y52dXV139RQsuvPTf/0T6S9972kv+MdWP7f/xc8HgACbW2I9XXQNPD7ER4PgaamsHl2vYyxtbWVsrKy4KCI2+PjeNbXE77Hg8yRJ3Oz83PArSSup1Ao7hDaoBM6/zJePB4PY2NjFBUVkZaWFvc6Ukr6+/vJzc3dJoDt7e1MTk4Gx78lk50ibLh2DcuXvwxCILOyMP3TP2H+2tcA8D/2GP5f/mUMc3MYlpbwffKTBD7wgYgbpsbXXiPjwQcpe897aPje92jz+6n8/d9He+wxlp95hqXR0bjnpR75HLkQogL4EPBl4NlkrKlQ3K/sjFi9Xi9dXV3U1dUlZUNOb/ix2Wx47kSu8eL1egkEAtTU1ODz+YJ13D6fD7vdTnd3NxkZGZSXl5OdnX0gOXjjtWtbP9wRSZmTg+lnP8P37LNgseD74hfxff7zYDDAHrlvQ3s7aZ/4BMLl2vr9D/8Qc3U1gXe8gzSbjazbt1n6q7/i2oc/TGFhYcyCHggEDiT3DsmLyF8Cfg+I+MmEEE8JIdqFEO0xu5cpFPcpuuieOHGCvLxd7qUxo6dBqqqqyMnJSSjFMTc3h9/vD1v+aDabqaysDEbp+mi3qamphAdSwPYcuSwo2OYvIdxu5M6/lcWyS8R33fN3vhMUcQDhcmEYH4f0dBACo81Gid1Oa0sLmZmZuN1uurq6sNvtUZUxJqt2PRwJC7kQ4leAeSllx17HSSlfllK2Silbj8okbYXiKKP7p1RUVFBcXJzwenorv97wk0iu2uFwMDIyQnp6+p75er3b8ty5czQ2NqJpGgMDA6yuruJwOOK6/s5z/B/60NbUn4UFxMICGI14t6bpvHUfMzOYX34Z89e+hqG/P/y6qanInZ9FiK38OsD6OrKgAKPJRHFxMenp6duGSg8ODrKehFx6PCQjzn878PCdOXKpQLYQ4m+llJ9MwtoKxX2JlJKbN29SUFBAeXl5xGNiifCGh4cxmUzBCT/xCrnL5aKnp4empibdIzsqzGYzNpuN3NxcRkdHmZqawuVyxTWAeRuZmbi/+U2MP//51mZmSwuyrCz4tpiaIvVTn0I4HFv38bd/S9bTT0Nb27ZlfE8+ieUb30A6nQhNQ6al4fvIRxCTk2A0Is1mfE8/DbzVnh86VHqvMsaDNhdLWMillF8EvggghHgI+F0l4gpF/OgbiGlpaRGHMOsiHK2QT0xMsLGxQUNDQ/CceITc7/fT3d3NuXPn4t50FUKQkpJCXV3dNufCrKwsysvLycrK2vNzhf3c6ekEfumXwh5v+s53EA4H8s63GrG6SsU//AM89tj2dauq2PjpT7F87WsIpxPfI48QePe7EcPDiI0NtMpKuOM/s/MedpYxzs7Obitj1L3h9/vv9eSTT/LP//zPFBcX09Oz5Va7vLzMo48+yvj4OMPDw/8KPCKlXAk9T7XoKxRHjJGREQKBAKdOnYp4TDRDIXTm5+ex2+1cuHBhe7VHDGvAW/n1mpoacnNzoz5vL/SSv7a2NkpLS5mYmKCjo2Nvf/EYEZubEFIjL00mDG532GNlbS2eP/oj3F//OoH3vGerEubECbT6+qCIw96GWRaLhaqqqm1VPF/60pfQNI25ubk97/WJJ57gBz/4wbbXXnzxRd7znvcwNDQE8GPguZ3nJVXIpZSvSyl/JZlrKhT3GwUFBfv6p0QrwqurqwwPD9PU1LSr4SeWiFxKSV9fH/n5+ZSWloY9JpGNPCEEeXl5nD9/nosXL+L3++ns7OTWrVusra3tupdYruXXI/X1ddjcRLhcLLzznVGda/zZz0j9zd8k9T/8B8wvvQQbG1Hfg74/UF9fz7PPPsvm5iaPPvooG3fWCMeDDz5Ifn7+ttdefUplpu4AACAASURBVPVVHn/8cf3XbwG/uvM81dmpUBwx8vPz962CiEbINzY26O3tpbm5OeKEn2iFXJ/sngxDrf2uq0e0NpuNlZUVJiYmcLvdwbxzrGhtbXi+8hXML7+M8PnwfvzjzNfUULXffQ4NYf6rv0IrKYGUFIxdXfDtb+N76qmYLWxzcnLIzs7m9ddfj/n+5+bmQt0o7cCuP4IScoXiGBLNAObr169z4cKFiA0/0Qq53W5naWmJ5ubmQ/ccz8/PJz8/P5h37uzsJCUlJeaN0cBDDxF46KG3XtBrz/fAMDmJFCJYn66VlmK4k7eO9VtBstrzpZRSCLHrP5rKkSsUx5C9RFivPd9vAHM0Qr66usrY2BgNDQ1JsQWIFz1Kb2trIzc3l7W1Ndrb25mZmUlaLn0nMjsboWnBGnWxvo68UzqtadqhCXlJSUlwbqcQwgrM7zxGCblCcQyJFJHrtec2m23fCT/7Cfnm5ia9vb0Je7EkEyEEWVlZFBUVceHCBbxeL52dnQwMDCS9hltraiLQ2ophYgIxNbU1xPlOrlpKeWjTgR5++GG+9a1v6b8+Dry68xiVWlEojhixTgnS0Rt+8vPzKQupo45lDR2fz0d3dzfnz59PyIvlINAfPikpKVRXV1NVVcXy8jLj4+N4vd5gLj0RN0cAjEZ8Tz9NYHgY3G40mw3uVOsclJCHG7783HPP8cgjj/CNb3wD4L3AIzvPU0KuUBxDwonw8PAwRqMxYu35TiJF5HqZ4YkTJ8jJyYn6ng6yBX0nO2u4CwoKKCgowOPxMDs7S0dHBzk5OZSXl5OZmRn/hQwGtNOnd718UKmVcMOXAX784x/rP7433PtKyBWKY4gQYpuQT05O4nQ6aWxsjN6rO4yQ61F9UVFRXBUid5udUbo+07SsrIzi4uLEo/Q7xLPZeZBe7UrIFYpjSOi8zfn5eWZnZ2Oe8BNOyEdHRzEajVRV7VecF369wyDaGu7QKH1mZoaOjg5yc3OTMlg61vLDg7SwBbXZqVAcS/TUit7w09jYGHO0uVPIZ2dnWV1d5ezZszGJ8mGmVOIhJSWFmpoa2trayM/PZ3R0lM3NTWZnZ+MevhzrZ3a73Qc2HQiUkCsUR45oNztdLhe9vb00NTXF9bU9VMhXVlYYHx+PucwwXgFPxEQq3nN1PxS9tt7tdtPR0cHg4OCe3ZbhiCciP0ghV6kVheIYomkaY2NjNDc3x11VoovwxsYGfX19tLS0xOVAqEenmqZF9a0gGdF7ImtIKTEajdTU1ARz6cPDwwQCAaxWa1S59KMWkSshVyiOGX6/n5mZGaxW654NP9Gg151fuHAhrhyuHtWPjo4yOTlJQUEBFRUViVWKHDChImwwGIKuhW63OzimTi/hzMjICLuGisgVCkXc6MKbl5eX8OxOTdNwuVw0NDTE/UAQQjA3N4fD4eBtb3sbq6urjI6O4vf7g5Uiye4ITTQnH+n81NRUamtrqa6uZmlpKRill5WVUVRUtC1Kj6dq5SA3O5WQKxTHBN2BMC8vD4vFEvdGnb5WT08PZrM5oelDfr+f8fFx2tra0DQtWCnidruZnp7m2rVrFBYWUlZWdmQai/YTYYPBQFFREUVFRbhcLmZnZ3dF6SoiVygUexJJZEZGRjAYDNTU1DA7O5vQhPqRkRHMZnNCtc1ut5vNzU0uXbqEyWTC6/UG30tNTeXEiRPU1NSwsLDArVu3MJlMlJeXH6igRUMs0bQ+ASg0Stc0jZSUlJiapZJlmhUJJeQKxTFgamqK9fX1YMNPrEMhQpmenmZtbY2mpibeeOONuNYIBAJ0d3eTkZGxZ6RtMBgoKSmhpKQEp9PJ9PQ0KysrGAwGfD5fXB4uB5Va2YudUfqtW7dYWVnB5XJRVla2b5rL7XYnvJ+xF0rIFYojzsLCAtPT07S2tiY0pg22xoZNTk5uWytW9HmiFRUVzM7ORn0fmZmZnDlzhrW1NQYGBoIPgoqKin3Hu+28fiIk+iBIS0sjPz+ftLQ0hBAMDg4ipQzm0sOlXFRqRaG4j3E4HAwNDdHa2rptsy2eiNzpdHLr1q24ywx1RkZGSE1NpaKiArvdHvP5RqOR9PR06uvrcTgcwcERZWVlUZtdHXZEvhO91LKwsJDi4mJcLhczMzOMj49TUFCwK0pXQq5Q3Kdsbm4Gp9XvzGXHKuRer5cbN25w8eLFhKonZmdncTgcNDc3B1+LJ0LWxTQ3N5fc3NxtZld5eXmUl5cnXJWz37WTuUZaWlpwT2BxcZGBgQEAysvLKSwsjKtqpbq6mqysLIxGIyaTifb29ojHKiFXKI4gXq+X7u5uLly4EFbQYhHyQCBAV1cXp06dIisrK+57cjgcwQqVRFI84URUN7uy2WwsLS0F0xUVFRUUFBRsS1ckI0eeKJHuwWAwUFxcTHFxMZubm8zMzNDZ2UlXVxcXLlyI+TqvvfYahYWF+x6nWvQViiOGLrynT5+OuEEWrZDr+Ww9fxsvbrebnp4eGhoaEkrL7Ie+qdjY2MiZM2dYXV2lvb2dsbExPB5PUq6RrNTKfuWH6enpnDx5koceeojMzEz+7M/+jG9/+9sJXTcSKiJXKI4YBoOBU6dO7ZqmHkq0kfDQ0BBpaWlUVlZGPGY/YdMrVOrq6nZ9O4h30zUa0tPTOXXqFIFAgPn5eXp6ekhJSSE1NTVhS9hkp1b2IjU1lby8PL74xS9y6dKlqK8hhOCXfumXEELwG7/xGzz11FMRj1VCrlAcMQwGAwUFBXsKZDQR+dTUFBsbGzQ2NkY8RhfiSKIUWqES7sFyGK6HRqMRq9WK1WplbW2NoaEh3G43BoOB0tLSmL8hHFZEHopeRx7LdX/+859TXl7O/Pw873vf+zh79iwPPvhg2GNVakWhOIbsJ+RLS0tMT09z8eLFPcVjv4h6eHg4WKESiYOKyMORnZ2N1WqloqICTdPo7Oykv78fp9MZ9RoHsdm5H/FsdpaXlwNQXFzMRz7yEa5evRrxWCXkCsUxZC8hdzqd9Pf3R+VRvpeQz87Osra2xpkzZ6I6P94GpVjR3QttNhttbW0UFRUxOjpKZ2cndrt93/u4WxF5LEK+sbERHCa9sbHBD3/4Q86fPx/xeJVaUSiOIZGE3OPxcP36dRoaGqKqW44k5OEqVPZC07TgP4PBkLSRapEIrZqJ1d/lbkXksdSRz83N8ZGPfATY8rP59V//dT7wgQ9EPF4JuUJxDAknwPqm5NmzZ6O2kQ23jl6h0tTUtG/+Wfch1zQNs9mMpmkEAoGgoOt2AofBXv4u+fn5QeG9GxF5rA1BtbW1XL9+PerjlZArFEeQ/XLXOyNyKSU3btygvLycgoKCqK+zcx39YVBfXx91Q04gEAgKttFoxGg0BsVcf89oNO6qBY8XKeWeIhrO32VkZITS0lKsVuuxyZHHghJyheIYslOABwcHg74lsRD6wNAfBhUVFeTl5e17rpQSs9nM+Ph4sAtRvzddaAOBAH6/P/i/RqPxUDdHdX8Xv9+P3W6nu7sbk8lEampqQoJ+1GxsE/7OI4SoFEK8JoToE0L0CiE+l4wbUygUkQkVoImJCVwuF6dOnYprHV1Yh4eHSU9Pj+phIKVE0zROnjyJ1WpleHiYjo4O5ubmtj1gjEYjKSkpWCwWTCZTMErXz4+HeB4EJpOJiooKWltbycvLY21tjY6ODmZmZuLydY/1IaBv0B4UyYjI/cDnpZSdQogsoEMI8a9Syr4krK1QKPZgYWGB2dnZuN0MdSGfmZlhfX2dpqamqM4L3djUR6Vtbm4yOTnJ6OgopaWllJeXBxt39CjdaDQyOjpKTk4OgUCAQCCwK+0S7X3HgxCC9PR0SkpKsFqtcfu7xBqRHzQJ34mUclZK2Xnn53XgFlCe6LpHlcuX7/YdKBRbBAIBhoaGaGpqijvaE0KwtrbG7du3960519GjaSHEtuPT09M5c+YMDzzwAGazma6uLnp6elhbWwseMzs7i9/vp6amBovFgsFgIBAI4PP5gnn1g0aP6HV/l9bWVnJzcxkcHKSrq4uFhYWo7uMwmqGiJak5ciFENdAEXAnz3lPAUwA2my2Zlz1UXnhBibni7uN2u3G5XLz97W9PqF1d0zSGh4dpbW2NqkNSSonf798l4qEYjUYqKiooLy9nZWWFsbExvF4vBQUFzM/PB0sahRBYLJZgukVPcewXpSfDNCt07dChEZubm0xPTzM2NkZRURFlZWUJ57YPY08gaUIuhMgE/j/g/5FSru18X0r5MvAyQGtr6+HtdigUx5C9hMrv99Pd3U1qamrEKe/R4Pf7WVlZ4dSpU1HN05RSBsU2GiEVQpCfn09+fj5ra2t0dnYGN0crKiqCAhmadtGvkUjaJZrPEen+I/m7lJeXk5ubm9AD5CAj+KT8hYQQZrZE/NtSyu8mY82jxOXLIMTWP3jrZxWZKw4bvbLEZrMl5EKoe6hkZmZGZW2rp1P2K/0Lh6ZpDA0Ncf78ed72treRlpbG9evXuXnzJqurq8GIVRdzi8WyLe3i9Xq3pV2SMSFoP3R/l5aWFmw2G3a7nfb2dqampvD7/Qld/yBIOCIXW4+ZbwC3pJR/nPgtHT0uX35LtIWAQ6yeUii20d/fT3Z2NmVlZYyPj8e9jl6holeQ7Ie+uRlPVDk0NER+fn7QV7usrAyr1YrD4WBycpLBwUEqKiq2TQcyGAy70i56JH3YMzuzs7PJzs7G5/MxOztLZ2cnHo8Hp9MZVePVYaRWkhGRvx34d8C7hRDdd/59MAnrKhSKEG7fvo3P5+PEiRPB1+IRCb1C5fTp01HZ0OrR8F558UjMzs7idruprq7e9ro+HejChQs0NDTgcrm4evVq0NlQx2AwYDabg9a1Qgg8Hs+2nHo8xPMgMJvNQX8Xk8kUtb9LvEOmYyHhiFxK+XPg6GzfHjDPP3+370BxPzI/P8/c3FzYAcyxiNLq6iq3b9/etuG4l5BHqlCJhrW1NSYmJmhpadnz3JSUlGBr/dzcHDdv3sRisVBZWUleXl7wXKPRiNPpxOl0cuLEiW3fEkKbkPYjnvRQKEIITCYTFy9ejMrfJVaflXhQnZ0xovLiisMgVPjW1taClSU7qy1iqWd2uVz09vbS3NwczK/v5aKobzzGI+Jer5fe3l4uXrwYdS7fYDBs8x2fnJxkaGiI8vJyrFYrgUCAvr4+GhsbSUlJCQq5Hp37/X5MJtO+f49ktOjrROPvchhCfnQq2o8wSrwVdwu3283NmzdpbGxMaACzXulSX1+/LWKMFJHrIh6P6GmaRk9PDydPnoy7qiY7O5tz587R1NSEz+fj6tWrXLlyhcrKyuD9GwwGTCYTKSkpmM3moMfLfjXpyRRyHd3fpbm5mdraWhYXF7l27RoTExOsr6/HLOQ/+MEPOHPmDCdPnuTFF1/c//rx3vj9xAsv3O07UNyP+P1+urq6IhpYRTtmTa9QsdlsuzxUwq0RWqESj+CNjIyQk5OT0IxQHYvFQk1NDcXFxWRlZWG32+nq6mJpaWnbfYdWuxiNRjRNw+/34/f7dwn6QQ9v1v1dmpubMRgMfOELX2B4eJgrV3a114QlEAjw27/923z/+9+nr6+PV155hb6+vRvl7ykhV5Gz4l5B0zSuX79OdXV1RAOraCPyoaEh0tPTgxNnQgkn5KG551gFz263s7GxQW1tbUzn7cXS0hKrq6tcvHiR1tZWTp48ydzcHFeuXGFiYmJbOaC+OaqLOhCM0kPLFxMR8mjTWbq/y/PPP09lZSXf//73o1r/6tWrnDx5ktraWiwWC5/4xCd49dVX9zznnhLyZEbOqnZccTeZnp4mNzcXq9Ua8ZhohHxmZgan08np06fDvr9TyBPZ3FxfX2d8fJxz584lLXXhdrsZHBzk/PnzQfHMysqivr6elpYWNE3j2rVr9Pf3s7GxETxPr0nXDbt2WgEkaqMby+fz+XxUVFRwOUrxmJ6e3jYsu6Kigunp6T3PUZudEVC144q7SWVlJT6fb89jIgq5z4f43vfwXL2K3+Wi4T/+x4jCs9PGdr/2+0j4fD56e3u5cOFC0krt9Fz72bNnw+aYzWYz1dXVVFVVsbi4yMDAALD1tyssLAx+htCadL/fj9PpJCcnB5/PF1fn6FGzsIV7ICJXkbPiXiQaMY0k5IZXX0X7wQ+YX1qi0mzG8hd/AcvLe64Ra/t9KFJKenp6qK2tTcgyYCfDw8MUFhbu640uhKCoqIjm5mbOnDnD0tISb775JuPj49sehgaDgYWFBSwWC4WFhbui9GiJNSJ3u90xCXl5eTmTk5PB36empsKmxUI5NkIeSZgvX96KlvWIWf85mUKuascVR5GIFSdvvMG0yUSRzYapuBi8XsTExJ5rxNt+D1ubm1lZWRQXF8d8biTm5ubY3NykqqoqpvMyMjI4e/YsbW1tGAwGOjo66OvrY319HafTyeTkJHV1dbusADRN22UFEIl4Bi/HIuRtbW0MDQ0Fzcb+7u/+jocffnjPc45NauVuug6q6F5xFAkXkUspWVhfJz8zk/SxMcTaGuj/wiCECNZhxyPi8/PzrK+v09jYGNdnCMfGxgZjY2P7NhLthclkwmazUVlZyfLyMkNDQ6yurlJbW7ttzZ1WAHuNp9M56MHLJpOJP//zP+f9738/gUCAJ598knPnzu19TtSrHwNU5Ky4nwgn5ENDQ6R+6EOUfuMbCI8HUlKQVitiaAjZ1gZhXA4dDgelpaUxe5o7nU5GR0cTEtydBAIBenp6qK+vT0quXQgRtM/NycnB6/Vy5cqViIMv9HvYOZ4u9G8Tq5B7vd6Y53V+8IMf5IMfjN7p5EgL+eXL2ytR9L/d88+Hj5JV5Ky4V4hGKHYK+czMDBsbG5x6//uR/f1bL6amIquqYH4eVle3CbmUkqKiIrxeL+3t7RQUFGxruNkLn89HT08P58+fT6qPSH9/P+Xl5WRnZydtzbm5ObxeL2fPng1+A5mdnaWrq4uMjAxsNtu26+nCHdo1GmoFcNCplXg48kJ+UJUjoWsrFMeR0Bz5ysrKWx4qBgOUlCCzsiA1FQKBrf/z6D/PzCA1jUBJCSaTierqamw2GwsLC/T09GCxWKiqqiInJyfsA0VKSW9vLzU1NVG5/0XL9PQ0Usp9N/ZiweVy7UrTRBp8oTsw6iId6pMeKuper/dAUyvxcKSF/CBRk34Uxx09OnS5XPT19W3zUNHe/W4MP/zhloBrGrK1FdLTMfz3/47o60MTAmGzEfit34KMjGCLeUlJCQ6Hg9u3b+PxeKisrKS4uHhbBDo6OkpGRgYlJSVJ+yzr6+tMTU3FPXs0HKHli+G+NYQOvnC73UxOTjI+Pk5xcXHYwRf68Oi1tTUsFgs+nw+DwRCM1CPh8Xii8nxPhGMj5Cr/rVBsx2Aw4PP56O7u5ty5c9tTIjYb2iOPgMOxlU4pLET86EeInh606uqtPO/EBIYf/hDtIx/Ztm5OTk7Q2W9ycpKxsTGsVivl5eWsrq7icDiiHtIcDaE16MmcND8yMkJhYSG5ubn7HpuamsqpU6eora1lbm6O69evk5qais1m2/bNZGNjg+XlZZqbm3f50UTaHI0nRx4rx778MNY1VM254l5BCMHk5CRVVVXhxSo7Gyor4c5AB2G3o6Wnv7VZl5mJYWYm4vq6uLW1tWE0Grl27Rq9vb27Kj8SQUpJX18f1dXVSa1BX1xcZH19fZcP+n4YjUbKyspoa2vDZrMxOTnJtWvXmJmZwev10tfXx7lz5zCZTNusAEJr0nduQCv3wyRzGDXnCsVhoTe3lJWVbb0gJSwuwtwchKmF1qqrwelEaBpoGmJ1FS1kSEUkTCYTVqsVg8FAbW0tIyMjdHV1sby8nPD0m4mJCVJSUigtLU1onVA8Hg9DQ0MJWQWEG3zxf//v/91lkxvNeLrD6Ow8NqkVheJ+Yj8Bmp6exuv1vuXFEghg+G//DcNPfwpCIOvrCfze78Ed10QpJf7WVky3b2P8+c8hEEBrbUV797v3vRd9c7Oqqgqr1YrNZmN9fZ2JiQmGhoaorKyktLQ05jr01dVV5ufnaWlpiem8/e61p6eH06dPJ008U1JSyMnJweFwYLVagxvCOwdfhBtPB1umX4nMV42G+1bIVc5dcVxZWVlhYmICm82G1+sFQPzkJxhee22r1FAIRG8vhr//e7SPfxzx3e8iurowlZWhffSjGIaGEIODGK9eRebkQGYmZGWhtbZCGPEbHx8nNTV1m4FXVlYW586dw+PxMDU1xZUrVyguLqaysnKXb3o4vF4vt27dorGxMaFpPeHuNScnh4KCgqSt6fV6GRoaorm5mZSUlIiDL0LnjeqRutfr5Sc/+Qm//Mu/nLT7CYc4jMGgO2ltbZXt7e2Hfl3F/YEQokNK2XqXLp+0/0N5vd5dqQuXy0VnZyfNzc2sra2xvr7OyZMnMbz8MuLf/g30FIXDARsbMDyMYXh4azPIZIKcHLRLl5DV1YjbtzG8+SZaQwOkpqKdPo3/P/0nCBHixcVFbt++TVNT056CGwgEsNvtTE1NkZmZSVVVVcTSRCklXV1d2Gy24EDmZLC6usrQ0BAtLS1JezhIKbl+/Trl5eVh/dW9Xi/T09PY7Xby8/OprKzc5h3/5S9/GYvFwvPJixzDflW7byNyheK4oU/50StUnE7nWx7bVVUYXK6g8IveXsTsLKysgNcb3NkX6+tbYl9ejuGNN7Zy5qOjaA88sBWp9/Yi71SkbG5uRi2MRqOR8vJyysrKgi3xUkpsNhsFBQXbUkWjo6Pk5OQkVcR9Pt+BRPjT09OkpKREHJKhD76oqqoKjnozGAwUFRUxMzPD66+/zuuvv560+4nEfbXZqVAcV6SU3LhxY1uFSmhnp3z3u9F+4RcQU1NbBlnr62hlZVubmoDQNMSd3X1DTw+iuxs2N8FshtRUDDdvIj0exJ1UTSAQ4ObNm9TX10eVKtHRW+Kbmpo4ffo08/PzXLlyhampKQKBAIuLizgcjqQOntBz+CdOnIiqKzVaNjY2mJqaiujlHopeh9/S0sLJkyf57ne/yyOPPMJDDz20rx1xMlAReZSoTlDF3WRwcJCsrKy3KlTY4X5oNqN9/vNoMzPgdGJ8/nnE4OCuMi1pMoHBgOHmza0uT00Do3FL1DUN7cSJoDBWVlaSk5MT9z1nZmZSX18fTD+8+eab+Hy+pHqzAExOTpKamppU90VN0+jr66O+vj7m2vbMzExmZmb43Oc+R3p6+rYJRgeFisijRM3tVNwtpqen2dzc5OTJk9teNxgMaH4/6BGfEFBejjx9mkBd3Za/yp2OQsmdWZOpqcEcuCwshIwMxNgYmM0EPvtZKCzk9u3b28saE0Rv+TebzVRWVnLr1i1u3rzJWgRHxlhYX1/Hbrdz6tSpJNzpW4yOjlJYWBiX58vPfvYzent7ee6553jmmWeS6hsTCRWRKxRHFCEEy8vLTE5Ohm1dT/npT7F985uYMjO3SgmfegqZlrZlw/q2t2F6/XWk243MysIwM4M0GiE7e+s1qxVhsSDT0xEOB4F3vQvtwQdZWlpicXGR5ubmpH6WwcFBSkpKqKqqora2ltXVVcbGxvD5fNhsNoqKimKO0v1+/4F0hK6urrK6uhpXWeTa2hpf+MIXePXVV5N6T/uhIvI9UJ2girvJ5uYmfX19NDY27qpDFrdukfa3f4s/OxtZWYnh6lXE//yfwQERorgYabOhvfe9aO9731Zbvs1G4O1vR9bWgsGAtFrB79863uHAPT/P4OAgFy5cSOqGod1ux+PxYLPZtu5dCPLy8mhoaKC+vp6VlRXefPPNXYOU96O/vx+bzZbUjlC/309/f39czURSSr74xS/yO7/zOzF3lCaKEvI9UJ2girtJb28v586dC+vTIcbGwGBAM5u3GoBKShA3bgTtVuX58wTe/naYnEQsLUFGBv5f/3W0Bx9ENjQgNjcRw8NbJYnFxWipqcx9/evU1dUltQtxY2OD8fFx6uvrwwpjeno6Z86cobW1FSkl165dY2BgAJfLtee6M3esBZKV/tHp7++nqqoqrk3TH/zgBywtLfHpT386qfcUDUc6taI2GBX3M83NzRGjQpmbiwidBr+2hlZV9dasTyEIPPnkVuemy4UYGsLQ2QluN7K0FP9nPoPxe9+DkhIC9fXY19YocjpJj8JgKlr0IRHnzp3b17PcbDZTVVVFZWXlvna6GxsbTExM0Nqa3FYBu92OlHJb41O0LC4ucvnyZf71X/81qd9mouVIC/lRspp917vu9h0o7jd0H+xwyAceQDY3Y/npT8HnQ0tLI/CpT20XfiGQNTVbx9fXo126hHC5kEVFCLt9K29us7GytoZ5bY2spiaiH0G8N1JKbt26RUVFRUwWrjvtdCcmJhgcHKSyspKSkpJgRU19fX1S297dbjdjY2NxPRyklDzzzDP8/u//flI9Y2LhSAv5UeLf/u1u34FCEYLJhP+ZZ5g4cYLMkycJ2GwYQqbNi5ERxOoqMj8/KOaUlgbbTuWJEwR+8Rfx/su/4F1bo6SujsCv/mrSbm96ehohREJDInJycrhw4cI2r3Cj0UhRUVFSK0H0h8OZM2fimnb093//96SmpvLxj388afcUK0dOyGMd76ZQ3Kvst9kmTCbWKysJ7NicNPzLv2D80Y+26sM1jcCHPoS28yulEGw+9BA3heDimTMEysq2teYnwtraGjMzM0kzw9LtdDMzMxkfH2dubg632520jc7bt2+TnZ1Nfn5+zOfOzMzwX//rf+X1119Pam18rCQlmSOE+IAQYkAIMSyEeC6RtY7SBqOqWlEcdVJTU+ns7GRhYWErX768jPEnP0HabFv/Kiowfv/74HRuOy8QCHDjxg1OPvAAKdXVSRNxn89HX18f58+fT2r5ncvlYnx8nNbWVi5dukRhYSH9/f0J2+mur68zPz/PiSjsfHeiaRq//du/zVe+8pW45MiTiAAAGYBJREFUHgLJJGEhF0IYga8BvwzUA78mhKhPdN3DYi9RPkoPFYUiFH0yzYULF6irq2NxcZErV65gv30baTBsReMAJhNSCPB4gudKKenv78dqtZIXko5JFD1FUVtbu804KlE0TaO3t5e6ujrMZjNCCIqKioLt8LOzs1y9epXp6emIewrhCAQCwUER8WxQ/vVf/zW1tbW8//3vj/ncZJOMiPwBYFhKOSql9AJ/B3w4CeseitWs6thUHDeklMEyQyEEGRkZ1NXV0dzczEZqKhNOJ47BQQJu95ZxVmkphFSjTE1NAVBRUZHU+7p9+zbp6elJbZWHrS7L/Pz8sFOQdDvdxsZG3G43V65cYWRkJGjvuxdDQ0OUlZXFlZ4ZHR3lm9/8Jl/96lfvakpFJxlCXg5Mhvw+dee1bQghnhJCtAsh2hcWFqJa+LAi32iuo/zLFUcBfU4k7M6hWywWTtTVUfKlL+EvK8N+6xazGRlsfOITwQh9dXWV2dlZzp49m1QBWl5eZnFxcZeNQKIsLS3hcDio0TdsI5CSksKJEyd44IEHSE1Npauri97eXtbX18Mev7i4iMvliuth5vf7+exnP8uf/dmfJbUZKRES9iMXQnwc+ICU8t/f+f3fAZeklE9HOudu+5Hv3FDVede74BAcJxUHzL3iRx4IBLZ1OuqRuKZpUaUCNE1jfn6eiYkJ0tPTsVqtDA4O0tTUlNRhwB6Ph87OzgNbVx/oEAtSSpaXl5mYmEBKSWVlJYWFhQgh8Hq9dHR00NLSEpOzo85LL72Ew+HgD/7gD2I+NwkcmB/5NFAZ8nvFndeOLKGNRnpQIuVbPysUR41QEY82kjYYDJSWllJSUsLS0hI3btwgPT2dzc1NUlJSkhKRa5oWHK2WTBHX8+2nTp2Kq9NUt9MtKCgINhCNjIxQXl4e/OYQj4j39vbyj//4j/zsZz+L+dyDJBmplWvAKSFEjRDCAnwC+KckrHug6BUpOvrPaiNTcRQJFfFYBVgIwcLCArW1tdTX1zMzM8O1a9ew2+0xbQ6GY2RkhLy8vKSOVoOtfHtWVlZShk+E7iHoqZrV1VXcbndM63i9Xp5++mm+/vWvH/gw5VhJWMillH7gaeBfgFvA30spexNdNxlEU5Gys7z2hRdUiaHiaLFzczNWpqam8Pv92Gw2srKyOH/+PBcvXsThcHD16lUmJyeDefdYWFhYwOl07pu/jpXV1VUWFhbiKgncC5/Ph9vt5u1vfzuZmZncuHGDmzdv4nA4ojr/xRdf5OGHH6bpzgSlWJicnOQXf/EXqa+v59y5c/zJn/wJsLW38L73vY9Tp07xvve9j5WVlZjXhnt8ZqcQb5UORnNctMcrjjb3Uo7c5/Ph9/vjFnGHw8HAwAAtLS1h67p9Ph9TU1PY7faYhie7XC66u7vjzjNHwufz0dHRQUNDQ1Kn/WiaRkdHB2fOnAl2hUopWV1dZWJiAp/PR2VlJcXFxWH/zteuXeM//+f/zGuvvRaXNcDs7Cyzs7M0Nzezvr5OS0sL//iP/8j/+B//g/z8fJ577jlefPFFVlZW9su9q5mdkVAVKYqjSqQKlWjweDxBG9xIzTlmszk4c3J2dpbOzk5ycnKoqqqKWAuuaVpcY+D2Q0pJX18fNTU1SRVxCD8oQrfTzcvLY3Nzk8nJSUZHRykrK6O8vDwo2JubmzzzzDO88sorcfu7WK3WoBlXVlYWdXV1TE9P8+qrrwZnej7++OM89NBDcW2i3nM2tvF0Y+rvKUFXHCX++I//mO985ztx5bF1sT19+nRUomgwGCgvLw92Tfb19XH9+vWwaYeBgQFKS0sTGgMXjunpacxmMyUlJUldd2VlBYfDsadHuG6n29bWBhC0011YWOD555/n05/+NGfOnEnK/YyPj9PV1cWlS5eYm5sLCnxpaSlzc3NxrXlPCnloeiSWbkyVF1ccJR599FG6u7t58MEH+cu//Mt9PbpDGRwcpLCwMOZNSL1rsrW1laqqKsbHx2lvbw9aAMzOzgbTEMlkfX2d6enppImljs/nY2BgIKIf+k5MJhNVVVVcunSJ3NxcHnvsMV599dW48uLhcDqdfOxjH+Oll17aZfwVb/oM7kEhVyjuFWw2Gy+99BI/+tGPWFlZ4V3vehdf/epXWV1d3fO8mZkZvF4vVVVVCV0/NzeXhoYG6urqWFhY4I033mB4eDjpzUShrfLJHo82MDAQ16AIg8FAamoq6+vr/MVf/AU3btxI+F58Ph8f+9jHeOyxx/joRz8KQElJCbOzs8BWHj3erth7Ssh3plVAVaAojj+FhYVcvnyZK1eukJOTwwc+8AG+9KUvYbfbdx27trbG5ORkXKPKIpGRkcHp06cBKCgooKOjg/HxcXz60OcE6e/vp6KigszMzKSsp5PIoAgpJb/3e7/Hs88+y6/8yq/w2c9+NqF7kVLymc98hrq6Op599tng6w8//DDf+ta3APjWt77Fhz8cn7vJPVu1oipQ7l/ulaqVSPh8Pl555RX+9E//lIaGhv+/vXsPqrpcFzj+fdCtm5SOAnuJcpEMRIzIErULGamkpuKp3VaZKTW2qaWl7mYn2ZRdZhg3hzx6sjxxxj1hU5JTR3PYHo9YpuUcZaND5GUrmCg3FyqopBHBes8fLNZGXUthXVi39zPDuK78Hn/os17e3/s+Dy+99BKxsbGWHYv33HOPU4tWKaU4fPgwISEhDBo0iNbWVqqrq6muriYkJISoqCi7NwPV1tZy/vx5EhISnDrKb2pqsqyqsafGeEFBAfn5+Xz++edO6fjz3Xff8fDDD1/TDzUrK4sxY8YwY8YMzpw5w+DBg9m8efOtKilaPUk6kWs+x9cTeTuTyURBQQE5OTmEhIRw9epV1q5da2ly7CyVlZU0NjYyfPi1RU2vLwEQHR3dpVH11atXKS0tJSkpyandfpRSHDp0iCFDhthV3bGuro5p06axa9cup194dQKridynplY60itQNF8XEBBAWloae/bsoXfv3hiNRhYvXszu3bsd3rHZ7tKlS9TW1lq9CNleAmDUqFGWOi6drQ/evrU/Pj7eqUkc/tkowp4kbjKZWLp0KW+//bYnJnGbfDaR63lxzZ+kpqZy6NAhVq9eTX5+PhMmTGDLli127dhs19zczNGjR7n77rtvehGyva7JfffdR0xMjKUEgNFotJnQy8rKGDBggNOXMDrSKAIgPz+f22+/nX91Ytu77uCzUyua//KXqZWbOXXqFO+++y779u1j3rx5pKend2keWylFSUkJERER/O53v+vy8X/++WfOnDlDQ0MD4eHhDBo0yPJhcO7cOaqqqhgxYoTTV78UFxeTkJBgV3nZqqoq/vCHP7Bnzx6rtc89hO9PrXSsatgdx7J2W9M8wR133MG6devYuXMnZ8+e5ZFHHmHNmjVcvny5U+8/deoUffv2tSuJAwQGBhIXF8fIkSNpaWmhqKiIkydP0tjYSHl5uVNX1bRzpFFEe9u2nJwcT07iNvnUiLxjSVpX63gxVV9Y9Sx6RH6jxsZGcnNz+eijj5g0aRIvvPCCzTngCxcuUFFRwb333uuUFRvQlihramo4ceIEwcHBDB061Kkra86fP09lZaXdo/zc3FzKy8t57733PKLjz014x4hcj241zfmCgoJ4+eWXOXjwIPHx8Tz11FMsXbqUH3/88ZrXNTU1ceLECRISEpyWxKHtwugvv/xCVFQUgwYN4ujRo5SWlna68uDNNDc3U1ZWZvcov6ysjI0bN5Kdne3pSdwmj0vkXe2haWsTkCs2Atmq4+LKY2qaM/Xq1YuMjAyKioqYPHkyCxYsYO7cuZSWltLU1GSpEOjsetv19fVcvHiRO++8E4PBQFJSElFRUVRUVHDw4EFLCYCuai+0ZW+jiJaWFhYtWsQHH3zg1N8QupvHTa04Mk1xs6kVZ8+f66kVz6WnVjrPZDKxd+9eVq1aRV1dHcnJyWRlZTl1NN6+UclWK7grV65w+vRpGhsbiYyMJCwsrNPHr6qqorGxkfj4eLtiy8nJoampiaysLLve7waeO7ViT8XCrurqSF/T/EFAQAApKSksXryY22+/nbq6OiZOnMi2bdscWrrYrr1lW0xMjM1VM3369GH48OGMGDGCq1evcuDAASoqKq7pV2rNlStXqKqqspQP6KrS0lK2b9/OSh/YdOIR9civ76Fp7+i2O38eHY/lA/8OND+XnJzMQw89RP/+/Tl58iQ5OTlkZ2czf/58Zs6cafdUS/uuz86sfunduzcxMTFER0db1qKHhoYSGRl5w4eAyWTiyJEjDB8+3K5CW7/88gsvvvgiGzZs8Li2bfbwqamV6735pvWR+MqVei7bl+mpFec4e/Ysa9eupaCggKeffpq5c+cSFBTU6fe3dydKSkqya6rGZDJhNBqprKykT58+DB482FICoLy8nB49etjdZm7lypWEhoayfPlyu97vRt5Ra8VVa8H1PLb/0IncuS5dusSHH37Ixx9/zNSpU1m4cOEtR9gtLS0UFxeTmJjo8EVEpRT19fWcPn3asou0rq6OkSNH2rXKZP/+/bz11lt89dVXTi8P0A28I5G7ik7k/kMnctdoampi48aNrF+/nvvvv58lS5ZYLdCllOKHH37AYDAQFhbm1BgaGhooKSmxFOmy1WPTlp9++omJEyeyefNmYmNjnRpbN/Hci53dQc9ja5pjfvvb3zJ//nyKi4sZN24cGRkZzJs3jyNHjlyzdLCmpoYePXo4PYlDWzu4YcOGkZiYyMWLFzlw4ACVlZWdujCrlOL1119n/vz5DiXxjIwMDAYDCQkJlsfq6+tJTU0lNjaW1NRUGhoa7P7+9vCbEbnmP/SIvHuYTCa+/vprsrOz6dmzJ3/6058IDg7m3LlzPPjgg07v9nP27FnOnTvH3XffbXns119/paqqirNnz2IwGIiMjLS5nvyrr77igw8+4G9/+5tDyyv37t1L3759mT17NocPHwbglVdeITg4mMzMTFatWkVDQ4NdTZQ7wb+nVjT/oRN591JKcfDgQbKysti/fz+vvPIK8+bNc+pa9Fs1ijCZTNTW1lJZWUm/fv2Iioq6Zm6+oaGByZMns337diIiIhyOp6KigqlTp1oSeVxcHN988w0DBw6ktraWlJQUjh8/7vBxrPDvqRVN01xDREhKSmLYsGEsXLiQo0ePkpKSwqeffuqUdnDta9Hj4uJsdvsJCAggPDycMWPGEBwcfE0JAKUUf/7zn1m+fLlTkrg1RqPR0lIuLCwMo9HokuPY4nWXbDVN80yvvvoqffv2RUSoqalhzZo1PPzww8yePZs5c+bYVZUQutYoQkQwGAwYDAYuXrxISUmJJa709HS7jt9VItLtNVv0iFzTNKcICgqyJLBBgwaRnZ3Nnj17aG5uZty4cWRlZXHhwoUufc/Lly/b3SiiX79+DBs2DKUUQ4YMYe3atV3+Hp01YMAAamtrgbY+pAaDwWXHskYnck3TXKZ///6sWLGCoqIiwsPDmTZtGsuXL6eqquqW721tbeXYsWPcdddddm8oWrJkCe+88w4bN25k2bJl9vwVOiUtLY28vDwA8vLymD59usuOZY1O5JrmRXbs2EFcXBwxMTGsWrXK3eF0WmBgIM8//zzFxcU8+OCDPPPMMyxYsIBjx47dtB1ceHi43VMyn3zyCaGhoUybNs2R0G+Qnp7OAw88wPHjx4mIiGDDhg1kZmZSWFhIbGwsu3btIjMz06nHvBW9akXzOb66aqW1tZWhQ4dSWFhIREQEo0aNYtOmTTd0t/cGJpOJwsJCsrOz6dOnD8uWLWP06NGWqRlH28GdOXOGmTNnsnfvXqf3BXUz569aEZF/E5F/iEipiGwREe/rkaRpXqKoqIiYmBiGDBlCr169mDVrFl9++aW7w7JLQEAAEydOZNeuXbz22mu89957TJkyhZ07d1JTU0NJSYndjSJaW1t54YUXWLNmja8lcZscnVopBBKUUonACeBVx0PSNM2a6upqIiMjLfcjIiKorq52Y0SOExHGjBnDF198wfr169m6dSsTJkxg3759dq/8yM3NJTExkZSUFOcG68EcWn6olNrZ4e5+4CnHwtE0zR+JCPHx8aSmpqKUoqWlhbFjx/Lss8/y9NNPd7rw1vHjx9m0aRPffvut17Zts4cz15FnAJ/ZelJE5gPzAauFdjRNu7nw8HAqKyst96uqqggPD3djRM735JNP8uSTTxIYGMiFCxd4//33efTRR/n973/Pc889d9O15L/++iuLFy9m/fr1BAYGdmPU7nfLqRUR2SUih618Te/wmteAFuATW99HKZWrlEpSSiV1psi8pmnXGjVqFGVlZZw6dYrm5mby8/NJS0tzd1hOFRgYaEnCISEhvPHGGxw4cICQkBAef/xxVqxYQU1NjdX3rl69mvHjxzNq1KjuDNkj3HJErpSacLPnRWQuMBUYr9yxBEbT/ETPnj1Zt24dEydOpLW1lYyMDO666y53h+Vyt912Gy+++CILFy7ks88+Iz09nYSEBJYsWUJsbCwiQklJCYWFhXzzzTfuDtctHFp+KCKTgNXAI0qpc519n15+qLmSry4/dJbo6GiCgoLo0aMHPXv2xNv+L5pMJrZv305OTg79+/dn0aJFZGZmkpeX5w8fbFYn/h2dI18H9AYKzRcW9iulFjr4PTVNc7Hdu3cTGhrq7jDsEhAQwNSpU5kyZQr79u1j6dKlJCYm+kMSt8nRVSsxzgrEFVzVNk7TNPcTEZKTkykuLra5O9Rf+PQWfWuNlzXN34kIjz32GCNHjiQ3N9fd4TiFPy01tEaXsdU0P/Pdd98RHh5OXV0dqampDBs2jLFjx7o7LM0BPjcif/PNtkbL7R/Q7bf1FIumtWlfe24wGHjiiScoKipyc0Sao3wykSvV9gX/vK0TuabBlStXaGxstNzeuXPnNU2ENe/kc4lc0zTbjEYjycnJ3HPPPYwePZopU6YwadIkd4flVt5aGrgjn54jX7nS3RFommcZMmQI33//fZffl5GRQUFBAQaDwdJwuL6+npkzZ1JRUUF0dDSbN2/uVDs2T9La2sqiRYuuKQ2clpbmdaWBfXpErqdTNM055s6dy44dO655bNWqVYwfP56ysjLGjx/vlaNZXykN7NOJXNM05xg7dizBwcHXPPbll18yZ84cAObMmcPWrVvdEZpDfKU0sE7kmqbZxWg0MnDgQADCwsIwGo1ujsh/6USuaZrDRMQrN+X4Smlgncg1TbPLgAEDqK2tBaC2thaDweDmiLrOV0oD60SuaZpd0tLSyMvLAyAvL4/p06ff4h2ep2Np4Pj4eGbMmOGVxbccKmNr90FFzgGnu/3AbUKB8246ti2eGBN4ZlydiWmwUkp3L3EiEdkEpNB2/o3ASmArsBmIou3/8wylVL27YvRnbknk7iQixW6sVW2VJ8YEnhmXJ8akdY6I/JW2JjR1SqkE82NvAs8B7f0MViiltrsnQu+lp1Y0TesuHwHWtpH+u1JqhPlLJ3E76ESuaVq3UErtBfTUiwv4YyL3xALMnhgTeGZcnhiT5pjFIlIqIn8VEe/a4+8h/G6OXNM09xGRaKCgwxz5ANouXivgHWCgUirDbQF6KX8ckWua5iGUUkalVKtSygT8FzDa3TF5I53INU1zGxEZ2OHuE8Bhd8XizXw2kYvIJBE5LiLlIpJp5fneIvKZ+fkD5l/5XBlPpIjsFpGjInJERJZYeU2KiFwSkRLz1xuujMl8zAoR+cF8vGIrz4uI/If5PJWKyH3dEFNch3NQIiKXRWTpda/p9nOlOca8Fv3/gDgRqRKRPwLZ5n9/pcCjwDK3BumlfLIeuYj0AN4HUoEq4O8isk0pdbTDy/4INCilYkRkFvAXYKYLw2oBXlZKHRKRIOCgiBReFxPAt0qpqS6Mw5pHlVK2NtlMBmLNX2OA9eY/XUYpdRwYAZafZTWwxcpL3XGuNDsppdKtPLyh2wPxQb46Ih8NlCulflRKNQP5wPX7h6cDeebbnwPjxYVVf5RStUqpQ+bbjcAxwBuq80wHNqo2+4F+1/067GrjgZNKKXftBNY0j+eriTwcqOxwv4obk6blNUqpFuASENIdwZmnce4FDlh5+gER+V5E/kdEuqPogwJ2ishBEZlv5fnOnEtXmgVssvFcd58rTfNIPjm14slEpC/wBbBUKXX5uqcP0VYn5CcReZy2WhaxLg4pWSlVLSIGoFBE/mHeuOF2ItILSANetfK0O86VpnkkXx2RVwORHe5HmB+z+hoR6Qn8C3DBlUGJyG9oS+KfKKX++/rnlVKXlVI/mW9vB34jIqGujEkpVW3+s462eejrl3915ly6ymTgkFLqho4F7jhXmuapfDWR/x2IFZE7zKO6WcC2616zDZhjvv0U8LVy4e4o8/z7BuCYUmq1jdeEtc/Ti8ho2n4+LvtwEZE+5guviEgf4DFuXP61DZhtXr1yP3BJKVXrqpiuk46NaZXuPlea5sl8cmpFKdUiIouB/wV6AH9VSh0RkbeBYqXUNtqS6sciUk5b/YdZLg7rIeAZ4AcRKTE/toK2EqAopf6Ttg+U50WkBfgZmOXKDxdgALDFnA97Ap8qpXaIyMIOMW0HHgfKgavAsy6Mx8L8wZIKLOjwWMe4uvtcaZrH0lv0NU3TvJyvTq1omqb5DZ3INU3TvJxO5JqmaV5OJ3JN0zQvpxO5pmmal9OJXNM0zcvpRK5pmubl/h9cNht6HE7AZwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_samples = 20  # nb samples\n",
    "n_noise = 10  # nb of samples (noise)\n",
    "\n",
    "p = ot.unif(n_samples + n_noise)\n",
    "q = ot.unif(n_samples + n_noise)\n",
    "\n",
    "mu_s = np.array([0, 0])\n",
    "cov_s = np.array([[1, 0], [0, 1]])\n",
    "\n",
    "mu_t = np.array([0, 0, 0])\n",
    "cov_t = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])\n",
    "\n",
    "\n",
    "xs = ot.datasets.make_2D_samples_gauss(n_samples, mu_s, cov_s)\n",
    "xs = np.concatenate((xs, ((np.random.rand(n_noise, 2) + 1) * 4)), axis=0)\n",
    "P = sp.linalg.sqrtm(cov_t)\n",
    "xt = np.random.randn(n_samples, 3).dot(P) + mu_t\n",
    "xt = np.concatenate((xt, ((np.random.rand(n_noise, 3) + 1) * 10)), axis=0)\n",
    "\n",
    "fig = pl.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax1.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "ax2 = fig.add_subplot(122, projection='3d')\n",
    "ax2.scatter(xt[:, 0], xt[:, 1], xt[:, 2], color='r')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute partial Gromov-Wasserstein plans and distance\n",
    "-----------------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----m = 1\n",
      "Wasserstein distance (m = 1): 62.612867557378095\n",
      "Entropic Wasserstein distance (m = 1): 64.09799387131392\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEtCAYAAADHtl7HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZRdVZnn8d+PANJJgBADIfJiVFAHcRraNDI2KI5oo42Drm5RVMTVSnBaZmTUEZuZaUK3jtgjvg2Oq0EYEBRlfEF0aNR2QMQXJLhQXhUai05CXoEAAVGSPPPH2QU3ZdXZe9e9t+6tyvezVq3cOvvUPs89N/ep5+576rmOCAEAAKDcDoMOAAAAYLqhgAIAAKhEAQUAAFCJAgoAAKASBRQAAEAlCigAAIBKFFCYUWzfZvuoQccBYOrY3t/2JtuzBh3LdGH7Lba/M+g4pjMKqBnI9l/b/scx2+6aYNubpja6OrbD9gGl+0fECyLi2j6GBGyXbI/Y/k0qVEa/zi382Wttv7NfsUXEv0TE3IjYUvoztndM9+HFHdveknLO2G139jrmXrO9zPalpftHxBci4lX9jGmmo4Cama6T9JLRV2O2F0naSdKhY7YdkPYdGNs7DvL4AKq8NhUqo1+n9mLSQeSBiNgs6ceSXtqx+aWS7hxn20DzpESuHEYUUDPTjWoKpkPS90dKukbSL8ds++eIuM/2p2yvsP2w7ZtsHzk6ke3DbC9PY2ttfzxt38X2pbbvt73R9o22F6ax3W1fYHu17VW2P9RRuL3d9g9tf8L2/ZKW2T7A9vdtP2R7g+0vp31Hk9bP0yvFN6btx9q+OR33R7b/dUe8I7aPTreX2b7c9udtP5Le3lvSjxMObM/S8/p62x+z/aDtX9t+dRr7sJp8c27nqlVa6Xm37bsk3ZW2vSTlkofSvy/pOMa1tj9i+6cpH33D9vw0tjjNt2P6fr7t/237vhTPFROEfp22LZaOlPTRcbZdl+Y9zPaPU+5Zbftc2zunMae8ti7Fd4vtg9PYa2zfnvLQKtvv77hfuXx2uu1fSHo0rZqdnuZ4xPYvbb/C9jGSzpD0xnSOf55+PpeLr+84Vth+l5t3Jjba/oxtF/0H2F5FBF8z8EtNwfSf0u1zJf2lpA+P2XZhuv1WSU+XtKOk90laI2mXNPZjSSem23MlHZ5unyLpm5JmS5ol6UWSdktjX5f0D5LmSNpL0k8lnZLG3i5ps6T/kI73B5Iuk/Rf1BT0u0g6ouN+hKQDOr4/VNI6SS9Oxz1J0oikp6XxEUlHp9vLJD0u6TVp349I+smgHxu++JqOX53PrXHG3i7pCUknp+fav5d0nySn8WslvXPMz4Sk70qan/LAfEkPSjox5YYT0vdP75hjlaSDU275qqRL09jiNN+O6fv/K+nLkvZQ82LyZRPE/TJJD6Tcs0DSvSmnre3YFpL2T/u/SNLhKb7Fku6QdFoa+1NJN0maJ8mS/pWkRWlstaQj0+09JP1Rul2Sz26WtF86R8+TtELSMzru93PS7WWj56Pj/uVy8fVjHo9vpfj3l7Re0jGD/n83zF+sQM1c39dTr6KOlPSD9NW57fuSFBGXRsT9EbE5Is6R9DQ1T1SpSYoH2F4QEZsi4icd25+uprjZEhE3RcTDaRXqNWqSyqMRsU7SJyR1Xmt1X0T8z3S836S5nqkmKTweEddrYksl/UNE3JCOe7Gk36pJauO5PiKuiubaiEsk/WHmvAGY2BVpdWL06+SOsXsj4vz0XLtY0iJJCzPzfSQiHkh54M8k3RURl6TccJmat9Ne27H/JRFxa0Q8Kum/STreYy4cd3N5wqslvSsiHoyIJyLi+xMc/wY1BdML1eTE6yPiMUm/7tg2EhH/Ikkpz/0kxTeipjh5WZrrCUm7Snq+msLxjohY3TF2kO3dUkw/S9tL8tmnI2JFOkdb1OTng2zvFBEjEfHP492xwlw81tkRsTHd32v01DsWGAcF1Mx1naQj0hL3nhFxl6Qfqbk2ar6aV3Gjy9Lvt31HWjbfKGl3Na+8JOkdkp4r6c60pH5s2n6JpG9L+lJaJv972zupKYR2krR6NMmqSTJ7dcS2YkysH1Dziu2n6W22v2y5X8+U9L7OJK7m1dkzJth/TcftxyTtYq4lACbrdRExr+Pr/I6xJ59rqQiRmlXrNp254BlqVoA63Stpnwn2v1dNrlmgbe0n6YGIeDBzbEXE42pWZV6avn6Qhq7v2Pbk9U+2n2v7W7bX2H5Y0n8fPX5E/D81K/ufkbTO9nm2d0s/+udqipl73Vyu8G/S9pJ89uR9joi7JZ2mZrVpne0v2Z4o95Xk4rHG5svc47ddo4CauX6sphA6WdIPJSkiHlazrH6ymlWgX7u53ukDko6XtEdEzJP0kJqCRhFxV0ScoOZJ91FJX7E9J72qOysiDpL0EknHSnqbmif7byUt6Eiyu0XECzpii85AI2JNRJwcEc9Q89bg//LEf3m3QtKHxyTx2enVKoDhFAXb71PzS7/T/mrethu135ixJyRtGPMzKyTNtz2vMLbR66BGV+qlp1brn7z+KfmsmlWxAyNiNzXXHT15nVBEfDoiXiTpIDUvPP9z2n5jRBynJo9eIenyjlhz+WxsvvxiRByh5lyFmrz8e/upLBejCxRQM1Ra7l0u6b16KilIzSur9+qppLCrmmuS1kva0fbfSBp91STbb7W9Z0RslbQxbd5q++W2X5iWzx9Wk8i2piXr70g6x/Zutnew/Rzbo8vcv8f2G2zvm759UE0i2Jq+Xyvp2R27ny/pXbZfnC7anGP7z2zvWnWCAEylsc/j8Vwl6bm235wuln6jmkLkWx37vNX2QbZnS/pbSV+JMa0LUg76RzUvxPawvZPtzovCx7pO0svVFGe3p20/lHSUmrewOguoXdXku022n6/mWi9Jku0/TnlpJ0mPqrn+cqvtnd20Qtg9Ip5IPz+a36ryme3n2f63tp+W5v+Nts2Vi23v0HEeqnIx6lBAzWzfV/OKp/Oaoh+kbaNJ4duSrpb0KzVL4o9r22XyYyTdZnuTpE9JelMqzvaW9BU1yeCOdKxL0s+8TdLOapLRg2m/RS1x/rGkG9IxrpT0noi4J40tk3RxWoI+PiKWq1lBOzfNfbeaiyEB9N83vW0fqK8X/tynJP2Fm7+I+/R4O0TE/WpWst8n6X41K+PHRkTnCtMlki5S+kMXSf9xguOdqOZF3Z1qLtI+rSW2H6lZrb8hIl1N3RxzvaR16fKHUe+X9GZJj6gpfr7cMbZb2vagmlx6v6T/0RHPSHrb712S3pKOU5vPnibpbDWrbmvU5PK/TmP/J/17v+3Ra6xqczEqjP6FBAAAQ8v2tWr+yuxzg44FkFiBAgAAqEYBBQAAUIm38AAAACqxAgUAAFCJAgoAAKBSVx2Z3XyA4afUfIbP5yLi7Pb9Z0fzMTuTt0irW8dX8xeawJBZvSEi9hx0FOOpyWGesyA0b/HEk5VcDZHLuJsL5si97N2aGS8xK7+LtmTGS16e9yLWXsg9dr34SN2peNxW39eDSbCtifPXpAuo1EDxM5JeKWmlpBttXxkRt0/8U/PUfPTP5C3VWa3jZ3U5P4BeO2vsx3MMheocNm+x9FfLJ57w8YKDjv3QkbE2ZsalpvtRm00Fc+SUvM7NHScXZ8kcJb+hckVnyRy5xy53X0oK39yHovTicfvQsh5Mgm1NnL+6eQvvMEl3R8Q9EfE7SV+SdFwX8wHAVCKHAZi0bgqofbRtx+qV2vZDHwFgmJHDAExa3y8it73U9nLby5sPdwaA6WGb/PXo+kGHA2CIdFNArdK2n4y9r7b91GxJUkScFxFLImKJNLuLwwFAT2Vz2Db5a85QXgcPYEC6KaBulHSg7WfZ3lnSm9R8ECwATAfkMACTNum/wouIzbZPlfRtNX/0emFE3NazyACgj8hhALrRVR+oiLhK0lU9iqXIWTqz6znOzLZC6P4YAIZfVQ7borI2A21yfy6/pmCO0zLjnyuYI9OmYPdT84E89Mm923dYWRBHbp9DCubIndO7C+ZYnBnPxVnSxuDYzHjJY19yHEwZOpEDAABUooACAACoRAEFAABQiQIKAACgEgUUAABAJQooAACAShRQAAAAlSigAAAAKnXVSLPWIq3W0pYmllPVwJJGmQCqzVJ7A8peNDnct2CfSzPjcwvmyMT60IcyTTKl/G+Pgim0oGCfnNz9XVIwR+6xyzQeLTLSZQwYOqxAAQAAVKKAAgAAqEQBBQAAUIkCCgAAoBIFFAAAQCUKKAAAgEoUUAAAAJWmtA/Uai3SWVo6lYfsmzNb+llJ9JoCZpyQ9HjLeEn/pZHMeElfpF26HC+xsQdxrOnBHCXnNPdbbEPBHLmeVbnzUfKbNNfnqeS+YqiwAgUAAFCJAgoAAKASBRQAAEAlCigAAIBKFFAAAACVKKAAAAAqUUABAABUooACAACoNKWNNGeS6dIoM9fwU5o+9wUYqFB7M8Rco0Qp3ziyrVHnqAMy4yt7EEfuGFK+UWbJ+cg1qCxpLpk7zqaCOXLnvWSOnH0z471oxokpxQoUAABAJQooAACAShRQAAAAlSigAAAAKlFAAQAAVKKAAgAAqEQBBQAAUIk+UDPcTOrxlOtpNZPuK4aQ1Z4xS3o4LciMl/QbymXtkjhyetF/qWSOXB+okr5HveiN1ItzltOLx40+UEOlqwLK9oikRyRtkbQ5Ipb0IigAmArkMACT1YsVqJdHxIYezAMAg0AOA1CNa6AAAAAqdVtAhaTv2L7J9tJeBAQAU4gcBmBSun0L74iIWGV7L0nftX1nRFzXuUNKSikx7d7l4QCgp1pz2Db5a7f9BxQigGHU1QpURKxK/66T9HVJh42zz3kRsaS5OHN2N4cDgJ7K5bBt8tecPQcRIoAhNekCyvYc27uO3pb0Kkm39iowAOgnchiAbnTzFt5CSV+3PTrPFyPi6p5EBQD9Rw4DMGmTLqAi4h5Jf9jDWDBGrnGktH01j9ye7iv6rzqH7aD25pAlTTBzzSVLMvJPMuPzCubIxXptwRy54+xdMEdun1yjTSl/zp5fMEeuicW+mfGSBpe5Y9Akc9qhjQEAAEAlCigAAIBKFFAAAACVKKAAAAAqUUABAABUooACAACoRAEFAABQqdvPwkMfTae+R7meVdPpvgDj2irp8Zbxkv5Ld2fGS3on5ZTEkes5tKZgjlysJXPkYi35DZXbp+0xG7WgyzlK4szNUfK48Rt7qLACBQAAUIkCCgAAoBIFFAAAQCUKKAAAgEoUUAAAAJUooAAAACpRQAEAAFSigAIAAKhEWy70RK5RJo02Me1Z7Rkz15xSkuZmxncpmGNjZrwkjlxTx9wxSuYoaQyZ26ekGWfut1hJHLn7m3tcSn6T5h6XkoafGCqsQAEAAFSigAIAAKhEAQUAAFCJAgoAAKASBRQAAEAlCigAAIBKFFAAAACV6AM1QNtTb6SZdF+wndoqaVPLeK7HkzI/L5X1cDo6M76yYI5cb6S/KJgjd5ySOEYy4/sWzJE7Z8sL5licGR8pmCMnd1960UsKU4oVKAAAgEoUUAAAAJUooAAAACpRQAEAAFSigAIAAKhEAQUAAFCJAgoAAKASBRQAAEClbOsu2xdKOlbSuog4OG2bL+nLatqPjUg6PiIe7F+Y5XLNKaXhaeo4LHEAM1nPctgOam+WWdLkcHFm/PGCOe4u2Ccnd5wNBXPkfnv0oglmLxpHLi7YJ3c+FmTGS+LMHaPkscdQKVmBukjSMWO2fVDS9yLiQEnfS98DwDC6SOQwAD2WLaAi4jpJD4zZfJyki9PtiyW9rsdxAUBPkMMA9MNkr4FaGBGr0+01khb2KB4AmArkMABd6foi8ogISTHRuO2ltpfbXi491u3hAKCn2nLYNvnr0fVTHBmAYTbZAmqt7UWSlP5dN9GOEXFeRCyJiCXS7EkeDgB6qiiHbZO/5uw5pQECGG6TLaCulHRSun2SpG/0JhwAmBLkMABdyRZQti+T9GNJz7O90vY7JJ0t6ZW275J0dPoeAIYOOQxAP2T7QEXECRMMvaLHsfTEsPRWmk79qICZrGc5bIukjS3jJX18cv2ESno8/dfMgT62S36Oee3Dzz7ztuwU93z4Be07rMmHoTsz40sK5tiUGS85p7meVW2Pu1TWB+pNmfGSOHvRFws9QydyAACAShRQAAAAlSigAAAAKlFAAQAAVKKAAgAAqEQBBQAAUIkCCgAAoBIFFAAAQKVsI01MzlQ1ycw17KRZJ9AjO6u94eKGgjlyPS4PKZjj7ZlJSppPZppc3nNkpkmmJB3e5bgkHZUZX14wx9zM+LsK5vhJZjx3X0oaXP5TwT45/MYeKqxAAQAAVKKAAgAAqEQBBQAAUIkCCgAAoBIFFAAAQCUKKAAAgEoUUAAAAJXoKjHN0efpKbmeWBLnC114Qu39k+YVzJHra1QyxzGZ8ZKeRHtnxhcXzJFzfcE+ub5YJecj91uspP9SW38vSbqzyxhKjpE7Fxg6rEABAABUooACAACoRAEFAABQiQIKAACgEgUUAABAJQooAACAShRQAAAAlSigAAAAKg1VI00aIaIb/N9AX1ntzQ5Lsmlun7kFc4xkxhcXzPF4ZjzXOFKSDsiMl9yXXKPMkqagObmmoSVycZY89m1NWEuOgaHDChQAAEAlCigAAIBKFFAAAACVKKAAAAAqUUABAABUooACAACoRAEFAABQKdu9wvaFko6VtC4iDk7blkk6WdL6tNsZEXFVt8HQx2f7RQ8w9EvPctgOau8DleutJOV7/bTNPyrXT6ik71Eu1twxJGnfzPiCgjly+4wUzJGTi1PK399cT6uSPlC5c76pYI6h6tyIkhWoiyQdM872T0TEIemr6+IJAPrkIpHDAPRYtoCKiOskPTAFsQBAz5HDAPRDN9dAnWr7F7YvtL1HzyICgKlBDgMwaZMtoD4r6TmSDpG0WtI5E+1oe6nt5baXS49N8nAA0FNFOWyb/PXo+vF2AbCdmlQBFRFrI2JLRGyVdL6kw1r2PS8ilkTEEmn2ZOMEgJ4pzWHb5K85e05tkACG2qQKKNuLOr59vaRbexMOAPQfOQxAt0raGFwm6ShJC2yvlHSmpKNsHyIp1Pyh6Sl9jBEAJo0cBqAfsgVURJwwzuYL+hALhtRU9GiixxP6pWc5bJba+zjdXTDHUZnxkv5LG7ocl6SNmfGSnkS5OY4omCP3G+jOgjk2Z8aPKpjj5sz44Znxkh5guThXFsxR0icMU4ZO5AAAAJUooAAAACpRQAEAAFSigAIAAKhEAQUAAFCJAgoAAKASBRQAAEAlCigAAIBK2Uaa26OpaBw5nWxP9xWY0BNqb3S5d8EcV2fGS+Y4ODO+b8EcJcfpdo5rC+Zoa0wqSXML5sj9FiuJY3FmfKTLGKT8fcnFgKHDChQAAEAlCigAAIBKFFAAAACVKKAAAAAqUUABAABUooACAACoRAEFAABQiT5Q46Dv0cyV6/HFY48JhaTHW8ZLsmlun40FcxyeGR8pmCPXk+iYgjlyx9lUMMeGzPjigjk2Z8ZHCubIHaet/5dU9tjn+nOVzJG7r5hSrEABAABUooACAACoRAEFAABQiQIKAACgEgUUAABAJQooAACAShRQAAAAlSigAAAAKtFIE9tVc8mZdF8wxWZJmtcyXtLksO3npXyDS0m6OjP+/II5cg0sry2Y4+DM+OKCOXL3N9fAUsr/FltSMEfufJSc05y7M+O7FMzBb+yhwgoUAABAJQooAACAShRQAAAAlSigAAAAKlFAAQAAVKKAAgAAqEQBBQAAUImuEqA3UodcTyyJ87XdWnWf9MFlg45ieFwx6AAw1jmxLrvPCu3XOr5Zs7JznHvBB1rHL3/Ha7NzHKKbW8dn67HsHLO0pXV80b3rs3N89JnvbR0/3ROPZVegbO9n+xrbt9u+zfZ70vb5tr9r+6707x7ZSAFgCpG/APRLyVt4myW9LyIOknS4pHfbPkjSByV9LyIOlPS99D0ADBPyF4C+yBZQEbE6In6Wbj8i6Q5J+0g6TtLFabeLJb2uX0ECwGSQvwD0S9U1ULYXSzpU0g2SFkbE6jS0RtLCCX5mqaSlzXe7Ty5KAOgS+QtALxX/FZ7tuZK+Kum0iHi4cywiQlKM93MRcV5ELImIJdLsroIFgMkgfwHotaICyvZOapLPFyLia2nzWtuL0vgiSfnL/wFgipG/APRDyV/hWdIFku6IiI93DF0p6aR0+yRJ3+h9eAAweeQvAP1Scg3Un0g6UdIttkcbN5wh6WxJl9t+h6R7JR3fnxABYNLIXwD6IltARcT1kiZqJfWK3oYDDNZUNcnMNeykWWdvkL+wPXms4Dq91+iq1vG9tDY7x1HvPL11/LJ3XJKd4/hdv5k5SHYKaUn7cFzc0gUzOeue9qagbfgoFwAAgEoUUAAAAJUooAAAACpRQAEAAFSigAIAAKhEAQUAAFCJAgoAAKBS1YcJA+iNXJ+nXJ+okjkAbF+2aFZ2ny/oza3jvynoJXX5aSe1jl+kN2bn+Mo1J7bvUPDZ3fH09nEfPu5HXG7j7/T+/IEmwAoUAABAJQooAACAShRQAAAAlSigAAAAKlFAAQAAVKKAAgAAqEQBBQAAUIkCCgAAoBKNNIEh1IsmmTTjHEY7ZcafmJIoMDP9Wouz+9yhg1rHN2hB/kBXtw+v+MR++TmOz4y/IT+FX5HZIROnJK3404JYJ8AKFAAAQCUKKAAAgEoUUAAAAJUooAAAACpRQAEAAFSigAIAAKhEAQUAAFCJPlBAjw1L/yV6PA0j+jyhf47W97L7vE2XtI7P1mPZOV548y2t45/RX2Xn8CnRvsPR2Sk0/5BVreNxsLNzXKo/bx0/r2WMFSgAAIBKFFAAAACVKKAAAAAqUUABAABUooACAACoRAEFAABQiQIKAACgEgUUAABApWwjTdv7Sfq8pIWSQtJ5EfEp28sknSxpfdr1jIi4ql+BAsMi1yiTBpbDg/yF7clGzcvus59WtI7vqkeyc8x5w9bW8RVX7pedQ1/JjG/MT/HAyn3ad7gmP8e6T+6V32kCJZ3IN0t6X0T8zPaukm6y/d009omI+Nikjw4A/UX+AtAX2QIqIlZLWp1uP2L7DkmZsg8ABo/8BaBfqq6Bsr1Y0qGSbkibTrX9C9sX2t6jx7EBQM+QvwD0UnEBZXuupK9KOi0iHpb0WUnPkXSImld450zwc0ttL7e9XAUfUggAvUb+AtBrRQWU7Z3UJJ8vRMTXJCki1kbElojYKul8SYeN97MRcV5ELImIJdLsXsUNAEXIXwD6IVtA2bakCyTdEREf79i+qGO310u6tffhAcDkkb8A9EvJX+H9iaQTJd1i++a07QxJJ9g+RM2fBo9IOqUvEQLA5JG/APRFyV/hXS/J4wzRMwXbJfo8bSvfF2twyF/YnizU2uw+j2hu6/jvtHN2jn+68ujW8UN1c+u4JOmi9uEDX/Dz7BTP069ax68+7mXZOfbSuuw+E6ETOQAAQCUKKAAAgEoUUAAAAJUooAAAACpRQAEAAFSigAIAAKhEAQUAAFCJAgoAAKBSSSdyADNUvglmvmlofp9BttIEth/X6cjsPr/S81rHN2pedo6fvry9QeUnr8k39o+PjdfftsPh2SmajwJv4WsjO8Vpp38ks8dXJxxhBQoAAKASBRQAAEAlCigAAIBKFFAAAACVKKAAAAAqUUABAABUooACAACoRB8oYDuW6+GU6xNVMgeAqfEsjWT3eb2uaB1fqLXZOV5w7z2t46fok9k5fH2mR9Mu2Smy+8RNmV5Tkj6i0woOND5WoAAAACpRQAEAAFSigAIAAKhEAQUAAFCJAgoAAKASBRQAAEAlCigAAIBKFFAAAACVaKQJTFO5Jpe9aHBJk8xOO0j6g5bxzQVzzM+MP1wwR8lxhsGuBfvk7stvehBH22NWepxe/Krs/+N2pf5ddp+//e3ftI4/NLJ3/kC/bh9+q87PTnHm3e1NLpc9Nx+GXts+7PWZZp2SXqVvZPaYuCkoK1AAAACVKKAAAAAqUUABAABUooACAACoRAEFAABQiQIKAACgEgUUAABAJUfk+yT07GD2ekn3dmxaIGnDlAUwedMlTmn6xEqcvTessT4zIvYcdBDdGid/ScN7zscizt6aLnFK0yfWYY1zwvw1pQXU7x3cXh4RSwYWQKHpEqc0fWIlzt6bTrHOFNPlnBNnb02XOKXpE+t0ibMTb+EBAABUooACAACoNOgC6rwBH7/UdIlTmj6xEmfvTadYZ4rpcs6Js7emS5zS9Il1usT5pIFeAwUAADAdDXoFCgAAYNoZWAFl+xjbv7R9t+0PDiqOHNsjtm+xfbPt5YOOp5PtC22vs31rx7b5tr9r+6707x6DjDHFNF6cy2yvSuf1ZtuvGWSMKab9bF9j+3bbt9l+T9o+VOe0Jc6hO6czFfmre+Sv3iJ/Tb2BvIVne5akX0l6paSVkm6UdEJE3D7lwWTYHpG0JCKGrj+F7ZdK2iTp8xFxcNr295IeiIizU2LfIyJOH8I4l0naFBEfG2RsnWwvkrQoIn5me1dJN0l6naS3a4jOaUucx2vIzulMRP7qDfJXb5G/pt6gVqAOk3R3RNwTEb+T9CVJxw0olmkrIq6T9MCYzcdJujjdvljNf8yBmiDOoRMRqyPiZ+n2I5LukLSPhuyctsSJqUH+6gHyV2+Rv6beoAqofSSt6Ph+pYb3BIak79i+yfbSQQdTYGFErE6310haOMhgMk61/Yu0RD7wpfpOthdLOlTSDRriczomTmmIz+kMQv7qn6F9ro1jaJ9r5K+pwUXkeUdExB9JerWkd6fl3Gkhmvdnh/XPLD8r6TmSDpG0WtI5gw3nKbbnSvqqpNMi4uHOsWE6p+PEObTnFAND/uqPoX2ukb+mzqAKqFWS9uv4ft+0behExKr07zpJX1ezfD/M1qb3mEffa1434HjGFRFrI2JLRGyVdL6G5Lza3knNk/oLEfG1tHnozul4cQ7rOZ2ByF/9M3TPtfEM63ON/DW1BlVA3SjpQNvPsr2zpDdJunJAsUzI9px0kZtsz5H0Kkm3tv/UwF0p6aR0+yRJ3xhgLBMafZNsvHIAAADeSURBVEInr9cQnFfblnSBpDsi4uMdQ0N1TieKcxjP6QxF/uqfoXquTWQYn2vkr6k3sEaa6U8UPylplqQLI+LDAwmkhe1nq3nVJkk7SvriMMVp+zJJR6n5FOu1ks6UdIWkyyXtr+aT44+PiIFeADlBnEepWaoNSSOSTul4n34gbB8h6QeSbpG0NW0+Q83780NzTlviPEFDdk5nKvJX98hfvUX+mnp0IgcAAKjEReQAAACVKKAAAAAqUUABAABUooACAACoRAEFAABQiQIKAACgEgUUAABAJQooAACASv8fCi5jD+vsRMUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----m = 2/3\n",
      "Partial Wasserstein distance (m = 2/3): 0.252736149616858\n",
      "Entropic partial Wasserstein distance (m = 2/3): 1.407282181449262\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEtCAYAAADHtl7HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZRdZZXn8d82CYGEFJAg6YQXE+moMGAHkiZOGyVtGwUEwdWK2kKTGSXiyEKmdUaHnlkkttrYA74tHBWUAUXFF0TRRhsaBQytkYSVIUhiB6GQhJBggiQkISZhzx/nFN6UVWc/T91z38L3s1ZWqu4+d59dp+o+d99Tp/Y1dxcAAADSvaDTBQAAAPQaGigAAIBMNFAAAACZaKAAAAAy0UABAABkooECAADIRAP1PGBml5jZFxO3vdbMPtLqmnqNmb3DzG7tdB1ArzOzo8zsaTMb1elaBmOtbN7zaa2kgeoCZtZvZjvKRWVD+cA8cIS55pnZ2sbb3P1j7v6uJmscXdY3p+G2d5iZD3Hb6mb21Q5mtsjMrk/d3t2/6u6va2VNQF0GrSkD/65MvO8dZtbUelHF3X/j7ge6+55W7SMFa2Ua1srh0UB1jzPc/UBJJ0qaLel/5iYws9G1V1Vy992Sfibp1Q03v1rS6iFuu6tVdaRq5bEAesQZZaMy8O/COpLuC48t1so/2Be+n51CA9Vl3H2dpB9KOk6SzOw/mdkqM9tqZg+Z2bsHth14BWVmHzSzxyV9vbzv1IZXnVMHv4Iws2+Z2eNm9pSZ3WVm/yGxvLu09wLwKkkfH+K2u8r9nGRmPzOz35nZejO70sz2K2NmZp80s41mtsXMVprZwNd8mpk9UH7N68zsAw21n25mK8qc/2ZmL2+I9ZfH4j5J28pXgh8sc2w1s1+Z2V+Z2SmSLpH01vIY/b/y/geZ2ZfKWteZ2UcGfs1gZgvMbEnDvtzMLjCzNWUtnzUzSzyOQMcM/Cyb2eVm9qSZPWxmp5axj6p4DF/ZeNaq/Hl/r5mtkbSmvO0vzOyech25x8z+omEfd5jZP5rZL8rH9/fMbGIZm1bmG11+PtHM/q+ZPVbW892Kuu8u15GnzGy1mf1VQ5y1krWyvdydfx3+J6lf0mvLj4+U9EtJ/1B+/gZJR0sySSdL2i7pxDI2T9JuFQ/MsZIOKG9bOyj/IknXN3z+nyVNKO/zKUkrGmLXSvrIMHWeLGmzisb7UEmPSBonaUPDbS7pqHL7WZJeIWm0pGmSVkm6uIy9XtJySQeXX9sxkqaUsfWSXlV+fEjD13uCpI2S5kgaJem88tiNbTiOK8pjeICkl0p6VNLUMj5N0tFDHZPytpskfUHSeEmHSfqFpHeXsQWSljRs65J+UNZ/lKQnJJ3S6Z8l/vHPfe81ZYjYAkm7JJ1fPo7eI+kxSVbG75D0rkH3cUm3SZpYPrYmSnpS0rnl4/vt5eeTGnKsU/FCcLykGwceb+Xj0CWNLj//Z0nfKB/rYySdXFH3bkn/tdzurZKekjSxjLNWsla29R9noLrHd83sd5KWSLpT0sckyd3/2d1/7YU7Jd2q4pXLgGclXeruO919R8qO3P0ad9/q7jtVPDj+zMwOSrjrUhWLwPFlDUvcfbukhxtu63f335T7We7uP3f33e7er+IBd3KZa5eKhellKhbuVe6+viF2rJn1ufuT7n5veftCSV9w96Xuvsfdr5O0U8XCM+Az7v5oeSz2qFj4jjWzMe7e7+6/HuoLM7PJkk5TsWhtc/eNkj4p6W0Vx+Myd/9d+fX+RNLM+BACbfPd8hX/wL/zG2KPuPvVXlyHdJ2kKZImB/n+0d03l4+tN0ha4+5fKR/fX1fxK6ozGrb/irvf7+7bJP0vSWfboAvHzWyKpFMlXVA+1neV69xwNkr6VLndNyT9qqyFtZK1su1ooLrHWe5+sLu/yN3/y8AD3MxONbOfm9nmssE6TcWrlwFPuPszqTsxs1FmdpmZ/drMtqh4JaJBOYdU7ucXKk5Dv1rST8vQkobbnvudvpm9xMx+UJ4C36KiKTy0zPVjSVdK+qykjWZ2lZn1lXf96/LrfMTM7jSz/1je/iJJ7298UlDxCmpqQ5mPNtT7oKSLVSx8G83sBjNr3LbRi1S8ql3fkPsLKl5dDefxho+3SxrRhf9AiwysKQP/rm6IPfezWz6xS/HP76MNH09VcVal0SOSDh9m+0dUPL4GrzNHStrs7k8G+x6wzr04rdGQd6rEWsla2X40UF3MzMaqOPV9uaTJ7n6wpFtUnMYd4IPuNvjzwf5G0pmSXivpIBWnajUoZ5WB3+2/Sn9YFH7acFvjRZGfU/GqdIa796n4Xfpz+3H3z7j7LEnHSnqJpP9W3n6Pu5+p4gH5XUnfLO/yqKSPDnpSGFe++n0ubWOx7v41d5+r4kHvKk7h/9F2Ze6dkg5tyN3n7qnXPAD7iuHWkMbbH1PxmGp0lIpf2w04clBsl6TfDrrPo5ImmtnBibUdPuj6maMkPcZayVrZCTRQ3W0/FadVn5C024oLPaM/D90gaVLFaeYJKn74N6k4xfyxzJrukvSXKhbHB8rb7lZxPcFM7b0oTJC0RdLTZvYyFddaSJLM7M/NbI6ZjZG0TdIzkp41s/2s+PPeg9x9V3n/Z8u7XS3pgvJ+ZmbjzewNZjZhqELN7KVm9ppycX1G0o6GXBskTTOzF0hSeUr8VklXmFmfmb3AzI42s5OHyg3swzZIenGwzS2SXmJmf1NegPxWFU/uP2jY5hwzO9bMxkn6sKRv+6DRBeXj7oeS/o+ZHWJmY8ys8ULrwQ6TdFG53VtUXA90i1grWSs7gAaqi7n7VkkXqXhV8aSKV0Q3B/dZreIvTB4qT68OPg37ZRWnvdepeFD/PLOsf1PxamzpwKl0d/+tioVro7uvadj2A2XNW1U8oL/REOsrb3uyrGeTpP9dxs6V1F+eyr5A0jvK/SxTceHrleX9HlRxweJwxkq6TMWr3sdVLL7/o4x9q/x/k5kNXDfwtyoW4gfK/N9WcW0I0Iu+b3vPgbop8X6flvRmK/4i7jNDbeDumySdLun9Kh67/13S6eVaMOArKi60flzS/irWsqGcq+Ls1GoV1zhdXFHbUkkzVDymPyrpze6+ibWStbITBv7qAgCAWpjZHSr+citpqndizgUq/jpwbl05gWZwBgoAACATDRQAAEAmfoUHAACQiTNQAAAAmWigAAAAMjX1LsxWvNHgp1W8184X3f2y6u3HefF2OK0zRevDbdY///7aEuig9b919xd2uoqh5KxhB5v5cKOZU40bUx3ftivOMSqIPxvEU4xNGBX5THD1R8qTy+6kalov+nLruNAlOltRx/dtNc9tLTD8+jXiBsqK9zT6rKT5ktZKusfMbnb3B4a/18Eq3qKndRZqcbjN4hbXAKDR4sFv+dEVctewqSoGGzVjVtBGLn0szjHkJMQGSW/yFpixX7zNqp3V8YkJ+9mcVE3rRU+EdTR6BwTxOr5vr+C5rQWGX7+a+RXeSZIedPeH3P33km5QMfYeAHoBaxiAEWumgTpce79Z5Frt/UaSANDNWMMAjFhT10ClMLOFeu73dsO95RAAdJ/G9etPOlwLgO7SzBmoddr73baP0N7vxC1Jcver3H22u88u3o8RALpCuIY1rl+HtLU0AN2umQbqHkkzzGy6me0n6W0K3rwRALoIaxiAERvxr/DcfbeZXSjpX1T8Ze017v7L2ioDgBZiDQPQjKaugXL3WyTdUlMttVisSztdgiTp0qRxCt1RK/B8lbOGjTtQmjWrYoPDEpI8Ux2e87KEHHuCeLPDqiRpU7zJnLHBBtsS9rN/SjGB8UE84WsJj1k0byE6Fik5Un5+It+qIQeSMYkcAAAgEw0UAABAJhooAACATDRQAAAAmWigAAAAMtFAAQAAZKKBAgAAyEQDBQAAkKnlbyb8fMWQzL1Fg0U5Xuh2O5+WHr5z+PjGhBxzXlMdf+DHcY4jgqGNq3YmFBJIGeh53+rq+PET4xwro+GSCXYE8QkJOaIy+oL4roR9TBtVHX8wGpAqnrC7DWegAAAAMtFAAQAAZKKBAgAAyEQDBQAAkIkGCgAAIBMNFAAAQCYaKAAAgExdNVYimhUkMS+oV/F9Q68bO06afmzFBssSkryvOnzsowk5jqsOT7wpThEu/BfEOY7/cHXc3hLneHlUa8Iz1IbHquOTZ8Q57ltTHX/5YdXxTZvifUw6rTq++/txjtHBLCklzJJCfTgDBQAAkIkGCgAAIBMNFAAAQCYaKAAAgEw0UAAAAJlooAAAADLRQAEAAGSigQIAAMjUVYM0GbaIVouGtfIziGGZpIpBhtODYYuSpGhQ5viEHLurwzNS6ohsjDexScEG2xL2MzGlmGqTow3GxjmOj+oI4pMS9hENuZx8ZEKOSMogVtSGM1AAAACZaKAAAAAy0UABAABkooECAADIRAMFAACQiQYKAAAgEw0UAABApq6aAwVUqWOGE3OeMFL+jLRr9fDxHc/EOfpuqY5vXxPnGBfMaNryVJwj0nd7vM32x6rj45Y2nyPF7mAu1uiE4xHmCOYrRfeXpL4V1fEtm+Ico3nG7ipNfTvMrF/SVhUjwna7++w6igKAdmANAzBSdfSzf+nuv60hDwB0AmsYgGxcAwUAAJCp2QbKJd1qZsvNbGEdBQFAG7GGARiRZn+FN9fd15nZYZJuM7PV7n5X4wblolQuTAc1uTsAqFXlGta4fh1lnSoRQDdq6gyUu68r/98o6SZJJw2xzVXuPru4OHNcM7sDgFpFa1jj+nUoFzwAaDDiJcHMxpvZhIGPJb1O0v11FQYArcQaBqAZzfwKb7Kkm8xsIM/X3P1HtVQFAK3HGgZgxEbcQLn7Q5L+rMZa0KWiAZZSewZUMgQTdcpdw7buke6oGMq4ISHHOSur4/+6Lc4xOdgmYRZn6I3B0EdJ+uHO6vgxCYXUUWs0w/KAoE5J2tJkjoQ5mjo8GBq6LiGHEr4WtA+/1QcAAMhEAwUAAJCJBgoAACATDRQAAEAmGigAAIBMNFAAAACZaKAAAAAyNfteeHgeSJm/FM2KYoYTel3fWGn+UcPH70gZavTp6vAZ74pT2Ozq+OZb4xzRwt/34TjHSR+sjk8/J84x+frq+OhRcY5Ve6rjcxLegvXmivlekjQ3uP/aeBeadXJ1fOmdcQ6esLsLZ6AAAAAy0UABAABkooECAADIRAMFAACQiQYKAAAgEw0UAABAJhooAACATDRQAAAAmZjLtY+LBlxK9Qy5ZFAm9nU7d0oPVwzLXJWQY97d1fF7N8c5TlxWHU+Z5zkmiM//eZxjZRCfHtQpSfdFGwRDMiVpXRAfEwzJlKQHg3hfEN8Q70IvX1EdT/n54Qm7u3AGCgAAIBMNFAAAQCYaKAAAgEw0UAAAAJlooAAAADLRQAEAAGSigQIAAMjEWIkRiuYrdctcpG6pA+h1YydI0/98+Ph7UgYwzawOz5qTkOPF1eGLliTkiMyON3lj9PW+Ms4xf/+kaqpFc54Oi1PMebTJHNvifUTf+wXBnKgU56b8DKI2nIECAADIRAMFAACQiQYKAAAgEw0UAABAJhooAACATDRQAAAAmWigAAAAMtFAAQAAZDJ3r97A7BpJp0va6O7HlbdNlPQNSdMk9Us6292fjHY21cwXVsQZ+gjsixYvd/eE0YytUdcadpyZf7MivjmhlrnBMMXlCcMUJwfx/oQ6InOnxtssfaw6/qej4hz9e9LqqbIjiB/Qhhy7E/YRzeJcl5BjTBB/Bc+hLTD8+pVyBupaSacMuu1Dkm539xmSbi8/B4BudK1YwwDULGyg3P0u/fGLqzMlXVd+fJ2ks2quCwBqwRoGoBVGeg3UZHdfX378uOKzygDQTVjDADSl6YvIvbiIatgLqcxsoZktM7Nl25vdGQDUrGoNa1y/Uq5xAvD8MdIGaoOZTZGk8v+Nw23o7le5+2x3nz1uhDsDgJolrWGN69fEtpYHoNuNtIG6WdJ55cfnSfpePeUAQFuwhgFoSthAmdnXJf1M0kvNbK2ZvVPSZZLmm9kaSa8tPweArsMaBqAVwjlQte7MprpUNQmqPS7V4so486iAOnV2DlRdZh9kvmxuxQbHJST5cRCfnpBj2AsmSicn5IgGF92SkCP6jq5KyDEjYZvIpCD+m4QcwXwurQniByXsY2UQPyEhR8Cu4Lmrfs3NgQIAAEADGigAAIBMNFAAAACZaKAAAAAy0UABAABkooECAADIRAMFAACQiQYKAAAg0+hOF9AJvTIoMxr4KfXO1wL0uk1bpGsrBkxuSBg++Z6x1fFvLItzRO/J9+CdcY7IGQnb/GhFdfz4hByr7k6pptqOIN6XkGPzt6rjE4L7R3NJJenwIN4fDVlF1+EMFAAAQCYaKAAAgEw0UAAAAJlooAAAADLRQAEAAGSigQIAAMhEAwUAAJDpeTkHqle0a8ZTNG+KWVOANGm8tGDm8PFNP49z9H2zOn7+3ycUckx1eEMw00iSRo+qjk+6Is5x7D8FG7wpzjE/YXZWZMtj1fG+6XGOh1dXx6cHOXZtjvcx5jXV8e23JuQInrE/9FScA/XhDBQAAEAmGigAAIBMNFAAAACZaKAAAAAy0UABAABkooECAADIRAMFAACQiQYKAAAgE4M0O6hbBlgyKBNIsEtSxdDGlXviFPN2Vsc9GAopSXZYdfzeOIVGB7XO3xbneDiodXrCcMkND8fbRNYE8Tnr4xzRMRsX1Lk23oVmbayOr0w45jxhdxfOQAEAAGSigQIAAMhEAwUAAJCJBgoAACATDRQAAEAmGigAAIBMNFAAAACZwrESZnaNpNMlbXT348rbFkk6X9IT5WaXuPstrSpyX1XH/KVumSUFdKu61rAdv5fuq5gHlDJ/ad6t1fE7E2YnnXh3dTxpDlQQn397nCMoQ9PvjHMsiTcJrQviO56Kc6xssoZgxJMk6YjggC1NyDEmpRi0TcoZqGslnTLE7Z9095nlP5onAN3qWrGGAahZ2EC5+12SEl4XAUD3YQ0D0ArNXAN1oZndZ2bXmNkhtVUEAO3BGgZgxEbaQH1O0tGSZkpaL+mK4TY0s4VmtszMlknbR7g7AKhV0hrWuH492c7qAHS9ETVQ7r7B3fe4+7OSrpZ0UsW2V7n7bHefLY0baZ0AUJvUNaxx/eIUFYBGI2qgzGxKw6dvknR/PeUAQOuxhgFoVsoYg69LmifpUDNbK+lSSfPMbKYkl9Qv6d0trBEARow1DEArmLu3b2c21aWFTeVg7hHQaxYvL36F39uOMvMPVMQnJuQ4IIinzBPaFcRT6ojsrqGOHQk56phrFB3TlDqiYxb9CWdUgyRtCeJ9CTmiY76Q578WGH79YhI5AABAJhooAACATDRQAAAAmWigAAAAMtFAAQAAZKKBAgAAyEQDBQAAkIkGCgAAIFM4ibzbtGNQZjSss111AOgeh+0vXTRt+PiS1XGOud8JNrggoZCZ1eHbbo1TRAMs510e51hbNVVU0hEL4hybvlIdH5PwDLVqZ3X8pITJot8PJmWePao6vm5PvI+Xv6Y6vvzHcY46Bo+iPpyBAgAAyEQDBQAAkIkGCgAAIBMNFAAAQCYaKAAAgEw0UAAAAJlooAAAADL13ByodmDGU29ifhdaaoykqcOH5wbziCRJTwXx6Qk5JlWH56fkiCR8LUccF2wwNs4x6YSkairNiY5pwhyoNz4abHBYdXjStngf0fdt1oyEHME8KiXMIkN9OAMFAACQiQYKAAAgEw0UAABAJhooAACATDRQAAAAmWigAAAAMtFAAQAAZKKBAgAAyMQgTXSFOoZgMiQTrbR9q7T8x8PH70zI8Xe3V8fvWBrnOH5ZdfyLe+IcY4L4390S57j+/ur4ORvjHDcmbBNZF8SPScjxiyD+sseq45sT9vHGNdXx7yTk4Am7u3AGCgAAIBMNFAAAQCYaKAAAgEw0UAAAAJlooAAAADLRQAEAAGSigQIAAMhk7t6+ndlUlxa2bX8AusHi5e4+u9NVNMsmzHbNrhjCdEec40RfUhm/9/i5cZLXBvFP1bCm9++Ot5kZTJP6SMJ+UraJPB7EX5aQY/X26vifjGuuBkm6IIh/PiFHaFEdSbCX4dev8AyUmR1pZj8xswfM7Jdm9r7y9olmdpuZrSn/P6TusgGgGaxfAFol5Vd4uyW9392PlfQKSe81s2MlfUjS7e4+Q9Lt5ecA0E1YvwC0RNhAuft6d7+3/HirpFWSDpd0pqTrys2uk3RWq4oEgJFg/QLQKlkXkZvZNEknSFoqabK7ry9Dj0uaPMx9FprZMjNbJgW/ZwaAFml6/dr1RFvqBNAbkhsoMztQ0o2SLnb3LY0xL65EH/LKRXe/yt1nFxdhBRfiAUAL1LJ+jXlhGyoF0CuSGigzG6Ni8fmquw+8afQGM5tSxqdIquF9tQGgXqxfAFoh5a/wTNKXJK1y9080hG6WdF758XmSvld/eQAwcqxfAFpldMI2r5R0rqSVZraivO0SSZdJ+qaZvVPSI5LObk2JADBirF8AWoJBmgBabB8ZpGnTXbq0YovNCVmiPu37CTkmBvF1CTkiJydsc3cQn5GQoz9hm8iOIH5AQo4tQbyvyRqk4o8/q9TxfYu+DuRrYpAmAAAA9kYDBQAAkIkGCgAAIBMNFAAAQCYaKAAAgEw0UAAAAJlooAAAADKlDNIEABw6STprwfDxG+IU52y9ujJ+/anviZOcFcQvjFNEK/+Ld/wyTPHQ8RdVb3BlQh0fCOIpz1Crg/grEnL8KIhHU8zWJuwj+lovT8gRHY+1ixKSoC6cgQIAAMhEAwUAAJCJBgoAACATDRQAAEAmGigAAIBMNFAAAACZaKAAAAAy0UABAABkYpAmAKQ4RNKbK+IJwxTfpS9Wxq8/6/wwx4HnPFEZf/rbL4wLCVb+v9WXwxSLTvl4ZfzEk5eEOe49fW64TeiIID4vIcfvgvgpQfzBeBdjFmypjO9a1hcniZ6xr49ToD6cgQIAAMhEAwUAAJCJBgoAACATDRQAAEAmGigAAIBMNFAAAACZaKAAAAAymbu3b2c21aWFbdsfgG6weLm7z+50Fc0yO8alayu2WBMnmXlOdXzFbQmV/GkQvzMhRzBQ6MCgTkl6+sZgg/kJddyRsE2ker6SNDkhx7ogPjGI70jYx4lB/N6EHJFVNeTA3oZfvzgDBQAAkIkGCgAAIBMNFAAAQCYaKAAAgEw0UAAAAJlooAAAADLRQAEAAGSigQIAAMgUTFOTzOxISV9WMY3MJV3l7p82s0WSzpf0RLnpJe5+S6sKBYBc9a5fv5e0tiK+OS5oRbTBhjhHKBosmeDplDqirzdlqGPCMQulDLGMNHvMUmroD+J1HAu0U9hASdot6f3ufq+ZTZC03MwGxuV+0t0vb115ANAU1i8ALRE2UO6+XtL68uOtZrZK0uGtLgwAmsX6BaBVsq6BMrNpkk6QtLS86UIzu8/MrjGzQ2quDQBqw/oFoE7JDZSZHSjpRkkXu/sWSZ+TdLSkmSpe4V0xzP0WmtkyM1smba+hZADIU8/6VcO1RQD2GUkNlJmNUbH4fNXdvyNJ7r7B3fe4+7OSrpZ00lD3dfer3H128W7G4+qqGwCS1Ld+9bWvaABdL2ygzMwkfUnSKnf/RMPtUxo2e5Ok++svDwBGjvULQKuk/BXeKyWdK2mlmQ38Ee4lkt5uZjNV/Glwv6R3t6RCABg51i8ALZHyV3hLJNkQIWY+Aehq9a5f4yXNqYivjFNcEMQ/f2JCHccG8ZSZRGOqw3MnxymWHFMdP7DqWJWeTnkNH4nmJ6X80WU0s2pakzVI0msTtmlWHXPEkIpJ5AAAAJlooAAAADLRQAEAAGSigQIAAMhEAwUAAJCJBgoAACATDRQAAEAmGigAAIBMdUwxA4B931H7SZccMXz88xWx0lWfO7cyvvDgr8R1vC2In/PGOEew8p/z06vDFNe/5fzK+MQb1oU5Nr9rVvUGKc9Q0ZvwzEvIcUMwnDTK0Z+wj8uC+Afmxzn2D+L/endCIagLZ6AAAAAy0UABAABkooECAADIRAMFAACQiQYKAAAgEw0UAABAJhooAACATObu7duZTXVpYdv2B6AbLF7u7rM7XUWzzF7s0ocrttickOU9QfyrCTkmB/E1CTkiZyRs86MgfkxCjjpq3RHE+xJyRN+7CU3WIEnTgnh/Qo5Iys8g8gy/fnEGCgAAIBMNFAAAQCYaKAAAgEw0UAAAAJlooAAAADLRQAEAAGSigQIAAMhEAwUAAJBpdKcLAIBeMGbWeE1eNmfY+NqPzwhz+OutMm4vjAcy/vXh366M37j4ojBHtPKv+PuXhClm/su/V8YXvf6DYY5FSz8ebhNaFsRP3xXn+NSY6vjbgvvfH+/i9Hd+qzL+gy8kfN/2D+ILFsU5UBvOQAEAAGSigQIAAMhEAwUAAJCJBgoAACATDRQAAEAmGigAAIBMNFAAAACZzN3btzOzJyQ90nDToZJ+27YCRq5X6pR6p1bqrF+31void39hp4to1hDrl9S9x3ww6qxXr9Qp9U6t3VrnsOtXWxuoP9q52TJ3n92xAhL1Sp1S79RKnfXrpVr3Fb1yzKmzXr1Sp9Q7tfZKnY34FR4AAEAmGigAAIBMnW6grurw/lP1Sp1S79RKnfXrpVr3Fb1yzKmzXr1Sp9Q7tfZKnc/p6DVQAAAAvajTZ6AAAAB6TscaKDM7xcx+ZWYPmtmHOlVHxMz6zWylma0ws2WdrqeRmV1jZhvN7P6G2yaa2W1mtqb8/5BO1ljWNFSdi8xsXXlcV5jZaZ2ssazpSDP7iZk9YGa/NLP3lbd31TGtqLPrjum+ivWreaxf9WL9ar+O/ArPzEZJ+ndJ8yWtlXSPpLe7+wNtLyZgZv2SZrt7182nMLNXS3pa0pfd/bjytn+StNndLysX9kPc/YNdWOciSU+7++WdrK2RmU2RNMXd7zWzCZKWSzpL0gJ10TGtqPNsddkx3RexftWD9aterF/t16kzUCdJetDdH3L330u6QdKZHaqlZ7n7XZI2D7r5TEnXlR9fp+IHs6OGqbPruPt6d7+3/HirpFWSDleXHdOKOtEerF81YP2qF205l98AAAIQSURBVOtX+3WqgTpc0qMNn69V9x5Al3SrmS03s4WdLibBZHdfX378uKTJnSwmcKGZ3VeeIu/4qfpGZjZN0gmSlqqLj+mgOqUuPqb7ENav1unax9oQuvaxxvrVHlxEHpvr7idKOlXSe8vTuT3Bi9/PduufWX5O0tGSZkpaL+mKzpbzB2Z2oKQbJV3s7lsaY910TIeos2uPKTqG9as1uvaxxvrVPp1qoNZJOrLh8yPK27qOu68r/98o6SYVp++72Ybyd8wDv2ve2OF6huTuG9x9j7s/K+lqdclxNbMxKh7UX3X375Q3d90xHarObj2m+yDWr9bpusfaULr1scb61V6daqDukTTDzKab2X6S3ibp5g7VMiwzG19e5CYzGy/pdZLur75Xx90s6bzy4/Mkfa+DtQxr4AFdepO64LiamUn6kqRV7v6JhlBXHdPh6uzGY7qPYv1qna56rA2nGx9rrF/t17FBmuWfKH5K0ihJ17j7RztSSAUze7GKV22SNFrS17qpTjP7uqR5Kt7FeoOkSyV9V9I3JR2l4p3jz3b3jl4AOUyd81ScqnVJ/ZLe3fB7+o4ws7mSfipppaRny5svUfH7+a45phV1vl1ddkz3VaxfzWP9qhfrV/sxiRwAACATF5EDAABkooECAADIRAMFAACQiQYKAAAgEw0UAABAJhooAACATDRQAAAAmWigAAAAMv1/vxOFT+Vr1yMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "C1 = sp.spatial.distance.cdist(xs, xs)\n",
    "C2 = sp.spatial.distance.cdist(xt, xt)\n",
    "\n",
    "# transport 100% of the mass\n",
    "print('-----m = 1')\n",
    "m = 1\n",
    "res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m, log=True)\n",
    "res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C2, p, q, 10,\n",
    "                                                          m=m, log=True)\n",
    "\n",
    "print('Wasserstein distance (m = 1): ' + str(log0['partial_gw_dist']))\n",
    "print('Entropic Wasserstein distance (m = 1): ' + str(log['partial_gw_dist']))\n",
    "\n",
    "pl.figure(1, (10, 5))\n",
    "pl.title(\"mass to be transported m = 1\")\n",
    "pl.subplot(1, 2, 1)\n",
    "pl.imshow(res0, cmap='jet')\n",
    "pl.title('Wasserstein')\n",
    "pl.subplot(1, 2, 2)\n",
    "pl.imshow(res, cmap='jet')\n",
    "pl.title('Entropic Wasserstein')\n",
    "pl.show()\n",
    "\n",
    "# transport 2/3 of the mass\n",
    "print('-----m = 2/3')\n",
    "m = 2 / 3\n",
    "res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m, log=True)\n",
    "res, log = ot.partial.entropic_partial_gromov_wasserstein(C1, C2, p, q, 10,\n",
    "                                                          m=m, log=True)\n",
    "\n",
    "print('Partial Wasserstein distance (m = 2/3): ' +\n",
    "      str(log0['partial_gw_dist']))\n",
    "print('Entropic partial Wasserstein distance (m = 2/3): ' +\n",
    "      str(log['partial_gw_dist']))\n",
    "\n",
    "pl.figure(1, (10, 5))\n",
    "pl.title(\"mass to be transported m = 2/3\")\n",
    "pl.subplot(1, 2, 1)\n",
    "pl.imshow(res0, cmap='jet')\n",
    "pl.title('Partial Wasserstein')\n",
    "pl.subplot(1, 2, 2)\n",
    "pl.imshow(res, cmap='jet')\n",
    "pl.title('Entropic partial Wasserstein')\n",
    "pl.show()"
   ]
  }
 ],
 "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": 0
}