summaryrefslogtreecommitdiff
path: root/notebooks/plot_OT_1D.ipynb
blob: 93d156d746dbdb5dac52426c8399a3caa3cd793b (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
{
 "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 <remi.flamary@unice.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "import ot.plot\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": "iVBORw0KGgoAAAANSUhEUgAAAZQAAADFCAYAAABzYARGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4U2X2wPHvoWW1yFZEsSCgiGylQEUERHBQdFBRB0cEFVdcBmVUEHQQEZfB0UH9CeIwMo47KCqDK4iAICrSYgsUZEepiiJiEdm6vL8/ToKltCVt09ws5/M8edIkNzcnN2nOfXdxzmGMMcZUVBWvAzDGGBMdLKEYY4wJCksoxhhjgsISijHGmKCwhGKMMSYoLKEYY4wJCksoxhhjgsISijHGmKCwhGKMMSYo4r0OoKjExETXrFkzr8MwxhgDpKen/+ScaxjItmGXUJo1a0ZaWprXYRhjjAFE5OtAt7UqL2OMMUFhCcUYY0xQBJRQRORcEVkrIhtEZHQxj1cXkRm+x5eKSDPf/VVF5HkRWSkia0Tk7uCGb4wxJlwcsQ1FROKAycDZQDawTERmO+dWF9rsOmCnc+4kERkIPAJcBlwKVHfOtReRWsBqEXnVObcl2G/EGFN+ubm5ZGdns2/fPq9DMR6pUaMGSUlJVK1atdz7CKRRvguwwTm3CUBEpgP9gcIJpT8wzvf3TGCSiAjggKNEJB6oCRwAdpU7WhOQOXNg9WoYPhyqWKWmCUB2dja1a9emWbNm6L+uiSXOOXbs2EF2djbNmzcv934C+bk5Htha6Ha2775it3HO5QE5QAM0ufwGfA98AzzmnPu56AuIyFARSRORtO3bt5f5TZjfPfMM/PGPcMcdMGgQ7N/vdUQmEuzbt48GDRpYMolRIkKDBg0qXEKt7PPXLkA+0BhoDtwpIi2KbuScm+qcS3XOpTZsGFB3Z1OEczBmDNx8M5x3Hjz4IMyYoX/n5HgdnYkElkxiWzA+/0CqvL4FmhS6neS7r7htsn3VW3WAHcAg4APnXC7wo4gsAVKBTRUN3Bzqtttg0iS49lr4178gPh6aNtXbPXvCkiWQkOB1lMaYaBZICWUZ0FJEmotINWAgMLvINrOBIb6/BwDznS5W/w1wFoCIHAV0Bb4KRuDmd19/DZMnw9Ch8OyzmkwArrwSXn8dVqyAV17xNkZjjuShhx6ibdu2JCcnk5KSwtKlS70O6TAJvrOy7777jgEDBpS43S+//MLTTz9d6r66desGwMKFCzn//PPLFMesWbNYvfr3ZuyxY8cyb968Mu2jMhwxofjaRIYBc4A1wGvOuSwRGS8iF/o2mwY0EJENwB2Av2vxZCBBRLLQxPScc25FsN9ErHvuOb2++24oWmrt3x/atoVp00IflzGB+uyzz3jnnXdYvnw5K1asYN68eTRp0uTITzyCvLy8IER3uMaNGzNz5swSHy8tofhj+vTTT8v9+kUTyvjx4+nTp0+59xcsAU294px7D3ivyH1jC/29D+0iXPR5u4u73wRPfr4mlD59oLgp0ETguuu0kX7VKmjXLuQhmgjz179CRkZw95mSAk88UfLj33//PYmJiVSvXh2AxMTEg4999NFHjBgxgry8PE499VSmTJlC9erVD07TlJiYSFpaGiNGjGDhwoWMGzeOjRs3smnTJpo2bcpLL73EqFGj+OCDD6hSpQo33HADt956K+np6dxxxx3s3r2bxMRE/vvf/3LccccdEtfmzZsZNGgQu3fvpn///gfv37JlC+effz6rVq0iKyuLa665hgMHDlBQUMAbb7zBvffey8aNG0lJSeHss8+mX79+3HvvvdSrV4+vvvqKdevWkZCQwO7duwHYtWsX/fr1Y8OGDfTu3Zunn36aKlWqHLLNzJkzeeeddxg6dCizZ8/m448/5sEHH+SNN97ggQce4Pzzz2fAgAGlHq8hQ4bw9ttvk5uby+uvv84pp5wSrI8YsJHyEe+jj+CbbzRplOTKK6FqVSulmPB1zjnnsHXrVk4++WRuueUWPv74Y0B7n1199dXMmDGDlStXkpeXx5QpU464v9WrVzNv3jxeffVVpk6dypYtW8jIyGDFihUMHjyY3Nxcbr31VmbOnEl6ejrXXnstf/vb3w7bz/Dhw7n55ptZuXLlYcnG75lnnmH48OFkZGSQlpZGUlISEyZM4MQTTyQjI4NHH30UgOXLl/Pkk0+ybt26w/bxxRdf8NRTT7F69Wo2btzIm2++WeJ769atGxdeeCGPPvooGRkZnHjiiQcfO9LxSkxMZPny5dx888089thjRzyOZRV2k0Oaspk2DerXh4suKnmbxESt+nrxRZgwAXwngcYUq7SSRGVJSEggPT2dxYsXs2DBAi677DImTJhAx44dad68OSeffDIAQ4YMYfLkyfz1r38tdX8XXnghNWvWBGDevHncdNNNxPsaF+vXr8+qVatYtWoVZ599NgD5+fnFJowlS5bwxhtvAHDllVcyatSow7Y5/fTTeeihh8jOzuaSSy6hZcuWxcbUpUuXEsd4dOnShRYttAPs5ZdfzieffFJqG01J1q5dW+rxuuSSSwDo3LlzqUmrvKyEEsF27IBZs7QEcqQkcf31uv3sot0pjAkTcXFx9OrVi/vvv59JkyYd/CEvSXx8PAUFBQCHjZ846qijSn2uc462bduSkZFBRkYGK1euZO7cucVue6TutIMGDWL27NnUrFmTP/7xj8yfP7/Y7UqLqehr+G8Xvj8Ysxj4qxTj4uIqpX3JEkoEe+klOHCg9Oouvz59oEkTq/Yy4Wnt2rWsX7/+4O2MjAxOOOEEWrVqxZYtW9iwYQMAL774ImeeeSagS12kp6cDlJp8zj77bP71r38d/AH9+eefadWqFdu3b+ezzz4DdOqZrKysw57bvXt3pk+fDsDLL79c7P43bdpEixYtuO222+jfvz8rVqygdu3a/PrrrwG//y+++ILNmzdTUFDAjBkz6NGjBwCNGjVizZo1FBQU8NZbbx3cvqT9l3a8QsESSoRyTrsIn3oqtG9/5O3j4uCaa2DuXG1zMSac7N69myFDhtCmTRuSk5NZvXo148aNo0aNGjz33HNceumltG/fnipVqnDTTTcBcN999zF8+HBSU1OJi4srcd/XX389TZs2JTk5mQ4dOvDKK69QrVo1Zs6cyahRo+jQoQMpKSnF9rp68sknmTx5Mu3bt+fbb4sOv1OvvfYa7dq1IyUlhVWrVnHVVVfRoEEDunfvTrt27Rg5cuQR3/+pp57KsGHDaN26Nc2bN+fiiy8GYMKECZx//vl069btkCq5gQMH8uijj9KxY0c2btx48P7SjlcoiA4XCR+pqanOFtg6spUrITkZpkyBQL8vW7ZA8+bajlJMVbCJYWvWrKF169Zeh2E8Vtz3QETSnXOpgTzfSigRyl9N269f4M9p1kzHpCxYUCkhGWNinCWUCLVgAZx4oraLlEXv3vDJJ5CbWzlxGWNilyWUCFRQAIsWQa9eZX9ur17w229gtYrGmGCzhBKBMjNh587yJRR/hw+r9jLGBJsllAi0cKFelyehJCZqrzD/PowxJlgsoUSgBQvgpJMgKal8z+/VS6ezP3AgqGEZY2KcJZQIk5+v7Se9e5d/H717w549sGxZ8OIypiJ27NhBSkoKKSkpHHvssRx//PEHbx+opDOf5cuX88EHHwS0bY8ePcjwzZjZt2/fUgctTpw4sdRR7ddccw1r164lLy+PunXrVijmt9566+BcYeHA5vKKMBkZugJjeaq7/Hr21OsFC6B796CEZUyFNGjQ4OAP9rhx40hISGDEiBEBPz8/P7/UwY3FWb58OatWreLcc88t0/PmzJlT6uMTJ07k2muvpUaNGsXG+ZxvvYnyTH1SNGb/AMhwYQklwlSk/cSvQQMdFLlwoS4bbMwhvJi/vhQXXHAB3333Hfv27eP222/n+uuvJy8vj8TERK6++mrmz5/Pv/71L7Zv387IkSNJSEigW7dubN26lVmzZrF7926GDRvG6tWryc3NPbh2yPjx49m7dy8LFy5kzJgxh0zGuGfPHoYMGcKqVato06bNISWOpKQkVq1aRVxcHH/+85/57rvvyM/PZ9y4cWzdupUff/yRM844g0aNGvHBBx8cFufIkSOZNGkS7XxrSdx222189NFHNG7cmOnTp9OgQQN69OjBpEmTSElJYdu2bfTo0YOVK1ceFvMvv/zCqlWreOKJJ9i8eTPXXnstO3bsoFGjRjz33HMkJSVxxRVX0KBBA5YtW8a2bdv45z//WWmJyKq8IsyCBXDyydC4ccX207u3tqPs3x+cuIypLM8//zzp6eksW7aMiRMnsnPnTgBycnLo2bMnK1asoEOHDtxyyy3MnTuXtLQ0tm3bdvD548eP59xzz+WLL75g/vz53HnnnYgIY8eOZfDgwWRkZBw2s++kSZOoV68ea9asYcyYMXz55ZeHxfXee+/RrFkzMjMzD85cfPvtt3PMMcewePHigysoFo7z9NNPP2QfOTk5dO/enaysLE4//XQeeOCBEo9DzZo1S435lltu4frrr2fFihVceumlh8zI/OOPP7JkyRJmzZrF3XffHeCRLzsroUSQvDxYvBgGDqz4vnr1giefhC++gDPOqPj+TBTxYv76Ujz++OPM9k2TnZ2dfXDhqmrVqh080169ejWtWrXihBNOAHQK+BdeeAGAuXPn8v777zNhwgRAZ+395ggT2i1atIi77roLgI4dO9K2bdvDtklOTmb06NGMHj2aCy64gO4l1B8XjrOo+Ph4Lr1U1yC84oorGDRoUKlxlWbp0qW88847AFx11VXce++9Bx+76KKLEBGSk5NLnJMsGKyEEkG+/BJ27apYg7xfz566mqONRzHhbN68eSxatIjPP/+czMxMkpOTD1Y/1axZ84hTy4NOVT9r1qyDU9V/8803B9cLqYjWrVuTlpZG27ZtGT16NA8//HCx2wUaJ/w+XX1pU/OXR/VC61tU5vyNllAiyKJFeh2M2ajr19d2FP8+jQlHOTk51K9fn5o1a5KVlcWyEromtmnThrVr17J161acc8yYMePgY3379uWpp546eNtffVXaFPM9e/bklVdeASAzM7PYqe2//fZbEhISuPLKK7nzzjtZvnz5EfdbVF5e3sGFrl555ZWD09YXnpq/8Nr1pe27a9euvPbaawC89NJL9PT3vgkhSygRJC0NmjaFElYiLbOuXXWfYTbhtDEH9evXjz179tCmTRvGjBnDaaedVux2tWrVYtKkSfTp04fU1FTq1q1LnTp1AJ3m/rfffqN9+/a0bduWcePGAXDWWWeRmZlJx44dD/nRBhg2bBg7duygdevWPPDAA3Ts2PGw18zMzOTUU08lJSWFhx9+mHvuuQeAoUOH0qdPH/r06XPE91enTh0WL15M27Zt+eSTTxjj6yUzcuRInnzySTp16nSwzehIMU+ePJmpU6eSnJzMjBkzePzxx4/4+sFm09dHkJYttVRxhIXsAvbss3DDDbB+vQ6UNLErGqav3717NwkJCTjnuPHGG2nfvj233nqr12FFFJu+Pkbs3AkbNkBqQB9rYPz7svxtosGUKVNISUmhTZs27N27lxtuuMHrkGKO9fKKEL7q2aAmlLZtdS36tLTg9BwzxksjR44MaHVEU3mshBIh/KWIzp2Dt8+qVXW8mZVQDFRu7x8T/oLx+VtCiRBpadCihfbOCqbUVEhP1zVWTOyqUaMGO3bssKQSo5xz7Nixo9jpYsrCqrwiRFoadOkS/P2mpsLkybBuHZxySvD3byJDUlIS2dnZbN++3etQjEdq1KhBUnmnMPexhBIBfvoJtmyBm28O/r79VWhpaZZQYlnVqlVp3ry512GYCGdVXhHAN74pqA3yfq1bQ82a1o5ijKm4gBKKiJwrImtFZIOIjC7m8eoiMsP3+FIRaVbosWQR+UxEskRkpYhUrJIuBvkTSqdOwd93fDx07Pj7axhjTHkdMaGISBwwGTgPaANcLiJtimx2HbDTOXcS8DjwiO+58cBLwE3OubZALyA3aNHHiLQ0HdRYxrV4Apaaqt2S8/MrZ//GmNgQSAmlC7DBObfJOXcAmA70L7JNf+B5398zgT+IznJ2DrDCOZcJ4Jzb4Zyzn60ySkurnOouv9RUXcHxq68q7zWMMdEvkIRyPLC10O1s333FbuOcywNygAbAyYATkTkislxE7iruBURkqIikiUia9TI51A8/wNatlZ9QwNpRjDEVU9mN8vFAD2Cw7/piEflD0Y2cc1Odc6nOudSGDRtWckiRpTIb5P1OPhkSEiyhGGMqJpCE8i3QpNDtJN99xW7jazepA+xASzOLnHM/Oef2AO8BldC0HL3S0nTdkmImOw2auDht8LeEYoypiEASyjKgpYg0F5FqwEBgdpFtZgNDfH8PAOY7HXI7B2gvIrV8ieZMYHVwQo8NaWnQqhXUrl25r5OaqsuI5+VV7usYY6LXEROKr01kGJoc1gCvOeeyRGS8iFzo22wa0EBENgB3AKN9z90JTESTUgaw3Dn3bvDfRvTKyKic7sJFdewI+/bB2rWV/1rGmOgU0Eh559x7aHVV4fvGFvp7H3BpCc99Ce06bMro55+1Qb5Dh8p/Lf9rZGbqLMTGGFNWNlI+jK1YodehSCinnALVqmlCMcaY8rCEEsYyMvQ6JaXyX6tqVS2Z+F/TGGPKyhJKGMvMhEaN9BIKHTpYCcUYU36WUMJYZmZoqrv8OnTQgZTbtoXuNY0x0cMSSpjKzYWsrNBUd/n5X8tKKcaY8rCEEqa++goOHAh9CQUsoRhjyscSSpjy/6iHMqHUqwdNmlhCMcaUjyWUMJWRAdWr6yj5UEpJsZ5expjysYQSpjIzoV07XQArlDp00NHy+/aF9nWNMZHPEkoYci70Pbz8OnTQhbayskL/2saYyGYJJQxt2wbbt4e2h5ef9fQyxpSXJZQw5G/D8KKE0qKFro1i7SjGmLKyhBKG/KWD5OTQv3aVKtC+vZVQjDFlZwklDGVmQrNmULeuN6+fkqIxOOfN6xtjIpMllDCUkeFNdZdfhw6QkwNff+1dDMaYyGMJJczs3Qvr1nmfUMCqvYwxZWMJJcysWgUFBd708PJr317XsbeEYowpC0soYcbLHl5+Rx0FLVtaTy9jTNlYQgkzmZlQu7Y2ynvJ1kYxxpSVJZQwk5mp3YWrePzJdOgAmzbBrl3exmGMiRyWUMJIQYEmFC/bT/z8MfjXtTfGmCOxhBJGtmyBX3/1tv3Ez3p6GWPKyhJKGPFiDZSSHH881K9vCcUYEzhLKGEkM1PbTtq18zoS7TZsDfPGmLKwhBJGMjK0u26tWl5Hojp0gJUrdTp7Y4w5EksoYSRcGuT9UlJ05P769V5HYoyJBJZQwsQvv2ijfDi0n/hZw7wxpiwsoYQJf/fccEoorVvrEsQ2Yt4YE4iAEoqInCsia0Vkg4iMLubx6iIyw/f4UhFpVuTxpiKyW0RGBCfs6BNOPbz8qlfXpGIlFGNMII6YUEQkDpgMnAe0AS4XkTZFNrsO2OmcOwl4HHikyOMTgfcrHm70ysyExERo3NjrSA7lXxvFGGOOJD6AbboAG5xzmwBEZDrQH1hdaJv+wDjf3zOBSSIizjknIhcBm4HfghZ1FPKvgSLidSSH6tABXnxR17hv2NDraExA9u2DtWth61a95OTomUqTJnDCCdC8efh90UxUCCShHA9sLXQ7GzitpG2cc3kikgM0EJF9wCjgbKDE6i4RGQoMBWjatGnAwUeLvDydtv4vf/E6ksMVbpjv08fbWEwp9u6FDz6A11+Ht9+G3btL3rZFC/jzn/WSkmLJxQRNZTfKjwMed86V8u0G59xU51yqcy61YQyeBq9bB/v3h1f7iZ/19Apz+/fD449DUhJccgnMnQuDBsGMGfD555Cdrcll3Tr46COYMkUHOz36KHTqBF27wqJFXr8LEyUCKaF8CzQpdDvJd19x22SLSDxQB9iBlmQGiMg/gLpAgYjsc85NqnDkUeTLL/U6nMag+DVsqLUl/hhNmHAOpk+He+7R/ubnnAN33glnnaVd84pq2VIvZ50FN90EP/0Er70GDz8MZ54JF1wAjzyivTCMKadASijLgJYi0lxEqgEDgdlFtpkNDPH9PQCY79QZzrlmzrlmwBPAw5ZMDpeeDjVqQJuiXR3CROfOGqMJEz//DBddpCWROnVgzhy9nHNO8cmkOImJcMstWnL5+9/h44/1jObppzVZGVMOR0wozrk8YBgwB1gDvOacyxKR8SJyoW+zaWibyQbgDuCwrsWmZOnpWrUU6G9BqHXurG28v/7qdSSGzz6Djh3h/fdh4kRYvlwTSXnVqgWjR+t0CH36aEPepZfqSFtjyiignzDn3HvAe0XuG1vo733ApUfYx7hyxBf1Cgq0OunKK72OpGSdO+tJa2Ym9OjhdTQxbMoUuO027a21ZAmcemrw9n3MMdqYP3Ei3H23Jqp337UqMFMmNlLeYxs26Jl/p05eR1Iyf2xW7eUR52DcOK2iOu88PQMJZjLxq1IFRoyAxYu119gZZ8AXXwT/dUzUsoTiMf+PdOfO3sZRmsaN4dhjLaF4Ij8fhg2D+++Ha6+FN9/UdpPK1LWrloDq1NFG/A8/rNzXM1HDEorH0tN1ipO2bb2OpHTWMO+B/Hy46iptKL/rLnj22dA1tLVooUnlpJOgXz+YNSs0r2simiUUjy1fDsnJULWq15GUrnNn+Oor+M3mOwgN5+Dmm+GVV7Rr7yOPhH4A4rHHwsKF+uFfdhnMmxfa1zcRxxKKh5zThBLO1V1+nTtrBwIb4BgCzmmJ5N//hr/9TRvJvVK3Lrz3HrRqpV2VP/vMu1hM2LOE4qGNG3WapUhJKGDVXiHx97/DY49pF94HHvA6GqhXT0fgH3cc/PGPdlZhSmQJxUOR0CDv17gxNGpkCaXSvfCClkquuAL+7//CZ56tY4/VKq+EBE0q3xadLMMYSyieSk+HatXCv0Ee9HetUydLKJXqk0/g+uu1Z9V//qPdeMPJCSdo9deuXXDhhdagZg4TZt/Y2LJ8ObRvr0klEnTuDKtXw549XkcShTZtgosv1qnlZ84M314a7dvrHGIZGdoDraDA64hMGLGE4pFIapD38zfM+5crNkGSk6OTM+bnwzvvaJtFOOvXD/75Tx0TM2aM19GYMGIJxSObN8POnZGXUMCqvYKqoEDn3Vm3Dt54Q2cEjgTDh8PQodqB4PXXvY7GhAlLKB6JpAZ5v6Qknc7eEkoQTZjw+xxavXt7HU3gROCpp+D00+Gaa2DNGq8jMmHAEopH0tO1mrxdO68jCZyIjZgPqrlztcpo0CCdXiXSVKumpZOjjtL2n127vI7IeMwSikf8s5BXr+51JGVz2mm6XLFNZV9BW7bA5ZdrF7+pU8One3BZHX+8rg65YYOWVGwtlZhmCcUDubk6iWu3bl5HUnbdumm1/9KlXkcSwQ4c0PXc8/K0Yfuoo7yOqGJ69dKpYd58U5cjNjHLEooHMjJg3z7o3t3rSMqua1c9mf70U68jiWCjRsGyZfDcc5HTCH8kd9yh1V6jRtnZRgyzhOIB/49xJJZQjj5ahyJYQimnWbPgiSd0oaxLLvE6muARgWnTtOfGZZdpF0YTcyyheODTT3XQcePGXkdSPt26aRuQjWkroy1btJ2hc2f4xz+8jib46tXT9pTvvrP2lBhlCSXEnNNlJiKxdOLXrZt26MnK8jqSCJKbCwMHahZ+7bXI640RqC5dNFn+7386F5mJKZZQQmzrVp1XL9ITCli1V5mMGaNtC9Om6eJV0Wz4cJ3r6667dLliEzMsoYSY/0c4Ehvk/Vq00JmHLaEEaO5cPWu/8UYYMMDraCqfiE5uecwx2p5ifcxjhiWUEPv0U+0l2r6915GUn4iWUiyhBGDbNp1apV272OpS26ABvPyyLvoTiYM2TblYQgmxJUt0cGColgavLN266Vi2H37wOpIwVlCgM/L++qvO0FuzptcRhVbPnjB2rK7x8uKLXkdjQsASSgjt3q2L3UVy+4mf/z3YirCleOwx+PBD7SYcCYveVIYxY+DMM+GWW2D9eq+jMZXMEkoILVumM5RHQ0Lp1EmncrJqrxIsXaorL156Kdxwg9fReCcuDl56Sb8sAwfC/v1eR2QqkSWUEPL/+J5+urdxBEONGpCaagmlWL/8oj+exx8f2fN0BUtSks4KsHw53H2319GYSmQJJYQ+/VRrPurW9TqS4OjWDdLS7KTzEM5pb66tW+HVV6Pnw66oCy+EW2/Vjgnvvut1NKaSWEIJkdxcWLwYevTwOpLg6dFDk8nnn3sdSRiZNk0HLj74YHQURYPpH/+AlBQYMkQHY5moE1BCEZFzRWStiGwQkdHFPF5dRGb4Hl8qIs18958tIukistJ3fVZww48cn3+unX3OOcfrSIKnd2/trTZnjteRhImVK/UsvE8fHdRnDlWjhvZ227dP14DJy/M6IhNkR0woIhIHTAbOA9oAl4tImyKbXQfsdM6dBDwOPOK7/yfgAudce2AIELN9B+fM0fbJP/zB60iC5+ij9STcEgrahe/Pf9YqrpdegipW+C9Wq1bwzDOwaBHcf7/X0ZggC+Rb3wXY4Jzb5Jw7AEwH+hfZpj/wvO/vmcAfREScc186577z3Z8F1BSRKJ3EqHRz5+rU73XqeB1JcPXtq22t27d7HYnH/vIXWLtWB/M1auR1NOHtiivg2mvhoYdg3jyvozFBFEhCOR7YWuh2tu++YrdxzuUBOUCDItv8CVjunDusCVdEhopImoikbY/CX6afftLG62iq7vLzv6cPP/Q2Dk/99786eG/sWDgrZmt1y+app6BNGxg8GL7/3utoTJCEpFwuIm3RarAbi3vcOTfVOZfqnEtt2LBhKEIKqXnztPNP375eRxJ8nTrpLBsxW+21YoUO2uvVC+691+toIketWtp5Yfdu7WJt7SlRIZCE8i3QpNDtJN99xW4jIvFAHWCH73YS8BZwlXNuY0UDjkRz5+pSEampXkcSfHFx2gY9d24MLn+RkwN/+pO2m7z6qh4ME7g2bXSczqJFcM89XkdjgiCQhLIMaCkizUWkGjAQmF1km9nZ9NAgAAAQpklEQVRoozvAAGC+c86JSF3gXWC0c25JsIKOJM7p2XufPtH7e9O3r86BuHKl15GEkHNw9dWwebOeaR97rNcRRabBg7WE9+ijuia9iWhHTCi+NpFhwBxgDfCacy5LRMaLyIW+zaYBDURkA3AH4O9aPAw4CRgrIhm+yzFBfxdhLCtLF7CLxuouP387SkxVez32mC7n++ij0TW4yAsTJ+rCXFdfDevWeR2NqQBxYVZPkZqa6tLS0rwOI2j++U8YMQK++QaaNDny9pGqXTs9SY+JTjvz5sG558LFF2vpJNanVgmGb77RBrlGjXTQVu3aXkdkfEQk3TkXUIW9dZavZHPnQuvW0Z1MQEtgixfDnj1eR1LJNm7U8SannKKLSFkyCY6mTTU5r12r68cUFHgdkSkHSyiVaO9ebW+M5uouv7594cABWLjQ60gq0a+/Qv/+mkRmz7az6GA76yyd6+t//4P77vM6GlMOllAq0bvv6iwT/fp5HUnlO+MM/X194w2vI6kkBQV65vzVV3omHe3rwntl2DAd9Pjgg/D6615HY8rIEkolevllbVfo3dvrSCpfzZrapDBzpibRqHPPPXrmPHFidM2fE25E4OmndSrrIUPgiy+8jsiUgSWUSrJzJ7z3no7ZitbuwkUNHgy7dun7jirPPAOPPAI33aSTP5rKVb06vPWWno1dcIF2zTYRwRJKJXnjDW1TGDTI60hC56yz4JhjtGQWNd55R+fp6tdPpwuxRvjQOOYYeP99XffhvPPg55+9jsgEwBJKJXnlFWjZMjpHx5ckPl5LZO++q4sWRrz0dLjsMl3DY/p0fYMmdFq10mrGzZu1M0RU1qVGF0soleDbb7W306BBsXdCO2iQLroV8YOeV6/WsSYNG2opJSHB64hi0xlnwPPPwyefaHft3FyvIzKlsIRSCaZP15k5Yqm6y69LFzjxRC2hRaxNm+Dss7VE8uGHcNxxXkcU2wYOhMmT4e234aqrID/f64hMCSyhVIKXX9aqrpNP9jqS0BPRRDp/vk45E3Gys7UX1759mkxatvQ6IgM639cjj+jZ2o032sDHMGUJJcjWrIEvv9QeT7Fq8GAtoc2Y4XUkZeRPJjt26MRk7dp5HZEp7K67YMwYmDZNx6tYUgk7llCC7JlntKbkssu8jsQ7rVrBqafqsYiY2olNm7S+/vvvtXdRLPWmiCTjx2timTIFrrsugr5gscESShD9+CP8+986oDrWq91HjtSJY996y+tIAvDVV9Czpw6imT8funf3OiJTEhGYMAHGjdOVMgcPtob6MGIJJYiefFKr3keN8joS711yibYhPfxwmC+8lZYGZ56pP0oLF1rJJBKI6Fxfjz6q9aoXX6wrPxrPWUIJkpwcmDQJBgzQKp9YFxcHo0dre1LYrpPy1ltaMqlVS2fxbN/e64hMWYwYofWq77+vn+O3RReSNaFmCSVInn5aa0zuvtvrSMLH4MGQlKSllLDinC6Q9ac/QXKyrr9hZwGR6cYbdZzQ+vVw2mmQkeF1RDHNEkoQ7Nmjs26fey507Oh1NOGjWjVtS1m8WC9h4bffdGXAkSO1OLlggS7qZCLXeefpwEcRXT0zogdBRTZLKEHw7LOwfbtOSGsOdf31kJgIDz3kdSTo6PcuXeDFF7VRd/p0nSbZRL4OHXRm4o4dtWh84402VYsHLKFU0Lffavtgr17a69QcqlYtLQzMmePhWinOwXPPaV/mn37SZTTvuw+q2Nc/qhx3nJY4R42CqVOha1cdGGZCxv6jKsA5uOEGnbtq6lSvowlft9+uy4XffLOW5EJq61adKfjaazWhfPkl9OkT4iBMyMTHa7fid97RgaopKfD3v0NenteRxQRLKBXwn/9oB5MJE2yGjtJUrarz++XkaFIJSTfiggLN8m3bwscfa5/u+fOhceMQvLjxXL9+kJWl66ncc4+WVqzBvtJZQimnb77RM+9evXQWCFO6du3g/vu12qvSp2T5+GMdT3LjjXq9ciXcdptVccWaRo10CdHXX9d/2E6dtErhhx+8jixq2X9YORw4oB2FnNNSiv1OBWbECO3Z+Ze/wJYtlfACWVnaFbhXL20reeUV+OgjW/891g0YAGvXwvDhOrr+pJO0l8iuXV5HFnXsp7CMDhzQZRkWLNAF/Jo39zqiyBEfr1VfBQXQu3cQk0p6ug7Nb9dOG9wffFB/QC6/PPYWpDHFq1dP+/ZnZenSomPGwAknaOcMWw0yaCyhlIE/mfzvf5pMrr7a64giT6tWOiv8L79UMKkcOKB1Z717a7XWggUwdqzu8G9/s+7Apngnn6z/wMuW6Xdn/Hho2hRuukk7bJgKsYQSoD17Dk0m1m5SfqmpMG+eJpVevWDDhgCf6JzOvTVyJDRpogsvbdmi62R8/bU20jRoUImRm6iRmqrLiq5cqf/YL7ygbSynnabTXlg7S7mIC7OZ+1JTU11aWprXYRxi3jwYOlSXtp40SdsATMWlp+vCiPv3wwMPaLv5Ycu2HzgAn36q3elef10/hPh47cVz8826A2vEMhW1c6cOeJ06VavFqlTRSUMvvhj69tVunDFafSoi6c65gGZNtYRSim3bdG6u//5XS8r//rfOQWeCJztbF+N7+209aZzy+D5SJV2TyOLFWpW1e7cmkT/8Qc8mL7oI6tf3OnQTjZzThPL661qlunat3t+smZ68dOumlxhKMEFPKCJyLvAkEAc865ybUOTx6sALQGdgB3CZc26L77G7geuAfOA251ypc896nVDy87Vd99//1h8553Q9n7FjoUYNz8KKPvn52pVz/Xrcqiy2zF7BriUraJ23kmro+hYFLU6iSt+z4ZxztCH16KM9DtrEnI0b9Qdhzhxd3iAnR+9v0EAHTSYn6yzVrVtrkqlfP+oSTVATiojEAeuAs4FsYBlwuXNudaFtbgGSnXM3ichA4GLn3GUi0gZ4FegCNAbmASc750pcZi2UCaWgQBfo27BBJ5xdskRPjHfsgIYNYcgQreqyQYtlsG+fVh/88ot23f3xR718/70WR7Zu1USyefOhCyMdeyy5bZLJlI68sP50ZnzTlZzqjTj1VF3vqnt3/Z9t2lQnnTQm5AoKdCqXzz6DpUshMxNWrYK9e3/fpm5d7frZpIlekpJ0PMwxx+ilfn3dpm5dXeMhAgQ7oZwOjHPO9fXdvhvAOff3QtvM8W3zmYjEA9uAhsDowtsW3q6k16tIQtm+9mcyH36X/HwOXnJz4UAuHNgP+/bDb7t1wtldu/T3LrfQjAzHHavJI7kDdO5UTH1+ZSrtcyj82JH+Lnxd9FJQ8Pt14UvhA5afr9NU5OXpwfNf9u/Xy759etm7V3sq7NkDv/76+2X//uLfQ5UqOtdSUpL+o514oh7sE0+ENm30n63QW1m6VGsdliyB5ct/zz1VqugujjtO/zfr14c6dbRTl/9Stape4uP1fzYuTp8ncvjFr6S/jSmNFOST8MNGjt62jto/bODoH9aTsH0ztXZmU2vHVqrv+aXE5+bWSCC3Rm1ya9Qmr0YCedVqkV+tJvlVa5JftQb5VatTULUG+fHVcHFVKfBdXJU4CuLicVXicBKn11XicFIFqlTBSRWcCPivq1al69Tryv8ey5BQAvnJPB7YWuh2NnBaSds45/JEJAdo4Lv/8yLPPb6YgIcCQwGaNm0aSNzF+in9a/q8cFW5n8823yVcplr3ioj+GletqsWBqlWhevVDL7VqQe3amghq1/79Ureu9vmvV0+nGfafmSUmBpyhRXSmjK5d9fbevZpU1q/Xgs3mzdoJZ/t2reLOydFt9u4N89UhTRSKA072XQ5Xi99oyHaO4Uca8QP12HnwcvS+XdTe9ytHs4sEdlOTvdTkV2rxA9XZT3X2U5P9VOMAVck9eB1PiRU8xdpDTahAQimLUJ6Dl8g5NxWYClpCKe9+TrqwDT99voH4eA5eqlWLoE5ApZ0aB3I67f+78LX/UvgUvfApu/9v/6m8/3YYqVnz92qv0jinncJyc38vYOXnH1oYK1xgK/y84v42puKO8l2alfmZBcBe3+UQ/lqG/HykIP/368L383uthAicUNG3EaBAEsq3QJNCt5N89xW3TbavyqsO2jgfyHODpmpCdRJPO7Gydm/CnMjvBShjopegJaPwa4MJ5FR0GdBSRJqLSDVgIDC7yDazgSG+vwcA8502zswGBopIdRFpDrQEvghO6MYYY8LJEUsovjaRYcAcNCX+xzmXJSLjgTTn3GxgGvCiiGwAfkaTDr7tXgNWA3nAX0rr4WWMMSZy2cBGY4wxJSpLL6/wan01xhgTsSyhGGOMCYqwq/ISke3A1xXcTSLwUxDCiWR2DJQdBzsGfnYcyncMTnDONQxkw7BLKMEgImmB1vlFKzsGyo6DHQM/Ow6VfwysyssYY0xQWEIxxhgTFNGaUKZ6HUAYsGOg7DjYMfCz41DJxyAq21CMMcaEXrSWUIwxxoSYJRRjjDFBEVUJRUTOFZG1IrJBREZ7HU+oiEgTEVkgIqtFJEtEhvvury8iH4rIet91Pa9jrWwiEiciX4rIO77bzUVkqe87McM3wWlUE5G6IjJTRL4SkTUicnqsfRdE5Hbf/8IqEXlVRGrEwndBRP4jIj+KyKpC9xX72Yv6P9/xWCEinSr6+lGTUHxLFU8GzgPaAJf7liCOBXnAnc65NkBX4C++9z4a+Mg51xL4yHc72g0H1hS6/QjwuHPuJGAnEJqVhrz1JPCBc+4UoAN6PGLmuyAixwO3AanOuXbopLYDiY3vwn+Bc4vcV9Jnfx46A3xLdIHDKRV98ahJKOi69Rucc5uccweA6UB/j2MKCefc98655b6/f0V/QI5H3//zvs2eBy7yJsLQEJEkoB/wrO+2AGcBM32bxMIxqAP0RGcAxzl3wDn3CzH2XUBnUq/pW5+pFvA9MfBdcM4tQmd8L6ykz74/8IJTnwN1ReS4irx+NCWU4pYqPmy54WgnIs2AjsBSoJFz7nvfQ9uARh6FFSpPAHehi92BLkP9i3Muz3c7Fr4TzYHtwHO+qr9nReQoYui74Jz7FngM+AZNJDlAOrH3XfAr6bMP+m9mNCWUmCciCcAbwF+dc7sKP+Zb8Cxq+4iLyPnAj865dK9j8Vg80AmY4pzrCPxGkeqtGPgu1EPPvpsDjdE1eItWA8Wkyv7soymhhHS54XAjIlXRZPKyc+5N390/+IuwvusfvYovBLoDF4rIFrS68yy0LaGur9oDYuM7kQ1kO+eW+m7PRBNMLH0X+gCbnXPbnXO5wJvo9yPWvgt+JX32Qf/NjKaEEshSxVHJ11YwDVjjnJtY6KHCSzMPAf4X6thCxTl3t3MuyTnXDP3s5zvnBgML0GWpIcqPAYBzbhuwVURa+e76A7piasx8F9Cqrq4iUsv3v+E/BjH1XSikpM9+NnCVr7dXVyCnUNVYuUTVSHkR+SNaj+5fqvghj0MKCRHpASwGVvJ7+8E9aDvKa0BTdEmAPzvnijbYRR0R6QWMcM6dLyIt0BJLfeBL4Arn3H4v46tsIpKCdkyoBmwCrkFPHmPmuyAi9wOXoT0gvwSuR9sHovq7ICKvAr3Qaep/AO4DZlHMZ+9LtpPQ6sA9wDXOuQotlxtVCcUYY4x3oqnKyxhjjIcsoRhjjAkKSyjGGGOCwhKKMcaYoLCEYowxJigsoRhjjAkKSyjGGGOC4v8B1ZcTbJpyYBQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9e5b81a780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYHWWZ9/Hv3d1ZCAkkIWELS7OKDGuMEAQZJGwimBnhRRQlw5ZB2SagAoZhYBzBeZXNKEhYFDQgGFYB8VVkUwmSAMomISYhgJCEfYtZup/3j7tOztadPt2nllN9fp/rqut0naquejih7/511VPPYyEEREQkeS1ZN0BEpFmo4IqIpEQFV0QkJSq4IiIpUcEVEUmJCq6ISEpUcEWahJl90syez7odzUwFV5qemX3RzGab2ftm9qqZ/crM9qrzmAvNbL+42ljD+YKZbb2mfUIID4cQPtLH4y80sxVmNqri/Seic7f35bjNRgVXmpqZnQ5cClwAbABsBlwOTMyyXXEzs7YYDrMA+ELJMXcEhsRw3KahgitNy8zWBf4bOCmEcGsI4YMQwsoQwi9DCF+P9hlkZpea2d+j5VIzGxRtG2Vmd5nZ22b2ppk9bGYtZvZTvHD/MkrN3+ji3PuY2ctm9g0zWxIl638xs4PNbG50vG+W7L+bmT0SnetVM/uBmQ2Mtj0U7fbn6HyfLzn+mWb2GvDjwnvR92wVnWNstL6xmS01s33W8JH9FDi6ZH0ScH2fPvwmpYIrzWwPYDBw2xr2mQqMB3YBdgZ2A86Jtp0BvAyMxtPxN4EQQvgysAg4NIQwNITwf7s59obR+ccA5wJXAV8CPgZ8EvhPM9si2rcDmAKMito9AfgqfsK9o312js53U8nxRwKbA5NLTxxC+BtwJvAzMxsC/Bi4LoTwwBo+i1nAOmb2UTNrBY4EfraG/aWCCq40s/WA10MIq9awz1HAf4cQloQQlgLnA1+Otq0ENgI2j5Lxw6F3g5OsBL4dQlgJ/BwvppeFEN4LITwDPIsXeUIIc0IIs0IIq0IIC4ErgX/u4fidwH+FEJaHEJZVbgwhXAXMAx6N/jum1tDmQsrdH3gOeKWG75GICq40szeAUT1c39wYeLFk/cXoPYDv4gXr/5nZfDM7q7fnDyF0RF8XCuLiku3LgKEAZrZtdPniNTN7F7/mXHYDqwtLQwj/6GGfq4AdgGkhhOU1tPmnwBeBf0OXE3pNBVea2SPAcuBf1rDP3/E/yQs2i94jSqJnhBC2BD4LnG5mE6L94h6G7wrgr8A2IYR18MsX1sP3rLENZjYUv2F4DXCemY3sqREhhBfxm2cHA7fW0G4poYIrTSuE8A5+7fSH0Q2rIWY2wMw+bWaF6643AueY2eioS9S5RNctzewQM9vazAx4B7/O2hl932JgyxibOwx4F3jfzLYDvlKxvS/nuwyYHUI4Hrgb+FGN33ccsG8I4YNenq/pqeBKUwshXAScjt8IWwq8BJwM3B7t8j/AbOAvwFPA49F7ANsAvwXex9Py5SGE+6NtF+KF+m0z+1oMTf0a/qf8e/hlgJsqtp8HXBed74ieDmZmE4GDKBbu04GxZnZUT98bQvhbCGF2L9ouEdMA5CIi6VDCFRFJiQquiEhKVHBFRFKigisikpI4BrSQfmDUqFGhvb0962aI5MqcOXNeDyGMrnV/FVwBoL29ndmz1dNHpDfM7MWe9yrSJQURkZSo4Io0i2XL4LHHYN48UP/7TKjgivR3TzwB++4L66wDu+0G22wDI0fC6afDe+9l3bqmooIr0l91dsJZZ8G4cfDMM/D1r8Mtt8BVV8GnPw2XXgrbbQf339/zsSQWumkm0h91dsJXvgLTp8Nxx8F3vwsjRhS3H388nHYaHHMMHHww3HUXTJjQ/fEkFkq4Iv3RKad4sf3mNz3Rlhbbgt13hwcf9EsMhxwCDz1UvY/ESgVXpL+5/nq4/HI44wz4n/8BW8OwuaNHw333wWabwec/D0uWpNfOJqSCK9KfPP88fPWrsPfe8J3vrLnYFoweDb/4Bbz1Fkya5JcjJBEquCL9RUcHHHUUDB4MM2ZAWy9u0ey0E1xyCdx7L/zwh8m1scmp4Ir0F9Onw5w5XjA32aT333/iiXDAAXDOObB4cc/7S6+p4Ir0B6+/DlOnen/bI3qc8KFrZjBtmj8gcVZv58OUWqjgivQHU6f6QwzTptV23bY7227rN9t+8hN45JHYmidOBVck7+bOhauv9ptl229f//GmToUNN4Szz9YjwDFTwRXJu/PP9xtlU6fGc7yhQ73/7oMPwu9+F88xBVDBFcm3Z56BG2/0Bx3WXz++455wgt94O+ccpdwYqeCK5Nn553si/frX4z3u4MFebGfN8q5iEgsVXJG8euEFmDkTTj4Z1lsv/uMfcwxsuqk/QCGxUMEVyauLLoKBA+HUU5M5/sCBMGWKj7Hw6KPJnKPJqOCK5NHixd51a9Ik71GQlOOPh+HDfbQxqZsKrkgeTZsGK1Z4n9kkDRvm3c1uvdUvYUhdVHBF8mbZMrjiCpg40R9USNopp8CAAfD97yd/rn5OBVckb264Ad580wcQT8OGG8KRR/oljHffTeec/ZQKrkiehOCXE3bcEf75n9M77ymnwPvve9GVPlPBFcmThx+GP//ZC2A9Yyb01rhxsMceXuw1Xm6fqeCK5Mm0aT5dzlFHpX/uU0/1Kdb1IESfqeCK5MXf/w633QbHHgtDhqR//sMO8+u5V1yR/rn7CRVckby45hqf1eHf/z2b8w8Y4DMA3303vPhiNm3IORVckTxYtcpndNh/f59lNyuTJ/u146uuyq4NOaaCK5IH99wDL78MX/lKtu3YbDP4zGd8/N0VK7JtSw6p4IrkwY9+BBtvDIccknVLfO6zxYvh9tuzbknuqOCKNLqFC71nwPHH+3XUrB14IGy+uV/ikF5RwRVpdFdf7ddNjzsu65a41lYv/vfd593EpGYquCKNbOVKuPZa+PSn/fppozj2WC+8unnWKyq4Io3s7rvh1Vez6wrWnY03hkMPhR//WDfPekEFV6SRXXkljBnjCbfRTJ4MS5fq5lkvqOCKNKqFC+HXv/brpW1tWbem2gEH6OZZL6ngijSqws2yY4/NuiVd082zXlPBFWlEjXqzrJJunvWKCq5II7rrrsa8WVZJN896RQVXpBE18s2ySrp5VjMVXJFGM3++3yw74YTGvFlW6YADoL1dwzbWQAVXpNFceWXxhlQetLb6pY8HHoC//jXr1jQ0FVyRRrJ8ud8s++xn/ZJCXhx7rI/z8KMfZd2ShqaCK9JIbr0VXn89+2EYe2v99X1GiOuugw8/zLo1DUsFV6SRXH45bLUVTJiQdUt678QT4e234cYbs25Jw1LBFWkUTz4Jv/89nHQStOTwR3PvvWGHHeAHP/Dp3KVKDv9VRfqpH/zAJ4c85pisW9I3ZnDyyf6L449/zLo1DUkFV6QRvPEGzJgBX/4yDB+edWv67ktfgnXX9V8eUkUFV6QRXHst/OMffjkhz9Ze23sszJzpT8pJGRVckaytWgXTpsE++8COO2bdmvqddJJP53755Vm3pOGo4IpkbeZMeOklOP30rFsSj622gokT/ckzdREro4IrkqUQ4KKLYJttfPrx/uKMM/y69PXXZ92ShqKCK5KlP/wBZs+GKVPy2RWsO3vuCR//OFxyCXR2Zt2ahtGP/oVFcui734WRI+Hoo7NuSbzM/BLJ3Llw551Zt6ZhqOCKZOWpp7wYnXqq393vbw4/HLbcEi68UA9CRFRwRbJy4YUwdCicckrWLUlGWxuceSb86U8+DY+o4IpkYt48uOkmH6Rm5MisW5OcSZN8VogLLsi6JQ1BBVckCxde6MMZTpmSdUuSNWgQfO1rcP/9Pk5Ek1PBFUnb3Lk+jOGJJ8JGG2XdmuRNngwbbADnnNP013JVcEXSdt55nvzOPjvrlqRj7bVh6lR48MGmv5argiuSpqeegp//HE47zVNfs5g82ad7nzq1qVOuCq5Ims48E9ZZx69rNpNBg+Dcc73Hwi23ZN2azKjgiqTlV7/y5T//s3/3TOjOpEk+OM/Xv+4jozUhFVyRNKxc6U9ebb11/+1325O2Nrj0Uli40B/5bUIquCJp+OEPfQrxiy6CgQOzbk129t3XRxL79rfhlVeybk3qVHBFkrZokXeJOuggOPTQrFuTvYsv9gFtTjqp6W6gqeCKJCmE4pTnP/qRD+rS7LbcEs4/H+64w6eFbyIquCJJuvFGuOce/xN6882zbk3jmDIFdt3VJ518882sW5MaFVyRpCxY4Ol2jz28sEhRW5vP4/bGG3D88U1zaUEFVyQJq1bBUUf51zfcAK2t2banEe2yi48pcdttcOWVWbcmFSq4Ikk4+2x45BEvJO3tWbemcU2ZAgce6K9z5mTdmsSp4IrE7brr4Hvf88sJRx6ZdWsaW0uLz3u2/vreXayfT62ugisSp4cf9nEDJkyAyy7LujX5sP76PvPF22970X3//axblBgVXJG4/OlPPvPuFlvAzTf7eLdSm513hhkz/LLCZz8Ly5Zl3aJEqOCKxGH2bL8WOXq0D0HYjGMl1GviRL8c88AD8LnPwYcfZt2i2KngitTrV7+CffaBddf1YjtmTNYtyq8vfQmuugp+/Wu/LPP661m3KFYquCJ9FYIPxnLoobDttt4rQT0S6nfccTBzJjz5JIwf76/9hAquSF8sWeJ/Ak+ZAocc4rMZNMN0OWn53Ofgd7/za7njx8O0aT7+Qs6p4Ir0RkeH/8m73XZw773eE+G222DYsKxb1v/ssYen2333hVNPhb32gj//OetW1UUFV6QWHR3e82Cnnbzb1047+Q//qadqQJokjR4Nd98NP/kJvPCCj7/wxS/Cc89l3bI+UcEVWZMFC3zgmS23hM9/3t+7+Waf9vujH822bc3CzGeLeP55n6Lozjth++1h//393+KDD7JuYc0sNMmgEbJm48aNC7Nnz866Gdl7913vT3v//X7J4PHH/f0JE+CrX/XrthoXIVtLl8L06b4sWgRrrQUHHAD77Qd77+3FuK0tlaaY2ZwQwria91fBFWiigrt8uQ8HuHgxvPaa/8AuWODp6ZlnYO5c36+tDXbbDf71X+Gww/xhBmksHR3+ZN/MmT4E5oIF/v5aa/kln49+FLbZxofF3HRTnyV59GjvvhfTL00VXOmTsoJ7++1w113Jn7S7//cK75e+drd0dPjd61Wr/OuVK31ZvtyXZct8ee89X7p6gqmtzeca++hHYexY+PjH4ROf0I2wvFmwAP74R3jsMZ+O/rnnuh6bwcz/bYcNg7XXhiFDYPBgn1l44EB/QrCtzZfWVh/vYe+9/S+cqkP1ruCmk7slX+bN8z+n09DdDafC+6WvXS2FH4jCD8eAAb4MGgRDh3raWWut4g/YiBH+FNj663vi2Wwz2HhjXSboD7bYwpfCsJjg13cXLfL50xYv9gcp3nwT3nnHfwF/8IE/0faPf/gv6Pfe81/Yq1b50tnpS0yDxyvhCtBElxREYtTbhKteCiIiKVHBFRFJiS4pCABmthR4seStUUAjjxzSyO1r5LZBY7evkdsG1e3bPIQwutZvVsGVLpnZ7N5cm0pbI7evkdsGjd2+Rm4b1N8+XVIQEUmJCq6ISEpUcKU707NuQA8auX2N3DZo7PY1ctugzvbpGq6ISEqUcEVEUqKCKyKSEhVcqWJmB5nZ82Y2z8zOyrgtm5rZ/Wb2rJk9Y2anRe+PNLPfmNkL0euIDNvYamZPmNld0foWZvZo9PndZGYDM2zbcDObaWZ/NbPnzGyPBvvspkT/rk+b2Y1mNjjLz8/MrjWzJWb2dMl7XX5e5r4ftfMvZja2p+Or4EoZM2sFfgh8Gtge+IKZbZ9hk1YBZ4QQtgfGAydF7TkLuC+EsA1wX7SeldOA0ikI/he4JISwNfAWcFwmrXKXAfeGELYDdsbb2RCfnZmNAU4FxoUQdgBagSPJ9vP7CXBQxXvdfV6fBraJlsnAFT0ePYSgRcvqBdgD+HXJ+tnA2Vm3q6Q9dwD7A88DG0XvbQQ8n1F7Nol+CPcF7gIMfxKpravPM+W2rQssILo5XvJ+o3x2Y4CXgJH4yIV3AQdm/fkB7cDTPX1ewJXAF7rar7tFCVcqFX4ICl6O3sucmbUDuwKPAhuEEAqDnb4GbJBRsy4FvgEUppRdD3g7hLAqWs/y89sCWAr8OLrkcbWZrU2DfHYhhFeA7wGLgFeBd4A5NM7nV9Dd59XrnxUVXMkFMxsK3AL8Rwjh3dJtweNF6v0bzewQYEkIYU7a565RGzAWuCKEsCvwARWXD7L67ACia6ET8V8MGwNrU/3nfEOp9/NSwZVKrwCblqxvEr2XGTMbgBfbGSGEW6O3F5vZRtH2jYAlGTRtT+CzZrYQ+Dl+WeEyYLiZFQb3z/Lzexl4OYTwaLQ+Ey/AjfDZAewHLAghLA0hrARuxT/TRvn8Crr7vHr9s6KCK5UeA7aJ7hQPxG9i3JlVY8zMgGuA50IIF5dsuhOYFH09Cb+2m6oQwtkhhE1CCO345/S7EMJRwP3A4Vm2LWrfa8BLZvaR6K0JwLM0wGcXWQSMN7Mh0b9zoX0N8fmV6O7zuhM4OuqtMB54p+TSQ9eyuFiupbEX4GBgLvA3YGrGbdkL/xPuL8CT0XIwfq30PuAF4LfAyIzbuQ9wV/T1lsCfgHnAL4BBGbZrF2B29PndDoxopM8OOB/4K/A08FNgUJafH3Ajfj15Jf4XwnHdfV74DdIfRj8nT+G9LdZ4/Loe7TWzg/A/oVqBq0MI3+nzwURE+rk+F9yov+ZcvIvOy/ifol8IITwbX/NERPqPembt3Q2YF0KYD2BmP8fvOHZbcEeNGhXa29vrOKXU6+23fWLSTTctf3/OnDmvh2jk+v1b/k9ffwtXrLdUrHazveJ9KxynpaX8uNF6cXthRt+K47S0rv569b6FWXkrZwNuLXyvv4aWinO3lrchrH6/u/XoNfq+0Gpdb49eO9sqv6/wSvl6dJrOiu2V64X9itvLj9PZVr696rVwnsr92kI3xw3l26NXovdpC1j0tbV5z7nW1ug1Wh8wwHuADWjtAGBgm78Oai28+va12lYCMDh6Xbt1hb/fGq23LQdgSIu/P6z1HwAMjV7XaVkWbV8erfv7w1a/+nGGWYjW/UMY2jIYgJYNX+hmiuna1VNwu+qDtnvlTmY2GX8Kg8022wzNDJutb3wDpk2Dyn8GM3ux6+8QkbjUU3BrEkKYTjSG5Lhx4zQWZMZWroQBAxI6eOHyVCE9huhZgCiBhs4o6bRUbO8sT6iFy1zWGW0vpMhovZAqrfCoQUvFcejAbysUE5t1eFpanXQLOjrLVi3quBMof7+QdAttKjwzZFSuF6xuXLSdiu3R1qh7f2fVT2Jhz/LvbInWO7tZr1T4RDqj/Vqi/Tq73LuLdnXTnuJxuz5v6PJr/66Oqr1jLkP1Hi5KunQWnr2IEnKdh4X6uoU1XH9N6dmKFTAws6FURJpbPb8LVvfXxAvtkcAXY2mVJObDD2HttRM+SUMk3UKOUtItnl1Jt1cqkm4cCbfPTQohrDKzk4Ff4/9XXxtCeCaGNkmC3n0Xhg7NuhUizamu3wEhhHuAe2Jqi6TgjTdgvfVSOlmWSbfsei4o6ZaePd2kW3njJrdJNwZ6tLfJvPIKbLhh1q0QaU6J91KQxtHRAQsXwuGH97hrvLJIul32XAAl3dKzK+mmTQm3iTz9NKxaBTvskHVLRJpTA9R8ScuDD/rr3ntn1IAUk+6a++iCkm7p2ZV006KE20R+9jPYfvvqx3pFJB1KuE1i1ix47DF/rDdzKSTd2p5GAyXd0rMnk3Rr6aNbvt5/k64SbhNYsQK+8hXYYAM4+uisWyPSvJRw+7kQYOpUePJJuP12WKeWx2VWJ8+Eh75IMun2atwFUNItPXu8Sbc3T6OVr/e/pKuE24+F4KODfe97nnAnTsy6RSLNTQm3n3rnHTj1VLj+ejjpJPj+9/twkDwn3T6NMAZKuqVnjyvp9n7chfL1/pN0lXD7odtv994IP/sZnHuu3yhr0b+0SOaUcPuRRx6BCy+EX/4SdtrJC+/HP96HAxVmUFidNPOXdOsbSxeUdEvPXm/S7fsIY+Xr+U+6yj0519EBt90Ge+4Jn/gE/P73cMEFPqNDn4qtiCRGCTennnsObroJZsyAefOgvd2v0x5zTIzDL+Y46cYzawQo6ZaevW9JN46xdMvX85t0VXBzZO5cL7I33+zjIpj5Y7rf/jZ87nPQpn9NkYamH9EG9v778PDD8Nvfwm9+A0895UV2r738Rthhh8FGG8V/3kJyLCTJXCbdWOdHAyXd0rP3LunGOWtE+Xr+kq4KbgNZuRIefdQL7H33+eO4q1b5HGR77gmXXupDK44Zk3VLRaQvVHAztGSJF9jC8sgj8MEHHug+9jE44wzYbz8vtmutlX77cp10E5kJGJR0S8/eZEk3Bo3Xon5q+XJ/vHbWLC+us2bBggW+rbUVdtwRJk2CCRNgn31g5MhMmysiCVDBTcD778Nf/uIF9okn/PUvf/FBZMAvCYwf74/b7r67p9nEZ9LtjYokm8ukm8hMwKCkW3H81a1R0q1F47QkpxYvLi+sTzwBL7xQrA0jR8Kuu8Jpp3lx3X132GSTbNssItlQwa3RsmXw7LPeU6B0ee214j7t7V5cjzoKdtnFv95kk2Igy43CtdDO/CbdRGYCLj2Okm758Ve3Rkl3TbJvQYPp6ID586sL67x5q38WGTzYxyo48MBiYd15Zxg+PNu2i0hja+qCu2RJdWF95hn48EPfbgZbbeU3tI480l933BG23ro6wPRLOU66icwEDEq6fUi69YylW37M/Cfdpii4H37ohbSyuC5ZUtxn9GgvpiecUCys//RPDXYzS0Ryrd8V3Dff9BtXjz9efJ07txiGBg/2acI/85liYd1xR59+RpxZefrLZdJNYCZg3x6dU0m3bL1SadJNYybg0hZXrzdO0s1twQ0B/v736uK6aFFxn003hbFjyy8HbLVVk1wOEJGGk5uC29nplwUeesjHF3joIXj11eL2bbeFPfbw2Q123dWXUaOya2+utVSkqjwm3QRmAgYl3b4k3SRmAl7TORs56TZswV25EubMKRbXP/wB3nrLt40Z409jjR/vCXbnnWHYsEybKyLSo4YquCH4ANrXX+9DEL77rr+/7bY+/OAnP+nDEba357Bva54UPtw8J90kZgIGJd0+JN2+9NEta1M3bclj0u3xDGa2KXA9sAHe5ukhhMvMbCRwE9AOLASOCCG81ZdGzJ8PP/2pF9r5871nwGGHwaGH+lCEG27Yl6OKiDSWWkr6KuCMEMLjZjYMmGNmvwH+DbgvhPAdMzsLOAs4s7cNuOgi+NrXPJBMmADnnedpVt2xMlSR5nKZdJOYCRiUdJV069LjkUMIrwKvRl+/Z2bPAWOAicA+0W7XAQ/Qy4I7fboX28MOg0su8V4FIiL9Va9KuZm1A7sCjwIbRMUY4DX8kkPNFiyAE0+EcePghht8kG1pDKv74SrpKukq6caqpeddnJkNBW4B/iOE8G7ptuD/N3b5U2Nmk81stpnNXrp06er3N9kEPvUp7z97zz19a7yISJ7UVMLNbABebGeEEG6N3l5sZhuFEF41s42AJV19bwhhOjAdYNy4cauL8oABcPvtsP/+fknhwAN9AO6JE/1pMMlQlOoKaU1JFyXdOpJuEjMBl++Xn6TbY8I1//vyGuC5EMLFJZvuBCZFX08C7ujtyYcNg3vvhTPP9LENjjzSeyRMnuz9bjt7+uRFRHKkltK9J/Bl4CkzezJ675vAd4Cbzew44EXgiL40YPhwuOAC+Na34IEH4LrrYMYMuOoqGDHCu4Xtvbf3wR071pOxJKyQQKN0pqSLkm4dSTeJsXTL2tRNW+JPuvWrpZfC76n+/AsmxNWQ1lbvFjZhAlx+OdxxB9x/vz9p9stf+j5Dhvjju4UHIHbbTd3HRCQ/GupJs4KhQ33WhKOO8vXXXvMn0ArjKJx/fvFhoo98xMdNGDu2OIaCJmCsU2VCVdJV0i1s7UPSTXLWiLI2ddOWRkq6DVlwK224IRx+uC8Ab78Nf/wjPPaY93L4/e/hxhuL+2++ebEAjx3rszJsvLEeBxaRbOWi4FYaPhwOPtiXgtdf9+JbWB5/3HtBFILMiBHFIRp32slfd9hBg950qaWQqqLf7Uq6SrrlR+1V0k1jfrSyNnXTlnqTbhxyWXC7MmqUdzHbf//ie++9B3/+sy9PPeVTlV9/vb9f0N5ePhD5jjv6YDlNfXOuUFxQ4VXhjaPw9v7hiOKZG6nw1q/fFNyuDBvmvRz22qv4Xgjw4ovV0+386lewKvqHGjgQttuuuhDncgZeEWkY/brgdsXMU217u49GVrB8OTz/vKfgQhF+8EHvolYwfLhfhqi8NLHOOmn/VySr+Ghv4R0lXSXdepJu3x8DLp45+6Qbh6YruN0ZNMgL6E47lb//1lvw9NPlaXjGjOJYveCz+BamSy+8brih0rCIlFPB7cGIEd7v95OfLL4XArz0kqfhJ5/05fHHYebM4j7rr18swLvsAh/7mBfmXBThKD0V0pSSLkq6dSXd+ge8KZ4530lXBbcPzGCzzXw55JDi+++840X4iSe8CD/xBFx8sU8XBN4/ePfdfRk/3h/cGDEim/8GEUmfCm6M1l23Og2vWAHPPguzZ8Ojj8KsWT5+RCFAbbutF9/dd4dPfMIvabS0dH381BRSnZKukq6SbqxUcBM2cGDxssLxx/t7777rBXjWLC/C997r3dXAu7d96lOw337+mPOWW+bkMoSI9EgFNwPrrAP77usLFLuqPfQQ3HefL7/4hW9rby+OMXHQQSldgqis8Eq6Srp1JN0kp2D37flJuiq4DaC0q9rRR/vP9vPPF4vvLbfANdf4wxgHHABHHOHjBq+7btYtF5HeUMFtQGb+4MV228FJJ0FHh48bccstPn383Xf7pYqDDvLie9hhMQ/aHqWdynSkpIuSbh+SbpLT9RTPnI+km/XtGalBa6vfWPvud2HhQnjkES/Ec+bAl77kg/V861vwxhtZt1Q4Q9TjAAAOdElEQVRE1kQFN2fMvPhefDEsWgS//a338T33XJ/1+OST/XpwXVpafGmNFrPypbUVWlsxM38qraWwtEZLtG4tvkTrxf2j4xeOF62v3r66HRXHKXwGLVZMkP5G2fZU7jKGUJ6kQ2cxzeJJd3USL93eGXxZfZjgabez05fCcaP14vZoqTxOZ0e0+Prq/Ts6fCkcr7B0dPoSHd86A9ZZcv7C9mh/6+z0tNsRoKOr9Wjp6IyWgEXbyrZ3BlpW+VL8nsJCtETrnZ7sWzoCLSXbK9cL+xW3+1I4TssqT7HF41cshfNU7rfKfKk4bhxUcHOspcVvpt1zT3GKounTYfvt4dJL/edNRBqHCm4/scMOcO218MILsM8+MGWK9+t95pneHyu0GKHF8p10C8dMmpJuj0m3NOXmOenGQQW3n9l8c7jrLrjhBpg/34vuH/6QdatEBNRLoV8ygy98wYel3G8/70p2551++aEmLeV3sK3y93Ieei+k2XOh9PjqvRAp9l7o2whjftTS9ex7L9RPCbcf23RTf5hiyy19eqJXXsm6RSLNTQm3n9tgA7jtNh+j4YQTvA9vj5c2W8t/D+cy6WbRR7f0+Eq6kU76Nu5C6Z6NknTrp4TbBLbeGi64wGe1uP/+rFsj0rxUcJvEiSf68JCXX17Dzqt7IUS9E6LeA3nqvVCyUt1HV70X0u29UEMf3Tz0XoiDCm6TGDwYJk2CO+6ADz7IujUizUkFt4kccIBPlDlr1pr3Cy0tfo0wx0m3pqfRlHRTSbq9eRqtkZNuHFRwm8gee/jrn/6UbTtEmpV6KTSRddf1Xgt/+1sPO7ZGd8oLd5wLd8rz1HuhkUYYKz1+E/ZeSGYm4OrvTKv3Qj2UcJtMe7sPeiMi6VPCbTKjRsGrr655n9UpqbJvZZ6SbiOOpVt6/CZKusnMj1a6Z36SrhJukxkxAt56K+tWiDSnmhOumbUCs4FXQgiHmNkWwM+B9YA5wJdDCCuSaabEZcgQWLash50qr+HmMek28qwRpcdvgqSb7EzApXs2ftLtTcI9DXiuZP1/gUtCCFsDbwHHxdguScjgwTUUXBFJRE0J18w2AT4DfBs43bwT5L7AF6NdrgPOA65IoI0SowEDYOXKNe8TovRUzA35S7q5mB+t9Pj9OOkmMRMwpJ9041Brwr0U+AbFT2Q94O0QQuGBt5eBMV19o5lNNrPZZjZ76dKldTVW6tfW5g8/iEj6eky4ZnYIsCSEMMfM9untCUII04HpAOPGjUs4LkhPWlqKgbA7obX4Ox7ymXRzNRNw6fH7Y9JNYCZgSD/pxqGWSwp7Ap81s4OBwcA6wGXAcDNri1LuJoBGW82BWgquiCSjx4IbQjgbOBsgSrhfCyEcZWa/AA7HeypMAu5IsJ0SE7Oew1poLfzGz3HSjWPWiJLtSrpRM5o46cahnmOdid9Am4df070mniZJkmopuCKSjF49aRZCeAB4IPp6PrBb/E2SJNU0QNbqfrgFOUy6cc6PVrJdSTdqRi+Sbpzzo/l2KrZHWxNOunHQk2YiIinRWApNppaEW90PtyA/STeRmYBLtivpRs2oIekmMROwb6die7S1gZOuEq6ISEqUcKVKIeEW5DHpJjITsL9Rtl1JN2qGkm5NlHBFRFKihCtVOtui3/AVjwDnKunGMWsEKOnGkXRjmDWifL0g3aQbByVcEZGUKOFKlcI13EIiyGXSjXN+NFDSVdKNhRKuiEhKlHClSnEsBZfPpJvATMCgpNuXpJvATMDl6wXJJt04KOGKiKRECVeqhNWBJM9JN4GZgEuPo6Qbbe856SYxE7Dvn27SjYMSrohISpRwpUrxGm7Xv+HzkHQTmQkYlHT7kHSTmAnYt6eddOunhCsikhIlXKlSTIf5TbpJzAQMSrp9S7rxzwRc2qb0km79lHBFRFKihCtVOlsr+x/mMOkmMBNwaZuUdHuTdBOYHw0ySLr1U8IVEUmJEq5UKfTDrR4tSUlXSbcPSTfJmYAhV0lXCVdEJCVKuFKl8hpuLpNuEjMBg5JuX5JuAjMBQ/pJNw5KuCIiKVHClSrdXcNV0lXS7UvSTWJ+NEg/6cZBCVdEJCVKuFIlVPwazmPSTWImYD9ndA4lXX+tIekmORMw5CvpKuGKiKRECVeqhNau389T0k1kJmBQ0u1L0k1ifjRIPenGQQlXRCQlSrhSpTMKEt39NlbSRUm3N0k3yZmAIbWkGwclXBGRlNSUcM1sOHA1sAMeQI4FngduAtqBhcARIYS3EmmlpKow40NnlDVzmXQTmQkYlHR7n3QTmQkYMki69as14V4G3BtC2A7YGXgOOAu4L4SwDXBftC4iIt3oMeGa2brA3sC/AYQQVgArzGwisE+023XAA8CZSTRS0rU6eZLfpJvMTMCgpNv7pJvETMCQRdKtXy0JdwtgKfBjM3vCzK42s7WBDUIIr0b7vAZs0NU3m9lkM5ttZrOXLl0aT6tFRHKolmu4bcBY4JQQwqNmdhkVlw9CCMHMuvxVG0KYDkwHGDduXMK/jiUOxbEUXB6TbhIzAYOSbp+SbgIzAfv+URNSSrpxqOVILwMvhxAejdZn4gV4sZltBBC9LomtVSIi/VCPCTeE8JqZvWRmHwkhPA9MAJ6NlknAd6LXOxJtqaSm8kmzPCbdJGYCBiXdPiXdRGYChrSTbhxqffDhFGCGmQ0E5gPH4P+332xmxwEvAkfE1ioRkX6opoIbQngSGNfFpgnxNkcaQfdjKbg8JN1E5kcDJd2C3iTdRGYChrSTbhz0pJmISEo0loJUqRwPt1Iekm6ss0aAkm53aki6icwEXHqclJJuHJRwRURSooQrVbpLlFX7Ra+NmXR7fhoNlHRjs4akm8hMwJB60o2DEq6ISEqUcKVK1ZNmuUy6PT+NVraupBuPrpJuAjMB+/bonGkl3Rgo4YqIpEQJV6qEtvJkmcekm8RMwL6/km5NSpNuAjMBQz6TrhKuiEhKlHClSvFJs/wm3VhmjQAl3XqF0LtxF6Bhk24clHBFRFKihCtVqmftVdJV0q1DEjMBQ+pJNw5KuCIiKVHClSqhtTJhrt4SvTZ+0k1kJmBQ0q2Hkq4SrohIWpRwpUqxH67LY9JNYiZgUNKNRV6TbgyUcEVEUqKEK1UqE26Bkq6SbqyaMOkq4YqIpEQJV6pFvRS6yzb5SLr1zxpRvl1JNzF5SboxUMIVEUmJEq5Uq7iGm8+kG9/8aOXblXQT0+BJNw4quFLFurmkkK/C2/ehHUu/S4VXhbf42ddPlxRERFKihCtVrK2YNSGfSTfJKdjLtyvpJqbhkm79lHBFRFKihCtVWlv9N3zx93r+km6yU7AXW6ek20xJt35KuCIiKVHClSqtbeUpK49JN9kp2EvPrqTr602QdGOghCsikpKaEq6ZTQGOx3/BPwUcA2wE/BxYD5gDfDmEsCKhdkqKBgwoxLvy/z3ylHSTnIK99LuUdJso6cagx4RrZmOAU4FxIYQd8H/ZI4H/BS4JIWwNvAUcF1+zRET6n1qv4bYBa5nZSmAI8CqwL/DFaPt1wHnAFXE3UNI3oLWy36GSrpKukm4cejxSCOEV4HvAIrzQvoNfQng7hFD43/BlYExX329mk81stpnNXrp0aTytFhHJoR4TrpmNACYCWwBvA78ADqr1BCGE6cB0gHHjxiX4607iMrCtuydrlHSVdPuQdJNMud7I4rnKzh1v0o1DLVl5P2BBCGFpCGElcCuwJzDczAo/gZsAr8TWKhGRfqiWa7iLgPFmNgRYBkwAZgP3A4fjPRUmAXck1UhJ16Cqa7iVGj/p9mbchdJ1Jd0Ekm4a13NLj59U0o1BLddwHwVmAo/jXcJa8EsEZwKnm9k8vGvYNbG1SkSkH6qpl0II4b+A/6p4ez6wW+wtkswNau0hUq7WuEk3yYkpQUm3V0k3zZ4LpcePO+nGQE+aiYikRGMpSJW12lb28jsaL+mmMQU7KOnWlHSz6KNbevyYkm4clHBFRFKihCtVBvc64RY0TtLtywhjvl1J188bnSeGpJvp02ilx68z6cZBCVdEJCVKuFJl7dZ6B31T0lXSLabZhhh3ofT4fU26MVDCFRFJiRKuVFmrta/XcCsp6SrpFvvh5j7pxkAJV0QkJUq4UmXttuUxH1FJt6mTbiOOpVt6/FqTbgyUcEVEUqKEK1WGtES9FGL/vyO9pJvETMC+XUnXzxudp5ak28izRpQev4ekGwclXBGRlCjhSpVhrf8ofyOHSbeesXTLj9ndOZV0/bzRedaQdHMxP1rp8btJunFQwhURSYkSrlQZWplwC3KVdOufNaL8mN2dU0nXzxudp4ukm6uZgEuPX5l0Y6CEKyKSEiVcqbJOy7I175CLpFv/rBG+rqQLdSbdJGYC9gORqMqkGwMlXBGRlCjhSpUhLTU+adbASbeep9HK9lPS9f3qSbpJzAQM6SfdGCjhioikRAlXqqzT0k0vhe40YNJNYiZgX1fShd4l3URmAi7ZnlrSjYESrohISpRwpcqw3ibcAiXdLs6rpJvITMD+Rtn2PCRdJVwRkZQo4UqVYS11zvigpNvFeZs46SYwE7Cv5i/pKuGKiKRECVeqDLNCwshv0k1y1oiy/ZR0fb81Jd0EZgKGfCZdJVwRkZQo4UqVYS3R/xadUZxS0u2Wkm5kjUk3gZmAIZdJVwlXRCQlSrhSZWjL4OirqD9uDpNunPOjgZJufUk3gZmAS4+To6SrhNtkxo+H007LuhUizclCitXezJYCL6Z2QumNzUMIo7NuhEh/lmrBFRFpZrqkICKSEhVcEZGUqOCKiKREBVdEJCUquCIiKVHBFRFJiQquiEhKVHBFRFKigisikhIVXBGRlKjgioikRAVXRCQlKrgiIilRwRURSYkKrohISlRwRURSooIrIpISFVwRkZSo4IqIpEQFV0QkJSq4IiIpUcEVEUnJ/weCS3jrpPMW6wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9e597f0160>"
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcVNWZ//HPQzeL7CDNvrQsxhhBxY6CGsaA4hINGXGUDI6MG+OKg8aoITFxMnH5uYcEFZcEE1xxDW4/g7hFQGkRFZVFQARlFUQQge4+88dzO91CA71U31vL9/163Vd13bpd9XQ1/eXUueeeYyEERESk/jVIugARkVyhwBURiYkCV0QkJgpcEZGYKHBFRGKiwBURiYkCVyRLmNkPzGx+0nXIrilwRSJm9p9m9p6ZfW1mK83sDjNrHT12p5ltirZtZra90v3nYqgtmFnv3R0TQngthPCdOrzGCDObZWabzWx19PUFZmbR42ZmN5jZumi7ofwxqR4FrghgZpcBNwCXA62AAUAP4EUzaxRCOC+E0DyE0By4Fni4/H4I4fjkKndmll/H778MuB24EegIdADOA44AGkWHjQZ+AhwI9ANOAv6rLq+baxS4kvPMrCVwDXBxCOH5EML2EMJS4FSgEDi9Fs95lJktN7OfR63Fz83sJ2Z2gpktMLMvzOwXlY4/1MxmmNmG6Ng/mFmj6LFXo8PmRi3q0yo9/xVmthL4U/m+6Ht6Ra/RP7rf2czWmNlRVdTaCvgf4IIQwpQQwlfBzQkhjAwhbI0OHQXcHEJYHkJYAdwM/GdN35tcpsAVgcOBJsDjlXeGEDYBzwLH1PJ5O0bP2wW4GrgbD+9DgB8AvzKzfaJjS4GxQDtgIDAEuCCqY1B0zIFRi/rhSs/fFm+Jj96h9o+BK4C/mllT4E/ApBDCy1XUORBoDDy1h5/ne8DcSvfnRvukmhS4Ih5ya0MIJVU89nn0eG1sB34XQtgOPBQ9z+1RC3Ie8AH+8ZwQQnEIYWYIoSRqXd8F/Msenr8M+HUIYWsIYcuOD4YQ7gYWAbOATsC4XTzPTj+/mb0Rtba3mFl54DcHvqz0fV8CzdWPW30KXBFYC7TbRT9op+jx2lgXQiiNvi4PxFWVHt+Chxhmtq+ZTY1O1m3E+4n3FPRrQgjf7OGYu4EDgPGVugZ2qpMdfv4QwuEhhNbRY+U5sQloWen7WgKbgmbAqjYFrgjMALYCJ1feaWbNgeOBaTHUcAfwEdAnhNAS+AWwp5bjboMuqv824F7gN2bWdheHlv/8w/bwevOIWuSRA6N9Uk0KXMl5IYQv8ZNm483sODNraGaFwCPAcuAvMZTRAtgIbDKz/YDzd3h8FdCzhs95OzA7hHAO8AxwZ1UHhRA24D//BDM7xcxamFkDMzsIaFbp0PuBS82si5l1Bi4D/lzDmnJanYaSiGSLEML/M7N1wE1ALzz8ngRG7uajeCr9DJgI/ByYAzwMDK70+G+ASWa2F36CbPXunszMhgHHAX2jXZcC75jZyBDC5B2Pj37+FdHr3w9sBhbjJ97eiA67Cw/996L790T7pJpM3S8iIvFQl4KISEwUuCIiMVHgiojERIErIhITjVIQANq1axcKCwuTLkMkoxQXF68NIRRU93gFrgBQWFjI7Nmzky5DJKOY2Sc1OV5dCiIiMVHgiuSKLVvgrbdg0SLQ+PtEKHBFst2cOTB4MLRsCYceCn36QNu2cOml8NVXSVeXUxS4ItmqrAyuvBKKimDePLj8cnjsMbj7bjj+eLjtNthvP5g+PelKc4ZOmolko7IyOP98mDgRzj4bbrwR2rSpePycc+CSS+DMM+GEE2DqVBgyJLl6c4RauCLZ6OKLPWx/8Qtv0VYO23KHHQavvOJdDCeeCK++uvMxklIKXJFsc//9MGECXHYZ/O//wu4WZCgogGnToHt3OO00WL3bScikjhS4Itlk/ny44AIYNAiuv373YVuuoAAefRTWr4dRo7w7QuqFAlckW5SWwsiR0KQJTJ4M+TU4RdOvH9x6Kzz/PPzxj/VXY45T4Ipki4kTobjYA7Nr15p//3nnwdCh8MtfwqpVez5eakyBK5IN1q6FceN8vO2pp9buOcxg/Hi/QOLKK1NbnwAKXJHsMG6cX8Qwfnz1+m13Zd99/WTbn/8MM2akrDxxClyRTLdgAdxzj58s23//uj/fuHHQsSNcdZUuAU4xBa5IprvmGj9RNm5cap6veXMfv/vKK/DSS6l5TgEUuCKZbd48ePBBv9ChffvUPe+55/qJt1/+Uq3cFFLgimSya67xFunll6f2eZs08bCdOdOHiklKKHBFMtXChTBlClx0Eey9d+qf/8wzoVs3v4BCUkKBK5Kpbr4ZGjWCMWPq5/kbNYKxY32OhVmz6uc1cowCVyQTrVrlQ7dGjfIRBfXlnHOgdWufbUzqTIErkonGj4dt23zMbH1q0cKHmz3+uHdhSJ0ocEUyzZYtcMcdMGyYX6hQ3y6+GBo2hN//vv5fK8spcEUyzQMPwBdf+ATicejYEUaM8C6MjRvjec0spcAVySQheHdC377wL/8S3+tefDFs2uShK7WmwBXJJK+9BnPnegDWZc6EmioqgoEDPew1X26tKXBFMsn48b5czsiR8b/2mDG+xLouhKg1Ba5IpvjsM3jiCTjrLGjaNP7XHz7c+3PvuCP+184SClyRTHHvvb6qw3/9VzKv37ChrwD8zDPwySfJ1JDhFLgimaCkxFd0OOYYX2U3KaNHe9/x3XcnV0MGU+CKZIJnn4Xly+H885Oto3t3+NGPfP7dbduSrSUDKXBFMsGdd0LnznDiiUlX4mufrVoFTz6ZdCUZR4Erku6WLvWRAeec4/2oSTv2WOjRw7s4pEYUuCLp7p57vN/07LOTrsTl5Xn4T5vmw8Sk2hS4Iuls+3a47z44/njvP00XZ53lwauTZzWiwBVJZ888A59/ntxQsF3p3BlOOgn+9CedPKsBBa5IOrvrLujSxVu46Wb0aFizRifPakCBK5Kuli6FF17w/tL8/KSr2dnQoTp5VkMKXJF0VX6y7Kyzkq6kajp5VmMKXJF0lK4ny3akk2c1osAVSUdTp6bnybId6eRZjShwRdJROp8s25FOnlWbAlck3Sxe7CfLzj03PU+W7WjoUCgs1LSN1aDAFUk3d91VcUIqE+TledfHyy/DRx8lXU1aU+CKpJOtW/1k2Y9/7F0KmeKss3yehzvvTLqStKbAFUknjz8Oa9cmPw1jTbVv7ytCTJoEX3+ddDVpS4Erkk4mTIBevWDIkKQrqbnzzoMNG+DBB5OuJG0pcEXSxTvvwOuvw4UXQoMM/NMcNAgOOAD+8Adfzl12koG/VZEs9Yc/+OKQZ56ZdCW1YwYXXeT/cbzxRtLVpCUFrkg6WLcOJk+G//gPaN066Wpq7/TToVUr/89DdqLAFUkH990H33zj3QmZrFkzH7EwZYpfKSffosAVSVpJCYwfD0cdBX37Jl1N3V14oS/nPmFC0pWkHQWuSNKmTIFPP4VLL026ktTo1QuGDfMrzzRE7FsUuCJJCgFuvhn69PHlx7PFZZd5v/T99yddSVpR4Iok6R//gNmzYezYzBwKtitHHAHf/z7ceiuUlSVdTdrIot+wSAa68UZo2xbOOCPpSlLLzLtIFiyAp59Oupq0ocAVScp773kYjRnjZ/ezzSmnQM+ecN11uhAiosAVScp110Hz5nDxxUlXUj/y8+GKK+DNN30ZHlHgiiRi0SJ4+GGfpKZt26SrqT+jRvmqENdem3QlaUGBK5KE667z6QzHjk26kvrVuDH87GcwfbrPE5HjFLgicVuwwKcxPO886NQp6Wrq3+jR0KED/PKXOd+Xq8AVidtvfuMtv6uuSrqSeDRrBuPGwSuv5HxfrgJXJE7vvQcPPQSXXOKtvlwxerQv9z5uXE63chW4InG64gpo2dL7NXNJ48Zw9dU+YuGxx5KuJjEKXJG4PPecb7/6VXaPTNiVUaN8cp7LL/eZ0XKQAlckDtu3+5VXvXtn77jbPcnPh9tug6VL/ZLfHKTAFYnDH//oS4jffDM0apR0NckZPNhnEvvd72DFiqSriZ0CV6S+LVvmQ6KOOw5OOinpapJ3yy0+oc2FF+bcCTQFrkh9CqFiyfM77/RJXXJdz55wzTXw1FO+LHwOUeCK1KcHH4Rnn/WP0D16JF1N+hg7Fg4+2Bed/OKLpKuJjQJXpL4sWeKt24EDPVikQn6+r+O2bh2cc07OdC0ocEXqQ0kJjBzpXz/wAOTlJVtPOjroIJ9T4okn4K67kq4mFgpckfpw1VUwY4YHSWFh0tWkr7Fj4dhj/ba4OOlq6p0CVyTVJk2Cm27y7oQRI5KuJr01aODrnrVv78PFsnxpdQWuSCq99prPGzBkCNx+e9LVZIb27X3liw0bPHQ3bUq6onqjwBVJlTff9JV399kHHnnE57uV6jnwQJg82bsVfvxj2LIl6YrqhQJXJBVmz/a+yIICn4IwF+dKqKthw7w75uWX4eST4euvk64o5RS4InX13HNw1FHQqpWHbZcuSVeUuU4/He6+G154wbtl1q5NuqKUUuCK1FYIPhnLSSfBvvv6qASNSKi7s8+GKVPgnXdgwAC/zRIKXJHaWL3aPwKPHQsnnuirGeTCcjlxOflkeOkl78sdMADGj/f5FzKcAlekJkpL/SPvfvvB88/7SIQnnoAWLZKuLPsMHOit28GDYcwYOPJImDs36arqRIErUh2lpT7yoF8/H/bVr5//8Y8Zowlp6lNBATzzDPz5z7Bwoc+/8O//Dh9+mHRltaLAFdmdJUt84pmePeG003zfI4/4st/f/W6yteUKM18tYv58X6Lo6adh//3hmGP8d7F5c9IVVpuFHJk0QnavqKgozJ49O+kykrdxo4+nnT7duwzeftv3DxkCF1zg/baaFyFZa9bAxIm+LVsGe+0FQ4fC0UfDoEEexvn5sZRiZsUhhKJqH6/AFcihwN261acDXLUKVq70P9glS7z1NG8eLFjgx+Xnw6GHwr/+Kwwf7hczSHopLfUr+6ZM8Skwlyzx/Xvt5V0+3/0u9Onj02J26+arJBcU+PC9FP2nqcCVWvlW4D75JEydWv8vuqt/e+X7K9/uaist9bPXJSX+9fbtvm3d6tuWLb599ZVvVV3BlJ/va41997vQvz98//tw+OE6EZZpliyBN96At97y5eg//LDquRnM/HfbogU0awZNm0KTJr6ycKNGfoVgfr5veXk+38OgQf4JZ6enqlngxtPulsyyaJF/nI7Drk44le+vfFvVVv4HUf7H0bChb40bQ/Pm3trZa6+KP7A2bfwqsPbtvcXTvTt07qxugmywzz6+lU+LCd6/u2yZr5+2apVfSPHFF/Dll/4f8ObNfkXbN9/4f9BffeX/YZeU+FZW5luKJo9XC1eAHOpSEEmhmrZwNUpBRCQmClwRkZioS0EAMLM1wCeVdrUD0nnmkHSuL51rg/SuL51rg53r6xFCKKjuNytwpUpmNrsmfVNxS+f60rk2SO/60rk2qHt96lIQEYmJAldEJCYKXNmViUkXsAfpXF861wbpXV861wZ1rE99uCIiMVELV0QkJgpcEZGYKHBlJ2Z2nJnNN7NFZnZlwrV0M7PpZvaBmc0zs0ui/W3N7EUzWxjdtkmwxjwzm2NmU6P7+5jZrOj9e9jMGiVYW2szm2JmH5nZh2Y2MM3eu7HR7/V9M3vQzJok+f6Z2X1mttrM3q+0r8r3y9zvozrfNbP+e3p+Ba58i5nlAX8Ejgf2B35qZvsnWFIJcFkIYX9gAHBhVM+VwLQQQh9gWnQ/KZcAlZcguAG4NYTQG1gPnJ1IVe524PkQwn7AgXidafHemVkXYAxQFEI4AMgDRpDs+/dn4Lgd9u3q/Toe6BNto4E79vjsIQRt2v65AQOBFyrdvwq4Kum6KtXzFHAMMB/oFO3rBMxPqJ6u0R/hYGAqYPiVSPlVvZ8x19YKWEJ0crzS/nR577oAnwJt8ZkLpwLHJv3+AYXA+3t6v4C7gJ9WddyuNrVwZUflfwTllkf7EmdmhcDBwCygQwihfLLTlUCHhMq6Dfg5UL6k7N7AhhBCSXQ/yfdvH2AN8Keoy+MeM2tGmrx3IYQVwE3AMuBz4EugmPR5/8rt6v2q8d+KAlcygpk1Bx4D/juEsLHyY8GbF7GPbzSzE4HVIYTiuF+7mvKB/sAdIYSDgc3s0H2Q1HsHEPWFDsP/Y+gMNGPnj/Nppa7vlwJXdrQC6FbpftdoX2LMrCEetpNDCI9Hu1eZWafo8U7A6gRKOwL4sZktBR7CuxVuB1qbWfnk/km+f8uB5SGEWdH9KXgAp8N7B3A0sCSEsCaEsB14HH9P0+X9K7er96vGfysKXNnRW0Cf6ExxI/wkxtNJFWNmBtwLfBhCuKXSQ08Do6KvR+F9u7EKIVwVQugaQijE36eXQggjgenAKUnWFtW3EvjUzL4T7RoCfEAavHeRZcAAM2sa/Z7L60uL96+SXb1fTwNnRKMVBgBfVup6qFoSneXa0nsDTgAWAB8D4xKu5Uj8I9y7wDvRdgLeVzoNWAj8HWibcJ1HAVOjr3sCbwKLgEeBxgnWdRAwO3r/ngTapNN7B1wDfAS8D/wFaJzk+wc8iPcnb8c/IZy9q/cLP0H6x+jv5D18tMVun79Ol/aa2XH4R6g84J4QwvW1fjIRkSxX68CNxmsuwIfoLMc/iv40hPBB6soTEckedVm191BgUQhhMYCZPYSfcdxl4LZr1y4UFhbW4SWlrjZs8IVJu3X79v7i4uK1IZq5/pgG/6YZjUR28GLZo7tYYrr66hK4VY1BO2zHg8xsNH4VBt27d0crwybr5z+H8eNhx1+DmX1S9XeISKrU+yiFEMLEEEJRCKGooKDaS/9IPdm+HRo2TLoKkdxUl8BNu/GasmfbtkGjxKZSEcltdQnctBqvKdXz9dfQrFnSVYjkplr34YYQSszsIuAFfFjYfSGEeSmrTOrFxo3QvHnSVYjkprqcNCOE8CzwbIpqkRisWwd77510FSK5SZf25pgVK6Bjx6SrEMlNCtwcUloKS5dCr15JV5IgM99EEqDAzSHvvw8lJXDAAUlXIpKb6tSHK5nllVf8dtCgZOtIVPml7OWtXIvaHKHs24+L1AO1cHPIX/8K+++/82W9IhIPtXBzxMyZ8NZbflmvUNGSDaV+2yAPAGvkfxJhe7TCS1lp3JVJFlMLNwds2wbnnw8dOsAZZyRdjUjuUgs3y4UA48bBO+/Ak09Cy5ZJV5SmopZs2Oq3lu9/GraXX5YXtm332+3bEihOsoVauFksBJ8d7KabvIU7bFjSFYnkNrVws9SXX8KYMXD//XDhhfD73yddUWYJJSXfum3QpInftvLL9MLXWwAo+/rrBKqTTKUWbhZ68kkfjfDXv8LVV/uJsgb6TYskTi3cLDJjBlx3Hfztb9Cvnwfv97+fdFXZoeybb/yL6DYv6gzP7+Fj7MLGTQCUrl8ff3GSMdTuyXClpfDEE3DEEXD44fD663Dttb6ig8JWJL2ohZuhPvwQHn4YJk+GRYugsND7ac88U9MvxqF040b/IrrNi1YzadBvP79dswGAks9Xxl+cpC0FbgZZsMBD9pFHfF4EM79M93e/g5NPhnz9NkXSmv5E09imTfDaa/D3v8OLL8J773nIHnmknwgbPhw6dUq6SgEoXbPGv4huGxR29/0/7A9Ao+Xe4i1duDj+4iRtKHDTyPbtMGuWB+y0aX45bkmJr0F2xBFw221wyinQpUvSlYpIbShwE7R6tQds+TZjBmze7K3YQw6Byy6Do4/2sN1rr6SrlZooWboMgLzolu99B4CvRgwAoOXCrwAIxVqVKpcocGOydatfXjtzpofrzJmwZIk/lpcHffvCqFEwZAgcdRS0bZtouSJSDxS49WDTJnj3XQ/YOXP89t13fRIZ8C6BAQP8ctvDDvPWrFbSzW6l8+YD0CJq0JYdfiAAq8YeDkC7uVsByH+pOP7iJDYK3DpaterbwTpnDixcWDH7X9u2cPDBcMklHq6HHQZduyZbs4gkQ4FbTVu2wAcf+EiBytvKSsMsCws9XEeOhIMO8q+7dtUSWrIze2MuAB3f8PtbT/CrVD65cSAAnWb4ChRNH58Vf3FSbxS4OygthcWLdw7WRYugLFqFpUkTn6vg2GMrgvXAA6F162RrF5H0ltOBu3r1zsE6bx6UTwBl5ivc9u0LI0b4bd++0Lu3n+gSSZXGz74FQK9n/f7Gf/fRDEse6gdAq//vnfxt75sRf3GSMjkRuF9/7UG6Y7iuXl1xTEGBh+m551YE6/e+p5NZIpI6WRe4X3zhJ67efrvidsGCipNYTZr4MuE/+lFFsPbt68vPiKSLlg/MjG79/uoLfDRDs1d9zoaPn+gDQMdb34i/OKm1jA3cEOCzz3YO12XLKo7p1g369/92d0CvXuoOEJFkZEzglpV5t8Crr/r8Aq++Cp9/XvH4vvvCwIG+usHBB/vWrl1y9YqkUvsJ3pLdPMHvl17pLdwT5/n8u3c89CMAuv1WLd50lraBu307FBdXhOs//gHlczt36eJXYw0Y4C3YAw+EFi0SLVdEZI/SKnBD8Am077/fpyAsn3J03319+sEf/MCnIyws1NhWyW1drveW7NTr2wBQ+ls/SXHDEh+3e+rD/w3APldqVEM62WPgmlk34H6gAxCAiSGE282sLfAwUAgsBU4NIdRqfZHFi+Evf/GgXbzYRwYMHw4nneRTEXbsWJtnFRFJL9Vp4ZYAl4UQ3jazFkCxmb0I/CcwLYRwvZldCVwJXFHTAm6+GX72M2+xDhkCv/mNt2Y1HEuk+gp/5S3ZK351GABlN/r+x5b7aIeDohZvr8tmxl+c/NMeAzeE8DnwefT1V2b2IdAFGAYcFR02CXiZGgbuxIketsOHw623+qgCEZFsVaM+XDMrBA4GZgEdojAGWIl3OVTbkiVw3nlQVAQPPOCTbItIavS63Fu8wy/3K9bCrb7/hc/e8ccfOg+A3peqxRunaq/aa2bNgceA/w4hbKz8WAgh4P27VX3faDObbWaz15QvQ4JP6vLDH/r42WefrV3xIiKZpFqBa2YN8bCdHEJ4PNq9ysw6RY93AlZX9b0hhIkhhKIQQlFBtLIpQMOG8OST3sIdPhxOOMEXSPzmmzr9PCJShd5jZ9J77EyO7XwQx3Y+CAtgwft4H1s+k49vHMjH0UxlUn/2GLhmZsC9wIchhFsqPfQ0MCr6ehTwVE1fvEULeP55uOIKn9tgxAgfkTB6tI+7LZ+dS0QkG1gIVfYEVBxgdiTwGvAeUB6Bv8D7cR8BugOf4MPCvtjdcxUVFYXZs2dX+VhpKbz8MkyaBI895hPOtGnjw8IGDfIxuP37e8tYUs/MikMIRQDHNPi33f+jkKyz5Hpv3T5y2m0AnPZXH9VQPvpB4MWyR+s8+r86oxReB3b1QkPqWkC5vDwfFjZkCEyYAE89BdOn+5Vmf/ubH9O0qV++W34BxKGHaviYiGSOPbZwU2l3LdzdWbnSr0Arn0dh7ly/Ks0MvvMdnzehf/+KORS0AGPNqYUrlS272mcnu+C0ZwC46y8+V0P5FW65KJYWbjro2BFOOcU3gA0b4I034K23fJTD66/Dgw9WHN+jR0UA9+/vqzJ07qzLgUUkWRkRuDtq3dpHNZxwQsW+tWs9fMu3t9/2URDlDfg2bSqmaOzXz28POECT3ohUpfv/RHM1/I/P1ZA31veXz8e79MHeABTcoT7emsiILoXa+uor736YO9dHQbz7Lrz/vu8vV1j47YnI+/b1yXJy7eScuhSkJr44y0+yfTl0MwB7P90UqJg4PRvlTJdCbbVo4aMcjjyyYl8I8MknOy+389xzUFLixzRqBPvtt3MQawVeEamLrG7h1sTWrTB/vreCKwfx8uUVx7Ru7d0QO3ZNtGyZXN2pohau1MXm4T5pzsoBPrS/63RvvZQvjpkN1MJNocaNPUD79fv2/vXrvRuicghPnlwxVy/4Kr7ly6WX33bsqNawiHybWri1EAJ8+qm3ht95x7c5c3wu33Lt21cE8EEHwSGHeDCnawirhSupVDL4EADWHtiYDrO8n9femJtkSXWmFm5CzKB7d99OPLFi/5dfegjPmVMRwrfc4ssFgY8PPuww3wYM8As32rRJ5mcQkfgpcFOoVSu/Cu4HP6jYt20bfPABzJ4Ns2bBzJk+f0T5B4t99/XwPewwOPxw79JoUO053ETSU/5LxQB0fAnskO8BsHGETxXZ+j1fGKZ03vxkikuQAreeNWpU0a1wzjm+b+NGD+CZMz2En3/elxcCX2n4hz+Eo4/2y5x79kzfbggRqRkFbgJatoTBg32DiqFqr74K06b59uij/lhhYcUcE8cdpy4IyTyheB4ALYqjHX16AlD6w/4ANF7s82SXfPJp7LXFTYGbBsw8WAsL4YwzPIDnz68I38ceg3vv9Ysxhg6FU0+FYcO8C0NEMocCNw2Z+YUX++0HF17oU1e+9ZYH7yOPwDPPeFfFccd5+A4fDk2aJF21SPWULvThPHkLox1dOgPQoN9+ANjn6/y4SivEZAudnskAeXl+Yu3GG2HpUpgxw4O4uBhOP90n6/ntb2HduqQrFZHdUeBmGDMP31tugWXL4O9/9zG+V1/tqx5fdJH3B4tkipIVn1Gy4jPK3v2Isnc/8mvsS0rI79GN/B7daNCiBQ2yZJYpBW4Ga9DAT6Y9+2zFEkUTJ8L++8Ntt3lXhIikDwVuljjgALjvPli4EI46CsaO9XG98+YlXZlIzZSuX0/p+vWUfPKpj1woK4OyMvLa7U1eu71p0KQJDTL0pIUCN8v06AFTp8IDD/ilxocf7gtyikjyFLhZyAx++lOfhL1jRx9KNm1a0lWJ1E7Z5s2Ubd5M6dp1lK5dRwiBEAINmjWjQbNmWH4+lp8ZA64UuFmsWze/mKJnT1+eaMWKpCsSyW0K3CzXoQM88YTP93vuuRVzOIhkqrB1K2Hr1n+2fMtZ48ZY48b+ES9Nr4dX4OaA3r3h2mt9VYvp05OuRiR3KXBzxHnn+fSQEyYkXYlIaoWSEt+ilu8/NcjzLY0ocHNEkyYwahQ89RRU+hQFYQepAAAGdElEQVQmIjFS4OaQoUP9Ip6Z2buwqoifqAgBykp9K+/TTYO+XQVuDhnoK1vz5pvJ1iGSqzJj8JqkRKtWPmrh44+TrkQkRmk0NEct3BxTWOiT3ohI/BS4OaZdO03jKJIUBW6OadMG1q9PugqR3FTtwDWzPDObY2ZTo/v7mNksM1tkZg+bWaP6K1NSpWlT2LIl6SpEclNNWriXAB9Wun8DcGsIoTewHjg7lYVJ/WjSRIErkpRqBa6ZdQV+BNwT3TdgMDAlOmQS8JP6KFBSq2FD2L496SpEclN1W7i3AT8HyqL7ewMbQggl0f3lQJeqvtHMRpvZbDObvSYLF4XLNPn5fvGDiMRvj4FrZicCq0MIxXs6tiohhIkhhKIQQlFBQUFtnkJSqEEDn0BfROJXnQsfjgB+bGYnAE2AlsDtQGszy49auV0BzbaaARS4IsnZYws3hHBVCKFrCKEQGAG8FEIYCUwHTokOGwU8VW9VSsqYpdWFNyI5pS7jcK8ALjWzRXif7r2pKUnqkwJXJDk1mkshhPAy8HL09WLg0NSXJPUpTSfCF8kJutJMRCQmCtwcoxauSHIUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMalW4JpZazObYmYfmdmHZjbQzNqa2YtmtjC6bVPfxYqIZLLqtnBvB54PIewHHAh8CFwJTAsh9AGmRfdFRGQX9hi4ZtYKGATcCxBC2BZC2AAMAyZFh00CflJfRYqIZIPqtHD3AdYAfzKzOWZ2j5k1AzqEED6PjlkJdKjqm81stJnNNrPZa9asSU3VIiIZqDqBmw/0B+4IIRwMbGaH7oMQQgBCVd8cQpgYQigKIRQVFBTUtV4RkYxVncBdDiwPIcyK7k/BA3iVmXUCiG5X10+JIiLZYY+BG0JYCXxqZt+Jdg0BPgCeBkZF+0YBT9VLhSIiWSK/msddDEw2s0bAYuBMPKwfMbOzgU+AU+unRBGR7FCtwA0hvAMUVfHQkNSWIyKSvXSlmYhITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxqVbgmtlYM5tnZu+b2YNm1sTM9jGzWWa2yMweNrNG9V2siEgm22PgmlkXYAxQFEI4AMgDRgA3ALeGEHoD64Gz67NQEZFMV90uhXxgLzPLB5oCnwODgSnR45OAn6S+PBGR7LHHwA0hrABuApbhQfslUAxsCCGURIctB7pU9f1mNtrMZpvZ7DVr1qSmahGRDFSdLoU2wDBgH6Az0Aw4rrovEEKYGEIoCiEUFRQU1LpQEZFMV50uhaOBJSGENSGE7cDjwBFA66iLAaArsKKeahQRyQrVCdxlwAAza2pmBgwBPgCmA6dEx4wCnqqfEkVEskN1+nBn4SfH3gbei75nInAFcKmZLQL2Bu6txzpFRDJe/p4PgRDCr4Ff77B7MXBoyisSEclSutJMRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQI3xwwYAJdcknQVIrnJQgjxvZjZGuCT2F5QaqJHCKEg6SJEslmsgSsiksvUpSAiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITP4PrQ161dhEqnEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9e5978a588>"
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XecVfWd//HXZxq9ykjXERsxgkpQscQ1Yk00ZKNrdM1KbKyJUZeYRA3ZbPztRpONNSQWLAkmhGiwBks2MdgFBbEhggiIImXonWnf3x+fc5nCDNPPueX9fDzO43LvPXPvhwvzns98z/d8j4UQEBGR9peXdAEiIrlCgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrkjCzOyLZragHV73AjP7vybu+y0ze7m5z0nzKHAl40WB8K6ZbTOzlWZ2l5n1jJ6728y2RFuZmZXXuP9MDLUFMztgT/uEEF4KIRzcwtc/3sxeNbONZrbOzF4xsyOj150SQji1Ja8r7UOBKxnNzK4BfgH8AOgBjAL2Bf5mZkUhhMtDCF1DCF2BG4GHUvdDCGckV7kzs4JWfG13YDowEegNDARuAHa2TXVtrzV/32ygwJWMFQXODcCVIYRnQwjlIYSlwLlACfDNFrzmiWb2qZn90MxWm9kKM/uamX3ZzBZGXeSPaux/lJm9ZmYbon1/bWZF0XMvRru9HXXU36jx+tea2Urgt6nHoq/ZP3qPEdH9AWZWamYn1lPuQQAhhKkhhMoQwvYQwv+FEN6JvrbWUEDUbV9uZh9G9f7GzKyBz+GXZvaymfWo8djNZrbezJaY2Rk1Hh9gZk9GdS8ys8tqPPdTM5tmZn8ws03At6LHHjazB81ss5nNM7ORzfynykgKXMlkxwIdgUdrPhhC2AI8DZzSwtftF73uQOAnwL14eH8B+CLwn2a2X7RvJTAe6AMcA4wGvhPVcUK0z2FRR/1QjdfvjXfi4+rU/hFwLfAHM+sM/BaYHEJ4vp46FwKVZjbZzM4ws15N+LudCRwJDMd/MJ1W80kzyzOze6PnTw0hbIyeOhpYEP09/xe4v0ZY/wn4FBgAnAPcaGYn1XjZMcA0oCcwJXrsq9HX9QSeBH7dhNozngJXMlkfYE0IoaKe51ZEz7dEOfCzEEI5Hgp9gDtCCJtDCPOA94HDAEIIc0IIM0MIFVF3fQ/wT428fhXwXyGEnSGE7XWfDCHcCywCZgH9gQn1vUgIYRNwPBDwHwqlUafZdw/v/fMQwoYQwjJgBnB4jecKgan4D4OzQgjbajz3cQjh3hBCJTA5qquvmQ0GjgOuDSHsCCG8BdwHXFjja18LITweQqiq8fd9OYTwdPR6vyf6PLOdAlcy2RqgTwPjgv2j51tibRQEAKmAWFXj+e1AVwAzO8jMpkcH6zbh48SNBX1pCGFHI/vcCxwKTAwhNDgmG0KYH0L4VghhULT/AOD2Pbzuyhp/3pb6e0QOwLvRG0IIZQ19XY0g7hq937oQwuYa+36M/3aQ8kkT6uiYC+O7ClzJZK/hB4i+XvNBM+sKnAE8F0MNdwEfAAeGELoDPwLqHRetYY9L9EX13w7cD/zUzHo3pZAQwgfA7/DgbYn5wEXAM2bW1FkTnwG9zaxbjcf2AZbXLK2F9WQdBa5krGh88QZgopmdbmaFZlYCPIyPKf4+hjK6AZuALWY2FPh2nedXAUOa+Zp3ALNDCJcCTwF317eTmQ01s2vMbFB0fzBwPjCzme+3SwhhKv5D4+9mtn8T9v8EeBW4ycw6mtlw4BLgDy2tIZspcCWjhRD+Fw+Im/Hgm4X/Cjt6T7+Kt6HvA/8KbMaHAR6q8/xPgcnRrIBzG3sxMxsDnE51cH8PGGFmF9Sz+2b8YNYsM9uKB+17wDUt+HvsEkKYDPw/4B/RD7DGnI/PCvkMeAwfn/57a2rIVqYFyEVE4qEOV0QkJgpcEZGYKHBFRGKiwBURiUnWTzSWpunTp08oKSlJugyRjDJnzpw1IYTipu6vwBUASkpKmD17dtJliGQUM/u4OftrSEFEJCYKXJFcsX07vPEGLFoEmn+fCAWuSLabOxdOOgm6d4ejjoIDD4TeveF734PNmxv/emkzClyRbFVVBdddByNHwrx58IMfwCOPwL33whlnwO23w9ChMGNG0pXmDB00E8lGVVXw7W/DpElwySXwy19Crxrrk196KVx9NVx0EXz5yzB9OowenVy9OUIdrkg2uvJKD9sf/cg72l71XAzi6KPhhRd8iOHMM+HFF3ffR9qUAlck2zz4INx5J1xzDfzP/0D9ly1zxcXw3HOwzz7wjW/A6tXx1ZmDFLgi2WTBAvjOd+CEE+DnP99z2KYUF8Of/wzr18PYsT4cIe1CgSuSLSor4YILoGNHmDIFCppxiGb4cLjtNnj2WfjNb9qvxhynwBXJFpMmwZw5HpiDBjX/6y+/HE49FX78Y1i1qvH9pdkUuCLZYM0amDDB59ue2+iFJepnBhMn+gkS113XtvUJoMAVyQ4TJvhJDBMnNm3ctiEHHeQH2373O3jttTYrT5wCVyTTLVwI993nB8sOOaT1rzdhAvTrB9dfr1OA25gCVyTT3XCDHyibMKFtXq9rV5+/+8IL8I9/tM1rCqDAFcls8+bB1Kl+osPee7fd6152mR94+/GP1eW2IQWuSCa74QbvSH/wg7Z93Y4dPWxnzvSpYtImFLgimerDD2HaNPjud2Gvvdr+9S+6CAYP9hMopE0ocEUy1S23QFERXHVV+7x+URGMH+9rLMya1T7vkWMUuCKZaNUqn7o1dqzPKGgvl14KPXv6amPSagpckUw0cSKUlfmc2fbUrZtPN3v0UR/CkFZR4Ipkmu3b4a67YMwYP1GhvV15JRQWwq9+1f7vleUUuCKZ5o9/hHXrfAHxOPTrB+ed50MYmzbF855ZSoErkklC8OGEYcPgn/4pvve98krYssVDV1pMgSuSSV56Cd5+2wOwNWsmNNfIkXDMMR72Wi+3xRS4Iplk4kS/XM4FF8T/3ldd5ZdY14kQLabAFckUn30Gjz0GF18MnTvH//5nn+3juXfdFf97ZwkFrkimuP9+v6rDv/97Mu9fWOhXAH7qKfj442RqyHAKXJFMUFHhV3Q45RS/ym5Sxo3zseN7702uhgymwBXJBE8/DZ9+Ct/+drJ17LMPfOUrvv5uWVmytWQgBa5IJrj7bhgwAM48M+lK/Npnq1bB448nXUnGUeCKpLulS31mwKWX+jhq0k47Dfbd14c4pFkUuCLp7r77fNz0kkuSrsTl53v4P/ecTxOTJlPgiqSz8nJ44AE44wwfP00XF1/swauDZ82iwBVJZ089BStWJDcVrCEDBsBZZ8Fvf6uDZ82gwBVJZ/fcAwMHeoebbsaNg9JSHTxrBgWuSLpauhT++lcfLy0oSLqa3Z16qg6eNZMCVyRdpQ6WXXxx0pXUTwfPmk2BK5KO0vVgWV06eNYsClyRdDR9enoeLKtLB8+aRYErko7S+WBZXTp41mQKXJF0s3ixHyy77LL0PFhW16mnQkmJlm1sAgWuSLq5557qA1KZID/fhz6efx4++CDpatKaAlcknezc6QfLvvpVH1LIFBdf7Os83H130pWkNQWuSDp59FFYsyb5ZRiba++9/YoQkyfDtm1JV5O2FLgi6eTOO2H//WH06KQrab7LL4cNG2Dq1KQrSVsKXJF08dZb8PLLcMUVkJeB35onnACHHgq//rVfzl12k4H/qiJZ6te/9otDXnRR0pW0jBl897v+g+PVV5OuJi0pcEXSwdq1MGUK/Nu/Qc+eSVfTct/8JvTo4T88ZDcKXJF08MADsGOHDydksi5dfMbCtGl+ppzUosAVSVpFBUycCCeeCMOGJV1N611xhV/O/c47k64k7ShwRZI2bRp88gl873tJV9I29t8fxozxM880RawWBa5IkkKAW26BAw/0y49ni2uu8XHpBx9MupK0osAVSdIrr8Ds2TB+fGZOBWvIccfBkUfCbbdBVVXS1aSNLPoXFslAv/wl9O4NF16YdCVty8yHSBYuhCefTLqatKHAFUnKu+96GF11lR/dzzbnnANDhsBNN+lEiIgCVyQpN90EXbvClVcmXUn7KCiAa6+F11/3y/CIAlckEYsWwUMP+SI1vXsnXU37GTvWrwpx441JV5IWFLgiSbjpJl/OcPz4pCtpXx06wPe/DzNm+DoROU6BKxK3hQt9GcPLL4f+/ZOupv2NGwd9+8KPf5zzY7kKXJG4/fSn3vldf33SlcSjSxeYMAFeeCHnx3IVuCJxevdd+NOf4OqrvevLFePG+eXeJ0zI6S5XgSsSp2uvhe7dfVwzl3ToAD/5ic9YeOSRpKtJjAJXJC7PPOPbf/5nds9MaMjYsb44zw9+4Cuj5SAFrkgcysv9zKsDDsjeebeNKSiA22+HpUv9lN8cpMAVicNvfuOXEL/lFigqSrqa5Jx0kq8k9rOfwfLlSVcTOwWuSHtbtsynRJ1+Opx1VtLVJO/WW31BmyuuyLkDaApckfYUQvUlz+++2xd1yXVDhsANN8ATT/hl4XOIAlekPU2dCk8/7b9C77tv0tWkj/Hj4Ygj/KKT69YlXU1sFLgi7WXJEu9ujznGg0WqFRT4ddzWroVLL82ZoQUFrkh7qKiACy7wP//xj5Cfn2w96ejww31Nicceg3vuSbqaWChwRdrD9dfDa695kJSUJF1N+ho/Hk47zW/nzEm6mnanwBVpa5Mnw803+3DCeeclXU16y8vz657tvbdPF8vyS6srcEXa0ksv+boBo0fDHXckXU1m2Htvv/LFhg0eulu2JF1Ru1HgirSV11/3K+/utx88/LCvdytNc9hhMGWKDyt89auwfXvSFbULBa5IW5g928cii4t9CcJcXCuhtcaM8eGY55+Hr38dtm1LuqI2p8AVaa1nnoETT4QePTxsBw5MuqLM9c1vwr33wl//6sMya9YkXVGbUuCKtFQIvhjLWWfBQQf5rATNSGi9Sy6BadPgrbdg1Ci/zRIKXJGWWL3afwUePx7OPNOvZpALl8uJy9e/Dv/4h4/ljhoFEyf6+gsZToEr0hyVlf4r79Ch8OyzPhPhscegW7ekK8s+xxzj3e1JJ8FVV8Hxx8PbbyddVasocEWaorLSZx4MH+7TvoYP92/+q67SgjTtqbgYnnoKfvc7+PBDX3/hX/8V5s9PurIWUeCK7MmSJb7wzJAh8I1v+GMPP+yX/f7c55KtLVeY+dUiFizwSxQ9+SQccgiccor/W2zdmnSFTWYhRxaNkD0bOXJkmD17dtJlJG/TJp9PO2OGDxm8+aY/Pno0fOc7Pm6rdRGSVVoKkyb5tmwZdOoEp54KJ58MJ5zgYVxQEEspZjYnhDCyyfsrcAVyKHB37vTlAFetgpUr/Rt2yRLvnubNg4ULfb+CAjjqKPjnf4azz/aTGSS9VFb6mX3TpvkSmEuW+OOdOvmQz+c+Bwce6MtiDh7sV0kuLvbpe230Q1OBKy1SK3AffxymT2//N23o/17q8Zq3DW2VlX70uqLC/1xe7tvOnb5t3+7b5s2+1XcGU0GBX2vsc5+DESPgyCPh2GN1ICzTLFkCr74Kb7zhl6OfP7/+tRnM/N+2Wzfo0gU6d4aOHf3KwkVFfoZgQYFv+fm+3sMJJ/hvOLu9VPMCN56+WzLLokX+63QcGjrglHq85m19W+obIvXNUVjoW4cO0LWrdzudOlV/g/Xq5WeB7b23dzz77AMDBmiYIBvst59vqWUxwcd3ly3z66etWuUnUqxbBxs3+g/grVv9jLYdO/wH9ObN/gO7osK3qirf2mjxeHW4AuTQkIJIG2puh6tZCiIiMVHgiojEREMKAoCZlQIf13ioD5DOK4ekc33pXBukd33pXBvsXt++IYTipn6xAlfqZWazmzM2Fbd0ri+da4P0ri+da4PW16chBRGRmChwRURiosCVhkxKuoBGpHN96VwbpHd96VwbtLI+jeGKiMREHa6ISEwUuCIiMVHgym7M7HQzW2Bmi8zsuoRrGWxmM8zsfTObZ2ZXR4/3NrO/mdmH0W2vBGvMN7O5ZjY9ur+fmc2KPr+HzKwowdp6mtk0M/vAzOab2TFp9tmNj/5d3zOzqWbWMcnPz8weMLPVZvZejcfq/bzM/Sqq8x0zG9HY6ytwpRYzywd+A5wBHAKcb2aHJFhSBXBNCOEQYBRwRVTPdcBzIYQDgeei+0m5Gqh5CYJfALeFEA4A1gOXJFKVuwN4NoQwFDgMrzMtPjszGwhcBYwMIRwK5APnkezn9zvg9DqPNfR5nQEcGG3jgLsaffUQgjZtuzbgGOCvNe5fD1yfdF016nkCOAVYAPSPHusPLEionkHRN+FJwHTA8DORCur7PGOurQewhOjgeI3H0+WzGwh8AvTGVy6cDpyW9OcHlADvNfZ5AfcA59e3X0ObOlypK/VNkPJp9FjizKwEOAKYBfQNIaQWO10J9E2orNuBHwKpS8ruBWwIIVRE95P8/PYDSoHfRkMe95lZF9LkswshLAduBpYBK4CNwBzS5/NLaejzavb3igJXMoKZdQUeAf4jhLCp5nPB24vY5zea2ZnA6hDCnLjfu4kKgBHAXSGEI4Ct1Bk+SOqzA4jGQsfgPxgGAF3Y/df5tNLaz0uBK3UtBwbXuD8oeiwxZlaIh+2UEMKj0cOrzKx/9Hx/YHUCpR0HfNXMlgJ/wocV7gB6mllqcf8kP79PgU9DCLOi+9PwAE6Hzw7gZGBJCKE0hFAOPIp/puny+aU09Hk1+3tFgSt1vQEcGB0pLsIPYjyZVDFmZsD9wPwQwq01nnoSGBv9eSw+thurEML1IYRBIYQS/HP6RwjhAmAGcE6StUX1rQQ+MbODo4dGA++TBp9dZBkwysw6R//OqfrS4vOroaHP60ngwmi2wihgY42hh/olMViuLb034MvAQuAjYELCtRyP/wr3DvBWtH0ZHyt9DvgQ+DvQO+E6TwSmR38eArwOLAL+DHRIsK7DgdnR5/c40CudPjvgBuAD4D3g90CHJD8/YCo+nlyO/4ZwSUOfF36A9DfR98m7+GyLPb5+q07tNbPT8V+h8oH7Qgg/b/GLiYhkuRYHbjRfcyE+RedT/FfR80MI77ddeSIi2aM1V+09ClgUQlgMYGZ/wo84Nhi4ffr0CSUlJa14S2mtDRv8wqSDB9d+fM6cOWtCtHL9KXn/ohWNROr4W9WfG7jEdNO1JnDrm4N2dN2dzGwcfhYG++yzD7oybLJ++EOYOBHq/jOY2cf1f4WItJV2n6UQQpgUQhgZQhhZXNzkS/9IOykvh8LCpKsQyU2tCdy0m68pjSsrg6LEllIRyW2tCdy0mq8pTbNtG3TpknQVIrmpxWO4IYQKM/su8Fd8WtgDIYR5bVaZtItNm6Br16SrEMlNrTloRgjhaeDpNqpFYrB2Ley1V9JViOQmndqbY5Yvh379kq5CJDcpcHNIZSUsXQr77590JQkwa94m0g4UuDnkvfegogIOPTTpSkRyU6vGcCWzvPCC355wQrJ1pAVrpNdINbmhqvbjrVh7REQdbg75wx/gkEN2P61XROKhDjdHzJwJb7zhp/XmpFRnmhqfTXWuUadreVbrPqn7VXU62ujrQurx1Ouo85UmUIebA8rK4Nvfhr594cILk65GJHepw81yIcCECfDWW/D449C9e9IVJaxup9sASz1fVP+3iEWvEyor/YGqOvfV+Uo91OFmsRB8dbCbb/YOd8yYpCsSyW3qcLPUxo1w1VXw4INwxRXwq18lXVGaSXWeoTK66x2t5VP7frS7FUbfKvnRDlEHnOp0iTrbUBl1tuXlte6r8xVQh5uVHn/cZyP84Q/wk5/4gbI8/UuLJE4dbhZ57TW46Sb4y19g+HAP3iOPTLqqDFFnTHbXNNzU89FPLIs6XEstKlwQfQtFsxosNXuhosK/PtXp7izzx1P3U8/v6nzV8eYC9T0ZrrISHnsMjjsOjj0WXn4ZbrzRr+igsBVJL+pwM9T8+fDQQzBlCixaBCUlPk570UVafrFVGph9sJu8qNONOtzQIVrVPRrrDflRx1vpX29l3tnaDu90w/bt0e0Ovy1LjflqrDebKXAzyMKFHrIPP+zrIpj5abo/+xl8/evVv92KSHrSt2ga27IFXnoJ/v53+Nvf4N13PWSPP94PhJ19NvTvn3SVWaruLIayqONMdbxVqQ7Ub1NnqoWow63q7B1vVQfvhENqjLfCvz5/u3e0eZu9w83bGnW8W7f5baoD1lhvVlHgppHycpg1ywP2uef8dNyKCr8G2XHHwe23wznnwMCBSVcqIi2hwE3Q6tUesKnttddg61bvYr/wBbjmGjj5ZA/bTp2SrjbHpcZ2K2qPtealOs9ovm1qlkJeNE+3qtD/4cq7+bdaRaeoE87rAEB+uV9grnCLv07R+p3+9Ru2+utt2uL7b/H7u8Z6ozrU8WYWBW5Mdu7002tnzvRwnTkTlizx5/LzYdgwGDsWRo+GE0+E3r0TLVdE2oECtx1s2QLvvOMBO3eu377zji8iAz4kMGqUn2579NHezepKuhmizthu1c6os011uNE/cl4037awvCLa36eOVBV5Z7u9u3e6FV18ZmaIZj3k7/D5vR03+H+ITqW9ACha42O7+Ws3+f5bvPOtSs1ySL1PVWUb/CWlvShwW2nVqtrBOncufPhh9fdl795wxBFw9dUerkcfDYMGJVuziCRDgdtE27fD++/7TIGa28qV1fuUlHi4XnABHH64/3nQIF0iK6ulxnbLo/m1qTUVotkFedEZZkXRbf7OHgBYpY/tbin0znZHsf8n2RbNOtkc/acp2Oodccc1fttlpS/31vkz72wLV23099sYdb7RLIeqstQYr+bzphMFbh2VlbB48e7BumhR9Uygjh19rYLTTqsO1sMOg549k61dRNJbTgfu6tW7B+u8ebDNmwTM/Aq3w4bBeef57bBhcMAB1YtGidRSVXveblVqbYad0eyDaMy161b/6VywtRsAm8t83u6Wwd7Z7ugfjckO8s5488F+d8Mm36/jZz7G221Z5+h2LwA6rPBON3/NBq9ja2p2Q+0OXB1vMnIicLdt8yCtG66rV1fvU1zsYXrZZdXB+vnP62CWiLSdrAvcdev8wNWbb1bfLlxY/QO9Y0e/TPhXvlIdrMOG+eVnRNrMrnm7tc8Us2g2QarT7bTJO9CiDT4bocNG/wm/cZt/a27Z31+u90DvWAcPXAFAxVCf3bB0nc8fXPuxd8rdFhcD0HOxd9CdPtkMQP7q9QBUbY7m9UYdtzreeGVs4IYAn322e7guW1a9z+DBMGJE7eGA/ffXcICIJCNjAreqyocFXnzR1xd48UVYsaL6+YMOgmOO8asbHHGEb336JFevSC11ZjNU1h3bjWYX9Fjnsxg6rfaOd2NpRwA2HOxjtIsO9nm6xw1aDMDJxfMB2HaAz2J49YghAMxb7Od/d1ngHXCvhf66XZZ6h7ur490QzXJQxxuLtA3c8nKYM6c6XF95Bdb7/xEGDvSzsUaN8g72sMOgW7dEyxURaVRaBW4IvoD2gw/6EoSb/IArBx3kyw9+8Yu+HGFJiea2SoZLzWbY6beVu8Z2fZWwwk0+9tpnjXeoXVZ4h7r+M7/9v88PA2D5UB+rPbffGwDcsO+7AKwd5GPB04cfDsCzHx4CQME8//peC/zMt26LvZPOX6mONw6NBq6ZDQYeBPriVxyZFEK4w8x6Aw8BJcBS4NwQwvqWFLF4Mfz+9x60ixf7zICzz4azzvKlCPv1a8mrioikFwuN/OQys/5A/xDCm2bWDZgDfA34FrAuhPBzM7sO6BVCuHZPrzVy5Mgwe/bsWo/dcgt8//vesY4eDRde6N2spmPFy8zmhBBGApyS9y9qZ5KS+tXNfBZCXkcfm7Vu3pGGvt7xbt3PzzhbN9R7pu3DvDP+ytD3ABi71ysAHFzo84E/izrVqRv8uksPfTjCX+8df53eH/h+3T6KZjWo493N36r+3OrfqxvtcEMIK4AV0Z83m9l8YCAwBjgx2m0y8Dywx8Cta9IkD9uzz4bbbvNZBSIi2arRDrfWzmYlwIvAocCyEELP6HED1qfuN6Rmh7tkiU/R+sIX/IBYUVHL/gLSNtThpqmo4911teAO3vHmdfejxJX9fPbClv29A143NLrCxGHeqZ530BwALuz5OgD7Ffp+S8p9tsKDG44CYOqCL/jrvuOv2/sD72S7fhTNali5FoCqaGw5Fzvetuhwm3zVXjPrCjwC/EcIYVPN54Kndr2fuJmNM7PZZja7tLR01+ODBsGXvuTzZ59+umXFi4hkkiZ1uGZWCEwH/hpCuDV6bAFwYghhRTTO+3wI4eA9vU7dMdzNm+GUU+CNN3whmLFjYcwYPxtM4qUON0M01PH29NkHlf181sHm/b1TXRedkWbDvUf6xoFvAnB+T5/VMCB6nQXlvt/ktccB8NQHhwLQ6V1f1az3Bz6LossSfx1btQ6AEJ25VrXDO95sXp0slg43Gi64H5ifCtvIk8DY6M9jgSea++bdusGzz8K11/raBued5zMSxo3zYYbU6lwiItmgKbMUjgdeAt4FUhH4I2AW8DCwD/AxPi1s3Z5eq75ZCimVlfD88zB5MjzyiC8406uXTws74QSfgztiBBQWNuevJ02lDjdD1e14o4vfVY/xRh3vEB+7XX+w71fxeV/D4fQD3wfgzJ5vAbBXvj/+7k5fJX/aSh/bnTffj2j3eN+Ps/da6GfMdVrmsxgojTre6Npru9bjzaIrUMQ1S+FloKE3Gt3aAlLy831a2OjRcOed8MQTMGOGn2n2l7/4Pp07++m7qRMgjjpK08dEJHM0a5ZCa+2pw92TlSv9DLTUOgpvv+1DRGZw8MG+bsKIEdVrKOgCjM2nDjdLNDbGu7d3vFv3izreA32/rQd7x/r5IcsBOLa3r9XQOd/HZudv9UtRvPKpr9VQtsDn7/Zc4G/b4yNf/axouc/fDeu9863ati1rrrcWS4ebDvr1g3PO8Q1gwwZ49VU/2DZ3rofx1KnV+++7b3UAjxjhV2UYMECnA4tIsjKiw22KNWs8fFPbm2/Wvphjr17VSzQOH+63hx6qRW9S1OFmqYY63ujMtarjtFm6AAAMPElEQVSo490+2L8RNgzxgySbh/jhmk77+rzbkt4+Rltg/vgnm7xjXrfcp953/ch7tx6Lo/m7H/tYbv7K9YRo7u6uKwxXpK63lln/zXKmw22KPn18itkpp1Q/tnmzDz+8/bbPgnjnHV+vYfPm6n1KSmovRD5smC+Wo4NzkhUaWAg9deKCRWHYeYUfDOkcLVy+s78PGWzex28/2scDdscAf52C7mW1blMLpVd09kjZ2cMDvHuPDnSIXjtvbXTZn9Qi6Dl42Z+sCdz6dOvmsxyOP776sRDg4493v9zOM89A9H+SoiIYOnT3INYVeEWkNbI6cOtj5l1tSYmvRpaycycsWOBdcCqEX3gBpkyp3qdnTx+GqDs00b173H8LkRZqoOMlOnEhL/r1r8Nqvzhlx6XRwbG+3uFuG+BnJW3t58/vSC3y38Vft6Kz324d6J1JZVERXbt619y5iw9n5EeLqoeN0VBDdNXWbDm4tic5F7gN6dDBA3T48NqPr18P771XuxueMqV6rV7wq/imLpeeuu3XT92wiNSmwG1Er14+7/eLX6x+LAT45BPvht96y7c334Rp06r32Xvv6gA+/HBfpOeAAxTCkmZS46bBu8qqHVF3GZ24YNGJDAXrfLpXj0+ihcv38s53Z3Gq0/WDHjt6Rp1tdHp+VZGxrU+0oE5+1DV38PuFHXzFqvyNfrvrApfRe2fqwbU9UeC2gBnss49vZ55Z/fjGjR7Cc+d6CM+dC7fe6pcLAp8ffPTRvo0a5Sdu9OqVzN9BROKnwG1DPXrs3g2XlcH778Ps2TBrFsyc6etHpH5oH3SQh+/RR8Oxx/qQRl6T13ATaSd1LwEUzSiwrb7Qua33GQcdV/oMhI7do2lmPf1+WS8fry3vmk9loXe9VdHVsst6Rmuxms9kKCj0GMqLpgaFaBw57PBvhF2nCWfBwjgK3HZWVFQ9rHDppf7Ypk0ewDNnegg/+6xPVwOf3valL8HJJ/tpzkOGaBhCJFsocBPQvTucdJJvUD1V7cUX4bnnfPvzn/25kpLqNSZOP11DEJKQOpd5T42vpi7zbhv9KHLeGl88p1OXaLy2SyequvmAbmWnqIMtiDqI6De5qq7e8Vrw7tjyoucLfPZCXuqEiSwY21XgpoGaU9UuvND/Hy1YUB2+jzwC99/vJ2Oceiqce66vG9yjR9KVi0hzKHDTkJmfeDF0KFxxhS9d+cYbHrwPPwxPPeVDFaef7uF79tlatF1i1tB83jqzG6yokLxO/p8zP1o6MnSMxnCLoo43v/aYWYhmL1hVnQ42dXCjzPffNW83g8Z2dXgmA+Tn+4G1X/4Sli6F117zIJ4zB775TV+s57//G9auTbpSEdkTBW6GMfPwvfVWWLYM/v53n+P7k5/4VY+/+10fDxaJVQi+VVVCVSWhvIxQXkbVtm1UbdjoW+kaqkrXEFaWElaWwqo1sGoNeaUbyCvdgG3a6tvOMmxn2a6XtsJC34pSW5FvhQW+5ef74jxmaX+EWYGbwfLy/GDa009XX6Jo0iQ45BC4/XYfihCR9KHAzRKHHgoPPOBLUp54Iowf7/N6581LujLJaSEQKioIFRVU7dzp29Ztvm3c7Nv6Db5FnXDYtMW3bdt8Ky8npM4eAiw/z7eCAqygwMfcoi433TtdBW6W2XdfmD4d/vhHWLzYQ/eVV5KuSkRAsxSykhmcf74vS3nyyT6V7MknffhBJDF11m0IqbPZKqJuNDXDIVosnbzai6c31LVa9Hiw6JLw+am3qTOmlgazGNThZrHBg/1kiiFD/PJEy5cnXZFIblPgZrm+feGxx3y938suS4sf8iK1NTDDIZT5VrVjJ1U7dhJ2RltZee0tBOq9VJjl+bbrfvJjuwrcHHDAAXDjjX5Vixkzkq5GJHcpcHPE5Zf78pB33pl0JSJNVLfzjWY7hIpy3yorfSuv8C11vyoQap6lVrfTTVB6VCHtrmNHGDsWnngCtm5NuhqR3KTAzSGnnuoXypw5M+lKRFqhTue7awtVe95SEhzLVeDmkGOO8dvXX0+2DpFcpXm4OaRHD5+18NFHSVci0g4yYAqOOtwcU1Lii96ISPwUuDmmTx8t4yiSFAVujunVC9avT7oKkdzU5MA1s3wzm2tm06P7+5nZLDNbZGYPmVlR+5UpbaVzZ9i+PekqRHJTczrcq4H5Ne7/ArgthHAAsB64pC0Lk/bRsaMCVyQpTQpcMxsEfAW4L7pvwEnAtGiXycDX2qNAaVuFhVBjaVERiVFTO9zbgR8CqdnDewEbQgjRVdz4FBhY3xea2Tgzm21ms0tLS1tVrLReQYGf/CAi8Ws0cM3sTGB1CGFOS94ghDAphDAyhDCyuLi4JS8hbSgvD6qqGt9PRNpeU058OA74qpl9GegIdAfuAHqaWUHU5Q4CtNpqBlDgiiSn0Q43hHB9CGFQCKEEOA/4RwjhAmAGcE6021jgiXarUtqMWUackCOSlVozD/da4Htmtggf072/bUqS9qTAFUlOs9ZSCCE8Dzwf/XkxcFTblyTtKU0vZiqSE3SmmYhITBS4OUYdrkhyFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFpUuCaWU8zm2ZmH5jZfDM7xsx6m9nfzOzD6LZXexcrIpLJmtrh3gE8G0IYChwGzAeuA54LIRwIPBfdFxGRBjQauGbWAzgBuB8ghFAWQtgAjAEmR7tNBr7WXkWKiGSDpnS4+wGlwG/NbK6Z3WdmXYC+IYQV0T4rgb71fbGZjTOz2WY2u7S0tG2qFhHJQE0J3AJgBHBXCOEIYCt1hg9CCAEI9X1xCGFSCGFkCGFkcXFxa+sVEclYTQncT4FPQwizovvT8ABeZWb9AaLb1e1ToohIdmg0cEMIK4FPzOzg6KHRwPvAk8DY6LGxwBPtUqGISJYoaOJ+VwJTzKwIWAxchIf1w2Z2CfAxcG77lCgikh2aFLghhLeAkfU8NbptyxERyV4600xEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmDQpcM1svJnNM7P3zGyqmXU0s/3MbJaZLTKzh8ysqL2LFRHJZI0GrpkNBK4CRoYQDgXygfOAXwC3hRAOANYDl7RnoSIima6pQwoFQCczKwA6AyuAk4Bp0fOTga+1fXkiItmj0cANISwHbgaW4UG7EZgDbAghVES7fQoMrO/rzWycmc02s9mlpaVtU7WISAZqypBCL2AMsB8wAOgCnN7UNwghTAohjAwhjCwuLm5xoSIima4pQwonA0tCCKUhhHLgUeA4oGc0xAAwCFjeTjWKiGSFpgTuMmCUmXU2MwNGA+8DM4Bzon3GAk+0T4kiItmhKWO4s/CDY28C70ZfMwm4FviemS0C9gLub8c6RUQyXkHju0AI4b+A/6rz8GLgqDavSEQkS+lMMxGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYKXBGRmChwRURiosAVEYmJAldEJCYK3BwzahRcfXXSVYjkJgshxPdmZqXAx7G9oTTHviGE4qSLEMlmsQauiEgu05CCiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjFR4IqIxESBKyISEwWuiEhMFLgiIjH5/xjX3J71f+OyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9e59652358>"
      ]
     },
     "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 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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}