{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 1D optimal transport\n", "\n", "\n", "This example illustrates the computation of EMD and Sinkhorn transport plans\n", "and their visualization.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Author: Remi Flamary \n", "#\n", "# License: MIT License\n", "\n", "import numpy as np\n", "import matplotlib.pylab as pl\n", "import ot\n", "from ot.datasets import get_1D_gauss as gauss" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generate data\n", "-------------\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#%% parameters\n", "\n", "n = 100 # nb bins\n", "\n", "# bin positions\n", "x = np.arange(n, dtype=np.float64)\n", "\n", "# Gaussian distributions\n", "a = gauss(n, m=20, s=5) # m= mean, s= std\n", "b = gauss(n, m=60, s=10)\n", "\n", "# loss matrix\n", "M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))\n", "M /= M.max()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot distributions and loss matrix\n", "----------------------------------\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAADFCAYAAABzYARGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4U2X2wPHvoWW1yFZEsSCgiGylQEUERHBQdFBRB0cE\nFVdcBmVUEHQQEZfB0UH9CeIwMo47KCqDK4iAICrSYgsUZEepiiJiEdm6vL8/ToKltCVt09ws5/M8\nedIkNzcnN2nOfXdxzmGMMcZUVBWvAzDGGBMdLKEYY4wJCksoxhhjgsISijHGmKCwhGKMMSYoLKEY\nY4wJCksoxhhjgsISijHGmKCwhGKMMSYo4r0OoKjExETXrFkzr8MwxhgDpKen/+ScaxjItmGXUJo1\na0ZaWprXYRhjjAFE5OtAt7UqL2OMMUFhCcUYY0xQBJRQRORcEVkrIhtEZHQxj1cXkRm+x5eKSDPf\n/VVF5HkRWSkia0Tk7uCGb4wxJlwcsQ1FROKAycDZQDawTERmO+dWF9rsOmCnc+4kERkIPAJcBlwK\nVHfOtReRWsBqEXnVObcl2G/EGFN+ubm5ZGdns2/fPq9DMR6pUaMGSUlJVK1atdz7CKRRvguwwTm3\nCUBEpgP9gcIJpT8wzvf3TGCSiAjggKNEJB6oCRwAdpU7WhOQOXNg9WoYPhyqWKWmCUB2dja1a9em\nWbNm6L+uiSXOOXbs2EF2djbNmzcv934C+bk5Htha6Ha2775it3HO5QE5QAM0ufwGfA98AzzmnPu5\n6AuIyFARSRORtO3bt5f5TZjfPfMM/PGPcMcdMGgQ7N/vdUQmEuzbt48GDRpYMolRIkKDBg0qXEKt\n7PPXLkA+0BhoDtwpIi2KbuScm+qcS3XOpTZsGFB3Z1OEczBmDNx8M5x3Hjz4IMyYoX/n5HgdnYkE\nlkxiWzA+/0CqvL4FmhS6neS7r7htsn3VW3WAHcAg4APnXC7wo4gsAVKBTRUN3Bzqtttg0iS49lr4\n178gPh6aNtXbPXvCkiWQkOB1lMaYaBZICWUZ0FJEmotINWAgMLvINrOBIb6/BwDznS5W/w1wFoCI\nHAV0Bb4KRuDmd19/DZMnw9Ch8OyzmkwArrwSXn8dVqyAV17xNkZjjuShhx6ibdu2JCcnk5KSwtKl\nS70O6TAJvrOy7777jgEDBpS43S+//MLTTz9d6r66desGwMKFCzn//PPLFMesWbNYvfr3ZuyxY8cy\nb968Mu2jMhwxofjaRIYBc4A1wGvOuSwRGS8iF/o2mwY0EJENwB2Av2vxZCBBRLLQxPScc25FsN9E\nrHvuOb2++24oWmrt3x/atoVp00IflzGB+uyzz3jnnXdYvnw5K1asYN68eTRp0uTITzyCvLy8IER3\nuMaNGzNz5swSHy8tofhj+vTTT8v9+kUTyvjx4+nTp0+59xcsAU294px7D3ivyH1jC/29D+0iXPR5\nu4u73wRPfr4mlD59oLgp0ETguuu0kX7VKmjXLuQhmgjz179CRkZw95mSAk88UfLj33//PYmJiVSv\nXh2AxMTEg4999NFHjBgxgry8PE499VSmTJlC9erVD07TlJiYSFpaGiNGjGDhwoWMGzeOjRs3smnT\nJpo2bcpLL73EqFGj+OCDD6hSpQo33HADt956K+np6dxxxx3s3r2bxMRE/vvf/3LccccdEtfmzZsZ\nNGgQu3fvpn///gfv37JlC+effz6rVq0iKyuLa665hgMHDlBQUMAbb7zBvffey8aNG0lJSeHss8+m\nX79+3HvvvdSrV4+vvvqKdevWkZCQwO7duwHYtWsX/fr1Y8OGDfTu3Zunn36aKlWqHLLNzJkzeeed\ndxg6dCizZ8/m448/5sEHH+SNN97ggQce4Pzzz2fAgAGlHq8hQ4bw9ttvk5uby+uvv84pp5wSrI8Y\nsJHyEe+jj+CbbzRplOTKK6FqVSulmPB1zjnnsHXrVk4++WRuueUWPv74Y0B7n1199dXMmDGDlStX\nkpeXx5QpU464v9WrVzNv3jxeffVVpk6dypYtW8jIyGDFihUMHjyY3Nxcbr31VmbOnEl6ejrXXnst\nf/vb3w7bz/Dhw7n55ptZuXLlYcnG75lnnmH48OFkZGSQlpZGUlISEyZM4MQTTyQjI4NHH30UgOXL\nl/Pkk0+ybt26w/bxxRdf8NRTT7F69Wo2btzIm2++WeJ769atGxdeeCGPPvooGRkZnHjiiQcfO9Lx\nSkxMZPny5dx888089thjRzyOZRV2k0Oaspk2DerXh4suKnmbxESt+nrxRZgwAXwngcYUq7SSRGVJ\nSEggPT2dxYsXs2DBAi677DImTJhAx44dad68OSeffDIAQ4YMYfLkyfz1r38tdX8XXnghNWvWBGDe\nvHncdNNNxPsaF+vXr8+qVatYtWoVZ599NgD5+fnFJowlS5bwxhtvAHDllVcyatSow7Y5/fTTeeih\nh8jOzuaSSy6hZcuWxcbUpUuXEsd4dOnShRYttAPs5ZdfzieffFJqG01J1q5dW+rxuuSSSwDo3Llz\nqUmrvKyEEsF27IBZs7QEcqQkcf31uv3sot0pjAkTcXFx9OrVi/vvv59JkyYd/CEvSXx8PAUFBQCH\njZ846qijSn2uc462bduSkZFBRkYGK1euZO7cucVue6TutIMGDWL27NnUrFmTP/7xj8yfP7/Y7UqL\nqehr+G8Xvj8Ysxj4qxTj4uIqpX3JEkoEe+klOHCg9Oouvz59oEkTq/Yy4Wnt2rWsX7/+4O2MjAxO\nOOEEWrVqxZYtW9iwYQMAL774ImeeeSagS12kp6cDlJp8zj77bP71r38d/AH9+eefadWqFdu3b+ez\nzz4DdOqZrKysw57bvXt3pk+fDsDLL79c7P43bdpEixYtuO222+jfvz8rVqygdu3a/PrrrwG//y++\n+ILNmzdTUFDAjBkz6NGjBwCNGjVizZo1FBQU8NZbbx3cvqT9l3a8QsESSoRyTrsIn3oqtG9/5O3j\n4uCaa2DuXG1zMSac7N69myFDhtCmTRuSk5NZvXo148aNo0aNGjz33HNceumltG/fnipVqnDTTTcB\ncN999zF8+HBSU1OJi4srcd/XX389TZs2JTk5mQ4dOvDKK69QrVo1Zs6cyahRo+jQoQMpKSnF9rp6\n8sknmTx5Mu3bt+fbb4sOv1OvvfYa7dq1IyUlhVWrVnHVVVfRoEEDunfvTrt27Rg5cuQR3/+pp57K\nsGHDaN26Nc2bN+fiiy8GYMKECZx//vl069btkCq5gQMH8uijj9KxY0c2btx48P7SjlcoiA4XCR+p\nqanOFtg6spUrITkZpkyBQL8vW7ZA8+bajlJMVbCJYWvWrKF169Zeh2E8Vtz3QETSnXOpgTzfSigR\nyl9N269f4M9p1kzHpCxYUCkhGWNinCWUCLVgAZx4oraLlEXv3vDJJ5CbWzlxGWNilyWUCFRQAIsW\nQa9eZX9ur17w229gtYrGmGCzhBKBMjNh587yJRR/hw+r9jLGBJsllAi0cKFelyehJCZqrzD/Powx\nJlgsoUSgBQvgpJMgKal8z+/VS6ezP3AgqGEZY2KcJZQIk5+v7Se9e5d/H717w549sGxZ8OIypiJ2\n7NhBSkoKKSkpHHvssRx//PEHbx+opDOf5cuX88EHHwS0bY8ePcjwzZjZt2/fUgctTpw4sdRR7ddc\ncw1r164lLy+PunXrVijmt9566+BcYeHA5vKKMBkZugJjeaq7/Hr21OsFC6B796CEZUyFNGjQ4OAP\n9rhx40hISGDEiBEBPz8/P7/UwY3FWb58OatWreLcc88t0/PmzJlT6uMTJ07k2muvpUaNGsXG+Zxv\nvYnyTH1SNGb/AMhwYQklwlSk/cSvQQMdFLlwoS4bbMwhvJi/vhQXXHAB3333Hfv27eP222/n+uuv\nJy8vj8TERK6++mrmz5/Pv/71L7Zv387IkSNJSEigW7dubN26lVmzZrF7926GDRvG6tWryc3NPbh2\nyPjx49m7dy8LFy5kzJgxh0zGuGfPHoYMGcKqVato06bNISWOpKQkVq1aRVxcHH/+85/57rvvyM/P\nZ9y4cWzdupUff/yRM844g0aNGvHBBx8cFufIkSOZNGkS7XxrSdx222189NFHNG7cmOnTp9OgQQN6\n9OjBpEmTSElJYdu2bfTo0YOVK1ceFvMvv/zCqlWreOKJJ9i8eTPXXnstO3bsoFGjRjz33HMkJSVx\nxRVX0KBBA5YtW8a2bdv45z//WWmJyKq8IsyCBXDyydC4ccX207u3tqPs3x+cuIypLM8//zzp6eks\nW7aMiRMnsnPnTgBycnLo2bMnK1asoEOHDtxyyy3MnTuXtLQ0tm3bdvD548eP59xzz+WLL75g/vz5\n3HnnnYgIY8eOZfDgwWRkZBw2s++kSZOoV68ea9asYcyYMXz55ZeHxfXee+/RrFkzMjMzD85cfPvt\nt3PMMcewePHigysoFo7z9NNPP2QfOTk5dO/enaysLE4//XQeeOCBEo9DzZo1S435lltu4frrr2fF\nihVceumlh8zI/OOPP7JkyRJmzZrF3XffHeCRLzsroUSQvDxYvBgGDqz4vnr1giefhC++gDPOqPj+\nTBTxYv76Ujz++OPM9k2TnZ2dfXDhqmrVqh080169ejWtWrXihBNOAHQK+BdeeAGAuXPn8v777zNh\nwgRAZ+395ggT2i1atIi77roLgI4dO9K2bdvDtklOTmb06NGMHj2aCy64gO4l1B8XjrOo+Ph4Lr1U\n1yC84oorGDRoUKlxlWbp0qW88847AFx11VXce++9Bx+76KKLEBGSk5NLnJMsGKyEEkG+/BJ27apY\ng7xfz566mqONRzHhbN68eSxatIjPP/+czMxMkpOTD1Y/1axZ84hTy4NOVT9r1qyDU9V/8803B9cL\nqYjWrVuTlpZG27ZtGT16NA8//HCx2wUaJ/w+XX1pU/OXR/VC61tU5vyNllAiyKJFeh2M2ajr19d2\nFP8+jQlHOTk51K9fn5o1a5KVlcWyEromtmnThrVr17J161acc8yYMePgY3379uWpp546eNtffVXa\nFPM9e/bklVdeASAzM7PYqe2//fZbEhISuPLKK7nzzjtZvnz5EfdbVF5e3sGFrl555ZWD09YXnpq/\n8Nr1pe27a9euvPbaawC89NJL9PT3vgkhSygRJC0NmjaFElYiLbOuXXWfYTbhtDEH9evXjz179tCm\nTRvGjBnDaaedVux2tWrVYtKkSfTp04fU1FTq1q1LnTp1AJ3m/rfffqN9+/a0bduWcePGAXDWWWeR\nmZlJx44dD/nRBhg2bBg7duygdevWPPDAA3Ts2PGw18zMzOTUU08lJSWFhx9+mHvuuQeAoUOH0qdP\nH/r06XPE91enTh0WL15M27Zt+eSTTxjj6yUzcuRInnzySTp16nSwzehIMU+ePJmpU6eSnJzMjBkz\nePzxx4/4+sFm09dHkJYttVRxhIXsAvbss3DDDbB+vQ6UNLErGqav3717NwkJCTjnuPHGG2nfvj23\n3nqr12FFFJu+Pkbs3AkbNkBqQB9rYPz7svxtosGUKVNISUmhTZs27N27lxtuuMHrkGKO9fKKEL7q\n2aAmlLZtdS36tLTg9BwzxksjR44MaHVEU3mshBIh/KWIzp2Dt8+qVXW8mZVQDFRu7x8T/oLx+VtC\niRBpadCihfbOCqbUVEhP1zVWTOyqUaMGO3bssKQSo5xz7Nixo9jpYsrCqrwiRFoadOkS/P2mpsLk\nybBuHZxySvD3byJDUlIS2dnZbN++3etQjEdq1KhBUnmnMPexhBIBfvoJtmyBm28O/r79VWhpaZZQ\nYlnVqlVp3ry512GYCGdVXhHAN74pqA3yfq1bQ82a1o5ijKm4gBKKiJwrImtFZIOIjC7m8eoiMsP3\n+FIRaVbosWQR+UxEskRkpYhUrJIuBvkTSqdOwd93fDx07Pj7axhjTHkdMaGISBwwGTgPaANcLiJt\nimx2HbDTOXcS8DjwiO+58cBLwE3OubZALyA3aNHHiLQ0HdRYxrV4Apaaqt2S8/MrZ//GmNgQSAml\nC7DBObfJOXcAmA70L7JNf+B5398zgT+IznJ2DrDCOZcJ4Jzb4Zyzn60ySkurnOouv9RUXcHxq68q\n7zWMMdEvkIRyPLC10O1s333FbuOcywNygAbAyYATkTkislxE7iruBURkqIikiUia9TI51A8/wNat\nlZ9QwNpRjDEVU9mN8vFAD2Cw7/piEflD0Y2cc1Odc6nOudSGDRtWckiRpTIb5P1OPhkSEiyhGGMq\nJpCE8i3QpNDtJN99xW7jazepA+xASzOLnHM/Oef2AO8BldC0HL3S0nTdkmImOw2auDht8LeEYoyp\niEASyjKgpYg0F5FqwEBgdpFtZgNDfH8PAOY7HXI7B2gvIrV8ieZMYHVwQo8NaWnQqhXUrl25r5Oa\nqsuI5+VV7usYY6LXEROKr01kGJoc1gCvOeeyRGS8iFzo22wa0EBENgB3AKN9z90JTESTUgaw3Dn3\nbvDfRvTKyKic7sJFdewI+/bB2rWV/1rGmOgU0Eh559x7aHVV4fvGFvp7H3BpCc99Ce06bMro55+1\nQb5Dh8p/Lf9rZGbqLMTGGFNWNlI+jK1YodehSCinnALVqmlCMcaY8rCEEsYyMvQ6JaXyX6tqVS2Z\n+F/TGGPKyhJKGMvMhEaN9BIKHTpYCcUYU36WUMJYZmZoqrv8OnTQgZTbtoXuNY0x0cMSSpjKzYWs\nrNBUd/n5X8tKKcaY8rCEEqa++goOHAh9CQUsoRhjyscSSpjy/6iHMqHUqwdNmlhCMcaUjyWUMJWR\nAdWr6yj5UEpJsZ5expjysYQSpjIzoV07XQArlDp00NHy+/aF9nWNMZHPEkoYci70Pbz8OnTQhbay\nskL/2saYyGYJJQxt2wbbt4e2h5ef9fQyxpSXJZQw5G/D8KKE0qKFro1i7SjGmLKyhBKG/KWD5OTQ\nv3aVKtC+vZVQjDFlZwklDGVmQrNmULeuN6+fkqIxOOfN6xtjIpMllDCUkeFNdZdfhw6QkwNff+1d\nDMaYyGMJJczs3Qvr1nmfUMCqvYwxZWMJJcysWgUFBd708PJr317XsbeEYowpC0soYcbLHl5+Rx0F\nLVtaTy9jTNlYQgkzmZlQu7Y2ynvJ1kYxxpSVJZQwk5mp3YWrePzJdOgAmzbBrl3exmGMiRyWUMJI\nQYEmFC/bT/z8MfjXtTfGmCOxhBJGtmyBX3/1tv3Ez3p6GWPKyhJKGPFiDZSSHH881K9vCcUYEzhL\nKGEkM1PbTtq18zoS7TZsDfPGmLKwhBJGMjK0u26tWl5Hojp0gJUrdTp7Y4w5EksoYSRcGuT9UlJ0\n5P769V5HYoyJBJZQwsQvv2ijfDi0n/hZw7wxpiwsoYQJf/fccEoorVvrEsQ2Yt4YE4iAEoqInCsi\na0Vkg4iMLubx6iIyw/f4UhFpVuTxpiKyW0RGBCfs6BNOPbz8qlfXpGIlFGNMII6YUEQkDpgMnAe0\nAS4XkTZFNrsO2OmcOwl4HHikyOMTgfcrHm70ysyExERo3NjrSA7lXxvFGGOOJD6AbboAG5xzmwBE\nZDrQH1hdaJv+wDjf3zOBSSIizjknIhcBm4HfghZ1FPKvgSLidSSH6tABXnxR17hv2NDraExA9u2D\ntWth61a95OTomUqTJnDCCdC8efh90UxUCCShHA9sLXQ7GzitpG2cc3kikgM0EJF9wCjgbKDE6i4R\nGQoMBWjatGnAwUeLvDydtv4vf/E6ksMVbpjv08fbWEwp9u6FDz6A11+Ht9+G3btL3rZFC/jzn/WS\nkmLJxQRNZTfKjwMed86V8u0G59xU51yqcy61YQyeBq9bB/v3h1f7iZ/19Apz+/fD449DUhJccgnM\nnQuDBsGMGfD555Cdrcll3Tr46COYMkUHOz36KHTqBF27wqJFXr8LEyUCKaF8CzQpdDvJd19x22SL\nSDxQB9iBlmQGiMg/gLpAgYjsc85NqnDkUeTLL/U6nMag+DVsqLUl/hhNmHAOpk+He+7R/ubnnAN3\n3glnnaVd84pq2VIvZ50FN90EP/0Er70GDz8MZ54JF1wAjzyivTCMKadASijLgJYi0lxEqgEDgdlF\ntpkNDPH9PQCY79QZzrlmzrlmwBPAw5ZMDpeeDjVqQJuiXR3CROfOGqMJEz//DBddpCWROnVgzhy9\nnHNO8cmkOImJcMstWnL5+9/h44/1jObppzVZGVMOR0wozrk8YBgwB1gDvOacyxKR8SJyoW+zaWib\nyQbgDuCwrsWmZOnpWrUU6G9BqHXurG28v/7qdSSGzz6Djh3h/fdh4kRYvlwTSXnVqgWjR+t0CH36\naEPepZfqSFtjyiignzDn3HvAe0XuG1vo733ApUfYx7hyxBf1Cgq0OunKK72OpGSdO+tJa2Ym9Ojh\ndTQxbMoUuO027a21ZAmcemrw9n3MMdqYP3Ei3H23Jqp337UqMFMmNlLeYxs26Jl/p05eR1Iyf2xW\n7eUR52DcOK2iOu88PQMJZjLxq1IFRoyAxYu119gZZ8AXXwT/dUzUsoTiMf+PdOfO3sZRmsaN4dhj\nLaF4Ij8fhg2D+++Ha6+FN9/UdpPK1LWrloDq1NFG/A8/rNzXM1HDEorH0tN1ipO2bb2OpHTWMO+B\n/Hy46iptKL/rLnj22dA1tLVooUnlpJOgXz+YNSs0r2simiUUjy1fDsnJULWq15GUrnNn+Oor+M3m\nOwgN5+Dmm+GVV7Rr7yOPhH4A4rHHwsKF+uFfdhnMmxfa1zcRxxKKh5zThBLO1V1+nTtrBwIb4BgC\nzmmJ5N//hr/9TRvJvVK3Lrz3HrRqpV2VP/vMu1hM2LOE4qGNG3WapUhJKGDVXiHx97/DY49pF94H\nHvA6GqhXT0fgH3cc/PGPdlZhSmQJxUOR0CDv17gxNGpkCaXSvfCClkquuAL+7//CZ56tY4/VKq+E\nBE0q3xadLMMYSyieSk+HatXCv0Ee9HetUydLKJXqk0/g+uu1Z9V//qPdeMPJCSdo9deuXXDhhdag\nZg4TZt/Y2LJ8ObRvr0klEnTuDKtXw549XkcShTZtgosv1qnlZ84M314a7dvrHGIZGdoDraDA64hM\nGLGE4pFIapD38zfM+5crNkGSk6OTM+bnwzvvaJtFOOvXD/75Tx0TM2aM19GYMGIJxSObN8POnZGX\nUMCqvYKqoEDn3Vm3Dt54Q2cEjgTDh8PQodqB4PXXvY7GhAlLKB6JpAZ5v6Qknc7eEkoQTZjw+xxa\nvXt7HU3gROCpp+D00+Gaa2DNGq8jMmHAEopH0tO1mrxdO68jCZyIjZgPqrlztcpo0CCdXiXSVKum\npZOjjtL2n127vI7IeMwSikf8s5BXr+51JGVz2mm6XLFNZV9BW7bA5ZdrF7+pU8One3BZHX+8rg65\nYYOWVGwtlZhmCcUDubk6iWu3bl5HUnbdumm1/9KlXkcSwQ4c0PXc8/K0Yfuoo7yOqGJ69dKpYd58\nU5cjNjHLEooHMjJg3z7o3t3rSMqua1c9mf70U68jiWCjRsGyZfDcc5HTCH8kd9yh1V6jRtnZRgyz\nhOIB/49xJJZQjj5ahyJYQimnWbPgiSd0oaxLLvE6muARgWnTtOfGZZdpF0YTcyyheODTT3XQcePG\nXkdSPt26aRuQjWkroy1btJ2hc2f4xz+8jib46tXT9pTvvrP2lBhlCSXEnNNlJiKxdOLXrZt26MnK\n8jqSCJKbCwMHahZ+7bXI640RqC5dNFn+7386F5mJKZZQQmzrVp1XL9ITCli1V5mMGaNtC9Om6eJV\n0Wz4cJ3r6667dLliEzMsoYSY/0c4Ehvk/Vq00JmHLaEEaO5cPWu/8UYYMMDraCqfiE5uecwx2p5i\nfcxjhiWUEPv0U+0l2r6915GUn4iWUiyhBGDbNp1apV272OpS26ABvPyyLvoTiYM2TblYQgmxJUt0\ncGColgavLN266Vi2H37wOpIwVlCgM/L++qvO0FuzptcRhVbPnjB2rK7x8uKLXkdjQsASSgjt3q2L\n3UVy+4mf/z3YirCleOwx+PBD7SYcCYveVIYxY+DMM+GWW2D9eq+jMZXMEkoILVumM5RHQ0Lp1Emn\ncrJqrxIsXaorL156Kdxwg9fReCcuDl56Sb8sAwfC/v1eR2QqkSWUEPL/+J5+urdxBEONGpCaagml\nWL/8oj+exx8f2fN0BUtSks4KsHw53H2319GYSmQJJYQ+/VRrPurW9TqS4OjWDdLS7KTzEM5pb66t\nW+HVV6Pnw66oCy+EW2/Vjgnvvut1NKaSWEIJkdxcWLwYevTwOpLg6dFDk8nnn3sdSRiZNk0HLj74\nYHQURYPpH/+AlBQYMkQHY5moE1BCEZFzRWStiGwQkdHFPF5dRGb4Hl8qIs18958tIukistJ3fVZw\nw48cn3+unX3OOcfrSIKnd2/trTZnjteRhImVK/UsvE8fHdRnDlWjhvZ227dP14DJy/M6IhNkR0wo\nIhIHTAbOA9oAl4tImyKbXQfsdM6dBDwOPOK7/yfgAudce2AIELN9B+fM0fbJP/zB60iC5+ij9STc\nEgrahe/Pf9YqrpdegipW+C9Wq1bwzDOwaBHcf7/X0ZggC+Rb3wXY4Jzb5Jw7AEwH+hfZpj/wvO/v\nmcAfREScc186577z3Z8F1BSRKJ3EqHRz5+rU73XqeB1JcPXtq22t27d7HYnH/vIXWLtWB/M1auR1\nNOHtiivg2mvhoYdg3jyvozFBFEhCOR7YWuh2tu++YrdxzuUBOUCDItv8CVjunDusCVdEhopImoik\nbY/CX6afftLG62iq7vLzv6cPP/Q2Dk/99786eG/sWDgrZmt1y+app6BNGxg8GL7/3utoTJCEpFwu\nIm3RarAbi3vcOTfVOZfqnEtt2LBhKEIKqXnztPNP375eRxJ8nTrpLBsxW+21YoUO2uvVC+691+to\nIketWtp5Yfdu7WJt7SlRIZCE8i3QpNDtJN99xW4jIvFAHWCH73YS8BZwlXNuY0UDjkRz5+pSEamp\nXkcSfHFx2gY9d24MLn+RkwN/+pO2m7z6qh4ME7g2bXSczqJFcM89XkdjgiCQhLIMaCkizUWkGjAQ\nmF1km9nZ9NAgAAAQpklEQVRoozvAAGC+c86JSF3gXWC0c25JsIKOJM7p2XufPtH7e9O3r86BuHKl\n15GEkHNw9dWwebOeaR97rNcRRabBg7WE9+ijuia9iWhHTCi+NpFhwBxgDfCacy5LRMaLyIW+zaYB\nDURkA3AH4O9aPAw4CRgrIhm+yzFBfxdhLCtLF7CLxuouP387SkxVez32mC7n++ij0TW4yAsTJ+rC\nXFdfDevWeR2NqQBxYVZPkZqa6tLS0rwOI2j++U8YMQK++QaaNDny9pGqXTs9SY+JTjvz5sG558LF\nF2vpJNanVgmGb77RBrlGjXTQVu3aXkdkfEQk3TkXUIW9dZavZHPnQuvW0Z1MQEtgixfDnj1eR1LJ\nNm7U8SannKKLSFkyCY6mTTU5r12r68cUFHgdkSkHSyiVaO9ebW+M5uouv7594cABWLjQ60gq0a+/\nQv/+mkRmz7az6GA76yyd6+t//4P77vM6GlMOllAq0bvv6iwT/fp5HUnlO+MM/X194w2vI6kkBQV6\n5vzVV3omHe3rwntl2DAd9Pjgg/D6615HY8rIEkolevllbVfo3dvrSCpfzZrapDBzpibRqHPPPXrm\nPHFidM2fE25E4OmndSrrIUPgiy+8jsiUgSWUSrJzJ7z3no7ZitbuwkUNHgy7dun7jirPPAOPPAI3\n3aSTP5rKVb06vPWWno1dcIF2zTYRwRJKJXnjDW1TGDTI60hC56yz4JhjtGQWNd55R+fp6tdPpwux\nRvjQOOYYeP99XffhvPPg55+9jsgEwBJKJXnlFWjZMjpHx5ckPl5LZO++q4sWRrz0dLjsMl3DY/p0\nfYMmdFq10mrGzZu1M0RU1qVGF0soleDbb7W306BBsXdCO2iQLroV8YOeV6/WsSYNG2opJSHB64hi\n0xlnwPPPwyefaHft3FyvIzKlsIRSCaZP15k5Yqm6y69LFzjxRC2hRaxNm+Dss7VE8uGHcNxxXkcU\n2wYOhMmT4e234aqrID/f64hMCSyhVIKXX9aqrpNP9jqS0BPRRDp/vk45E3Gys7UX1759mkxatvQ6\nIgM639cjj+jZ2o032sDHMGUJJcjWrIEvv9QeT7Fq8GAtoc2Y4XUkZeRPJjt26MRk7dp5HZEp7K67\nYMwYmDZNx6tYUgk7llCC7JlntKbkssu8jsQ7rVrBqafqsYiY2olNm7S+/vvvtXdRLPWmiCTjx2ti\nmTIFrrsugr5gscESShD9+CP8+986oDrWq91HjtSJY996y+tIAvDVV9Czpw6imT8funf3OiJTEhGY\nMAHGjdOVMgcPtob6MGIJJYiefFKr3keN8joS711yibYhPfxwmC+8lZYGZ56pP0oLF1rJJBKI6Fxf\njz6q9aoXX6wrPxrPWUIJkpwcmDQJBgzQKp9YFxcHo0dre1LYrpPy1ltaMqlVS2fxbN/e64hMWYwY\nofWq77+vn+O3RReSNaFmCSVInn5aa0zuvtvrSMLH4MGQlKSllLDinC6Q9ac/QXKyrr9hZwGR6cYb\ndZzQ+vVw2mmQkeF1RDHNEkoQ7Nmjs26fey507Oh1NOGjWjVtS1m8WC9h4bffdGXAkSO1OLlggS7q\nZCLXeefpwEcRXT0zogdBRTZLKEHw7LOwfbtOSGsOdf31kJgIDz3kdSTo6PcuXeDFF7VRd/p0nSbZ\nRL4OHXRm4o4dtWh84402VYsHLKFU0Lffavtgr17a69QcqlYtLQzMmePhWinOwXPPaV/mn37SZTTv\nuw+q2Nc/qhx3nJY4R42CqVOha1cdGGZCxv6jKsA5uOEGnbtq6lSvowlft9+uy4XffLOW5EJq61ad\nKfjaazWhfPkl9OkT4iBMyMTHa7fid97RgaopKfD3v0NenteRxQRLKBXwn/9oB5MJE2yGjtJUrarz\n++XkaFIJSTfiggLN8m3bwscfa5/u+fOhceMQvLjxXL9+kJWl66ncc4+WVqzBvtJZQimnb77RM+9e\nvXQWCFO6du3g/vu12qvSp2T5+GMdT3LjjXq9ciXcdptVccWaRo10CdHXX9d/2E6dtErhhx+8jixq\n2X9YORw4oB2FnNNSiv1OBWbECO3Z+Ze/wJYtlfACWVnaFbhXL20reeUV+OgjW/891g0YAGvXwvDh\nOrr+pJO0l8iuXV5HFnXsp7CMDhzQZRkWLNAF/Jo39zqiyBEfr1VfBQXQu3cQk0p6ug7Nb9dOG9wf\nfFB/QC6/PPYWpDHFq1dP+/ZnZenSomPGwAknaOcMWw0yaCyhlIE/mfzvf5pMrr7a64giT6tWOiv8\nL79UMKkcOKB1Z717a7XWggUwdqzu8G9/s+7Apngnn6z/wMuW6Xdn/Hho2hRuukk7bJgKsYQSoD17\nDk0m1m5SfqmpMG+eJpVevWDDhgCf6JzOvTVyJDRpogsvbdmi62R8/bU20jRoUImRm6iRmqrLiq5c\nqf/YL7ygbSynnabTXlg7S7mIC7OZ+1JTU11aWprXYRxi3jwYOlSXtp40SdsATMWlp+vCiPv3wwMP\naLv5Ycu2HzgAn36q3elef10/hPh47cVz8826A2vEMhW1c6cOeJ06VavFqlTRSUMvvhj69tVunDFa\nfSoi6c65gGZNtYRSim3bdG6u//5XS8r//rfOQWeCJztbF+N7+209aZzy+D5SJV2TyOLFWpW1e7cm\nkT/8Qc8mL7oI6tf3OnQTjZzThPL661qlunat3t+smZ68dOumlxhKMEFPKCJyLvAkEAc865ybUOTx\n6sALQGdgB3CZc26L77G7geuAfOA251ypc896nVDy87Vd99//1h8553Q9n7FjoUYNz8KKPvn52pVz\n/Xrcqiy2zF7BriUraJ23kmro+hYFLU6iSt+z4ZxztCH16KM9DtrEnI0b9Qdhzhxd3iAnR+9v0EAH\nTSYn6yzVrVtrkqlfP+oSTVATiojEAeuAs4FsYBlwuXNudaFtbgGSnXM3ichA4GLn3GUi0gZ4FegC\nNAbmASc750pcZi2UCaWgQBfo27BBJ5xdskRPjHfsgIYNYcgQreqyQYtlsG+fVh/88ot23f3xR718\n/70WR7Zu1USyefOhCyMdeyy5bZLJlI68sP50ZnzTlZzqjTj1VF3vqnt3/Z9t2lQnnTQm5AoKdCqX\nzz6DpUshMxNWrYK9e3/fpm5d7frZpIlekpJ0PMwxx+ilfn3dpm5dXeMhAgQ7oZwOjHPO9fXdvhvA\nOff3QtvM8W3zmYjEA9uAhsDowtsW3q6k16tIQtm+9mcyH36X/HwOXnJz4UAuHNgP+/bDb7t1wtld\nu/T3LrfQjAzHHavJI7kDdO5UTH1+ZSrtcyj82JH+Lnxd9FJQ8Pt14UvhA5afr9NU5OXpwfNf9u/X\ny759etm7V3sq7NkDv/76+2X//uLfQ5UqOtdSUpL+o514oh7sE0+ENm30n63QW1m6VGsdliyB5ct/\nzz1VqugujjtO/zfr14c6dbRTl/9Stape4uP1fzYuTp8ncvjFr6S/jSmNFOST8MNGjt62jto/bODo\nH9aTsH0ztXZmU2vHVqrv+aXE5+bWSCC3Rm1ya9Qmr0YCedVqkV+tJvlVa5JftQb5VatTULUG+fHV\ncHFVKfBdXJU4CuLicVXicBKn11XicFIFqlTBSRWcCPivq1al69Tryv8ey5BQAvnJPB7YWuh2NnBa\nSds45/JEJAdo4Lv/8yLPPb6YgIcCQwGaNm0aSNzF+in9a/q8cFW5n8823yVcplr3ioj+GletqsWB\nqlWhevVDL7VqQe3amghq1/79Ureu9vmvV0+nGfafmSUmBpyhRXSmjK5d9fbevZpU1q/Xgs3mzdoJ\nZ/t2reLOydFt9u4N89UhTRSKA072XQ5Xi99oyHaO4Uca8QP12HnwcvS+XdTe9ytHs4sEdlOTvdTk\nV2rxA9XZT3X2U5P9VOMAVck9eB1PiRU8xdpDTahAQimLUJ6Dl8g5NxWYClpCKe9+TrqwDT99voH4\neA5eqlWLoE5ApZ0aB3I67f+78LX/UvgUvfApu/9v/6m8/3YYqVnz92qv0jinncJyc38vYOXnH1oY\nK1xgK/y84v42puKO8l2alfmZBcBe3+UQ/lqG/HykIP/368L383uthAicUNG3EaBAEsq3QJNCt5N8\n9xW3TbavyqsO2jgfyHODpmpCdRJPO7Gydm/CnMjvBShjopegJaPwa4MJ5FR0GdBSRJqLSDVgIDC7\nyDazgSG+vwcA8502zswGBopIdRFpDrQEvghO6MYYY8LJEUsovjaRYcAcNCX+xzmXJSLjgTTn3Gxg\nGvCiiGwAfkaTDr7tXgNWA3nAX0rr4WWMMSZy2cBGY4wxJSpLL6/wan01xhgTsSyhGGOMCYqwq/IS\nke3A1xXcTSLwUxDCiWR2DJQdBzsGfnYcyncMTnDONQxkw7BLKMEgImmB1vlFKzsGyo6DHQM/Ow6V\nfwysyssYY0xQWEIxxhgTFNGaUKZ6HUAYsGOg7DjYMfCz41DJxyAq21CMMcaEXrSWUIwxxoSYJRRj\njDFBEVUJRUTOFZG1IrJBREZ7HU+oiEgTEVkgIqtFJEtEhvvury8iH4rIet91Pa9jrWwiEiciX4rI\nO77bzUVkqe87McM3wWlUE5G6IjJTRL4SkTUicnqsfRdE5Hbf/8IqEXlVRGrEwndBRP4jIj+KyKpC\n9xX72Yv6P9/xWCEinSr6+lGTUHxLFU8GzgPaAJf7liCOBXnAnc65NkBX4C++9z4a+Mg51xL4yHc7\n2g0H1hS6/QjwuHPuJGAnEJqVhrz1JPCBc+4UoAN6PGLmuyAixwO3AanOuXbopLYDiY3vwn+Bc4vc\nV9Jnfx46A3xLdIHDKRV98ahJKOi69Rucc5uccweA6UB/j2MKCefc98655b6/f0V/QI5H3//zvs2e\nBy7yJsLQEJEkoB/wrO+2AGcBM32bxMIxqAP0RGcAxzl3wDn3CzH2XUBnUq/pW5+pFvA9MfBdcM4t\nQmd8L6ykz74/8IJTnwN1ReS4irx+NCWU4pYqPmy54WgnIs2AjsBSoJFz7nvfQ9uARh6FFSpPAHeh\ni92BLkP9i3Muz3c7Fr4TzYHtwHO+qr9nReQoYui74Jz7FngM+AZNJDlAOrH3XfAr6bMP+m9mNCWU\nmCciCcAbwF+dc7sKP+Zb8Cxq+4iLyPnAj865dK9j8Vg80AmY4pzrCPxGkeqtGPgu1EPPvpsDjdE1\neItWA8Wkyv7soymhhHS54XAjIlXRZPKyc+5N390/+IuwvusfvYovBLoDF4rIFrS68yy0LaGur9oD\nYuM7kQ1kO+eW+m7PRBNMLH0X+gCbnXPbnXO5wJvo9yPWvgt+JX32Qf/NjKaEEshSxVHJ11YwDVjj\nnJtY6KHCSzMPAf4X6thCxTl3t3MuyTnXDP3s5zvnBgML0GWpIcqPAYBzbhuwVURa+e76A7piasx8\nF9Cqrq4iUsv3v+E/BjH1XSikpM9+NnCVr7dXVyCnUNVYuUTVSHkR+SNaj+5fqvghj0MKCRHpASwG\nVvJ7+8E9aDvKa0BTdEmAPzvnijbYRR0R6QWMcM6dLyIt0BJLfeBL4Arn3H4v46tsIpKCdkyoBmwC\nrkFPHmPmuyAi9wOXoT0gvwSuR9sHovq7ICKvAr3Qaep/AO4DZlHMZ+9LtpPQ6sA9wDXOuQotlxtV\nCcUYY4x3oqnKyxhjjIcsoRhjjAkKSyjGGGOCwhKKMcaYoLCEYowxJigsoRhjjAkKSyjGGGOC4v8B\n1ZcTbJpyYBQAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYHWWZ9/Hv3d1ZCAkkIWELS7OKDGuMEAQZJGwimBnh\nRRQlw5ZB2SagAoZhYBzBeZXNKEhYFDQgGFYB8VVkUwmSAMomISYhgJCEfYtZup/3j7tOztadPt2n\nllN9fp/rqut0naquejih7/511VPPYyEEREQkeS1ZN0BEpFmo4IqIpEQFV0QkJSq4IiIpUcEVEUmJ\nCq6ISEpUcEWahJl90syez7odzUwFV5qemX3RzGab2ftm9qqZ/crM9qrzmAvNbL+42ljD+YKZbb2m\nfUIID4cQPtLH4y80sxVmNqri/Seic7f35bjNRgVXmpqZnQ5cClwAbABsBlwOTMyyXXEzs7YYDrMA\n+ELJMXcEhsRw3KahgitNy8zWBf4bOCmEcGsI4YMQwsoQwi9DCF+P9hlkZpea2d+j5VIzGxRtG2Vm\nd5nZ22b2ppk9bGYtZvZTvHD/MkrN3+ji3PuY2ctm9g0zWxIl638xs4PNbG50vG+W7L+bmT0SnetV\nM/uBmQ2Mtj0U7fbn6HyfLzn+mWb2GvDjwnvR92wVnWNstL6xmS01s33W8JH9FDi6ZH0ScH2fPvwm\npYIrzWwPYDBw2xr2mQqMB3YBdgZ2A86Jtp0BvAyMxtPxN4EQQvgysAg4NIQwNITwf7s59obR+ccA\n5wJXAV8CPgZ8EvhPM9si2rcDmAKMito9AfgqfsK9o312js53U8nxRwKbA5NLTxxC+BtwJvAzMxsC\n/Bi4LoTwwBo+i1nAOmb2UTNrBY4EfraG/aWCCq40s/WA10MIq9awz1HAf4cQloQQlgLnA1+Otq0E\nNgI2j5Lxw6F3g5OsBL4dQlgJ/BwvppeFEN4LITwDPIsXeUIIc0IIs0IIq0IIC4ErgX/u4fidwH+F\nEJaHEJZVbgwhXAXMAx6N/jum1tDmQsrdH3gOeKWG75GICq40szeAUT1c39wYeLFk/cXoPYDv4gXr\n/5nZfDM7q7fnDyF0RF8XCuLiku3LgKEAZrZtdPniNTN7F7/mXHYDqwtLQwj/6GGfq4AdgGkhhOU1\ntPmnwBeBf0OXE3pNBVea2SPAcuBf1rDP3/E/yQs2i94jSqJnhBC2BD4LnG5mE6L94h6G7wrgr8A2\nIYR18MsX1sP3rLENZjYUv2F4DXCemY3sqREhhBfxm2cHA7fW0G4poYIrTSuE8A5+7fSH0Q2rIWY2\nwMw+bWaF6643AueY2eioS9S5RNctzewQM9vazAx4B7/O2hl932JgyxibOwx4F3jfzLYDvlKxvS/n\nuwyYHUI4Hrgb+FGN33ccsG8I4YNenq/pqeBKUwshXAScjt8IWwq8BJwM3B7t8j/AbOAvwFPA49F7\nANsAvwXex9Py5SGE+6NtF+KF+m0z+1oMTf0a/qf8e/hlgJsqtp8HXBed74ieDmZmE4GDKBbu04Gx\nZnZUT98bQvhbCGF2L9ouEdMA5CIi6VDCFRFJiQquiEhKVHBFRFKigisikpI4BrSQfmDUqFGhvb09\n62aI5MqcOXNeDyGMrnV/FVwBoL29ndmz1dNHpDfM7MWe9yrSJQURkZSo4Io0i2XL4LHHYN48UP/7\nTKjgivR3TzwB++4L66wDu+0G22wDI0fC6afDe+9l3bqmooIr0l91dsJZZ8G4cfDMM/D1r8Mtt8BV\nV8GnPw2XXgrbbQf339/zsSQWumkm0h91dsJXvgLTp8Nxx8F3vwsjRhS3H388nHYaHHMMHHww3HUX\nTJjQ/fEkFkq4Iv3RKad4sf3mNz3Rlhbbgt13hwcf9EsMhxwCDz1UvY/ESgVXpL+5/nq4/HI44wz4\nn/8BW8OwuaNHw333wWabwec/D0uWpNfOJqSCK9KfPP88fPWrsPfe8J3vrLnYFoweDb/4Bbz1Fkya\n5JcjJBEquCL9RUcHHHUUDB4MM2ZAWy9u0ey0E1xyCdx7L/zwh8m1scmp4Ir0F9Onw5w5XjA32aT3\n33/iiXDAAXDOObB4cc/7S6+p4Ir0B6+/DlOnen/bI3qc8KFrZjBtmj8gcVZv58OUWqjgivQHU6f6\nQwzTptV23bY7227rN9t+8hN45JHYmidOBVck7+bOhauv9ptl229f//GmToUNN4Szz9YjwDFTwRXJ\nu/PP9xtlU6fGc7yhQ73/7oMPwu9+F88xBVDBFcm3Z56BG2/0Bx3WXz++455wgt94O+ccpdwYqeCK\n5Nn553si/frX4z3u4MFebGfN8q5iEgsVXJG8euEFmDkTTj4Z1lsv/uMfcwxsuqk/QCGxUMEVyauL\nLoKBA+HUU5M5/sCBMGWKj7Hw6KPJnKPJqOCK5NHixd51a9Ik71GQlOOPh+HDfbQxqZsKrkgeTZsG\nK1Z4n9kkDRvm3c1uvdUvYUhdVHBF8mbZMrjiCpg40R9USNopp8CAAfD97yd/rn5OBVckb264Ad58\n0wcQT8OGG8KRR/oljHffTeec/ZQKrkiehOCXE3bcEf75n9M77ymnwPvve9GVPlPBFcmThx+GP//Z\nC2A9Yyb01rhxsMceXuw1Xm6fqeCK5Mm0aT5dzlFHpX/uU0/1Kdb1IESfqeCK5MXf/w633QbHHgtD\nhqR//sMO8+u5V1yR/rn7CRVckby45hqf1eHf/z2b8w8Y4DMA3303vPhiNm3IORVckTxYtcpndNh/\nf59lNyuTJ/u146uuyq4NOaaCK5IH99wDL78MX/lKtu3YbDP4zGd8/N0VK7JtSw6p4IrkwY9+BBtv\nDIccknVLfO6zxYvh9tuzbknuqOCKNLqFC71nwPHH+3XUrB14IGy+uV/ikF5RwRVpdFdf7ddNjzsu\n65a41lYv/vfd593EpGYquCKNbOVKuPZa+PSn/fppozj2WC+8unnWKyq4Io3s7rvh1Vez6wrWnY03\nhkMPhR//WDfPekEFV6SRXXkljBnjCbfRTJ4MS5fq5lkvqOCKNKqFC+HXv/brpW1tWbem2gEH6OZZ\nL6ngijSqws2yY4/NuiVd082zXlPBFWlEjXqzrJJunvWKCq5II7rrrsa8WVZJN896RQVXpBE18s2y\nSrp5VjMVXJFGM3++3yw74YTGvFlW6YADoL1dwzbWQAVXpNFceWXxhlQetLb6pY8HHoC//jXr1jQ0\nFVyRRrJ8ud8s++xn/ZJCXhx7rI/z8KMfZd2ShqaCK9JIbr0VXn89+2EYe2v99X1GiOuugw8/zLo1\nDUsFV6SRXH45bLUVTJiQdUt678QT4e234cYbs25Jw1LBFWkUTz4Jv/89nHQStOTwR3PvvWGHHeAH\nP/Dp3KVKDv9VRfqpH/zAJ4c85pisW9I3ZnDyyf6L449/zLo1DUkFV6QRvPEGzJgBX/4yDB+edWv6\n7ktfgnXX9V8eUkUFV6QRXHst/OMffjkhz9Ze23sszJzpT8pJGRVckaytWgXTpsE++8COO2bdmvqd\ndJJP53755Vm3pOGo4IpkbeZMeOklOP30rFsSj622gokT/ckzdREro4IrkqUQ4KKLYJttfPrx/uKM\nM/y69PXXZ92ShqKCK5KlP/wBZs+GKVPy2RWsO3vuCR//OFxyCXR2Zt2ahtGP/oVFcui734WRI+Ho\no7NuSbzM/BLJ3Llw551Zt6ZhqOCKZOWpp7wYnXqq393vbw4/HLbcEi68UA9CRFRwRbJy4YUwdCic\nckrWLUlGWxuceSb86U8+DY+o4IpkYt48uOkmH6Rm5MisW5OcSZN8VogLLsi6JQ1BBVckCxde6MMZ\nTpmSdUuSNWgQfO1rcP/9Pk5Ek1PBFUnb3Lk+jOGJJ8JGG2XdmuRNngwbbADnnNP013JVcEXSdt55\nnvzOPjvrlqRj7bVh6lR48MGmv5argiuSpqeegp//HE47zVNfs5g82ad7nzq1qVOuCq5Ims48E9ZZ\nx69rNpNBg+Dcc73Hwi23ZN2azKjgiqTlV7/y5T//s3/3TOjOpEk+OM/Xv+4jozUhFVyRNKxc6U9e\nbb11/+1325O2Nrj0Uli40B/5bUIquCJp+OEPfQrxiy6CgQOzbk129t3XRxL79rfhlVeybk3qVHBF\nkrZokXeJOuggOPTQrFuTvYsv9gFtTjqp6W6gqeCKJCmE4pTnP/qRD+rS7LbcEs4/H+64w6eFbyIq\nuCJJuvFGuOce/xN6882zbk3jmDIFdt3VJ518882sW5MaFVyRpCxY4Ol2jz28sEhRW5vP4/bGG3D8\n8U1zaUEFVyQJq1bBUUf51zfcAK2t2banEe2yi48pcdttcOWVWbcmFSq4Ikk4+2x45BEvJO3tWbem\ncU2ZAgce6K9z5mTdmsSp4IrE7brr4Hvf88sJRx6ZdWsaW0uLz3u2/vreXayfT62ugisSp4cf9nED\nJkyAyy7LujX5sP76PvPF22970X3//axblBgVXJG4/OlPPvPuFlvAzTf7eLdSm513hhkz/LLCZz8L\ny5Zl3aJEqOCKxGH2bL8WOXq0D0HYjGMl1GviRL8c88AD8LnPwYcfZt2i2KngitTrV7+CffaBddf1\nYjtmTNYtyq8vfQmuugp+/Wu/LPP661m3KFYquCJ9FYIPxnLoobDttt4rQT0S6nfccTBzJjz5JIwf\n76/9hAquSF8sWeJ/Ak+ZAocc4rMZNMN0OWn53Ofgd7/za7njx8O0aT7+Qs6p4Ir0RkeH/8m73XZw\n773eE+G222DYsKxb1v/ssYen2333hVNPhb32gj//OetW1UUFV6QWHR3e82Cnnbzb1047+Q//qadq\nQJokjR4Nd98NP/kJvPCCj7/wxS/Cc89l3bI+UcEVWZMFC3zgmS23hM9/3t+7+Waf9vujH822bc3C\nzGeLeP55n6Lozjth++1h//393+KDD7JuYc0sNMmgEbJm48aNC7Nnz866Gdl7913vT3v//X7J4PHH\n/f0JE+CrX/XrthoXIVtLl8L06b4sWgRrrQUHHAD77Qd77+3FuK0tlaaY2ZwQwria91fBFWiigrt8\nuQ8HuHgxvPaa/8AuWODp6ZlnYO5c36+tDXbbDf71X+Gww/xhBmksHR3+ZN/MmT4E5oIF/v5aa/kl\nn49+FLbZxofF3HRTnyV59GjvvhfTL00VXOmTsoJ7++1w113Jn7S7//cK75e+drd0dPjd61Wr/OuV\nK31ZvtyXZct8ee89X7p6gqmtzeca++hHYexY+PjH4ROf0I2wvFmwAP74R3jsMZ+O/rnnuh6bwcz/\nbYcNg7XXhiFDYPBgn1l44EB/QrCtzZfWVh/vYe+9/S+cqkP1ruCmk7slX+bN8z+n09DdDafC+6Wv\nXS2FH4jCD8eAAb4MGgRDh3raWWut4g/YiBH+FNj663vi2Wwz2HhjXSboD7bYwpfCsJjg13cXLfL5\n0xYv9gcp3nwT3nnHfwF/8IE/0faPf/gv6Pfe81/Yq1b50tnpS0yDxyvhCtBElxREYtTbhKteCiIi\nKVHBFRFJiS4pCABmthR4seStUUAjjxzSyO1r5LZBY7evkdsG1e3bPIQwutZvVsGVLpnZ7N5cm0pb\nI7evkdsGjd2+Rm4b1N8+XVIQEUmJCq6ISEpUcKU707NuQA8auX2N3DZo7PY1ctugzvbpGq6ISEqU\ncEVEUqKCKyKSEhVcqWJmB5nZ82Y2z8zOyrgtm5rZ/Wb2rJk9Y2anRe+PNLPfmNkL0euIDNvYamZP\nmNld0foWZvZo9PndZGYDM2zbcDObaWZ/NbPnzGyPBvvspkT/rk+b2Y1mNjjLz8/MrjWzJWb2dMl7\nXX5e5r4ftfMvZja2p+Or4EoZM2sFfgh8Gtge+IKZbZ9hk1YBZ4QQtgfGAydF7TkLuC+EsA1wX7Se\nldOA0ikI/he4JISwNfAWcFwmrXKXAfeGELYDdsbb2RCfnZmNAU4FxoUQdgBagSPJ9vP7CXBQxXvd\nfV6fBraJlsnAFT0ePYSgRcvqBdgD+HXJ+tnA2Vm3q6Q9dwD7A88DG0XvbQQ8n1F7Nol+CPcF7gIM\nfxKpravPM+W2rQssILo5XvJ+o3x2Y4CXgJH4yIV3AQdm/fkB7cDTPX1ewJXAF7rar7tFCVcqFX4I\nCl6O3sucmbUDuwKPAhuEEAqDnb4GbJBRsy4FvgEUppRdD3g7hLAqWs/y89sCWAr8OLrkcbWZrU2D\nfHYhhFeA7wGLgFeBd4A5NM7nV9Dd59XrnxUVXMkFMxsK3AL8Rwjh3dJtweNF6v0bzewQYEkIYU7a\n565RGzAWuCKEsCvwARWXD7L67ACia6ET8V8MGwNrU/3nfEOp9/NSwZVKrwCblqxvEr2XGTMbgBfb\nGSGEW6O3F5vZRtH2jYAlGTRtT+CzZrYQ+Dl+WeEyYLiZFQb3z/Lzexl4OYTwaLQ+Ey/AjfDZAewH\nLAghLA0hrARuxT/TRvn8Crr7vHr9s6KCK5UeA7aJ7hQPxG9i3JlVY8zMgGuA50IIF5dsuhOYFH09\nCb+2m6oQwtkhhE1CCO345/S7EMJRwP3A4Vm2LWrfa8BLZvaR6K0JwLM0wGcXWQSMN7Mh0b9zoX0N\n8fmV6O7zuhM4OuqtMB54p+TSQ9eyuFiupbEX4GBgLvA3YGrGbdkL/xPuL8CT0XIwfq30PuAF4LfA\nyIzbuQ9wV/T1lsCfgHnAL4BBGbZrF2B29PndDoxopM8OOB/4K/A08FNgUJafH3Ajfj15Jf4XwnHd\nfV74DdIfRj8nT+G9LdZ4/Loe7TWzg/A/oVqBq0MI3+nzwURE+rk+F9yov+ZcvIvOy/ifol8IITwb\nX/NERPqPembt3Q2YF0KYD2BmP8fvOHZbcEeNGhXa29vrOKXU6+23fWLSTTctf3/OnDmvh2jk+v1b\n/k9ffwtXrLdUrHazveJ9KxynpaX8uNF6cXthRt+K47S0rv569b6FWXkrZwNuLXyvv4aWinO3lrch\nrH6/u/XoNfq+0Gpdb49eO9sqv6/wSvl6dJrOiu2V64X9itvLj9PZVr696rVwnsr92kI3xw3l26NX\novdpC1j0tbV5z7nW1ug1Wh8wwHuADWjtAGBgm78Oai28+va12lYCMDh6Xbt1hb/fGq23LQdgSIu/\nP6z1HwAMjV7XaVkWbV8erfv7w1a/+nGGWYjW/UMY2jIYgJYNX+hmiuna1VNwu+qDtnvlTmY2GX8K\ng8022wzNDJutb3wDpk2Dyn8GM3ux6+8QkbjUU3BrEkKYTjSG5Lhx4zQWZMZWroQBAxI6eOHyVCE9\nhuhZgCiBhs4o6bRUbO8sT6iFy1zWGW0vpMhovZAqrfCoQUvFcejAbysUE5t1eFpanXQLOjrLVi3q\nuBMof7+QdAttKjwzZFSuF6xuXLSdiu3R1qh7f2fVT2Jhz/LvbInWO7tZr1T4RDqj/Vqi/Tq73LuL\ndnXTnuJxuz5v6PJr/66Oqr1jLkP1Hi5KunQWnr2IEnKdh4X6uoU1XH9N6dmKFTAws6FURJpbPb8L\nVvfXxAvtkcAXY2mVJObDD2HttRM+SUMk3UKOUtItnl1Jt1cqkm4cCbfPTQohrDKzk4Ff4/9XXxtC\neCaGNkmC3n0Xhg7NuhUizamu3wEhhHuAe2Jqi6TgjTdgvfVSOlmWSbfsei4o6ZaePd2kW3njJrdJ\nNwZ6tLfJvPIKbLhh1q0QaU6J91KQxtHRAQsXwuGH97hrvLJIul32XAAl3dKzK+mmTQm3iTz9NKxa\nBTvskHVLRJpTA9R8ScuDD/rr3ntn1IAUk+6a++iCkm7p2ZV006KE20R+9jPYfvvqx3pFJB1KuE1i\n1ix47DF/rDdzKSTd2p5GAyXd0rMnk3Rr6aNbvt5/k64SbhNYsQK+8hXYYAM4+uisWyPSvJRw+7kQ\nYOpUePJJuP12WKeWx2VWJ8+Eh75IMun2atwFUNItPXu8Sbc3T6OVr/e/pKuE24+F4KODfe97nnAn\nTsy6RSLNTQm3n3rnHTj1VLj+ejjpJPj+9/twkDwn3T6NMAZKuqVnjyvp9n7chfL1/pN0lXD7odtv\n994IP/sZnHuu3yhr0b+0SOaUcPuRRx6BCy+EX/4SdtrJC+/HP96HAxVmUFidNPOXdOsbSxeUdEvP\nXm/S7fsIY+Xr+U+6yj0519EBt90Ge+4Jn/gE/P73cMEFPqNDn4qtiCRGCTennnsObroJZsyAefOg\nvd2v0x5zTIzDL+Y46cYzawQo6ZaevW9JN46xdMvX85t0VXBzZO5cL7I33+zjIpj5Y7rf/jZ87nPQ\npn9NkYamH9EG9v778PDD8Nvfwm9+A0895UV2r738Rthhh8FGG8V/3kJyLCTJXCbdWOdHAyXd0rP3\nLunGOWtE+Xr+kq4KbgNZuRIefdQL7H33+eO4q1b5HGR77gmXXupDK44Zk3VLRaQvVHAztGSJF9jC\n8sgj8MEHHug+9jE44wzYbz8vtmutlX77cp10E5kJGJR0S8/eZEk3Bo3Xon5q+XJ/vHbWLC+us2bB\nggW+rbUVdtwRJk2CCRNgn31g5MhMmysiCVDBTcD778Nf/uIF9okn/PUvf/FBZMAvCYwf74/b7r67\np9nEZ9LtjYokm8ukm8hMwKCkW3H81a1R0q1F47QkpxYvLi+sTzwBL7xQrA0jR8Kuu8Jpp3lx3X13\n2GSTbNssItlQwa3RsmXw7LPeU6B0ee214j7t7V5cjzoKdtnFv95kk2Igy43CtdDO/CbdRGYCLj2O\nkm758Ve3Rkl3TbJvQYPp6ID586sL67x5q38WGTzYxyo48MBiYd15Zxg+PNu2i0hja+qCu2RJdWF9\n5hn48EPfbgZbbeU3tI480l933BG23ro6wPRLOU66icwEDEq6fUi69YylW37M/Cfdpii4H37ohbSy\nuC5ZUtxn9GgvpiecUCys//RPDXYzS0Ryrd8V3Dff9BtXjz9efJ07txiGBg/2acI/85liYd1xR59+\nRpxZefrLZdJNYCZg3x6dU0m3bL1SadJNYybg0hZXrzdO0s1twQ0B/v736uK6aFFxn003hbFjyy8H\nbLVVk1wOEJGGk5uC29nplwUeesjHF3joIXj11eL2bbeFPfbw2Q123dWXUaOya2+utVSkqjwm3QRm\nAgYl3b4k3SRmAl7TORs56TZswV25EubMKRbXP/wB3nrLt40Z409jjR/vCXbnnWHYsEybKyLSo4Yq\nuCH4ANrXX+9DEL77rr+/7bY+/OAnP+nDEba357Bva54UPtw8J90kZgIGJd0+JN2+9NEta1M3bclj\n0u3xDGa2KXA9sAHe5ukhhMvMbCRwE9AOLASOCCG81ZdGzJ8PP/2pF9r5871nwGGHwaGH+lCEG27Y\nl6OKiDSWWkr6KuCMEMLjZjYMmGNmvwH+DbgvhPAdMzsLOAs4s7cNuOgi+NrXPJBMmADnnedpVt2x\nMlSR5nKZdJOYCRiUdJV069LjkUMIrwKvRl+/Z2bPAWOAicA+0W7XAQ/Qy4I7fboX28MOg0su8V4F\nIiL9Va9KuZm1A7sCjwIbRMUY4DX8kkPNFiyAE0+EcePghht8kG1pDKv74SrpKukq6caqpeddnJkN\nBW4B/iOE8G7ptuD/N3b5U2Nmk81stpnNXrp06er3N9kEPvUp7z97zz19a7yISJ7UVMLNbABebGeE\nEG6N3l5sZhuFEF41s42AJV19bwhhOjAdYNy4cauL8oABcPvtsP/+fknhwAN9AO6JE/1pMMlQlOoK\naU1JFyXdOpJuEjMBl++Xn6TbY8I1//vyGuC5EMLFJZvuBCZFX08C7ujtyYcNg3vvhTPP9LENjjzS\neyRMnuz9bjt7+uRFRHKkltK9J/Bl4CkzezJ675vAd4Cbzew44EXgiL40YPhwuOAC+Na34IEH4Lrr\nYMYMuOoqGDHCu4Xtvbf3wR071pOxJKyQQKN0pqSLkm4dSTeJsXTL2tRNW+JPuvWrpZfC76n+/Asm\nxNWQ1lbvFjZhAlx+OdxxB9x/vz9p9stf+j5Dhvjju4UHIHbbTd3HRCQ/GupJs4KhQ33WhKOO8vXX\nXvMn0ArjKJx/fvFhoo98xMdNGDu2OIaCJmCsU2VCVdJV0i1s7UPSTXLWiLI2ddOWRkq6DVlwK224\nIRx+uC8Ab78Nf/wjPPaY93L4/e/hxhuL+2++ebEAjx3rszJsvLEeBxaRbOWi4FYaPhwOPtiXgtdf\n9+JbWB5/3HtBFILMiBHFIRp32slfd9hBg950qaWQqqLf7Uq6SrrlR+1V0k1jfrSyNnXTlnqTbhxy\nWXC7MmqUdzHbf//ie++9B3/+sy9PPeVTlV9/vb9f0N5ePhD5jjv6YDlNfXOuUFxQ4VXhjaPw9v7h\niOKZG6nw1q/fFNyuDBvmvRz22qv4Xgjw4ovV0+386lewKvqHGjgQttuuuhDncgZeEWkY/brgdsXM\nU217u49GVrB8OTz/vKfgQhF+8EHvolYwfLhfhqi8NLHOOmn/VySr+Ghv4R0lXSXdepJu3x8DLp45\n+6Qbh6YruN0ZNMgL6E47lb//1lvw9NPlaXjGjOJYveCz+BamSy+8brih0rCIlFPB7cGIEd7v95Of\nLL4XArz0kqfhJ5/05fHHYebM4j7rr18swLvsAh/7mBfmXBThKD0V0pSSLkq6dSXd+ge8KZ4530lX\nBbcPzGCzzXw55JDi+++840X4iSe8CD/xBFx8sU8XBN4/ePfdfRk/3h/cGDEim/8GEUmfCm6M1l23\nOg2vWAHPPguzZ8Ojj8KsWT5+RCFAbbutF9/dd4dPfMIvabS0dH381BRSnZKukq6SbqxUcBM2cGDx\nssLxx/t7777rBXjWLC/C997r3dXAu7d96lOw337+mPOWW+bkMoSI9EgFNwPrrAP77usLFLuqPfQQ\n3HefL7/4hW9rby+OMXHQQSldgqis8Eq6Srp1JN0kp2D37flJuiq4DaC0q9rRR/vP9vPPF4vvLbfA\nNdf4wxgHHABHHOHjBq+7btYtF5HeUMFtQGb+4MV228FJJ0FHh48bccstPn383Xf7pYqDDvLie9hh\nMQ/aHqWdynSkpIuSbh+SbpLT9RTPnI+km/XtGalBa6vfWPvud2HhQnjkES/Ec+bAl77kg/V861vw\nxhtZt1Q4Q9TjAAAOdElEQVRE1kQFN2fMvPhefDEsWgS//a338T33XJ/1+OST/XpwXVpafGmNFrPy\npbUVWlsxM38qraWwtEZLtG4tvkTrxf2j4xeOF62v3r66HRXHKXwGLVZMkP5G2fZU7jKGUJ6kQ2cx\nzeJJd3USL93eGXxZfZjgabez05fCcaP14vZoqTxOZ0e0+Prq/Ts6fCkcr7B0dPoSHd86A9ZZcv7C\n9mh/6+z0tNsRoKOr9Wjp6IyWgEXbyrZ3BlpW+VL8nsJCtETrnZ7sWzoCLSXbK9cL+xW3+1I4Tssq\nT7HF41cshfNU7rfKfKk4bhxUcHOspcVvpt1zT3GKounTYfvt4dJL/edNRBqHCm4/scMOcO218MIL\nsM8+MGWK9+t95pneHyu0GKHF8p10C8dMmpJuj0m3NOXmOenGQQW3n9l8c7jrLrjhBpg/34vuH/6Q\ndatEBNRLoV8ygy98wYel3G8/70p2551++aEmLeV3sK3y93Ieei+k2XOh9PjqvRAp9l7o2whjftTS\n9ex7L9RPCbcf23RTf5hiyy19eqJXXsm6RSLNTQm3n9tgA7jtNh+j4YQTvA9vj5c2W8t/D+cy6WbR\nR7f0+Eq6kU76Nu5C6Z6NknTrp4TbBLbeGi64wGe1uP/+rFsj0rxUcJvEiSf68JCXX17Dzqt7IUS9\nE6LeA3nqvVCyUt1HV70X0u29UEMf3Tz0XoiDCm6TGDwYJk2CO+6ADz7IujUizUkFt4kccIBPlDlr\n1pr3Cy0tfo0wx0m3pqfRlHRTSbq9eRqtkZNuHFRwm8gee/jrn/6UbTtEmpV6KTSRddf1Xgt/+1sP\nO7ZGd8oLd5wLd8rz1HuhkUYYKz1+E/ZeSGYm4OrvTKv3Qj2UcJtMe7sPeiMi6VPCbTKjRsGrr655\nn9UpqbJvZZ6SbiOOpVt6/CZKusnMj1a6Z36SrhJukxkxAt56K+tWiDSnmhOumbUCs4FXQgiHmNkW\nwM+B9YA5wJdDCCuSaabEZcgQWLash50qr+HmMek28qwRpcdvgqSb7EzApXs2ftLtTcI9DXiuZP1/\ngUtCCFsDbwHHxdguScjgwTUUXBFJRE0J18w2AT4DfBs43bwT5L7AF6NdrgPOA65IoI0SowEDYOXK\nNe8TovRUzA35S7q5mB+t9Pj9OOkmMRMwpJ9041Brwr0U+AbFT2Q94O0QQuGBt5eBMV19o5lNNrPZ\nZjZ76dKldTVW6tfW5g8/iEj6eky4ZnYIsCSEMMfM9untCUII04HpAOPGjUs4LkhPWlqKgbA7obX4\nOx7ymXRzNRNw6fH7Y9JNYCZgSD/pxqGWSwp7Ap81s4OBwcA6wGXAcDNri1LuJoBGW82BWgquiCSj\nx4IbQjgbOBsgSrhfCyEcZWa/AA7HeypMAu5IsJ0SE7Oew1poLfzGz3HSjWPWiJLtSrpRM5o46cah\nnmOdid9Am4df070mniZJkmopuCKSjF49aRZCeAB4IPp6PrBb/E2SJNU0QNbqfrgFOUy6cc6PVrJd\nSTdqRi+Sbpzzo/l2KrZHWxNOunHQk2YiIinRWApNppaEW90PtyA/STeRmYBLtivpRs2oIekmMROw\nb6die7S1gZOuEq6ISEqUcKVKIeEW5DHpJjITsL9Rtl1JN2qGkm5NlHBFRFKihCtVOtui3/AVjwDn\nKunGMWsEKOnGkXRjmDWifL0g3aQbByVcEZGUKOFKlcI13EIiyGXSjXN+NFDSVdKNhRKuiEhKlHCl\nSnEsBZfPpJvATMCgpNuXpJvATMDl6wXJJt04KOGKiKRECVeqhNWBJM9JN4GZgEuPo6Qbbe856SYx\nE7Dvn27SjYMSrohISpRwpUrxGm7Xv+HzkHQTmQkYlHT7kHSTmAnYt6eddOunhCsikhIlXKlSTIf5\nTbpJzAQMSrp9S7rxzwRc2qb0km79lHBFRFKihCtVOlsr+x/mMOkmMBNwaZuUdHuTdBOYHw0ySLr1\nU8IVEUmJEq5UKfTDrR4tSUlXSbcPSTfJmYAhV0lXCVdEJCVKuFKl8hpuLpNuEjMBg5JuX5JuAjMB\nQ/pJNw5KuCIiKVHClSrdXcNV0lXS7UvSTWJ+NEg/6cZBCVdEJCVKuFIlVPwazmPSTWImYD9ndA4l\nXX+tIekmORMw5CvpKuGKiKRECVeqhNau389T0k1kJmBQ0u1L0k1ifjRIPenGQQlXRCQlSrhSpTMK\nEt39NlbSRUm3N0k3yZmAIbWkGwclXBGRlNSUcM1sOHA1sAMeQI4FngduAtqBhcARIYS3EmmlpKow\n40NnlDVzmXQTmQkYlHR7n3QTmQkYMki69as14V4G3BtC2A7YGXgOOAu4L4SwDXBftC4iIt3oMeGa\n2brA3sC/AYQQVgArzGwisE+023XAA8CZSTRS0rU6eZLfpJvMTMCgpNv7pJvETMCQRdKtXy0Jdwtg\nKfBjM3vCzK42s7WBDUIIr0b7vAZs0NU3m9lkM5ttZrOXLl0aT6tFRHKolmu4bcBY4JQQwqNmdhkV\nlw9CCMHMuvxVG0KYDkwHGDduXMK/jiUOxbEUXB6TbhIzAYOSbp+SbgIzAfv+URNSSrpxqOVILwMv\nhxAejdZn4gV4sZltBBC9LomtVSIi/VCPCTeE8JqZvWRmHwkhPA9MAJ6NlknAd6LXOxJtqaSm8kmz\nPCbdJGYCBiXdPiXdRGYChrSTbhxqffDhFGCGmQ0E5gPH4P+332xmxwEvAkfE1ioRkX6opoIbQngS\nGNfFpgnxNkcaQfdjKbg8JN1E5kcDJd2C3iTdRGYChrSTbhz0pJmISEo0loJUqRwPt1Iekm6ss0aA\nkm53aki6icwEXHqclJJuHJRwRURSooQrVbpLlFX7Ra+NmXR7fhoNlHRjs4akm8hMwJB60o2DEq6I\nSEqUcKVK1ZNmuUy6PT+NVraupBuPrpJuAjMB+/bonGkl3Rgo4YqIpEQJV6qEtvJkmcekm8RMwL6/\nkm5NSpNuAjMBQz6TrhKuiEhKlHClSvFJs/wm3VhmjQAl3XqF0LtxF6Bhk24clHBFRFKihCtVqmft\nVdJV0q1DEjMBQ+pJNw5KuCIiKVHClSqhtTJhrt4SvTZ+0k1kJmBQ0q2Hkq4SrohIWpRwpUqxH67L\nY9JNYiZgUNKNRV6TbgyUcEVEUqKEK1UqE26Bkq6SbqyaMOkq4YqIpEQJV6pFvRS6yzb5SLr1zxpR\nvl1JNzF5SboxUMIVEUmJEq5Uq7iGm8+kG9/8aOXblXQT0+BJNw4quFLFurmkkK/C2/ehHUu/S4VX\nhbf42ddPlxRERFKihCtVrK2YNSGfSTfJKdjLtyvpJqbhkm79lHBFRFKihCtVWlv9N3zx93r+km6y\nU7AXW6ek20xJt35KuCIiKVHClSqtbeUpK49JN9kp2EvPrqTr602QdGOghCsikpKaEq6ZTQGOx3/B\nPwUcA2wE/BxYD5gDfDmEsCKhdkqKBgwoxLvy/z3ylHSTnIK99LuUdJso6cagx4RrZmOAU4FxIYQd\n8H/ZI4H/BS4JIWwNvAUcF1+zRET6n1qv4bYBa5nZSmAI8CqwL/DFaPt1wHnAFXE3UNI3oLWy36GS\nrpKukm4cejxSCOEV4HvAIrzQvoNfQng7hFD43/BlYExX329mk81stpnNXrp0aTytFhHJoR4TrpmN\nACYCWwBvA78ADqr1BCGE6cB0gHHjxiX4607iMrCtuydrlHSVdPuQdJNMud7I4rnKzh1v0o1DLVl5\nP2BBCGFpCGElcCuwJzDczAo/gZsAr8TWKhGRfqiWa7iLgPFmNgRYBkwAZgP3A4fjPRUmAXck1UhJ\n16Cqa7iVGj/p9mbchdJ1Jd0Ekm4a13NLj59U0o1BLddwHwVmAo/jXcJa8EsEZwKnm9k8vGvYNbG1\nSkSkH6qpl0II4b+A/6p4ez6wW+wtkswNau0hUq7WuEk3yYkpQUm3V0k3zZ4LpcePO+nGQE+aiYik\nRGMpSJW12lb28jsaL+mmMQU7KOnWlHSz6KNbevyYkm4clHBFRFKihCtVBvc64RY0TtLtywhjvl1J\n188bnSeGpJvp02ilx68z6cZBCVdEJCVKuFJl7dZ6B31T0lXSLabZhhh3ofT4fU26MVDCFRFJiRKu\nVFmrta/XcCsp6SrpFvvh5j7pxkAJV0QkJUq4UmXttuUxH1FJt6mTbiOOpVt6/FqTbgyUcEVEUqKE\nK1WGtES9FGL/vyO9pJvETMC+XUnXzxudp5ak28izRpQev4ekGwclXBGRlCjhSpVhrf8ofyOHSbee\nsXTLj9ndOZV0/bzRedaQdHMxP1rp8btJunFQwhURSYkSrlQZWplwC3KVdOufNaL8mN2dU0nXzxud\np4ukm6uZgEuPX5l0Y6CEKyKSEiVcqbJOy7I175CLpFv/rBG+rqQLdSbdJGYC9gORqMqkGwMlXBGR\nlCjhSpUhLTU+adbASbeep9HK9lPS9f3qSbpJzAQM6SfdGCjhioikRAlXqqzT0k0vhe40YNJNYiZg\nX1fShd4l3URmAi7ZnlrSjYESrohISpRwpcqw3ibcAiXdLs6rpJvITMD+Rtn2PCRdJVwRkZQo4UqV\nYS11zvigpNvFeZs46SYwE7Cv5i/pKuGKiKRECVeqDLNCwshv0k1y1oiy/ZR0fb81Jd0EZgKGfCZd\nJVwRkZQo4UqVYS3R/xadUZxS0u2Wkm5kjUk3gZmAIZdJVwlXRCQlSrhSZWjL4OirqD9uDpNunPOj\ngZJufUk3gZmAS4+To6SrhNtkxo+H007LuhUizclCitXezJYCL6Z2QumNzUMIo7NuhEh/lmrBFRFp\nZrqkICKSEhVcEZGUqOCKiKREBVdEJCUquCIiKVHBFRFJiQquiEhKVHBFRFKigisikhIVXBGRlKjg\nioikRAVXRCQlKrgiIilRwRURSYkKrohISlRwRURSooIrIpISFVwRkZSo4IqIpEQFV0QkJSq4IiIp\nUcEVEUnJ/weCS3jrpPMW6wAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#%% plot the distributions\n", "\n", "pl.figure(1, figsize=(6.4, 3))\n", "pl.plot(x, a, 'b', label='Source distribution')\n", "pl.plot(x, b, 'r', label='Target distribution')\n", "pl.legend()\n", "\n", "#%% plot distributions and loss matrix\n", "\n", "pl.figure(2, figsize=(5, 5))\n", "ot.plot.plot1D_mat(a, b, M, 'Cost matrix M')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve EMD\n", "---------\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcVNWZ//HPQzeL7CDNvrQsxhhBxY6CGsaA4hINGXGU\nDI6MG+OKg8aoITFxMnH5uYcEFZcEE1xxDW4/g7hFQGkRFZVFQARlFUQQge4+88dzO91CA71U31vL\n9/163Vd13bpd9XQ1/eXUueeeYyEERESk/jVIugARkVyhwBURiYkCV0QkJgpcEZGYKHBFRGKiwBUR\niYkCVyRLmNkPzGx+0nXIrilwRSJm9p9m9p6ZfW1mK83sDjNrHT12p5ltirZtZra90v3nYqgtmFnv\n3R0TQngthPCdOrzGCDObZWabzWx19PUFZmbR42ZmN5jZumi7ofwxqR4FrghgZpcBNwCXA62AAUAP\n4EUzaxRCOC+E0DyE0By4Fni4/H4I4fjkKndmll/H778MuB24EegIdADOA44AGkWHjQZ+AhwI9ANO\nAv6rLq+baxS4kvPMrCVwDXBxCOH5EML2EMJS4FSgEDi9Fs95lJktN7OfR63Fz83sJ2Z2gpktMLMv\nzOwXlY4/1MxmmNmG6Ng/mFmj6LFXo8PmRi3q0yo9/xVmthL4U/m+6Ht6Ra/RP7rf2czWmNlRVdTa\nCvgf4IIQwpQQwlfBzQkhjAwhbI0OHQXcHEJYHkJYAdwM/GdN35tcpsAVgcOBJsDjlXeGEDYBzwLH\n1PJ5O0bP2wW4GrgbD+9DgB8AvzKzfaJjS4GxQDtgIDAEuCCqY1B0zIFRi/rhSs/fFm+Jj96h9o+B\nK4C/mllT4E/ApBDCy1XUORBoDDy1h5/ne8DcSvfnRvukmhS4Ih5ya0MIJVU89nn0eG1sB34XQtgO\nPBQ9z+1RC3Ie8AH+8ZwQQnEIYWYIoSRqXd8F/Msenr8M+HUIYWsIYcuOD4YQ7gYWAbOATsC4XTzP\nTj+/mb0Rtba3mFl54DcHvqz0fV8CzdWPW30KXBFYC7TbRT9op+jx2lgXQiiNvi4PxFWVHt+Chxhm\ntq+ZTY1O1m3E+4n3FPRrQgjf7OGYu4EDgPGVugZ2qpMdfv4QwuEhhNbRY+U5sQloWen7WgKbgmbA\nqjYFrgjMALYCJ1feaWbNgeOBaTHUcAfwEdAnhNAS+AWwp5bjboMuqv824F7gN2bWdheHlv/8w/bw\nevOIWuSRA6N9Uk0KXMl5IYQv8ZNm483sODNraGaFwCPAcuAvMZTRAtgIbDKz/YDzd3h8FdCzhs95\nOzA7hHAO8AxwZ1UHhRA24D//BDM7xcxamFkDMzsIaFbp0PuBS82si5l1Bi4D/lzDmnJanYaSiGSL\nEML/M7N1wE1ALzz8ngRG7uajeCr9DJgI/ByYAzwMDK70+G+ASWa2F36CbPXunszMhgHHAX2jXZcC\n75jZyBDC5B2Pj37+FdHr3w9sBhbjJ97eiA67Cw/996L790T7pJpM3S8iIvFQl4KISEwUuCIiMVHg\niojERIErIhITjVIQANq1axcKCwuTLkMkoxQXF68NIRRU93gFrgBQWFjI7Nmzky5DJKOY2Sc1OV5d\nCiIiMVHgiuSKLVvgrbdg0SLQ+PtEKHBFst2cOTB4MLRsCYceCn36QNu2cOml8NVXSVeXUxS4Itmq\nrAyuvBKKimDePLj8cnjsMbj7bjj+eLjtNthvP5g+PelKc4ZOmolko7IyOP98mDgRzj4bbrwR2rSp\nePycc+CSS+DMM+GEE2DqVBgyJLl6c4RauCLZ6OKLPWx/8Qtv0VYO23KHHQavvOJdDCeeCK++uvMx\nklIKXJFsc//9MGECXHYZ/O//wu4WZCgogGnToHt3OO00WL3bScikjhS4Itlk/ny44AIYNAiuv373\nYVuuoAAefRTWr4dRo7w7QuqFAlckW5SWwsiR0KQJTJ4M+TU4RdOvH9x6Kzz/PPzxj/VXY45T4Ipk\ni4kTobjYA7Nr15p//3nnwdCh8MtfwqpVez5eakyBK5IN1q6FceN8vO2pp9buOcxg/Hi/QOLKK1Nb\nnwAKXJHsMG6cX8Qwfnz1+m13Zd99/WTbn/8MM2akrDxxClyRTLdgAdxzj58s23//uj/fuHHQsSNc\ndZUuAU4xBa5IprvmGj9RNm5cap6veXMfv/vKK/DSS6l5TgEUuCKZbd48ePBBv9ChffvUPe+55/qJ\nt1/+Uq3cFFLgimSya67xFunll6f2eZs08bCdOdOHiklKKHBFMtXChTBlClx0Eey9d+qf/8wzoVs3\nv4BCUkKBK5Kpbr4ZGjWCMWPq5/kbNYKxY32OhVmz6uc1cowCVyQTrVrlQ7dGjfIRBfXlnHOgdWuf\nbUzqTIErkonGj4dt23zMbH1q0cKHmz3+uHdhSJ0ocEUyzZYtcMcdMGyYX6hQ3y6+GBo2hN//vv5f\nK8spcEUyzQMPwBdf+ATicejYEUaM8C6MjRvjec0spcAVySQheHdC377wL/8S3+tefDFs2uShK7Wm\nwBXJJK+9BnPnegDWZc6EmioqgoEDPew1X26tKXBFMsn48b5czsiR8b/2mDG+xLouhKg1Ba5Ipvjs\nM3jiCTjrLGjaNP7XHz7c+3PvuCP+184SClyRTHHvvb6qw3/9VzKv37ChrwD8zDPwySfJ1JDhFLgi\nmaCkxFd0OOYYX2U3KaNHe9/x3XcnV0MGU+CKZIJnn4Xly+H885Oto3t3+NGPfP7dbduSrSUDKXBF\nMsGdd0LnznDiiUlX4mufrVoFTz6ZdCUZR4Erku6WLvWRAeec4/2oSTv2WOjRw7s4pEYUuCLp7p57\nvN/07LOTrsTl5Xn4T5vmw8Sk2hS4Iuls+3a47z44/njvP00XZ53lwauTZzWiwBVJZ888A59/ntxQ\nsF3p3BlOOgn+9CedPKsBBa5IOrvrLujSxVu46Wb0aFizRifPakCBK5Kuli6FF17w/tL8/KSr2dnQ\noTp5VkMKXJF0VX6y7Kyzkq6kajp5VmMKXJF0lK4ny3akk2c1osAVSUdTp6bnybId6eRZjShwRdJR\nOp8s25FOnlWbAlck3Sxe7CfLzj03PU+W7WjoUCgs1LSN1aDAFUk3d91VcUIqE+TledfHyy/DRx8l\nXU1aU+CKpJOtW/1k2Y9/7F0KmeKss3yehzvvTLqStKbAFUknjz8Oa9cmPw1jTbVv7ytCTJoEX3+d\ndDVpS4Erkk4mTIBevWDIkKQrqbnzzoMNG+DBB5OuJG0pcEXSxTvvwOuvw4UXQoMM/NMcNAgOOAD+\n8Adfzl12koG/VZEs9Yc/+OKQZ56ZdCW1YwYXXeT/cbzxRtLVpCUFrkg6WLcOJk+G//gPaN066Wpq\n7/TToVUr/89DdqLAFUkH990H33zj3QmZrFkzH7EwZYpfKSffosAVSVpJCYwfD0cdBX37Jl1N3V14\noS/nPmFC0pWkHQWuSNKmTIFPP4VLL026ktTo1QuGDfMrzzRE7FsUuCJJCgFuvhn69PHlx7PFZZd5\nv/T99yddSVpR4Iok6R//gNmzYezYzBwKtitHHAHf/z7ceiuUlSVdTdrIot+wSAa68UZo2xbOOCPp\nSlLLzLtIFiyAp59Oupq0ocAVScp773kYjRnjZ/ezzSmnQM+ecN11uhAiosAVScp110Hz5nDxxUlX\nUj/y8+GKK+DNN30ZHlHgiiRi0SJ4+GGfpKZt26SrqT+jRvmqENdem3QlaUGBK5KE667z6QzHjk26\nkvrVuDH87GcwfbrPE5HjFLgicVuwwKcxPO886NQp6Wrq3+jR0KED/PKXOd+Xq8AVidtvfuMtv6uu\nSrqSeDRrBuPGwSuv5HxfrgJXJE7vvQcPPQSXXOKtvlwxerQv9z5uXE63chW4InG64gpo2dL7NXNJ\n48Zw9dU+YuGxx5KuJjEKXJG4PPecb7/6VXaPTNiVUaN8cp7LL/eZ0XKQAlckDtu3+5VXvXtn77jb\nPcnPh9tug6VL/ZLfHKTAFYnDH//oS4jffDM0apR0NckZPNhnEvvd72DFiqSriZ0CV6S+LVvmQ6KO\nOw5OOinpapJ3yy0+oc2FF+bcCTQFrkh9CqFiyfM77/RJXXJdz55wzTXw1FO+LHwOUeCK1KcHH4Rn\nn/WP0D16JF1N+hg7Fg4+2Bed/OKLpKuJjQJXpL4sWeKt24EDPVikQn6+r+O2bh2cc07OdC0ocEXq\nQ0kJjBzpXz/wAOTlJVtPOjroIJ9T4okn4K67kq4mFgpckfpw1VUwY4YHSWFh0tWkr7Fj4dhj/ba4\nOOlq6p0CVyTVJk2Cm27y7oQRI5KuJr01aODrnrVv78PFsnxpdQWuSCq99prPGzBkCNx+e9LVZIb2\n7X3liw0bPHQ3bUq6onqjwBVJlTff9JV399kHHnnE57uV6jnwQJg82bsVfvxj2LIl6YrqhQJXJBVm\nz/a+yIICn4IwF+dKqKthw7w75uWX4eST4euvk64o5RS4InX13HNw1FHQqpWHbZcuSVeUuU4/He6+\nG154wbtl1q5NuqKUUuCK1FYIPhnLSSfBvvv6qASNSKi7s8+GKVPgnXdgwAC/zRIKXJHaWL3aPwKP\nHQsnnuirGeTCcjlxOflkeOkl78sdMADGj/f5FzKcAlekJkpL/SPvfvvB88/7SIQnnoAWLZKuLPsM\nHOit28GDYcwYOPJImDs36arqRIErUh2lpT7yoF8/H/bVr5//8Y8Zowlp6lNBATzzDPz5z7Bwoc+/\n8O//Dh9+mHRltaLAFdmdJUt84pmePeG003zfI4/4st/f/W6yteUKM18tYv58X6Lo6adh//3hmGP8\nd7F5c9IVVpuFHJk0QnavqKgozJ49O+kykrdxo4+nnT7duwzeftv3DxkCF1zg/baaFyFZa9bAxIm+\nLVsGe+0FQ4fC0UfDoEEexvn5sZRiZsUhhKJqH6/AFcihwN261acDXLUKVq70P9glS7z1NG8eLFjg\nx+Xnw6GHwr/+Kwwf7hczSHopLfUr+6ZM8Skwlyzx/Xvt5V0+3/0u9Onj02J26+arJBcU+PC9FP2n\nqcCVWvlW4D75JEydWv8vuqt/e+X7K9/uaist9bPXJSX+9fbtvm3d6tuWLb599ZVvVV3BlJ/va419\n97vQvz98//tw+OE6EZZpliyBN96At97y5eg//LDquRnM/HfbogU0awZNm0KTJr6ycKNGfoVgfr5v\neXk+38OgQf4JZ6enqlngxtPulsyyaJF/nI7Drk44le+vfFvVVv4HUf7H0bChb40bQ/Pm3trZa6+K\nP7A2bfwqsPbtvcXTvTt07qxugmywzz6+lU+LCd6/u2yZr5+2apVfSPHFF/Dll/4f8ObNfkXbN9/4\nf9BffeX/YZeU+FZW5luKJo9XC1eAHOpSEEmhmrZwNUpBRCQmClwRkZioS0EAMLM1wCeVdrUD0nnm\nkHSuL51rg/SuL51rg53r6xFCKKjuNytwpUpmNrsmfVNxS+f60rk2SO/60rk2qHt96lIQEYmJAldE\nJCYKXNmViUkXsAfpXF861wbpXV861wZ1rE99uCIiMVELV0QkJgpcEZGYKHBlJ2Z2nJnNN7NFZnZl\nwrV0M7PpZvaBmc0zs0ui/W3N7EUzWxjdtkmwxjwzm2NmU6P7+5jZrOj9e9jMGiVYW2szm2JmH5nZ\nh2Y2MM3eu7HR7/V9M3vQzJok+f6Z2X1mttrM3q+0r8r3y9zvozrfNbP+e3p+Ba58i5nlAX8Ejgf2\nB35qZvsnWFIJcFkIYX9gAHBhVM+VwLQQQh9gWnQ/KZcAlZcguAG4NYTQG1gPnJ1IVe524PkQwn7A\ngXidafHemVkXYAxQFEI4AMgDRpDs+/dn4Lgd9u3q/Toe6BNto4E79vjsIQRt2v65AQOBFyrdvwq4\nKum6KtXzFHAMMB/oFO3rBMxPqJ6u0R/hYGAqYPiVSPlVvZ8x19YKWEJ0crzS/nR577oAnwJt8ZkL\npwLHJv3+AYXA+3t6v4C7gJ9WddyuNrVwZUflfwTllkf7EmdmhcDBwCygQwihfLLTlUCHhMq6Dfg5\nUL6k7N7AhhBCSXQ/yfdvH2AN8Keoy+MeM2tGmrx3IYQVwE3AMuBz4EugmPR5/8rt6v2q8d+KAlcy\ngpk1Bx4D/juEsLHyY8GbF7GPbzSzE4HVIYTiuF+7mvKB/sAdIYSDgc3s0H2Q1HsHEPWFDsP/Y+gM\nNGPnj/Nppa7vlwJXdrQC6FbpftdoX2LMrCEetpNDCI9Hu1eZWafo8U7A6gRKOwL4sZktBR7CuxVu\nB1qbWfnk/km+f8uB5SGEWdH9KXgAp8N7B3A0sCSEsCaEsB14HH9P0+X9K7er96vGfysKXNnRW0Cf\n6ExxI/wkxtNJFWNmBtwLfBhCuKXSQ08Do6KvR+F9u7EKIVwVQugaQijE36eXQggjgenAKUnWFtW3\nEvjUzL4T7RoCfEAavHeRZcAAM2sa/Z7L60uL96+SXb1fTwNnRKMVBgBfVup6qFoSneXa0nsDTgAW\nAB8D4xKu5Uj8I9y7wDvRdgLeVzoNWAj8HWibcJ1HAVOjr3sCbwKLgEeBxgnWdRAwO3r/ngTapNN7\nB1wDfAS8D/wFaJzk+wc8iPcnb8c/IZy9q/cLP0H6x+jv5D18tMVun79Ol/aa2XH4R6g84J4QwvW1\nfjIRkSxX68CNxmsuwIfoLMc/iv40hPBB6soTEckedVm191BgUQhhMYCZPYSfcdxl4LZr1y4UFhbW\n4SWlrjZs8IVJu3X79v7i4uK1IZq5/pgG/6YZjUR28GLZo7tYYrr66hK4VY1BO2zHg8xsNH4VBt27\nd0crwybr5z+H8eNhx1+DmX1S9XeISKrU+yiFEMLEEEJRCKGooKDaS/9IPdm+HRo2TLoKkdxUl8BN\nu/GasmfbtkGjxKZSEcltdQnctBqvKdXz9dfQrFnSVYjkplr34YYQSszsIuAFfFjYfSGEeSmrTOrF\nxo3QvHnSVYjkprqcNCOE8CzwbIpqkRisWwd77510FSK5SZf25pgVK6Bjx6SrEMlNCtwcUloKS5dC\nr15JV5IgM99EEqDAzSHvvw8lJXDAAUlXIpKb6tSHK5nllVf8dtCgZOtIVPml7OWtXIvaHKHs24+L\n1AO1cHPIX/8K+++/82W9IhIPtXBzxMyZ8NZbflmvUNGSDaV+2yAPAGvkfxJhe7TCS1lp3JVJFlML\nNwds2wbnnw8dOsAZZyRdjUjuUgs3y4UA48bBO+/Ak09Cy5ZJV5SmopZs2Oq3lu9/GraXX5YXtm33\n2+3bEihOsoVauFksBJ8d7KabvIU7bFjSFYnkNrVws9SXX8KYMXD//XDhhfD73yddUWYJJSXfum3Q\npInftvLL9MLXWwAo+/rrBKqTTKUWbhZ68kkfjfDXv8LVV/uJsgb6TYskTi3cLDJjBlx3Hfztb9Cv\nnwfv97+fdFXZoeybb/yL6DYv6gzP7+Fj7MLGTQCUrl8ff3GSMdTuyXClpfDEE3DEEXD44fD663Dt\ntb6ig8JWJL2ohZuhPvwQHn4YJk+GRYugsND7ac88U9MvxqF040b/IrrNi1YzadBvP79dswGAks9X\nxl+cpC0FbgZZsMBD9pFHfF4EM79M93e/g5NPhnz9NkXSmv5E09imTfDaa/D3v8OLL8J773nIHnmk\nnwgbPhw6dUq6SgEoXbPGv4huGxR29/0/7A9Ao+Xe4i1duDj+4iRtKHDTyPbtMGuWB+y0aX45bkmJ\nr0F2xBFw221wyinQpUvSlYpIbShwE7R6tQds+TZjBmze7K3YQw6Byy6Do4/2sN1rr6SrlZooWboM\ngLzolu99B4CvRgwAoOXCrwAIxVqVKpcocGOydatfXjtzpofrzJmwZIk/lpcHffvCqFEwZAgcdRS0\nbZtouSJSDxS49WDTJnj3XQ/YOXP89t13fRIZ8C6BAQP8ctvDDvPWrFbSzW6l8+YD0CJq0JYdfiAA\nq8YeDkC7uVsByH+pOP7iJDYK3DpaterbwTpnDixcWDH7X9u2cPDBcMklHq6HHQZduyZbs4gkQ4Fb\nTVu2wAcf+EiBytvKSsMsCws9XEeOhIMO8q+7dtUSWrIze2MuAB3f8PtbT/CrVD65cSAAnWb4ChRN\nH58Vf3FSbxS4OygthcWLdw7WRYugLFqFpUkTn6vg2GMrgvXAA6F162RrF5H0ltOBu3r1zsE6bx6U\nTwBl5ivc9u0LI0b4bd++0Lu3n+gSSZXGz74FQK9n/f7Gf/fRDEse6gdAq//vnfxt75sRf3GSMjkR\nuF9/7UG6Y7iuXl1xTEGBh+m551YE6/e+p5NZIpI6WRe4X3zhJ67efrvidsGCipNYTZr4MuE/+lFF\nsPbt68vPiKSLlg/MjG79/uoLfDRDs1d9zoaPn+gDQMdb34i/OKm1jA3cEOCzz3YO12XLKo7p1g36\n9/92d0CvXuoOEJFkZEzglpV5t8Crr/r8Aq++Cp9/XvH4vvvCwIG+usHBB/vWrl1y9YqkUvsJ3pLd\nPMHvl17pLdwT5/n8u3c89CMAuv1WLd50lraBu307FBdXhOs//gHlczt36eJXYw0Y4C3YAw+EFi0S\nLVdEZI/SKnBD8Am077/fpyAsn3J03319+sEf/MCnIyws1NhWyW1drveW7NTr2wBQ+ls/SXHDEh+3\ne+rD/w3APldqVEM62WPgmlk34H6gAxCAiSGE282sLfAwUAgsBU4NIdRqfZHFi+Evf/GgXbzYRwYM\nHw4nneRTEXbsWJtnFRFJL9Vp4ZYAl4UQ3jazFkCxmb0I/CcwLYRwvZldCVwJXFHTAm6+GX72M2+x\nDhkCv/mNt2Y1HEuk+gp/5S3ZK351GABlN/r+x5b7aIeDohZvr8tmxl+c/NMeAzeE8DnwefT1V2b2\nIdAFGAYcFR02CXiZGgbuxIketsOHw623+qgCEZFsVaM+XDMrBA4GZgEdojAGWIl3OVTbkiVw3nlQ\nVAQPPOCTbItIavS63Fu8wy/3K9bCrb7/hc/e8ccfOg+A3peqxRunaq/aa2bNgceA/w4hbKz8WAgh\n4P27VX3faDObbWaz15QvQ4JP6vLDH/r42WefrV3xIiKZpFqBa2YN8bCdHEJ4PNq9ysw6RY93AlZX\n9b0hhIkhhKIQQlFBtLIpQMOG8OST3sIdPhxOOMEXSPzmmzr9PCJShd5jZ9J77EyO7XwQx3Y+CAtg\nwft4H1s+k49vHMjH0UxlUn/2GLhmZsC9wIchhFsqPfQ0MCr6ehTwVE1fvEULeP55uOIKn9tgxAgf\nkTB6tI+7LZ+dS0QkG1gIVfYEVBxgdiTwGvAeUB6Bv8D7cR8BugOf4MPCvtjdcxUVFYXZs2dX+Vhp\nKbz8MkyaBI895hPOtGnjw8IGDfIxuP37e8tYUs/MikMIRQDHNPi33f+jkKyz5Hpv3T5y2m0AnPZX\nH9VQPvpB4MWyR+s8+r86oxReB3b1QkPqWkC5vDwfFjZkCEyYAE89BdOn+5Vmf/ubH9O0qV++W34B\nxKGHaviYiGSOPbZwU2l3LdzdWbnSr0Arn0dh7ly/Ks0MvvMdnzehf/+KORS0AGPNqYUrlS272mcn\nu+C0ZwC46y8+V0P5FW65KJYWbjro2BFOOcU3gA0b4I034K23fJTD66/Dgw9WHN+jR0UA9+/vqzJ0\n7qzLgUUkWRkRuDtq3dpHNZxwQsW+tWs9fMu3t9/2URDlDfg2bSqmaOzXz28POECT3ohUpfv/RHM1\n/I/P1ZA31veXz8e79MHeABTcoT7emsiILoXa+uor736YO9dHQbz7Lrz/vu8vV1j47YnI+/b1yXJy\n7eScuhSkJr44y0+yfTl0MwB7P90UqJg4PRvlTJdCbbVo4aMcjjyyYl8I8MknOy+389xzUFLixzRq\nBPvtt3MQawVeEamLrG7h1sTWrTB/vreCKwfx8uUVx7Ru7d0QO3ZNtGyZXN2pohau1MXm4T5pzsoB\nPrS/63RvvZQvjpkN1MJNocaNPUD79fv2/vXrvRuicghPnlwxVy/4Kr7ly6WX33bsqNawiHybWri1\nEAJ8+qm3ht95x7c5c3wu33Lt21cE8EEHwSGHeDCnawirhSupVDL4EADWHtiYDrO8n9femJtkSXWm\nFm5CzKB7d99OPLFi/5dfegjPmVMRwrfc4ssFgY8PPuww3wYM8As32rRJ5mcQkfgpcFOoVSu/Cu4H\nP6jYt20bfPABzJ4Ns2bBzJk+f0T5B4t99/XwPewwOPxw79JoUO053ETSU/5LxQB0fAnskO8BsHGE\nTxXZ+j1fGKZ03vxkikuQAreeNWpU0a1wzjm+b+NGD+CZMz2En3/elxcCX2n4hz+Eo4/2y5x79kzf\nbggRqRkFbgJatoTBg32DiqFqr74K06b59uij/lhhYcUcE8cdpy4IyTyheB4ALYqjHX16AlD6w/4A\nNF7s82SXfPJp7LXFTYGbBsw8WAsL4YwzPIDnz68I38ceg3vv9Ysxhg6FU0+FYcO8C0NEMocCNw2Z\n+YUX++0HF17oU1e+9ZYH7yOPwDPPeFfFccd5+A4fDk2aJF21SPWULvThPHkLox1dOgPQoN9+ANjn\n6/y4SivEZAudnskAeXl+Yu3GG2HpUpgxw4O4uBhOP90n6/ntb2HduqQrFZHdUeBmGDMP31tugWXL\n4O9/9zG+V1/tqx5fdJH3B4tkipIVn1Gy4jPK3v2Isnc/8mvsS0rI79GN/B7daNCiBQ2yZJYpBW4G\na9DAT6Y9+2zFEkUTJ8L++8Ntt3lXhIikDwVuljjgALjvPli4EI46CsaO9XG98+YlXZlIzZSuX0/p\n+vWUfPKpj1woK4OyMvLa7U1eu71p0KQJDTL0pIUCN8v06AFTp8IDD/ilxocf7gtyikjyFLhZyAx+\n+lOfhL1jRx9KNm1a0lWJ1E7Z5s2Ubd5M6dp1lK5dRwiBEAINmjWjQbNmWH4+lp8ZA64UuFmsWze/\nmKJnT1+eaMWKpCsSyW0K3CzXoQM88YTP93vuuRVzOIhkqrB1K2Hr1n+2fMtZ48ZY48b+ES9Nr4dX\n4OaA3r3h2mt9VYvp05OuRiR3KXBzxHnn+fSQEyYkXYlIaoWSEt+ilu8/NcjzLY0ocHNEkyYwahQ8\n9RRU+hQFYQepAAAGdElEQVQmIjFS4OaQoUP9Ip6Z2buwqoifqAgBykp9K+/TTYO+XQVuDhnoK1vz\n5pvJ1iGSqzJj8JqkRKtWPmrh44+TrkQkRmk0NEct3BxTWOiT3ohI/BS4OaZdO03jKJIUBW6OadMG\n1q9PugqR3FTtwDWzPDObY2ZTo/v7mNksM1tkZg+bWaP6K1NSpWlT2LIl6SpEclNNWriXAB9Wun8D\ncGsIoTewHjg7lYVJ/WjSRIErkpRqBa6ZdQV+BNwT3TdgMDAlOmQS8JP6KFBSq2FD2L496SpEclN1\nW7i3AT8HyqL7ewMbQggl0f3lQJeqvtHMRpvZbDObvSYLF4XLNPn5fvGDiMRvj4FrZicCq0MIxXs6\ntiohhIkhhKIQQlFBQUFtnkJSqEEDn0BfROJXnQsfjgB+bGYnAE2AlsDtQGszy49auV0BzbaaARS4\nIsnZYws3hHBVCKFrCKEQGAG8FEIYCUwHTokOGwU8VW9VSsqYpdWFNyI5pS7jcK8ALjWzRXif7r2p\nKUnqkwJXJDk1mkshhPAy8HL09WLg0NSXJPUpTSfCF8kJutJMRCQmCtwcoxauSHIUuCIiMVHgiojE\nRIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIi\nMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6I\nSEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMalW4JpZazObYmYfmdmHZjbQzNqa2YtmtjC6bVPf\nxYqIZLLqtnBvB54PIewHHAh8CFwJTAsh9AGmRfdFRGQX9hi4ZtYKGATcCxBC2BZC2AAMAyZFh00C\nflJfRYqIZIPqtHD3AdYAfzKzOWZ2j5k1AzqEED6PjlkJdKjqm81stJnNNrPZa9asSU3VIiIZqDqB\nmw/0B+4IIRwMbGaH7oMQQgBCVd8cQpgYQigKIRQVFBTUtV4RkYxVncBdDiwPIcyK7k/BA3iVmXUC\niG5X10+JIiLZYY+BG0JYCXxqZt+Jdg0BPgCeBkZF+0YBT9VLhSIiWSK/msddDEw2s0bAYuBMPKwf\nMbOzgU+AU+unRBGR7FCtwA0hvAMUVfHQkNSWIyKSvXSlmYhITBS4IiIxUeCKiMREgSsiEhMFrohI\nTBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsi\nEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCK\niMREgSsiEhMFrohITBS4IiIxqVbgmtlYM5tnZu+b2YNm1sTM9jGzWWa2yMweNrNG9V2siEgm22Pg\nmlkXYAxQFEI4AMgDRgA3ALeGEHoD64Gz67NQEZFMV90uhXxgLzPLB5oCnwODgSnR45OAn6S+PBGR\n7LHHwA0hrABuApbhQfslUAxsCCGURIctB7pU9f1mNtrMZpvZ7DVr1qSmahGRDFSdLoU2wDBgH6Az\n0Aw4rrovEEKYGEIoCiEUFRQU1LpQEZFMV50uhaOBJSGENSGE7cDjwBFA66iLAaArsKKeahQRyQrV\nCdxlwAAza2pmBgwBPgCmA6dEx4wCnqqfEkVEskN1+nBn4SfH3gbei75nInAFcKmZLQL2Bu6txzpF\nRDJe/p4PgRDCr4Ff77B7MXBoyisSEclSutJMRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgo\ncEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQm\nClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJ\niQJXRCQmClwRkZgocEVEYqLAFRGJiQI3xwwYAJdcknQVIrnJQgjxvZjZGuCT2F5QaqJHCKEg6SJE\nslmsgSsiksvUpSAiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMF\nrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMRE\ngSsiEhMFrohITP4PrQ161dhEqnEAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#%% EMD\n", "\n", "G0 = ot.emd(a, b, M)\n", "\n", "pl.figure(3, figsize=(5, 5))\n", "ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve Sinkhorn\n", "--------------\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "It. |Err \n", "-------------------\n", " 0|8.187970e-02|\n", " 10|3.460174e-02|\n", " 20|6.633335e-03|\n", " 30|9.797798e-04|\n", " 40|1.389606e-04|\n", " 50|1.959016e-05|\n", " 60|2.759079e-06|\n", " 70|3.885166e-07|\n", " 80|5.470605e-08|\n", " 90|7.702918e-09|\n", " 100|1.084609e-09|\n", " 110|1.527180e-10|\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVfWd//HXZxq9ykjXERsxgkpQscQ1Yk00ZKNrdM1K\nbKyJUZeYRA3ZbPztRpONNSQWLAkmhGiwBks2MdgFBbEhggiIImXonWnf3x+fc5nCDNPPueX9fDzO\n43LvPXPvhwvzns98z/d8j4UQEBGR9peXdAEiIrlCgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMF\nrkjCzOyLZragHV73AjP7vybu+y0ze7m5z0nzKHAl40WB8K6ZbTOzlWZ2l5n1jJ6728y2RFuZmZXX\nuP9MDLUFMztgT/uEEF4KIRzcwtc/3sxeNbONZrbOzF4xsyOj150SQji1Ja8r7UOBKxnNzK4BfgH8\nAOgBjAL2Bf5mZkUhhMtDCF1DCF2BG4GHUvdDCGckV7kzs4JWfG13YDowEegNDARuAHa2TXVtrzV/\n32ygwJWMFQXODcCVIYRnQwjlIYSlwLlACfDNFrzmiWb2qZn90MxWm9kKM/uamX3ZzBZGXeSPaux/\nlJm9ZmYbon1/bWZF0XMvRru9HXXU36jx+tea2Urgt6nHoq/ZP3qPEdH9AWZWamYn1lPuQQAhhKkh\nhMoQwvYQwv+FEN6JvrbWUEDUbV9uZh9G9f7GzKyBz+GXZvaymfWo8djNZrbezJaY2Rk1Hh9gZk9G\ndS8ys8tqPPdTM5tmZn8ws03At6LHHjazB81ss5nNM7ORzfynykgKXMlkxwIdgUdrPhhC2AI8DZzS\nwtftF73uQOAnwL14eH8B+CLwn2a2X7RvJTAe6AMcA4wGvhPVcUK0z2FRR/1QjdfvjXfi4+rU/hFw\nLfAHM+sM/BaYHEJ4vp46FwKVZjbZzM4ws15N+LudCRwJDMd/MJ1W80kzyzOze6PnTw0hbIyeOhpY\nEP09/xe4v0ZY/wn4FBgAnAPcaGYn1XjZMcA0oCcwJXrsq9HX9QSeBH7dhNozngJXMlkfYE0IoaKe\n51ZEz7dEOfCzEEI5Hgp9gDtCCJtDCPOA94HDAEIIc0IIM0MIFVF3fQ/wT428fhXwXyGEnSGE7XWf\nDCHcCywCZgH9gQn1vUgIYRNwPBDwHwqlUafZdw/v/fMQwoYQwjJgBnB4jecKgan4D4OzQgjbajz3\ncQjh3hBCJTA5qquvmQ0GjgOuDSHsCCG8BdwHXFjja18LITweQqiq8fd9OYTwdPR6vyf6PLOdAlcy\n2RqgTwPjgv2j51tibRQEAKmAWFXj+e1AVwAzO8jMpkcH6zbh48SNBX1pCGFHI/vcCxwKTAwhNDgm\nG0KYH0L4VghhULT/AOD2Pbzuyhp/3pb6e0QOwLvRG0IIZQ19XY0g7hq937oQwuYa+36M/3aQ8kkT\n6uiYC+O7ClzJZK/hB4i+XvNBM+sKnAE8F0MNdwEfAAeGELoDPwLqHRetYY9L9EX13w7cD/zUzHo3\npZAQwgfA7/DgbYn5wEXAM2bW1FkTnwG9zaxbjcf2AZbXLK2F9WQdBa5krGh88QZgopmdbmaFZlYC\nPIyPKf4+hjK6AZuALWY2FPh2nedXAUOa+Zp3ALNDCJcCTwF317eTmQ01s2vMbFB0fzBwPjCzme+3\nSwhhKv5D4+9mtn8T9v8EeBW4ycw6mtlw4BLgDy2tIZspcCWjhRD+Fw+Im/Hgm4X/Cjt6T7+Kt6Hv\nA/8KbMaHAR6q8/xPgcnRrIBzG3sxMxsDnE51cH8PGGFmF9Sz+2b8YNYsM9uKB+17wDUt+HvsEkKY\nDPw/4B/RD7DGnI/PCvkMeAwfn/57a2rIVqYFyEVE4qEOV0QkJgpcEZGYKHBFRGKiwBURiUnWTzSW\npunTp08oKSlJugyRjDJnzpw1IYTipu6vwBUASkpKmD17dtJliGQUM/u4OftrSEFEJCYKXJFcsX07\nvPEGLFoEmn+fCAWuSLabOxdOOgm6d4ejjoIDD4TeveF734PNmxv/emkzClyRbFVVBdddByNHwrx5\n8IMfwCOPwL33whlnwO23w9ChMGNG0pXmDB00E8lGVVXw7W/DpElwySXwy19Crxrrk196KVx9NVx0\nEXz5yzB9OowenVy9OUIdrkg2uvJKD9sf/cg72l71XAzi6KPhhRd8iOHMM+HFF3ffR9qUAlck2zz4\nINx5J1xzDfzP/0D9ly1zxcXw3HOwzz7wjW/A6tXx1ZmDFLgi2WTBAvjOd+CEE+DnP99z2KYUF8Of\n/wzr18PYsT4cIe1CgSuSLSor4YILoGNHmDIFCppxiGb4cLjtNnj2WfjNb9qvxhynwBXJFpMmwZw5\nHpiDBjX/6y+/HE49FX78Y1i1qvH9pdkUuCLZYM0amDDB59ue2+iFJepnBhMn+gkS113XtvUJoMAV\nyQ4TJvhJDBMnNm3ctiEHHeQH2373O3jttTYrT5wCVyTTLVwI993nB8sOOaT1rzdhAvTrB9dfr1OA\n25gCVyTT3XCDHyibMKFtXq9rV5+/+8IL8I9/tM1rCqDAFcls8+bB1Kl+osPee7fd6152mR94+/GP\n1eW2IQWuSCa74QbvSH/wg7Z93Y4dPWxnzvSpYtImFLgimerDD2HaNPjud2Gvvdr+9S+6CAYP9hMo\npE0ocEUy1S23QFERXHVV+7x+URGMH+9rLMya1T7vkWMUuCKZaNUqn7o1dqzPKGgvl14KPXv6amPS\nagpckUw0cSKUlfmc2fbUrZtPN3v0UR/CkFZR4Ipkmu3b4a67YMwYP1GhvV15JRQWwq9+1f7vleUU\nuCKZ5o9/hHXrfAHxOPTrB+ed50MYmzbF855ZSoErkklC8OGEYcPgn/4pvve98krYssVDV1pMgSuS\nSV56Cd5+2wOwNWsmNNfIkXDMMR72Wi+3xRS4Iplk4kS/XM4FF8T/3ldd5ZdY14kQLabAFckUn30G\njz0GF18MnTvH//5nn+3juXfdFf97ZwkFrkimuP9+v6rDv/97Mu9fWOhXAH7qKfj442RqyHAKXJFM\nUFHhV3Q45RS/ym5Sxo3zseN7702uhgymwBXJBE8/DZ9+Ct/+drJ17LMPfOUrvv5uWVmytWQgBa5I\nJrj7bhgwAM48M+lK/Npnq1bB448nXUnGUeCKpLulS31mwKWX+jhq0k47Dfbd14c4pFkUuCLp7r77\nfNz0kkuSrsTl53v4P/ecTxOTJlPgiqSz8nJ44AE44wwfP00XF1/swauDZ82iwBVJZ089BStWJDcV\nrCEDBsBZZ8Fvf6uDZ82gwBVJZ/fcAwMHeoebbsaNg9JSHTxrBgWuSLpauhT++lcfLy0oSLqa3Z16\nqg6eNZMCVyRdpQ6WXXxx0pXUTwfPmk2BK5KO0vVgWV06eNYsClyRdDR9enoeLKtLB8+aRYErko7S\n+WBZXTp41mQKXJF0s3ixHyy77LL0PFhW16mnQkmJlm1sAgWuSLq5557qA1KZID/fhz6efx4++CDp\natKaAlcknezc6QfLvvpVH1LIFBdf7Os83H130pWkNQWuSDp59FFYsyb5ZRiba++9/YoQkyfDtm1J\nV5O2FLgi6eTOO2H//WH06KQrab7LL4cNG2Dq1KQrSVsKXJF08dZb8PLLcMUVkJeB35onnACHHgq/\n/rVfzl12k4H/qiJZ6te/9otDXnRR0pW0jBl897v+g+PVV5OuJi0pcEXSwdq1MGUK/Nu/Qc+eSVfT\nct/8JvTo4T88ZDcKXJF08MADsGOHDydksi5dfMbCtGl+ppzUosAVSVpFBUycCCeeCMOGJV1N611x\nhV/O/c47k64k7ShwRZI2bRp88gl873tJV9I29t8fxozxM880RawWBa5IkkKAW26BAw/0y49ni2uu\n8XHpBx9MupK0osAVSdIrr8Ds2TB+fGZOBWvIccfBkUfCbbdBVVXS1aSNLPoXFslAv/wl9O4NF16Y\ndCVty8yHSBYuhCefTLqatKHAFUnKu+96GF11lR/dzzbnnANDhsBNN+lEiIgCVyQpN90EXbvClVcm\nXUn7KCiAa6+F11/3y/CIAlckEYsWwUMP+SI1vXsnXU37GTvWrwpx441JV5IWFLgiSbjpJl/OcPz4\npCtpXx06wPe/DzNm+DoROU6BKxK3hQt9GcPLL4f+/ZOupv2NGwd9+8KPf5zzY7kKXJG4/fSn3vld\nf33SlcSjSxeYMAFeeCHnx3IVuCJxevdd+NOf4OqrvevLFePG+eXeJ0zI6S5XgSsSp2uvhe7dfVwz\nl3ToAD/5ic9YeOSRpKtJjAJXJC7PPOPbf/5nds9MaMjYsb44zw9+4Cuj5SAFrkgcysv9zKsDDsje\nebeNKSiA22+HpUv9lN8cpMAVicNvfuOXEL/lFigqSrqa5Jx0kq8k9rOfwfLlSVcTOwWuSHtbtsyn\nRJ1+Opx1VtLVJO/WW31BmyuuyLkDaApckfYUQvUlz+++2xd1yXVDhsANN8ATT/hl4XOIAlekPU2d\nCk8/7b9C77tv0tWkj/Hj4Ygj/KKT69YlXU1sFLgi7WXJEu9ujznGg0WqFRT4ddzWroVLL82ZoQUF\nrkh7qKiACy7wP//xj5Cfn2w96ejww31Nicceg3vuSbqaWChwRdrD9dfDa695kJSUJF1N+ho/Hk47\nzW/nzEm6mnanwBVpa5Mnw803+3DCeeclXU16y8vz657tvbdPF8vyS6srcEXa0ksv+boBo0fDHXck\nXU1m2Htvv/LFhg0eulu2JF1Ru1HgirSV11/3K+/utx88/LCvdytNc9hhMGWKDyt89auwfXvSFbUL\nBa5IW5g928cii4t9CcJcXCuhtcaM8eGY55+Hr38dtm1LuqI2p8AVaa1nnoETT4QePTxsBw5MuqLM\n9c1vwr33wl//6sMya9YkXVGbUuCKtFQIvhjLWWfBQQf5rATNSGi9Sy6BadPgrbdg1Ci/zRIKXJGW\nWL3afwUePx7OPNOvZpALl8uJy9e/Dv/4h4/ljhoFEyf6+gsZToEr0hyVlf4r79Ch8OyzPhPhsceg\nW7ekK8s+xxzj3e1JJ8FVV8Hxx8PbbyddVasocEWaorLSZx4MH+7TvoYP92/+q67SgjTtqbgYnnoK\nfvc7+PBDX3/hX/8V5s9PurIWUeCK7MmSJb7wzJAh8I1v+GMPP+yX/f7c55KtLVeY+dUiFizwSxQ9\n+SQccgiccor/W2zdmnSFTWYhRxaNkD0bOXJkmD17dtJlJG/TJp9PO2OGDxm8+aY/Pno0fOc7Pm6r\ndRGSVVoKkyb5tmwZdOoEp54KJ58MJ5zgYVxQEEspZjYnhDCyyfsrcAVyKHB37vTlAFetgpUr/Rt2\nyRLvnubNg4ULfb+CAjjqKPjnf4azz/aTGSS9VFb6mX3TpvkSmEuW+OOdOvmQz+c+Bwce6MtiDh7s\nV0kuLvbpe230Q1OBKy1SK3AffxymT2//N23o/17q8Zq3DW2VlX70uqLC/1xe7tvOnb5t3+7b5s2+\n1XcGU0GBX2vsc5+DESPgyCPh2GN1ICzTLFkCr74Kb7zhl6OfP7/+tRnM/N+2Wzfo0gU6d4aOHf3K\nwkVFfoZgQYFv+fm+3sMJJ/hvOLu9VPMCN56+WzLLokX+63QcGjrglHq85m19W+obIvXNUVjoW4cO\n0LWrdzudOlV/g/Xq5WeB7b23dzz77AMDBmiYIBvst59vqWUxwcd3ly3z66etWuUnUqxbBxs3+g/g\nrVv9jLYdO/wH9ObN/gO7osK3qirf2mjxeHW4AuTQkIJIG2puh6tZCiIiMVHgiojEREMKAoCZlQIf\n13ioD5DOK4ekc33pXBukd33pXBvsXt++IYTipn6xAlfqZWazmzM2Fbd0ri+da4P0ri+da4PW16ch\nBRGRmChwRURiosCVhkxKuoBGpHN96VwbpHd96VwbtLI+jeGKiMREHa6ISEwUuCIiMVHgym7M7HQz\nW2Bmi8zsuoRrGWxmM8zsfTObZ2ZXR4/3NrO/mdmH0W2vBGvMN7O5ZjY9ur+fmc2KPr+HzKwowdp6\nmtk0M/vAzOab2TFp9tmNj/5d3zOzqWbWMcnPz8weMLPVZvZejcfq/bzM/Sqq8x0zG9HY6ytwpRYz\nywd+A5wBHAKcb2aHJFhSBXBNCOEQYBRwRVTPdcBzIYQDgeei+0m5Gqh5CYJfALeFEA4A1gOXJFKV\nuwN4NoQwFDgMrzMtPjszGwhcBYwMIRwK5APnkezn9zvg9DqPNfR5nQEcGG3jgLsaffUQgjZtuzbg\nGOCvNe5fD1yfdF016nkCOAVYAPSPHusPLEionkHRN+FJwHTA8DORCur7PGOurQewhOjgeI3H0+Wz\nGwh8AvTGVy6cDpyW9OcHlADvNfZ5AfcA59e3X0ObOlypK/VNkPJp9FjizKwEOAKYBfQNIaQWO10J\n9E2orNuBHwKpS8ruBWwIIVRE95P8/PYDSoHfRkMe95lZF9LkswshLAduBpYBK4CNwBzS5/NLaejz\navb3igJXMoKZdQUeAf4jhLCp5nPB24vY5zea2ZnA6hDCnLjfu4kKgBHAXSGEI4Ct1Bk+SOqzA4jG\nQsfgPxgGAF3Y/df5tNLaz0uBK3UtBwbXuD8oeiwxZlaIh+2UEMKj0cOrzKx/9Hx/YHUCpR0HfNXM\nlgJ/wocV7gB6mllqcf8kP79PgU9DCLOi+9PwAE6Hzw7gZGBJCKE0hFAOPIp/puny+aU09Hk1+3tF\ngSt1vQEcGB0pLsIPYjyZVDFmZsD9wPwQwq01nnoSGBv9eSw+thurEML1IYRBIYQS/HP6RwjhAmAG\ncE6StUX1rQQ+MbODo4dGA++TBp9dZBkwysw6R//OqfrS4vOroaHP60ngwmi2wihgY42hh/olMViu\nLb034MvAQuAjYELCtRyP/wr3DvBWtH0ZHyt9DvgQ+DvQO+E6TwSmR38eArwOLAL+DHRIsK7DgdnR\n5/c40CudPjvgBuAD4D3g90CHJD8/YCo+nlyO/4ZwSUOfF36A9DfR98m7+GyLPb5+q07tNbPT8V+h\n8oH7Qgg/b/GLiYhkuRYHbjRfcyE+RedT/FfR80MI77ddeSIi2aM1V+09ClgUQlgMYGZ/wo84Nhi4\nffr0CSUlJa14S2mtDRv8wqSDB9d+fM6cOWtCtHL9KXn/ohWNROr4W9WfG7jEdNO1JnDrm4N2dN2d\nzGwcfhYG++yzD7oybLJ++EOYOBHq/jOY2cf1f4WItJV2n6UQQpgUQhgZQhhZXNzkS/9IOykvh8LC\npKsQyU2tCdy0m68pjSsrg6LEllIRyW2tCdy0mq8pTbNtG3TpknQVIrmpxWO4IYQKM/su8Fd8WtgD\nIYR5bVaZtItNm6Br16SrEMlNrTloRgjhaeDpNqpFYrB2Ley1V9JViOQmndqbY5Yvh379kq5CJDcp\ncHNIZSUsXQr77590JQkwa94m0g4UuDnkvfegogIOPTTpSkRyU6vGcCWzvPCC355wQrJ1pAVrpNdI\nNbmhqvbjrVh7REQdbg75wx/gkEN2P61XROKhDjdHzJwJb7zhp/XmpFRnmhqfTXWuUadreVbrPqn7\nVXU62ujrQurx1Ouo85UmUIebA8rK4Nvfhr594cILk65GJHepw81yIcCECfDWW/D449C9e9IVJaxu\np9sASz1fVP+3iEWvEyor/YGqOvfV+Uo91OFmsRB8dbCbb/YOd8yYpCsSyW3qcLPUxo1w1VXw4INw\nxRXwq18lXVGaSXWeoTK66x2t5VP7frS7FUbfKvnRDlEHnOp0iTrbUBl1tuXlte6r8xVQh5uVHn/c\nZyP84Q/wk5/4gbI8/UuLJE4dbhZ57TW46Sb4y19g+HAP3iOPTLqqDFFnTHbXNNzU89FPLIs6XEst\nKlwQfQtFsxosNXuhosK/PtXp7izzx1P3U8/v6nzV8eYC9T0ZrrISHnsMjjsOjj0WXn4ZbrzRr+ig\nsBVJL+pwM9T8+fDQQzBlCixaBCUlPk570UVafrFVGph9sJu8qNONOtzQIVrVPRrrDflRx1vpX29l\n3tnaDu90w/bt0e0Ovy1LjflqrDebKXAzyMKFHrIPP+zrIpj5abo/+xl8/evVv92KSHrSt2ga27IF\nXnoJ/v53+Nvf4N13PWSPP94PhJ19NvTvn3SVWaruLIayqONMdbxVqQ7Ub1NnqoWow63q7B1vVQfv\nhENqjLfCvz5/u3e0eZu9w83bGnW8W7f5baoD1lhvVlHgppHycpg1ywP2uef8dNyKCr8G2XHHwe23\nwznnwMCBSVcqIi2hwE3Q6tUesKnttddg61bvYr/wBbjmGjj5ZA/bTp2SrjbHpcZ2K2qPtealOs9o\nvm1qlkJeNE+3qtD/4cq7+bdaRaeoE87rAEB+uV9grnCLv07R+p3+9Ru2+utt2uL7b/H7u8Z6ozrU\n8WYWBW5Mdu7002tnzvRwnTkTlizx5/LzYdgwGDsWRo+GE0+E3r0TLVdE2oECtx1s2QLvvOMBO3eu\n377zji8iAz4kMGqUn2579NHezepKuhmizthu1c6os011uNE/cl4037awvCLa36eOVBV5Z7u9u3e6\nFV18ZmaIZj3k7/D5vR03+H+ITqW9ACha42O7+Ws3+f5bvPOtSs1ySL1PVWUb/CWlvShwW2nVqtrB\nOncufPhh9fdl795wxBFw9dUerkcfDYMGJVuziCRDgdtE27fD++/7TIGa28qV1fuUlHi4XnABHH64\n/3nQIF0iK6ulxnbLo/m1qTUVotkFedEZZkXRbf7OHgBYpY/tbin0znZHsf8n2RbNOtkc/acp2Ood\nccc1fttlpS/31vkz72wLV23099sYdb7RLIeqstQYr+bzphMFbh2VlbB48e7BumhR9Uygjh19rYLT\nTqsO1sMOg549k61dRNJbTgfu6tW7B+u8ebDNmwTM/Aq3w4bBeef57bBhcMAB1YtGidRSVXveblVq\nbYad0eyDaMy161b/6VywtRsAm8t83u6Wwd7Z7ugfjckO8s5488F+d8Mm36/jZz7G221Z5+h2LwA6\nrPBON3/NBq9ja2p2Q+0OXB1vMnIicLdt8yCtG66rV1fvU1zsYXrZZdXB+vnP62CWiLSdrAvcdev8\nwNWbb1bfLlxY/QO9Y0e/TPhXvlIdrMOG+eVnRNrMrnm7tc8Us2g2QarT7bTJO9CiDT4bocNG/wm/\ncZt/a27Z31+u90DvWAcPXAFAxVCf3bB0nc8fXPuxd8rdFhcD0HOxd9CdPtkMQP7q9QBUbY7m9UYd\ntzreeGVs4IYAn322e7guW1a9z+DBMGJE7eGA/ffXcICIJCNjAreqyocFXnzR1xd48UVYsaL6+YMO\ngmOO8asbHHGEb336JFevSC11ZjNU1h3bjWYX9Fjnsxg6rfaOd2NpRwA2HOxjtIsO9nm6xw1aDMDJ\nxfMB2HaAz2J49YghAMxb7Od/d1ngHXCvhf66XZZ6h7ur490QzXJQxxuLtA3c8nKYM6c6XF95Bdb7\n/xEGDvSzsUaN8g72sMOgW7dEyxURaVRaBW4IvoD2gw/6EoSb/IArBx3kyw9+8Yu+HGFJiea2SoZL\nzWbY6beVu8Z2fZWwwk0+9tpnjXeoXVZ4h7r+M7/9v88PA2D5UB+rPbffGwDcsO+7AKwd5GPB04cf\nDsCzHx4CQME8//peC/zMt26LvZPOX6mONw6NBq6ZDQYeBPriVxyZFEK4w8x6Aw8BJcBS4NwQwvqW\nFLF4Mfz+9x60ixf7zICzz4azzvKlCPv1a8mrioikFwuN/OQys/5A/xDCm2bWDZgDfA34FrAuhPBz\nM7sO6BVCuHZPrzVy5Mgwe/bsWo/dcgt8//vesY4eDRde6N2spmPFy8zmhBBGApyS9y9qZ5KS+tXN\nfBZCXkcfm7Vu3pGGvt7xbt3PzzhbN9R7pu3DvDP+ytD3ABi71ysAHFzo84E/izrVqRv8uksPfTjC\nX+8df53eH/h+3T6KZjWo493N36r+3OrfqxvtcEMIK4AV0Z83m9l8YCAwBjgx2m0y8Dywx8Cta9Ik\nD9uzz4bbbvNZBSIi2arRDrfWzmYlwIvAocCyEELP6HED1qfuN6Rmh7tkiU/R+sIX/IBYUVHL/gLS\nNtThpqmo4911teAO3vHmdfejxJX9fPbClv29A143NLrCxGHeqZ530BwALuz5OgD7Ffp+S8p9tsKD\nG44CYOqCL/jrvuOv2/sD72S7fhTNali5FoCqaGw5Fzvetuhwm3zVXjPrCjwC/EcIYVPN54Kndr2f\nuJmNM7PZZja7tLR01+ODBsGXvuTzZ59+umXFi4hkkiZ1uGZWCEwH/hpCuDV6bAFwYghhRTTO+3wI\n4eA9vU7dMdzNm+GUU+CNN3whmLFjYcwYPxtM4qUON0M01PH29NkHlf181sHm/b1TXRedkWbDvUf6\nxoFvAnB+T5/VMCB6nQXlvt/ktccB8NQHhwLQ6V1f1az3Bz6LossSfx1btQ6AEJ25VrXDO95sXp0s\nlg43Gi64H5ifCtvIk8DY6M9jgSea++bdusGzz8K11/raBued5zMSxo3zYYbU6lwiItmgKbMUjgde\nAt4FUhH4I2AW8DCwD/AxPi1s3Z5eq75ZCimVlfD88zB5MjzyiC8406uXTws74QSfgztiBBQWNuev\nJ02lDjdD1e14o4vfVY/xRh3vEB+7XX+w71fxeV/D4fQD3wfgzJ5vAbBXvj/+7k5fJX/aSh/bnTff\nj2j3eN+Ps/da6GfMdVrmsxgojTre6Npru9bjzaIrUMQ1S+FloKE3Gt3aAlLy831a2OjRcOed8MQT\nMGOGn2n2l7/4Pp07++m7qRMgjjpK08dEJHM0a5ZCa+2pw92TlSv9DLTUOgpvv+1DRGZw8MG+bsKI\nEdVrKOgCjM2nDjdLNDbGu7d3vFv3izreA32/rQd7x/r5IcsBOLa3r9XQOd/HZudv9UtRvPKpr9VQ\ntsDn7/Zc4G/b4yNf/axouc/fDeu9863ati1rrrcWS4ebDvr1g3PO8Q1gwwZ49VU/2DZ3rofx1KnV\n+++7b3UAjxjhV2UYMECnA4tIsjKiw22KNWs8fFPbm2/Wvphjr17VSzQOH+63hx6qRW9S1OFmqYY6\n3ujMtarjtFm6AAAMPElEQVSo490+2L8RNgzxgySbh/jhmk77+rzbkt4+Rltg/vgnm7xjXrfcp953\n/ch7tx6Lo/m7H/tYbv7K9YRo7u6uKwxXpK63lln/zXKmw22KPn18itkpp1Q/tnmzDz+8/bbPgnjn\nHV+vYfPm6n1KSmovRD5smC+Wo4NzkhUaWAg9deKCRWHYeYUfDOkcLVy+s78PGWzex28/2scDdscA\nf52C7mW1blMLpVd09kjZ2cMDvHuPDnSIXjtvbXTZn9Qi6Dl42Z+sCdz6dOvmsxyOP776sRDg4493\nv9zOM89A9H+SoiIYOnT3INYVeEWkNbI6cOtj5l1tSYmvRpaycycsWOBdcCqEX3gBpkyp3qdnTx+G\nqDs00b173H8LkRZqoOMlOnEhL/r1r8Nqvzhlx6XRwbG+3uFuG+BnJW3t58/vSC3y38Vft6Kz324d\n6J1JZVERXbt619y5iw9n5EeLqoeN0VBDdNXWbDm4tic5F7gN6dDBA3T48NqPr18P771XuxueMqV6\nrV7wq/imLpeeuu3XT92wiNSmwG1Er14+7/eLX6x+LAT45BPvht96y7c334Rp06r32Xvv6gA+/HBf\npOeAAxTCkmZS46bBu8qqHVF3GZ24YNGJDAXrfLpXj0+ihcv38s53Z3Gq0/WDHjt6Rp1tdHp+VZGx\nrU+0oE5+1DV38PuFHXzFqvyNfrvrApfRe2fqwbU9UeC2gBnss49vZ55Z/fjGjR7Cc+d6CM+dC7fe\n6pcLAp8ffPTRvo0a5Sdu9OqVzN9BROKnwG1DPXrs3g2XlcH778Ps2TBrFsyc6etHpH5oH3SQh+/R\nR8Oxx/qQRl6T13ATaSd1LwEUzSiwrb7Qua33GQcdV/oMhI7do2lmPf1+WS8fry3vmk9loXe9VdHV\nsst6Rmuxms9kKCj0GMqLpgaFaBw57PBvhF2nCWfBwjgK3HZWVFQ9rHDppf7Ypk0ewDNnegg/+6xP\nVwOf3valL8HJJ/tpzkOGaBhCJFsocBPQvTucdJJvUD1V7cUX4bnnfPvzn/25kpLqNSZOP11DEJKQ\nOpd5T42vpi7zbhv9KHLeGl88p1OXaLy2SyequvmAbmWnqIMtiDqI6De5qq7e8Vrw7tjyoucLfPZC\nXuqEiSwY21XgpoGaU9UuvND/Hy1YUB2+jzwC99/vJ2Oceiqce66vG9yjR9KVi0hzKHDTkJmfeDF0\nKFxxhS9d+cYbHrwPPwxPPeVDFaef7uF79tlatF1i1tB83jqzG6yokLxO/p8zP1o6MnSMxnCLoo43\nv/aYWYhmL1hVnQ42dXCjzPffNW83g8Z2dXgmA+Tn+4G1X/4Sli6F117zIJ4zB775TV+s57//G9au\nTbpSEdkTBW6GMfPwvfVWWLYM/v53n+P7k5/4VY+/+10fDxaJVQi+VVVCVSWhvIxQXkbVtm1Ubdjo\nW+kaqkrXEFaWElaWwqo1sGoNeaUbyCvdgG3a6tvOMmxn2a6XtsJC34pSW5FvhQW+5ef74jxmaX+E\nWYGbwfLy/GDa009XX6Jo0iQ45BC4/XYfihCR9KHAzRKHHgoPPOBLUp54Iowf7/N6581LujLJaSEQ\nKioIFRVU7dzp29Ztvm3c7Nv6Db5FnXDYtMW3bdt8Ky8npM4eAiw/z7eCAqygwMfcoi433TtdBW6W\n2XdfmD4d/vhHWLzYQ/eVV5KuSkRAsxSykhmcf74vS3nyyT6V7MknffhBJDF11m0IqbPZKqJuNDXD\nIVosnbzai6c31LVa9Hiw6JLw+am3qTOmlgazGNThZrHBg/1kiiFD/PJEy5cnXZFIblPgZrm+feGx\nx3y938suS4sf8iK1NTDDIZT5VrVjJ1U7dhJ2RltZee0tBOq9VJjl+bbrfvJjuwrcHHDAAXDjjX5V\nixkzkq5GJHcpcHPE5Zf78pB33pl0JSJNVLfzjWY7hIpy3yorfSuv8C11vyoQap6lVrfTTVB6VCHt\nrmNHGDsWnngCtm5NuhqR3KTAzSGnnuoXypw5M+lKRFqhTue7awtVe95SEhzLVeDmkGOO8dvXX0+2\nDpFcpXm4OaRHD5+18NFHSVci0g4yYAqOOtwcU1Lii96ISPwUuDmmTx8t4yiSFAVujunVC9avT7oK\nkdzU5MA1s3wzm2tm06P7+5nZLDNbZGYPmVlR+5UpbaVzZ9i+PekqRHJTczrcq4H5Ne7/ArgthHAA\nsB64pC0Lk/bRsaMCVyQpTQpcMxsEfAW4L7pvwEnAtGiXycDX2qNAaVuFhVBjaVERiVFTO9zbgR8C\nqdnDewEbQgjRVdz4FBhY3xea2Tgzm21ms0tLS1tVrLReQYGf/CAi8Ws0cM3sTGB1CGFOS94ghDAp\nhDAyhDCyuLi4JS8hbSgvD6qqGt9PRNpeU058OA74qpl9GegIdAfuAHqaWUHU5Q4CtNpqBlDgiiSn\n0Q43hHB9CGFQCKEEOA/4RwjhAmAGcE6021jgiXarUtqMWUackCOSlVozD/da4Htmtggf072/bUqS\n9qTAFUlOs9ZSCCE8Dzwf/XkxcFTblyTtKU0vZiqSE3SmmYhITBS4OUYdrkhyFLgiIjFR4IqIxESB\nKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR\n4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhM\nFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFpUuCaWU8zm2ZmH5jZfDM7xsx6m9nfzOzD6LZXexcr\nIpLJmtrh3gE8G0IYChwGzAeuA54LIRwIPBfdFxGRBjQauGbWAzgBuB8ghFAWQtgAjAEmR7tNBr7W\nXkWKiGSDpnS4+wGlwG/NbK6Z3WdmXYC+IYQV0T4rgb71fbGZjTOz2WY2u7S0tG2qFhHJQE0J3AJg\nBHBXCOEIYCt1hg9CCAEI9X1xCGFSCGFkCGFkcXFxa+sVEclYTQncT4FPQwizovvT8ABeZWb9AaLb\n1e1ToohIdmg0cEMIK4FPzOzg6KHRwPvAk8DY6LGxwBPtUqGISJYoaOJ+VwJTzKwIWAxchIf1w2Z2\nCfAxcG77lCgikh2aFLghhLeAkfU8NbptyxERyV4600xEJCYKXBGRmChwRURiosAVEYmJAldEJCYK\nXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJ\nAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURi\nosAVEYmJAldEJCYKXBGRmDQpcM1svJnNM7P3zGyqmXU0s/3MbJaZLTKzh8ysqL2LFRHJZI0GrpkN\nBK4CRoYQDgXygfOAXwC3hRAOANYDl7RnoSIima6pQwoFQCczKwA6AyuAk4Bp0fOTga+1fXkiItmj\n0cANISwHbgaW4UG7EZgDbAghVES7fQoMrO/rzWycmc02s9mlpaVtU7WISAZqypBCL2AMsB8wAOgC\nnN7UNwghTAohjAwhjCwuLm5xoSIima4pQwonA0tCCKUhhHLgUeA4oGc0xAAwCFjeTjWKiGSFpgTu\nMmCUmXU2MwNGA+8DM4Bzon3GAk+0T4kiItmhKWO4s/CDY28C70ZfMwm4FviemS0C9gLub8c6RUQy\nXkHju0AI4b+A/6rz8GLgqDavSEQkS+lMMxGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAV\nEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChw\nRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYK\nXBGRmChwRURiosAVEYmJAldEJCYK3BwzahRcfXXSVYjkJgshxPdmZqXAx7G9oTTHviGE4qSLEMlm\nsQauiEgu05CCiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgi\nIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWu\niEhMFLgiIjH5/xjX3J71f+OyAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#%% Sinkhorn\n", "\n", "lambd = 1e-3\n", "Gs = ot.sinkhorn(a, b, M, lambd, verbose=True)\n", "\n", "pl.figure(4, figsize=(5, 5))\n", "ot.plot.plot1D_mat(a, b, Gs, 'OT matrix Sinkhorn')\n", "\n", "pl.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.12" } }, "nbformat": 4, "nbformat_minor": 0 }