summaryrefslogtreecommitdiff
path: root/notebooks/plot_OT_1D.ipynb
blob: 9b2b44682225d47267509f1297d06fd69295b928 (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
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
{
 "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 make_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": [
    "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": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f8f3d878710>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAADCCAYAAABnjpSEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU5fXA8e8hYQmCIIsbi4AsJSEhkAgIiGBBQCy4CyKLSwVbFNtKReWHFJdqtaKtiKXuKyhuWBcUhYKgYAKJbCKrEuoCCAGEQJb398eZiSFkmSST3JnJ+TzPfWa7M/fMZDLnvrs45zDGGGPKqobXARhjjAlPlkCMMcaUiyUQY4wx5WIJxBhjTLlYAjHGGFMulkCMMcaUS7TXARTWpEkT16pVK6/DMMYYA6Smpu52zjUt6rGQSyCtWrUiJSXF6zCMMcYAIvJNcY9ZFZYxxphyCSiBiMggEdkoIptFZHIRj9cWkbm+x1eISCvf/TVF5DkRWSMiG0Tk9uCGb4wxxiulJhARiQJmAoOBWGCEiMQW2u06YK9zri0wA3jAd//lQG3nXDyQBIzzJxdjjDHhLZA2kG7AZufcVgARmQMMA9YX2GcYMM13fR7wmIgI4IATRCQaiAGOAvuDE7oxJliys7PJyMggKyvL61CMR+rUqUPz5s2pWbNmwM8JJIE0A3YUuJ0BdC9uH+dcjohkAo3RZDIM+A6oC/zBOfdTwNGZclmwANavh4kToYa1cpkAZGRkUL9+fVq1aoWe+5nqxDnHnj17yMjIoHXr1gE/r7J/XroBucDpQGvgTyLSpvBOInKDiKSISMquXbsqOaTI9sQTcMEF8Mc/wlVXwZEjXkdkwkFWVhaNGze25FFNiQiNGzcucwk0kASyE2hR4HZz331F7uOrrmoA7AGuAj5wzmU7534ElgHJhQ/gnJvtnEt2ziU3bVpkd2NTCudgyhS48UYYPBjuuQfmztXrmZleR2fCgSWP6q08f/9AEsgXQDsRaS0itYDhwPxC+8wHxviuXwZ84nShkW+B83zBnQD0AL4qc5SmVDffDPfeC9deC2+9BXfeCc8/D0uXQp8+cPCg1xEaU7J7772XuLg4EhISSExMZMWKFV6HdJx69eoB8L///Y/LLrus2P327dvH448/XuJr9ezZE4DFixdz4YUXlimOt956i/Xrf2mGnjp1KgsXLizTawSFc67UDbgA+BrYAtzpu286MNR3vQ7wGrAZWAm08d1fz3f/OrTRfVJpx0pKSnKmbLZvd07EuRtucC4v79jH3nzTOXDuX//yJjYTHtavX+/p8ZcvX+569OjhsrKynHPO7dq1y+3cubPCr5udnV3h1yjohBNOCGi/bdu2ubi4uCIfKxzTokWL3JAhQ8oUx5gxY9xrr71WpucEoqjvAZDiivm9DqgNxDn3nnOuvXPuTOfcvb77pjrn5vuuZznnLnfOtXXOdXO+HlvOuYO+++Occ7HOuQeDkvXMMZ55Ri9vvx0Kl0KHDYO4OHjqqaqPy5hAfffddzRp0oTatWsD0KRJE04//XQAPv74Y7p06UJ8fDzXXnstR3wNe61atWL37t0ApKSk0LdvXwCmTZvGqFGj6NWrF6NGjSI3N5dbb72VTp06kZCQwD//+U8AUlNTOffcc0lKSmLgwIF89913x8W1bds2zj77bOLj45kyZUr+/du3b6dTp04ArFu3jm7dupGYmEhCQgKbNm1i8uTJbNmyhcTERCZNmsTixYs555xzGDp0KLGxOgrCX5oB2L9/P0OGDKFDhw6MHz+evLy84/aZN28eY8eOZfny5cyfP59JkyaRmJjIli1bGDt2LPPmzSv187rrrrvo2rUr8fHxfPVVxSuDQm4qE1M2ubmaQPr3h6KmEBOB667TRvW1a8H3nTemWLfcAmlpwX3NxER45JHiHz///POZPn067du3p3///lx55ZWce+65ZGVlMXbsWD7++GPat2/P6NGjmTVrFrfcckuJx1u/fj2ffvopMTExzJo1i+3bt5OWlkZ0dDQ//fQT2dnZ3HTTTbz99ts0bdqUuXPncuedd/L0008f8zoTJ07kxhtvZPTo0cycObPIYz3xxBNMnDiRkSNHcvToUXJzc7n//vtZu3Ytab4PcvHixaxatYq1a9cW2ctp5cqVrF+/njPOOINBgwbxxhtvFFtF1rNnT4YOHcqFF1543D6lfV5NmjRh1apVPP744zz00EM8+eSTJX6OpbFOnmHu44/h2281SRRn1CioWdNKISZ01atXj9TUVGbPnk3Tpk258sorefbZZ9m4cSOtW7emffv2AIwZM4YlS5aU+npDhw4lJiYGgIULFzJu3Diio/V8uVGjRmzcuJG1a9cyYMAAEhMTueeee8jIyDjudZYtW8aIESMAGDVqVJHHOvvss7nvvvt44IEH+Oabb/KPW1i3bt2K7SLbrVs32rRpQ1RUFCNGjODTTz8t9T0WpbTP65JLLgEgKSmJ7du3l+sYBVkJJMw99RQ0agQXXVT8Pk2aaFXWCy/A/feDr5bAmCKVVFKoTFFRUfTt25e+ffsSHx/Pc889R5cuXYrdPzo6Or+qp3D30xNOOKHEYznniIuL47PPPis1rtJ6J1111VV0796dd999lwsuuIB//etftGlz3GiFEmMqfAz/7YL3B2OQp7+KMCoqipycnAq/npVAwtiePdrjatSo0pPC9dfr/vML958zJgRs3LiRTZs25d9OS0vjjDPOoEOHDmzfvp3NmzcD8MILL3DuuecCWqefmpoKwOuvv17saw8YMIB//etf+T+YP/30Ex06dGDXrl35CSQ7O5t169Yd99xevXoxZ84cAF566aUiX3/r1q20adOGm2++mWHDhvHll19Sv359Dhw4EPD7X7lyJdu2bSMvL4+5c+fSu3dvAE455RQ2bNhAXl4eb775Zv7+xb1+SZ9XZbAEEsZefBGOHi25+sqvf39o0cKqsUxoOnjwIGPGjCE2NpaEhATWr1/PtGnTqFOnDs888wyXX3458fHx1KhRg/HjxwNw1113MXHiRJKTk4mKiir2ta+//npatmxJQkICnTt35uWXX6ZWrVrMmzeP2267jc6dO5OYmMjy5cuPe+6jjz7KzJkziY+PZ+fOwsPf1KuvvkqnTp1ITExk7dq1jB49msaNG9OrVy86derEpEmTSn3/Z511FhMmTKBjx460bt2aiy++GID777+fCy+8kJ49e3Laaafl7z98+HAefPBBunTpwpYtW/LvL+nzqgyivbRCR3JysrP1QErnHCQkQEwMrFwZ2HPuugvuvhu2b4eWLSs1PBNmNmzYQMeOHb0Ow3isqO+BiKQ6544bAA5WAglba9fqdu21gT/nmms08bzySuXFZYypPiyBhKlPPtHLIUMCf06rVjomZNGiSgnJGFPNWAIJU4sWwZlnartGWfTrB59+CtnZlROXMab6sAQShvLyYMkS8A28LZO+feHnn8GamYwxFWUJJAylp8PeveVLIP4efVaNZYypKEsgYWjxYr0sTwJp0gTi4395DWOMKS9LIGFo0SJo2xaaNy/f8/v2hWXLdAyJMaFgz549JCYmkpiYyKmnnkqzZs3ybx+tpC/qqlWr+OCDDwLat3fv3vnzWg0cOLDEQYIPP/xwiaPGr7nmGjZu3EhOTg4NGzasUMxvvvkmDz7o3Ry1lkDCTG6utn/061f+1+jXDw4dgi++CF5cxlRE48aNSUtLIy0tjfHjx/OHP/wh/3atWrVKfX5ubm6Zj1mWBFLQggULqF+/frGPl5RAcnNzeeaZZ+jQoUOZjwvHx3zxxRcHNFCxslgCCTNpabrCYHmqr/z69NFLawcx4eA3v/kNSUlJxMXF5c8e6z97v+WWW0hISGDlypXMnz+fDh06kJSUxE033cRFvgniDh48yNixY+nWrRtdunThnXfe4fDhw0yfPp2XXnqJxMTE/KnQ/Q4dOsTll19Ox44dufTSS49JCM2bN2ffvn0cOHCAwYMH07lzZzp16sS8efOYMWMGP/74I+eccw79+/cvMs6CpRmAm2++mbi4OAYMGMCePXuAY0s833//PW3bti0y5ieffDJ/pt1t27bRr18/EhISGDBgQP7kkFdffTUTJ06kZ8+etGnT5pgpUSrKJlMMMxVp//Br3FhHsS9erMvgGnMML+ZzL8Fzzz1Ho0aNOHToEMnJyVx66aXUr1+fzMxM+vTpwyOPPMKhQ4do3749y5Yto2XLllxxxRX5z58+fTqDBg3i2WefZe/evXTv3p0vv/ySqVOnsnbtWh4pIq7HHnuMk046iQ0bNrB69WqSk48fiP3ee+/RqlUr3n//fQAyMzNp0KABf//731m6dCkNGzYkJyfnmDgLy8zMpFevXvzjH/9g6tSp3H333UXuBxATE3NczAWnY//d737H9ddfz8iRI5k9eza33HJLfmL88ccfWbZsGWvWrOGKK67InyqloqwEEmYWLYL27cG31k659eun7SC+tWaMCVkzZsygc+fOnH322WRkZOTP/VSrVq38H8L169fToUMHzjjjDEQkfwp2gA8//JB7772XxMRE+vXrR1ZWFt9++22Jx1yyZAlXX301AF26dCEuLu64fRISEvjggw+YPHkyy5Yto0GDBkW+VsE4C4uOjubyyy8HtKRQ3mncAVasWMHw4cMBGD16NEuXLs1/7KKLLkJESEhIKHZOr/KwEkgYycnRNc5935EK6dsXHn1U59E655yKv56JIF7N516EhQsXsmTJEj7//HNiYmLo3bt3fnVSTExMqVOtg07d/tZbb3HmmWcec38g64qUpGPHjqSkpPDee+8xefJkBg8ezB133HHcfoHGCb9M317SVPXlUbvAdN3BnP/QSiBhZPVq2L+/Yg3ofn366GqF1g5iQllmZiaNGjUiJiaGdevW8UUxPT9iY2PZuHEjO3bswDnH3Llz8x8bOHBg/jK2AKtXrwaKnxIdoE+fPrz88ssApKenFznV+86dO6lXrx6jRo3iT3/6E6tWrSr1dQvLycnhjTfeAODll1/On8a94FT1BdtnSnrtHj168OqrrwLw4osv0sff2FmJLIGEEf8JUzCm92/USNtBKngSZkylGjJkCIcOHSI2NpYpU6bQvXv3IverW7cujz32GP379yc5OZmGDRvmVynddddd/Pzzz8THxxMXF8e0adMAOO+880hPT6dLly7HNaJPmDCBPXv20LFjR+6+++4iF7ZKT0/nrLPOIjExkfvuuy+/9HHDDTfQv39/+vfvX+r7a9CgAUuXLiUuLo5PP/00f931SZMm8eijj9K1a1f27t2bv39JMc+cOZPZs2eTkJDA3LlzmTFjRqnHryibzj2MjBgBy5fDN98E5/XGj4c5c3RUe4AlbBOhImE694MHD1KvXj2cc4wbN474+Hhuuukmr8MKKzadewRLSYEiOoOUW3KydgkusB6NMWFr1qxZJCYmEhsby+HDh/ntb3/rdUgRzxrRw8TevbB5c9nW/yiNPxmlpOjIdmPC2aRJkzwdVFcdWQkkTPja54JaAomL07XUrcbQGFMelkDChP9HPikpeK9Zs6aO77IEYiC43TtN+CnP398SSJhISYE2bbT3VDAlJ0Nqqq4xYqqvOnXqsGfPHksi1ZRzjj179lCnTp0yPc/aQMJESgp06xb8101Ohpkz4euv4Ve/Cv7rm/DQvHlzMjIy2LVrl9ehGI/UqVOH5mWc4tsSSBjYvRu2b4cbbwz+a/urxFJSLIFUZzVr1qR169Zeh2HCTEBVWCIySEQ2ishmEZlcxOO1RWSu7/EVItKqwGMJIvKZiKwTkTUiUrYyksE3IDWoDeh+HTtCTIy1gxhjyq7UBCIiUcBMYDAQC4wQkdhCu10H7HXOtQVmAA/4nhsNvAiMd87FAX2B7KBFX034E0jXrsF/7eho6NLll2MYY0ygAimBdAM2O+e2OueOAnOAYYX2GQY857s+D/i16Kxg5wNfOufSAZxze5xzZV/5pZpLSYF27aCMi5cFLDlZuwmXY00eY0w1FkgCaQbsKHA7w3dfkfs453KATKAx0B5wIrJARFaJyJ8rHnL1E+wR6IUlJ+sKhV99VXnHMMZEnsruxhsN9AZG+i4vFpFfF95JRG4QkRQRSbFeIMf64QfYsaPyEwhYO4gxpmwCSSA7gRYFbjf33VfkPr52jwbAHrS0ssQ5t9s5dwh4DziuJt85N9s5l+ycS27atGnZ30UEq8wGdL/27aFePUsgxpiyCSSBfAG0E5HWIlILGA7ML7TPfGCM7/plwCdORyQtAOJFpK4vsZwLrA9O6NVDSorOlFvEbNJBExWlDfSWQIwxZVFqAvG1aUxAk8EG4FXn3DoRmS4iQ327PQU0FpHNwB+Byb7n7gUeRpNQGrDKOfdu8N9G5EpJgQ4doH79yj1OcrIug52TU7nHMcZEjoAGEjrn3kOrnwreN7XA9Szg8mKe+yLaldeUQ1pa1Sw526ULZGXBxo06yaIxxpTG5sIKYT/9pA3onTtX/rH8x0hPr/xjGWMigyWQEPbll3pZFQnkV7+CWrUsgRhjAmcJJISlpellYmLlH6tmTa268h/TGGNKYwkkhKWnwymn6FYVOne2EogxJnCWQEJYenrVVF/5de6sAxe//77qjmmMCV+WQEJUdjasW1c11Vd+/mNZKcQYEwhLICHqq6/g6NGqL4GAJRBjTGAsgYQo/494VSaQk06CFi0sgRhjAmMJJESlpUHt2joKvSolJlpPLGNMYCyBhKj0dOjUSRd8qkqdO+to9Kysqj2uMSb8WAIJQc5VfQ8sv86ddWGpdeuq/tjGmPBiCSQEff897NpVtT2w/KwnljEmUJZAQpC/DcKLEkibNro2iLWDGGNKYwkkBPnP/hMSqv7YNWpAfLyVQIwxpbMEEoLS06FVK2jY0JvjJyZqDM55c3xjTHiwBBKC0tK8qb7y69wZMjPhm2+8i8EYE/osgYSYw4fh66+9TyBg1VjGmJJZAgkxa9dCXp43PbD84uN1HXZLIMaYklgCCTFe9sDyO+EEaNfOemIZY0pmCSTEpKdD/fraiO4lWxvEGFMaSyAhJj1du+/W8Pgv07kzbN0K+/d7G4cxJnRZAgkheXmaQLxs//Dzx+Bfl90YYwqzBBJCtm+HAwe8bf/ws55YxpjSWAIJIV6sAVKcZs2gUSNLIMaY4lkCCSHp6dr20amT15FoN15rSDfGlMQSSAhJS9Pus3Xreh2J6twZ1qzR6d2NMaYwSyAhJFQa0P0SE3Vk/KZNXkdijAlFlkBCxL592ogeCu0fftaQbowpiSWQEOHvLhtKCaRjR11S10akG2OKElACEZFBIrJRRDaLyOQiHq8tInN9j68QkVaFHm8pIgdF5NbghB15QqkHll/t2ppErARijClKqQlERKKAmcBgIBYYISKxhXa7DtjrnGsLzAAeKPT4w8D7FQ83cqWnQ5MmcPrpXkdyLP/aIMYYU1h0APt0AzY757YCiMgcYBiwvsA+w4BpvuvzgMdERJxzTkQuArYBPwct6gjkXwNExOtIjtW5M7zwgq7R3rSp19GYgGRlwcaNsGOHbpmZembSogWccQa0bh16XzQTlgJJIM2AHQVuZwDdi9vHOZcjIplAYxHJAm4DBgDFVl+JyA3ADQAtW7YMOPhIkZOj07j//vdeR3K8gg3p/ft7G4spweHD8MEH8Npr8M47cPBg8fu2aQNXXKFbYqIlE1Nuld2IPg2Y4Zwr4dsMzrnZzrlk51xy02p4mvv113DkSGi1f/hZT6wQd+QIzJgBzZvDJZfAhx/CVVfB3Lnw+eeQkaHJ5Ouv4eOPYdYsHWz04IPQtSv06AFLlnj9LkyYCqQEshNoUeB2c999Re2TISLRQANgD1pSuUxE/gY0BPJEJMs591iFI48gq1frZSiNAfFr2lRrP/wxmhDhHMyZA3fcof2/zz8f/vQnOO887TpXWLt2up13HowfD7t3w6uvwn33wbnnwm9+Aw88oL0mjAlQICWQL4B2ItJaRGoBw4H5hfaZD4zxXb8M+MSpc5xzrZxzrYBHgPsseRwvNRXq1IHYwl0TQkRSksZoQsRPP8FFF2lJo0EDWLBAt/PPLzp5FKVJE/jd77Rk8te/wn//q2cwjz+uycmYAJSaQJxzOcAEYAGwAXjVObdORKaLyFDfbk+hbR6bgT8Cx3X1NcVLTdWqokD/96taUpK2yR444HUkhs8+gy5d4P334eGHYdUqTRzlVbcuTJ6s0w30768NcZdfriNbjSlFQD9Zzrn3gPcK3Te1wPUs4PJSXmNaOeKLeHl5Wj00apTXkRQvKUlPStPToXdvr6OpxmbNgptv1t5Uy5bBWWcF77VPPlkb3x9+GG6/XRPTu+9alZYpkY1E99jmzXpm37Wr15EUzx+bVWN5xDmYNk2rnAYP1jOOYCYPvxo14NZbYelS7dV1zjmwcmXwj2MihiUQj/l/lJOSvI2jJKefDqeeagnEE7m5MGEC/OUvcO218MYb2u5RmXr00BJOgwba6P7RR5V7PBO2LIF4LDVVpwyJi/M6kpJZQ7oHcnNh9Ght2P7zn+HJJ6uuoaxNG00ibdvCkCHw1ltVc1wTViyBeGzVKkhIgJo1vY6kZElJ8NVX8LPNJ1A1nIMbb4SXX9autg88UPUD/k49FRYv1j/+lVfCwoVVe3wT8iyBeMg5TSChXH3ll5SkDf42oLAKOKcljn//G+68Uxu1vdKwIbz3HnTooF2HP/vMu1hMyLEE4qEtW3SaonBJIGDVWFXir3+Fhx7SLrV33+11NHDSSTrC/bTT4IIL7CzC5LME4qFwaED3O/10OOUUSyCV7vnntdRx9dXwj3+EzjxVp56qVVj16mkS2Vl4MgpTHVkC8VBqKtSqFfoN6KC/Y127WgKpVJ9+Ctdfrz2fnn5au9WGkjPO0Oqs/fth6FBrEDOWQLy0ahXEx2sSCQdJSbB+PRw65HUkEWjrVrj4Yp1qfd680O1VER+vc3ClpWkPsbw8ryMyHrIE4pFwakD38zek+5ffNUGSmamTGebmwn/+o20OoWzIEPj733VMypQpXkdjPGQJxCPbtsHeveGXQMCqsYIqL0/nsfn6a3j9dZ0xNxxMnAg33KAN/q+95nU0xiOWQDwSTg3ofs2b6/TulkCC6P77f5mDql8/r6MJnAj8859w9tlwzTWwYYPXERkPWALxSGqqVnN36uR1JIETsRHpQfXhh1oFdNVVOl1JuKlVS0sfJ5yg7Tf793sdkalilkA84p+Vu3ZtryMpm+7ddfldm9q9grZvhxEjtAve7Nmh0123rJo109UPN2/WkoitJVKtWALxQHa2TnLas6fXkZRdz55abb9ihdeRhLGjR3U98pwcbYg+4QSvI6qYvn11qpU33tDldU21YQnEA2lpkJUFvXp5HUnZ9eihJ8vLl3sdSRi77Tb44gt45pnwaTQvzR//qNVYt91mZxfViCUQD/h/fMOxBHLiiToUwBJIOb31FjzyiC4MdcklXkcTPCLw1FPa0+LKK7WLoYl4lkA8sHy5Duo9/XSvIymfnj21DcfGkJXR9u3aTpCUBH/7m9fRBN9JJ2l7yP/+Z+0h1YQlkCrmnC6zEI6lD7+ePbXDzbp1XkcSRrKzYfhwzbqvvhp+vScC1a2bJse339a5vExEswRSxXbs0Hnowj2BgFVjlcmUKdo28NRTulhTJJs4UefK+vOfdfldE7EsgVQx/49uODag+7VpozPzWgIJ0Icf6ln5uHFw2WVeR1P5RHQyyJNP1vYQ6/MdsSyBVLHly7XXZny815GUn4iWQiyBBOD773Wqkk6dqlcX18aN4aWXdNGbcBwkaQJiCaSKLVumg/GqamnrytKzp44d++EHryMJYXl5OmPtgQM6g21MjNcRVa0+fWDqVF3j5IUXvI7GVAJLIFXo4EFdzC2c2z/8/O/BVjgtwUMPwUcfabfdcFj0pTJMmQLnngu/+x1s2uR1NCbILIFUoS++0Bm7IyGBdO2qUyFZNVYxVqzQlQUvvxx++1uvo/FOVBS8+KJ+WYYPhyNHvI7IBJElkCrk/7E9+2xv4wiGOnUgOdkSSJH27dMfy2bNwnueq2Bp3lxH3a9aBbff7nU0JogsgVSh5cu1JqNhQ68jCY6ePSElxU4qj+Gc9rbasQNeeSVy/tgVNXQo3HSTdiR4912vozFBYgmkimRnw9Kl0Lu315EET+/emjw+/9zrSELIU0/pQMF77omMomYw/e1vkJgIY8boYCgT9gJKICIySEQ2ishmEZlcxOO1RWSu7/EVItLKd/8AEUkVkTW+y/OCG374+Pxz7Yxz/vleRxI8/fppb7IFC7yOJESsWaNn2f376yA6c6w6dbQ3WlaWroGSk+N1RKaCSk0gIhIFzAQGA7HACBGJLbTbdcBe51xbYAbwgO/+3cBvnHPxwBig2vblW7BA2xN//WuvIwmeE0/Uk2xLIGgXuyuu0CqrF1+EGla4L1KHDvDEE7BkCfzlL15HYyookG95N2Czc26rc+4oMAcYVmifYcBzvuvzgF+LiDjnVjvn/ue7fx0QIyIROglQyT78UKdCb9DA60iCa+BAbRvdtcvrSDz2+9/Dxo06eO6UU7yOJrRdfTVcey3cey8sXOh1NKYCAkkgzYAdBW5n+O4rch/nXA6QCTQutM+lwCrn3HFNriJyg4ikiEjKrgj8Jdq9WxubI6n6ys//nj76yNs4PPXsszpYbupUOK/a1tKWzT//CbGxMHIkfPed19GYcqqScraIxKHVWuOKetw5N9s5l+ycS27atGlVhFSlFi7UzjkDB3odSfB17aqzVlTbaqwvv9RBcn37wv/9n9fRhI+6dbWzwcGD2uXZ2kPCUiAJZCfQosDt5r77itxHRKKBBsAe3+3mwJvAaOfclooGHI4+/FCXSkhO9jqS4IuK0jbjDz+shss/ZGbCpZdqu8crr+iHYQIXG6vjZJYsgTvu8DoaUw6BJJAvgHYi0lpEagHDgfmF9pmPNpIDXAZ84pxzItIQeBeY7JxbFqygw4lzenbev3/k/r4MHKhzBq5Z43UkVcg5GDsWtm3TM+lTT/U6ovA0cqSW4B58UNdUN2Gl1ATia9OYACwANgCvOsOjHI8AABBVSURBVOfWich0ERnq2+0poLGIbAb+CPi7+k4A2gJTRSTNt50c9HcRwtat0wXaIrH6ys/fDlKtqrEeekiXp33wwcga3OOFhx/WhajGjoWvv/Y6GlMG4kKs3iE5OdmlpKR4HUbQ/P3vcOut8O230KJF6fuHq06d9CS8WnSqWbgQBg2Ciy/W0kd1n6okGL79VhvUTjlFB03Vr+91RMZHRFKdc0VWwFtn9Ur24YfQsWNkJw/QEtbSpXDokNeRVLItW3S8x69+pYsmWfIIjpYtNRlv3Kjrp+TleR2RCYAlkEp0+LC2D0Zy9ZXfwIFw9CgsXux1JJXowAEYNkyTxvz5dpYcbOedp3Nlvf023HWX19GYAFgCqUTvvquzNgwZ4nUkle+cc/T39PXXvY6kkuTl6ZnxV1/pmXKkr2vulQkTdJDhPffAa695HY0phSWQSvTSS9ou0K+f15FUvpgYbRKYN0+TZsS54w49M3744ciajybUiMDjj+tUz2PGwMqVXkdkSmAJpJLs3QvvvadjpCK1+25hI0fC/v36viPKE0/AAw/A+PE6WaKpXLVrw5tv6tnXb36jXaVNSLIEUklef13bBK66yutIqs5558HJJ2vJK2L85z86z9WQITr9hjWaV42TT4b339d1EAYPhp9+8joiUwRLIJXk5ZehXbvIHH1enOhoLXG9+64uyhf2UlPhyit1DYs5c/QNmqrToYNWG27bpp0XIrJuNLxZAqkEO3dqb6Srrqp+J6xXXaWLTIX9oOL163WsR9OmWgqpV8/riKqnc86B556DTz/V7tPZ2V5HZAqwBFIJ5szRmS6qU/WVX7ducOaZWgILW1u3woABWuL46CM47TSvI6rehg+HmTPhnXdg9GjIzfU6IuNjCaQSvPSSVl21b+91JFVPRBPnJ5/oFC5hJyNDe1llZWnyaNfO64gM6HxZDzygZ2fjxtlAwxBhCSTINmyA1au1R1J1NXKklsDmzvU6kjLyJ489e3Rir06dvI7IFPTnP8OUKbru/IQJlkRCgCWQIHviCa35uPJKryPxTocOcNZZ+lmETW3D1q1a3/7dd9r7pzr1fggn06drIpk1C667Loy+YJHJEkgQ/fgj/PvfOmC5ulebT5qkE6u++abXkQTgq6+gTx8dxPLJJ9Crl9cRmeKIwP33w7RpuhLkyJHWsO4hSyBB9OijWnV+221eR+K9Sy7RNqD77gvxhaZSUuDcc/VHaPFiK3mEAxGdK+vBB7We9OKLdWVDU+UsgQRJZiY89hhcdplW4VR3UVEwebK2B4XsOiFvvqklj7p1ddbL+HivIzJlceutWk/6/vv6d9xZeKFUU9ksgQTJ449rDcjtt3sdSegYORKaN9dSSEhxTheEuvRSSEjQ9Scs64enceN0nM6mTdC9O6SleR1RtWIJJAgOHdJZqAcNgi5dvI4mdNSqpW0hS5fqFhJ+/llXvps0SYuLixbpIkYmfA0erAMNRXR1yLAehBReLIEEwZNPwq5dOmGrOdb110OTJnDvvV5Hgo4u79YNXnhBG2HnzNFphE3469xZZ+7t0kWLvuPG2dQnVcASSAXt3KnteX37ai9Qc6y6dfVkf8ECD9cKcQ6eeUb7Fu/erctE3nUX1LCvf0Q57TQtUd52G8yeDT166MAsU2nsP6gCnIPf/lbnfpo92+toQtcf/qDLXd94o5bUqtSOHTqT7rXXagJZvRr696/iIEyViY7Wbr7/+Y8ODE1MhL/+FXJyvI4sIlkCqYCnn9YOIPffbzNelKRmTZ0PLzNTk0iVdOvNy9OsHhcH//2v9rH+5BM4/fQqOLjx3JAhsG6dridyxx1aGrEG9qCzBFJO336rZ9Z9++qsCqZknTrBX/6i1ViVPsXJf/+r4znGjdPLNWvg5putyqq6OeUUXSLztdf0H7ZrV60y+OEHryOLGPYfVQ5Hj2pHHue0FGK/S4G59Vbtafn738P27ZVwgHXrtGtu377a1vHyy/Dxx7Z+eXV32WWwcSNMnKij19u21V4d+/d7HVnYs5++Mjp6VJclWLRIF6hr3drriMJHdLRWZeXl6TrxQUsiqak69L1TJ20gv+ce/cEYMaL6LchiinbSSdrXft06XTpzyhQ44wztTGGrHZabJZAy8CePt9/W5DF2rNcRhZ8OHXSW9H37KphEjh7VurB+/bSaatEimDpVX/DOO617rila+/b6D/zFF/rdmT4dWrbU9e5Xr/Y6urBjCSRAhw4dmzys3aP8kpNh4UJNIn37wubNAT7ROZ27atIkaNFCFxravl3XifjmG21kady4EiM3ESM5WZfNXLNG/7Gff17bSLp312klrJ0kIOJCbKa75ORkl5KS4nUYx1i4EG64QZdmfuwxrcM3FZeaqgv/HTkCd9+t7dzHLTt+9CgsX67d3V57Tf8I0dHay+bGG/UFrBHKVNTevTrAdPZsreaqUUMn2bz4Yhg4ULtZVtPqUBFJdc4VOcuoJZASfP+9zm317LNa8v33v3XONhM8GRm62Nw77+hJ4awZWSRLqiaNpUu1aurgQU0av/61ni1edBE0auR16CYSOacJ5LXXtIp040a9v1UrPVnp2VO3apRQKpxARGQQ8CgQBTzpnLu/0OO1geeBJGAPcKVzbrvvsduB64Bc4GbnXIlzs3qdQHJztR323//WHzXndP2aqVOhTh3Pwoo8ubnatXLTJtzadWyf/yX7l31Jx5w11ELXd8hr05YaAwfA+edrw+eJJ3octKl2tmzRH4QFC3S6/8xMvb9xYx2kmJCgszh37KhJpVGjiEssFUogIhIFfA0MADKAL4ARzrn1Bfb5HZDgnBsvIsOBi51zV4pILPAK0A04HVgItHfOFbuMWFUmkLw8XYBu82adkHXZMj3x3bMHmjaFMWO06soGCZZBVpZWB+zbp11pf/xRt+++0+LGjh2aOLZtO3YhoFNPJTs2gXTpwvObzmbutz3IrH0KZ52l6zv16qX/oy1b6iSNxlS5vDydGuWzz2DFCkhPh7Vr4fDhX/Zp2FC7ZrZooVvz5joe5eSTdWvUSPdp2FDXPAgDFU0gZwPTnHMDfbdvB3DO/bXAPgt8+3wmItHA90BTYHLBfQvuV9zxKpJAdm38ifT73iU3l/wtOxuOZsPRI5B1BH4+qBOy7t+vv2/ZBWY4OO1UTRYJnSGpaxH18ZWppL9DwcdKu17wsvCWl/fLZcGt4AeWm6vTPuTk6Ifn344c0S0rS7fDh7VnwaFDcODAL9uRI0W/hxo1dK6i5s31H+vMM/XDPvNMiI3Vf64Cb2XFCq1FWLYMVq36JdfUqKEvcdpp+r/YqBE0aKCdrvxbzZq6RUfr/2hUlD5P5PjNr7jrxpRE8nKp98MWTvz+a+r/sJkTf9hEvV3bqLs3g7p7dlD70L5in5tdpx7ZdeqTXac+OXXqkVOrLrm1YsitGUNuzTrk1qxNXs065EbXwkXVJM+3uRpR5EVF42pE4SRKL2tE4aQG1KiBkxo4EfBf1qxJj9nXlf89lpBAAvmJbAbsKHA7A+he3D7OuRwRyQQa++7/vNBzmxUR4A3ADQAtW7YMIKSi7U79hv7Pjy738/net4XK1ONeEdFf35o19XS/Zk2oXfvYrW5dqF9ff/jr1/9la9hQ+9yfdJJOw+s/82rSJOCMLKIzT/ToobcPH9YksmmTFly2bdNOMrt2aRV1Zqbuc/hwiK9+aCJQFNDetx2vLj/TlF2czI+cwg+cxN787cSs/dTPOsCJ7KceB4nhMDEcoC4/UJsj1OYIMRyhFkepSXb+ZTRlWwf+EDFQgQRSkqo8xy6Wc242MBu0BFLe12k7NJbdn28mOpr8rVatMOqkU9KpbyCny/7rBS/9W8FT8IKn5P7r/lN1/+0QEhPzSzVWSZzTTlvZ2b8UoHJzjy1sFSyQFXxeUdeNqbgTfFurMj8zDzjs247hr0XIzUXycn+5LHg/v9Q6iMAZFX0bxQgkgewEWhS43dx3X1H7ZPiqsBqgjemBPDdoatarTZPuZ1bWy5sQJ/JLAcmYyCVoycf7NpRATjW/ANqJSGsRqQUMB+YX2mc+MMZ3/TLgE6eNK/OB4SJSW0RaA+2AlcEJ3RhjjJdKLYH42jQmAAvQlPe0c26diEwHUpxz84GngBdEZDPwE5pk8O33KrAeyAF+X1IPLGOMMeHDBhIaY4wpVkm9sEKrtdQYY0zYsARijDGmXEKuCktEdgHfVPBlmgC7gxBOOLPPQNnnYJ+Bn30O5fsMznDONS3qgZBLIMEgIinF1dlVF/YZKPsc7DPws88h+J+BVWEZY4wpF0sgxhhjyiVSE8hsrwMIAfYZKPsc7DPws88hyJ9BRLaBGGOMqXyRWgIxxhhTySIqgYjIIBHZKCKbRWSy1/FUFRFpISKLRGS9iKwTkYm++xuJyEcissl3eZLXsVY2EYkSkdUi8h/f7dYissL3nZjrm88toolIQxGZJyJficgGETm7un0XROQPvv+FtSLyiojUqQ7fBRF5WkR+FJG1Be4r8m8v6h++z+NLEela1uNFTALxrZw4ExgMxAIjfCsiVgc5wJ+cc7FAD+D3vvc+GfjYOdcO+Nh3O9JNBDYUuP0AMMM51xbYiy6vHOkeBT5wzv0K6Ix+HtXmuyAizYCbgWTnXCd0Dr/hVI/vwrPAoEL3Ffe3H4xOcNsOXY9pVlkPFjEJBF02d7Nzbqtz7igwBxjmcUxVwjn3nXNule/6AfQHoxn6/p/z7fYccJE3EVYNEWkODAGe9N0W4Dxgnm+X6vAZNAD6oBOc4pw76pzbRzX7LqATxcb4lpeoC3xHNfguOOeWoBPaFlTc334Y8LxTnwMNReS0shwvkhJIUSsnHrf6YaQTkVZAF2AFcIpz7jvfQ98Dp3gUVlV5BPgzuhYP6KqY+5xz/oWLq8N3ojWwC3jGV5X3pIicQDX6LjjndgIPAd+iiSMTSKX6fRf8ivvbV/g3M5ISSLUnIvWA14FbnHP7Cz7mW58lYrvciciFwI/OuVSvY/FYNNAVmOWc6wL8TKHqqmrwXTgJPbtuDZyOLglYuFqnWgr23z6SEkiVrn4YakSkJpo8XnLOveG7+wd/kdR3+aNX8VWBXsBQEdmOVl+eh7YFNPRVY0D1+E5kABnOuRW+2/PQhFKdvgv9gW3OuV3OuWzgDfT7Ud2+C37F/e0r/JsZSQkkkJUTI5Kvrv8pYINz7uECDxVcKXIM8HZVx1ZVnHO3O+eaO+daoX/7T5xzI4FF6CqZEOGfAYBz7ntgh4h08N31a3RBt2rzXUCrrnqISF3f/4b/M6hW34UCivvbzwdG+3pj9QAyC1R1BSSiBhKKyAVoPbh/5cR7PQ6pSohIb2ApsIZf6v/vQNtBXgVaojMcX+GcK9zAFnFEpC9wq3PuQhFpg5ZIGgGrgaudc0e8jK+yiUgi2pGgFrAVuAY9Waw23wUR+QtwJdpDcTVwPVq/H9HfBRF5BeiLzrr7A3AX8BZF/O19yfUxtHrvEHCNc65Mq/lFVAIxxhhTdSKpCssYY0wVsgRijDGmXCyBGGOMKRdLIMYYY8rFEogxxphysQRijDGmXCyBGGOMKRdLIMYYY8rl/wED42j/NaledQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 460.8x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZhcVZmH36+6sydsIUJYQsIORsPSIrhgBFzABXRQNCObYhgUFRBZXEbAZVRUQEF2AXGBAREC4ygCMqAgmCCbIKAsEQwQVpElS/eZP75bXVXdXem6dbdafu/z3KdTp86990sl/fXb557zHQshIIQQIn9KRQcghBDdihKwEEIUhBKwEEIUhBKwEEIUhBKwEEIUhBKwEEIUhBKwEF2Cmb3ZzO4rOg5RQQlYdD1mNs/MFprZv8xsiZn9r5m9KeE1Hzaz3dKKsYH7BTPbdFV9Qgg3hhC2aPL6D5vZcjNbe0j7n6J7z2zmut2OErDoaszsCOBk4OvAOsAM4AfAnkXGlTZm1pvCZR4CPlx1zdcAE1O4bteiBCy6FjNbHTgB+GQI4bIQwoshhBUhhCtDCJ+L+owzs5PN7B/RcbKZjYveW9vMrjKz58zsGTO70cxKZnYhnsivjKz6qBHuPdfMHjWzo8zsyci89zKzPczs/uh6n6/qv4OZ3Rzda4mZnWpmY6P3boi63RHdb5+q6x9tZo8D55XbonM2ie6xXfR6PTNbamZzV/GRXQjsV/V6f+BHTX34AlACFt3NTsB44Ber6PMFYEdgG2AOsAPwxei9zwKPAtNwe/48EEII+wKLgfeEECaHEL5V59rrRvdfH/hP4GzgI8D2wJuBL5nZrKhvP3A4sHYU967AJ/Ab7hz1mRPd7+Kq668FbATMr75xCOFvwNHAj81sInAecEEI4fpVfBZ/AFYzs63MrAf4EPDjVfQXo6AELLqZqcBTIYSVq+jz78AJIYQnQwhLgeOBfaP3VgDTgY0ic74xxCuusgL4WghhBXARnlxPCSG8EEL4M3APnvQJISwKIfwhhLAyhPAwcCbwllGuPwB8OYSwLITw8tA3QwhnA38Fbon+Hl9oIOayBb8NuBd4rIFzRB2UgEU38zSw9ijjo+sBj1S9fiRqAzgRT2BXm9mDZnZM3PuHEPqjP5cT5BNV778MTAYws82j4Y7Hzeyf+Jh1zQOxEVgaQnhllD5nA7OB74cQljUQ84XAPOAANPyQGCVg0c3cDCwD9lpFn3/gv8KXmRG1EZnqZ0MIGwPvBY4ws12jfmmXGTwd+AuwWQhhNXy4w0Y5Z5UxmNlk/AHkucBxZrbWaEGEEB7BH8btAVzWQNxiFSgBi64lhPA8PvZ6WvQAbKKZjTGz3c2sPG77M+CLZjYtmoL1n0Tjnmb2bjPb1MwMeB4fpx2IznsC2DjFcKcA/wT+ZWZbAocMeb+Z+50CLAwhHAT8D3BGg+d9DNglhPBizPuJISgBi64mhPAd4Aj8wdpS4O/AocDlUZevAguBO4G7gNuiNoDNgGuAf+E2/YMQwm+j9/4LT9zPmdmRKYR6JP6r/wv4sMHFQ94/Drggut8HR7uYme0JvJNKIj8C2M7M/n20c0MIfwshLIwRu6iDqSC7EEIUgwxYCCEKQglYCCEKQglYCCEKQglYCCEKIo0CHaIDWHvttcPMmTOLDkOItmLRokVPhRCmNXu+ErAAYObMmSxcqJlFQsTBzB4ZvVd9NAQhhBAFoQQsRLfx4oug+f8tgRKwEN3AAw/AfvvBJpvA5MkwdSq8/e1w2WVKxgWiBCxEJ7NyJRx3HMyeDZdfDttuCyecAHvvDQ89BP/2b/Cud8HixUVH2pXoIZwQncrKlbDvvnDRRfDhD8N3vgPTp9e+f+qp8KUvwZvfDNdfD7Nm1b2cSB8ZsBCdSH+/DzlcdBF84xvw05/WJl+A3l447DC44QZ44QWYO9etWOSGErAQncjxx8PPfubJ9+ijV913223h2ms9Ce+5J7w8bPMMkRFKwEJ0GtddB1/9Khx44OjJt8y227ol33UXHHFEtvGJQZSAhegkli6Fj3wENt8cvv/9eOe+853wuc/BGWfApZdmE5+oQQlYiE7iqKPgqafg4oth0qT453/1q9DXB4ceCs8/n358ogYlYCE6hZtvhvPPh8MPhzlzmrvG2LFw+unw5JM+jiwyRQlYiE6gv9+tdb31fFpZEvr6YP58+N734O6704lPjIgSsBCdwPnnw223+VzfyZOTX+9rX4PVV3ebFpmhBCxEu7Nsma9u22EH2GefdK45dSp84QtwzTXwf/+XzjXFMJSAhWh3zj3XlxJ/5Stglt51DznEF2986UuqF5ERSsBCtDMvv+zDBW96E7ztbelee8IE+Pzn4cYb3YRF6igBC9HOnHMO/OMf6dtvmY9/HDbcEL78ZVlwBigBC9GurFgB3/622+/cudncY9w4X013883w+99nc48uRglYiHblkkt87Peoo7K9z4EH+kO5b30r2/t0IUrAQrQjIXhC3Gorr+ebJRMn+hzjK6+Ee+/N9l5dhhKwEO3Ib34Dd9zhtRtKOXwbf/KT/lDuxBOzv1cXoQQsRDty0kmw7rowb14+95s2zYcifvITX6YsUkEJWIh24/774Ve/8nm648bld99PfQqWL4ezzsrvnh2OErAQ7capp8KYMV6vIU+23NI38jz9dJ+BIRKjBCxEO/HCC173YZ99fAgibz79aZ93fNll+d+7A1ECFqKduOACT8Kf+lQx9999d9/a/nvfK+b+HYYSsBDtQgjwgx/A617nhXeKoFSCT3wCbrrJZ2GIRCgBC9Eu3HCDz8M95JBi4zjgABg/3rcuEolQAhaiXTjjDFhjjfRKTjbLWmt5DD/+sQ+HiKZRAhaiHXjiCfj5z90+J04sOhq38H/9y+cFi6ZRAhaiHfjhD33q18EHFx2Js8MOsM02PiVNVdKaRglYiFZnYADOPtsrnm25ZdHROGbwH/8Bd94Jt9xSdDRtixKwEK3ONdfAQw+1jv2WmTcPJk3SyrgEKAEL0eqcdRasvTa8731FR1LLlCmehC+6CJ57ruho2hIlYCFamccfhyuu8IdvedZ9aJSDD/ZtkfQwrimUgIVoZc47D1au9K2BWpHtt4fttoMzz9TDuCZQAhaiVal++Lb55kVHU5/58+Guu/QwrgmUgIVoVVr14dtQ9DCuaZSAhWhVzjyzNR++DaX6YdzzzxcdTVuhBCxEK7JkCSxY0LoP34Yyf74exjWBErAQrUirP3wbSl+fHsY1gRKwEK1Gf78/fHvrW1v74dtQ5s/3lXF/+EPRkbQNSsBCtBpXXw0PP+xLfduJefNg8mS3YNEQSsBCtBpnnAHrrAN77VV0JPGYMgX23RcuvhieeaboaNoCJWAhWonFi+Gqq+BjH4OxY4uOJj4HHwyvvOJbJ4lRUQIWopU45xx/iNUuD9+GMmcO7LSTW7wexo2KErAQrcLy5f7wbffdYebMoqNpnkMOgfvvh2uvLTqSlkcJWIhW4bLLvPjOoYcWHUkyPvABX0By6qlFR9LyKAEL0Sqceipsuim84x1FR5KM8eN9StqVV/psDlEXJWAhWoE//Ql+/3v45Cd96/d2pzyFTjsnr5IO+JcWogM49VTfbPOAA4qOJB023NCn0Z1zji9RFiOiBCxE0Sxd6jUU9t3Xt53vFD71KXj6adWHWAVKwEIUzemnw7JlcNhhRUeSLm95i++cfNJJmpJWByVgIYrklVfgtNPgXe9qnR2P08IMjjgC7rkHfv3roqNpSZSAhSiSn/wEnnzSE1Unss8+sN568N3vFh1JS6IELERRDAx4YtpmG6981omMHetjwb/5DdxxR9HRtBxKwEIUxVVX+a/nn/2s/7reqRx8sFdJ++Y3i46k5VACFqIIQoCvfQ1mzYIPfajoaLJlzTV9efLFF8Nf/1p0NC2FErAQRXDddXDrrXD00dDbW3Q02XP44TBmDHzrW0VH0lIoAQtRBF//OkyfDvvvX3Qk+TB9Onz0o3D++fDoo0VH0zIoAQuRNzfc4AZ85JFeN6FbOOooH3r5r/8qOpKWQQlYiDwJAb74RTfCQw4pOpp8mTnTC82ffbaK9EQoAQuRJ1dfDTfe6El4woSio8mfL37Riw0df3zRkbQESsBC5EXZfjfaCA46qOhoimGDDeATn4Af/Qj+8peioykcJWAh8uKii2DhQjjuuPbc7y0tjj0WJk3yGSBdjhKwEHnw0kuecLbdFvbbr+hoimXaNPjCF2DBArjmmqKjKRQlYCHy4Nvfhr//HU4+uTMKriflM5/xRSiHHw4rVxYdTWHof4IQWbN4sS/D3Xtv2HnnoqNpDcaPhxNPhLvvhjPPLDqawlACFiJLQqhMNzvxxGJjaTXe/37YbTcfE37ssaKjKQQlYCGy5OKL4Ze/9LoP7bzVfBaY+Z5xK1f6zIguLNquBCxEVjz1FHz60/C613lJRjGcTTbxOcELFsAllxQdTe4oAQuRBSF47YPnn/eNKXt6io6odTn8cOjr852U//73oqPJFSVgIbLgtNPgyiu9+tdrX1t0NK1Nby/89KewfDl85CPQ3190RLmhBCxE2tx2mxfa2WMPH4IQo7PZZv5D64Yb4IQTio4mN5SAhUiTJUtgzz3hVa+C887r7J0u0ma//bw85wknwKWXFh1NLnRBJWghcuKVV+B974NnnoHf/96TsGic8qyI++/3ZLzxxrDddkVHlSkyYCHSYPlyX2hxyy1w4YW+0aaIz/jx8Itf+HLl3Xfv+II9SsBCJGXlSpg3D/7nf+D0032BgWieddbxsp1msOuu8Le/FR1RZigBC5GEF1/0hPvzn8NJJ/lUKpGcLbbwQj3LlsGb3wy33150RJmgBCxEszzxBOyyi5vvaafBYYcVHVFnMXs2XH+9z6HeeWf49a+Ljih1lICFaIZrroE5c+Cuu3zM8hOfKDqizmT2bLj5Zq+ctvvu8J//2VHV05SAhYjDCy/AEUfA298Oa63lD93e+96io+psNtjAZ5Xstx985Sswdy78+c9FR5UKSsBCNEJ/P/z4x7DVVj7We/DB8Mc/wmteU3Rk3cHkyb6l/YUXwr33+iyTI4/0ehttjBKwEKvipZd8QcWrXw377utze2++2Wc7TJpUdHTdx0c+Avfd5/8W3/2uD00ccww88kjRkTWFErAQQ1m5En77W/jkJ2H99b2oztixvjpr4ULYcceiI+xu1l4bfvhDL+b+rnd5neWNN/Y//+hH8NxzRUfYMBa6sAanGE5fX19YuHBh0WEUwwsv+DSnP/7RaxFcf71XMZs40ZcV/8d/+FQoLStuTRYvhrPPhgsu8GpqPT1eAvStb4XXv94rra23Xib/fma2KITQ1/T5SsACOiwB9/fDyy/78MELL/jx7LO+RPjJJ3362OLF/mvr/ffDo49Wzp01yyf/v/OdfmiYoX0IwX+ILlgA110Ht95aqay2+uo+t3jWLJgxA6ZP9wUfU6fCmmv6+1Om+L/3xIleoa2BhK0ELFKhb+LEsHDTTfO74dD/d+XX1e0hDD8GBvzo768cK1f6sWKFT9wfbZqSmX8DzpjhVbi23NJLRvb1wbrrpvv3FMXx0ktwxx2waJEvab7vPv+hu3ix/z9ZFWYwbpwPPY0Z4wm5p6dybLwxXHtt4gSsYjzCGTcO8kzAMNwwyq+r281qj56eytfqY8wYP8aP97/LhAl+TJnix5pr+rSxV73K6wz06r9+xzNxIuy0kx/VhODjxE88AU8/7b8Z/fOf/pvSSy/5b0+vvOJJevly/8G+cmXlB/7AQGqFlvS/UDibbAKXXVZ0FEJkj5n/QF5zzaIj0SwIIYQoCiVgIYQoCD2EEwCY2QvAfUXHMQprA62+9KkdYoT2iLMdYtwihDCl2ZM1BizK3JfkaW4emNlCxZgO7RBnu8SY5HwNQQghREEoAQshREEoAYsyZxUdQAMoxvRohzg7PkY9hBNCiIKQAQshREEoAQshREEoAXc5ZvZOM7vPzP5qZscUHU8ZM9vQzH5rZveY2Z/N7DNR+1pm9hszeyD6Wvh6UjPrMbM/mdlV0etZZnZL9JlebGZjC45vDTO71Mz+Ymb3mtlOrfY5mtnh0b/z3Wb2MzMb3wqfo5n90MyeNLO7q9pG/OzM+V4U751mtt1o11cC7mLMrAc4Ddgd2Br4sJltXWxUg6wEPhtC2BrYEfhkFNsxwLUhhM2Aa6PXRfMZ4N6q198ETgohbAo8C3yskKgqnAL8KoSwJTAHj7VlPkczWx/4NNAXQpgN9AAfojU+x/OBdw5pq/fZ7Q5sFh3zgdNHvXoIQUeXHsBOwK+rXh8LHFt0XHVivQJ4G75ab3rUNh1fQFJkXBtE34S7AFcBhq/e6h3pMy4gvtWBh4geuFe1t8znCKwP/B1YC18cdhXwjlb5HIGZwN2jfXbAmcCHR+pX75ABdzfl//hlHo3aWgozmwlsC9wCrBNCWBK99TiwTkFhlTkZOAoYiF5PBZ4LIZSLEhf9mc4ClgLnRcMk55jZJFrocwwhPAZ8G1gMLAGeBxbRWp9jNfU+u9jfT0rAoqUxs8nAz4HDQgj/rH4vuGYUNo/SzN4NPBlCWFRUDA3QC2wHnB5C2BZ4kSHDDS3wOa4J7In/sFgPmMTwX/tbkqSfnRJwd/MYsGHV6w2itpbAzMbgyfcnIYRyseInzGx69P504Mmi4gPeCLzXzB4GLsKHIU4B1jCzcp2Voj/TR4FHQwi3RK8vxRNyK32OuwEPhRCWhhBWAJfhn20rfY7V1PvsYn8/KQF3N38ENoueNo/FH3wsKDgmwJ8oA+cC94YQvlv11gJg/+jP++Njw4UQQjg2hLBBCGEm/tldF0L4d+C3wN5Rt6JjfBz4u5ltETXtCtxDC32O+NDDjmY2Mfp3L8fYMp/jEOp9dguA/aLZEDsCz1cNVYxMUQPvOlrjAPYA7gf+Bnyh6Hiq4noT/qvdncDt0bEHPsZ6LfAAcA2wVtGxRvHOBa6K/rwxcCvwV+ASYFzBsW0DLIw+y8uBNVvtcwSOB/4C3A1cCIxrhc8R+Bk+Lr0C/23iY/U+O/wB7GnR99Jd+KyOVV5fS5GFEKIgEg1BtOokfiGEaAeaNuBoEv/9+NzMR/HxxA+HEO5JLzwhhOhckuyIsQPw1xDCgwBmdhE+laRuAl577bXDzJkzE9xSJGXpUnj2Wdh889r2NW+bFf9iQ7eVB7DSkJdDt56P3h/SbuVrlUq1145eV96v3rJ+yLVKPTWvB8/p6am95mB7+Xz/GkpDYuipjSWUqv5uPUPaBl9HX6NzQ4+N/H70daB36Hnlr9SeD4ToVgND+gx9Xe5Xeb/2WgO9te8P+1r11xzWtzfUuXaofT/6StRO9Nqi19br06Z7eqKvveVp1DBmjE/9HdPTD8DYXv86rqf81d+f0LsCgPHR10k9y729J3rduwyAiSVvn9LzCgCTo6+rlV4evOfE0rKozd+bMvjVrzXFQvTaP5DJpfEAlNZ9YIRvgsZJkoBHmnT8+qGdzGw+viyPGTNmsHBhoh08REI+9zk47TQY+s/wttIH4l+s+renclIL0TdSlBzDQPQNVxry/kBt8iz/JmYD0fvlxBa9Lic6q3yfQmnIteiPvvZE50Sh9Eft5URcpn+g5qVFI3KB2vZyIh6MDSivYbKob+V1mXLf8jUZ8n70brTMYGDYd2K5ZxjWVoraBuq8Hkr50xmI+pWifgMj9q4TX524Ktce+d5Df7+uvPYz+xmJlHdKS+NyUSJmoLwuJEriSS+b8PxRCSGcFULoCyH0TZs2LevbiVFYuRLGjCk6CiEEJPvZ0NKT+MXILFsGY7OoKVW24ZxM2PtEf8jZhKvjy8+Eh58tE45JBiZcpAG37CR+UZ+XXoKJE4uOQggBCX4mhBBWmtmhwK9x1fhhCOHPqUUmMuFf/4LJkzO8QU4mDCOMC+dlwjBsXDh7E67u3bkmDCPZcAubcEIShRJC+CXwy1QiEbnw9NMwdWrRUQghIPUfLaLVefxxePWrc7hR1iYM9WdIyIRHfD2U1jThytltYcIJUTGeLmJgAB56CGY1MeVXCJE+LfAzQOTFgw/6LIihizAyJSsThtHnCmdkwjD6XOG0TRgamSucrglXR1mPtE24ti0nE87oko0gA+4ifvc7//qGNxQbhxDCkQF3EQsWwKteBVttVcDNUzZhiLFqLmUThsZXzaVlwhBn1Vw6JuxtjY0Ly4SbQwbcJSxZ4gl4//1rn2MJIYpDBtwlHH+8S+jHP16ng1ltbYesSMmE/VIx60fIhEd8Pez6VX+OO0MiqQlXRz/8deeZsFyoC/i//4Mzz4TDDoPNNis6GiFEGSXgDufPf4YPfAA23hhOOKHoaIQQ1WgIooO5+27YZRfo7YVf/hImTRrlhMFhgTYYioDmS1kmHYqApktZNjsUUdsn6pHxUEQlCg1FZIUMuAMZGIAf/AB23NGT7/XXwxZbjHqaECJnZMAdxkMPwUEHwXXXwTveAeecAxts0MCJVqqy0DYwYUhe1L2NTNj7MKRP1EMmPIT2MWEZcIfwwAM+w2GLLeCPf4Szz4b//d8Gk68QohBkwG3O7bfDN74Bl1ziO118/ONwzDGw4YajnzuM8h5r7WDC1X3yNmFIbXujxk0Y0tveqHVNeFX370QTVgJuQx59FC66CH76U/jTn2DKFN/r7bDDYN11i45OCNEoSsBtwjPPwKWXetK94QYXxte9Dk46CQ44ANZYI9n1rWSDltkOJuxNKW/02aAJQwqlLGObMKS/0eeqTbi6bShZmXDttTvfhJWAW5SVK33n4l//Gq6+Gm65Bfr7fYz3uONg3jzYdNOioxRCJEEJuIV45BFPtldfDddcA8895yLY1+fjuu9/P2y7bUUO06ZslO1gwt6U7kafjZqwnxOFk5MJe1uZfEy4XlvNPQYjSsuEK3G1hQknpPUi6hJWroS77oKbbvLj5pt9Chn4zIX3v9+nke26q7YQEqJTUQLOiWefhT/8oZJwb7kFXnzR35s+Hd74RvjMZ+Dtb4ctt8zOcutSNQ+4LUy4Kq78TRiy2vK+ngnXtpXJ1oT9DtlubzTchIfH1ckm3DqRdBDPPeezExYtgttu86/33+/v9fTAnDlw4IFeGP0Nb4AZMwpIuEKIwlECTsgzz1SSbPnr3/5WeX/DDWH77WG//TzZvu51GW8L3ywlq5heXBOG7G14qAnXxJG3CUNmG33WMWE/J1kpy7gm7FfMZ6PPmq2X6sTViSZcfARtQnlDyzvuqD0efrjSZ+ZMT7Yf/ah/3W47mDatqIiFEK2OEvAIvPiiPyCrTrR33QUvvODvl0q+seXrXw+HHOKJdrvtYK21io07MWVzjGvCkN+4cPX1s9ryfjQTrr5WTibsfZqrpNa8CVfaOmPL+3omDEWlwq5OwCH4qrKhVvvAA5Xv89VWg9e+1ocQ5syBbbaBV78aJk4sNnYhRPvTNQl42TK4557hyfaZZyp9Zs3yJDtvnn+dM8eHFbrhAZmZVTa8jGvCVX1aYoZE1iYMqW9vNJoJV8eXnwkPP1smnC4dmYBfesmL1Cxa5KvJbrsN/vIXn3sLMGECvOY18G//Vkm0r32t264QQuRF2yfgl192ky0n24UL3XTLMrPOOv5A7L3vrSTbTTcdubRrV1MqVQwrrgl7Y02fTjZh7xP9IS8Thsw2+qxvwtW9O9eEobgZEm2XgP/xD99k8vrrfTHD3Xd7jQTwGQd9ffC+9/nX7beH9dbrjiEEIUT70fIJeMkST7blo7ygYbXVfMudd7/bE21fny/hVbJtErNB44trwt7UwqvmUjZhiF8/QiZcHUGrmXDl7LxNuOUScAi+o8OPfuQFae67z9tXWw123hnmz4e5c302goYRhBDtTMsk4KeeggsvhB/+0IcVJkzwHX0POqiScHtbJtoOpFSqGF5ME4Y2qR+RlglDdvvM1TFhiF8/IqkJQ3a7a9Qz4eoo65G2Cde25WvCo17VzDYEfgSsg8d5VgjhFDNbC7gYmAk8DHwwhPBs3AAWL4Yjj4TLL4cVK2CHHeDMM2GffWD11eNeTQgh2odG0vpK4LMhhNvMbAqwyMx+AxwAXBtC+IaZHQMcAxwd5+ZPPQVve5uP8x56qC/hnT077l9BpIGZVZ7yxzVhiF8/op1NGFLbXaNRE4b49SOSmjCkt7tGoybsbc1WUms/Ex71aiGEJcCS6M8vmNm9wPrAnsDcqNsFwPXESMDLlvkDtMWLfaz3jW+MGbkQQrQ5sdK5mc0EtgVuAdaJkjPA4/gQxUjnzAfmA8yYMWOw/fHHfc5uX59XCBMFU7JBe4ttwtB8JbU2NGFoon6ETLjm9bDrV/05qx2XV7Xbc1GV1Eqjd3HMbDLwc+CwEMI/q98L/j9zxO+kEMJZIYS+EELftKrSYBttBGed5XN5DzjAhyOEEKKbaCgBm9kYPPn+JIRwWdT8hJlNj96fDjwZ9+Yf/Sh8/evws5/5gokPfhB+9avKwgohhOhkGpkFYcC5wL0hhO9WvbUA2B/4RvT1imYCOPZYeM974NxzfRraJZf4gooDDoC99tJ839yw0uCvy7GHIqD5UpZtOBThl0pa1D3mUARkuOX9yEMRtX2iHhkPRVSiaJehiGQ0YsBvBPYFdjGz26NjDzzxvs3MHgB2i143xezZcNJJ8NhjnoBnz4avfc3Hh6dO9QT9ne94UR3ZsRCiU2hkFsTvGHlcHmDXNIMZNw723tuP8hLk3/7Wv151lfdZfXVfETd3rs+cmDMHxo9PM4oupXpLopgmDI0vW+4IE4bkRd3bwIS9D0P6RD1kwqnQsmvLpk+HD3/YD3A7Lhfhuf56uPJKb+/tdWMu14Po6/NSk+PGFRW5EEI0Rssm4KGsv74XSp83z18/9hjcemulBOUvfuHjyABjxngSLifk7bbzXSxkyqug1MPgz/WYJux94hXwaWsThvS2N2rUhCGzLe/rmzAkL+re+ia8qvs3UsoyCW2TgIey/vpedvJ97/PXIcAjj1QS8qJF8N//7VPdwP/vbrFFpSZweXuhddct7u8ghOhu2jYBD8XMtw+aOdPHkMGT8oMP+sO78hZEv/udT3sr86pX1SblOaJHqlUAABYfSURBVHNgyy3doruKklG2r9gmDE2XsmxLE67uk5MJQ/wCPslNGNLb3qgxE65uG0pWJlx77WZKWTZPxyTgkTCDTTbx4wMfqLQ/8wzceWft3nDf/74vjwYYOxa23np4Yp46tZi/hxCiM+noBFyPtdbyWRRz51baVqzwYu/VSflXv4ILLqj0WX/9ytBFOSlvtlltlcJ2xYvxlF/FNGFoupRlO5qwNzVZyrJJE/ZzonvnZMLeViYfE67XVnOPwYjSMuFKXM0U8ElCVybgkRgzxh/UvfrVlQd9AE88MXwn5auvrmzwOXmyJ+Ttt/eHfdtv72PNql0shBgNpYlRWGcdePvb/ShT3uL+9tvhT3/yB35nn+27MYMXk99mm0pC3n572GqrFh9X7ukZNKu4JuznxCzg08Ym7E0Ji7rHNmHIasv7eiZc21YmWxP2OzRXyrJ5Ex4eV14mrATcBOPGwbbb+nHggd7W3+/bJy1a5A/9Fi3y4YvTTvP3J0zwYvNveIMfO+2kMWUhuh0l4JTo6fEHd1tvDfvu620DA/DAA56Mb70VbroJTjyxMnyxxRaVhPyGN/jsi8LGk80GTSquCUMT9SOSmrBfrPG/XzPUM+GquPIzYchso886JuznxKsfkdSE/YrJirrHNeHavvFXzSVBCThDSiVPsltsURlXfukln6d8001+LFgA553n702dCrvtVhny2GCD4mIXQmSPEnDOTJzotSx23tlfh+CWfPPNXvfi6qvh4ov9va23riTjnXeGSZMyDKzaLOOaMDRfSa1ZE66OOW8Tro4jLxOuvlZOJux94tWPSG7ClbZ22vK+WTpgAlV7Ywabbw777w/nn+9LrO+8E779bTfgM86APfbwqXN77unJufywTwjR3siAWwwzr2PxmtfAZz8LL78MN94Iv/yll+pcsMCnvu21lw9r7LZbSrMrekrDt79p1ISh6UpqTZtwVZ/cTbjmnjmZMCSvpBbThKvjy8+Eh5/dySYsA25xJkzwIYiTT/YNTK+7Dj70IS/PuccevjjkuOPg6aeLjlQIERcl4Daipwfe+lafc/z443D55T6d7fjjfY+9I47wIYymKJXcfHpKbnfVR0+PzxM2w8zc4krmFdQGj6jNSn5EryvnlPwoXzN6Pfj+YBxDrhNhJautxeCNtUZcvnYehFBrxGGgZnw6DISalXOD7w8EPwYvE9yGBwb8KF83el1+3/tEx9BrDfRHh78e7N/f70f5muWjf8CP6B42ELCBqhj6q47oHBsYcBvuD9A/0uvo6B+IjoBF79W8Hx2llX5Uzqs+iI7o9YD/JlDqD5Sq3h/6utyv8r4f5euUVrrhVq4/wlG+19C+K82PIddOihJwmzJunI8JX3EF3HWXV4X73vdg1ixPxBonFqL1UQLuAGbP9v30HnjAH+addJKvxLvppsavEUpVlhrbhC13E66x4S4w4Wobzs2EBwow4YH2MuGkKAF3ELNm+fDEddd5caE3vcnHh/NYvSuEiI9mQXQgb32rT2U79FAfH37pJfjmN0eRwlJp8Gm4Df25PMrsCIhfPyLp7AhoYK5wC9WPSDw7AhLXFI47OwLi149IOjsCktcUjjs7ojrKeqy6klrzKAF3KFOm+LziSZN8+fPYsfDVrxYdlRCiGiXgDsbMiwG98gp8/euw++6+k/SI9Ay3nkZNGBqYK5y2CUP8+hHtbMKQoJJacyYMo88VTtuEoflKas2asLclqaTWPBoD7nDMfHbEjBlw0EGwfHnREQkhysiAu4DJk+GUU3z13OWXwwc/OEInsxoLhhgmDE1XUmvahKuv1QUmDCOMC8uEa143a8Lep/lKakmQAXcJ73mPL9Y488yiIxFClFEC7hJKJdhvP6+49vzzRUcjhAAl4K7iLW/x34xvvnn4e6F6cUR5IUa0SCKUrLGFGg0tW05poQY0vmy5AxZq+KUaW7ac2kKNOMuWU1qoEW/ZcjoLNZparDFQNeSTACXgLmKHHfzrbbcVG4cQwtFDuC5iyhTfZPShh0Z4s8cqD0zKD3safSgHzZeybPKhHIy+WKOjHsrB6Is1snooB6lveV/voZz3YUifqEfGD+Wqo2imlGUzyIC7jI028rKWQojikQF3GVOnwtKlw9tDqTTcUlrYhKvj6woThuRF3eOaMGS25X19E4ZGly13ggnLgLuMNdeEZ58tOgohBMQwYDPrARYCj4UQ3m1ms4CLgKnAImDfEILWWbU4EyfWqRVcPQYc04T9nMaWLadmwhC7gE9bm3B1n5xMGBpfrJGeCUPcAj5JTbi6bShZm3AcA/4McG/V628CJ4UQNgWeBT6WUkwiQyZM8H3mhBDF05ABm9kGwLuArwFHmE/I3AWYF3W5ADgOOD2DGEWKjBkDK0dYRhlKVuUN8UzY+zS4bDktE4amS1m2owl7U5OlLJs0YT8nundOJuxtZfIx4XptNfcYjGh4KcskNGrAJwNHUflEpgLPhRDK38qPAuunEpHIlN5eL9YuhCieUQ3YzN4NPBlCWGRmc+PewMzmA/MBZsyYETtAkS49PbXDqGVCTwkGLTZqa2ET9nNWPVe4k0zYmxIWdY9twpDVlvf1TLi2rUy2Jux3aL6UZRIaMeA3Au81s4fxh267AKcAa5hZ+a+6ATDifrwhhLNCCH0hhL5p06alELJIQqnkq06FEMUzqgGHEI4FjgWIDPjIEMK/m9klwN54Ut4fuCLDOEVKmI0sdKHHqP75DjFMGJovZdmkCfv9o3t1gwlXxZWfCUNa2xs1asJ+TmOr5tIyYb9i86Usk5DkOkfjD+T+io8Jn5tOSCJL6iVgIUT+xFoJF0K4Hrg++vODwA7phySypG4RsB6r2aDFacyEoYlVc0lNGJIXdW/WhP1iZMpQE66JIycTrr5WTibsfRpdNZeWCVfamlk1lwSthBNCiIJQLYguo54B184DLtOYCde05WXCkP5Gn42aMOQ3Llx9/ay2vK9nwpC8klpME66OLz8THn52XiYsAxZCiIKQAQvADbhMXBP2tsbmCqdnwpDZlvejmXBVn5aYIZGRCXuf6A95mTBkttFnfROu7p2vCcuAhRCiIGTAAoCBXhu21XajJux94q2aS27CkNmW96OZsDfW9OlEE4Y4q+a604STIgMWQoiCkAELwMeAy0YQ14Rr++Rjwh7z4Jvlk6KwsjVhb2rhVXNpmTBkt89cHROGxlfNpWXCkKySWhJkwEIIURAyYAEQ1YJw4pswJK2kFteEIcaquZRNGGKsmmtnE4b095kbxYSh8VVzaZkwJKuklgQZsBBCFIQMWAAQemCoGzRuwtB0JbVmTRgy23F5VBOGxlfNtbEJQ4xVc11qwkmRAQshREHIgAVQHgMe+af7aCYM1J0h0ZEmXH2tDjZhv1Q2Oy7XNWHIcMflkU24tk/UI1YlteaRAQshREEoAQshREFoCEIA5V8NV/2god5QxEhnZj4UAdlteT/aUAS09/ZGjQ5FQAYbfbbeUIT3YUifqEdDpSybRwYshBAFIQMWAAz0WNXkcpmwvz+yCUPjy5bb2oQhwy3v65gwZLblfX0ThmQFfJpHBiyEEAUhAxaAL8QYXmqvMROG+AV8kpqwn9PYsuW0Tbg6vo424eo+OZkwNL5YIz0ThnRKWcZHBiyEEAUhAxZA7RhwXBOGZgr4JDNh79PgYo20TRjS3+izBU3Ym9LZ3qhRE/ZzonvnZMLeViaeCSdFBiyEEAUhAxbAyGPAMuE6JgzZbXnfQibsTRlteV/XhCGrLe/rmXBtW5nGTDgpMmAhhCgIGbAAap9Cxzfh6rboGlmbMGS45f2qTdjPaXDVXDubcFVc+ZkwZLbRZx0T9nMaWzWXtgnLgIUQoiBkwAIgKsheS+MmDHFXzSU1YWhi1VxKJuz3ju6Vlwn7xciUoSZcE0dOJlx9rZxM2Ps0umqu1oSTIgMWQoiCkAELAAZ66v80Ht2EodFVc2mZcE1b3iYMqW/0OaoJQ37jwtXXz2rL+3omDKlvbzSaCVfH10wltSTIgIUQoiBkwAKA0GMMRB4a14Sr2/IyYW9rbK5w6iYMmW15X9eEq/q0xAyJjEzY+0R/yMuEIWElteZpyIDNbA0zu9TM/mJm95rZTma2lpn9xsweiL6umUpEQgjRJTRqwKcAvwoh7G1mY4GJwOeBa0MI3zCzY4BjgKMzilNkzEAvlKKf73FNeOS2bE3Y+8RbNZeeCUNmW97XM2FvrOnTiSYMcVbNtYIJJ2NUAzaz1YGdgXMBQgjLQwjPAXsCF0TdLgD2SikmIYToChox4FnAUuA8M5sDLAI+A6wTQlgS9XkcWGekk81sPjAfYMaMGYkDFtngtSCcuCbsfZqrH9GsCdf2yduEIbMt7+uYsDe18Kq5tEwYsttnro4JQ+Or5kaaj56ERq7SC2wHnB5C2BZ4ER9uGCT4pzniv34I4awQQl8IoW/atGlJ4xVCiI6hEQN+FHg0hHBL9PpSPAE/YWbTQwhLzGw68GRWQYrsqV4JF9eEvU+ySmrxTbgSad4mDI2vmkvLhCHGqrl2NmFIf5+5UUwYGl81N9J89CSMasAhhMeBv5vZFlHTrsA9wAJg/6htf+CKVCISQoguodFZEJ8CfhLNgHgQOBBP3v9tZh8DHgE+mE2IIg9GrgXhtKYJV0eSrwlDjFVzaZkwNL5qro1NGGKsmmsBE05KQwk4hHA70DfCW7umEoUQQnQhWgkngKFP+GtpRROuPjN3E4bMdlyua8LV1+pgE/ZLZbPjcl0ThkSV1JKQjkcLIYSIjRKwEEIUhIYgBBAtRR5lq+16QxHeVu+cbIYiRjqzo4cioL23N2p0KAJS3N4o+6GIpMiAhRCiIGTAAhiyFDmmCXtbc6Us29KEIbst7+uYMDS+bLmtTRgy2OhzFBOGZAV8EiADFkKIgpABCwBCb4DBsV2nUROG5ktZNmvCI8XXySZcHV9Hm3B1n5xMGOIX8EnLhGXAQghREDJgAZSXItdaZ6MmXNM3JxOG9Df6bNSE/Zx4BXwSmzCkt71RC5uwN6WzvVGjJuznRPduopRlEmTAQghREDJgATBkW3qZMNQ3Ye8Tr4BPYhOG9Df6bEET9qaMtryva8KQpIBPEmTAQghREDJgAUDoCVX2OdgafW1FE65ui67RwSbs58Qs4NOOJlwVV34mDGmUsmwGGbAQQhSEDFgA5XnATnuY8PC4Bq+RtQlDhlvej2zCfv/oXnmZsF+MTBlqwjVx5GTC1ddqpn5EAmTAQghREDJgAdQacJlGTRiar6TWvAlX4sjbhKGJVXNJTRiSF3WPa8KQ37hw9fWz2vK+nglDOpXUmkAGLIQQBSEDFk5PoJ7jjGbC3sfJy4Sr2/I24Zq2vEwY0tveqFETrurTEjMkMjJh7xP9oZlKagmQAQshREHIgIVTNQYc34Qh7gyJpCY8cls+Juxt8epHJDdhyGzL+3om7I01fTrRhCF+/YhVbWIbBxmwEEIUhAxYAGAjjAE3bsLVvfMxYe+T7u4ajZqw94m3ai65CUNmW97XMWFvauFVc2mZMCSqpJYEGbAQQhSEDFgAYL0DlH8exzVhGGlcOFsTrr5/3iZc2ycfE4b49SOSmjC0Sf2IpCYMCSupNY8MWAghCkIGLADo6Rmo+pkez4QhvfoR7WHClUjzMmGIXz8isQlD/PoRbWjC0ET9iIF03FUGLIQQBSEDFgD09FbMSyZce92RKhGnVVO4YROGzHZcrmvC1dfqYBP2SyWopJYAGbAQQhSEErAQQhREQ0MQZnY4cBD+W9hdwIHAdOAiYCqwCNg3hLA8ozhFxowZs5Kh/x1aeSii9tr17p3NUET1mR09FAHJi7q3w1AEJCtlmYBRDdjM1gc+DfSFEGbj/8ofAr4JnBRC2BR4FvhYOiEJIUR30OhDuF5ggpmtACYCS4BdgHnR+xcAxwGnpx2gyIcxPdUTy+OZcL02yM6EvS2t7Y3imfBIZ2ZuwpDdlvd1TBiaKODTjiYMiUpZJmHUq4QQHgO+DSzGE+/z+JDDcyGE8n/NR4H1RzrfzOab2UIzW7h06dJUghZCiE5gVAM2szWBPYFZwHPAJcA7G71BCOEs4CyAvr6+HH4MimYY2zvS0srWNWG/VrKi7jLh6F51TLg6vo424eo+TZSyTEIjHr0b8FAIYWkIYQVwGfBGYA0zK3+HbgA8lkpEQgjRJTQyBrwY2NHMJgIvA7sCC4HfAnvjMyH2B67IKkiRPeN6VlVcZNUmDPEL+CQ14Zq+OZvwSPFlbcJ+TrwCPolNGJIXdW8DE/amBKUsE9DIGPAtwKXAbfgUtBI+pHA0cISZ/RWfinZuKhEJIUSX0NAsiBDCl4EvD2l+ENgh9YhEIYzraUA565owNFvKsh1NGNLb3qhRE/Y+zZWybNqEIb3tjVrYhL0pQSnLBGglnBBCFISK8QgAJvSuiNF7+H+bZlfNyYRb14T9nJgFfNrRhKviaqaUZRJkwEIIURAyYAHA+FgGXKZbTbi6LbpG1iYMGW55P7IJ+/2je+Vlwn4xMmWoCdfEEc+EkyIDFkKIgpABCwAm9SQtZJdOJbVGTRjSL+reuAkPj2vwGhmZMKRQSS2uCUPyou5xTRjyGxeuvn6SSmoJkAELIURByIAFABN6mhkDHol8TNj7RFfO3YSr48jHhGva8jJhSH+jz9FMuKpPS8yQaKSSWgJkwEIIURAyYAHApN5lKV8xaxOGZiupJTXh6ra8TNjbElZSi23CkNmW9/VM2Btr+rS0CSdEBiyEEAUhAxYATCwtz+h/Q1YmXN07XxMeuS1bE/Y+zVVSa96EIbMt7+uYsDe18Kq5ESqpJUEGLIQQBSEDFgBM6Xml8qINTLhyZnXvfEzY+6S7u8ZoJlzbJx8Thvj1I5KaMLRJ/QhLx11lwEIIURAyYAHA5GoDLtPCJgzp7zPXqAnX3jsvE65EmpcJQ/OV1Jo2YYhfP6LoSmoJkAELIURByIAFAKuVXq7/pky4gXtnbcLVkeRkwpDZjst1Tbj6Wu1gwgmRAQshREEoAQshREFoCEIAMLHUwFJkDUU0cO9shiKqz+zooQhor+2NEiIDFkKIgpABCwBWK40wDa0eLWDC9dogexOuvXa9e6drwiOdmbkJQ3Zb3tcxYUhQyrIIE06IDFgIIQpCBiwAmBLHgMt0qQl7W9obfcqEB/9+7bjlfZPIgIUQoiBkwAKAKaWEWxKl/j9p1SYMSUpZJjNhv1Za2xs1ZsIjxZe1Cfs5SYu6xzRhSG97ozYwYRmwEEIUhAxYADDFAiS1YMjRhCG97Y3imXBN35xMGLLb8r6eCXuftLY3atCEIbst71vQhGXAQghREDJgAcCUUi8MRGrVBiYM2W15LxMuzoT9nGRF3dvJhGXAQghREDJgAcDk0nggmgssE47ej67aEiZc3RZdI2sThgy3vB/ZhP3+0b3yMmG/GEUgAxZCiIKwkGPmN7OlwCO53VDEYaMQwrSigxCim8g1AQshhKigIQghhCgIJWAhhCgIJWAhhCgIJWAhhCgIJWAhhCgIJWAhhCgIJWAhhCgIJWAhhCgIJWAhhCiI/wenKlFkVpAEFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "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": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXyU1dn/8c8Vwr7JJiKLICCoVASidWktrhVt61KtS921uFdw72Krtvpo1aJWRXDXqlCXR9H6uCLVX0E0URBcUVEE2ZRFBAWSnN8f1z0mspiQZObM8n2/XvdrMvfcM3NlIN+cnPvc51gIARERybyi2AWIiBQqBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFskTZvZjM3svdh1SewpgkYSZnWBmM8xslZktMLPRZrZZ8titZvZVsq0xs7XV7v9fBmoLZtbn+44JIbwcQuhXj/c40symmtlKM1uUfH2GmVnyuJnZ1Wb2RbJdnXpM6kYBLAKY2XnA1cAFQFtgF2Ar4DkzaxJCOC2E0CqE0Aq4Ehifuh9CGBavcmdmxfV8/nnADcA1wBZAZ+A0YHegSXLYcOBgYCCwA/Bz4NT6vG+hUwBLwTOzNsBlwNkhhKdDCGtDCB8DvwJ6AsfU4TWHmtlcM7swaU3ON7ODzewAM3vfzJaY2e+rHb+zmU0xs2XJsTeZWZPksZeSw6YnLe4jqr3+RWa2ALgrtS95Tu/kPQYn97c0s8VmNnQDtbYFLgfOCCE8HEJYEdwbIYRfhxBWJ4ceD1wXQpgbQpgHXAecsKmfjVRRAIvAbkAz4NHqO0MIXwFPAfvW8XW3SF63K/An4DY8zIcAPwYuMbNeybEVwEigI7ArsDdwRlLHHskxA5MW9/hqr98eb6kPX6f2D4GLgH+aWQvgLuCeEMKkDdS5K9AUeLyG72d7YHq1+9OTfVJHCmARD73PQwjlG3hsfvJ4XawFrgghrAXGJa9zQ9LCfAt4G/9znhBCWQjhlRBCedL6HgP8pIbXrwT+HEJYHUL4et0HQwi3AR8AU4EuwB828jrrff9mNjlpjX9tZqlfAK2A5dWetxxopX7gulMAi8DnQMeN9KN2SR6viy9CCBXJ16mAXFjt8a/xUMPMtjGzJ5OTf1/i/cw1Bf/iEMI3NRxzGzAA+Ee1roT16mSd7z+EsFsIYbPksVROfAW0qfa8NsBXQTN61ZkCWASmAKuBQ6vvNLNWwDDghQzUMBp4F+gbQmgD/B6oqWX5vcGX1H89cAdwqZm138ihqe//oBre7y2SFntiYLJP6kgBLAUvhLAcPwn3DzPb38wam1lP4F/AXOC+DJTRGvgS+MrM+gOnr/P4QmDrTXzNG4DSEMIpwL+BWzd0UAhhGf7932Jmh5lZazMrMrMdgZbVDr0XONfMuprZlsB5wN2bWJNUU6+hKyL5IoTwNzP7ArgW6I2H4WPAr7/nT/eGdD4wFrgQeAMYD+xV7fFLgXvMrDl+wm3R972YmR0E7A/8INl1LjDNzH4dQrh/3eOT739e8v73AiuBj/ATeZOTw8bgvwRmJPdvT/ZJHZm6b0RE4lAXhIhIJApgEZFIFMAiIpEogEVEItEoCAGgY8eOoWfPnrHLEMkpZWVln4cQOtX1+QpgAaBnz56UlpbGLkMkp5jZJ/V5vrogREQiUQCLFJqVK0Hj/7OCAlikEMyaBccdB717Q6tW0KED7LcfPPqowjgiBbBIPisvh0svhQED4LHHYNAguPxyOOwwmD0bfvlLOPBAmDMndqUFSSfhRPJVeTkceyyMGwdHHQXXXQddunz38ZtugksugR//GCZNgl69Nvpy0vDUAhbJRxUV3uUwbhxcdRU88MB3wxeguBhGjICXXoIVK2DoUG8VS8YogEXy0WWXwYMPevhedNH3HztoELzwgofwQQfB1+striFpogAWyTcTJ8Jf/wonnlhz+KYMGuSt5Bkz4Nxz01uffEsBLJJPFi+GY46BbbaBf/xj0567//5wwQVw663w8MPpqU++QwEskk8uvBA+/xzGj4eWLWs+fl1//SuUlMBZZ8Hy5TUfL/WiABbJF1OmwN13w8iRMHBgjYdvUJMmMHo0LFrk/ciSVgpgkXxQUeGt1i239GFl9VFSAsOHw403wsyZDVOfbJACWCQf3H03vP66j/Vt1ar+r3fFFdC2rbemJW0UwCK5bvVqv7pt553hiCMa5jU7dIA//AGefx7+85+GeU1ZjwJYJNfdcYdfSvyXv4BZw73u6af7xRuXXKL5ItJEASySy77+2rsLfvQj2Hffhn3t5s3h97+Hl1/2lrA0OAWwSC67/Xb47LOGb/2m/OY30L07/PnPagWngQJYJFetXQvXXuut36FD0/MeTZv61XRTpsB//5ue9yhgCmCRXPXQQ973e+GF6X2fE0/0k3J/+1t636cAKYBFclEIHojbbuvz+aZTixY+xviJJ+Cdd9L7XgVGASySi557DqZP97kbijLwY3zmmX5S7ppr0v9eBUQBLJKLRo2CLbaAo4/OzPt16uRdEfff75cpS4NQAIvkmvffh6ef9nG6TZtm7n3PPhvWrIGxYzP3nnlOASySa266CRo39vkaMql/f1/Ic/RoH4Eh9aYAFsklK1b4vA9HHOFdEJn229/6uONHH838e+chBbBILrnnHg/hs8+O8/7DhvnS9jfeGOf984wCWCRXhAC33AI77eQT78RQVARnnAGTJ/soDKkXBbBIrnjpJR+He/rpces44QRo1syXLpJ6UQCL5Ipbb4XNNmu4KSfrqn17r+Gf//TuEKkzBbBILli4EB55xFufLVrErsZb4V995eOCpc4UwCK54M47fejXqafGrsTtvDPsuKMPSdMsaXWmABbJdpWVcNttPuNZ//6xq3FmcNpp8OabMHVq7GpylgJYJNs9/zzMnp09rd+Uo4+Gli11ZVw9KIBFst3YsdCxIxxySOxKvqt1aw/hceNg2bLY1eQkBbBINluwAB5/3E++ZXLeh9o69VRfFkkn4+pEASySze66C8rLfWmgbDRkCAweDGPG6GRcHSiARbJV9ZNv22wTu5qNGz4cZszQybg6UACLZKtsPfm2Lp2MqzMFsEi2GjMmO0++rav6ybjly2NXk1MUwCLZaP58mDAhe0++rWv4cJ2MqwMFsEg2yvaTb+sqKdHJuDpQAItkm4oKP/m2557ZffJtXcOH+5Vxr7wSu5KcoQAWyTbPPgsff+yX+uaSo4+GVq28FSy1ogAWyTa33gqdO8PBB8euZNO0bg3HHgvjx8OSJbGryQkKYJFsMmcOPPkknHwyNGkSu5pNd+qp8M03vnSS1EgBLJJNbr/dT2Llysm3dQ0cCLvu6q14nYyrkQJYJFusWeMn34YNg549Y1dTd6efDu+/Dy+8ELuSrKcAFskWjz7qk++cdVbsSurn8MP9ApKbbopdSdZTAItki5tugj594Kc/jV1J/TRr5kPSnnjCR3PIRimARbLBG2/Af/8LZ57pS7/nutQQOq2c/L3y4F9aJA/cdJMvtnnCCbEraRjdu/swuttv90uUZYMUwCKxLV7scygce6wvO58vzj4bvvhC80N8DwWwSGyjR8Pq1TBiROxKGtZPfuIrJ48apSFpG6EAFonpm2/g5pvhwAOzZ8XjhmIG554Lb78NzzwTu5qspAAWien++2HRIg+qfHTEEbDllvD3v8euJCspgEViqaz0YNpxR5/5LB81aeJ9wc89B9Onx64m6yiARWJ58kn/8/y88/zP9Xx16qk+S9rVV8euJOsogEViCAGuuAJ69YIjj4xdTXq1a+eXJ48fDx98ELuarKIAFolh4kR49VW46CIoLo5dTfqNHAmNG8Pf/ha7kqyiABaJ4coroUsXOP742JVkRpcucNJJcPfdMHdu7GqyhgJYJNNeeslbwOef7/MmFIoLL/Sul//5n9iVZA0FsEgmhQB//KO3CE8/PXY1mdWzp080f9ttmqQnoQAWyaRnn4WXX/YQbt48djWZ98c/+mRDl10Wu5KsoAAWyZRU63erreCUU2JXE0e3bnDGGXDvvfDuu7GriU4BLJIp48ZBaSlcemlurvfWUH73O2jZ0keAFDgFsEgmrFrlgTNoEBx3XOxq4urUCf7wB5gwAZ5/PnY1USmARTLh2mvh00/h+uvzY8L1+jrnHL8IZeRIKC+PXU00+p8gkm5z5vhluIcdBnvsEbua7NCsGVxzDcycCWPGxK4mGgWwSDqFUDXc7Jpr4taSbQ49FPbZx/uE582LXU0UCmCRdBo/Hp56yud9yOWl5tPBzNeMKy/3kREFOGm7AlgkXT7/HH77W9hpJ5+SUdbXu7ePCZ4wAR56KHY1GacAFkmHEHzug+XLfWHKRo1iV5S9Ro6EkhJfSfnTT2NXk1EKYJF0uPlmeOIJn/1rhx1iV5PdiovhgQdgzRo45hioqIhdUcYogEUa2uuv+0Q7BxzgXRBSs759/ZfWSy/B5ZfHriZjFMAiDWn+fDjoINh8c7jrrvxe6aKhHXecT895+eXw8MOxq8mIApgJWiRDvvkGDjkEliyB//7XQ1hqLzUq4v33PYy33hoGD45dVVqpBSzSENas8Qstpk6F++7zhTZl0zVrBv/7v3658rBheT9hjwJYpL7Ky+Hoo+Hf/4bRo/0CA6m7zp192k4z2Htv+PDD2BWljQJYpD5WrvTAfeQRGDXKh1JJ/fXr5xP1rF4NP/4xTJsWu6K0UACL1NXChbDXXt7yvflmGDEidkX5ZcAAmDTJx1DvsQc880zsihqcAlikLp5/HgYOhBkzvM/yjDNiV5SfBgyAKVN85rRhw+BPf8qr2dMUwCKbYsUKOPdc2G8/aN/eT7r94hexq8pv3br5qJLjjoO//AWGDoW33opdVYNQAIvURkUF/POfsO223td76qnw2mvwgx/ErqwwtGrlS9rfdx+8846PMjn/fJ9vI4cpgEW+z6pVfkHF9tvDscf62N4pU3y0Q8uWsasrPMccA++95/8Wf/+7d01cfDF88knsyupEASyyrvJyePFFOPNM6NrVJ9Vp0sSvziothV12iV1hYevYEe680ydzP/BAn2d5663963vvhWXLYldYaxYKcA5OWV9JSUkoLS2NXUYcK1b4MKfXXvO5CCZN8lnMWrTwy4pPO82HQumy4uw0Zw7cdhvcc4/PptaokU8Buuee8MMf+kxrW26Zln8/MysLIZTU+fkKYIE8C+CKCvj6a+8+WLHCt6VL/RLhRYt8+NicOf5n6/vvw9y5Vc/t1csH/++/v2/qZsgdIfgv0QkTYOJEePXVqpnV2rb1scW9ekGPHtCli1/w0aEDtGvnj7du7f/eLVr4DG21CGwFsDSIkhYtQmmfPpl7w3X/36XuV98fwvpbZaVvFRVVW3m5b2vX+sD9moYpmfkPYI8ePgtX//4+ZWRJCWyxRcN+nxLPqlUwfTqUlfklze+9579058zx/yffxwyaNvWup8aNPZAbNaratt4aXnih3gGsyXjENW0KmQxgWL+Fkbpffb/Zd7dGjapuq2+NG/vWrJl/L82b+9a6tW/t2vmwsc0393kGivVfP++1aAG77upbdSF4P/HChfDFF/6X0Zdf+l9Kq1b5X0/ffOMhvWaN/2IvL6/6hV9Z2WATLel/objeveHRR2NXIZJ+Zv4LuV272JVoFISISCwKYBGRSHQSTgAwsxXAe7HrqEFHINsvfcqFGiE36syFGvuFEFrX9cnqA5aU9+pzNjcTzKxUNTaMXKgzV2qsz/PVBSEiEokCWEQkEgWwpIyNXUAtqMaGkwt15n2NOgknIhKJWsAiIpEogEVEIlEAFzgz29/M3jOzD8zs4tj1pJhZdzN70czeNrO3zOycZH97M3vOzGYlt9GvJzWzRmb2hpk9mdzvZWZTk890vJk1iVzfZmb2sJm9a2bvmNmu2fY5mtnI5N95ppk9aGbNsuFzNLM7zWyRmc2stm+Dn525G5N63zSzwTW9vgK4gJlZI+BmYBiwHXCUmW0Xt6pvlQPnhRC2A3YBzkxquxh4IYTQF3ghuR/bOcA71e5fDYwKIfQBlgInR6mqyg3A0yGE/sBAvNas+RzNrCvwW6AkhDAAaAQcSXZ8jncD+6+zb2Of3TCgb7INB0bX+OohBG0FugG7As9Uu/874Hex69pIrY8D++JX63VJ9nXBLyCJWVe35IdwL+BJwPCrt4o39BlHqK8tMJvkhHu1/VnzOQJdgU+B9vjFYU8CP82WzxHoCcys6bMDxgBHbei4jW1qARe21H/8lLnJvqxiZj2BQcBUoHMIYX7y0AKgc6SyUq4HLgQqk/sdgGUhhNSkxLE/017AYuCupJvkdjNrSRZ9jiGEecC1wBxgPrAcKCO7PsfqNvbZbfLPkwJYspqZtQIeAUaEEL6s/ljwZka0cZRm9jNgUQihLFYNtVAMDAZGhxAGAStZp7shCz7HdsBB+C+LLYGWrP9nf1aq72enAC5s84Du1e53S/ZlBTNrjIfv/SGE1GTFC82sS/J4F2BRrPqA3YFfmNnHwDi8G+IGYDMzS82zEvsznQvMDSFMTe4/jAdyNn2O+wCzQwiLQwhrgUfxzzabPsfqNvbZbfLPkwK4sL0G9E3ONjfBT3xMiFwT4GeUgTuAd0IIf6/20ATg+OTr4/G+4ShCCL8LIXQLIfTEP7uJIYRfAy8ChyWHxa5xAfCpmfVLdu0NvE0WfY5418MuZtYi+XdP1Zg1n+M6NvbZTQCOS0ZD7AIsr9ZVsWGxOt61ZccGHAC8D3wI/CF2PdXq+hH+p92bwLRkOwDvY30BmAU8D7SPXWtS71DgyeTrrYFXgQ+Ah4CmkWvbEShNPsvHgHbZ9jkClwHvAjOB+4Cm2fA5Ag/i/dJr8b8mTt7YZ4efgL05+VmagY/q+N7X16XIIiKR1KsLIlsH8YuI5II6t4CTQfzv42Mz5+L9iUeFEN5uuPJERPJXfVbE2Bn4IITwEYCZjcOHkmw0gDt27Bh69uxZj7eU+lq8GJYuhW22+e7+dq/3ilOQSA57rvIhq8/z6xPAGxp0/MN1DzKz4fhlefTo0YPS0nqt4CH1dMEFcPPNsO4/w75Fh8cpSKSApX0YWghhbAihJIRQ0qlTp3S/ndSgvBwaN45dhYhA/QI4qwfxy4atXg1Nos7NJSIp9QngrB3ELxu3ahW0aBG7ChGBevQBhxDKzews4Bl8+rg7QwhvNVhlkhZffQWtWsWuQkSgfifhCCE8BTzVQLVIBnzxBXToELsKEQHNBVFwFiyAzTePXUUkZr6JZAkFcAGprITZs6GXhvyKZIV6dUFIbvnoIx8Fse5FGAUjddVnqhVsSfsjVH73cZEMUQu4gPy//+e3u+0Wtw4RcWoBF5AJE7z/d9ttY1cSWaqlGyr8tqgRANbEfxzC2mQVnMqKTFcmBUYt4AIxf74H8PHHQ5H+1UWyglrABeKyy7zh95vfxK4kCyUt3bDab63YfyyseUvfv2at365dE6E4yWdqCxWA//wHxoyBESOgb9/Y1YhIigI4z731Fhx+OGy9NVx+eexqRKQ6dUHksZkzYa+9oLgYnnoKWraMXVFuCOXl37ktatbMb9tWXUIYVn0NQOWqVRmuTvKJWsB5qLISbrkFdtnFw3fSJOjXr8aniUiGqQWcZ2bPhlNOgYkT4ac/hdtvh27dYleV2yq/+ca/SN0Cjdq0AaB4K5+RNXz5FQAVS5dmtjjJaWoB54lZs3yEQ79+8NprcNtt8H//p/AVyWZqAee4adPgqqvgoYd8pYvf/AYuvhi6d6/5uVJ3FV9+6V8kt42S1V6Kdujvt4uXAVA+f0Hmi5OcoQDOQXPnwrhx8MAD8MYb0Lq1r/U2YgRssUXs6kSkthTAOWLJEnj4YQ/dl17yiyp22glGjYITToDNNotdYWGrWLzYv0hui3r28P17DgagyVxvEVfM+ijzxUnWUgBnqfJyX7n4mWfg2Wdh6lSoqPA+3ksvhaOPhj59YlcpIvWhAM4in3ziYfvss/D887Bsmc+cWFLi/bqHHgqDBmlO8VxQ/vEcABolt2zv4wBXHLkLAG1mrQAglGkVr0KmAI6kvBxmzIDJk32bMsWHkIGPXDj0UB9GtvfeWkJIJF8pgDNk6VJ45ZWqwJ06FVau9Me6dIHdd4dzzoH99oP+/dXKzTcVb70HQOukwVu520AAFo70yZk7Tl8NQPHEsswXJ9EogNNg2TIfnVBWBq+/7rfvv++PNWoEAwfCiSf6xOi77QY9eihwRQqRArieliypCtnU7YcfVj3evTsMGQLHHedhu9NOWhZewCZPB2CLyX5/9QE7AfDJNbsC0GWKL5PU4tGpmS9OMkYBXEupBS2nT//u9vHHVcf07Olhe9JJfjt4MCTj80VE1qMA3oCVK/0EWfWgnTEDVviJa4qKfGHLH/4QTj/dg3bwYGjfPm7dkruaPvUaAL2f8vtfHu2jJWaP2+HbY9o+69PZtb9zSmaLk7Qp6AAOwa8qW7dVO2tW1bJhbdrADjt4F8LAgbDjjrD99tCiRdzaRST3FUwAr14Nb7+9ftguWVJ1TK9eHrJHH+23Awd6t4JOkEmmtXngleS2at+iM3zERMuXvF/rw//15U22GDU5s8VJg8nLAF61yiepKSvzq8lefx3efdfH3gI0bw4/+AH88pdVQbvDDt7aFRHJlJwP4K+/9pZsKmxLS72lW+knkenc2U+I/eIXVWHbp48PBxPJJZvf4i3dlbf4/YqLvQX8s7d8DuLR4w4EoPtf1CLOFTkXwJ995otMTprkFzPMnOlzJICPOCgpgUMO8dshQ2DLLdWFICLZKesDeP58D9vUlrqgoU0bX3LnZz/zoC0p8Ut4FbZSKLpe5S3dJ69qB0DFX/zM8dWzfezwr8aPAKDXxRo1ka2yLoBD8BUd7r3XJ6R5z6/gpE0b2GMPGD4chg710QjqRhCRXJY1Afz553DffXDnnd6t0Ly5r+h7yilVgVucNdWKZJ+el3hL96JLfghA5TW+/5G5PqJix6RF3Pu8VzJfnGxQjZFmZt2Be4HOQADGhhBuMLP2wHigJ/Ax8KsQwiavSDhnDpx/Pjz2GKxdCzvvDGPGwBFHQNu2m/pqIiK5ozZtynLgvBDC62bWGigzs+eAE4AXQghXmdnFwMXARZvy5p9/Dvvu6/28Z53ll/AOGLCp34KIbEjvC7xF/MsL/Kq6MMr3P/PZNH983GkA9DlXLeJYagzgEMJ8YH7y9QozewfoChwEDE0OuweYxCYE8OrVfgJtzhzv6919902sXEQkx21Sr6qZ9QQGAVOBzkk4AyzAuyg29JzhwHCAHj16fLt/wQIfs1tS4jOEiUh69RnpLd2fjtwRALvO93/bR/zgSKCq5SzpV1TbA82sFfAIMCKE8GX1x0IIAe8fXk8IYWwIoSSEUNKp2tRgW20FY8f6WN4TTvDuCBGRQlKrADazxnj43h9CeDTZvdDMuiSPdwEWbeqbn3QSXHklPPigXzDxq1/B009XXVghIpLPLIQNNlyrDjAzvI93SQhhRLX91wBfVDsJ1z6EcOH3vVZJSUkoLS1db//MmXDHHT4M7Ysv/IKKE06Agw/WeN9M2bfo8NglSGSzr/LJ4P91xPXf7jvin/4jnxriJt/1XOVD9br0qzYt4N2BY4G9zGxash0AXAXsa2azgH2S+3UyYACMGgXz5sFDD/n9K67w/uEOHeDnP4frrvNJddQ6FpF8UWMLuCFtrAW8IalLkF980W9nzfL9bdv6FXFDh/rIiYEDoVmzdFVcONQClg2Z8yefAvOMI/4NwJj7fMKf1GXQha6+LeCsvbasSxc46ijfwFvHqUl4Jk2CJ57w/cXF3mJOzQdRUuJTTTZtGqtyEZHaydoWcE3mzYNXX62agrK0tGpy9caNPYRTgTx4sK9ioZbyxqkFLLWxYKS3iHsf4n+SfvxgHwA6jS7MPuK8bQHXpGtXn3bykEP8fgjwySdVYVxWBv/6lw91Az+R169f1ZzAqeWFttgi3vcgIoUtZ1vAtRECfPSRn7yrvgzRp59WHbP55t8N5YEDoX9/b0UXErWApS6WnOQjJ5bvtxKADhN8scTUkkr5rmBbwLVhBr17+3Z4tXxZsgTefPO7ofyPf/jl0QBNmsB2260fzB06xPk+RCQ/5XULeFOsXeuTvVcP5WnTYOHCqmO6dq3qukiFct++vkx9rlMLWBrCyl/6VJgLdvEfim4v+kKMTZ96LVpN6aQWcANp3NhP1G2/va+KnLJw4forKT/7bNUCn61aeSAPGeIn+4YM8b5mzV0sIjVRC7gOUkvcT5sGb7zhJ/ymTfPVmMEnk99xx6pAHjIEtt02u/uV1QKWdCjfawgAnw/0caGdp3pfsU2eHq2mhqQWcARNm8KgQb6deKLvq6jw5ZPKyvykX1kZ3HMP3HyzP968uU82v9tuvu26q/qURQqdWsBpVFnpV/CVlfmY5cmTvcWc6r7o168qkHfbzUdfxOpPVgtYMsGGbA/Al31bA7DZDF9Ep+Kt96LVVB9qAWexoiIP2X79qvqVV63yccqTJ/s2YQLcdZc/1qED7LMP7Lefb926xatdRNJPAZxhLVr4XBZ77OH3Q/BW8pQpPu/Fs8/C+PH+2HbbVYXxHntAy5bx6hZpCKHsLQBalyU7+m4NQMWeg789pulHiwEo/+RT8l0eDKDKbWawzTZw/PFw991+ifWbb8K113oL+NZb4YADoH17OOggD+fUyT4RyW3qA85yX38NL78MTz3lU3V+9pkPfTv4YO/W2GefhhldoT5gyRbFXbcEoLJDGwBs/hcAVCxeHK2mjcnEfMASUfPm3gVx/fW+gOnEiXDkkfDkk94y7toVLr3UJ7IXkdyiAM4hjRrBnnvCbbf5oqaPPebD2S67zNfYO/dc78IQyWXl8z6jfN5nVL75LpVvvuvDhsrLKd6qO8VbdaeodWuKWreOXWaDUADnqKZNvU/48cdhxgyfFe7GG6FXLw9i9ROLZD8FcB4YMMDX05s1y0/mjRrlV+JN1qIFkgcqli6lYulSyj/51EdGVFZCZSWNOnagUccOFDVrRlGOTvatAM4jvXp598TEiT650I9+5P3DGTzPKiKbQOOA89Cee/pQtrPO8iO5VDkAAAbpSURBVP7hVavg6qt9yJtIrqtc6fNJkNxasv5YUTJQPiTzyobUJadZTAGcp1q39nHFLVvCNdf4HMd//WvsqkSkOgVwHjPzyYC++QauvBKGDfOVpEXyybct3uTWkrlgUy3jsGZNcmD29cWpDzjPmfnoiB494JRTIPV/UUTiUwAXgFat4IYb4N13feywSD4L5eW+rV79basYgKJGvmURBXCB+PnP/WKNMWNiVyIiKQrgAlFUBMcd5zOuLV8euxoRAQVwQfnJT/w8xJQpsSsRyaAQfKus8M1s/S0SBXAB2Xlnv3399bh1iIjTMLQC0ro1dO4Ms2fHrkQkoiwajqYWcIHZaiuf1lJE4lMAF5gOHWDJkthViAgogAtOu3awdGnsKkQENiGAzayRmb1hZk8m93uZ2VQz+8DMxptZk/SVKQ2lRQvNFSySLTalBXwO8E61+1cDo0IIfYClwMkNWZikR/Pmvs6ciMRXqwA2s27AgcDtyX0D9gIeTg65Bzg4HQVKw2rc2Fd4EZH4atsCvh64EKhM7ncAloUQUj/Kc4GuDVybpEFxsU/WLiLx1RjAZvYzYFEIoawub2Bmw82s1MxKF2fhstKFplEjX9FFROKrTQt4d+AXZvYxMA7vergB2MzMUhdydAM2uB5vCGFsCKEkhFDSqVOnBihZ6qOoCCoqYlchIlCLAA4h/C6E0C2E0BM4EpgYQvg18CJwWHLY8cDjaatSGoxZVl0IJFLQ6jMO+CLgXDP7AO8TvqNhSpJ0UgCLZI9NmgsihDAJmJR8/RGwc8OXJOmkhTlFsoeuhBMRiUQBXGDUAhbJHgpgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkkloFsJltZmYPm9m7ZvaOme1qZu3N7Dkzm5Xctkt3sSIi+aS2LeAbgKdDCP2BgcA7wMXACyGEvsALyX0REamlGgPYzNoCewB3AIQQ1oQQlgEHAfckh90DHJyuIkVE8lFtWsC9gMXAXWb2hpndbmYtgc4hhPnJMQuAzht6spkNN7NSMytdvHhxw1QtIpIHahPAxcBgYHQIYRCwknW6G0IIAQgbenIIYWwIoSSEUNKpU6f61isikjdqE8BzgbkhhKnJ/YfxQF5oZl0AkttF6SlRRCQ/1RjAIYQFwKdm1i/ZtTfwNjABOD7ZdzzweFoqFBHJU8W1PO5s4H4zawJ8BJyIh/e/zOxk4BPgV+kpUUQkP9UqgEMI04CSDTy0d8OWIyJSOHQlnIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkkloFsJmNNLO3zGymmT1oZs3MrJeZTTWzD8xsvJk1SXexIiL5pMYANrOuwG+BkhDCAKARcCRwNTAqhNAHWAqcnM5CRUTyTW27IIqB5mZWDLQA5gN7AQ8nj98DHNzw5YmI5K8aAziEMA+4FpiDB+9yoAxYFkIoTw6bC3Td0PPNbLiZlZpZ6eLFixumahGRPFCbLoh2wEFAL2BLoCWwf23fIIQwNoRQEkIo6dSpU50LFRHJN7XpgtgHmB1CWBxCWAs8CuwObJZ0SQB0A+alqUYRkbxUmwCeA+xiZi3MzIC9gbeBF4HDkmOOBx5PT4kiIvmpNn3AU/GTba8DM5LnjAUuAs41sw+ADsAdaaxTRCTvFNd8CIQQ/gz8eZ3dHwE7N3hFIiIFQlfCiYhEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRicRCCJl7M7PFwCcZe0PZFFuFEDrFLkKkkGQ0gEVEpIq6IEREIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQi+f9WLJycozuqKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "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": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Err         \n",
      "-------------------\n",
      "    0|2.861463e-01|\n",
      "   10|1.860154e-01|\n",
      "   20|8.144529e-02|\n",
      "   30|3.130143e-02|\n",
      "   40|1.178815e-02|\n",
      "   50|4.426078e-03|\n",
      "   60|1.661047e-03|\n",
      "   70|6.233110e-04|\n",
      "   80|2.338932e-04|\n",
      "   90|8.776627e-05|\n",
      "  100|3.293340e-05|\n",
      "  110|1.235791e-05|\n",
      "  120|4.637176e-06|\n",
      "  130|1.740051e-06|\n",
      "  140|6.529356e-07|\n",
      "  150|2.450071e-07|\n",
      "  160|9.193632e-08|\n",
      "  170|3.449812e-08|\n",
      "  180|1.294505e-08|\n",
      "  190|4.857493e-09|\n",
      "It.  |Err         \n",
      "-------------------\n",
      "  200|1.822723e-09|\n",
      "  210|6.839572e-10|\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZyVZfnH8c81C8sAsousAi6gogiO5lJGaiZaqWW5lHthVpZLqW2/1MpfpqWW5oIbmglp/hTN3CVLCQU3UAQUFEH2fWeW+/fH9ZyZwzjDDLPdZ/m+X6/zGs72nGsOM9+5zv3cz/1YCAEREWl9BbELEBHJVwpgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwSmZl9xsxmtcB2v2FmTzfwsWeZ2X929D5pGgWwZL0kIKab2UYzW2xmt5hZl+S+W81sfXLZamZladf/2Qq1BTPbfXuPCSH8O4QwpJHb/7SZvWxma8xspZm9ZGYHJtu9P4RwdGO2K61DASxZzcwuAa4Bfgx0Bg4GdgWeMbM2IYTvhBA6hhA6AlcDE1LXQwij41XuzKyoCc/dCXgc+BPQDegLXAlsaZ7qml9Tvt9cpACWrJUE0JXABSGEJ0MIZSGED4CvAwOBbzZim6PMbIGZXWpmS81skZmdYGbHmtnspMv8adrjDzKzyWa2OnnsTWbWJrnvxeRhbyYd98lp27/MzBYDd6duS56zW/IaI5PrfcxsmZmNqqXcPQFCCA+EECpCCJtCCE+HEN5KnrvN0EHSjX/HzOYk9d5sZlbH+3Ctmf3HzDqn3Xadma0ys3lmNjrt9j5mNjGp+z0z+3bafVeY2UNm9hczWwucldz2NzO718zWmdnbZla6g/9VOUEBLNnsUKAd8HD6jSGE9cATwOcbud1dku32Bf4HGIuH+QHAZ4BfmNmg5LEVwEVAD+AQ4Ejgu0kdhyePGZ503BPStt8N79TH1Kj9feAy4C9mVgLcDYwLIUyqpc7ZQIWZjTOz0WbWtQHf2xeBA4H98D9UX0i/08wKzGxscv/RIYQ1yV2fAmYl3+fvgDvTwns8sADoA5wEXG1mR6Rt9njgIaALcH9y25eT53UBJgI3NaD2nKMAlmzWA1geQiiv5b5Fyf2NUQb8JoRQhodED+DGEMK6EMLbwDvAcIAQwrQQwn9DCOVJ930b8Nl6tl8J/DKEsCWEsKnmnSGEscB7wBSgN/Cz2jYSQlgLfBoI+B+JZUkn2ms7r/3bEMLqEMJ84AVg/7T7ioEH8D8OXwohbEy778MQwtgQQgUwLqmrl5n1Bw4DLgshbA4hvAHcAZyR9tzJIYRHQgiVad/vf0IITyTbu4/k/cw3CmDJZsuBHnWMK/ZO7m+MFUkwAKQCY0na/ZuAjgBmtqeZPZ7s/FuLjzPXF/zLQgib63nMWGAY8KcQQp1juiGEmSGEs0II/ZLH9wFu2M52F6f9e2Pq+0jsjnerV4YQttb1vLRg7pi83soQwrq0x36If3pI+agBdbTLx/FhBbBks8n4DqevpN9oZh2B0cBzrVDDLcC7wB4hhJ2AnwK1jqum2e4ShEn9NwB3AleYWbeGFBJCeBe4Bw/ixpgJnA3808waOivjY6CbmXVKu20AsDC9tEbWk/MUwJK1kvHJK4E/mdkxZlZsZgOBv+Fjkve1QhmdgLXAejMbCpxf4/4lwOAd3OaNwNQQwreAfwC31vYgMxtqZpeYWb/ken/gVOC/O/h6VUIID+B/RJ41s90a8PiPgJeB/zWzdma2H3Au8JfG1pBPFMCS1UIIv8MD4zo8CKfgH3mP3N5H92b0I+A0YB0+bDChxv1XAOOSWQdfr29jZnY8cAzVQX4xMNLMvlHLw9fhO8emmNkGPHhnAJc04vuoEkIYB1wFPJ/8QavPqfisk4+B/8PHt59tSg35wrQgu4hIHOqARUQiUQCLiESiABYRiUQBLCISSd5NfJba9ejRIwwcODB2GSJZZdq0actDCD0b+3wFsAAwcOBApk6dGrsMkaxiZh825fkaghARiUQBLJJvNmwAzf/PCApgkXwwZw6ccQbstht07Ajdu8PRR8PDDyuMI1IAi+Sy8nK44goYNgweeQRGjICrroKTToJ58+CrX4XjjoP582NXmpe0E04kV5WXw+mnw/jxcOqp8PvfQ+/e295/003wi1/AZz4DkybBoEF1bk6anzpgkVxUUeFDDuPHw29/C3/967bhC1BUBBdeCC++COvWwahR3hVLq1EAi+SiK6+EBx7w8L3ssu0/dsQIeO45D+Hjj4dNnzhJh7QQBbBIrnn+efj1r+Hss+sP35QRI7xLnj4dLr64ZeuTKgpgkVyybBl885uw557wpz/t2HOPOQZ+/GO49VZ46KGWqU+2oQAWySWXXgrLl8OECdChw44//9e/htJS+P73Yc2a+h8vTaIAFskVkyfDPffARRfB8EaeZLhNG7jlFli61MeRpUUpgEVyQUWFd619+vi0sqYoLYUxY+CPf4QZM5qnPqmVAlgkF9xzD7z2ms/17dix3ofX6ze/gc6dvZuWFqMAFsl2W7b40W0HHQQnn9w82+zeHX72M3j2WfjXv5pnm/IJCmCRbHfnnX4o8a9+BWbNt93zz/eDN37xC60X0UIUwCLZbNMmHy749Kfh859v3m23bw8//Sn8+9/eCUuzUwCLZLM77oCPP27+7jfl29+G/v3hl79UF9wCFMAi2aqsDK67zrvfUaNa5jXatvWj6SZPhpdeapnXyGMKYJFs9eCDPvZ76aUt+zpnn+075X73u5Z9nTykABbJRiF4IO61l6/n25JKSnyO8WOPwcyZLftaeUYBLJKNnnkG3nzT124oaIVf4+99z3fKXXtty79WHlEAi2Sj66+HXXaB005rndfr2dOHIu6/3w9TlmahABbJNrNnw5NP+jzdtm1b73UvuAC2boXbb2+918xxCmCRbHPTTVBc7Os1tKahQ/1Enrfc4jMwpMkUwCLZZN06X/fh5JN9CKK1/eAHPu/44Ydb/7VzkAJYJJuMG+chfMEFcV5/9Gg/tf0f/xjn9XOMAlgkW4QAf/4zHHigL7wTQ0EBfPe78PLLPgtDmkQBLJItXnzR5+Gef37cOs46C9q181MXSZMogEWyxa23QpcuzbfkZGN16+Y1/OUvPhwijaYAFskGS5bA3//u3WdJSexqvAtfv97nBUujKYBFssFdd/nUr/POi12JO+gg2H9/n5KmVdIaTQEskukqK2HsWF/xbOjQ2NU4M/jOd+Ctt2DKlNjVZC0FsEime/ZZmDcvc7rflNNOgw4ddGRcEyiARTLd7bdDjx5w4omxK9lWp04ewuPHw+rVsavJSgpgkUy2eDE8+qjvfGvNdR8a6rzz/LRI2hnXKApgkUx2991QXu6nBspEBxwAI0fCbbdpZ1wjKIBFMlX6zrc994xdTd3GjIHp07UzrhEUwCKZKlN3vtWknXGNpgAWyVS33ZaZO99qSt8Zt2ZN7GqyigJYJBMtWgQTJ2buzreaxozRzrhGUACLZKJM3/lWU2mpdsY1ggJYJNNUVPjOt899LrN3vtU0ZowfGfff/8auJGsogEUyzdNPwwcf+KG+2eS006BjR++CpUEUwCKZ5tZboVcvOOGE2JXsmE6d4PTTYcIEWLkydjVZQQEskknmz4fHH4dzz4U2bWJXs+POOw82b/ZTJ0m9FMAimeSOO3wnVrbsfKtp+HA45BDv4rUzrl4KYJFMsXWr73wbPRoGDoxdTeOdfz7Mng3PPRe7koynABbJFA8/7IvvfP/7sStpmq99zQ8guemm2JVkPAWwSKa46SbYfXf4whdiV9I07dr5lLTHHvPZHFInBbBIJnj9dXjpJfje9/zU79kuNYVOZ07erhz4nxbJATfd5CfbPOus2JU0j/79fRrdHXf4IcpSKwWwSGzLlvkaCqef7qedzxUXXAArVmh9iO1QAIvEdsstsGULXHhh7Eqa12c/62dOvv56TUmrgwJYJKbNm+Hmm+G44zLnjMfNxQwuvhjeeQeeeip2NRlJASwS0/33w9KlHlS56OSToU8f+MMfYleSkRTAIrFUVnow7b+/r3yWi9q08bHgZ56BN9+MXU3GUQCLxPL44/7x/JJL/ON6rjrvPF8l7ZprYleScRTAIjGEAL/5DQwaBKecErualtW1qx+ePGECvPde7GoyigJYJIbnn4dXXoHLLoOiotjVtLyLLoLiYvjd72JXklEUwCIxXH019O4NZ54Zu5LW0bs3nHMO3HMPLFgQu5qMoQAWaW0vvugd8I9+5Osm5ItLL/Whl//939iVZAwFsEhrCgF+/nPvCM8/P3Y1rWvgQF9ofuxYLdKTUACLtKann4Z//9tDuH372NW0vp//3BcbuvLK2JVkBAWwSGtJdb+77grf+lbsauLo1w+++1249154993Y1USnABZpLePHw9SpcMUV2Xm+t+byk59Ahw4+AyTPKYBFWsPGjR44I0bAGWfEriaunj3hZz+DiRPh2WdjVxOVAlikNVx3HXz0EdxwQ24suN5UP/yhH4Ry0UVQXh67mmj0kyDS0ubP98NwTzoJDj88djWZoV07uPZamDEDbrstdjXRKIBFWlII1dPNrr02bi2Z5itfgaOO8jHhhQtjVxOFAlikJU2YAE884es+ZPOp5luCmZ8zrrzcZ0bk4aLtCmCRlrJ8OfzgB3Dggb4ko3zSbrv5nOCJE+HBB2NX0+oUwCItIQRf+2DNGj8xZWFh7Ioy10UXQWmpn0n5o49iV9OqFMAiLeHmm+Gxx3z1r/32i11NZisqgr/+FbZuhW9+EyoqYlfUahTAIs3ttdd8oZ1jj/UhCKnfHnv4H60XX4SrropdTatRAIs0p0WL4PjjYeed4e67c/tMF83tjDN8ec6rroKHHopdTavIg5WgRVrJ5s1w4omwciW89JKHsDRcalbE7NkexoMHw8iRsatqUeqARZrD1q1+oMWUKXDffX6iTdlx7drB//2fH648enTOL9ijABZpqvJyOO00+Mc/4JZb/AADabxevXzZTjM48kh4//3YFbUYBbBIU2zY4IH797/D9df7VCppuiFDfKGeLVvgM5+BN96IXVGLUACLNNaSJXDEEd753nwzXHhh7Ipyy7BhMGmSz6E+/HB46qnYFTU7BbBIYzz7LAwfDtOn+5jld78bu6LcNGwYTJ7sK6eNHg3/8z85tXqaAlhkR6xbBxdfDEcfDd26+U63L385dlW5rV8/n1Vyxhnwq1/BqFHw9tuxq2oWCmCRhqiogL/8Bfbay8d6zzsPXn0V9t03dmX5oWNHP6X9fffBzJk+y+RHP/L1NrKYAlhkezZu9AMq9tkHTj/d5/ZOnuyzHTp0iF1d/vnmN2HWLP+/+MMffGji8svhww9jV9YoCmCRmsrL4YUX4Hvfg759fVGdNm386KypU+Hgg2NXmN969IC77vLF3I87ztdZHjzY/33vvbB6dewKG8xCHq7BKZ9UWloapk6dGruMONat82lOr77qaxFMmuSrmJWU+GHF3/mOT4XSYcWZaf58GDsWxo3z1dQKC30J0M99Dj71KV9prU+fFvn/M7NpIYTSRj9fASyQYwFcUQGbNvnwwbp1flm1yg8RXrrUp4/Nn+8fW2fPhgULqp87aJBP/j/mGL9omCF7hOB/RCdOhOefh1deqV5ZrXNnn1s8aBAMGAC9e/sBH927Q9eufn+nTv7/XVLiK7Q1ILAVwNIsSktKwtTdd2+9F6z5c5e6nn57CJ+8VFb6paKi+lJe7peyMp+4X980JTP/BRwwwFfhGjrUl4wsLYVddmne71Pi2bgR3nwTpk3zQ5pnzfI/uvPn+8/J9phB27Y+9FRc7IFcWFh9GTwYnnuuyQGsxXjEtW0LrRnA8MkOI3U9/XazbS+FhdVf0y/FxX5p186/l/bt/dKpk1+6dvVpYzvv7OsMFOlHP+eVlMAhh/glXQg+TrxkCaxY4Z+M1q71T0obN/qnp82bPaS3bvU/7OXl1X/wKyubbaEl/RSK2203ePjh2FWItDwz/4PctWvsSjQLQkQkFgWwiEgk2gknAJjZOmBW7Drq0QPI9EOfsqFGyI46s6HGISGETo19ssaAJWVWU/bmtgYzm6oam0c21JktNTbl+RqCEBGJRAEsIhKJAlhSbo9dQAOoxuaTDXXmfI3aCSciEok6YBGRSBTAIiKRKIDznJkdY2azzOw9M7s8dj0pZtbfzF4ws3fM7G0z+2Fyezcze8bM5iRfox9PamaFZva6mT2eXB9kZlOS93SCmbWJXF8XM3vIzN41s5lmdkimvY9mdlHy/zzDzB4ws3aZ8D6a2V1mttTMZqTdVut7Z+6PSb1vmdnI+ravAM5jZlYI3AyMBvYGTjWzveNWVaUcuCSEsDdwMPC9pLbLgedCCHsAzyXXY/shMDPt+jXA9SGE3YFVwLlRqqp2I/BkCGEoMByvNWPeRzPrC/wAKA0hDAMKgVPIjPfxHuCYGrfV9d6NBvZILmOAW+rdeghBlzy9AIcAT6Vd/wnwk9h11VHro8Dn8aP1eie39cYPIIlZV7/kl/AI4HHA8KO3imp7jyPU1xmYR7LDPe32jHkfgb7AR0A3/OCwx4EvZMr7CAwEZtT33gG3AafW9ri6LuqA81vqBz9lQXJbRjGzgcAIYArQK4SwKLlrMdArUlkpNwCXApXJ9e7A6hBCalHi2O/pIGAZcHcyTHKHmXUgg97HEMJC4DpgPrAIWANMI7Pex3R1vXc7/PukAJaMZmYdgb8DF4YQ1qbfF7zNiDaP0sy+CCwNIUyLVUMDFAEjgVtCCCOADdQYbsiA97ErcDz+x6IP0IFPfuzPSE197xTA+W0h0D/ter/ktoxgZsV4+N4fQkgtVrzEzHon9/cGlsaqDzgM+LKZfQCMx4chbgS6mFlqnZXY7+kCYEEIYUpy/SE8kDPpfTwKmBdCWBZCKAMext/bTHof09X13u3w75MCOL+9CuyR7G1ug+/4mBi5JsD3KAN3AjNDCH9Iu2sicGby7zPxseEoQgg/CSH0CyEMxN+750MI3wBeAE5KHha7xsXAR2Y2JLnpSOAdMuh9xIceDjazkuT/PVVjxryPNdT13k0EzkhmQxwMrEkbqqhdrIF3XTLjAhwLzAbeB34Wu560uj6Nf7R7C3gjuRyLj7E+B8wBngW6xa41qXcU8Hjy78HAK8B7wINA28i17Q9MTd7LR4CumfY+AlcC7wIzgPuAtpnwPgIP4OPSZfiniXPreu/wHbA3J79L0/FZHdvdvg5FFhGJpElDEJk6iV9EJBs0ugNOJvHPxudmLsDHE08NIbzTfOWJiOSuppwR4yDgvRDCXAAzG49PJakzgHv06BEGDhzYhJeUplq2DFatgj333Pb2rq8NilOQSBZ7pvJBa8rzmxLAtU06/lTNB5nZGPywPAYMGMDUqU06g4c00Y9/DDffDDX/Gz5f8LU4BYnksRafhhZCuD2EUBpCKO3Zs2dLv5zUo7wciotjVyEi0LQAzuhJ/FK7LVugTdS1uUQkpSkBnLGT+KVuGzdCSUnsKkQEmjAGHEIoN7PvA0/hy8fdFUJ4u9kqkxaxfj107Bi7ChGBpu2EI4TwBPBEM9UirWDFCujePXYVIgJaCyLvLF4MO+8cu4pWZrZjF5FWogDOI5WVMG8eDNKUX5GM0KQhCMkuc+f6LIiaB2HkHaun70g1waFy29u1boo0M3XAeeQ///Gvhx4atw4RceqA88jEiT7+u9desStpZanONTW+m+psk07YCmyb66SuV9boeNM64pC6L3WbumNpBHXAeWLRIg/gM8+EAv2vi2QEdcB54sorvUn79rdjVxJRzU64Dpa6v03dvx6WbCtUVPgNlTWuqzOWBlAvlAf+9S+47Ta48ELYY4/Y1YhIigI4x739NnztazB4MFx1VexqRCSdhiBy2IwZcMQRUFQETzwBHTrErihDpIYFQkVy1YccrJBtrycPt+Lk16SwsHobyTBFaiiCZOghVCRDD2Vl21zX0ITURh1wDqqshD//GQ4+2MN30iQYMqTep4lIK1MHnGPmzYNvfQuefx6+8AW44w7o1y92VRmuxg61quMwUvcn00YsrQO21KLKRcmvUDJ1zVLT08rLfRupTnjLVr89dT11f1VnrI44H6kDzhFz5vgMhyFD4NVXYexY+Oc/Fb4imUwdcJZ74w347W/hwQf9TBff/jZcfjn071//c6WGOqaW1arAu2FLOuDQNlnlPhkvDoVJR1zh27Ct3vnaZu+Ew6ZNydfN/nVrasxYY8X5RAGchRYsgPHj4a9/hddfh06d/FxvF14Iu+wSuzoRaSgFcJZYuRIeeshD98UXvTE68EC4/no46yzo0iV2hTmk5iyJrUk3mt4RV6Y6VP+aOpw5JB1wZYl3xJVtvVMOqTHict9G4SbveAvWeQdcsCHpiDds9K+pDjkZK/ZyNF6caxTAGaq83M9c/NRT8PTTMGWKz3QaMgSuuAJOOw123z12lSLSFArgDPLhhx62Tz8Nzz4Lq1f7dNPSUh/X/cpXYMQIrRne6lJjw+Vl1Tcl3WhBqitN5vumZkEUJP9JlcXtASjr5L9q5e2TTrmgLQCFZT45u3i9b6fNqi3+/NUbfHtr11e/5nq/rWq8OFWPOuKspQCOpLwcpk+Hl1/2y+TJPoUMfObCV77i08iOPFKnEBLJVQrgVrJqFfz3v9WBO2UKbPCGht694bDD4Ic/hKOPhqFD1eVmpPROMxkfrtySdL6pDnirz3IoSOb7FpclY7jBz4Ra2cY73007+X9weQefCRqSWRWFm31+cbvV3hm3W9616iXbLvPx4cIVa/056707rkzNpEi9VmVFE75JaU0K4BawerXPTpg2DV57zb/Onu33FRbC8OFw9tm+MPqhh8KAAQpckXykAG6ilSurQzb19f33q+/v3x8OOADOOMPD9sADdVr4nJIaHy5L5vem1oRIZi8UJEfAtUm+Fm7pDIBV+Njw+mLvfDf39L/AG3v7Ztclf5GLNrSteql2y/3fHRbvBEDJx975Fi9Z46+5JumMk5kUlVtTY8SaU5ypFMANlDqh5Ztvbnv54IPqxwwc6GF7zjn+deRI6NkzVsUikukUwLXYsMF3kKUH7fTpsG6d319Q4Ce2/NSn4PzzPWhHjoRu3eLWLRmgctu5w5Wp9SW2JLMbkvHajht84nbRhk4ArNvq84bX9/fOd3PvZDy339aqTa9PFlRavdYf2+5jHyfuNL8k+ep7a9su8k64cPlqr2VDavbEtl26OuL48jqAQ/Cjymp2tXPmVP9s7rQT7LefDyEMHw777w/77AMlJXFrF5HslzcBvGULvPPOJ8N25crqxwwa5CF72mn+dfhwH1bQDjLZYVVzh7dd9cySmQqpTrj9Wu9O26z22Q5t13hXu2aj/2qu3616k936ekfbv+8iAMqH+gyKD1b6R68VH3o33Wmuj3t1metddvuP/KNb4dJVAFSu89kTIenK1RHHk5MBvHGjL1IzbZofTfbaa/Duu1UrBNK+Pey7L3z1q9VBu99+3u2KiLSWrA/gTZu8k02F7dSp3ummDtXv1ct3iH35y9Vhu/vu257cQKTF1ZgtUVFzbDiZudB5pc+SaL/UO+I1y9pVbWL1EB/jfW+IzxU+rN9cAI7qOROAjbv7LImXRwwG4O25fQHoMMs75K5zfNsd5nkHXNURr05mUdTsiNPqlpaRdQH88cd+kslJk/xghhkzqs4GQ8+eftjuiSf61wMOgD59NIQgIpkp4wN40SIP29QldUDDTjv5KXe++EUP2tJSP4RXYStZITVbYot/ragaG/ZV0IrX+rhtj+XVU2s6LPIOdtXH/vXpvfcFYOFePtb79V1eBeDKXacDsKKfjyc/vt/+ADw5Z28Ait7253ed5RPSO831brtwcdIRJ68NaauyaZy4RWRcAIfgZ3S4915fkGbWLL99p53g8MNhzBgYNcpnI2gYQUSymYVW/ItWWloapk6dWut9y5fDfffBXXf5sEL79n5G31GjqgO3KOP+XOSOzxd8LXYJAtUf4cxnOBS0qz4Szjp5xxp6eVe8YZDvNV451H8xNu3r3epxQ2cAcGb3lwAYUuw7RD5OutgHVh8IwIQ5I317b/l2ur3rj+v0fnUHXNUV1zVOnOcd8TOVDzbpM3e9kWZm/YF7gV74eQpvDyHcaGbdgAnAQOAD4OshhFU7WsD8+fCjH8Ejj/j5Cg86CG67DU4+GTp33tGtiYhkj3o7YDPrDfQOIbxmZp2AacAJwFnAyhDCb83scqBrCOGy7W2rZge8fLmvArZokZ/J95xzYNiwJn5H0ijqgDNU2k6N1FmZra13xQU7+bzfil18dsT63bxDXjk0OQvHcO9kT9lzGgBndHkFgEHF/rh5ZT4b4t7VBwEwfvYBvv03O1W9Zrd3vdPt+H4yc2LxCqB6nDjfO+IW74BDCIuARcm/15nZTKAvcDwwKnnYOGASsN0ATrdli+9Amz/fx3oPO2wHKxcRyXI7NAZsZgOBF4FhwPwQQpfkdgNWpa7XeM4YYAzAgAEDDvjwww8BP/vDbrv57IUXX4Q2bZr4nUiTqAPOIklX/ImOuIuP2VXs4rMa1u3mnezK5Ig528/XiDh5j9cAOLWLz5rok2xnVpk/btyK6m7oH+/6R9L20331tm7v+myNDvN8W7bEDyUNydF1lZu9I86XFdia2gEXNPSBZtYR+DtwYQhhbfp9wVO81nc6hHB7CKE0hFDaM21psF13hdtv97m8Z53lwxEiIvmkQQFsZsV4+N4fQng4uXlJMj6cGideuqMvfs45cPXV8MADfsDE178OTz5ZfWCFiEgua8hOOMPHeFeGEC5Mu/1aYEXaTrhuIYRLt7etuqahzZgBd97p09BWrPADKs46C044QfN9W4uGILJYzSGJ9j5cUL2TLhmSGOw731YN8ceV7+MLAR2zxzsAfLHLGwB0L9xQtenpW/oB8NBi30H39sz+AHR+x3cfdZ3th1a3n+/T1FiWDEkkJxCtWhQ+R0+T1BpDEIcBpwNHmNkbyeVY4LfA581sDnBUcr1Rhg2D66+HhQvhwQf9+m9+4+PD3bvDl74Ev/+9L6qj7lhEckXGHIhRU8asHQIAAA0fSURBVOoQ5Bde8K9z5vjtnTv7EXGjRvnMieHDoV277WxIGkQdcA6pbyfdzt4RbxiUdMR7+OM2DPFudp/BC6s2dWg3X/CnpNB3rs3c4OdMemmBL/izdZYfxNElOWK18/u+zGabhX5IQFjlnXHlRl9sKNdOHNri09Bi6d0bTj3VL+DdcWoRnkmT4LHH/PaiIu+YU+tBlJb6UpNt29a1ZRGRzJCxHXB9Fi6EV16pXoJy6tTqxdWLiz2EU4E8cqSfxUKdct3UAeewujri5NDmyqQj3tTfx4xXDy6ueuq6wT6drP2ufuDFwG7+S1ZkfvtHa72rXrnQZ6B2fN97us5zkwM4PvSx4NQhzSE5gKMyWZA+lKdOHJqd09VytgOuT9++vuzkiSf69RB8bnEqjKdNg7/9zae6ge/IGzKkek3g1OmFdtkl3vcgIvktawO4JjM/fdDAgXDSSX5bCDB3ru+8S52C6D//8WlvKTvvvG0oDx8OQ4d6Fy2SE+o4PVLqMGJLutKSRb58Zcnc6uOptvT2Md51A/zr+wO8493cx7dVtNPWbb6mTqFUXuLRsqWzd9U7dfauu23yGgUrkhOGpk6PtLX65KP5dFhzzgRwbcz8aLvddoOvpX3CXrkS3npr23PD/elPfng0+FF5e+/9yWDu3j3O9yEiuSmnA7gu3bpVL3OZUlbmi72nh/KTT8K4cdWP6du3eugiFcp77OGnqRfJGnV0xCSHEResq16Osu1SP/13uw+S2Q69vAPe2Md3qGzYxe/f3CN5QgffdnmJf93Q14dIK5K1Bjp29O66pIN3xIXJKZfCmurXzNUZE7XJywCuTXGx76jbZx8/K3LKkiWfPJPy009Xn+CzY0cP5AMO8J19BxzgY81au1hE6qOYqEevXnD00X5JSZ3i/o034PXXfYff2LF+NmbwxeT33786kA84APbaS+PKkqFSY63BO83KzWkdZ3IkmyVHthWt9NkMnT9KTmfU3TvjLT1TnbD/kG/uknS+ycyjyjZ+fWOPZKnMwqSzbuvXi9tWr8ZVuMb/XVk1Puw1ZPuMidoogBuhbVsYMcIvZ5/tt1VU+OmTpk3znX7Tpvnwxc03+/3t2/ti84ce6pdDDtGYski+UwA3k8JC33G3995w+ul+W2WlH8E3bZrPWX75Zbj22urhiyFDqgP50EN99oXGkyWj1Dx5aDJbwTb46Y9slc9maLfYZze02ymZW9zFr2/t6mO9ZR29060o9k64MlnfZWuXpPO16kXgi4o9lgqSj4whGZMOm/2Xo2p9iRxY8lIB3IIKCjxkhwypHlfeuNHnKb/8sl8mToS77/b7uneHo46qHvLo1y9e7SLS8hTAraykxNeyOPxwvx6Cd8mTJ/u6F08/DRMm+H17710dxocfDh06xKtbBKieQVHmnXBqXNZSc4rX+FLhBct9Rbb2HZKx3g5+vbKTDwpXtE+626LkQLK0T36VHb0rtuA/8FaQPKbId7IUpI6iy4GxYX3gjcwM9twTzjwT7rnHD7F+6y247jrvgG+9FY491qfOHX+8h3NqZ5+IZDd1wBnGzNex2HdfuOQS2LQJ/v1veOIJX6pz4kSf+nbCCT6scdRRml0hEdU1p7jG7Alr4z+kBe29Ay5M1iwO7ZIx4DbVP8ShcNvlFUIyQ8Iqa3S4qR0mW/3xVfOGs2hsWB1whmvf3ocgbrjBT2D6/PNwyinw+OPeGfftC1dc4QvZi0h2ydrV0PLdli1+pN5dd3lX3KEDjBnjXXPfvju+Pa2GJi2q5opsqSOVko9vljYP2FIf6VKPSY0BpzrgVLedzMioGgtOrSeRWuuiFdaUaLWTckpmadvWx4QffRSmT/dV4f74Rxg0CC6+WOPEItlAAZwDhg3z8+nNmeM7866/3o/Ee/nl2JWJJEKAEAjl5YTyciq3bPHLho1+WbOu+rJqtV9Wr6Fy9RrC2vV+2bjRL2VlhLKyqk1bYYFfioq8sy4shMJCLLlgVtWBZxoFcA4ZNMgPiX7+eV9c6NOf9vHhLNgXIZKXNAsiB33ucz6V7fvfhyuv9OGIa67J2CZA8lGN9SdC2opnoTz5QU3NpEidFr1g23Hkun6gLbk9WEHy+NRL1VhVLQM6EwVwjurUyecVd+jghz+3aQO//nXsqkQknQI4h5n5YkCbN8PVV8Po0X4maZGMVld3XLOzTc2OsBojqQV1fNRLPS41Tzi9g47UDWsMOMeZ+eyIAQPgW9+CtDO/iEhkCuA80LEj3HgjvPsuPPJI7GpEGimZSUFlBVRWVM2oCOVlfqmo8EtZuV9S1ysDIf0oOiv4ZNccSWZUIS3uS1+CXXeF226LXYmIpCiA80RBAZxxhq+4tmZN7GpEBBTAeeWzn/VPcJMnx65EpBnVGJqouoTK7V/SRTpYQwGcRw46yL++9lrcOkTEaRpaHunUyU8yOm9e7EpEWkEGHGhRH3XAeWbXXX1ZSxGJTwGcZ7p3h5UrY1chIqAAzjtdu8KqVbGrEBHYgQA2s0Ize93MHk+uDzKzKWb2nplNMLM29W1D4isp0VrBIpliRzrgHwIz065fA1wfQtgdWAWc25yFScto397PMyci8TUogM2sH3AccEdy3YAjgIeSh4wDTmiJAqV5FRdXndFFRCJraAd8A3ApkJq93B1YHUJI/SovABpxJjJpbUVFvli7iMRXbwCb2ReBpSGEaY15ATMbY2ZTzWzqsmXLGrMJaUaFhVBZWf/jRKTlNaQDPgz4spl9AIzHhx5uBLqYWepAjn7AwtqeHEK4PYRQGkIo7dmzZzOULE1RUFB10lgRiazeAA4h/CSE0C+EMBA4BXg+hPAN4AXgpORhZwKPtliV0mzMsuIAIZG80JR5wJcBF5vZe/iY8J3NU5K0JAWwSObYobUgQgiTgEnJv+cCBzV/SdKSdGJOkcyhI+FERCJRAOcZdcAimUMBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRNKgADazLmb2kJm9a2YzzewQM+tmZs+Y2Zzka9eWLlZEJJc0tAO+EXgyhDAUGA7MBC4Hngsh7AE8l1wXEZEGqjeAzawzcDhwJ0AIYWsIYTVwPDAuedg44ISWKlJEJBc1pAMeBCwD7jaz183sDjPrAPQKISxKHrMY6FXbk81sjJlNNbOpy5Yta56qRURyQEMCuAgYCdwSQhgBbKDGcEMIIQChtieHEG4PIZSGEEp79uzZ1HpFRHJGQwJ4AbAghDAluf4QHshLzKw3QPJ1acuUKCKSm+oN4BDCYuAjMxuS3HQk8A4wETgzue1M4NEWqVBEJEcVNfBxFwD3m1kbYC5wNh7efzOzc4EPga+3TIkiIrmpQQEcQngDKK3lriObtxwRkfyhI+FERCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpE0KIDN7CIze9vMZpjZA2bWzswGmdkUM3vPzCaYWZuWLlZEJJfUG8Bm1hf4AVAaQhgGFAKnANcA14cQdgdWAee2ZKEiIrmmoUMQRUB7MysCSoBFwBHAQ8n944ATmr88EZHcVW8AhxAWAtcB8/HgXQNMA1aHEMqThy0A+tb2fDMbY2ZTzWzqsmXLmqdqEZEc0JAhiK7A8cAgoA/QATimoS8QQrg9hFAaQijt2bNnowsVEck1DRmCOAqYF0JYFkIoAx4GDgO6JEMSAP2AhS1Uo4hITmpIAM8HDjazEjMz4EjgHeAF4KTkMWcCj7ZMiSIiuakhY8BT8J1trwHTk+fcDlwGXGxm7wHdgTtbsE4RkZxTVP9DIITwS+CXNW6eCxzU7BWJiOQJHQknIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgshtN6LmS0DPmy1F5QdsWsIoWfsIkTySasGsIiIVNMQhIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgk/w+/asBPRYt0dwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}