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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZBc13Wff6e36enZgJnBvg24yaZoiaTGtGRZihbKlhjZpFOOizKsomynkLgsR45UJdFiVQikChXGduy4YpdciMSEKSGWVFoilouxRMaSt1gUhxRBEgQXCBwQOwYYALP3evLHua/vnenXM93Tr/u97j5fVVd333fffbd7Zs6cd1ZiZiiKoijtSyzsDSiKoiiNoYJcURSlzVFBriiK0uaoIFcURWlzVJAriqK0OYkwLjo6OspjY2NhXFpxePbZZy8z86ag1tOfa3QI+merRJtQBPnY2BgmJibCuLTiQESnglxPf67RIeifrRJt1LSiKIrS5qggVxRFaXNUkCuKorQ5gQhyIvp3RHSMiF4ior8konQQ6yqKoihr07AgJ6IdAP4tgHFmvg1AHMD9ja6rKIqi1EZQppUEgF4iSgDIADgX0Lodz4EDYe9AUZR2p2FBzsxnAfwRgDcBnAdwnZm/u3IeEe0nogkimpiammr0sh3DwYNh70BpOUeOAGNjQCwmz0eOhL0jpc0JwrSyEcC9APYC2A6gj4h+feU8Zj7MzOPMPL5pk+YpKF3KkSPA/v3AqVMAszzv36/CXGmIIEwrdwN4g5mnmDkP4JsAfjaAdTuWAwcAInkA9rWaWbqAhx4CFhaWjy0syLiirJMgBPmbAN5JRBkiIgAfBHA8gHU7lgMHRBnzenp4r1WQdwFvvlnfuKLUQBA28qcBfB3AcwBeNGsebnRdRelIdu+ub1xRaiCQqBVmfpiZf4KZb2PmjzNzNoh1u4GHH65+rJqGrpp7G3PoEJDJLB/LZGRcUdaJZnaGzGpCuVpEi0a6tDH79gGHDwN79ohjZM8eeb9vX9g7U9oYFeSKsh4aCSHctw+YnARKJXlWIa40iAryiFEtouV979NIl8jQaAihxpErAaOCPGJUi2h53/taG+lCRLuI6HtE9LKpo/Op5lypDWkkhFDjyJUmoIK8TQjBLl4A8BlmvhXAOwH8DhHd2vJdRJFGQgg1jlxpAirII0y1iJbVIl2CgpnPM/Nz5vUsJDdgR/Ov3AY0EkKoceRKE1BBHhGqmUj87OKthojGANwB4OnWXz2CNBJCqHHkShNQQR4R/EwnUcgAJaJ+AN8A8HvMPONzvPuKoTUSQqhx5EoTUEGuVIWIkhAhfoSZv+k3p2uLoa03hFDjyJUmoII8ROopntUKu7iLqZvzJQDHmfmPW3v1DkfjyJWAUUEeIvWYTkKIF383gI8D+AARPW8e97R8F+2CxoYrIZIIewNKNGHmfwAQgmu1DfFiw72wQi82HFBtW2kJqpFHhFaYTjQLtElobLgSMirII0IrhKwW22oSGhuuhIwK8g5ANe2Q0dhwJWRUkHcAq2na2lauBWhsuBIyKsg7nCgkFXU8GhuuhIwK8jZhpeBVTTtiaGy4EiIqyNuEleaT9WjarU4qUhSlNagg7yJUW68BTexR2hAV5BGmVvOJatoBoU0flDZFBXnIrKUl12I+UU07IKKe2KN3C0oVVJCHzGqhg5rA02KinNijdwvKKqggbxPUfNICopzYE/W7BSVUAhHkRLSBiL5ORK8Q0XEielcQ63Yqq9m+qx1TWsB6EntaZe6I8t2CEj7M3PADwGMA/pV5nQKwYbX573jHO1gRgPUdC+bamOAAfv7cST/XL3+Zec8eZiJ5/vKXV5+byXiuC3lkMqufs15GRpZfx3uMjPhOD/pnq49oPxrWyIloCMB7IU0IwMw5Zr7W6LqKEgr1JPaouUOJCEGYVvYCmALw34noR0T0RSLqWzmpK3s71sBqtm+1i0ecVpo7pqfrG1e6iiAEeQLAnQC+wMx3AJgH8ODKSdytvR3XYLXQQQ0rjDitdI5G2RGrhE4QgvwMgDPM/LR5/3WIYFeqoAK6Qwi66qHrOB0dlYfnRL3nHq2wqFSlYUHOzBcAnCait5ihDwJ4udF1OxmND+8Qgqx6uDJO/MoVebCJGX/sMeCBB7TCouJLUD07fxfAESJKATgJ4DcCWldZgReiqESEffuCEaZ+jlOXhQXgiSfEAasoKwgkjpyZnzf277cx833MfDWIdTuJRsvOevNUm+9QanGQasy4UgXN7GwRjTZ4UAHe4dTitFTHplIFFeRtgCfstYlEB+PnOHVRx6ayCirIQ6DW+HDPHLNSG3/4YW3X1nGsdJyOjMhDHZtKDRB79/otZHx8nCcmJlp+3XaGSIS397yS9ThBiehZZh4PYn+A/lyjRNA/WyXaqEbeZlTT5tWGHhG0ZrgSAirIm0xQ5g9PgKs5JcJozXAlJFSQN5mgNGU/Ad5oSKMSMFpESwkJFeRtzFohjSrQW4zWDFdCQgV5E4iKpqx28xajha2UkFBB3gQaTf5ZD1ryNgIEXURLUWpEBXmH4JpTonA30JUEWURLUeogqKJZShVarSm78eTVYs6VJhJUEa0mQ0QfBvCnAOIAvsjMj4S8JaUBVCNvMqoJK1GDiOIA/hzARwDcCuBjRHRruLtSGkE18g6gWlan2s2VKtwF4AQznwQAIvoKgHtRpY9AKt7LvYlBa68DJOEJADtDlCtUnuwdLxTLQ9zbI4cKpcr5BWeNRFzmpyrFFC1k7ZukjxgryvW4JwUAyA1anTVmthLL2+nsHXbuYNln2ZIZiy/5jDnrFVNmn/Zjg3j5fACImY/LcWdecfka7l4XL5+5zMwVLdZUkHcABw9WjzNXFB92ADjtvD8D4GfcCUS0H8B+AEgnBvGuG34TxY3WkVtMi+iI5aykSp6V/qF89bpdaOdWec5ZKZfdPQwAiGcdKbdiDQBgI8jnbrNyK5YTaZh57bLda7Hkbdqea/5Z5Edlzyfv67HXmJV5mfP2usW0OdcR5Lmhiu1haYtI3sHXrOhcGjF7umivP7dTxlIzdixhUgwWt9iLpKfkeH7AXiMxL8/ze+z3kzkr38Xx//jpU5W7UtNKS1HBqrQLbo/dVGKVqoxKJFBB3kKCjOvW6BSlgbouZwHsct7vNGNKm6KmlTZFo1O6HK+ui1cSwKvrAtQSNfMMgJuJaC9EgN8P4Neqzs7lwWfOI34lXR5K9BktPWZ1QU4b80VPD1bC5y6WX3tHKe/Yw0s+v8DT1wAAfZO9dszY1UuT1jIUGxFTDS2z4cvr5AUxTyRn7BrpK+b5qrXRZwflOZ5ztpSQNcjZWrFXTBzJWTvIcZnX46yX75PvJTlnz40bs1DBMbekr/CyNdyx3Ab73SZnsSqqkTeZIDXnVmvbRPRhInqViE4Q0YOtvbqyKg3UdWHmAoBPAvgOgOMAvsbMx5qwS6VFqEbeZILUnKs5NZsRneKEqH0I4gx7hogeZ2bfyAZlnRw5IsL3zTcllf/Qodri0Bus68LMTwB4oqbJMQKle0DJpD1/SaJGilNXymOF998OAOiJWe1yaaeouunidrvexSk5d+/O8hCVRJst9tlQjYS5xrIIleuimsZ223ORNWp00ToHczs3AgCSP3xF9tY3ave5IPprrt/uM7tRXiec/425IfljdTXyQr/sM7fBhpkUzM1JbsCuV/RuAJybhELRc2zaBfN9ZNazY4lFMnu2Gn5u0FnIB9XII0i9mneTNPVyiBoz5wB4IWpKUDRS9lbruigOKshbSK2as+sUDdGp6ReitmPlJCLaT0QTRDQxNTXV9E11FI2UvdW6LoqDmlZayHrt4lF2ajLzYQCHAWn1FvJ22otGzCOe+WU9Zpk64XQK+Z/cjZiTwFPIiJml0DdWHuufkH0Xd9q4755/PC4vdmwtjy387C0AgMzr1izjOU2TF20cee4tojfkh6xJp5gSE8nAa9fKY9mxEdmno5b2/JOYVGY/8lMAgMET1jSRXJBf04HT1mRTetMkODlmob6LTpaOYW6bjI0ct17R+S2yv/5zdiyWFxNR+pr9zlLXxbk7t8Oaj3qvyFjmst18vBybb8XzwGmfZCsH1cgjQgTDCTVErdk0ah7Ztw+YnARKJXlugxovSnMITJATUZyIfkREfxXUmt1ELaVvW5xyXw5RI6IUJETt8ZbuoNNpI/MIMSRE0DzIPMCwj1hsWTgi4Ix5v9DMoBLkUSyVHyiZhx8l+yBmkLMWzHt5oPxw973yAb/His8ln40rboHLa7jfRXlPdj2/6/o/zLnuNYomTd9dz3xn1QhSI/8UJJRJaRKt1M41RK0FaNlbJSACsZET0U4A/xzAIQCfDmLNbiYqxa7qClFT1keblL1Vok1Qzs7/AuCzAAaqTXCL8OzWEKlV0TT7LmK9ceQthLIFJE9NgedsmmKivx8AkBrqL4/xYB8AIHZt3p68bbMcO3uhPOSFWJecGHTPrEIbN5SHUqfleDJjM0o5aSoinrYVr1LFLfLCcVTSqGR79r0he752g62A5RXeKvRaZ2Z2o7xOLFkzx+KIGCxcU03OSLilEeuALfR6cd9ObLkpwpXrt0aPYkrOcWPCU3Mxcy07lvEKadmvFksbKx2vLg2bVojoowAuMfOzq81zi/Bs2lRRhVFRuo9G4sgVxSEIjfzdAH6JiO4BkAYwSERfZuZfD2BtRelcVosjj5JWTiS1wbdvKQ8VMyl7zBB745y8cErW0qBRYdO2/kphVMYSPnXE+YzV3GlIskKz2wedCfLUc95q6eRldrqOVlP/vGTCJPODrkNR9hzPORq0p1Wn7ecp9FVmU3rr+GWFxrOuhm+269RQ8bJGc87HWVySPRcydt7isMk8HbJ7Ts41ObOTmX+fmXcy8xgksuFvVIgrSg00mGavKB4aR64oYaFp9kpABJrZyczfB/D9INdUlEjTiLPy0KHlpWiBaMaRl4rg2Tlg0fY3iy+Ky5JT1umHDcZmMGudnV6nnuI5azJJTJuys3PO5/ZiyDO23CwvLgIAktP2ul7MddEpB5HoMZUj4o75xHQXis/KubGcjddPLJrnrJ1fTBnzSN6OUbGya5DXci3mJFrGTIJo3FkvseSZb+yYX4s5z7madwpuJczHjeccs1V+eaz5SlQjV5T10qizUuPIlYDQWiuKsl6CcFa2QRx5YTCNq79wCxKLVissh9c5oXTDL4squXCndYoO/FjKzuY+dEd5bG67aPEbXl+0F/GcgkV7jdk9xqHp+Pnyxik4sOOny2PTmxMV80aOzgAALr5Lwg43HbUqNJkmFr2nbLeGoQXZe6nfOlE39Mg+3QbT87tFsx88au8IlvaIZ7Pngr0TSc6b8r2X7d1E4qL0Ml26wZbUjS/Ivja8Ykvwlnrk86TmrIO4/4TTB9UH1cgVZb2os1KJCCrIFWW9RNxZSUSPEtElInrJGRsmoieJ6HXzvDHMPSrBoKYVRVkv99wDfOEL/uPR4H8A+DMA/9MZexDA/2XmR0z7vgcBfG61RRJzeQz/v3NgJ8Oy1FfZl9OLIx86WhlH3jtpMzGTb5GimolrjmnFw4kjHzktZpHFmzfb9YxPtOdFe9fT1y/mDo5bvZQWxQO56UdiHnnzF6wTNWl6ZvZtsFmkXhy5W5jKiw93WdwsZpmlIWs+ypqszMwFW552dpeMpWZsemZiQTJfF7fYdXsvyTm5IafMrukHOjtmrzu3TRzEOFqxJQCqkSvK+nmiShmaauMthpn/DsD0iuF7ATxmXj8G4L6WbkppCqqRK8p6aU8b+RZm9tTjCwC2+E1yayOl4wPlTMkyXm0UxzlJPSbU8Lp1IsY2S9OHcl9NWGcjryx5C4B6rdbvZYguK0tb8gnD83p7JhxxZsIiKW8yPJ0oSTbTSgmrBXvHXY3cPceOsXl2zjXrFf2u4Yx5r0vONsvnWmUecbMvTtjP6rcXF9XIFWW9RNxGvhbM7FW89jtWro2UivX6TVEihApyRVkvbdQYwuEiEW0DAPN8KeT9KAGgphVFWS8t7JsZII8DeADAI+b52zWdxYziKz8uv40Z00Vsi61kmt8jr+MjtipUwThFk0vWtBJ//QwAgFLWnsDGVMNbR8pjdE1K0KaPnbH7iJvCVAlboIo3VFbPLqVl7fi0rJGcsXPSU3ITkr5uzUULZr3UvBMr7+MATc2Y3qILdp7XKzR93e1pasriOvO8dZK2GjB6r5SWzXfHshutnu2u44cKckVphAgn9BDRXwJ4H4BRIjoD4GGIAP8aEf0WgFMAfjW8HSpBoYJcUToUZv5YlUMfrGudVAL53ZuQcJpIlHpFI18csfbznkuS5VoYtA7L5CVxfLqhi/kx0dyTx07Zi5jQQVq0mnt+p2jnnsMSAAr9omn3XLAO1ew20bbdDMze18ViNH/bVgBA/1mr0aanZb3MabtG6rp8jtiSvVY8K3t2HawLi6I5D560oZPZEdlT5qzN8mWSUMOe6+568nphs70T6TtrasEU7feTviTO27xTd2bgrFOgxQe1kSuKorQ5KsgVRVHaHDWtKIqyOsyILRXK2ZIAEDPlZOMZx2FpCl+5seWeU9I9N74oZgKvTK0saEwrI7ZiQCxXNOdas0Lcy950uhDFcsaL6KqlxikaX5RjfjHenLLir5CW+UnHjFLs8Zyd1mZTNPHjxbR1Tnox5aWU44A134Ubb84xuZ5XMlfmmQ5BPc5YImbOtXsuplbXuVUjVxRFaXNUI1cUZVUol0fszfPgonXcUb84PpNuxqfRquPz18pD3CcOO7exRCy+XV4MOc0rjSZM81ZLj01L6VbeYkMS47Oi2RdP2uzZVN6UqHX6h3qvey5KrN/Cu+wapaRozhyzOQBLQ7L3hNN3c2GTuUtwbjCWRuVNas7eiWRNnRR2Eqfmt5m+m07DiJjx4y5usmNUkvDMhS1Wp+aYrL2w1ZnHdl9+qEauKFHjyBFgbEwE49hY7Y0qlK5FNXJFiRJe1yGvYYXXdQiIbLy6Ej4qyBUlSgTRdShoEglgdBi0ZB2WPC97JMfcMvOuPQCA9JSNBV/aJGaCvkGnF+cxyRDlW2+w1zAO0uxmO6/3jMR50wXbjad4+YqM3fFWO2aevX6eALC4XeK4+54TE0xhYNjOnxOTRXbQGiSWjLkjaUPLkTWnuJmd+X5TxnbEnls0FX2zTrekggkL55jr7PTWcNbrI7Oe3XtyXsYKfXZsaaSypK5LpEwrBw6EvQNFCZn2rKiohEykNPKDB1WYK13O7t1iTvEbD4tCEZi+Di65zk7ReDltG0wMvCiasxtqmLwgx0tv2H9EdINo7nTaqddlnJ2ZC04In1d/Zbut5xLbIU0mSkdfKY8ldu8wC1utte+qqNalzRLO6DWTAGy3+x6nNoqnLrtd770yt66z0yu92ztlz81ukLH0NTuW7zfhjHP25Ljn7Cw6TSRMbZfeKbeuSsmca8fSl1evtRIpjVxRup72rKiohEzogvzAAflH6v0z9V6rZq50Jfv2AYcPA3v2yB/Cnj3yXh2dyio0bFohol2QnoBbIEXqDzPzn9Z6/oEDVmgTAbz6HYSidD4Rq6jIPUkU925FfM6aTIqmaFZuo1Ps6ZyYM/KjtmRs8vxVAEB8987yWG6zHE9OOV3ovIzNtL0bKWy1PTXL182IyOq5aaw8lt0uvT1dx2L6+FkAwOJOuVbmvFM0y5hA+k7Pl8dSg2IC8rJJASCek7FlHYoKss/+M9ahm1iU76LvrI2BLyXkc6Rmis68olnXKZp1Xr7TeN6Opa9I1mqh1zFbnSlgNYKwkRcAfIaZnyOiAQDPEtGTzPxyAGsriqIoa9CwIDf9/86b17NEdBzADgB1C/KHH250N4qiBA1lc4ifPAfqtaGBXtnVtFMHxcvKTJy7aE8eFmdj8ez58lDClLSlDU5mp4FnbdeFhHldvGmHMyaacOmUbTbRM2s067iT/Wj6d/aek2MX77LafW5I5hV7bBxgrl+0+ZhTLXZp1Gj4bmbnZlNjJme15exGEy7o3E3M7TINKGat9Tq+JHta3OyGKcp34YYXemGRs2P2uqXk6qI6UBs5EY0BuAPA0z7H9hPRBBFNTE1NrTwMYH12cbWlK4ovSSL6HhG9TETHiOhTAEBEw0T0JBG9bp43rrWQEn2IAzJKE1E/gL8FcIiZv7na3PHxcZ6YmAjoumpXXy9E9Cwzjwe1XpA/167hyJH6WsXVOJ+IXgDwG67JE8B9AD4BYJqZHyGiBwFsZObPrbbF3q27+IYHPo3cBjcOT57SV6wmubDFC6WzY3O7TSidE/4XK5iEnFGnToupMBhftPO8deZ32LA+r9WaS/ry8j0BwMzNck7/GzL/od+xZQ6end8LADh61Wr6OzJS12XeaWf/tgGxs5dg9/Rzfa8BAB699HPlsVv6JIxy4uqe8th7R14HAJzN2juB6byEbI4PTpbH/v7qzcuuBQAvzUotmvs3W334b2ZuBQD81zu/4vs3G0gcORElAXwDwJG1hLiiKIZ60/Hrm59n5ueACpPnvZD2bwDwGIDvA1hVkCvRp2HTChERgC8BOM7Mf9z4ltZGQxaVjmC1dPwg5htWmDy3GL8WAFyARJv5nVM2hRYX5/2mKBEiCI383QA+DuBFInrejH2emZ8IYG1fNGRR6QjqTcdfR/q+MXl+A8DvMfMMOdmPzMxE5PvXw8yHARwGgPSOXZwf4GW1P2CsHXnHFFLYKCFy2bwjVjaIczIPa7KIe1F6Azakjo1ppdBjdcu8cQ6Wep3u9DlTMjZp9xLLVuqjvMGE8GXEKXlL0maRXkyJuWO6r688dlNGjk8X7Nhb0vL/ruiYVm5KzgAAbsxcLo/t7RGf36U+G3Z5Q4+slyRrPhpJyj/Em1O2pO/JtGSt3tprTSvXi+JUvjlpr3G613Eg+xBE1Mo/AM4nVRSlNupNx69zfhWT50Ui2sbM54loG4BLvicrbUWkaq2sBw1ZDB4i+kMAvwggB+DHEKfZtdXPUurm0KHlNm9g9XT8euf7mzwfB/AAgEfM87fX2masAKQvE2I5tyGCPKevWM240C/iJH3Z0dL7TNf7q24YnjwvJmwvM6+dmhv+5zkxS05YYWqmsoVar1eHxLlhyI7K2r1TMviDRVtp8fnZXQCAV69tLo/NFkRzn8nZBKe4ue3IO00d4uYiR69bR+nVgoQdHru2rTzWH5dEn4tZG2LpXcN1nr50TRybPTF7d/KiGXu6d6w89tyMdaT6EXqKfqOoXbwpPAngNmZ+G4DXAPx+yPvpTOpNx69vfj/E5PkBInrePO6BCPAPEdHrAO4275U2p+01ciV4mPm7ztsfAPiVsPbS8dSbjl/7/Dlmrmby/GDtF1TaARXkylr8JoCvVjtIRPsB7AeA3WGWWu0k6o0tbzIck2YIXlMFwDrFYnnHjJIpmWc7VjJjeccs43WHL/a7ZWTNekuOs9NkW+YHK/2xbgvL+FJludlinzgZc0Mi4m5IWVfA5Yw4Ja/nrRllR1oshwMJm726u0eaWBQdw8WNSXFs7sxYS+ON6Smzrs0UHTN2Iddk4plgPOcoALzRN7rsWgBwxThh3T2/0WtL+frR9qYVZX3cfffdAPBWInppxeNebw4RPQSppVO1aSQzH2bmcWYe37Rp9V82pQa8WPFTpyQcy4sV176dyiqoRt6lPPXUUyCiY9UyO4noEwA+CuCDHFT6r7I2UWz1xuKgjPdUWmritiAiEgtyPGGLACK2YLrTLzmOUuPQdLVvL5wx7szzwhTjWWcsV7mHxGKlszM+H1u2vwuFofKxSznRyK9mbW2UdFw054WCdcCeT0uYYpHtPi+kxHl5KWtDDQcT4r2ddta7nJfjUzk7bzonxy/22L1cyYr2fTFfOebu2V3HDxXkSgVE9GEAnwXwz5h5Ya35SoBoqzdlHahpRfHjzwAMAHjSRDv8Rdgb6hpWiyFXlCqoRq5UwMw3hb2HrqX+WPHmQwAnAI5XFs1ymzmUTLalG+PNKTOWcHpXmp6Vy9YzKiUnnHMT3jXsPCbvXGcr8cpys6Wk6fdpDvXFbCMIzwGZjNmsy76E2GBKTqBPxpzjmlbSxi7UG7cB72kT/J5y1vPGepzAeO+cTMzao1LGpJNxbFTeOn3OvF5n/36oRq4oUUJbvSnrQDVyRYkaEWv1JpmdDCpWOhp7rjld55OiF2Yu2rFi2nSTn3VO8rR5txGEGUs4NyJe5/h8n9U3U2Ydcjqf9V4uLVsDALIbZe3UjAz+cN5mdr46K3XCTl+1JWZnTaOIXMGKxAJX6rlJo82/ctVmhS4WxUH6xtXh8tjGHsn8PL9gMzuvZyXcMebESb5xbaRibPK6rPPD/hvLY14GaDVUI1cURWlzVJAriqK0OZE2rbjlahVFCYdiCpjZCxQGnExM87KQtrrg4pg487xsSgAo7JFg8KWrtoxtwsSWF7ZaB55XxhZupmifmCyWNls7SjFj+m2mrCkiP+B5Su32lm6W2G4yvUU/tsF22/n71C0AgC3pneWxW/qktOzVvC1je1f/SbmWY2L5mfQ5AEB2h403v6lHSsxu7rExAh8YOg4AODtoO+ldN5md45mT5bENydsAAO8ZeLU89o8Z2d/9Q8+Ux7YlJZP0KfgTaY384MGwd6AoihJ9Iq2RK4oSAQjg5PJmDp5GzklnXsKE/DlSJWZCDAtJ1ykqr2NJq+F7ucMlJ67Qmwefc3nZmNmmc8MQM+GO3l6GnDDAobg0ePBCDgEbapiP280PxORuwtXIB2Kyv36vFi+AAZOC2u+EEHrnDsRsPZdSXNbZELceXe+cwZhdbyjunWs/kBuy6EfkNHJt46YoilIfkdPItY2boihKfUROkCuKEjFKQGKOUHKyLr06tm7CYfyyODRT1+y8hX6Jz07MVd78F+Yc8WMyKuPOvMScjOVGnIJbJpY9MWvHknPeGnY5z7nqxZu/nrdOx4umGNW5RVuUysvozJYc00rcOkM9NsflYicWbM9qry/nmSUbl362V673RtbGm3uFr4YTc+Wxc1nZw4uJXc6YrHOqYItwvbZkuw/5ETnTiou2cVOU9UNEaSL6IREdJaJjRHTQjO8loqeJ6AQRfQVjNhwAABOkSURBVJWIUmutpUSbSGvk9drFNVxRUZaRBfABZp4zjZj/gYj+D4BPA/gTZv6KKYj2WwC+UG2ReBYYnCwhO2T1PiqJ+tt72dYXuX6DOAL7zjmZnSkZS08765mys4tbHE+p11jCydjsPyvOvnjOiqnUda+ei52XuSR7cBtLFEyY4sbX5NjjV+8sH3vpqmi356Zt1uWwCRNcyNmFL26U426PzcsbpHnEP12wPTRP9kt25uQVm9k5l5c7kcuLttnEkska9crZAsALFyVjc2rUzjsxLc0m3Dotz15ZvWhapDXyetFwRUWxsODdxyfNgwF8AMDXzfhjAO4LYXtKgHSUIFcUZTlEFCei5wFcgjTV/jGAa8zs6b5nAOzwOW8/EU0Q0UQhO9+6DSvrItKmlVo4cGC5Ju6FLT78sJpZFIWZiwBuJ6INAL4F4CdqPO8wgMMAkNm8i/N9hEKfM6HkdQOyumDRtLt0kiNRzJjCV06XH6/0baGvsixuzOkAlO8z86wlAlQwZWwdyVUuqrXMtLK84Nb2Httj80KvOB1nMj3lseFeie1Oxe3Ylt4Z+ahOHPmO1FWZn1l05kklr+mM3ejWXrdKmLBQEFfEtvT18tipzMaK+Zd7xcyyw9nzqV5rtvEjEEFuOsr8KYA4gC8y8yNBrFsrXoiihisqij/MfI2IvgfgXQA2EFHCaOU7AZwNd3dKozQsyIkoDuDPAXwIcpv2DBE9zswvN7p2LRw8qJp3x+J1kz91CojHgWJR6nOH3FW+XSCiTQDyRoj3Qv5G/xOA7wH4FQBfAfAAgG+vtg7HgaVhWtbN3suipJLVoLMjxumYt9mZ+WFx2LEjaoqm92du2DpKPW2a8lb7jefkdW7IyQAlGSumXW2elu1JriuWo6WNct23Z06Vj80VResulOw+b+yXzvbX8larvqNf2uvlHfX/9rSs88KgDRe8oVfO9crZAsBb++R/43DShiR61317n23bN5UT7fttfafLY946b++1e75ubne+BX+CsJHfBeAEM59k5hzkl+PeNc4JnAMHNFyxo3C7yQMixIHqXeWPHAHGxoBYTJ616zwAbAPwPSJ6AcAzAJ5k5r8C8DkAnyaiEwBGAHwpxD0qARCEIN8B4LTzfk3nydTUVEMX9Evj94tYUU09YtQjbP26yXt4XeXddT2hz1xd2HcZzPwCM9/BzG9j5tuY+T+Y8ZPMfBcz38TM/5KZVy/koUSeljk7XefJ+Ph4Q5ZsvzR+v3osanaJEJ6w9YSzJ2wBfzPJWl3j3eN+Qt8T9mqCaRgqAqnrsJoTbC/MxKL9U05fEFNF75QTR94rIiY5V9ldKHndp0OQM6/nqgzmhuyYF1qdWLRj6Wmv5ZBdLjUla8ezMvjyktUtvQzMs7M2s3OpKPvMFq1ITJjema6z0yuQ9eqMzdj04sxPz1ozysn0JrmGk+15dUnMNj1OsPzpOXF2/ihu48TPzMs5Lw/YzNJjM83P7DwLYJfzvmnOE08o+2nk3nOtBbZUwLeY1YStH2t1jXePVxP6a/0zUJQOIQiN/BkANxPRXogAvx/ArwWwbgWehl1NI2e24YieqaVaOOJa2rpmiQZMvcLWr5u8x8qu8rt3W1u6y1r/DJSaKCWBpc1A3mks4YUBxpw+nl5jCU8LB4DCHinPWrhiw/o8rTu/yWYuetp0ftjR+uPi9Ctk3OvCXMMpY5sQfdTN7MzvEmtR/rpc9z2Z1+x8n16cP9EnzSEu5QfKYz878Lqs4Tg732mcnW+MbiqP7U5dkT2x3fud/TLPDSGc7ZeStndkJstj8wXZ3/tNIwoA+NvYWyr2XDRaf9OcnSaE6ZMAvgPgOICvMfOxRtddLwcOiED3whC91/UKZc0SDZhqQrXauNtNHpCoFcC/q/yhQyLcXVYKe0XpYALJ7GTmJ5j5Fma+kZkD/etZqz65F6lSS8SK1joPkfUI2337gMlJ+U9cKMjz5GSl3dsV+kT+wl5ROpjIZ3auVZ/ctZuvZKVwr2UtzRJtEp5QfeghMafs3h1cPLgXbx70ugoAIL4EbHi1hOwGx2Fp/nZ6r9hY8HhWTCGZC9YUMjMv8c891+wfm1f6dnHeFl30YsCdJjvoOy+Dczutvpma8daxe8lMFZbtSRCTxdAbsr8vTb23fOT4VSlBe8EpmvVSvzgTs3krEl8bFodmyTGZvLBR3IF/feony2NbByQrc/Kyzb48v1UcqVcWrfKykJXP++Lw9vLYyxe3yrmj9tzJ6WFzXfu5f3TFc9b6d+2MvCBvhEaqJ2qWaBPYty94AVtvNIyidCBtJciDSPjxhLUmD3UIGnrYdEopYGYshtxQpWaTG7IhhNkRr8Ss1SQXt4tGXO50D6t950YKzphpLDFv55VS8np+t9X6EzPGsekkheYzlWJswVyXirK/+0aeLR/bmb4BAPCjXhtstz0j9U/c7Myf6pfgu6Kj/XsOyOlttqDMT/adBwD8oGdveewDw68AAN7MjpTHvMYS7xz8cXksk5Dbk9sHzpTHjqVFY//l4YnymNdf9AcVn1Roq+qH9WjY1eZ6ppO11lJB3yZo6KGitJcgr4dGo07UJt4m1BsNoygdSFuZVtaLOjE7GL94cw09DBYG4jmbVQnA9ux0xmLZ5c8AEFsSXTHulLG18x09slzG1jluXlPeyewseNe1Y/FcZWant7a3vysF24FnpiDx3LN5G9vude3JFa2p6HLenuNxpSQmlev5tD3X1Pedc9a7XpT1rjk1eL1zpovWLOOdM+3UCJ41seVXivb61/K9FXtx6SiNvFp44fe/H+aulKbhRassLKweZ64oHU5HaeS1RJ1oNEqHsDJapVi0mrgK8UApJYGF7SUUhwo+x6xzsLhFVOhSjw0rjO2S7kJLs1ZbhckKTY/a5gwl4+zML9r1SilZh7YslceyvZV9oos95p+4q5HvlN+LRaMZ/6JTOnZrQhybgwm77t6eSwCA2ZLVfN+Rnqy41jvTcq0LIzbn8a094hQdSti7wl/ol+OTaRtWOFWQcMf3ZKyzM21SVX8mc6I89lxqDADw0cyUM+8oAOCLFTsSOkojV7qIemu3KEoHEwlB3gw7dbWoE41G6RA0WkVRykTCtNKMcrPV1lPnZoeghbJaBpWAxDyhlLDiwitQ5ZanLQya0rFOidmsMZXE5q0T0evL6ZpJ2Os05DhAEwsytjRrzS3xWbNOzNpREvOVjtQFc93eWTl2qmCvfzInGZtnljaWx5ImMN3rxANYE4wbR76zIGaUMzlrMukz3t3TS052Zo+8Pp23ceRXjUNzMm9L276ZlXnbkyPO2IjZ8+vlsVM5W6TLj0ho5IpSN1ooS1HKEIfg+RsfH+ePfnTCN9ZbQwJbBxE9y8zjQa03Pj7OExMTa08MCq2xUpUgf7ZENAXgFIBRAJeDWDNE2v0z7GHmCvU8NEHu/sFrJEk4tL0gV6oS9M/WrDkR9JqtphM+gx9qWlFWhYg+Q0RMRKPrXkQbIytKU4mEs1MjSaIJEe0C8PMA1h8KotUJFaXpREIjV5t4ZPkTAJ/FikrPdaHx3p3E4bA3EACd8BkqiIQgV6IHEd0L4CwzH21oIY337hiYue2FYCd8Bj8iYVpRQuMWInrJZ/whAJ+HmFVWhYj2A9gPALv9Yrg13ltRmk5kNHI1r4TCa8x828oHgJMA9gI4SkSTAHYCeI6Itq5cgJkPM/M4M49v2uSTtLBavHctTlB1lCrKmkRGkLsx5SrUw4WZX2Tmzcw8xsxjAM4AuJOZL9S9WLXGyIA4PU+dkthTzwnqCmrPUbraHKXpENGHiehVIjpBRA+GvZ9aIKJdRPQ9InqZiI4R0afM+DARPUlEr5vnjWut1Q5EIo4cWB5LrnHlraHWWGOjlY8z86qJFHXFkY+N+Ztc9uwBJidrn6P4ElQcORHFAbwG4EOQf+jPAPgYM7/c6NrNhIi2AdjGzM8R0QCAZwHcB+ATAKaZ+RHzT2kjM38uxK0GQqgaebX64Uq0MJp5sNlwtThB1VEaBe4CcIKZTzJzDsBXANwb8p7WhJnPM/Nz5vUsgOMAdkD2/piZ9hhEuLc9oQty5up1w4nUzNKx1NKiTdu4RYEdAE4778+YsbaBiMYA3AHgaQBbmPm8OXQBwJaQthUoDQlyIvpDInqFiF4gom8R0Ya1z6qOK9S91yrIO5Rail5pYSylQYioH8A3APweM8+4x1jsyh1hxG1UI38SwG3M/DaIHe3317uQZnd2GdWcoG62Zy1zlGZzFsAu5/1OMxZ5iCgJEeJHmPmbZviisZ97dvRLYe0vSBoS5Mz8XWb2+j/9APJDXheu5q1CvUvYt0+clqWSPPsJ6FrmrERDFoPkGQA3E9FeIkoBuB/A4yHvaU2IiAB8CcBxZv5j59DjAB4wrx8A8O1W760ZBJkQ9JsAvlrt4JqJIw5qTlHWjdZ2CRRmLhDRJwF8B0AcwKPMfGyN06LAuwF8HMCLRPS8Gfs8gEcAfI2IfgtSmvdXQ9pfoKwZfkhETwGoSAQB8BAzf9vMeQjAOIB/wTXEM2q502jQkWVsNWQRQHPK2CrRZU2NnJnvXu04EX0CwEcBfLAWIa4oTUVDFpUupNGolQ9DquP9EjMvrDU/aNQEo1SgIYtKF9Jo1MqfARgA8CQRPU9EfxHAnmrGr1Wc0uVoyKLShTQatXITM+9i5tvN498EtbFaUa1cWYaGLCpdSGSKZtXKyrT+gwc1A1RZwXpCFhWljWm7euQHDlih7QlzdbEqitLNtJ1GDlit3EPrsnQxWtNcUdpXkDPbDFCty9Kl1FKvXGuaK11AWwpyDxXcXU4tjZ21+bPSBbS1IAe0LktXozXNFQVABwhy1cq7mFqSf4aH/edUG1eUNqTtBbnSxWjyj6IAUEGutDO1JP9MT/ufW21cUdoQFeRK+DQSHrhW8o/WXlG6ABXkSrg0OzxQzS9KF6CCXAmXZocHau0VpQtouxR9pcNoRXjgvn0quJWORjVyJVzUhq0oDaOCXAmXem3YWjdFUSpQQa40D1fojo7KY6UArseGrXVTFMWXNZsvN4NINOlVmtt8eWU3+5VkMvU7HbWxcs1o8+XuQjVypTn4RaO4rCcyReumKIovKsiV5lCLcK1XAKtjVFF8UUGuNIdahGu9AliTexTFFxXkSnPwE7ou6xHAmtyjKL6oIFeaw0qhOzIij0YFsFtb5dAhsbNrKKLS5Whmp9I8mplRuTIqxgtF9K6rKF2EauRKe6It3BSlTCCCnIg+Q0RMRKNBrKd0AM3OwNRQREUp07AgJ6JdAH4egP4FKUIrMjA1FFFRygShkf8JgM8CaDhFVPtvdgitMHtUi4qZm1Onp9J1NCTIieheAGeZ+WgNc/cT0QQRTUxNTfnOOXiwkd0okaFVpWkPH5ZIGJcrV7T+itJ1rCnIiegpInrJ53EvgM8D+Pe1XIiZDzPzODOPb9q0qdF9K1GmVWaPffuA/v7KcXV6Kl3GmoKcme9m5ttWPgCcBLAXwFEimgSwE8BzRLS1ng0cOCChxUTy3nutZpZwIaLfJaJXiOgYEf1BXSe3MgNTnZ6Ksv44cmZ+EcBm770R5uPMfLmedQ4csEKbSHxjSrgQ0fsB3Avg7cycJaLNa52zDC+O+6GHRKDu3i1CvBnx3bt3+1dEVKen0kVoHLnix28DeISZswDAzJfqXmGt7vZBofVXFCU4Qc7MY/Vq4yt5+OGgdqM0yC0A3kNETxPR3xLRT1ebWIsTu6lo/RVFiVaKvtrFW8fdd98NAG8lopdWHHoI8nsxDOCdAH4awNeI6Ab26ULCzIcBHAaksURTN10Nba6sdDmREuRK63jqqadARMf8usgQ0W8D+KYR3D8kohKAUQAhqNyKoqyF2sgVP/43gPcDABHdAiAFoCGzmaIozUM1csWPRwE8aswuOQAP+JlVFEWJBirIlQqYOQfg18Peh6IotaGmFUVRlDaHwrhjJqIpAD5ZHHUziva33Yb5GfYwc2D1EgL4uXbCz7NWmv1ZA/3ZKtEmFEEeFEQ04Rd10U50wmcIim76LrrpsyrNR00riqIobY4KckVRlDan3QX54bA3EACd8BmCopu+i276rEqTaWsbuaIoitL+GrmiKErXo4JcURSlzWlbQU5EHyaiV4noBBE9GPZ+1gMRTRLRi0T0PBFNhL2fKEBEB4jorPlOnieie8LeU9B0wu+uEi3a0kZORHEArwH4EIAzAJ4B8DFmfjnUjdXJersqdTJEdADAHDP/Udh7aQad8rurRIt21cjvAnCCmU+auiBfgbQmU5Soo7+7SuC0qyDfAeC08/6MGWs3GMB3iehZItof9mYixCeJ6AUiepSINoa9mYDplN9dJUK0qyDvFH6Ome8E8BEAv0NE7w17Q62AiJ4iopd8HvcC+AKAGwHcDuA8gP8c6mYVpQ1o1zK2ZwHsct7vNGNtBTOfNc+XiOhbkNvuvwt3V82Hme+uZR4R/TcAf9Xk7bSajvjdVaJFu2rkzwC4mYj2ElEKwP0AHg95T3VBRH1ENOC9BvDzAFb2z+w6iGib8/aX0XnfSdv/7irRoy01cmYuENEnAXwHQBzAo8x8LORt1csWAN8iIkB+Dv+Lmf863C1Fgj8gotsh/oNJAP863O0ES4f87ioRoy3DDxVFURRLu5pWFEVRFIMKckVRlDZHBbmiKEqbo4JcURSlzVFBriiK0uaoIFcURWlzVJAriqK0Of8f9MFwsuzk7bAAAAAASUVORK5CYII=\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",
    "by transporting 50% of the mass\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.485157824314826\n",
      "Entropic partial Wasserstein distance (m = 0.5): 0.5048991597745197\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEtCAYAAADHtl7HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfa0lEQVR4nO3de7xdZX3n8e8XEi7hFgKYhptRikWKGiAiVASK1eE64AyDWMuEGTHQKVOZqgV1ZjjMiI0tl9YXrZWbpGq5CHIRoYVRkGIlkNAAgaBcDCUxF65ylRL4zR/rObBzstc++9l77bX3Oefzfr32K/s8z1rr96x9zv7lt9dez1qOCAEAAKB9G/R7AAAAAGMNBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyUUABAABkooCaAGx/0fZFbS57qe0v93pMY43tT9q+ud/jAMY62zvbftH2hv0ey0jkyu5NpFxJATUAbC+z/UpKKqvTG3PzDrd1kO3ljW0R8ZWIOLHLMU5K4/tAQ9snbUeTtoe6iVUH20O2v93u8hHxnYj4aC/HBFRlRE4Zfpzf5rq32e4qX7QSEf8aEZtHxOu9itEOcmV7yJXlKKAGx5ERsbmkvSTNlvQ/czdge1Llo0oiYq2kn0o6oKH5AEkPNWm7vVfjaFcvXwtgjDgyFSrDj1Oq2Oh4eG+RK98yHn6f/UIBNWAiYoWkmyTtIUm2/4vtpbZfsP2Y7ZOGlx3+BGX7NNurJF2W1t2+4VPn9iM/Qdj+ru1Vtn9l+3bbv93m8G7XugngQ5K+2qTt9hRnH9s/tf2c7ZW2z7e9Ueqz7fNsr7H9vO37bQ/v82G2H0z7vML25xrGfoTtxWmb/2z7vQ19y9JrcZ+kl9InwdPSNl6w/TPbH7Z9iKQvSvp4eo3uTetvZfviNNYVtr88/DWD7RNs39EQK2yfbPvhNJa/tu02X0egb4b/lm2fbftZ27+wfWjqO0vFe/j8xqNW6e/9j2w/LOnh1PY7tu9OeeRu27/TEOM2239m+670/r7O9rTUNzNtb1L6eZrtb9r+ZRrPtS3G/ZOUR35l+yHbH27oJ1eSK+sVETz6/JC0TNLvpec7SXpA0v9NPx8uaRdJlnSgpJcl7ZX6DpK0VsUbc2NJm6a25SO2PyTp2w0//1dJW6R1/lLS4oa+SyV9uWScB0p6RkXhva2kxyVNkbS6oS0k7ZyW31vSvpImSZopaamkU1Pfv5O0SNLUtG/vljQj9a2U9KH0fOuG/d1T0hpJH5C0oaQ56bXbuOF1XJxew00l/ZakJyRtn/pnStql2WuS2q6R9A1Jm0l6m6S7JJ2U+k6QdEfDsiHphjT+nSU9KemQfv8t8eARsW5OadJ3gqTXJH06vY/+UNIvJTn13ybpxBHrhKRbJE1L761pkp6VdHx6f38i/bxNwzZWqPgguJmkq4ffb+l9GJImpZ9/IOmK9F6fLOnAFuNeK+l/pOU+LulXkqalfnIlubLWB0egBse1tp+TdIekH0v6iiRFxA8i4tEo/FjSzSo+uQx7Q9IZEfFqRLzSTqCIuCQiXoiIV1W8Od5ne6s2Vl2gIgm8J43hjoh4WdIvGtqWRcS/pjiLIuLOiFgbEctUvOEOTNt6TUVi2k1F4l4aESsb+na3vWVEPBsR96T2uZK+ERELIuL1iJgv6VUViWfY1yLiifRavK4i8e1ue3JELIuIR5vtmO3pkg5TkbReiog1ks6TdFyL12NeRDyX9vdWSbNGfwmB2lybPvEPPz7d0Pd4RFwYxXlI8yXNkDR9lO39WUQ8k95bh0t6OCK+ld7fl6n4iurIhuW/FRFLIuIlSf9L0rEeceK47RmSDpV0cnqvv5byXJk1kv4yLXeFpJ+lsZAryZW1o4AaHEdHxNSIeHtE/LfhN7jtQ23fafuZVGAdpuLTy7AnI+LX7QaxvaHtebYftf28ik8iGrHNplKcu1Qchj5A0j+lrjsa2t78Tt/2u2zfkA6BP6+iKNw2betHks6X9NeS1ti+wPaWadX/mPbzcds/tr1fan+7pM82/qeg4hPU9g3DfKJhvI9IOlVF4ltj+3Lbjcs2eruKT7UrG7b9DRWfrsqsanj+sqSOTvwHemQ4pww/Lmzoe/NvN/3HLo3+9/tEw/PtVRxVafS4pB1Kln9cxftrZJ7ZSdIzEfHsKLGHrYgoDms0bHd7iVxJrqwfBdQAs72xikPfZ0uaHhFTJd2o4jDusBix2sifR/p9SUdJ+j1JW6k4VKsR22xl+Lv9D+mtpPBPDW2NJ0V+XcWn0l0jYksV36W/GScivhYRe0vaXdK7JH0+td8dEUepeENeK+nKtMoTks4a8Z/ClPTp983NNg42Iv4+IvZX8aYPFYfw11subftVSds2bHvLiGj3nAdgvCjLIY3tv1Txnmq0s4qv7YbtNKLvNUlPjVjnCUnTbE9tc2w7jDh/ZmdJvyRXkiv7gQJqsG2k4rDqk5LWujjRc7TpoaslbdPiMPMWKv74n1ZxiPkrmWO6XdLvqkiOD6a2n6g4n2CW1k0KW0h6XtKLtndTca6FJMn2+21/wPZkSS9J+rWkN2xv5GJ671YR8Vpa/4202oWSTk7r2fZmtg+3vUWzgdr+LdsHp+T6a0mvNGxrtaSZtjeQpHRI/GZJ59je0vYGtnexfWCzbQPj2GpJ7xxlmRslvcv276cTkD+u4j/3GxqW+QPbu9ueIun/SLoqRly6IL3vbpL0N7a3tj3ZduOJ1iO9TdIfp+X+k4rzgW4UuZJc2QcUUAMsIl6Q9McqPlU8q+IT0fWjrPOQihkmj6XDqyMPw/6disPeK1S8qe/MHNY/q/g0tmD4UHpEPKUica2JiIcblv1cGvMLKt7QVzT0bZnank3jeVrSX6S+4yUtS4eyT5b0yRRnoYoTX89P6z2i4oTFMhtLmqfiU+8qFcn3C6nvu+nfp20Pnzfwn1Uk4gfT9q9ScW4IMBZ93+teB+qaNtf7K0nHuJgR97VmC0TE05KOkPRZFe/dP5V0RMoFw76l4kTrVZI2UZHLmjlexdGph1Sc43Rqi7EtkLSrivf0WZKOiYinyZXkyn4YnnUBAEAlbN+mYuZWW1f1bnObJ6iYHbh/VdsEusERKAAAgEwUUAAAAJn4Cg8AACATR6AAAAAydVVA2T7ExT1zHrF9elWDAoA6kMMAdKrjr/BcXJL/55I+Imm5pLslfSIiHixfZ0oUt8NBoxla2bR95cSbFYpxaeVTEbFdv0cxUm4OszeL4nZjI73RpK1Tm7boa+vuI295T4v8cX/znNOJXfZ+qWn7o4ta3fFkbWaUDVv05f4fVuXvC+Nfef6a1MVW95H0SEQ8Jkm2L1dx1dbSAqoonuZ2EXJ8mqszm7afyWuFceHMkbf8GBSZOWxrSf+9SXtmYdPSu1v0Lc3b1Pe/VN4386y8bbVwzsKfNm0/2oe3WOuZzChbtuh7LXNbVf6+MP6V569uvsLbQeve62i51r0PEgAMMnIYgI51cwSqLbbn6s3DTu3cxBoABsO6+YvTDwC8pZsjUCu07s0id9S6N5KUJEXEBRExOyJmF7cTAoCBMGoOWzd/bVbr4AAMtm4KqLsl7Wr7HbY3knScRrn3EAAMEHIYgI51/BVeRKy1fYqkf1QxReKSiHigspFNIGfqjH4PoakzSk9uH8zxAjnyc9gb6v0JyI+Ud90x1Lx9/5L2Ck8Ub+Vo71fS80KFUZ4v71pecvWJpzZp3j5rqOvRAFKX50BFxI2SbqxoLABQK3IYgE5xJXIAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADI1PMrkTeaoZVN7/vGtPjBxO8FeMv0vaU5C71e+5+7kxuyl90xZr1rEb+l7HIFHZle0r66uhDLP1vet+O8zI21uMnyjteVdHBXHvQWR6AAAAAyUUABAABkooACAADIRAEFAACQiQIKAAAgU62z8FZqhs7U3Eq2xY1uAdRp9aJOZ9yt7+clefBdteWvCmfblcmeaddKq5s4L81sb2VOSfv8DraF8Y4jUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCp1ll4Vapyth0z+gDUqXy2Xav7tzW/T96j8Y2m7bv4pLxBTRCTn/qT0r7Xtj23xpFgrOMIFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGTqahae7WWSXpD0uqS1ETG7ikEBQB3IYQA6VcVlDH43Ip6qYDt9w+UK8pRd9kHitcSYNEA5rNVNc5s7WLeW9Fze3VDGqTu32be0b28dW+NIMNbxFR4AAECmbguokHSz7UW251YxIACoETkMQEe6/Qpv/4hYYfttkm6x/VBE3N64QEpKKTFt1WU4AKhUyxxG/gJQpqsjUBGxIv27RtI1kvZpsswFETG7ODlzSjfhAKBSo+Uw8heAMh0XULY3s73F8HNJH5W0pKqBAUAvkcMAdKObr/CmS7rG9vB2/j4i/qGTDTGra2zhd4JxopIcdm38tLTv6OP/sXnHt4dK1ngmN7wed9lsu8nlK+37pebtdw6VrLBpixHkzxyUppW05+9/rr3NTDtUo+MCKiIek/S+CscCALUhhwHoBpcxAAAAyEQBBQAAkIkCCgAAIBMFFAAAQKYq7oXXNWZ19VfZLEh+L0CjDdRsNtrR3q/FOkNZEfaIw0v7lvj9mTFapPfS2XZlWsy0O6ZkW1dd3WJ792fGb6VkRmHpmM+tMDYmMo5AAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwDcRkD9B43bAa6tY2kE5q0X9pinbVZEZb4By16W/U1cdxp5X2Xf62ko4Ob+V41VNJxZIuV8i5jcFPcVtp3qMtudNziZspABTgCBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkcEfUF8/Yhza0tHjAejb2bP5+5KCJm93sU3SJ/oQobrPp80/Y3fuMvahrBsSXtV9YS/eb4YdP2j/rDtcTPV56/OAIFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmUa9F57tSyQdIWlNROyR2qZJukLSTEnLJB0bEc/2bpjt455vGO/4O85TXQ6zmt9f7bXKxvrVeLq07zSX3dvuopL2VveCq27MOnGoeftFF7ZYaUV18fUfStpfKWm/qcLY+eqbbVemntl2ZQZ3tl2+do5AXSrpkBFtp0v6YUTsKumH6WcAGESXihwGoGKjFlARcbvWv0X3UZLmp+fzJR1d8bgAoBLkMAC90Ok5UNMjYmV6vkrS9IrGAwB1IIcB6Mqo50CNJiLCdunlzG3P1ZuX792q23AAUKlWOYz8BaBMp0egVtueIUnp3zVlC0bEBRExu7gU+pQOwwFApdrKYevmr81qHSCAwdZpAXW9pDnp+RxJ11UzHACoBTkMQFfauYzBZZIOkrSt7eWSzpA0T9KVtj8l6XGV352wdq2meI+9m7AC6NZYymGneZvSvjnx/abt8122RoWXKmjloqGSjr1arFThZQzmvbd5+7KS5f+2v5cxwPgxagEVEZ8o6Ro/F3MAMG6RwwD0AlciBwAAyEQBBQAAkIkCCgAAIBMFFAAAQKauL6TZa1XOnGO2HYDOhZrPbPtg6RpL49NN29/tE7Kjz/fq7HVK3TbUvP2gkvaO3NOib4eS9g5m550+lL8OUAGOQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMA38ZAy49MHFx82eMDT8p7enkcgW1KLtcwY4l7atabGttyTotrelgnbHl355rnr82mjrW8lf5ZTpa/e1PBByBAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwDPwtvvCibUSYxq6wMrwvGvskl7c1uStx/ez1xR9P2j+uK0nVO8zYdRDqypP17HWxrMI292XZlJvZMu1Y4AgUAAJCJAgoAACATBRQAAEAmCigAAIBMoxZQti+xvcb2koa2IdsrbC9Oj8N6O0wA6Aw5DEAvtDML71JJ50v6uxHt50XE2ZWPaIzIvU8bM8qAvrlUVeSwSdtLU4fWbz+xxTrzymbblc1qe6Tt4YzuN8u7zv6Dps33eKh5uzqZaXdQi77c2XZ/2qLvvJL2sv/eXsmMDTQ36hGoiLhd0jM1jAUAKkcOA9AL3ZwDdYrt+9Lh8a0rGxEA1IMcBqBjnRZQX5e0i6RZklZKOqdsQdtzbS+0vVB6ucNwAFCptnLYOvnrjSfrHB+AAddRARURqyPi9Yh4Q9KFkvZpsewFETE7ImZLUzodJwBUpt0ctk7+2mC7egcJYKB1VEDZntHw48ckLSlbFgAGDTkMQLdGnYVn+zIV0ym2tb1c0hmSDrI9S1JIWibppB6OEQA6Rg4D0AuOiPqCefuQ5tYWbyLKvbwC0HtnLiq+wh/bxlf+GsybHJO/MHjK8xdXIgcAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBM7dxMeMIqmxEiDe6skEEdFzD2bSNpzvrNk95RvsraW0o67i9pfz5zTC38w1B53yFls6/Lc16+97Toa77/pflr3lD5pk4/q+0RFfo70xDjB0egAAAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBP3wqsY93ICRuJeeGPFEfHbTdtvWHFk+Uo7zusgUj/vxfeBFn0Lqguz+VDz9hdL2gfWli36Kpw1OrC4Fx4AAEBlKKAAAAAyUUABAABkooACAADIRAEFAACQiQIKAAAgE5cxGGBj8WbGwPrGyWUMNp0dmrlw/Y6HhjrYWh3T+E9s0XdRhXGam/zUn5T2vbbtuZlb26u8a+a/b97+GyXL3zmUGRsTG5cxAAAAqAwFFAAAQCYKKAAAgEwUUAAAAJlGLaBs72T7VtsP2n7A9mdS+zTbt9h+OP27de+HCwDtI38B6JVRZ+HZniFpRkTcY3sLSYskHS3pBEnPRMQ826dL2joiTmu9LWbhARNP/2bhkb8AdKeLWXgRsTIi7knPX5C0VNIOko6SND8tNl9FUgKAgUH+AtArWedA2Z4paU9JCyRNj4iVqWuVpOmVjgwAKkT+AlCltgso25tLulrSqRHxfGNfFN8DNv0u0PZc2wttL5Re7mqwANAJ8heAqrVVQNmerCL5fCcivpeaV6fzC4bPM1jTbN2IuCAiZhffIU6pYswA0DbyF4BeaGcWniVdLGlpRDRef/96SXPS8zmSrqt+eADQOfIXgF6Z1MYyH5R0vKT7bS9ObV+UNE/SlbY/JelxScf2ZogA0DHyF4CeGLWAiog7JLmk+8PVDgcAqkP+AtArXIkcAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADINGoBZXsn27faftD2A7Y/k9qHbK+wvTg9Duv9cAGgfeQvAL0yqY1l1kr6bETcY3sLSYts35L6zouIs3s3PADoCvkLQE+MWkBFxEpJK9PzF2wvlbRDrwcGAN0ifwHolaxzoGzPlLSnpAWp6RTb99m+xPbWFY8NACpD/gJQpbYLKNubS7pa0qkR8bykr0vaRdIsFZ/wzilZb67thbYXSi9XMGQAyEP+AlC1tgoo25NVJJ/vRMT3JCkiVkfE6xHxhqQLJe3TbN2IuCAiZkfEbGlKVeMGgLaQvwD0Qjuz8CzpYklLI+LchvYZDYt9TNKS6ocHAJ0jfwHolXZm4X1Q0vGS7re9OLV9UdInbM+SFJKWSTqpJyMEgM6RvwD0RDuz8O6Q5CZdN1Y/HACoDvkLQK9wJXIAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADIRAEFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACDTqAWU7U1s32X7XtsP2D4ztb/D9gLbj9i+wvZGvR8uAOQhhwHohXaOQL0q6eCIeJ+kWZIOsb2vpK9KOi8iflPSs5I+1bthAkDHyGEAKjdqARWFF9OPk9MjJB0s6arUPl/S0T0ZIQB0gRwGoBfaOgfK9oa2F0taI+kWSY9Kei4i1qZFlkvaoTdDBIDukMMAVK2tAioiXo+IWZJ2lLSPpN3aDWB7ru2FthdKL3c4TADoXKc5jPwFoEzWLLyIeE7SrZL2kzTV9qTUtaOkFSXrXBARsyNitjSlq8ECQDdycxj5C0CZdmbhbWd7anq+qaSPSFqqIgkdkxabI+m6Xg0SADpFDgPQC5NGX0QzJM23vaGKguvKiLjB9oOSLrf9ZUn/IuniHo4TADpFDgNQuVELqIi4T9KeTdofU3EuAQAMLHIYgF7gSuQAAACZKKAAAAAyUUABAABkooACAADIRAEFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZKKAAgAAyEQBBQAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACDTqAWU7U1s32X7XtsP2D4ztV9q+xe2F6fHrN4PFwDaR/4C0CuT2ljmVUkHR8SLtidLusP2Tanv8xFxVe+GBwBdIX8B6IlRC6iICEkvph8np0f0clAAUAXyF4BeaescKNsb2l4saY2kWyJiQeo6y/Z9ts+zvXHPRgkAHSJ/AeiFtgqoiHg9ImZJ2lHSPrb3kPQFSbtJer+kaZJOa7au7bm2F9peKL1c0bABoD3kLwC9kDULLyKek3SrpEMiYmUUXpX0TUn7lKxzQUTMjojZ0pTuRwwAHSB/AahSO7PwtrM9NT3fVNJHJD1ke0Zqs6SjJS3p5UABIBf5C0CvtDMLb4ak+bY3VFFwXRkRN9j+ke3tJFnSYkkn93CcANAJ8heAnmhnFt59kvZs0n5wT0YEABUhfwHoFa5EDgAAkIkCCgAAIBMFFAAAQCYKKAAAgEwUUAAAAJkooAAAADJRQAEAAGSigAIAAMhEAQUAAJCJAgoAACATBRQAAEAmCigAAIBMFFAAAACZKKAAAAAyUUABAABkooACAADIRAEFAACQiQIKAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZKKAAgAAyOSIqC+Y/aSkx9OP20p6qrbg65vI8Sfyvvc7/kTc97dHxHY1x6wc+Yv4AxB7oscfqPxVawG1TmB7YUTM7kvwCR5/Iu97v+NP5H0fT/r9OhKf9/BEjN/vfR+Jr/AAAAAyUUABAABk6mcBdUEfY0/0+BN53/sdfyLv+3jS79eR+BMz9kSP3+99X0ffzoECAAAYq/gKDwAAIFNfCijbh9j+me1HbJ9ec+xltu+3vdj2whriXWJ7je0lDW3TbN9i++H079Y1xx+yvSK9BottH9aj2DvZvtX2g7YfsP2Z1F7L/reIX9f+b2L7Ltv3pvhnpvZ32F6Q/v6vsL1RjbEvtf2Lhn2fVXXs8a6f+SvFnzA5rJ/5K8XqWw6byPlrlPiDk8MiotaHpA0lPSrpnZI2knSvpN1rjL9M0rY1xjtA0l6SljS0/bmk09Pz0yV9teb4Q5I+V8O+z5C0V3q+haSfS9q9rv1vEb+u/bekzdPzyZIWSNpX0pWSjkvtfyvpD2uMfamkY3q97+P10e/8lcYwYXJYP/NXitW3HDaR89co8Qcmh/XjCNQ+kh6JiMci4t8kXS7pqD6MoxYRcbukZ0Y0HyVpfno+X9LRNcevRUSsjIh70vMXJC2VtINq2v8W8WsRhRfTj5PTIyQdLOmq1N6T/W8RG92ZUPlL6m8O62f+SvH7lsMmcv4aJf7A6EcBtYOkJxp+Xq4a/yhU/AJutr3I9twa4zaaHhEr0/NVkqb3YQyn2L4vHSLv2VeIw2zPlLSnik8Rte//iPhSTftve0PbiyWtkXSLiqMXz0XE2rRIz/7+R8aOiOF9Pyvt+3m2N+5F7HGs3/lLIodJNecvqb85bCLmr2bxBy2HTcSTyPePiL0kHSrpj2wf0M/BRHF8su6q+uuSdpE0S9JKSef0MpjtzSVdLenUiHi+sa+O/W8Sv7b9j4jXI2KWpB1VHL3YrVexRottew9JX0hjeL+kaZJOq2s8qMxEz2G15i+pvzlsouavZvEHLYf1o4BaIWmnhp93TG21iIgV6d81kq5R8UdRt9W2Z0hS+ndNncEjYnX6w3xD0oXq4Wtge7KKN/93IuJ7qbm2/W8Wv879HxYRz0m6VdJ+kqbanpS6ev733xD7kPS1QETEq5K+qf78/Y9lfc1fEjms7vdvP3MY+Wu9+AOVw/pRQN0tadd0Jv9Gko6TdH0dgW1vZnuL4eeSPippSeu1euJ6SXPS8zmSrqsz+PAbP/mYevQa2LakiyUtjYhzG7pq2f+y+DXu/3a2p6bnm0r6iIrzGG6VdExarCf7XxL7oYakbxXnLvTj738s61v+kshhUn3v3xSrbzlsIuevFvEHK4dVdTZ6zkPSYSpmFDwq6Us1xn2nilkz90p6oI7Yki5TcZj1NRXfF39K0jaSfijpYUn/T9K0muN/S9L9ku5TkQhm9Cj2/ioObd8naXF6HFbX/reIX9f+v1fSv6Q4SyT974a/w7skPSLpu5I2rjH2j9K+L5H0baVZLjyyXtu+5K+Gv50Jk8P6mb9S/L7lsImcv0aJPzA5jCuRAwAAZJqIJ5EDAAB0hQIKAAAgEwUUAABAJgooAACATBRQAAAAmSigAAAAMlFAAQAAZKKAAgAAyPT/ARoaCgOaN7OYAAAAAElFTkSuQmCC\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eXBkWX3n+zk3N22pfclM7aoqVZWkqtJW3TQm2iwGYzwDAzgG2zCYxWZswzB4cMy0Bzu6mnC8aRxuG7w9Gx4843mAA2Yc0zOOpjEGykMz7mpJpaW0L6V935WScs/z/lDd2ykpU8qUUiWp6nwiKkrLyXNPZnV/77nf81uElBKFQqFQnF+0016AQqFQKI6HEnKFQqE45yghVygUinOOEnKFQqE45yghVygUinOO+TQump+fLysqKk7j0orHgLa2tiUpZcEpXV6FgSlOEhHth6ci5BUVFbS2tp7GpRWPAUKI8dNeg0LxMEmKtSKE+C0hRI8QolsI8S0hREoy5lUoFArF4RxbyIUQxcCngWYpZR1gAn7xuPMqFAqFIj6SddhpBlKFEGYgDZhJ0rwKhUKhOIRjC7mUchr4Q2ACmAXWpZT/sHecEOITQohWIUTr4uLicS+rUCgUigckw1rJAd4DVAIuIF0I8aG946SUX5ZSNkspmwsKTiugQKFQKB49kmGt/AwwKqVclFIGgL8D3piEeRVnhFu3TnsFCoXiIJIh5BPAG4QQaUIIAbwN6EvCvIozwnPPnfYKFArFQSTDI78D/DfgLnDvwZxfPu68CsVBqKcEheJ1khK1IqV8Vkp5RUpZJ6X8N1JKXzLmVZwet26BEDt/4PWvz4qAqqcEheJ1xGk0lmhubpYqs/P8IASctf4jB61JCNEmpWx+uCsyOGOflOIRI2qKviqapTg3nPWnBIXitDiVWiuK88Wzz572Cna4det10T6LTwkKxWmhhFxxKGrHqziIcDiMx+NBCIHVasVkMiFEVAdAcUIoIVecS87KU8LjjpQSv99PKBQiHA7j9/sRQmCz2bBYLGiapkT9IaCEXHEuUU8Jp48u4lJKhBCGaEsp8Xg8eDwezGYzVqvVEHXFyaA+WYVCkTC6iIfDYTRNIzL6TQiB2WzGZDIZtsu9e/fY3NwkEAhwGpFyjzpqR65QKBJir4i73W7u3r2LzWbD4XBQUFCAxWJBCGHYKsvLy1RUVLC1tQWA1WpVfnoSUUKuUCjiJlLEhRBsbm7S1dXFjRs3EEKwsLBAR0cHqampOBwOcnNzDUtF/1ufQ/fTdVFXfvrRUUKuUCjiYq+IezweOjs7uX79OjabDSkl5eXllJWVsbm5ydzcHCMjI2RnZxMKhQwvXQiByWQy5vT5fPh8PjRN2yXqivhRQq5QKA5lr4j7fD46Ojqoq6vDbrfj9/uNsUII7HY7drsdKSUrKyvMz8/T0tJCQUEBDoeD1NRUY6zJZEJKiZQSr9eL1+vFZDKpQ9IEUEKuUCgORErJ4uIiNpsNm81GIBDg7t27XL16laysrANfK4QgLy+P1NRUGhoaWFxcZGBggHA4TFFREYWFhfv8dCkl4XCY7e1thBBYLBasVitms1lZLzFQQq5QKGIipSQQCDA9PY3D4cBkMtHW1kZ1dTU5OTnGOD3s8CDMZjNOpxOn04nP52N+fp7Ozk7jkDQvL8/wyfVwRv36fr+fQCCA2WwmMzNTHZLuQQm5QqGIipSSYDBIKBRCCEEoFOLu3btcuHCB/Pz8Y81ts9koKyvb5aePjo6SlZVFUVERWVlZ+/z0xcVF44BU0zQj6Uj//eOMEnKFQrEPXcSDwaCx2x4aGqKyspLCwsKkXisjI4OLFy8ipWR1dZWZmRkGBwfJz8/H4XCQlpZmjNU0DbPZrJKO9qCEXKFQ7CKaiC8uLlJUVITT6Tyx6wohyM3NJTc3l1AoxOLiIkNDQwSDQYqKiqImHel+eqSo22y2x85PV0KuUCgMool4Z2cnaWlpx7ZTEsFkMuFwOHA4HPj9fubn55mcnMRkMmGz2cjLyzN88shD0lAo9FgmHSkhVygUBpEiDtDd3U1WVhbBYPDUUuutViulpaWYTCY8Hg9ut5uxsTHsdjsOh4Ps7GxD0CNF/XFKOlJCrlAogP0i3tvbS2pqKlVVVQwNDZ2JGik2m42SkhKqqqpYW1tjbm6OoaEh8vLycDgcpKenA+w6JH0cko6UkCsUCoLBIIFAwBDxgYEBNE3j4sWLQHzhhSeNnhmqrycnJ4ecnBxCoRBLS0sMDw8TCAQoKiqiqKgIq9VqjI1MOpqeniYcDuNwOB6ZQ1Il5ArFY06kiAshjAPG2traXcJ5VjGZTIZ4+/1+FhYWuHfvHmazGYfDQX5+/i4/PRgMEg6HH6mkIyXkCsVjzF4RHx0dZXt7m+vXr+8TtNPekceD1WqlpKSEkpIStra2mJ+fZ3x8fJefHlk/fW/SkaZphqifp0NSJeQKxWNKKBQyYsOFEExMTLC2tmZUMozkrFkr8ZCenk5VVRWVlZWsr68bfrrFYiE7O9sYt9dP9/v9+Hy+XfVeznrS0fk2hhSPDaojUHLR27LNzc0hhGB6epqFhQVu3LgR1S8+C0Kur+Mor8nOzubKlSs0NTWRkpLC4uIira2tTExM4PP5do01mUy7ko42NjZwu934fD7C4XAy307SUEKuOBc899xpr+DRIRwO4/P5DFGcm5tjenqahoaGmId+58ViOAyTyUR6ejplZWXGTau7u5uOjg7m5uYIhULGWD2LNLLT0cbGxpnsdKSEXJEwand8fokUcSEEgUCA8fFxGhsbD7UPTlu4knV93aKxWCyUlJTQ1NREdXU1Ho+HtrY2ent7WVlZMa6n++kmkwlN04yko/X1dZaXl081xl5HCbkiYR7W7vjWLRBi5w+8/rW6kRyNyC73QgiWl5fx+/00NDRgNh98XHaerZW9RPPa09LSqKys5ObNmxQXF7O4uEhLSwvDw8Nsbm7uur4u6kIIOjs72dzcZGNjA6/XazTQeNiow07FmeXWrddFWwiQcvfPFPGjizjsiNHq6ioDAwOkpqYa8dYHcVaEPBkcdGgqhCArK4usrCzC4TDLy8uMjo7i9XopLCzE4XBgs9mMeXRRj2yKcRpJR0m5ihAiWwjx34QQ/UKIPiHEU8mYV3F2OCu7Y+WVJ44u4rqAra+v09fXR2NjY9xC87gIeSSaplFQUMC1a9eor6/HbDbT09NDe3s7s7OzuxKo9ENSPZzR6/XuOiQ96c8uWbeLLwEvSymvADeAviTNqzgj3Lq1syPW/3vUvz6qkCf6umefPdp1Hnf0cDp99+h2u+nu7qa+vp6UlJS45zlM+JaWllhdXT1RwUo0/DCZWCwWiouLaWxs5MqVK3i9Xjo7O9ne3mZ5edmIZtnrp+tJR+vr62xtbbG0tLSrLV6yOLaQCyGygKeBrwJIKf1SyrXjzqt4tElkZ33r1s74034aOG/sFfGtrS2j431kje9E5ovGzMwMk5OTLCws0NLSwv379/F4PMdd/olx3BtCamoqlZWV3Lhxg9TUVJaXl2lpaWFoaAi3273vkNRsNqNpGoFAgM9+9rO0t7cn660YJMMjrwQWgf9XCHEDaAP+vZRyK3KQEOITwCcAysrKknBZxWnxsHfH0bxyxcFENkvWNA2Px0NHRwfXrl0jIyNj39jDhC2WtbK4uMjk5KSRRBQOh3f15XQ6nRQUFCT1vR2XZO3spZRYLBaqq6sJh8OsrKwwPj6Ox+MxmkzrTz269eL3+w2PPZkkw1oxA43A/y2lbAC2gGf2DpJSfllK2SylbD5r/7CKxDiOnXIWfPZHnb0i7vV6aW9vp7a2lszMzF1j4xW0aEK+urrK0NDQrqgXvY54fX09NTU1+Hw+7t69i8fjObb1kkwBTtY8+hmDpmnk5+dTV1dHQ0MDVquVvr4+2tvbmZmZIRgMAuDz+RKytOIlGUI+BUxJKe88+P6/sSPsCsUukuGzK6/8YKSUzM7OGmGGfr+f9vZ2rly5sistXSfeQ8y9wud2u+nt7aWxsTFm1EtKSgoVFRXcvHkTi8XC3NwcLS0tjI6Onqr1ctI3BLPZjMvloqGhgatXrxr/Bi+//DJbW1uHxutPTk7ylre8hZqaGmpra/nSl74EwMrKCkKI7wshhh78bXS/PraQSynngEkhxOUHP3ob0HvceRWKaKjde2z0nfj4+LhRW/zu3btcvHiR3NzcqK/R7ZDDiBT87e1tw2uPZ3ept2W7evUqTU1NpKam0t/fb0R/RGZTHvb+kkGyhFx/4jkI/WbW3NxMWVkZMzMzvOc97+HHP/5xzNeYzWZeeOEFent7efXVV/nzP/9zent7ef755wF+IKW8BPyACOcjWVEr/w74hhCiC6gH/q8kzatIMmdFCNXOOrnstVN0Ea+srDzQo04krFBv0NDR0UFdXd0+rz0edOuloaHBiP5oa2ujr6+PtbW1Q9dy1qyVRKypmpoacnJyuH37Njdv3ow51ul00ti4Y2rY7XauXr3K9PQ0L774IsDXHwz7OvCv9NckRcillB0P/O/rUsp/JaVcTca8iuRzVuKwz8oN5VFAL8MaDod3dfcpLS2lqKjowNcmYq2EQiHa29u5fPkyWVlZx163Hv1x8+ZNHA4Hs7OzhvXi9XqPPX8sHuaOfC8+n4+MjIy4ffKxsTHa29t58sknmZ+fR0o5++BXc4Dxj6syOxWKc4wu4qFQyBDltbU1iouLcblch74+XiHXvffLly+Tl5cXc8xRqxPq3X6CwSCLi4v09e2kouhRL3r2ZDI4zUPTRKJWNjc3ef/7388Xv/jFfYfUUkophDA+EFVr5TFARYs8mugd73URB+jq6sJms1FYWBjXHPEIeTgcZmxsjIyMDBwOx7HXfRBmsxmn02lYL9vb24b1cpK79KNwlB05ENdrAoEA73//+/ngBz/I+973PgCKiooQQjgB/e8FY86EV6E4dyQ7K1Nx+ugivrfjvd1uN+qExIOeUn7QdXp7e0lLS8Nutydl7fGiN37WrRe3283w8DBjY2PHEvWzFsYYbd6Pf/zjXL16lf/wH/6D8fN3v/vdAL/y4NtfAV7Uf6eEXKE4Z0QT8b6+PqxWK1VVVXFHosDhO/LBwUFMJhNOp/PUaq3o1ktBQQFVVVVYrVZ6e3uj1hCPh9P0yOP5DH/yk5/wX//rf+WHP/wh9fX11NfX89JLL/HMM88AvF0IMQT8DPC8/hrlkT9mqGiR8000ER8cHASgurraSAuPV3QPEv3R0VF8Ph/Xrl3bVZ/7tJBSYjabyc/Px+Vysb29zdzcHG1tbWRmZuJ0OsnMzDxUpE9rRx7v5/emN70p5lgp5dui/VwJ+WOGslPON5EiLoRgZGQEv99PXV2dISqJCnm0sdPT06ysrNDQ0GBc67SFfC9paWlGT87V1VWmpqbY3t7eV252L6e5I4eT6bakhFyhOCfsFfGxsTHcbve+ZsnHtVYWFhaYmpqiubnZEKqz3OpNCEFubi65ubkEAwE2vvc95l58EX9pKZlvfCP5+fm7silP0yM/qc9RCblCcQ4IBoNG/WshBJOTk6ysrFBfX79PHBLZPe8du7KywsjICM3NzftSyU97R36ocEpJ6h/9EfaXXwZ2dszTCwu01teTnZ2N0+nEbrfvmsf8v/4Xts99Dra3Cbz//fh///fBYolrPZFx+/Gu/6RQQq5QnHGCwSCdnZ3U1tYihGBmZob5+fmYzZL1OtjxEGnDbGxs0NfXR1NTE5Y9YnYWrZW9aENDWL73PcKFhaBpCL+f0u98h7wPf5hVj4fJyUm2t7eNm2JqWxspv/qriAd1X6x//dcgBP7nnz/4Qg+QUh5aNyWSQCCw73NNFipqRXEgylM/XXTRWV9fRwjB/Pw8U1NT1NfXxxSRo1gr29vb3Lt3L2bDifMg5LjdSJMJ9JubxYIIhxEeD7m5udTW1lJfXw/sRPmsfu1rhogDCI8Hy9/9XdyXS9Ra8Xq9J1LCFpSQKw7hrKT07+VxuMGEQqFddsrCwgKjo6OHNktO9LDT6/UatcrT09NjjjtsnpPmMOEMV1WBzYZYW4NQCLGwsPOzHKNIIBaLBZvNRn19PZklJYT3fI4ygYYbkYedWn8/pu9/H621FWLcRH0+nxJyhSKSs3qDSRZ7O96Hw2GGh4dpbGw89PE8EWtFn/fKlSv70sAjORc78pwcfC+8gCwoQKyuEr5+Hd9/+S+vpzQ/QL8hiF//dcjORlosSCBks9Hz0Y8yMTGx045tawvtzh1M//zPsLy873L6PKbvfQ/r5z+P5a//GusLL2D5y7+M2v3kJIVceeSKfeit1XT0/w+effbx2AmfNuFwGJ/PZ4j42toaHo+Hp556Kqkd70OhEAsLCxQXF8cscxvJmRdyIHz5Mt6/+ZsDx+gCLIuK2P7nf8by9a/DxgbBf/kvKWlsZH5+nt5XX6Xqb/6GrK0tLDYbZrsd/zPPIIuLX79WOIwWCGD55jeRTidYrRAOY/rJTwi+853Iqqpd1z2p7kCgduSKKJzVlP7HoWbMXhHf2Nigp6eHzMzMuA/K4rFWwuEwHR0dRkr/YZyFHflJFM2SRUX4/+N/xP/7v0/4ySexWCyUlJTQtL5Ort/PZmEh01YrG8vLhL797X3zaMHgjpWi/9toGphMCJ9v33WVtaJQcHZvMMlEL4AlhGBzc9M4gLRYLAkdYB40VkpJd3c3ubm5ZGVlxV3G9qBxbrc7rnrixyVhL35jA9unP01aczOp73wnWktLXPMItxtTWho5OTk4nU6smZmsj4/T2trK5OSkUftd2O2EL19GTE+D349YWEDa7YRLS/fNqQ47FaeGSul/uGiahhCC7e1tOjs7uXHjBunp6QmHFMYaK6Wkv78fm81GZWXlkVu9ReJ2u+nu7jZauY2Pj+OLsiM9DWyf/Szm738fhEDMzJD6a79Gyvz8oa8LX7u2E9GyvY3w+0n1eMj/+Z83kq+6urpYWlpifWMD76c+RejmTdjeJlxejv8//2eI0nTjJK0V5ZErDuSs7nYf5RuMHkUS2YUnESE/SJzv379PMBikrq7OmDeRDkHR1trV1cW1a9ew2WyEw2Hm5+e5d+8eVqsVl8sVl/9+1OsfSDiM+ZVXkJmZO5aH2Qybm9j7++Ff/IuDX3rtGoFf+zXMf/d34PUSeN/7CL3tbVg0jZKSEkpKSujq6mJzc5OWuTlyf+7ncH7sYwd2TVKHnQrFHs7qDea46F3nr169usu7TkRwY4n+5OQk6+vru7JBE9mR7x0XCARob2+npqaGjIwMAoEAZrOZ4uJiiouL2dzcZGZmhpGREQKBAB6Ph9TU1Ljew0HrSGAwMjUVgsGdg8gHXlwozjWE3vQmQm96U8zfm81mysrKSEtLY2VlxSgyVlRURFFR0b6DaSXkCsVjwubmJpcvXyYnIvYZjl8/ZW5ujtnZWZqamnZlgx6l+TK8flhaWVlJTk5O1FKyGRkZVFdXEwqFuHPnDv39/QC4XK599U9OBCHw/97vYfu934OtLTCZCF2/zmpDQ1Km1+PINU0jPz+f/Px8/H4/8/PzRoMPp9NJbm4umqYpIVcoHhfy8/MJBAL7fn4cj3x5eZnR0dGo9VOO4pFLKbl37x4FBQVxdQwymUxYLBbq6+vxeDzMzs4yNjZGTk4OLpcroSbOiR52Bt/7XsLl5Zju3kXm5RF817uQXV0JzRGLaAlKVquV0tJSSktLcbvdzM3NMTIyQm5uLrOzs4cK+cc+9jH+/u//nsLCQrq7uwG4desWX/nKVygoKKCzs7MD+M9SypciX6eEXKE4BxzVI19fX2dgYCBq/ZS9Y+Odc3BwEJvNRkVFRfxv4AFpaWlcuHCByspKVlZWuH//PoFAAIfDQVFR0YEZq0eNiAk3NhJ+0JU+mUgpDyxja7fbsdvthMNhFhYW+Pa3v838/DzV1dV84AMfiPqaj3zkI3zqU5/iwx/+8K6f/9Zv/Ra//du/DVAf7XVKyBWKc8BRduSbm5t0d3fT0NAQcyeYiGUjpWR8fByv18v169fjXnvk6yNrput2hM/nY25ujrt372K323G5XHE1iHhYiOVlrF/4AlpfH+HKyp3EIJcr7uqHmqbhcDj44Ac/yObmJlV7EoUiefrppxkbG0t4jSr8UKE4ByS6I9crJl6/fp20A+qHJLIj9/l8LCwscO3ataTW4bbZbJSXlxu9Oaenp2ltbX09VT7OeU6EYBDbpz+N6X//b/D7MbW2kvLJT4LXm3AZW5/PR15eHjdv3kx4GX/2Z3/G9evXEUJ8TQiRs/f3SsgVMXlUI0POMrGEIREhD4VCrKysUFNTc2jD5HiFfH19HY/HQ319/ZG64sSD3puzpqbGuE5XVxf37t1jeXk57vef1DXNzaGNjSELCsBmQ+bnIxYX0UZHD7VW9uL3+6NWljyM3/iN32BkZISOjg6AWeCFvWOUkCti8qgXpjpPxCvkwWCQ3t5e0tPT90W+xJr3MCHf3Nykv7+fjIyME6unvRc9Vb65uZmKigqWlpZYXFxkZmYGr9f7UNYAQErKTgq+/tk/+FqmpDy0MrZFRUWYTCb9pvEV4Im9Y5SQKxTngHiEXA8JdLlccQvuYTtyn89HZ2fnkeyUZGG327l8+TL5+fnYbDZ6e3vp6OhgYWHhxHfpMj+fwPvfj1haQszPIxYXCf7MzyArKhIW8qPuyGdnZyO/fS/QvXeMOuxU7EJVPjybaJoWNSxRRw8JzM/Px+l0srCwENe8Bx12BoNB2tvbuXLlitEi7agko+iWEIL8/HwqKyvZ3t5mZmaG0dFRcnNzcblcO7XU/X7ExgYyKyvulm2HEfjsZwk3NqKNjBAuKyP09rcb/2Mk6pEftiP/pV/6JW7fvs3S0hIlJSU899xz3L59m46ODv1abwH+7d7XKSFX7OLWrdcFW4ioZZUVp8BBQiilpK+vj7S0NCoqKggGg8dO59d39+Xl5eTl5SGlPFCIH3Z1xLS0NC5evEhVVRXLy8sMDw+T0tfHlW98A1soBBkZ+J99lnBt7fEvJgSht76V0Fvfeqxp4qm18q1vfWvfzz7+8Y9HfvvuaK9T1opCcQ44yFoZGRlBSsnFixeB4zVfhp0bQ09PD3l5eTidTmNcPHOdJNHek6ZpFBQUcKOykmt/+7eEgWWrlS23G+13fxe5vZ30dYj797H8wR9w4a/+CtN//++wJ7ImFucis1MIYQJagWkp5cEVaRTngke5MNVZJdGolYmJCTY3N42qfAeNjXW9vQI5PDyM2Ww+UsLPSRPr8xHz85hCIVILC0llZ/frn5qi54c/JPvGDRwOBxaL5fj2zuIi1t//fbSJCezb21g2NxGhEMF//a8Pfe25EHLg3wN9QOx+UYpzhfLEzw7RxHl2dpb5+Xmampp2CVwiO+O9USuTk5Nsbm7uKqyVLE7SepF5eTtf+Hxgs2ENh7HZ7Vx54xuZ29qio6OD1NRU4wnjqJh+/GNMt28jLBbSAgFMS0tgtcYt5Ec57IyHpFgrQogS4OeB/ycZ8ykUjzPRBHSvkC8tLTE+Pk5DQ8Ox4rojd+QLCwvMzs7qiSdHnjPWdU6U3Fz8/+7fIdbWEAsLiPV1/J/9LJbcXEpLS2lubqasrIyFhQW2trYYGxs7Us107Sc/QYTDyPR0QikpSCHQRkfjeu1Ro1biIVk78i8C/xGImX0ghPgE8AmAsrKyJF1WoXg8iBTytbU1BgcHaW5uPrA2STzoUStra2sMDw9HLax1Vjgs3C/0jnfgvXEDMT+PdDiQhYXG74QQZGZmkpaWxvb2Nlarle7ubsxmMy6Xi7y8vPhuiLm5O7Hl29togQACCMVpQZ1pa0UI8S+ABSllmxDizbHGSSm/DHwZoLm5WcVCKBQJEFk/paenh8bGxrgaMR+GEIJAIJDUOU8TWVSELCqK/fsH2ZgulwuXy8XW1hYzMzPcv3+fvLw8XC7XgSUNgu95D+Z/+Afweglub2NNTyfwb/dFA0blTAs58FPAu4UQ7wJSgEwhxP8npfxQEuZWKBTsCLnf7zfavx23QYNOIBBgZWWFmzdvJm3Os8xenz49PZ1Lly4RDodZWlpicHCQcDiM0+mksLBw39NJuK4O75/8CeZvfpO1mRlMv/7rhH7qp+K69kl65McWcinl7wC/A/BgR/7bSsQViuQSDAZZXV2lubk5ofrdh805MDCA3W4nM/NkYxSSEWeeaCZlLKJZKJqmUVhYSGFhIV6vl9nZWdra2sjKysLlcu2qWRO+cQNPXR3D7e1kNzfHfd3z4JErFIoksVf0IuunZGdnJ+Ua4XCYzs5OHA4HGxsbx55venqa+fl5iouLyc3NPTMlaPcSz80kJSWFyspKKioqWFlZMUr36jXT9TDGRA+ZQ6HQsc80YpHUWaWUt4HbyZxToXicCYfDtLe3U1paytzcXFLmlFLS29tLdnY2DoeD9fX1Y823tLTEzMwMlZWVzM/PMzIyQmFhIU6n88Q84aOSyK5eCEFeXh55eXn4/X7m5ubo6OggLS2NwoiD1JO4dqKoHblCcUaRUtLZ2UlRURFOp5OZmZmkzHv//n2EEFRVVeHz+Y5lebjdbqMDEUB2djahUIj5+Xm6u7uxWq0UFxcnJYY8GUJ41DmsVitlZWWUlpaysbHB5OQkbreb8fFxHA7Hqd+wlJArFGcQPU3ebrdTVlZGOBxOWAyjidbU1BTr6+s0NDQghDiWd+3z+ejq6uLGjRtYrVajCYTJZDKiQtxuN9PT07jdbiYnJykpKTnVyJjj3gyEEGRlZWGz2QgGg5jNZrq7u7FYLLhcLqPRcrTrniRKyBWKM8jQ0BCapnHhwgUgsZZs+vi9orW4uMj09DTNzc3GzxOdVycUChmVETMyMmLOYbfbuXLlCh6PB5PJRFdXFykpKRQXF5Odnf3QvfRk2RvhcBiz2UxxcTHFxcVsbm4yMzPDyMiIUYEyWhjjSb1fVTRLoThjjI6O4vF4uHr16i7BTYS9maDr6+sMDQ3R0NCwK6TuKDtyKSVdXV2UlJSQp6fGH4IQAqfTSVNTE2VlZczNzcVs53bQdZNmrWxuYv3CF0h9//tJ+bVfQ+veV+I7obVkZGRQXV1tRNeUhTEAACAASURBVBUNDAzQ3t7O3NwcoVAIiO/f8GMf+xiFhYXU1dUZP1tZWeHtb387ly5dQgjxfdXqTaE4B9jt9mM3coisobK9vU13dzf19fX7bI2jCPng4CDp6emUlJTsm+sw9AzLq1evGuUFOjs76enpYW1t7cQtCF2ArV/8IqYf/QiZng6rq9h+93cRuxs4HEg4HI5qoZhMJoqKimhoaODKlStsb2/T2trKN77xjbiafXzkIx/h5Zdf3vWz559/nre97W0MDQ0B/AB4Zu/rlJArFGeMgoKCY/fF1C0Tv99PR0cHdXV1MR/1ExHPyclJPB4Ply5dSng9ezGbzUY7t5KSEmZmZmhtbWVycvLAJhrHRkpM/+f/IJ3OneYTWVng96MNDCQwxeFPB6mpqVRVVfHEE0+QnZ3NxMQEb3rTmw6s8fL000+Tm5u762cvvvgiv/Irv6J/+3XgX+19nfLIFYpHECEEwWCQ7u5uLl26RFZWVsxx8Qp5MBhkZmZml8eeCLGuox8gZmVlEQgEjDC/9PR0iouLyczMPLLFFG0NQghkRgZ4vZCWttM9RUrkAan5e4m1I4+GEIK3vOUtVFZW8vd///cJR7jMz89HVm2cA/bVIFBCrlA8gggh6O3tpbi4mIKCggPHxSPkbrcbr9fLG97whhMtqmWxWCgtLaWkpIT19XWmpqbY3t7G6XQaXvNxkFIiNI3Apz6F9fnnYXUVgFBDA+HGxsTmOUKbt+MmdEkppRBi3z+YEnKF4hFDSsnm5iZ5eXn7fOy9xCNGephhampq3E2dj4sQguzsbLKzswkEAszOzrKxscHQ0BBlZWXY7fYjPxUIIQg9/TRepxPT4CDSbif01FOQQNZlIjtyOF7BrKKiImZnZ3E6nQghnMC+hqzKI1cozhjHtQ9GH9THLi0tPfZaIsMMT6u8rcVioaysjKysLAoLC5mYmKCtrY3p6WmCwWBCc0XupOWlSwR//ucJPf10wo2aj7ojPwrvfve7+frXv65/+yvAi3vHKCFXKM4J8dggMzMzrK6ukpOzL0ItYfaGGR52/YfRnDk7O5u6ujquX79OMBjk7t279Pf3x10vJllx5Ccl5L/0S7/EU089xcDAACUlJXz1q1/lmWee4fvf/75+wPwzwPN7X6esFYXinKBHosTaGS8vLzMxMUFzczMDAwNHSvSJZG+YYbLEOBlYrVbKy8spKytjdXWV8fFxfD4fLpeLwsLCmMWpkpkQlKi1Ek9G67e+9a2oP//BD36gf/kz0X6vhFyhOCfoST7RhDyy5onZbN7XizNR9DDDGzduGD87TMhPOkszmggLIcjNzSU3Nxefz8fs7Cx3796NWn421hxHIRwOJzTPSZawBSXkCsW5IVoDZgCPx0NXVxf19fXG4/tRU+/h9WqGiYQZnoWytTabjYqKCsrLy1lZWWF0dJRAIGDs0k0mU9KeKBItY+v1ek+0sJYScoXijBFLFKMJeSAQoKOjg9raWtLT03fNcRQh13f2N2/e3LfzP4618jBtmcjysz6fj5mZGdra2sjOziY1NfVUPHK/36+EXKFQ7BdyPaLkwoUL++KTj2Kt7K1muJfT9siPYovYbDajScTy8rKxS09JSYnayi1e9KJZ8aJ25AqFAtgt5FJK7t27h8PhiNrkIFEh31vNMBqnLeTHQQhBfn4+oVCIzc1NPB4Pra2t5Obm4nK5dj3NxIPakSsUiiMRKeQDAwOkpaVRVlYWdWwi1kq81QzPs5DrSCmNuHR9lz48PEwoFDK89Hi876NErajDToVCYQj52NgYfr+fa9euHTg2XtH1+XwUFBQcmgV6XM7CTSByJ61pGgUFBRQUFODxeJiZmaGlpYW8vDxcLlfUImPR5okHtSNXKB4zDjrsXFpaYn19ncbGxgOFJN4d+eTkJOFwOK5qhgftyDc3N5mcnMThcES1Zk7jgDGROVJTU7lw4QKVlZUsLS0xODiIlJLi4mLy8/P37b4fZop+PCghVyjOCV6vF7fbzVNPPXWoiMQKVYxEDzNMTU2N6/qxhNzv99PV1UVxcTEjIyOEw+GYAnjaHHYz0DSNwsJCCgsL2d7eZmZmhtHRUfLz83G5XMZnlehNxev1kpmZeez1x0IJuUJxDnC73aysrFBdXR1XtMRhO/LIMMO2tra4hCmakIfDYTo7O41SuS6XC4/Hw/T0NKOjoxQUFFBcXBzfm3wIJBL/nZaWxsWLF6mqqmJxcZH+/n6EELhcLkKhkLJWFApF/Hi9Xrq6unA6nXGHvGmaFrM5w94ww0T89MhxUkr6+vrIz8+noKDAaNmWmprKxYsXqaysZGFhge7ubrxeL+vr68eK406WtZIomqZRVFREUVER29vbTE9PMz8/TygUIiUlJa4nGiXkCsVjTCAQoL29nZqaGtbW1uKORIllg0QLM4w3GmWviE5MTBAOh6moqIg63mQy4XQ6cTqdtLe3s7S0xNTUFE6nE4fDkVAcdjI5zs0gLS2NS5cuEQgESEtLo6+vD5PJhMvlIi8vL+ZuX0WtKBSPKeFwmI6ODiorK8nJyWFjYyNuIY+2y44VZhjvwWik4C8tLTE/P2+k8R92I7BarVRUVGC1WnfVQykuLo4Zt34SJKvWCkB+fj4VFRVsbW0xPT3N/fv3KSwsxOl07hNtddipUDxm6MJ47949CgoKcDgcQHwHmJFz7B17UNPkeHfkUkq2trYYGBigubk54cNMPYa7tLSUlZUV7t+/TzAYNDoZHTTfSUatHGee9PR0qqurCYVCLCws0NPTg8ViMXbpQggl5ArF48jg4KBRBErnIN97L3tFP1o1Q51EhDwQCNDd3c3169ePJUyR9VD0GO6xsTEKCgpwuVwnZkOcZBnbSCtpc3PT2KXb7XY8Hk/Cn1dFRQV2ux2TyYTZbKa1tTXmWCXkCsUZIxgMIoSgurp6188T3ZFH2iAHVTOMV8illAwODnLx4sV95WHjIdY19BjuiooKFhcX6enpwWq1UlxcTE5OTlIrKyZzR37Q00NGRgaXL18mFArx2muv8corrzAzM8OXvvQl6urq4r7Oj370I/Lz8w8dd2whF0KUAn/DTmdnCXxZSvml486rUDyumM1mLl++vE/4EhFyfazb7WZwcJCmpqaYBaLiFfL19XVycnKi1nY5jHjE02Qy4XA4cDgcuN1upqamGBkZMX4W7zwH8bDrkZtMJp566ikaGhr40Ic+lHBNl3hJRrR+EPislLIGeAPwSSFETRLmVSgUESS6Iw8Gg3R1dR1qg8Qj5HqEitPpjDnmuB2JIrHb7Vy9epX6+noA2tvbjYSo43Card4aGhqorKyM+zVCCN7xjnfQ1NTEl7/85QPHHntHLqWcBWYffO0WQvQBxUDvcedWKBSvk4iQSylZWlrixo0bh0aFHBZHvry8zOzs7IGP+OFw2PijaVrSGjVbLBZKS0spKSnh1VdfZXR0lGAwmFCBq0hOs9Vboh75K6+8QnFxMQsLC7z97W/nypUrPP3001HHJjV/VghRATQAd6L87hNCiFYhROvi4mIyL3vq3Lp12itQPA7EK+RSSoaHh0lLSzuwmqHOQTvyra0t+vv7qa+vjzlOSkk4HMZisWAymQiHwwQCAUKhUNJ26UIIzGYz169fp6amhq2tLVpaWhgZGcHr9cY9TzI7BJ10ZqeeEVtYWMh73/teXnvttZhjkybkQogM4L8Dn5FS7mtpLaX8spSyWUrZXFBQkKzLngmee+60V6B41IgmEvFmYA4ODpKWlhZ35EesOPJAIEBnZyfXrl3DZrNFXZOU0khXN5lMWCwWrFar0VYtGAwakTbJEtGUlBQuXLjAzZs3SU9Pp7e3l66uLlZWVo6U2HQUjmKtJBKJs7W1ZdhIW1tb/MM//MOBh6RJiVoRQljYEfFvSCn/LhlzKhSK3cSTuKOHGV66dInBwcG4542VPFRVVWUUe9o7ThfxvREcmqYZ34dCIQKBAF6vl1AoZIh+MoppaZq263B0enqa4eFhI3PUYrHse00yE4JOckc+Pz/Pe9/7XmAniumXf/mXeec73xlz/LE/TbHzbr4K9Ekp/+i480XjLFoXt26BEDt/4PWvz+JaFY8Gh1krepjhtWvXDIsjHqIJ+cDAAFlZWUa0yN5xup1ymDCaTCbm5+fJyMggMzOTcDhMKBQiGAwm/XD0ypUrNDQ0ANDR0UFfX9++w9FkCnkiSCkTOjeoqqqis7OTzs5Oenp6+NznPnfg+GRYKz8F/BvgrUKIjgd/3pWEeQ3OonVx6xZIufMHXv9aCbnipDhIyPVqhvX19ZhMpoS6+ewdOzU1hdfr5cKFCzHH6QebQogDhXF5eZnFxUUuX75s2C56jRV9t55MQdcPR5ubmykqKmJsbIy2tjbm5ubiuvGcV5IRtfIK8Oh9MgrFGSOWkEdWM9Qf3xONOdcFemVlhampKW7evBlT8PTd+GEi7vF4GBwcpLGx0bBSdNvFbDYTDocJBoOG5aJpWtJsFyEEubm55Obm4vV6je4/QFwJNueNs1X1ndd3tOfJunj22dNegeJxIJo4h0IhOjo69jVNPsqOfHt7m76+PmNXH22cLr6HiXgoFOLevXvU1NTE9IY1TcNqtRqHo7rtkuxdekpKClVVVdy8eROz2czo6ChdXV0sLy8/lPZzD+MaZ07IdRvlPFkXZ3FNivNNrKiVSIHTDySLi4v3hRkeJXmos7OTurq6A6MrQqFQzPVFrqunp4eSkhKysrIOvb6maYbtYrFY0DTtRARd0zRSUlKoqamhqqqKpaUlWlpamJiYiLuGzXE4SUtH1VpRKM4Je8U5VjVDSLzj/ejoKBUVFTGFV0pJSkoKExMTRtPiWMI0Pj6O1WrF5XLFfX3AsFT03bluuwSDQSOc8bjoc+i1UILBIHNzc3R0dJCRkUFxcfGJtmQ7Kc7EjvwwG0VZFwrFbnGODDOMNTZeVlZWSElJOTT9vqioiNraWlZXV3n11VcZGxvbt5NdXl5meXl5X8GvRIm0XcxmsxHlctxd+t7DTrPZTElJCc3NzTgcDsbHx2lra2N2dtZ4+jguD8NaORM78lu3XhdtIV63UyJ/r1A87ugCdFg1w0SYnp4mEAhQWloac4yeoSmEID093djJzszM0NbWRmZmJmVlZWiaxtDQ0K7DzeOiH44ODQ1RUlKCEGLXLv0opQBiVYDMyckhJycHn8/H9PQ0ra2t5OXlUVxcvKudW6LCHAgEosa0J5MzIeQKhSI+QqGQ0dThuPVMVldXmZiYOHAnHitCxWw2Gw0ilpeXGRgYYH19nQsXLiRdtGZmZgiFQpSVlRkHrnqki762yCSkg4gn/NBms1FVVUVFRQVLS0v09/ejaZpxFnGUrM6TbCoBZ8RaiUTZKApF9F2jz+czmkMcVxg8Hg+9vb1GhEqsGip6JmYs4dIbRJjNZqqqqvB4PLz66quMj48n5QDR7XYzOTlJTU2NsYbIw9HImHS/33+o7ZKICGuaRmFhIQ0NDVy4cIHl5WVaWloYHx9P6D08DCE/cztyZaMoFPvRwwxTUlKO3eMyGAzS0dFBbW0tqampMft7Rku/j8bY2BgpKSlGN6NI2yUrK4vS0tIjrTkQCNDT02Nkqu4lMiZd36EfFpN+1ISgyMPR6elptra26O3tpaSk5NDDUa/X+/jtyM8i6uaiOE0iwwyP23le7wVaXl5OdnY2EL2GSrxZkEtLS6ysrOw6dNVtlyeffJL8/HwGBwe5e/cui4uLcfvLUkp6e3uprKyMqxmDyWTaF5MeDAb3lQI4bman2WzG4XCQlZWF0+lkYmKC1tZWw/6JxlEqHyaKEvI4OIslAhSPD3o1w2hhhomil7eNDA3cK+Txpt9vbW0xPDzMtWvXou7ahRAUFBTQ2NjI5cuXWV5e5s6dO3HFbY+Pj5OSkkJRUVFC7y/SdrFarcDuUgDJauCsaRo5OTnU1dVx7do1fD4fra2tDA0Nsb29vWv8Y+mRR+PNbz7tFSgUp4MeZhgZznfUcLaZmRncbve+0MDIqorxpt8Hg0G6u7upra01BPMg0tPTuXLlCk1NTUgpaWtro7+/n62trX1jV1dXWVxcjBlaGQ96cwubzYbVajWSjPQD0uOEMO5tKmGz2aisrOTmzZtkZWUxMDBAZ2cnS0tLSCmPJOQvv/wyly9f5uLFizz//POHjj8XQv5P//Twr3meSgQoHj2EELuqGeqimmiij87a2hrj4+Ncv359n0Drc+r1ww8TcT1zs6ysLOEmzBaLhfLycp588kny8vLo7+/n7t27u0Svv78/5i7/KOgx6evr61F36YkSa1cfeTh68eJFVlZWaGlp4fvf/35C7yUUCvHJT36S7373u/T29vKtb32L3t6DG66dCyE/Dc5TiQDFo8feaoY6iaTew47oeL1eenp6uHHjRlSPPTKkT//+IEZHR0lNTT0wbPEwdNulqamJ6upqFhcXuXPnDq2trVy4cCGhJgzx4PF4uH//PnV1ddhstn2lABJJ/omn8XJ6ejrV1dU0NjaytrbGd7/7XZ555pm45n/ttde4ePEiVVVVWK1WfvEXf5EXX3zxwNecWSF/85uj74iPYrMo8VWcNzRNixpmmGhVw0AgQEdHBzU1NaSlpUUdpwt5PBEqi4uLrK2tcfHixfjeSBxkZGRw9epVsrOzSU1NZWRkhIGBgX1e81EJh8P09PRw5coVw2aJPBzVP1O/3x9Xe7p4Picds9nMm9/8Zj74wQ/ymc98Jq7XTE9P70rQKikpYXp6+uDrxDXzKXD79utfR8v2TITnnjuemKvYdsXDJj09Par3nOiOXC9elZOTE3OMxWJhfn6etLS0Axsab21tMTIyktTMTZ2FhQU8Ho/RGGJxcZG+vj40TaOsrIzc3NwjH1Lev3+f3NzcqJ+BbrvoTyT633rrumjvMxEhh9fbvEU26Ug2Z1bIzxJqR684KyQi5D6fD7vdfmC0SygUIjMzkxs3bjA1NcXo6ChOp5Pi4uJdGZrBYJB79+7FfbiZCNvb29y/f5+mpiZDrAsLCyksLGRzc5OJiQmGhoYoLi7G6XQmFIK5vLzM+vo6jY2NB47b255ub8GuSHsrHmslEr/fn5BVVFxczOTkpPH91NSU0Yg55vrjnv0U+emfTvw16rBScZ45KJMynsPOubk5QqHQvi4/kURGqKSlpVFdXc3NmzfRNM2IKtne3kZKSXd3NxUVFQkfbh5GKBSiu7ubmpqaqKn9GRkZ1NTU0NTURCgUoqWlJW7bxefzMTg4SF1dXULCGxntojfACAQCRkz6Safo37x5k6GhIUZHR/H7/fzt3/4t7373uw98zbnYkUfaLPFyWCEuheI8Es+OfH19ndHRUTIzM2OKfqz0+8gaKrq94fV6sdvtCcd0x0N/fz8ul+vQ7EiLxUJFRQVlZWUsLS0darvokTWXLl06cgy3vkuPbHqhlwI4SSE3m8382Z/9GT/7sz9LKBTiYx/7GLW1tQe/Ju7ZFQrFqXOYkHu9Xrq7u2loaGBwcPDAGioHeb1CCAoLC4GdFHxN02hpaaGkpASHw5EUj3x6ehop5aG2QSR6iF9hYaFRh0W3XVwul2GBjI2NYbfbk9LWbW97uvX1ddLT0wkEAnG1p9MtrkR417vexbveFX/r43NhrRyXeA8rle2iOOscJOShUIjOjg6uulykeb2IKIkviaTfb21tcf/+fRoaGqirq+PGjRtsb29z584d7t+/j9/vP/L7cLvdTE1NcfXq1SMfYtrtdmpqamhsbCQYDPLaa68xODjI/Pw8S0tLB9pKR2V1dRWfz0dJScmuUgAHxaQn6pEfhcdCyOMVaJWKrzjrxBJyKSXdXV1UTU6S//LLaN/+Nrm3byP3eMnxpt8HAgHu3btHXV2d4V3bbDYuXrzIE088gc1mo729nZ6eHjY3NxN6D3oxrLq6umOX4gWwWq1UVlby5JNPkp6eTnd3N5qmsba2ltSmDn6/n8HBQWpqajCZTLtKARzUnu6xrH6oUChiE0vI79+/T8bCAkVzc1BWBpqGpa0Nra0Nfu7ngPjT7/XDzcrKyqhVC00mk2FlrKysMDQ0hJSS0tJS8vPz48oKjbcYViIIIVhcXKSuro7U1NRdtovT6TzWTUNKSV9fHxcuXNglyoe1pxNCqKJZ0Ui2/aGiWxRnkVhiGE3I5+fnWV1dpSovD1JT4YG4hO12tMVFgLjT7wFGRkbiOtzUa5E3NDRw+fJllpaWuHPnDpOTkzEzJcfGxkhLSzuRg9PJyUmj0FZmZia1tbU0NDQQCAR47bXXGBoawuPxHGnumZkZLBaLcW4Qjb3t6fQD0uXl5WNXrTyMcyfkybY/VCq+4jyxV8g3NjYYGRnhxo0biNxccLvhgYia3W6CBQXIcJjw4CDaP/4j2quv7oyJwfz8PG63O2F/OT09natXr9LU1GT41UNDQ3i9XmPMysoKy8vLSc0K1dnY2GBubm5foa1I28Vut9Pd3U1nZycrKytx2y7b29tMTk5y+fLluMZHVmAE+OEPf0gwGEzsDSWIslYUijNKtJjxSCH3+Xzcu3eP+vp6LO3tmL71LeTsLKK1FVlVBTYbvitXSO/rQ/vBDxDZ2TA9DSMjhN73PtiTsr+5ucno6OiuxJxEsVgsVFZWUl5ezvz8PF1dXUZdlqGhIRoaGpKeFRoMBunt7Y3ZgAJ2PjeHw4HD4WBjY8NIMtKjcGK9Tk/vv3r1asLWjKZpvPDCC3z4wx/mfe97X8LvKxHOhZDfurV7J67/N/bss2rnrHi80Lv5hMNhOjo6uHz5Mum9vZifew6ZkYEwmxHt7Wg/+Qn5mZmE2tvhrW+F1FS01lbExgbSaiXc2IisqTHmDQQCdHd37zrcPO46nU4nDoeD1dVVurq6sNlsRgXCZIm57l2XlZXF7blnZmZSV1dnNFl+7bXXyM/Pp7S0dF90yejoKHl5eWRlZSW8tra2Nm7fvs3toyTCJMi5sFZO0v5QHrniPKFHR/T09OBwOCgYH8f0h38IAwOI3l5EWxtibg4sFgKFhZiGhzG/9BLmf/xHxOwsSIkYGUH77neNOfXDzaqqqmO3kduLXo63vLyc+vp61tfXuXPnDmNjY0np6TkzM4MQYlejjHjRmyw/+eSTZGRkcO/ePTo7O1ldXUVKydraGqurq1RWViY8t8fj4TOf+Qxf+cpXkt6MOhrnQshPkqPeJJTQK04DXRg1TaPMYkG89hrC7QaTCSwWxNbWjkfu82GTEn9qKlvT0wTm5pBWK/h84HKhTU7Cgzhw/XDzoIO8ozI/P8/29jYVFRWkpqbGLANwFDY3N5mcnOTKlSvHWqP+9HDz5k0qKiqMXfq9e/eOFOcupeS5557jQx/6EDURTz0nSVKEXAjxTiHEgBBiWAgRX9HdI3JWKhGqmHPFabC5ucnW1taOwHg8YDYjs7ORmZk7Iu31gs+HWF/H0ttLxtwctsZGtgoKmPP72bDbCV66BCkpYDId+XAzHra2thgdHaW2tjZqGYAnn3yS3Nxc+vr66OjoSOgAUn8qqa2tTWpESFZWlhG+aLfb6erqYnh4eNeh7WH8+Mc/pqenh09/+tNJW9dhHPsTEEKYgD8H3g5MAS1CiP8ppTy4pcUROcmd8Fm5SSgUsP+w0+12Mzc393qKfEYGBINQWQmhEHJubsc+0TRkKARbWwhA/MIvYJ+awt7VhdvnY7G7G9973oNtbe3Yh5ux0Ith1dbWxrQW9DIAhYWFxgHk8PBwXGUABgYGcLlcSS/iBTtPEQA3btxASmkc2qakpFBWVkZWVlbMz2tjY4P/9J/+Ey+++GJSkp3iJRk78ieAYSnlfSmlH/hb4D1JmPehE4+dovx0xWng9/u5d+8elZWVr4tIURHyTW9COp1gt8ODdm3S4UAWFYHVupOK39lJ6AMfIPyRj5DxgQ+Qe+sWoZ/6KTo6OkhNTcXn8yV1rfoBZElJSdxCqx9AxlMGYH5+nkAgkJRm1Hvxer2Mjo4alopuuzzxxBOUl5czOTlJS0sLMzMzUcsf/M7v/A6f/vSnqaioSPraDiIZzyTFwGTE91PAk3sHCSE+AXwCoKysLAmXffioioqK00CPULl06RKapu3ylOW1a8hLlxAdHYjnn0dYrYjNTQiHwe9HBAKYX3oJmZVF+KMf3fHJpWSxvZ26ujrMZjODg4MAlJeXH6uBg8709DRCiISKYenoZQAqKyuZm5ujvb2djIwMysvLycjIMGqXNzc3J/0pQkpJb28v1dXVUZ8isrKyuHbtGl6vl6mpKe7cuUNBQQElJSWkpKTw8ssvs7y8zEc/+tGkriseHlr4oZTyy8CXAZqbm5UEKhRxoItLYWEhBQUFrK6u7k/RX16G731vxyMPhcDj2fnbYkFWVSGdTkytrcif+zlkWRnDw8NkZWUZ2ZV5eXlsbm4yPj7O8PAwpaWlR65wuLGxwczMDE1NTcd639HKAITDYXw+X8za5cdlYmICu91Obm7ugeNSUlKMnppzc3N0dXXx9a9/nddee40f/OAHSY+Tj4dkXHEaKI34vuTBz5LGWbQulJ+ueBiMj48DO7tliJ6iL370I0RLC9r0NNjtyLQ0sFggMxN5/fpOyv6DcKy5uTm2traoqqraNUdGRga1tbXU19cb1kaiIYKBQIDe3t6kFcOC3WUAUlJSMJvN9Pf3H1gG4Ci43W7m5+cTOvTVNA2Xy0VzczOrq6tkZ2fz+c9/PmlrSoRkCHkLcEkIUSmEsAK/CPzPJMxrcBYjRM7izUXxaLG4uMjCwgI1NTWGjbBXyKXPh+zu3tmBS4ksLIScHGReHtJmg7U1WFlB1taykZbG+Pj4gR1zIiscappGa2srAwMDh9YoiYxFj9Xk+TgsLi7i9/u5efPmgWUAjkIoFKK3t5eampoj7aa/853vkJ6eziuvvMKf/MmfHGstR+XY1oqUMiiE+BTwPcAEfE1K2XPslSkUL9G85QAAGixJREFUjzlut5v6+vpd4qJ3vIcH1QwDAcweD8LnA7cbbWsLabEQeuMbEaurO9bKzZt43vEOegYHuXbtWlzheiaTyegUtLCwQHd3NzabjfLy8qhZjqOjo2RkZJxILLrX62V4eNiIrolVBkCPKEmU4eFhXC7XkZKhZmZm+OM//mNu376NEOLEi2PFIilXlVK+BLyUjLl0EknLjzyEVCgeFS5cuLDPPtA0jZTBQbRXXyUcCEB9/U7suMtFOD8f0d+PCIWgooLACy9AYSHhcJjujg4uXLiQcOlYIQRFRUUUFRWxtrZm2C1lZWUUFBQghGB5eZnV1VUaGhqS+faB18veXr58eV/T58gyAOvr64yPj+Pz+Yy1xbO7Xl5eZnt7m+rq6oTXFg6H+eQnP8kf/MEfHOqrnzQimYXX46W5uVm2trbGPf6wCBEVQaKIRAjRJqVsPqXLJ+2/RL22dSS+/n5W/uqvcDY2IjUNMTmJWFqCzEzEygoyPR2Znk74F34B0dmJdv8+c4EAW297G+VPPJGUdW1vbzMxMcHa2hpFRUXMzs7S1NR0IjW3R0ZGAOL2rj0eD5OTkywvL+N0OikuLo55MOr3+2lra6OxsfFIa//qV79Kb28vf/EXf5H0CJoDiHqhxzJFX+3eFecV0+QkIZttJ4zQbIb8/J0Dzqwswo2NUFwMZWVor7yC1t/PmqYRXl2l6sc/hq2tpKwhLS2NK1eu0NDQwNTUFOFwmMnJyaTHo6+srCRc6yTeMgB6rPvFixePJOL379/na1/7Gn/4h3/4MEU8JudCyKNFiBwnOecsHp4qFPGgZWYS2t5+vUmEx0P4iScIv+ENyNRUwtXVhN7xDsTkJJ68PJY3Nii6ehXh9e7s3JPI2NgYpaWlvPGNbyQtLY2Ojg56e3sTbv0WDb/fz8DAAHV1dUc6gDysDIDeKKKgoCDhuYPBIL/5m7/Jn/7pnya9y9FRORfWymEkaq0oK+bR5lGxVvSWYcbEUhLe3mbja19jc3gYW0oK2aWlmH75lyHykC8UQnzhC4y73RRXVWE1mxFTUwR/9VfB6UzK2ubm5pifn+f69evGjlRKyfLyMhMTEwghKC8vJycn50hFpzo6OigpKTmS0MbC7XYzPj6O2+0mGAzyhje84Ujx6F/84hdZX1/nC1/4QtLWlgCPt7Wi0usV5xm9BnnYasX+8Y/j+tSn0N73PjqvX6dzbIz19XVjbFgIei9cwAHY5ufRJicJNzWBw5GUtWxtbTE2NrYrLBJ2Dkbz8/NpbGzk4sWLzMzM0NLSwuzsbMwO89EYHx8nLS0tqSIOYLfbjRDD3NxcWltbY5YBiEVPTw//43/8j1OLF4/FI7EjTzRqJZ4duYqEOb88ijvyUCgUs3Hy2toa4+PjBINBysvLWVpawmazUZmSglhehrQ0ZGXl67uYY66ptbWVmpqauOqoeL1eJicnWVpaOvTwEWB9fZ2BgQGam5tPJENyZGQETdOorKwkFAoxNzfH1NTUrjIAsfD7/fzsz/4sf/mXf3kiETpxEvUf8ZEQ8kSJR8iV/XJ+edSEPN7GyVtbW/T19eF2u7l8+fKR0+xjoSf95OXlJdzIIRgMMjMzw/T0NLm5uZSVlZGamrprTCAQoK2tjevXr59IUtHa2tqueHQdKSUrKytMTEwgpaS0tJT8/Px9n/XnP/957HY7n/vc55K+tgSI+h/AuWj1lmxUer3ivBCviMOO8IdCIZ544glmZma4c+cOLpeL4uLipCSqTE1NYTKZjtSNRz98LCkpYXFxke7ublJSUigvLyczM9OIIqmoqDgREQ8Gg/T39+80qd7zOeplAPLy8tja2mJiYoKRkRGj1ovJZKKlpYWf/OQn/OhHP0r62pLBY7kjj8XeJCSdn/5peAht9xRJ4lHakXu9XqSUh+6s9ZjoGzduGEIYDAaZnp5mZmaGgoICysrK9iXVxItueTQ1NSWljoreSm18fJxQKERaWhrhcJja2tpjzx0N/UnCGedhbyAQYGpqirm5OW7fvs1LL73EN7/5TS5fvnwi60uAx/uwMx5itX37p3861WUpHlP+6I/+iO985zuHHhSGw2Hu3btHdXX1rt2s2WymvLycJ598krS0NNrb2+nr60u4tZrf76evr+/ALvWJIoQgJyeH+vp6SktLjU5FU1NTSS2GBTv1y8PhMI4EDnv1MgBPPPEEAwMDrK2t8dWvfjWp60omSsgVijPKBz7wATo6Onj66af5/9u796iqy3SB49+XjcpEIY1cdFBAvJEpEYrOmZkgA7Q1Kmo2pS5N5HhrhSnHyanjzMpyzUQdV8c1aa4UpqgMp1XHy0zqUUdbaiWIl7yggaLc5BIKqSQh8J4/YHOwuO3Nvvt81nIt9mb/fu+zdfns335/z/u8mzZtardxVW5ubsvUQFuMXfrGjBmDj48POTk5nDp16o5Kl/YY58UHDRr0kzltS2hoaGjpLx4REUFdXR1ZWVlcvHjRpGqS9tTW1pKfn2/W3pvQtG1bdXU1ubm5vPDCC92Ox1okkbcjOlrKFYV9BQYGsnbtWvbt20dVVRXR0dGsWbOG6urqltdcuXKFurq6lja3HVFK4evry+jRowkMDOTy5cscO3aMysrKdvfKzM/Px8vLy+KlgEbnz59nwIAB3HvvvfTs2ZOQkBDGjBmDh4cHx48fJycnhxozV6S27tNiTr34d999x0svvURqaioGg8FqfweWIHPkXSAVLM7FVebIf6ympobU1FTS0tKIi4tj7NixfP/990yfPt3sKY+ampqWRTKBgYH4+/u3zMdXVlZSUFBARESEVZahl5aWUllZ2W5bXa01lZWVFBYWYjAYCAoKwtvbu8uxGBt8DRkyxOTYtNYsXryY2NhY5s6da/LxViRz5EI4M09PT5YuXcqxY8cIDg5m0aJFHDhwgPz8/G6dc/jw4YSHh3Pz5k0yMzMpLCzk5s2b5OXlMXLkSKskceMHSEdTHsZvEKNGjSIkJITi4mKOHj1KWVlZp/cNbty4QUVFhUkbRbT22WefUVNTw5w5c8w6PjExET8/P0aMGNHy3KpVqwgICCA8PJzw8HB27rRcw9i7svzQVFKuKBxJjx49UErx97//ndraWpYsWYKvry/Jyck8/PDDZiXeXr16MWTIEAYOHNiyH6VxKzhLa2xs5OzZswwfPrzLZZFeXl4t+2UWFhZy6dKldksrjRtFPPjgg2bV0VdUVLB69Wr27dtndh1+QkICSUlJPPPMM3c8n5yczO9//3uzztkRSeRdIPPiwtEkJSW1/Dx58mQOHTpESkoKtbW1JCcnEx0dbXazqdraWgYNGkTPnj05ceIEXl5eBAUFWay+Ozc3l759++Ll5WXysR4eHgwdOrSltPLo0aP06dOHwMBAPDw8AMjLyzN7o4jGxkaWLVvGq6++2q0PsqioKC5fvmz28aaSqRUhnJxSiqioKP7xj3/w5ptvsmXLFmJjY9m6davJpXylpaUtN0/NrXTpSEVFBT/88AMDBgzo/MUdaF1a6eXlxalTpzh9+jSFhYXcunWL/v37m3XeLVu24OXlxdSpU7sVX3vWrVtHWFgYiYmJVFVVWey8ksiFcBFKKcLCwvjggw/IyMjgq6++IioqinfffbdL+1revHmTgoKCO5phmVPp0p5bt26Rn5//k2Zb3eHm5kbfvn2JjIzE39+fCxcuUF9fz9WrV02Or7i4mLfeeou//vWvVrkv8Oyzz3Lx4kVOnjxJv379WL58ucXOLYlcCBc0cOBA1q1bx549eygrKyM6Opq1a9dy/fr1Nl9fX1/P2bNnGTFiRLvz1t7e3jz00EOEhoZSUVFBVlZWlzsbNjY2cubMGUJDQ80qBeyKK1euMHLkSIYPH05FRQWZmZmUlJR06VuJcdu2NWvW4O3tbZX4/P39MRgMuLm5sWDBArKysix2bknkQrgwX19fVq9ezZdffkmvXr2Ii4vj5Zdfpry8vOU1WmtycnIIDAzs0rxye5UurXun/9jFixfx8fGxWpIsKSmhV69e+Pr6tsQXERFBbW0tWVlZnbarTU1NZdiwYcTGxlolPmiatjLaunXrHRUt3SV15MLluGoduSXU1dXx4Ycfsn79ekaNGsXzzz9PUVERfn5+PPDAA2ads7OeLlevXqWgoMDsiprO1NTUcPr0aSIjI9usp29oaKC0tJTi4mJ69+5NYGDgHTv75OXlMW/ePA4fPmyxG7ozZ87k888/p7KyEn9/f1555RU+//xzTp48iVKK4OBg3nnnnS73fmlF2tiKu4Mk8s41NDSwY8cOVq1axfXr13n//fe7vfCnsbGRsrIyioqKWipdDAYDx48fN3uD466MeezYMYYNG9ZpFcyPFxgFBwfj6enJxIkTWbNmDb/85S8tHp8VyIIgIUQTg8FAfHw8np6evPbaa/z5z3/miSee4ODBgybt5tNaWz1djhw5QkBAgFWSODS1EPDx8elSKWPrBUYDBw7kiy++YMyYMfj5+REZGWmV+GxFErkQdymDwcD+/fuZMWMGu3btIiUlhfT0dCZMmMCOHTvM7kJoTJjGRl5VVVVmV7p0pLq6murqaoKDg00+tnfv3gwYMAAvLy8CAwP55JNPLBaXPcjUinA5MrXSPRcvXmTNmjVkZmaycOFCnn76aZOvqI278URERODm5tZhTxdz1NfXk52dzUMPPWRWV8YffviB8ePHk5aWRlhYmNlx2IFMrQghOjdo0CA2bNjAzp07KSgoICoqirfeeosbN2506fjbt29z7ty5O5bIm1Pp0pHz588TFBRkdmvdv/zlLzz55JPOlsTbJYlcCNGmvn378tprr3H48GGUUsTGxrJ69Wq+/fbbdo8xto5tr3+5sadLZGQkjY2NHD16lAsXLpjUe7ysrAyttTkVHwAcOXKErKysbi3Iaasp1rVr14iLi2PIkCHExcVZdOVmZySRCyE61Lt3b1asWMHRo0cZOHAgU6dOZfny5RQWFv7ktUVFRXh4eODn59fhOd3d3QkODjZ596La2louXbpEaGioWe/l5s2bLF++nNTU1G7tY5qQkMDu3bvveC4lJYWYmBjy8vKIiYkhJSXF7PObyuESuTSoEsIxeXh4sHDhQrKzs3nsscdITExk/vz5nD17Fq01169fp6ysjKFDh3b5nKbsXtTdjSK01vzpT39i4cKFZvUoby0qKoqf//zndzy3ffv2lt7lc+fOZdu2bd0awxTdutmplPovYDJQB1wE5mmtqzs+quObnc64icOqVfIB5EjkZqdtNDY2sn//ft544w3c3Nz47rvv2Lp1q1ldDVszbspcX19PUFAQffr0QSnVrY0iAP71r3/x9ttv89lnn3XrRqvR5cuXmTRpEmfOnAGaWhgYd2/SWnP//fffsZuThVjlZudeYITWOgzIBV7q5vmc0iuv2DsCIWzPzc2N2NjYlikGd3d3fve737Fz506za9Gh7Z4uly5dory83OyNIqqqqli5ciWbNm2ySBLvjFLKKqtY29Otd6S13qO1Nt52PgKY1Tty1SrZH1MIZ6W1JiEhgUOHDpGWlsaePXt49NFH+eijj7h9+7bZ5zVWuowcObKlwqW4uNjk+natNS+88AJ/+MMfzG5v2xX+/v4t/VRKS0s7vU9gSZb8aEoEdrX3S6XUQqVUtlIq+8d3vVetappOMU6pGH925EQuHz5CNDEYDMyYMQOAoUOHsnHjRnbs2EFubi6PPPIIb7/9ttkbKAMUFhYSEhLC2LFjaWxsJCsry6RKl+3bt3P79m1mzZpldgxdER8fT3p6OgDp6elMmTLFquO11ukcuVJqH9C3jV+t1Fpvb37NSmA08ITuwqS7q82RO2PMrszV5sh3797N0qVLaWhoYP78+bz44ouWHsJqqqqq2LBhAxkZGUyZMoVFixbRp0+fLh9fWVlJUVER4eHhLVMVbfV0aa/ZVXl5OfHx8ezfvx9fX1+LvCdouynW1KlTeeqppygsLCQoKIiPP/74JzdELcA6TbOUUgnAIiBGa91x7VCzjhK5M944lETuWFwpkTc0NDB06FD27t1L//79iYyMJCMjg+HDh1tyGKu7desW7733Hu+88w6PPPIIS5Ys6XSao66ujmPHjjFq1Kg7uikaGZtgFRQU0LNnT4KCgujdu3fL7xsbG5k1axbz588nPj7e4u/JTix/s1Mp9TiwAojvahLvjLMlcZDNmYX1ZGVlMXjwYEJCQujZsyczZsxg+/bt9g7LZD/72c949tlnyc7O5le/+hVz5sxh0aJFnDt3rs3+K8Ye6YMHD24ziUPnuxdt3rwZHx8fJk+ebO23Z3fd3Xx5HdAL2Nv8teeI1npxt6NyMs744SOcQ0lJyR37W/bv35/MzEw7RtQ97u7uzJw5k6effpq9e/eyYsUKPD09SU5OZsyYMS3TJ603iugKb29vvL29qamp4fLlyyQkJHDlypWWVamurluJXGs92FKBCCHuHm5ubkyYMIHx48eTlZXF66+/zrVr11i2bBlBQUGUlJQQHR1t8nk9PT0JDQ3Fzc2NiIgIkpOT+eCDD6zwDhyLw63sFEL8v4CAAIqKiloeFxcXExAQYMeILEspxdixY/n000/ZsGEDW7du5be//S0nTpwwuxZ948aNjBgxgs2bN/P+++9bOGLH1N2pFSGEFUVGRpKXl8elS5cICAhgy5YtfPTRR/YOy+KUUjzwwAOEhYUREBBAVVUVUVFRzJs3j9mzZ3d5C7ZvvvmGjIwMDh06dFdMqRhJIhfCgbm7u7Nu3TomTJhAQ0MDiYmJPPjgg/YOy2oWL16Mh4cHBoOBq1evsn79esaNG8f06dNZsGAB999/f7vH3r59m6SkJDZs2GB2e9sfCw4O5r777sNgMODu7o6j7qMgG0sIl+NK5YcCvv/+e9LS0khNTWXcuHEkJSXxi1/84ieve/3119Fa8+qrr1ps7ODgYLKzs/Hx8bHYObtJNpYQ4m4WHBzMyJEjCQ8PZ/Roe33Ome6ee+5hyZIlZGdnM3r0aGbOnMlzzz1Hbm5uS+niyZMn2bt3L3/84x/tHK19SCK3EClBFM7gwIEDnDx50mGnCDrSo0cPZs+eTWZmJtOnT2fZsmXMnj2bL7/8kueff55Nmza1W3NuLqUU48ePZ9SoUWzcuNGi57YkmSO3kFdekWQuhC24ubkxadIkJk6cyBdffMGyZcsICwuzyr2Dw4cPExAQQEVFBXFxcYSGhhIVFWXxcbrLKa/IJWEKYTpnubrsKqUUv/nNb8jOziYtLc0qYxhLPf38/Jg2bRpZWVlWGae7nDKRO0r/b+mAKJzJ4cOHOX78OLt27WL9+vUcPHjQ3iFZjDVKDWtqalo2nK6pqWHPnj137NHpSGRqpRtaN/iSxlnC0bV1demI0wSOory8nGnTpgFQX1/PrFmzePzxx+0cVduc5opcrn6FMJ8zXV06ipCQEL7++mu+/vprzp49y8qVK+0dUruc5orc0a9+pQOicGTOdHUpTOeUC4IcMZELxyELgoQLc50FQXL1K4TtJCYm4ufnd8dUzLVr14iLi2PIkCHExcVRVVVlxwiFUyZymRcXwnYSEhLYvXv3Hc+lpKQQExNDXl4eMTExpKSk2Ck6AU6ayIUQthMVFfWTvSe3b9/O3LlzAZg7dy7btm2zR2iimSRyIYTJysvL6devHwB9+/alvLzczhGZbvfu3QwbNozBgwc7/TcKSeRCiG5RSjld7++Ghgaee+45du3aRU5ODhkZGeTk5Ng7LLNJIhdCmMzf35/S0lIASktL8fPzs3NEpnGVTa2NJJELIUwWHx9Peno6AOnp6UyZMsXOEZmmrU2tS0pK7BhR99iljlwp9S1QYIOhfIBKG4wj4zvW+EFa665tvy46pZTKAB6l6d+zHHgZ2AZ8DATS9H/5Ka31NXvFaCql1JPA41rr+c2P5wBjtdZJ9o3MPHZZ2Wmr/2RKqWw7LgyR8e08vrAMrfXMdn4VY8p5lFJ/AyYBFVrrEc3PrQIWAN82v+w/tdY7zQzVFCXAgFaP+zc/55RkakUIYSvvAW31BfhvrXV48x9bJHGAo8AQpdRApVRPYAaww0ZjW5zT9FoRQjg3rfVBpVSwveMA0FrXK6WSgP8FDMDftNZn7RyW2Vw9kdu7e76ML0TnkpRSzwDZwHKttU3W+zdf/dvqG4BV2eVmpxDi7tR8Rf7PVnPk/jTdENfAaqCf1jrRbgE6KZkjF0LYjda6XGvdoLVuBDYBY+wdkzOSRC6EsBulVL9WD6cBZ+wVizNzyUSulHpcKfWNUuqCUupFO4w/QCl1QCmVo5Q6q5RaausYmuMwKKVOKKX+aYexvZVSnyilziulziml/s3WMQjH0lyP/hUwTClVrJT6d+ANpdRppdQpYByQbNcgnZTLzZErpQxALhAHFNNUZjRTa22zRgrNVxn9tNbHlVL3AceAqbaMoTmO/wBGA15a60k2HjsdOKS1Tm0u77pHa11tyxiEuFu44hX5GOCC1jpfa10HbAFsun5Ya12qtT7e/PMN4BwQYMsYlFL9gYlAqi3HbR67NxAFpAForeskiQthPa6YyAOAolaPi7FxEm2t+S79w0CmjYdeC6wAGm08LsBAmlbqvds8tZOqlPK0QxxC3BVcMZE7DKXUvcCnwDKt9XUbjmtcBn3MVmP+iDsQAWzQWj8M1AA2v1chxN3CFRO5Q/RQUEr1oCmJb9Za/4+Nh/81EK+UukzT1NJjSqkPbTh+MVCstTZ+C/mEpsQuhLACV0zkdu+hoJq67KcB57TWb9pybACt9Uta6/5a62Ca3v9+rfVsG45fBhQppYY1PxUDOG/XfiEcnMst0XeQHgq/BuYAp5VSJ5ufs1VXN0exBNjc/GGaD8yzczxCuCyXKz8UQoi7jStOrQghxF1FErkQQjg5SeRCCOHkJJELIYSTk0QuhBBOThK5EEI4OUnkQgjh5P4PNZItA/fc2zIAAAAASUVORK5CYII=\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",
    "by transporting 100% and 2/3 of the mass\n",
    "-----------------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----m = 1\n",
      "Partial Wasserstein distance (m = 1): 63.419317539744505\n",
      "Entropic partial Wasserstein distance (m = 1): 64.89236221074341\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEtCAYAAADHtl7HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfZElEQVR4nO3de7RcZZnn8d9PLsGQQIBgCNe0iDIMPaKm0VEUehQbFRc4zaC0zcB0S9CRsenBHh2mZwHToPQ0grpwbLk1SCteQJCxcRqaERAvSGChIMHh0qFJyI1LJBFIk/DMH3sfrBzPqfd9T+2q2nXy/ax1Vursd9e7n9on+6mndu16yhEhAAAA5HvZsAMAAAAYNRRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUIgCagtg+3Tbl2Sue7nts/sd06ix/UHbNw47DmDU2d7b9nrbWw07lvHIlb3bknIlBVQL2F5q+7k6qayqD8xZU5zrMNvLOpdFxKci4kM9xrh1Hd8bO5Z90HZMsOyBXrY1CLbPtP23uetHxFci4p39jAloyricMvZzYeZ9b7HdU77oJiL+KSJmRcSmfm0jB7kyD7lychRQ7fHeiJgl6fWSFkr689IJbG/deFS1iNgo6UeS3tax+G2SHphg2W39iiNXP/cFMCLeWxcqYz+nNDHpdDi2yJW/Nh3+nsNCAdUyEbFc0nclHShJtv+D7SW219l+xPbJY+uOvYKy/QnbKyVdVd93945XnbuPfwVh+5u2V9r+pe3bbP/LzPBu0+YJ4K2S/nKCZbfV2znY9o9sr7W9wvaFtretx2z7AturbT9j+17bY4/53bbvrx/zctsf74j9SNv31HP+0Pa/6hhbWu+Ln0n6Vf1K8BP1HOts/8L2220fIel0Se+v99FP6/vvaPvSOtblts8ee5vB9om2b+/YVtj+sO0H61i+YNuZ+xEYmrH/y7bPs/207X+0/a567BxVx/CFnWet6v/vH7X9oKQH62Vvtn1nnUfutP3mjm3cYvvTtn9SH9/ftr1zPbagnm/r+vedbf+N7cfreK7rEvcP6jzyS9sP2H57xzi5klw5WBHBz5B/JC2V9I769l6Sfi7pL+rf3yNpX0mWdKikZyW9vh47TNJGVQfmDEkvr5ctGzf/mZL+tuP3P5I0u77PZyXd0zF2uaSzJ4nzUElPqSq850p6VNJMSas6loWkvev13yDpTZK2lrRA0hJJp9ZjvyfpLklz6sf2LyTNr8dWSHprfXunjsf7OkmrJb1R0laSTqj33YyO/XhPvQ9fLuk1kh6TtHs9vkDSvhPtk3rZtZK+JGl7Sa+Q9BNJJ9djJ0q6vWPdkPSdOv69Ja2RdMSw/y/xw0/E5jllgrETJb0g6aT6OPqIpMcluR6/RdKHxt0nJN0kaef62NpZ0tOSjq+P7+Pq33fpmGO5qheC20u6Zux4q4/DkLR1/fvfSfp6faxvI+nQLnFvlPSn9Xrvl/RLSTvX4+RKcuVAfzgD1R7X2V4r6XZJt0r6lCRFxN9FxMNRuVXSjapeuYx5UdIZEbEhIp7L2VBEXBYR6yJig6qD47W2d8y46x2qksBv1zHcHhHPSvrHjmVLI+Kf6u3cFRE/joiNEbFU1QF3aD3XC6oS0/6qEveSiFjRMXaA7R0i4umIuLtevkjSlyLijojYFBFXSNqgKvGM+XxEPFbvi02qEt8BtreJiKUR8fBED8z2PEnvVpW0fhURqyVdIOkDXfbHuRGxtn6835N0UHoXAgNzXf2Kf+znpI6xRyPi4qiuQ7pC0nxJ8xLzfToinqqPrfdIejAirqyP76tUvUX13o71r4yI+yLiV5L+u6RjPe7CcdvzJb1L0ofrY/2FOs9NZrWkz9brfV3SL+pYyJXkyoGjgGqPoyNiTkTsExH/cewAt/0u2z+2/VRdYL1b1auXMWsi4vncjdjeyva5th+2/YyqVyIaN+eE6u38RNVp6LdJ+n49dHvHspfe07f9atvfqU+BP6OqKJxbz/V/JV0o6QuSVtu+yPYO9V1/v36cj9q+1fa/rpfvI+m0zicFVa+gdu8I87GOeB+SdKqqxLfa9tdsd67baR9Vr2pXdMz9JVWvriazsuP2s5KmdOE/0CdjOWXs5+KOsZf+79ZP7FL6/+9jHbd3V3VWpdOjkvaYZP1HVR1f4/PMXpKeioinE9seszyiOq3RMe/uErmSXDl4FFAtZnuGqlPf50maFxFzJN2g6jTumBh3t/G/j/cHko6S9A5JO6o6Vatxc3Yz9t7+W/XrpPD9jmWdF0V+UdWr0v0iYgdV76W/tJ2I+HxEvEHSAZJeLenP6uV3RsRRqg7I6yR9o77LY5LOGfekMLN+9fvStJ3BRsRXI+IQVQd9qDqF/xvr1XNvkDS3Y+4dIiL3mgdgupgsh3Quf1zVMdVpb1Vv243Za9zYC5KeGHefxyTtbHtOZmx7jLt+Zm9Jj5MryZXDQAHVbtuqOq26RtJGVxd6pj4eukrSLl1OM89W9Z//SVWnmD9VGNNtkn5XVXK8v172A1XXExykzZPCbEnPSFpve39V11pIkmz/ju032t5G0q8kPS/pRdvbuvp4744R8UJ9/xfru10s6cP1/Wx7e9vvsT17okBtv8b2v6mT6/OSnuuYa5WkBbZfJkn1KfEbJX3G9g62X2Z7X9uHTjQ3MI2tkvTKxDo3SHq17T+oL0B+v6on9+90rPOHtg+wPVPS/5B0dYxrXVAfd9+V9L9s72R7G9udF1qP9wpJH6vX+3eqrge6QeRKcuUQUEC1WESsk/QxVa8qnlb1iuj6xH0eUPUJk0fq06vjT8N+WdVp7+WqDuofF4b1Q1Wvxu4YO5UeEU+oSlyrI+LBjnU/Xse8TtUB/fWOsR3qZU/X8Twp6a/qseMlLa1PZX9Y0gfr7SxWdeHrhfX9HlJ1weJkZkg6V9Wr3pWqku9/rce+Wf/7pO2x6wb+vapEfH89/9Wqrg0BRtH/9uZ9oK7NvN/nJB3j6hNxn59ohYh4UtKRkk5Tdez+F0lH1rlgzJWqLrReKWk7VblsIserOjv1gKprnE7tEtsdkvZTdUyfI+mYiHiSXEmuHIaxT10AANAI27eo+uRWVlfvzDlPVPXpwEOamhPoBWegAAAAClFAAQAAFOItPAAAgEKcgQIAAChEAQUAAFCop29hdvVFg59T9V07l0TEud3XnxnV1+GMvvla0XV8xZb3iU5gEiueiIhdhx3FREpymLebG5q9YPLJtsvY4MYex6V0G8ecqzIm7AbU4ZmMOWYkxl9oYI6clpXrM9bp1Z6JnboyI9DUKs8mxnOsf7yBSbC5yfPXlAsoV99p9AVJh0taJulO29dHxP2T32uOqq/oGX2LdFbX8bOmyeMEenfW+K/8aIXiHDZ7gXT04skn3D9jo+P7cJeOS+lCLefLSg5LjP9DxhyvSowva2COnGeo2zPW6dV5iWrwvG3Sc6Qeyz0ZcaQK7NvPzJgEZSbPX728hXewpIci4pGI+GdJX1PV9h4ARgE5DMCU9VJA7aHNvyxymTb/IkkAaDNyGIAp6/tF5LYX2V5se3Ezb/ICwGBslr+eXzPscAC0SC8F1HJt/m3be2rzb+KWJEXERRGxMCIWVt/HCACtkMxhm+Wv7Vp5HTyAIemlgLpT0n62f8v2tpI+oMSXNwJAi5DDAEzZlD+FFxEbbZ8i6e9VfQT4soj4eWORAUAfkcMA9KKnPlARcYOkGxqKZaScpTOGHQLGOSPRWkLi74bNFeWwbVS9yTeZ+zLmSLUgmJsxR+qj7Dn9qK5OjB+YMcfSxHhOy7/UPstpyfCmBuZIPRMek2hTcEjGNlL+MGOdVJyDaOmAl9CJHAAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCop0aabUQzxS0Xf1f01YuS1ncZz8mmsxLjTWTknDlymlympBp25jT0TO2PnDmeaGCObn9XKd3gNGcbqQaoyzLmQKtwBgoAAKAQBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoNNA+UPO1Qou69Glqoo8PvYAA9IXVPWPm9FZK9QJK9SNSIgZJej5jjlSsOXOkeh818ViaeIZam7FO6rHsnxjP2V9N/N3QKpyBAgAAKEQBBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAApRQAEAABQaaCPNFZqvs7RokJuckjO6NPscQ8NOYAsT6t4IM6dx5KzEeE5GTjXjTG1DklYmxudmzJF6vDmNRXOaXPYqJ45UE8uHEuMHZmwj9XfL+dsP9BkbKZyBAgAAKEQBBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAArRVWIC06nHEz2tgAHZLWOdVO+knF5Sqaz9RMYcezYwR6rfVE6Pp5yeVSmpfZYTx3aJ8SMT4/dlbKOJvxtapacCyvZSSeskbZK0MSIWNhEUAAwCOQzAVDVxBup3I4LaGcCoIocBKMY1UAAAAIV6LaBC0o2277Ld/i+5A4DNkcMATEmvb+EdEhHLbb9C0k22H4iI2zpXqJNSnZh27HFzANCorjlss/w1e+8hhQigjXo6AxURy+t/V0u6VtLBE6xzUUQsrC7OnNnL5gCgUakctln+mrnrMEIE0FJTLqBsb2979thtSe9U3oc5AWDoyGEAetHLW3jzJF1re2yer0bE/2kkKgDoP3IYgCmbcgEVEY9Iem2DsWAKUo0yaZIJTKw4h72ovEaX3aQaNjZhTsY6GxPjTcSZ0yQz9QyU8wyVapSZ81ieT4w/lBhvoiFoTsPP1N8NA0UbAwAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCIAgoAAKBQr9+Fhz5K9XiS6PMEDFS3jJnTo6ctfXxGJY6cOJvoJZVaJ9Unqok+UE3EiYHiDBQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACg0BbZlivVoLItzSnbEgcASVbvGbMtDSxTjSFzHmfqsTQxR47tGogjtc76HmPI0UQzTgwUZ6AAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEAUUAABAIQooAACAQltkHyj6K42mUenfBUzZoHonpbQljhxN9KPq1aD2V1v2OSRxBgoAAKAYBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQKNn+y/Zlko6UtDoiDqyX7Szp65IWSFoq6diIeLp/YQLNNMqkGeeWp7EcNkPSq7qML84IJtUIcbeMOVLWZqzzQGL8TRlzLEuMz8mYY2Vi/PmMOY7pcRuStH9i/NzE+CEZ29guMZ56HDlzXJgxBxqTcwbqcklHjFv2SUk3R8R+km6ufweANrpc5DAADUsWUBFxm6Snxi0+StIV9e0rJB3dcFwA0AhyGIB+mOo1UPMiYkV9e6WkeQ3FAwCDQA4D0JOeLyKPiJAUk43bXmR7se3F0rO9bg4AGtUth22Wv9avGXBkANpsqgXUKtvzJan+d/VkK0bERRGxMCIWSjOnuDkAaFRWDtssf83adaABAmi3qRZQ10s6ob59gqRvNxMOAAwEOQxAT5IFlO2rJP1I0mtsL7P9x6o+1Hm47QclvUPpD3kCwFCQwwD0Q7IPVEQcN8nQ2xuOBSNsVPortSUODE5jOWyjuvcUyul7NAipXkGSdGBiPNWvSmqmZ9WeDcyxNDGe00sqtU6qz9OsjG2k3NLAHBgoOpEDAAAUooACAAAoRAEFAABQiAIKAACgEAUUAABAIQooAACAQhRQAAAAhSigAAAACiUbaWL6a6IJJg0qMe2FujeYzGnYmGpymZORU00ucxppLkuM5zS4TD3enOaSaxPjOQ09U/usiSaXSxPjBzWwjZy/G8/YrcIZKAAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCIAgoAAKBQq7pKpPoRSfQb6gf2KdCAV2Wsk+q/9ETGHKl+Qesz5liYGL8vY47dEuOpxyql+03lPEOtTIynek1J0pzE+NmJ8csztpF6LA9kzJHTFwsDwxkoAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUIgCCgAAoBAFFAAAQKFWNdKkoeOWiyaqaL1N6t6kMqeB5azEeKpJppRupphqCimlm1zOzZjj+cR4qkmmlH4sOY0jU89iCzLmSP3tLkmM5+yvFJpkjhzOQAEAABSigAIAAChEAQUAAFCIAgoAAKAQBRQAAEAhCigAAIBCFFAAAACFkn2gbF8m6UhJqyPiwHrZmZJOkrSmXu30iLihX0FOV6neR1tS36Mt6bFisBrNYd169eT0X0pl3FRvpRw5/YRS/ajWNjBHThw5fa9SnmhgjtR+X5AYz+kBltofTcyBgco5A3W5pCMmWH5BRBxU/1A8AWiry0UOA9CwZAEVEbdJemoAsQBA48hhAPqhl2ugTrH9M9uX2d6psYgAYDDIYQCmbKoF1Bcl7SvpIEkrJH1mshVtL7K92PZi6dkpbg4AGpWVwzbLX8+tmWgVAFuoKRVQEbEqIjZFxIuSLpZ0cJd1L4qIhRGxUJo51TgBoDG5OWyz/PXyXQcbJIBWm1IBZXt+x6/vk3RfM+EAQP+RwwD0KqeNwVWSDpM01/YySWdIOsz2QZJC0lJJJ/cxRgCYMnIYgH5IFlARcdwEiy/tQyzTSqrHk0TvI2AQGsthVve+RTk9nJIZd0BSfZ5y4kw93pw5cnofpaT6UeXEkVpnWY8x5JjbwBwYKDqRAwAAFKKAAgAAKEQBBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAAq1pa1btlSDyrY0p2xLHADwG0Yu8wPtwxkoAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUIgCCgAAoNDIdQOhv9LgpXpvSfxdsAWwumfMtRlzzOpxXJI2JsZzsvrKxPieGXOkHu+cjDnWJ8ZTj1WS5ibGc/ZHajv3JcYXNrCN1ONA63AGCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFBo5BppYvCaaJKZasZJI0603iZ1bx55UMYcy3ocl6TtEuOp5pSSdHRi/LqMOfZPjD+UMceBifHUY5WkWxLjz2fMkdjObjc90nV85Udemd5G6tn2H9JTZD0WDAxnoAAAAApRQAEAABSigAIAAChEAQUAAFCIAgoAAKAQBRQAAEAhCigAAIBCjojBbcy7h7RoYNsDhiHV80ra0vpenXVXRCwcdhS9In+h7a6JnyTX+TP9Vdfxhz+eas4l6d7uw34qo664JDE+N9306sY9Du86/ke6LDnHBs3oOr7G+0yav5JnoGzvZft7tu+3/XPbf1Iv39n2TbYfrP/dKRkpAAwQ+QtAv+S8hbdR0mkRcYCkN0n6qO0DJH1S0s0RsZ+km+vfAaBNyF8A+iJZQEXEioi4u769TtISSXtIOkrSFfVqVyj95QAAMFDkLwD9UvRdeLYXSHqdpDskzYuIFfXQSknzJrnPIr104cCOU4sSAHpE/gLQpOxP4dmeJekaSadGxDOdY1FdiT7hVWMRcVFELKwuwprZU7AAMBXkLwBNyyqgbG+jKvl8JSK+VS9eZXt+PT5f0ur+hAgAU0f+AtAPOZ/Cs6RLJS2JiPM7hq6XdEJ9+wRJ324+PACYOvIXgH7JuQbqLZKOl3Sv7XvqZadLOlfSN2z/saRHJR3bnxABYMrIXwD6IllARcTtkjzJ8NubDQdov1SjzC2rSWa7kb+wJfmeDkuu81F9oev4Ln+5LDnHU8te0XX81H0+nZxjXz2cXCflah3TdfzjOi85xyZt1XX8tC5jfJULAABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUKjoy4QB0OdpvHRfLACDsLU29TzHLls9kV5pz+7Dv9BrklPspce6jj+b8d2Tc7S26/gGbZucY7bWJdeZDGegAAAAClFAAQAAFKKAAgAAKEQBBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIVopAkMQbr55Og060zHOl1aaW4laYcu4y9kzLGxgThSaXsQ22hqO014eWI85++yTWL8ucxYerFzzzMcoB8m11l015Vdxx9Z6OQcVyTGz1Ik5/ju2f+2+wqJZp2SdOUJx3QdP/4HV6cnSfYNnfzRcgYKAACgEAUUAABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQfKKBDqj+T1EyPplHq84QxmyQ9M+wglNfXaBS20ZQmYh1En6eUVT3P8KxmJtc59Q2f7jr+zTglOcfWmtF9hZ8mp9Bxr70svVLCc4keYB95y/nJOWZoQ9fxz3YZ4wwUAABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUIgCCgAAoFCykabtvSR9WdI8SSHpooj4nO0zJZ0kaU296ukRcUO/AgUGYTo1uBxUU9A2I39hS/JDvTm5zlv1/a7j/+nWS9IbWtp9+OYT0nE8qV26jm+lTck5LtGHuo4fo6uTc+Q0H51MTifyjZJOi4i7bc+WdJftm+qxCyLivClvHQD6i/wFoC+SBVRErJC0or69zvYSSXv0OzAA6BX5C0C/FF0DZXuBpNdJuqNedIrtn9m+zPZODccGAI0hfwFoUnYBZXuWpGsknRoRz0j6oqR9JR2k6hXeZya53yLbi20vlp5tIGQAKEP+AtC0rALK9jaqks9XIuJbkhQRqyJiU0S8KOliSQdPdN+IuCgiFkbEQvVwsRYATAX5C0A/JAso25Z0qaQlEXF+x/L5Hau9T9J9zYcHAFNH/gLQLzmfwnuLpOMl3Wv7nnrZ6ZKOs32Qqo8GL5V0cl8iBICpI38B6IucT+HdLskTDNEzBRiiVJ+n6d7jKQf5C1uSx7V7cp3ZWtd9hZxzsUu7D++rh5JTbKsNXcef1NzkHGs1p+t4Ti+pvfRYcp3J0IkcAACgEAUUAABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUCinEzmAFmqiUSbNOIHp4091QXKd/6zzu47Hkon6zo6zpPuwfyfSc/x1Yny355NTfHePI7qOn6SLk3Os2zA7scb8SUc4AwUAAFCIAgoAAKAQBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUog8UsAVL9XlK9YnKmQPAYHxeH0uukzpe/dsZPZx26j788F9M3jtpzGPau+v4Ws1JznGO/lvX8bP158k5Zs9Y13X897uMcQYKAACgEAUUAABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUIhGmkALtaWBJU0ygdFx60+PSK5zzGuv7r7CJzM2tLb78JfPXpmc4mPRfZKcRpq3//3hXcdP/r0vJefYoG2T60yGM1AAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEAUUAABAIUfE4DZmr5H0aMeiuZKeGFgAUzcqcUqjEytxNq+tse4TEbsOO4heTZC/pPbu8/GIs1mjEqc0OrG2Nc5J89dAC6jf2Li9OCIWDi2ATKMSpzQ6sRJn80Yp1uliVPY5cTZrVOKURifWUYmzE2/hAQAAFKKAAgAAKDTsAuqiIW8/16jEKY1OrMTZvFGKdboYlX1OnM0alTil0Yl1VOJ8yVCvgQIAABhFwz4DBQAAMHKGVkDZPsL2L2w/ZPuTw4ojxfZS2/favsf24mHH08n2ZbZX276vY9nOtm+y/WD9707DjLGOaaI4z7S9vN6v99h+9zBjrGPay/b3bN9v++e2/6Re3qp92iXO1u3T6Yr81TvyV7PIX4M3lLfwbG8l6f9JOlzSMkl3SjouIu4feDAJtpdKWhgRretPYfttktZL+nJEHFgv+5+SnoqIc+vEvlNEfKKFcZ4paX1EnDfM2DrZni9pfkTcbXu2pLskHS3pRLVon3aJ81i1bJ9OR+SvZpC/mkX+GrxhnYE6WNJDEfFIRPyzpK9JOmpIsYysiLhN0lPjFh8l6Yr69hWq/mMO1SRxtk5ErIiIu+vb6yQtkbSHWrZPu8SJwSB/NYD81Szy1+ANq4DaQ9JjHb8vU3t3YEi60fZdthcNO5gM8yJiRX17paR5wwwm4RTbP6tPkQ/9VH0n2wskvU7SHWrxPh0Xp9TifTqNkL/6p7XH2gRae6yRvwaDi8jTDomI10t6l6SP1qdzR0JU78+29WOWX5S0r6SDJK2Q9JnhhvNrtmdJukbSqRHxTOdYm/bpBHG2dp9iaMhf/dHaY438NTjDKqCWS9qr4/c962WtExHL639XS7pW1en7NltVv8c89l7z6iHHM6GIWBURmyLiRUkXqyX71fY2qg7qr0TEt+rFrdunE8XZ1n06DZG/+qd1x9pE2nqskb8Ga1gF1J2S9rP9W7a3lfQBSdcPKZZJ2d6+vshNtreX9E5J93W/19BdL+mE+vYJkr49xFgmNXZA196nFuxX25Z0qaQlEXF+x1Cr9ulkcbZxn05T5K/+adWxNpk2Hmvkr8EbWiPN+iOKn5W0laTLIuKcoQTShe1XqnrVJklbS/pqm+K0fZWkw1R9i/UqSWdIuk7SNyTtreqb44+NiKFeADlJnIepOlUbkpZKOrnjffqhsH2IpO9LulfSi/Xi01W9P9+afdolzuPUsn06XZG/ekf+ahb5a/DoRA4AAFCIi8gBAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEAUUAABAIQooAACAQhRQAAAAhf4/scjAzWIgBKQAAAAASUVORK5CYII=\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.17456327357887044\n",
      "Entropic partial Wasserstein distance (m = 2/3): 1.070213997054379\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAEtCAYAAADHtl7HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfbklEQVR4nO3df5RkdZnf8c8HGMYBZoARmDP8HGRRISyLOAGzILBBXVQ86GaDsEogWfmRSHbJsomGkxwgqy7ZKCoH4xEWFvyFuihIWNzAGgFxZWTgEEEGMwiNzDDMAAPMAAPMwJM/7m0smu76fr9dt6pu9bxf5/Tp7vu99b1P3er79FO3bj3liBAAAADybTHsAAAAAEYNBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooDaDNg+x/ZfZ657he1P9TumUWP7I7ZvHHYcwKizvaftZ21vOexYJiJX9m5zypUUUC1ge8z2hjqprK4PzO2mOddRtld0LouIz0TEx3qMcas6vkM7ln3Edkyy7P5etjUIts+z/fXc9SPiGxHxnn7GBDRlQk4Z/7o487Y32+4pX3QTEb+OiO0i4uV+bSMHuTIPuXJqFFDt8YGI2E7SwZIWS/ovpRPY3qrxqGoRsUnSTyUd0bH4CEn3T7Ls1n7Fkauf+wIYER+oC5XxrzObmHQmHFvkyt+YCY/nsFBAtUxErJT0A0kHSJLtf217me31th+0ffr4uuPPoGx/wvZjkq6qb7trx7POXSc+g7D9t7Yfs/2M7Vtt/5PM8G7VaxPAOyX990mW3Vpv5xDbP7X9tO1Vti+2vXU9Ztuft73G9jrb99gev8/vs31ffZ9X2v7zjtiPtX13Pec/2j6wY2ys3hc/l/Rc/UzwE/Uc623/0vbRto+RdI6kD9f76P/Wt9/e9mV1rCttf2r8ZQbbp9i+rWNbYfsM28vrWL5k25n7ERia8b9l25+1/ZTth2y/tx77tKpj+OLOs1b13/vHbS+XtLxe9ru276jzyB22f7djGzfb/kvbP6uP7+/bnl+PLarn26r+fb7tv7H9aB3PtV3i/kmdR56xfb/tozvGyZXkysGKCL6G/CVpTNK76p/3kPQLSX9R//5+SftIsqQjJT0v6eB67ChJm1QdmLMlzamXrZgw/3mSvt7x+7+RNLe+zRck3d0xdoWkT00R55GS1qoqvHeS9LCkbSSt7lgWkvas13+7pHdI2krSIknLJJ1Vj/2+pDsl7VDft/0kLazHVkl6Z/3zjh33922S1kg6VNKWkk6u993sjv14d70P50h6i6RHJO1ajy+StM9k+6Redo2kr0jaVtIukn4m6fR67BRJt3WsG5Kur+PfU9Ljko4Z9t8SX3xFvDanTDJ2iqSNkk6tj6N/K+lRSa7Hb5b0sQm3CUk3SZpfH1vzJT0l6aT6+D6x/v2NHXOsVPVEcFtJ3x0/3urjMCRtVf/+d5K+XR/rsyQd2SXuTZL+Q73ehyU9I2l+PU6uJFcO9IszUO1xre2nJd0m6RZJn5GkiPi7iPhVVG6RdKOqZy7jXpF0bkS8GBEbcjYUEZdHxPqIeFHVwfE7trfPuOkSVUngt+sYbouI5yU91LFsLCJ+XW/nzoi4PSI2RcSYqgPuyHqujaoS01tVJe5lEbGqY2x/2/Mi4qmIuKtefpqkr0TEkoh4OSKulPSiqsQz7qKIeKTeFy+rSnz7254VEWMR8avJ7pjtBZLepyppPRcRayR9XtIJXfbHBRHxdH1/fyTpoPQuBAbm2voZ//jXqR1jD0fEpVFdh3SlpIWSFiTm+8uIWFsfW++XtDwivlYf31epeonqAx3rfy0i7o2I5yT9V0nHe8KF47YXSnqvpDPqY31jneemskbSF+r1vi3pl3Us5Epy5cBRQLXHByNih4jYKyL+3fgBbvu9tm+3vbYusN6n6tnLuMcj4oXcjdje0vYFtn9le52qZyKaMOek6u38TNVp6CMk/bgeuq1j2auv6dt+s+3r61Pg61QVhTvVc/0fSRdL+pKkNbYvsT2vvum/qO/nw7Zvsf3P6uV7STq785+CqmdQu3aE+UhHvA9IOktV4ltj+1u2O9fttJeqZ7WrOub+iqpnV1N5rOPn5yVN68J/oE/Gc8r416UdY6/+7db/2KX03+8jHT/vquqsSqeHJe02xfoPqzq+JuaZPSStjYinEtsetzKiOq3RMe+uErmSXDl4FFAtZnu2qlPfn5W0ICJ2kHSDqtO442LCzSb+PtEfSTpO0rskba/qVK0mzNnN+Gv779RvksKPO5Z1XhT5ZVXPSveNiHmqXkt/dTsRcVFEvF3S/pLeLOk/1svviIjjVB2Q10r6Tn2TRyR9esI/hW3qZ7+vTtsZbER8MyIOV3XQh6pT+K9br577RUk7dcw9LyJyr3kAZoqpckjn8kdVHVOd9lT1st24PSaMbZT0xITbPCJpvu0dMmPbbcL1M3tKepRcSa4cBgqodtta1WnVxyVtcnWhZ+rtoaslvbHLaea5qv74n1R1ivkzhTHdKun3VCXH++plP1F1PcFBem1SmCtpnaRnbb9V1bUWkiTb/9T2obZnSXpO0guSXrG9tau3924fERvr279S3+xSSWfUt7PtbW2/3/bcyQK1/Rbb/7xOri9I2tAx12pJi2xvIUn1KfEbJX3O9jzbW9jex/aRk80NzGCrJb0psc4Nkt5s+4/qC5A/rOqf+/Ud63zU9v62t5H03yRdHRNaF9TH3Q8k/U/bO9qeZbvzQuuJdpH0J/V6/1LV9UA3iFxJrhwCCqgWi4j1kv5E1bOKp1Q9I7oucZv7Vb3D5MH69OrE07BfVXXae6Wqg/r2wrD+UdWzsSXjp9Ij4glViWtNRCzvWPfP65jXqzqgv90xNq9e9lQdz5OS/kc9dpKksfpU9hmSPlJvZ6mqC18vrm/3gKoLFqcyW9IFqp71PqYq+f7neuxv6+9P2h6/buBfqUrE99XzX63q2hBgFP0vv7YP1DWZt/uipD909Y64iyZbISKelHSspLNVHbv/SdKxdS4Y9zVVF1o/JukNqnLZZE5SdXbqflXXOJ3VJbYlkvZVdUx/WtIfRsST5Epy5TCMv+sCAIBG2L5Z1Tu3srp6Z855iqp3Bx7e1JxALzgDBQAAUIgCCgAAoBAv4QEAABTiDBQAAEAhCigAAIBCPX0Ks6sPGvyiqs/a+euIuKDb+tvY0a1b2qrN712QwGZg1RMRsfOwo5hMSQ7baSvHotndJsvYYLfb53olMT4rY44XE+M5/xlScTQxR46XEuM5+zzVn3xTYjynr3bqapnUY5Ixx53P8z+0eVPnr2kXUK4+0+hLkt4taYWkO2xfFxH3TXWbHVR9QM9Uzu86CmA0nT/xIz9aoTSHLZotLd2/y4Rbdhkbt295nK/zXGJ8qg/g6PRQYjzn095SRUe3D/YYty5jnZRHE+N7Z8yxPDG+JjF+WMY2UgXSgxlzvNx92Ev5H9q8qfNXLy/hHSLpgYh4MCJekvQtVW3vAWAUkMMATFsvBdRueu2HRa7Qaz9IEgDajBwGYNr6fhG57dNsL7W99Pn06gDQGp356/HUdTAANiu9FFAr9dpP295dr/0kbklSRFwSEYsjYvE2PWwMABqWzGGd+Wvnnt5yA2Cm6aWAukPSvrb3tr21pBOU+PBGAGgRchiAaZv2c6qI2GT7TEn/W9X7Ty6PiF80FhkA9BE5DEAvejopHRE3SLohd/1VWkirAkzbuTq/6/j5OndAkWCmKMlhzz0vLVk69fjGjDnGluTF1c2cxPjajDl2T4w/OaA4UnPktLRakBgfuyU9x28lxlPvLLgl1QZB6fu6KD0FWoZO5AAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCIAgoAAKAQBRQAAEAhCigAAIBCfLpTi6UaR0qbV/PIzem+on1mqXtDxQ0Zc9zTUBzdvO4DSSexKDGec1+aiCPVoDInjlQTzJyGnvMS43vv0X38mkfS25ibGJ+fnoJ/2C3DGSgAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEG0lWoy+R2Xom4V+2moLaf6c3uY4uZlQuvqDjKy+cVP38cNaEkcT9mtgjo3ruo//++0z5hjAfdVzA9gGXsUZKAAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCIAgoAAKAQBRQAAEAhGmmiFZpogkmTTPTTxlek1V0aFW7ImOMnifGcXovzE+MPZMxxcGJ8RcYc8xqIY5fE+KyMOVKNMpdkzHFUYvzARKAXrklvI/W4NdG8FIPFGSgAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEH2gpinVt4ieRGXYX2i72VtJe3dr5pORTfffo4FAXkyMp5orSdKvE+PbZcyRalqVanwkSV36akmSXs6YY3X34UP3zJjjoe7D657sPv5nR2ZsI3VfH82YI7XPM/pRoTk9FVC2xyStV/VnvikiFjcRFAAMAjkMwHQ1cQbq9yLiiQbmAYBhIIcBKMY1UAAAAIV6LaBC0o2277R9WhMBAcAAkcMATEuvL+EdHhErbe8i6Sbb90fErZ0r1EmpTkzb97g5AGhU1xzWmb/25Hw9gA49pYSIWFl/XyPpGkmHTLLOJRGxuLo4c5teNgcAjUrlsM78tTMFFIAO004Jtre1PXf8Z0nvkXRvU4EBQD+RwwD0opeX8BZIusb2+DzfjIi/byQqAOg/chiAaZt2ARURD0r6nQZjGSmpxo+pRps5cwDon9Ic9uImaXmXRoUbMuZYktMsMWFeYnwsY44DE+OJ3pSSpDkNxLEgY52U/RLjdz2SnuN1155McGDi8t2LbklvI/W4pWJA+/CqPgAAQCEKKAAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCo18/Ca5229F+ixxOAzdmsAW1n44C202+bMtaZcf+wRxxnoAAAAApRQAEAABSigAIAAChEAQUAAFCIAgoAAKAQBRQAAEAhCigAAIBCFFAAAACFZlxfLhpYvlZbGosCo+5lSeu6jOc0QpzTQByppD1vANvIWWduQ9tJSe33nH2+ITG++pnu4zn7PHVfUzGgfTgDBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAApRQAEAABSacX2gcqR6I82kvkgz6b4Aw7SFuvcUyunjk9Mrqtc5mugn1EScOXHMSow30Y9qYwNxzJvdfXzDi+ltpPpiNdEjDIPFGSgAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEAUUAABAoWSfMtuXSzpW0pqIOKBeNl/StyUtkjQm6fiIeKp/YTarLc0lN6eGnsCwNJXDXpD0QI+x7JcYz2kcmWpQ+dsZc6xMjKfilNINKnPmWJ+xTsraxPjBGXOk9sdYolHmYRnbSD1uYxlzoF1yzkBdIemYCcs+KemHEbGvpB/WvwNAG10hchiAhiULqIi4Va8v8o+TdGX985WSPthwXADQCHIYgH6Y7jVQCyJiVf3zY5IWNBQPAAwCOQxAT3q+iDwiQlJMNW77NNtLbS+Vnu91cwDQqG45rDN/rRtwXADabboF1GrbCyWp/r5mqhUj4pKIWBwRi6Vtprk5AGhUVg7rzF/zBhoegLabbgF1naST659PlvT9ZsIBgIEghwHoSbKAsn2VpJ9KeovtFbb/WNIFkt5te7mkd9W/A0DrkMMA9EOy7UhEnDjF0NENx7LZaUufJ/pRYSZrKodtK+nQ3sPpWar/0qyMOd44g+JY1MAcuyfG58zuPr4h0Scqx6Lep8CA0YkcAACgEAUUAABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUCjZSBPt1kQTTBplAmmbJK3uMr4hY45lDcQxJzHeLcZx+ybGn2wgjik/ILXD3Ix1UhYlxnP2+X6p8U3dx6/L2EZqf6ViQPtwBgoAAKAQBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRB+oETcqPZya6FcFDNOcnaUDj++ywhvScxx6RAOBrEuMp5o8SdLtifE9M+Z4poE4Us2iEv2XJEn3dh8+/LCMOW5JjN/dffjUszO2kWqu9ZOMOVL740sZc6AxnIECAAAoRAEFAABQiAIKAACgEAUUAABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFBq5Rpo0ZBxNPC4YdS8+Lj3UpVHhhow5bvlc73HMT4wvy5jj0MT4iow55iXGH8iYY5fE+KyMOfZLjC/JmOOoxPiBiUAvPDq9jdTjltPvc+T+Yc9wnIECAAAoRAEFAABQiAIKAACgEAUUAABAIQooAACAQhRQAAAAhSigAAAACiXbSti+XNKxktZExAH1svMknSrp8Xq1cyLihn4F2Yl+Qu1Efy60VVM5LNS911NOH6hUL6AccxLju2XMkUr8qR5POXGkejzlbCen79HGxHjOPk89dk8+2fs2Uvsr5+8H7ZJzBuoKScdMsvzzEXFQ/TWQ4gkApuEKkcMANCxZQEXErZLWDiAWAGgcOQxAP/RyDdSZtn9u+3LbOzYWEQAMBjkMwLRNt4D6sqR9JB0kaZWkKT/hyfZptpfaXio9P83NAUCjsnJYZ/7iFBaATtMqoCJidUS8HBGvSLpU0iFd1r0kIhZHxGJpm+nGCQCNyc1hnfmriQvAAcwc0yqgbC/s+PVDku5tJhwA6D9yGIBe5bQxuErSUZJ2sr1C0rmSjrJ9kKp39o5JOr2PMQLAtJHDAPRDsoCKiBMnWXxZH2JBS6V6PEn0eUJ7NZXDtlD3vkWpPj/SYPpA5UjFMStjjpweTSk5/aZSUvdlfQNzzNuu+/iCZ9LbSD1uc9NTZD0uGBw6kQMAABSigAIAAChEAQUAAFCIAgoAAKAQBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKNdELDTPcTGqSmWoKOpPuK9on1Qhx40CiaEZbmjqm4uCfHPqFM1AAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEAUUAABAIVpkTCLVK0iiX9Co4nFDv2zKWGdZA9uZ38A2Ur2TVmTMMS8x/kDGHLskxnN6Tc1NjOfsjwWpOLZNbOOZ9DZSce6WngItwxkoAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUIgCCgAAoBAFFAAAQKFWNdJsSwNLmi0CmGjrBdLuH+2ywuz0HPse1kAgqaaN+2bMcWdifGHGHOsaiOPRjHVSEp0y3744Y47bE+P3dh/+szMytpF63JZmzJHq1vpXGXOgMZyBAgAAKEQBBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAAq1qg8U/ZcAtNWdL71dfrhLs54XMia5qIFAtkuMP5Yxx1sbmCMVx4qMOXbIWCflgMT4X2TMcXiP47+fsY3U/kptQ8r4j31exiRoSvIMlO09bP/I9n22f2H7T+vl823fZHt5/X3H/ocLAPnIXwD6JeclvE2Szo6I/SW9Q9LHbe8v6ZOSfhgR+0r6Yf07ALQJ+QtAXyQLqIhYFRF31T+vV9U4fzdJx0m6sl7tSkkf7FeQADAd5C8A/VJ0EbntRZLeJmmJpAURsaoeekzSgiluc5rtpbaXSs/3ECoATF/P+evFxwcSJ4DRkF1A2d5O0nclnRURr/kYyYgISTHZ7SLikohYHBGLpW16ChYApqOR/DV75wFECmBUZBVQtmepSj7fiIjv1YtX215Yjy+UtKY/IQLA9JG/APRDzrvwLOkyScsi4sKOoesknVz/fLKk7zcfHgBMH/kLQL/k9IE6TNJJku6xfXe97BxJF0j6ju0/lvSwpOP7EyIATBv5C0BfuHr5f0Ab864hnTaw7QFog/PvrK6BHG32W0L6Spc1NmbMsjoxviljjjmJ8XWJcUmanxjfkDHHrAbiSN2XHG9MjKf2uTTFewg67JYYvytjG6nzFakYpPQ+vzljDpSZOn/xUS4AAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEAUUAABAoZxGmgAAbaHufYty0mmqj09qPGc7OXOkPpc0px9VKo6cHk85saaktpPzuKTmmNvANlL3dVD7C03hDBQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQiAIKAACgEI00ASDHm7aVPnPo1OPPZsxxfWI8p3/lTonxezPmOCoxPpYxxw4NxLF7xjop70iM/0PGHMemxjd2H//YR9PbSO2vVAyS9IbE+Ak/yJgETeEMFAAAQCEKKAAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCIPlAAkOPBl6UT1nVZIdErSJK0rIFA5iTG16SnWLpXYoW1DcTRxByz0lNcvVtihQfSc/z9ft3Hz1yQmODm9DZS9/Xq38qYI2N/YGA4AwUAAFCIAgoAAKAQBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAolGykaXsPSV+VtEBSSLokIr5o+zxJp0p6vF71nIi4oV+BAkCpZvPXOkk3dRnflBFRTnPJlFTa3pAxx+oG5mgijiYaQy5PjK/PmCO1P+YnxjOadSaNZaxD7+s2yXk0Nkk6OyLusj1X0p22x7PI5yPis/0LDwB6Qv4C0BfJAioiVklaVf+83vYySane+QAwdOQvAP1SdA2U7UWS3iZpSb3oTNs/t3257R0bjg0AGkP+AtCk7ALK9naSvivprIhYJ+nLkvaRdJCqZ3ifm+J2p9leanup9HwDIQNAmWbyV7cPEgawuckqoGzPUpV8vhER35OkiFgdES9HxCuSLpV0yGS3jYhLImJxRCyWtmkqbgDI0lz+mje4oAG0XrKAsm1Jl0laFhEXdixf2LHahyTd23x4ADB95C8A/ZLzLrzDJJ0k6R7bd9fLzpF0ou2DVL01eEzS6X2JEACmj/wFoC9y3oV3myRPMkTPJwCt1mz+CuX1NupmY4+3z5HTjyoVRxNx5sTRhNR2cu5L6nFNjTexv3LmGMTfD3LRiRwAAKAQBRQAAEAhCigAAIBCFFAAAACFKKAAAAAKUUABAAAUooACAAAoRAEFAABQKKcTOQBgx/nS0R+devyFjDlubiCO7RLjj2XMcUBifEUDceTMsUPGOimp+7I0Y47Dexz/VMY2UvsrtQ1JekNi/OrzMiZBUzgDBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUooAAAAApRQAEAABSiDxQA5HhqrXT113ucZG3vcTw7K7HChvQc987tfY6nG4jj6dS/oNQ2JN02J7HCuvQc/zAvMT4/McFYehtPJ8avT20DbcMZKAAAgEIUUAAAAIUooAAAAApRQAEAABSigAIAAChEAQUAAFCIAgoAAKAQBRQAAEAhGmkCQJaXJK3scY5NTQTSgIwmlyMj9W8sZ5+n9keqAWoT+3MmPSabB85AAQAAFKKAAgAAKEQBBQAAUIgCCgAAoBAFFAAAQCEKKAAAgEIUUAAAAIUcEYPbmP24pIc7Fu0k6YmBBTB9oxKnNDqxEmfz2hrrXhGx87CD6NUk+Utq7z6fiDibNSpxSqMTa1vjnDJ/DbSAet3G7aURsXhoAWQalTil0YmVOJs3SrHOFKOyz4mzWaMSpzQ6sY5KnJ14CQ8AAKAQBRQAAEChYRdQlwx5+7lGJU5pdGIlzuaNUqwzxajsc+Js1qjEKY1OrKMS56uGeg0UAADAKBr2GSgAAICRM7QCyvYxtn9p+wHbnxxWHCm2x2zfY/tu20uHHU8n25fbXmP73o5l823fZHt5/X3HYcZYxzRZnOfZXlnv17ttv2+YMdYx7WH7R7bvs/0L239aL2/VPu0SZ+v26UxF/uod+atZ5K/BG8pLeLa3lPT/JL1b0gpJd0g6MSLuG3gwCbbHJC2OiNb1p7B9hKRnJX01Ig6ol/2VpLURcUGd2HeMiE+0MM7zJD0bEZ8dZmydbC+UtDAi7rI9V9Kdkj4o6RS1aJ92ifN4tWyfzkTkr2aQv5pF/hq8YZ2BOkTSAxHxYES8JOlbko4bUiwjKyJulbR2wuLjJF1Z/3ylqj/MoZoiztaJiFURcVf983pJyyTtppbt0y5xYjDIXw0gfzWL/DV4wyqgdpP0SMfvK9TeHRiSbrR9p+3Thh1MhgURsar++TFJC4YZTMKZtn9enyIf+qn6TrYXSXqbpCVq8T6dEKfU4n06g5C/+qe1x9okWnuskb8Gg4vI0w6PiIMlvVfSx+vTuSMhqtdn2/o2yy9L2kfSQZJWSfrccMP5DdvbSfqupLMiYl3nWJv26SRxtnafYmjIX/3R2mON/DU4wyqgVkrao+P33etlrRMRK+vvayRdo+r0fZutrl9jHn+tec2Q45lURKyOiJcj4hVJl6ol+9X2LFUH9Tci4nv14tbt08nibOs+nYHIX/3TumNtMm091shfgzWsAuoOSfva3tv21pJOkHTdkGKZku1t64vcZHtbSe+RdG/3Ww3ddZJOrn8+WdL3hxjLlMYP6NqH1IL9atuSLpO0LCIu7Bhq1T6dKs427tMZivzVP6061qbSxmON/DV4Q2ukWb9F8QuStpR0eUR8eiiBdGH7TaqetUnSVpK+2aY4bV8l6ShVn2K9WtK5kq6V9B1Je6r65PjjI2KoF0BOEedRqk7VhqQxSad3vE4/FLYPl/RjSfdIeqVefI6q1+dbs0+7xHmiWrZPZyryV+/IX80ifw0encgBAAAKcRE5AABAIQooAACAQhRQAAAAhSigAAAAClFAAQAAFKKAAgAAKEQBBQAAUIgCCgAAoND/B6t0ptTWdhE1AAAAAElFTkSuQmCC\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",
    "print('-----m = 1')\n",
    "m = 1\n",
    "res0, log0 = ot.partial.partial_gromov_wasserstein(C1, C2, p, q, m=m,\n",
    "                                                   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 = 1): ' + str(log0['partial_gw_dist']))\n",
    "print('Entropic partial Wasserstein distance (m = 1): ' +\n",
    "      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('Partial Wasserstein')\n",
    "pl.subplot(1, 2, 2)\n",
    "pl.imshow(res, cmap='jet')\n",
    "pl.title('Entropic partial Wasserstein')\n",
    "pl.show()\n",
    "\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
}