summaryrefslogtreecommitdiff
path: root/notebooks/plot_fgw.ipynb
blob: b41f280473cb77a47af8d35da557adb4ba275ebc (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
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Plot Fused-gromov-Wasserstein\n",
    "\n",
    "\n",
    "This example illustrates the computation of FGW for 1D measures[18].\n",
    "\n",
    ".. [18] Vayer Titouan, Chapel Laetitia, Flamary R{'e}mi, Tavenard Romain\n",
    "      and Courty Nicolas\n",
    "    \"Optimal Transport for structured data with application on graphs\"\n",
    "    International Conference on Machine Learning (ICML). 2019.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Titouan Vayer <titouan.vayer@irisa.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import matplotlib.pyplot as pl\n",
    "import numpy as np\n",
    "import ot\n",
    "from ot.gromov import gromov_wasserstein, fused_gromov_wasserstein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%% parameters\n",
    "# We create two 1D random measures\n",
    "n = 20  # number of points in the first distribution\n",
    "n2 = 30  # number of points in the second distribution\n",
    "sig = 1  # std of first distribution\n",
    "sig2 = 0.1  # std of second distribution\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "phi = np.arange(n)[:, None]\n",
    "xs = phi + sig * np.random.randn(n, 1)\n",
    "ys = np.vstack((np.ones((n // 2, 1)), 0 * np.ones((n // 2, 1)))) + sig2 * np.random.randn(n, 1)\n",
    "\n",
    "phi2 = np.arange(n2)[:, None]\n",
    "xt = phi2 + sig * np.random.randn(n2, 1)\n",
    "yt = np.vstack((np.ones((n2 // 2, 1)), 0 * np.ones((n2 // 2, 1)))) + sig2 * np.random.randn(n2, 1)\n",
    "yt = yt[::-1, :]\n",
    "\n",
    "p = ot.unif(n)\n",
    "q = ot.unif(n2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot data\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHwCAYAAABZrD3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xec1OW5///XPWULC8vSBBSQrmCHVcFewBhb1IAmscQUNdF8j8nv5GhyjiknzZCck97QmKjJSTRgLLFEwRoRC2BHRRZBqYK7sGydnZnr98d8FmaXmdmZ3Sk7u+/n48EDZz73fD7Xgg+uudt1OzNDREREiouv0AGIiIhI5pTARfoA59wC55xl+VddoX8uEUnOaQhdpPg55yYCNd7LG8zsRxl8tgqoBmYAc4E5cZfnmtnSLMY5w7v/TmAdsMLMdmbr/iL9iRK4SB/hnJsHLPJezjSzVd28TxWwALgKWGxm87MU31Xef9YCFwPzgJ1mNiQb9xfpb5TARfoQ59wiYolxnZlN6uG95gGLzMxlKbaFZnZ13OsqYGJ3v2iI9HdK4CJ9jHOuBphIFnrPzrnriX0ZWJyFuFaa2cye3kdEYrSITaTvaU/a87xedLd5c+lH9zwkANY55xZk6V4i/Z564CJ9kNdzbk+Wk8xsXSHjAXDOzQGWkOWFcSL9lRK4SB/lnFtCbMX3qt4ydO0N7w8FJmj1uUjPaAhdJAnn3Dwv4SS6tsg5tzDfMWVoPrHtWjN6w9C11wNfAFQBtxQ4HJGipwQuktzXiSXADrzV0/MSXetNvB5u+3z49V4CzTvn3Azn3CIzW2pmNwNLic3PzyhEPCJ9hRK4SHIziCWbzqq935fkMZZu8eaab/ZeLvK+fOSNl6QfA66Me7t95OLifMYi0tcECh2ASG8U11tNlKTnwp7kmO79FtKxwlkmFmZSWa0zM7va+3kmEiv0Mre79+qGRcBNnea72//cJuYxDpE+RwlcJLFUSXoekFHxkfgCJgUyl1ip1TnOuRn5KJ7iVV6byN4RACA2tO+cg1hFNhHpJiVwkcTmkCBJt1cPA7rdIy4EM1vnnFtKrChLviqfzQeWplhtnnCBoIikRwlcpBMvSc8gcZJONbTea3n7woeaWdaGz72V7RNTVHurplPv2/tc+9D5qk7vtxedmUtstXr7n/XM7oxgeD9z+5eHKthTmEakT1ACF9lXe+J4McG1DkPrzrnr00kKhZwD9+a/vw5M6O49klhI18PgHyZ4r71W+57pCS/Gpd7w+lBiX5BmAhfRjbly59xK4Mr20QbvaNSsHMoi0lsogYvsq72Xmmjo9yJix2C2G5bODQs1B+71dhcBp2e7cEoa1d1WAIkOVLmavVXi4u/XHt8kYsVndhLrwXfoxXs961XJFhHG1W+P7+FXefGkfR+R3k6V2EQ6iTsMpMO52t6Q8QwAM5vrLdJa0ZtP0/J+lgXe/ut8P3sisJK4qmvtBWXM7IYUn1tC7BS0bsXsnDPiyrV6W9kW9fR0NpHeRglcJI43/10HtJ++Fd/LvMn7/RZiw+s7C5EY0+UdLVrb096/c25ifG87bkvapFSJ2Gs7g9jw/YvEetY1XU0HeAm4W/XbvS8NNfFHoHo97aOzda65SG+hIXSRjtrnqRemGFrt9YkgjQVm6d5nInADsWHv+DO8b/YSbcoE7o1OpB1De3W29uTtPW+ot4q+itjfz8QUXwKq2HfqYy6xaYT2Lx8r0riPSK+nSmwiHWVcpKW38ZLUVcDpWbjdIvZWTgOo9pL3HDqOTnSbc25OXM35izvd96K4nni1dy75xXGfneiNNLTrEJP3haAaWBGXvPe5j0gxUg9cpKM5JC6fWhSytWjN6+0uItZL3TPHH/fF5mo6JvaeqAWWetvIFnrPv4pYb3rPFIWZLW1frR732YnEitPM8Wqt73TOXdm+kM37dQOxZP2i92eS6D4iRUdz4CKeuPnvG4pxaNWLfyWx+Bd31T7FPa4iNm9dBfyo8zx33J/TkHwfCer1tm/oPD/unJuXyc+c7D4ixUQ9cJG9JhKbP+1W8usF2oeSh3o92K5UEdsGV0XsjO4Z7Lvn+q4En7uIvXu25+R5umGiNx++57neqEOmiXif+4gUGyVwEY83VDyk0HF0h5ew9yzAy9Jtk5VdnQksiSsrm09LE/S253VjxCTRfUSKiobQRSQj3vzxXGLbtQq6jc77ElGtXrT0R0rgIiIiRUjbyERERIqQEriIiEgRUgIXEREpQkrgIiIiRahg28iGDx9u48ePL9TjRUREeqWVK1fuMLMRXbUrWAIfP348K1as6LqhiIhIP+Kc25BOOw2hi4iIFCElcBERkSKkUqoFZhahoekedjX8lrbwe/h8g6isuITKis/g9w8tdHgiItJLKYEXkFmYrTs+TUvoOcyaAIhGGmiu/xWtDX9g+PC/Eyg5qMBRiohIb6Qh9AKqb7iDltble5J3lSthrH8g+/kDDHMhXO15RHd9DbPmAkcqIiK9jRJ4Ae1s+DVGLDkP85VS6SvB5xw+5/A7h8Ow5gew2k9jFilwtCIi0psogReIWRuRyBYAAjgqXBCfc/u0c4QgvAZan8xzhCIi0pspgRdMAPADMNAXZN/UHceasKb/y0dQIiJSJJTAC8Q5x4Cy0wGHHx8uQe+7g+iOvMQlIiLFQQm8gIZUXo9zZbRZhGjKc9kd+CfkLS4REen9tI2sgEpLpjNq+J/Z/uHnqaItRcsyXMXleYsrXbWt7/JS7V283/giAGMrjuGooRcztHR8YQMTEekH1AMvsPLS4xg7+jUi5Z/AEn2fcuVQfg4EZ+Q/uBRq6p9i8YYvsqZ+CU2RWpoitaypf5TFG75ATf3ThQ5PRKTPUwLvBZzzU1r1XXxVv4bANGKL2wLg2x8G/ieu8ntdz5HnUVO4lqVbbyJsrRjRPe8bUcLWytKtP6ApXFfACEVE+j4NofciruxUXNmpWLQRiIAb1KsSd7s3dj4ApJqzN1bvepDqYZfmKyQRkX4nowTuYtnkfOAjwAFAaYJmZmYfyUJs/ZbzVRQ6hJS2Nr9OxEJJr0csxNbm1/IYkYhI/5N2AnfOlQAPALG9T7EuWHz30OLelz4s6Cvvuo3ruo2IiHRfJnPg1wNzgB8Co4gl6+8A44DLgU3AnYD+5e7jplbOjUvQRufvbEFXztTBc/Mel4hIf5LJEPrFwEtm9l9A+9xs1Mw2An92zi0HXgauBX6a7UCl9xg/cDaVgQGURLdQ4WsBoA0/dZEB7I4OoiI4ggMrZhU4ShGRvi2THvhEYFncawOCe16Y1QAPAp/NTmjSW9W1PM9wXw0V/hacA+egxEUY4W9gXEkb54/5GT7nL3SYIiJ9WiYJPAw0xb1uAEZ0arOeWKKXPipqbbz6wVeIWss+9dt9zih3DTSEVhYkNhGR/iSTBL4JGBP3eg3QeZz0CEAbgPuw7U1PpDzaNGLNrN/1hzxGJCLSP2WSwJfRMWHfBxzunFvonPuIc+4m4AzgySzGJ71MU9t6ItbSRZsNeYpGRKT/ymQR21+B8c658Wa2nthCtfOBK4HPE1uV/i7wtWwHKb1HwFeJzwWJpuiFB3wD8xiRiEj/lHYCN7PHgcfjXjc652YDFwKTic1/32dmDdkOUnqPkRVnsKb2pqTXfa6UMYMuymNEIiL9U49KqZpZG3BXlmKRIlDiH8rYyst4v/4vRK25wzWHn6CvigMqlcBFRHJNtdAlY1OG/Dt+V86GXbcSO3gFjDYqSw7jsP3+l6BvUGEDFBHpB5TAJWPOOSYNuZbxgz9LbcsLRK2VQSUHMyA4rtChiYj0G0rg0m1+XzkjBpxc6DBERPolnQcuIiJShJTARUREipASuIiISBHSHHiRi1qEN+tXsaruaVojzYyrmMqsYXOpDA4pdGgiIpJDGSdw59whwCeBaUCFmZ3pvT8OqAYeN7OdWY1SEmoM7+Z3a79FXdt2QtFYedN1jat58oP7mD/2Cxw15MQCRygiIrmSUQJ3zn0T+CZ7h94t7nIQWARcB/wqK9FJSrev/xE7WrcQIbznvbC1AbD4/d8xsmwM+5dPKFR4IiKSQ2nPgTvnLgK+TaycajWwIP66dx74SuC8LMYnSWxteZ9NTes6JO94YWvjiQ/uzXNUIiKSL5ksYrsOqAHONbNVQKIjqVYDU7IRmKRWs/t1rMMASEeG8c7uV/MYkYiI5FMmCfxw4J9m1pqizRZgZM9CkuxxhQ5ARERyJJME7oBoF21GAKkSvGTJpEGH4lIkaIdjyqDD8hiRiIjkUyYJfC0wO9lF55wPOIHYMLrk2KiysYwZMAm/d5hIZwEX5NT9LshzVCIiki+ZJPC/ATOdc9cluX49sfnvv/Y4KknL5eP/g+Gl+1PqK9vzXsAFCboS5o+9hv3LxxcuOBERySlnlnwhVIeGzg0AngUOA54jNqR+LPBj4ERgFvAicKJ3TnhK1dXVtmLFim6G3beYGc99uJJ7Nz3CpuatBH0BZg2bwfkHnMnIshEpPxu1CG/tfplVtU/RGm3mwAFTOXbYHAapkIuISFFyzq00s+ou26WbwL2bVhHb4/0JOvbeDbgT+KKZ1adzLyXwGDPjdzV3sPzDlbRGQ3ve9+GjxBfkG4d8hckDxxcuQBERyat0E3hGtdDNbKeZXQqMBs4FrgAuAA4ws0vSTd6y14q6V/ZJ3gBRorREW/mft35L1LpaOygiIv1Nt2qhm9l24MEsx9Iv3b/p0X2Sd7zmSAuv7XqLI6qm5zEqERHp7XQaWYFtat6a8nrYwrzftDlP0YiISLFIuwfunLs5zaZmZld3M55+p9RfQmOkKel1v/NT7i9Lel1ERPqnTIbQP9/FdSO2Mt0AJfA0nTRiFg9sXkrYEtc0j1qUo4cekeeoRESkt8tkCH1Kkl9HA9cQK6N6FzA1yzH2aR8ddRpl/tKEVdVKfSXMGXkilcFBBYhMRER6s7R74N5pY8msdM49DLwKPELs0BNJQ1VJJd899Hr+5+3fsqO1FoiVQY1YhDkjT+LSAy8scIQiItIbZbQPvMubOXcHcLiZHdlVW+0D78jMqGncwIbGjZT6SjhyyCEMDFQUOiwREcmzdPeBd2sbWQrb0BB6tzjnmDxwvIq2iIhIWrK2jcw7zORUQMVcREREciyTbWTHpbjHWOCzwFHArVmIS0RERFLIZAj9GWJbxJJxxA47ub5HEfXAqtq13L5+Ka/vXA/A4VUT+PSEuRw5ZGKhQhIREcmJTBL4D0icwKNAHfCCmT2blai64e/vL+O3ax+kNbr3ILSVdWt5fdcGvjTlXD42JulR5iIiIkUnk21kN+YykJ7Y2lzHb9Y+QCi6bzGU1mgbv3znfmYPn8Z+ZVUFiE5ERCT7+kQt9Hs2Pkuq7XBmcP+m53Iex65QM/e+9zJ/rnmeZz+o0SliIiKSM9neRlYQa3Zvos0iSa+3WZg1uzfl7Plmxs/ffJzb1i4n4HyELULA56fCX8Ivjr2YI4aOzdmzRUSkf0qawJ1za7p5TzOzg7r52W6pDJan0WZAzp7/m7ee5I6a5whFw7QfDBqKRmgKh/jcsjtYfOrVjB84PGfPFxGR/ifVEPoAoLwbv3KXKZM4e/9jKPeXJL1e7i/ho6OPzsmzG8Ot3Lp2GS2RtoTXW6Nhfvf20zl5dlfC0SiPb1zL7W+v4IENb9IcThyjiIgUn6Q9cDMbk89AeqJ66BTGDxjJ2oYttHU61SvoAkwcOJoZQybl5NnLP1iH3yX/HhQx49HNb/LDmTl5fFJPba7hy8/cT1s0SjgaJeDzEcX4+lGncdlBM/IbjIiIZF2fWMTmcz5+OuMqjhk2lRJfgAH+Usr9pZT4AswefjD/e+SVOLfvaV/Z0Bxpo6ty8m0JVsfn0qrtG/nCU39nZ6iFxnCI1miYxnCI5nAbN616nLtrXstrPCIikn19YhEbwIBAGTcd8Rm2tezk1Z3vAnBk1URGlA3O6XMPqhxJNGV9GxhbMTSnMXS24KUnaYkk/tLQHGnjhy89wQUTD8WXoy81IiKSexkncOdcEJgJHACUJmpjZn/pYVzdNrKsirmjjsrb86YOHsn4gcNYs2tbwkRe7g/y+Skn5C2exrYQq7anXnHfHG7j9dqtHD5sdJ6iEhGRbMsogTvnLgd+DCRbUu2IVWsrWAIvhJ8dcxGfeOoWmsIhQtG929nK/UFOGDmZ88cdkbdYmiNt+J0jnGJQwOccTVrQJiJS1NKeA3fOnQH8EfgQ+BqxZP0P4FvAE97rxcBV2Q+zdxtXMZT7T7uWKybPZkTZICoCpUwfPJrvHnUePz16Pr4Ui9yybUhJOaX+1N/LWiNhJlUOy1NEIiKSC5n0wL9KrOb5LDOrd84tAFaZ2feA7znnrgZ+CfwsB3H2esPLBvLl6XP48vQ5BY3D7/Nx+UEz+f2bLyScBw84HyeMnsCI8ooCRCciItmSSddwJnC/mcWf973n82a2EHgO6LU10/uLLx12PNOGjKS8U0+81Odnv/KBLJh9VoEiExGRbMmkB14BbIl73QoM6tTmBeAzPQ1KeqbUH+DOuZewqOYVbn3zRbY172ZwSRmXTJ3BZVNnUFlSVugQRUSkhzJJ4FuBEXGvNwOdS6ZWZnhPyZESv59Lps7gkqkq2iIi0hdlMoS+mo4JexlwunNuNoBzbhpwkddOREREciiTBP4wcLxzrn3z8I+IbRl7xjm3BXiNWA/8+9kNUURERDrLJIEvBA4EagHM7HVgLrAEaCC2lewcM3sg20GKiIhIR2nPV5tZCNjU6b1lwJnZDkpERERSy6SQS2UuAxEREZH0ZTKEvsU593/OuTNcro72EhERkbRkksA3AZ8ktpjtfefcTd7KcxEREcmztBO4mU0Fjgd+DwwAbgBed84975y7xjk3JEcxioiISCcZnbJhZsvN7GpgNLHe+CPADOBXwGbn3GLn3DnZD1NERETideuYLDNrNbO7zOwsYAxwPbAGuBC4N4vxiYiISAI9PufSzLYBK4GXgDZix4qKiIhIDnW7brlzbjLwaeAyYCyxxP0ucEd2QhMREZFkMkrgzrnBwMXEEvcsYkm7AbgduM3Mns56hJJSU1uIe9a8yd/ffoPmSJijRx3AZw6fwfjBWlMoItKXpZ3AnXN3AecCpd5bTwK3AXebWVPWI5MuvV+/k4/f81caQiGawm0AvFO7g7vefI3/PvF0Lp52WIEjFBGRXMmkBz4fqCHW277DzN7LTUiSDjPj0w/ezY7mJqJme95vi0ZpI8q3nnmMQ0eM5JDh+xUwShERyZVMFrGdaGZTzOx7St6F9+KWTWxtaOiQvOOFIhEWvvRCnqMSEZF8yaSQy7JcBiKZeXHrJloj4aTXo2Y8t/n9PEYkIiL51ONtZFIYAefwdVGSPuDTX6+ISF+lf+GL1CnjJuB3yf/6gj4fH5kwJY8RiYhIPnV7H7gU1kHDRnDkyNGs2rqZUDSyz/Wgz89nD59ZgMhERArnvZ07uX/1W2xvbGTysKF8bPo0KsvKCh1WTiiBF7GFZ36MTz9wN2tqd9AcbsOAAcEgDsetZ13A2MrBhQ5RRCQvomZ849Gl3PvGaiJmhKNRygMBfvjUv/ju3NO58NBDCh1i1imBF7HBpWXcc+GnWLl1Mw/WvE1TWxszR+3PuZMPpjwYLHR4IiJ589N/LeP+1W/SGtk7Itkcji30/eaSxxhdOYjZ48YVKrycUAIvcs45qkcfQPXoAwodiohIQTS3tXHbqlV7EnZnLeEwP3vmWWZ/qp8ncOdcADgFmAYMNLObvPdLgIFAnVmSzckiIiJZtmrT5pSLegFe2ryFUCRCid+fp6hyL6NV6M65OcA6YueA/xz4XtzlmcB2YrXSRURE8iIcjXbZxgHRNNoVk7QTuHNuBvAAsV77fwB3xl83s+XAeuCCLMYnIiKS0iEj9yMU2Xc3Trz9Kysp62NrgzLpgX8TaAaqzewnwNsJ2rwIHJmNwERERNIxvKKC0yZNTDo8Xh4Mcs2sY/McVe5lksBPAO4xs80p2rwHjO5ZSCIiIpm56cwzmDBkCAPietkOKA8GOOfgg5h3WP/eRjaQ2Bx3KuWoupuIiOTZoNJS7rn8Eh5Z8w7/99Ir1DY3MWnoUK6onsExY8bguig9XYwySeCbgK6+whwJvNv9cERERLqnxO/n3GkHc+60gwsdSl5k0lt+BDjTOTc70UXn3BnA8cQWuomIiEgOZZLAfwDsApY6574PHAzgnPuI9/puYBvwk6xHKSIiIh2kPYRuZhudcx8B/gZ8HTBiawQe8n5fD1xoZl3Nk4uIiEgPZVSJzcxWOOemAh8DZgHDiPXKnyO2Qj2U/RBFRESks7QTuHNuf6DN62Hf7f0SERHpVULhCEtWv8N9r7xJc6iNGeP255PHHMGowYMKHVpWZdIDfx+4A/hMjmIRERHpkQ/qG/jU7++irqmZplAbAK9u3Mrty1/i++fP5ezD+84K9UwWse0EPshVICIiIj1hZlz1p3vYWr97T/IGCEUitIbD3HjfEtZs21HACLMrkwT+PHBUrgIRERHpiVc3buW92p1EookPxGyLRPjDMyvyHFXuZJLA/xs42Tl3RY5iERER6bYX129MeahJJGo8W/NeHiPKrUzmwE8HHgdudc59gdjBJVuJbSeLZ+1nhIuIiOSLzzkcjn3T0l59qaJqJgk8/uzvY7xfiRigBC4iInk1e9I4fvnE8qTngwd8jlMOmpjnqHInkwQ+N2dRiIiI9NC00ftx8KgRvL55G+HIvkk86PdzxXEzCxBZbmRSie2xXAYiIiLSU7+55GNc8cdFbKyr37MSvSwYwAE/vfgcJgwfUtgAsyijSmwiItJ/hSNRnntjPVs+rGfIoAGccNgEykqDXX8wj4YMKOeeL17GszUbeOC1t2gKtTFj7P5cMOMQBpeXFTq8rFICFxGRDtZs3s7L6zfj9/k47qADGT2kkuWvr+fG3z9EWzhCJGL4/T7MjOvmncS8U48odMgd+HyOE6aM54Qp4wsdSk5lUkq1jVRL+/YyMyvtfkgiIlII2+sbuO7W+1m79UMMw+GImnHEuNG8tXozraG4LVpenZSfLXqKstIA5xx3SGGC7scy6YE/T+IEXgVMBkqB14D6LMQlIiJ51BIKc9nP72Lbrt37FEJZsW4j+I1Eg+UtoTC/WPwvzpo1HZ+vD+3RKgKZLGI7Idk151wl8AugGjg3C3GJiEge/fPlt6lrbE5YxcwAAhD1gS/BDq2W1jbeem8b08ePynmcslcmldiSMrN64HPE/p6/n417iohI/vz9uddpjqsfnkg0yXo1n8/R3Jr6s5J9WUngAGYWAZ4ALsjWPUVEJD8aW0OpGzgHSUbIQ20RDhw1NPtBSUpZS+CeEqDvbLITEeknDh03En+qOqMGLkGZcb/PcfS0cQwfXJG74CShrCVw59wUYD5Qk617iohIflx60gwCAX/S6z7nGOA6Xi8J+hk+eCDfvOKMXIcnCWSyjezmFPcYC5zk/fcNWYhLRETyaMro4Vx75mx+88/ltLSF97zvc46SgJ8fXXYWO7bv5i9LX2L7zgYGDSjl4ycfzkWnHcmgAX2rQEqxcGbpbO0G51zi6vB7rQV+bGa3pHO/6upqW7Gi75zLKiLSF6xYu5HfP/YCr6zfjM/n46TpE/jc6UczedTwQofWbzjnVppZdVftMtkHPiXJ+1Ggzsx2ZnAvERHphaonj6F68phCh9HrtLS08dgTq/nHQy/T2NjKxAkjmHfh0Rx2SOH+rDLZB665bRER6XfqdjZy7XV/om5nEy0tse1ymzbX8fyL6zjv7KO45urTChJX2ovYnHM3O+fO6aLNWSnmykVERIrOf3//Pj7YvntP8gYwg9bWMP946GX+tWxNQeLKZBX654EZXbQ5ilhBFxERkaK3aVMdb761hUiC88UhNrT+p788m+eoYnKxDzzBTkEREZHi8+bbW/D7U6fKde9+kKdoOso0gSddsu6cCwInAtt6FJGIiEgvEQj4SFXfBsDny3ZfOD0pF7E55zoP7F/nnLssQVM/sB8wANAcuIhIgYRawzz9+GqW/+ttMDj2hKmcfPp0SsuSFDKXlGYceSDhcPKBZefgmOqJeYxor65WoQ9gb6/bgCBQnqBdBFgDPAb8d9aiExGRtG14dztf/eLthEJhmptitc1XPF/Dzb9Ywo9/fRkTJo8scITFp7KynDPnHsYjS1+ntTW8z/WSkgCXX3JcASLrIoGb2Z4Nbl4hl/81s+/kPCoREclIa0sbX/3i7eza2dTh/eamEM2E+I9r7+BP915HeXlJgSIsXv/vmjk0NLbyzLPvEIlEiUSilHkjGjd+7VymTinMMaqZFHKZC6zLVSAiItJ9Ty19I2EPsV0oFOHJJW/w0fOOymNUfUMg4OcbXz+P9zfW8sRTb1K/u5kJ40dw2inTKC8r3BeiTAq5PJbLQEREpPuWPf02Lc3JjwRtaQ6x7Mk3lcB7YOyYoVx+yfGFDmOPTHrgwJ7V5jOBA4DSRG3M7C89jEtERDJg0a7PtYim0UaKR0YJ3Dl3OfBjIFlVe0dssZsSuIhIHh17whReXvkuLc1tCa+XlQeZdeLUPEcluZRJKdUzgD8CHwJfI5as/wF8C3jCe70YuCr7YYqISCqnfeSwlOd5+/0+5px5eB4jklzLZPf5V4E6YJaZ/dh7b5WZfc/M5gBfBM4H3sxyjCIi0oXy8hIW/PIyBg4s27NCGqCsLEjFwFIW/PIyBlQknPWUIpXJEPpM4D4zq497b88XADNb6Jy7BLgROCtL8YmISJqmHDyaP993HUsffpVnn34LM5h94kHMPfsIKpS8+5xMEngFsCXudSswqFObF4DP9DQoERHpngEVpZw372jOm3d0oUORHMtkCH0rMCLu9WbgoE5tKunGynYRERHJTCbJdjUdE/YyYL5zbraZLXfOTQMu8tqJiEga1r2/g7sfeZn1Gz9k+JCBfGzO4Rw1fQyuqxM0pN/LJIE/DPzUOTfazLYAPwLmAc845z4g1jv3Ad/PfpgiIn2LmfHbv/yLRQ+/RDgcIRI1HPDMihoOPWh/fnzD+ZQENaApyWUyhL4QOBCoBTCz14mVV10CNBDbSnaOmT0VD/VUAAAgAElEQVSQ7SBFRPqaJcveYvE/X6I1FCbiFVgxoLm1jVff2sTPb3+yoPFJ75d2AjezkJltMrPWuPeWmdmZZjbFzOaa2UO5CVNEpG/5/aJnaUlSu7w1FOahJ9+gsak14XURyKwHLiIiWdDYHGLLtvqUbQIBH2+u25aniKQYdacW+iHAJ4FpQIWZnem9Pw6oBh43s51ZjVJEpA9xjljtylSs6ybSv2XUA3fOfRN4BfhP4AJic+DtgsAi4NKsRSci0gcNKCth7KiqlG0ikSjTJ4/OU0RSjDKphX4R8G3gcWI97QXx182sBlgJnJfF+ERE+qQrLz6estLEg6ClJQE+NvdwyuNKoop0lkkP/DqgBjjXzFYBLQnarAamZCMwEZG+7NRZU7n8/GMpCfoJ+GP/FPuco6w0wLFHjOfaS04qcITS22UyB344cFv8KvQEtgAjexaSiEj/cMXHZ3HGidO4d8krewq5nHPaYUyfPKrQoUkRyCSBOyDaRZsRxGqki4hIGvbfbzDXqLcNwPs1H7Bq2TuYGYcdPZFJ0/cvdEi9WiYJfC0wO9lF55wPOAGVUhURkQw01Dfz/f/3Z1avWo/Fatrg8zkOnDKSb/3uCoaO6HxulkBmc+B/A2Y6565Lcv16YvPff+1xVCIi0i+YGf95xS28vuJdQq1h2kKxX60tbdS8uZmvfvK3tIUSF7zp7zJJ4D8DXgN+4pxbBnwEwDn3Q+/194EXiZVcFRER6dIrz9Xw/rrthNsi+1yLhKPU7djNs0veKEBkvV8mpVSbgFOI9bCPBWYRmxe/3vvvO4EzzKwt+2GKiEhf9MT9L9HSFEp6vaUpxJK/r8hjRMUjo0psXoW1S51zXwGOAYYBu4DnzWxrDuITEZE+rKkh0Y7kjppTJPj+rFtn1ZnZduDBLMciIiL9zLSjDuTFp96mtSXx4G0g6GfakePyHFVxSDmE7py73Dl3eL6CERGR/mXux6tTFn33+RznXnpc/gIqIl3Ngd8GnB//hnPu0865x3MWkYiI9BuDBg/gaz/5FKVlQXz+vSnJ+RylZUGu/fYFjDxgSAEj7L26M4Q+Hjg5y3GIiEg/Nev06fz87i+x6PdPseKpt7Gocfixk5h/1SlMPWxMocPrtbo1By4iIpJNB04ZxVcXXFzoMIpKRseJioiISO+gHriI9HtrVtaw/P4VhENhph93EMecdRR+v7/QYYmklE4Ct5xHISJSAA07G7nxnJtY+/J6Qi0hLGqUDyqjfGA5P3zkRiYcqu1L0ns5s+T52TkXJfMEbmbW5ReD6upqW7FC1XVEpDDMjOuOv5G1q9YlrLU9cEgFt6/5JZXDdJCG5JdzbqWZVXfVLp05cJfhL82ri0iv9/aLa3n3tQ1JD8oItbTx0O8fy3NUIulL2VM2MyVjEemTnr3vRVpTlOgMNYd4/K//4hM3nJ+0jWTu9eXv8Jf/fYg3nnsHgENmTeZT/342h86eUuDIio8WsYlIvxRqCZFqChGgrUXHWGbTQ7c/zc3fWExr894vTi899RarX1jHld+Zx9lXnFTA6IqPetgi0i8dctzBlA8qS3rdH/Bx+MnT8xhR37Z9Uy0L/2tRh+TdrrU5xM03LmL7ptoCRFa8lMBFpF+afV41JWUlSa8HggEu/PLZeYyob3vwtqdTjniYGQ/+8ek8RlT8lMBFpF8KBAP88J83UjF4ACVlwT3v+/w+SgeU8KVffY4Dp6mMZ7asffW9pAsGAdpCYd559b08RlT8NAcuIv3W5KMm8Me3f8GDCx/liTuXEW6LcPhJ05j3/53LgdPHFjq8PmVQVUUabQbkIZK+QwlcRPq1IfsN5tJvzOfSb8wvdCh92pxPzOa5R16lpbE14fWyilLmflLHhmZCQ+giIpJzR518MAceNJpg6b79xmBpgHFTR3HUyQcXILLipQQuIiI55/P5+MHdX2bmqdMJlgYoqyilrKKUYGmAmadO56a/fwWfTykpExpCFxGRvBgwsIxv/ekaPthYy+vLY4VcDp09hf3GDC1wZMVJCVxERPJqvzFDOW3+sYUOo+hpvEJERKQIKYGLiIgUoZTHieb0wc5tBzYU5OEiIiK914FmNqKrRgVL4CIiItJ9GkIXEREpQkrgIiIiRUgJXEREpAgpgYuIiBQhJXAREZEipAQuIiJShJTARUREipASuIiISBFSAhcRESlCSuAiIiJFSAlcRESkCCmBixQx59wC55xl+VddoX8uEemaDjMRKWLOuYlAjffyBjP7UQafrQKqgRnAXGBO3OW5ZrY0a4HufeYM7zk7gXXACjPbme3niPQHSuAiRc45Nw9Y5L2caWarunmfKmABcBWw2MzmZynE9vtf5f1nLXAxMA/YaWZDsvkckf5CCVykD3DOLSKWENeZ2aQe3msesMjMXFaC23vfhWZ2ddzrKmBid79wiPR3SuAifYRzrgaYSBZ6z86564l9GVicleBi91xpZjOzdT+R/k6L2ET6jvakPc/rRXebN5d+dM9D6mCdc25Blu8p0m+pBy7Sh3g95/YkOcnM1hUynnjOuTnAEnK0QE6kv1ECF+ljnHNLiK30XtXbhqy9Yf6hwAStPhfpGQ2hi3icc1XOuZXOuTpvP/RVSdot8BaN9VbziW3TmtGbhqy9HvgCoAq4pcDhiBQ9JXARj5nt9HqsN3lv7TPM662cvp7YHuZeyevZts+HX+8lzoJxzs1wzi0ys6VmdjOxP9d53p5wEekmJXCRfR1NbH9yoiTd3iu/K4/xZMybY77Ze7nI++KRd16Sfgy4Mu7thd7vF+fyud5oiuYIpc/SHLhIJ14p0aWJtmJ5c7jrzGxuGvdZSMfqZplYmElVtRQxtG8tW5pOzNnmPb/Dz+J9magjB8ViOj37KmB+IX5ukXwIFDoAkd7EK01aRWy1dOdr84glw7QSQnzRkgKaS6zU6hzn3Ix8Fk3xEuhE9o4EALEhfuccxCqy5dJcEvw9ivQVSuAiHbX3mDvMf3u9xluAq3vT1qyumNk659xSYqMG+a54Np9Yzz/ZavOaJO9nyxz2rmcQ6XM0By7S0VwSz38/RiwZ3ZzgM72Wty98aIFGA6qBfb40eKMcxF9zzs1xzl3VedW8c25Rdxa7xT1jonNunrdzYGLKD4kUGfXARTqaA6yIf8PbMrYu0/naQs+Be6vPvw5M6Ml9Utz/emJFWVJNKXyY4L32mu1LvftUAVVmdrO36OyGTm2vTHCPruwZSfGG7NcRO/ClV+2LF+kJJXART+f5by+xPAbc1Z1kWsg5cO9nWQScnsOCKTfTaX67kxVAooNVrmZvtTiAajNb7H3h6NArJzYasjPuvSrgljS+TM0ldrxq+2d3Ejs2dY8M7iXSK2kVuojHW3S1kNg//jOIbSe7oZjmvNt5q78XFHLI3/sSsZK4qmvtQ+RmdkOC9guBmvYvS14P/+juJFhvJ8HM9r87bwHiLTq6VPoS9cBF9ppJrEDLXLrZ6+4NvCH/Hs/XO+cmJvry4iXDocQSZNJRBm8B3enALc65F4n1xmtS/LleRMch7rnsPec8k7jbh+TjY59LgsI8IsVMPXCRPsTr4c7paQ10r/d8Q+cE7Y1SLPWScx1ZqmnuPa8m/gxybz58Zvvq+bhtfItTjYp491oSfy66F+vpmd5LpDfTKnSRPsKbM74KOD0Lt1vE3opp8Wq95F1FbL1AttQSm6cG9nxRIC7hziHWg040lz0xvja9l5Dj77UAuDmde4kUE/XARfqAuPnmPb3Mbt6niljyrk41X+wlwYXxvdye8ua8IZZY5xIbBp/bqc1C77mdF7stIlZ1rX1l+zxiaxg+JLYQbp/phET3EikmSuAiRc5LuiuJDXkv7sE9riK27awK+FGihWZx7RcQS7A5WWnvHYm6pPN8uXNuZbLpAefcvEx+/lT3EikGWsQmUvzah4+HJjsCtZMqYJj3+1Biw8idi5x0dVhLVqucxc9Re18mqtl7olq8hOVXvRGITOeyc13KVSSnlMBFipiXsNuLliSas+6OdMquziBLq7q9hH1z3DNvAa7svDiu8z7xTuZlsmugi3uJFAUlcJEi5s3t5nWvt1fadGfnBNtdXqW0Gu/LyCRi89KJvhwkPJzE+wKQaTLWQSdS9DQHLiIZ6UmBlW4+b55Xqa3Hc9bZvJdIoakHLiKZupj8nvJ1sXNuKN2riZ7Le4kUlHrgIpKSN0T9LrEqaVXAomxuHxOR7lEPXETSsZTYYrlJ6EQvkV5BPXAREZEipFKqIiIiRahgQ+jDhw+38ePHF+rxIiIivdLKlSt3mNmIrtoVLIGPHz+eFStWFOrxIiIivZJzbkM67TSELiIiUoS0Cl1ERIqeRbZCy6NgjRCYCqUn41zfTnF9+6cTEZE+zSyM1X8bmu8FHNAGrhwogSG/xZX03SPfNYQuIiJFy+q/D833AyGgFYjGeuFWh9V9FgunNZ1clJTARUSkKFm0FpoXAy1JGrRijXk96yevlMBFRKQ4tf4L8KdoEIGWf+YrmrxTAhcRkeJkLUC0izZteQmlEJTARUSkOAUPAedStwlMyU8sBaAELiIiRckFDwXfGJKmMleOG3h1XmPKJyVwEREpWm7Ib8ENBko7XSiHsgugdG5B4soH7QMXEZGi5QLjYMQ/saa/QvPdYE0QmIqruBJKjselGGJvDr3Bh7v/SGu4hqB/FEMHXkpF6XEpP9ObKIGLiEhRc74huIHXwMBr0mpvZmzZ+S1qG/+KWQiIAFDf8hgVJUdz4Ig/4nMlOYw4OzSELiIieReNfEjj7l+w88NL2VX7BVpbHsUskpdn1zXe5SXvZtqTN4BZE42tz7N153fzEkdPZdwDd87NAj4DHAUMBnYBq4DbzOy57IYnIiJ9TUvzQ+yu+zfAaC/C0tb6OD7//lQNuxuff1jOnm1mfFD/My95J7hOC7WNdzJy8Nfw+ypyFkc2ZNQDd879FFgGXAlUA1O8368CljnnfpL1CEVEpM8It73N7p3/BjQTX0HNrJFIeD27aq/I6fMj0Z2EI9tStnEEaGl7PadxZEPaCdw5dw1wHbCBWAKfAgzyfr8KeA+4zjn3xRzEKSIifUBTw2/BQkmuthEOv0k4h8nTOR+WVsveP8OcSYTXAFuBajO71cxqzKzR+/33wNHANuDaXAQqIiLFL9T6FPHzzvuwCKHWZ3L2fL9vMCWBsV20ilBecljOYsiWTBL4JGCRmdUmumhmO4C7vXYiIiLdkF7/uCdGDv4PnCtPeM25coYO/Cw+V5bzOHoqkwReS+ystlRagA+7H46IiPRlJaUnkfIAEhcgWHpCTmOoGnAuIwZdi6MUaN8u5se5MirLP8Kowf+R0+dnSyar0O8DznPO/aeZhTtfdM4FgfO8diIiIvsYMPCLtDY/QOJh9ACBwMEEg4fmPI6Rg7/MkIqPU9vwf7SGawj4RzG04pOUl0zP+bOzxZmlN1zhnBsMPE5s29jXzOyFuGvHAj8ktqjtNDOr7+p+1dXVtmLFim4FLSIixaul6QF277yO2HC5N7DrKvD7RlI1/B58/uGFDK/gnHMrzay6q3ZJe+DOuTUJ3i4FjgSWO+dage3ACPYWod0IvAgclHHEIiLSL5QNOIeS0mNobvo/2lpfwPkqKCufR0nZHJxTgdB0pfqTGkDi1QSb4/7bR8c5b5/3ORERkaR8/v2oGPSV2LitdEvSBG5mY/IZiIiIiKSv9+9UFxERkX0ogYuIiBQhrRYQEZFexyzKhobHeaPuz+xu20jQV8GUyvM4qGoepf7BhQ6vV1ACFxHpA5raNrCxfjFN4fWUBUYzZtA8BpZMLXRY3WIW5aktX2dL0/OEzTutLNrA63W38/auuzlr7B+oCI4qcJSFpwQuIlLEzIy1dT/lvfo7MItghAE/m3YvYlTFWUwf/l2cS3+2tC60idW7HqMxXMvQkrEcMngu5YHK3P0ACbyz6342Nz1PxFo6vB+xENHITp7eeiMfHfv7vMbUGymBi4gUsc0N9/Be/Z+JWnyl6whRi7C18SHKA2OZOOQLXd7HLMqjW3/B6l2PYRYhSoSAK+Ff22/jtJFf5IghZ+Xuh+jkjbo/7ZO898RJlLrWd6gPvUdlybi8xdQbaRGbiEiRMjPW1f2aqDUnvB61FjbU/4GotXV5r2Xb/8Sbux6P9XK9MqdhCxGxEE9s+x3vNryY1diTMYvSEN6cso3DT13rO3mJpzdTAhcRKVKtke2EoqnPjzKL0hBKnezaoq2srP07YUt8XlXYWnlm+x3djjMzDl+qw04A5xwBX+LTxPqTjBO4c+4Q59z3nHN3O+f+Gff+OOfchc65quyGKCIiiUXTbJf6zIstzW91OU/+Qcta2qJdHUjZc845xgw8EXBJ25hFGVk+I+ex9HYZJXDn3DeBV4D/BC4A5sZdDgKLgEuzFp2IiCRV6t+PgK+rWqRGRXByFy0iuBQJE8DhMEt0glj2HT708/hdacJrflfGIUMuJ+Dr/ed151raCdw5dxHwbWInklUDC+Kvm1kNsJLYkaIiIpJjzvmYMPgL+Fzi4WSfK2NM5afw+xInw3b7lU4mbKGUbQYGhxPM07D1kNJJnL7/Tyj1VRFwA/ARJODK8bsSplV9gsOGXpGXOHq7TFahXwfUAOeaWatz7twEbVYDJ2clMhER6dLYyk/REHqLLY0PELUQsWF1h9+VMaRsFpOHXNflPcoDlUwddAJrdj9DJMGCt4Ar5dhhn8C51L30bBo5YAbzJj7I5qbnqA+9R4lvEGMHnqgiLnEySeCHA7eZJVnlELMFGNmzkEREJF3OOaaP+C4HVF7M+/V/obltA6WBUYyt/CRVpTPTTrpzR19HbWgjta0baduzqt0RdKUcXHkyh1d9NHc/RBI+52dMxfFQcXzen10MMkngjq5XTIxgz+nsIiKSL4NLD2XwiB90+/MlvnIuGf9z1u5+lpfq/kFTZCdDS8Yyc+gFHFB+SF5735KeTBL4WmB2sosutoTxBGLD6CIiUmR8zs/UyhOZWnlioUORNGSyCv1vwEznXLIJleuBKcBfexyViIiIpJRJD/xnwEXAT7wV6Q7AOfdD4ERgFvAisDDbQYqIiEhHaSdwM2tyzp0C/Ar4BHt779cTqxJwJ/BFszRq9omIiEiPZHSYiZntBC51zn0FOAYYBuwCnjezrTmIT0RERBLo1mlkZrYdeDDLsYiIiEiadJyoiEg/VBvaxit1y2iI7GJk6RiOqDqBUv++ldbMjJZoCwHnJ+grKUCkkkzaCdw5d3OaTc3Mru5mPCIikkNRi3DvpltYVfcUUTOihClxpfxj8x+ZN+ZajhhyvNcuymMfLOGfWx+mIbwbM2PiwElceMA8pg46qMA/hQA4s9Sn1Oxp6FxXRVyM2Mp0M7PUZ8EB1dXVtmLFirSeLSIi2fHwlj/z7I6HaUtQVDPoSvjcxG8wbsBB/KbmV6ze9TqhTjXSg66Ez034PNVDj8lXyP2Oc26lmVV31S6TfeBTkvw6GriGWBnVu4CpGUcrIiI51xpp5tkdDyVM3gBtFuLRrXfxys6XeLP+jX2Sd3ubP66/ldbI3ntELUpzpIWopXu8qWRDJtvIalJcXumcexh4FXiE2KEnIiLSi6xveguf86c8HvzdxjeoC5XQ2sXZ36t2rmDaoMO48/2HeOKD54hYBL/zc8qIY/nkuLOpKqnMcvTSWUbngadiZhuA+4AvZ+ueIiKSPZE0z/PeHtqW8nprtJX3mt7nKy/fxNKtzxKKthGxKKFoG49tW85XXv4Bta07sxGypJC1BO7ZhobQRUR6pQPKJxCxcMo2w0pGMTCQuvccdEFerltDfVsDETp+KYgQob6tkVvW/a3H8UpqWUvg3mEmpwL12bqniIhkz+DgMCYPPBx/ktnTElfKqft9nFNGnEpJii1jhrGuYTvRJAdURomysu4NGsNNWYlbEks7gTvnjkvy6yTn3CXAEuAo4P6cRSsiIj1y0dgvMbR0JCW+0rh3HSW+Uo6oOoEZQ05m9rDjGRocit/tu6GoxFfCsUOPx+8LpnxOwBdgh4bRcyqTQi7PkHLpAw54llhtdBER6YUGBAbxb1N+zKu7nuW5HY/QFGlgROkBnDjiHCZWxM79LvWX8vVp3+C29bfy2q7XCPoCmIFzcPbo8zh6yGwe3vJyyueEo2EGBSvy9FP1T5kk8B+QOIFHgTrgBTN7NitRiYhIzgR9Jcwccgozh5yStE1FoIJrJ/8bu9p2sbHpfYK+IBMrJhHwxdLGhIoxrGlYn/Tz4ysOYGjJ4CxHLvEy2UZ2Yy4DERGR3mdwcDCDB++biD83cT7ffP3ntEb33Ste4gvyuYnz93k/HI3w9PbX+Pv7z1Ab2s3+5cOYP/ZEjhl2MM65nMTfl6kWuoiIZGzqoPF885Br+eU7f6YutAu/8xGxKENKKvl/ky/joEETOrRvjrTylZW/Y0PTNpojsaS/uflDXt+1nuohU/nvwy/H77K9MapvUwIXEZFumV45md/M+BbrGzfyYWgXw0oGM75iTMLe9C/evpd1jVsIRTtuY2uJhHix9m0WbXiKT4w/NV+h9wlJE7hzbk0372lmpkr3IiL9gHOOCQPHMoGxSds0hJt5bNtL+yTvdq3RNu5870kuOvBkfOqFpy1VD3wAqVedi4hIL2JmPLjpZf5Q8xTvNX6I3+fjpP0O5qrJpzKlclTB4lrfsI2gCxAieRGZpkgrdaEGhpWqBGu6kiZwMxuTz0BERKT7zIwbX1nME9tW0xJpAyAajfDE1tUs276Gn828lGOGTypIbAGfH+uiPxg1I+jTrG4mNFYhItIHPP3BWzwZl7zbRTFaIm1c/9KdtEXTq4WebZMH7p+wKEy8cQP2ozI4IE8R9Q1K4CIifcAd656huVPyjheORnj6g7fyGNFeAZ+fyyecTlmS6m2lviCfn/zRPEdV/DIer3DOBYGZwAFAaaI2ZvaXHsYlIiIZWN+4I+X11kgbGxpSt8mleWNPoi7UwKL3/wVmtFmEUl8Qw7hmyrkcN3x6wWIrVhklcOfc5cCPgeHJmhBb+KYELiKSRwMDpdSFGpNeD/oCDAyW5TGijpxzXDX5bC4cewKPbX2JHa31HDBgGKePPIpBGjrvlrQTuHPuDOCPwNvA/wALiB1csgI4BTgNWAw8kvUoRUQkpfPHzuTmd56gNclWrSjG6aMOyXNU+xpeOpiLDzyl0GH0CZn0wL9KrOb5LDOrd84tAFaZ2feA7znnrgZ+CfwsB3GKiEgK88Ydw53rn6M21EjEOh7zWeYPcuHYoxlWOrBA0fVMbUsTi9a9wmu1WxlcUsYFEw5l5vDEBWPamRmPbH6T3771L9bWbyfg83HyqClce/BJHDR4ZB6jzx1nlt5Wb+fch8B9ZvZZ73UU+I6ZfTuuzdNAg5md1dX9qqurbcWKFd0KWkRE9rWteRc3vHQnb9dvIeCLrfoOR6NcMmE210ydU5RFUv6x/g2uf/5BHNASCeOA8kCQw4aO5tZTLmJAYN9zy82M77zyMPdseKXDwj4fjlJ/gF/PuojjRxZmS106nHMrzay6q3aZ9MArgC1xr1uBQZ3avAB8JoN7iohIlowsH8xtx13N+obtvF2/hVJ/kGOGTWRAIOF6417vlQ83c8PzD9Ia2TstYEBTuI2Xd2ziq8v/wW9O/Pg+n3tu+7v7JG+ITSM0R9r4t+cXs/zsf6fEX9z7zjOJfiswIu71ZqBzydTKDO8pIiJZNn7gCMYPHNF1w17u168v65C847VGIzy+aS1bmuoZPaBj9bY/vPNcyi11hvHo5rc4Z+yhWY033zIZT1lNx4S9DDjdOTcbwDk3DbjIayciItIjy7auT1m/LeDzsXzrhn3ef6f+g5T3bQyHqKnf3sPoCi+TBP4wcLxzbrT3+kfERjOecc5tAV4j1gP/fnZDFBGR/qir8qtYbFi8s4HB1FMGJT5/l22KQSYJfCFwIFALYGavA3OBJUAD8ARwjpk9kO0gRUSk/5kxPPWRHBEzqkfs22b++BmU+RNXfYtxnDmm+AvHpD1fbWYhYFOn95YBZ2Y7KBERiXltx1Zufu0FXtmxlYpgCfOnHMb8KYcyqKT4e5Bd+dKhx/PSjk0J57ODPj8zR4xh/KCh+1z7+IFHcuuaZ2mLRhJuqTtrzCEcMKAqZ3HnSybbyCrNrD5bD9Y2MhGR1H758rP86pXnCEUjRL1/q8v9AQaWlHLPOZcydtDgAkeYe7e++QL/8+qTRKJRwl4yrgiUsP+ASu6aexlVpeUJP7etuZ5/e34xb+3aig8fzsW21F00fgZfO/wMAr7eu6Uu3W1kmSTwRuBe4HZgiaX7wSSUwEVEklu+5T0+s2QxzeF9V2H7nGNK1TAeveCzBYgs/9bVf8gda1bsKeQyf+IRzBkzhaAv9QlnADW7d/Ba3WZKfQGO228Cg0sSJ/zeJBf7wDcBnwQ+AWxxzv0JuMPM3uxmjCIiksSvX1meMHlD7Ozs93bv4pXtWzhixOiEbfqSiZXD+Hb1R7r12UmDhjNpULLjOzJT19LMG9s/wO/zcdTIUZQFUs2z514mc+BTvS1jVwDzgRuA651zK4j1yv9qZnU5iVJEpJ95ZcfWlNej0Sirtm/uFwm80BpDIf7rqaU8XLOGEn+s1x8x43NHzOQrxxyHL0VJ11zKaBLAzJab2dXAaGK98UeAGcCvgM3OucXOuXOyH6aISP8S6KLsqXOOkjSGkKVnwtEon7z3bzxcs4bWSITdoRC7QyGa2tr4/csr+K8nlxQstm7N4ptZq5nd5dU8HwNcD6wBLiQ2Ty4iIj3wkQOn4E91WAfGKWMm5jGi/unRdWtZu7OW1khkn2vN4TB/f3s1G3btLEBk3Uzg8cxsG7ASeAloI3YmuIiI9MDVhx2TtFZ3qd/PnLGTOWBgZcLrkj1/fv1lmtqSl2WNmHH3W2/kMaK9up3AnXOTnXPfdc6tB5YClxXuKlAAACAASURBVAMbge9kKTYRkX5rwuCh3Db34wwKllLhLZbyO0eZP8AJ+4/nJyd1eeijZMH2pqaU18PRKNubGvMUTUcZHTzinBsMXAx8GphFrLfdQGwR221m9nTWIxQR6admjR7Hyk9dy8Pr1/Dajq0MKinlrPEHMXVIdlZVS9cmDRnC2roPkxZ1LQsEmDxkWF5japd2AnfO3QWcC7SX/3kSuA2428xSf0UREZFuKfUHOH/SdM6fVPylP4vRZ4+YyVPvbaA5nHgY3cy48ODC/N1kMoQ+n9he8G8BE8zsdDP7k5K3dLatoYEXN25kzY4d9LDej4hIQR09+gDOnXIQ5YF9+7tlgQDfOuFUhpQVpjhMJkPoJ3q1z0US2lRfz9cfeZQXN22i1O8nHI0yvKKCb512KqdO1GpZESk+zjkWnHoGR40cza9WPs+2xgbMjOnD9+Pfjz2eUw6cULjYCtVDUinVvuWDhgbOvuNP7GppIdLp/6myQICfn302cyZPKlB0IiI9Z2Y0trUR8LmcVmFLt5Rq763mLkXlV889T31r6z7JG6AlHObGJUv2HMYgIlKMnHMMLCkpeAnVdkrgkhV/X72acDSa9HpTuI2V/397dx7nRH0/fvz13iR7cSgIcgseqAioHCoqCkqp900p3kcrHvWobe23WtuirW3VCrWXZ621/hQFtYKiqBxSVFAUlYpHQRG5FHA5dzfZJO/fHzOBISTZZDfn8n4+HvNIMvP5zLwzm807M/OZz2flyqTLjTHGZMYSuGm2SDRKXYqODgAEYX0j91MaY4xJnyVw02y+sjJ2q6hIWSaqSte21muUMcZkiyVwkxXnH3rItlF6EulQXU3/Tp3yGJExxrRslsBNVow97DC6tmmTMIlX+f3cffJJSIGG3DPGmJYo4wQuIn4R+ZaIXCsiN3nml4tIe7Fv6V1Sm4oKnj3/PL7bvx9VgQCVfj+BsjKG9tyLp84dw8CuXQsdojHGtCgZ3QcuIt8CHga64fSDrqrqc5cdCcwFzlfViY2ty+4Db7lCkQg1dXW0Ki+ndXl5ocMxxpiSkvX7wEVkIPA8Tu9tNwI7JGlVfRNYBpyVUaSmxSn3+ejUurUlb2OMyaFMulL9JVAHDFbVVSLyqwRl3gYGZCUykzMNkQjzln1JTW0de7XbnUO6dbbr08YYU2IySeBDgWdVdVWKMssBG6S2iP37g8X85uXZRKIKKKrQvlUVE846mUO6dSl0eMYYY9KUSQJvDaxtpEwVu2DL9uXrN/DSok/YVB+kd6cOfLtvb6rKi6OrPa+p//2YX06bQX04vMP82g0NXPTYZCZdei7772njDBtjTCnIJIGvBPo2UuZQ4POmh9M86zfX8o/ZC/j32x+ypT5Ix7atueCYAYw56hAqApm81fSEI1F+/sx0Xv7wf0SiSjgapao8wK+nzuTu757CsAMKN0pNvEg0ym9fnr1T8o6pbwgzYfbr3Dv6jDxHZowxpikyOVqeDpzotjbfiYh8Gzgap6Fb3q2u2cTZdz/K43MXsrG2nkhUWbNhM39+6Q0u+uuT1DckTlzN8dsXZvHK4iUEw5Ft/YDXhRqoDTVww8Tn+Wj111nfZlMtWvVV0uQNoMBrS5bREInkLyhjjDFNlkkC/y2wEXhVRG4HDgQQkRPc108DXwHjsx5lGn4+cTobttbTENlxQI1gQ5ila9bz4Kvzs7q9DbV1PPPuh0l/GITCEe6bld1tNsfmYJCyNBqqBVMkeWOMMcUj7QSuqiuAE3CS9E3AaJx7wae5r9cCJ6lqY9fJs251zSY++GJ10uEqg+EIT7zxPtFo9oazfHPpcgIpug6NqjLn02VZ215z9WrfjlAjR9fV5QGq7dYvY4wpCRldGFbVBSKyP3AGMATYA+eofB5OC/VQ9kNs3Odrawj4fQTDyRNUfaiBLcEgbasqs7LNhkiUxjrBiaQYXjPferTbjb6dO/HeysQ/dCr8Pi4YfEhaR+nGGGMKL+0ELiJdgQb3CPtpdyoKrSvKkx59x0RVqcxiQ7b+3To1mqD371RcLbrvOuNERj38OJuDoR3G7q70+9m3Q3uuOPrwAkZnjDEmE5lcA/8SuDNXgTRHvx6dUyZnAY7avyfl/uwl8L07tqdv1074yxIfsVYF/FwxvLgSYo92uzF17IWcO+hgWleUUybCnq1bcc2xQ3ji4u9SFSi+W9+MMcYklnZf6CKyHnhIVf8vGxvOdl/oL7zzEeMmv5qwUVllwM9j147hgK4ds7Y9gLWbt3Lu/ROp2VpLnbtdcbf3ncP6838nDbMezowxxmQk3b7QMzkknU8Rd5N6yqA+hCIR7pzyGqpKVJ1k2qaqgrsuODnryRugY5tWTLnuIqa+9xGT3l7E5vog+3fuwCVHD2JQr25Z354xxhgTk8kR+BHAHOAKVX2kuRvO1WhkDeEIb/5vORu21tGtfVsG7t3NjoKNMcaUjFwcgY8AZgJ/F5ErcQYuWYPTB4iXqurvMlhvVgX8Po7tUzw9oBljjDG5kEkC/43n+eHulIgCBUvgxhhjzK4gkwQ+MmdRmJwKNoSZ9/EXbNxaT88923Hw3l3ssoIxxpS4tBO4qs7IZSAmN578z/vc8+//IAgaGz60TRV3XnYKfXt2LnR4xhhjmij7Q3S1QOFIlP+8v5TZC5cSjUY5om8vRh62f05GOMumyXM/YMKzc6gP7Xhr3cr1DVz+p8n8vxvPY+/O7QsUnTHGmOYo7gxUBFat28jYO55i09Z6aoMNALy2cCkTJs7mbz8ZxQF77VngCBNriES457m5OyXvmPpQmHunvcmdl52S58iMMcZkQ9o9sYlIg4iE0piCuQw4n8KRKGPvfIqva7ZsS94AtcEGNm6t58q7JrGltjjf7sIlK1N2LxtVZeb7Sxrtz90YY0xxyrQjl0Tf9rsD+wEVwCJgUxbiKgpz3/+MTVvrkybChnCEqa9/yLkjB+Y5ssZtqQ/RWDO1SDRKOBpNOaqaMcaY4pRJI7ahyZaJSFvgT8Bg4LQsxFUUXntvKbX1DUmX14fCvLrg06JM4Pt0br/T2Ojx2reutuRtjDElKpPBTJJS1U3A93CO0G/PxjqLQaSRBAg0OgpaKsFQmHCKIVCbo1en9uzbZY+kw4NWBPxccHzx/fAwxhiTnqwkcABVjQCzgLOytc5CO7J/L6orko/QVRHwMbR/Zr2+qSrTZixizJUP8u0xf2TE6Alc+dPHWPD+F80Ndye/v/RkWldV4Pft+GeuDPjp02NPLjjOErgxxpSqrCVwVznQLsvrLJgRg3pTHkh+irmsrIyzhh2c9vpUlT/c+woTHniVlWs2EI0q0ajy4aer+dntz/DCq4uyEfY2e3Xcnck3X8h3hh5Mm6oKfGVCl/Ztue6MoTx43SgCfjt9bowxpSrtwUwaXZFIb2AesFpV+zVWPleDmWTb/75cy9g7n6IhHNl2S1ZFwE9ZmTDhujMYfOBeaa/rg8Ur+PGtk6kPJr6uXl7u59m/X0nbNlVZid0YY0zpyfpgJiLyQIp19ACOdZ9nZbzwYtG7R0em3vl9pr3xEa8u+IRIVDm6/96cNaw/7dpUZ7SuSVPfIRhK3ihOgJdfW8yoUwc1M2pjjDEtXSa3kX2/keVLgLtU9aFmxFOUWldVMHrEoYwecWiz1rNsxXpSnfAIhsJ8vnx9s7ZhjDFm15BJAu+dZH4UqFHVDVmIp0Xbfbdq+DJ5gvb5ymjfrlUeIzLGGFOqMrkPfGkuA9kVnHnioXyyZA11Se4t9/nKOGH4QXmOyhhjTCnKpCvVB0Tk1EbKnJziWvkub9iQ3nTv0o5AgpbtlRV+Rh7Th+5d8teI//MV65k+9yNmvfVp0XYJa4wxJrFMr4GvAJ5PUWYATocuY5sTVEvl9/v48+1juPOv05n71pJtiTwSiXLOKQO5/Pxj8hLHmnWbuPmPU/lsxTp8ZQIiRMJRRp0wgKvHHENZmY0VbowxxS7bo5GVA7npWqyFaFVdwa03nk7Nhq18svQr/P4y+h3YjcoUHcZk06Yt9Xz/F49Ts7mWaHTHFnVPv7KQUEOYH118fF5iMcYY03SZduSStA21iASAY4CvmhXRLqLd7q0YMmgfBh/SK2/JG+DZGe+zpbZ+p+QNUB8M89zMD1i/YWve4jHGGNM0KY/AReTTuFnXi8iFCYr6gD2BasCugRexqbMWEWxIfpJERJj99v84Z2TzbpkzxhiTW42dQq9m+1G3AgEgUTdhEeBTYAZwa9aiM1mXanQ1cIZItQZtxhhT/FImcFXtHnsuIlHgblW9LedRmZzp2bU9NZtqky6vLPfTq2v7PEZkjDGmKTK5Bj4S+FeuAjH5cf6pg6lKcc3d7/Nx9MB98xiRMcaYpsikI5cZuQykJQrWNzBvxmK++XoTHbrsxhHH9aE8jw3WEjl6wD4cd0RvZs7/lPpgeNv8MhHKy3387ken7zT8qDHGmOKT8W1kbmvzQUA3oCJRGVV9vJlxlbyXn36be389BRGhIRQmUO7s6ut+fTbDT81PA7E1X21kydKvKS/3cXC/HlRWBhARbrniRA7v14t/Pjef5au/we/3ceygfbnkrCHs071DXmIzxhjTPBkNJyoiFwF3Acm+5QVQVW10oOlSGU60KV5/eRF3/eRJggkajFVUBrjlLxcy+NgDdlq2emUNy5evp1V1BX36dsPnb9qRcE3NVm6/43kWLfpye2cxUWX0qMO5+IKjraMWY4wpYrkYTvTbwD+AT4A/AHcAU4AFwHDgeGAyML0J8bYYqspDv38hYfIG57T6Q3e8sEMCX71qA7+/7d8s+fQrAgEfqorf72PsD0ZwwimHZLT9uroQV1/3L9au20wkEiXkuWXsyUnzqd0a5AdXjWjamzPGGFM0MjnE+wlQAwxR1bvcee+q6m9U9VvAVcCZwEdZjrGkrF6+npp1W1KWWbV8Peu/2gTAN+u3cM3lD/Px4lWEQmG2bg1SWxti06Y6/jz+JV54bmFG239x+iI2bNhKJBLdaVkwGOa55xeyfn3q+IwxxhS/TBL4IGCKqm5KVF9V7wfmAbdkKbaSVF8bavTUt99XRn1dCICnHn+T2q3BhD2jBYNhHvjbq4RC4Z2WJTN12ns7NE6LJyLMnvNx2uszxhhTnDJJ4K2A1Z7XQaBNXJm3gCOaG1Qp69xjD8IpejoDUIWOXXYDYPq0DwiHdz5a3l4Y3n3787S3v2VzfcrloVCYLVtSlzHGGFP8Mknga4COntergPiWWG3J/gApJaW6dQXDTjmUQHnidnzlFX5Gjhq87XayutpQyvVFVdm0qS7t7XfvnroTlqqqAN26WUctxhhT6jJJ4IvZMWG/DowQkSMBRKQPMNott0u74pbT6LLXHlRU7njPd0VVgB777sllPz5p27wOHeNPYuxIgK7ddh4jfP3azcyc9gGvTH2PVV9+s23+6HMOo7Iy+b3mgnDs0P3TfCfGGGOKVSZHyy8CE0Ski6quBu4ERgFzReRrnKPzMuD27IdZWlq1ruSep6/l5clv89y/3mDj+i2069CGMy8ZysizB+3QmcvZow/n4ftnEUxy3bpN2yr69t/Woy3B+gbG3/ocr89cjN/vQ4FoOMpBh/bg53eMZsgR+3LUkfvxxptLqI9rCV9R4efnPzuN8vJd+iSJMca0CGnfBy4i5ThJep2qBt15RwO/APYFlgETVHVaOutryfeBZyIUCvOTax9j6ZKvCHmSuAhUVAS444/nc1C/btvm33z1oyx694sdygL4Az66dm/PvU9eRZmvjGkvvc/jT85nzZoNiAiDBvbi0ouG0ufArnl7b8YYYzKX7n3gGXXkkk2WwLcLBcM89cQ8np30Flu31APCkKN7c8n3h9Frn+3NDj75cCU3fv8fSe8xr6ou50fjzuTYkX23zQuHI5SVlVnnLcYYUyKy3pGLyZ3yCj8XXDKU8y8+mvr6BsoD/oS3os184X0aUtxSVlcb4sVnFuyQwP3+RjvFM8YYU4Ka0hd6X+BcoA/QSlVPdOfvBQwGZqrqhqxGuYsQEaqqypMu37SxLuH94l6bN9ktYsYYsyvIqLNtEfkl8D5wM3AWzhCjMQFgEnBB1qIzO9jvwC47tWz38vnK6N2nSx4jMsYYUyhpJ3ARGQ2MA2biHGnf4V2uqkuBd4DTsxif8Rh5+qFOLzBJ+P1lnHnukDxGZIwxplAyOQK/HlgKnKaq7wKJztUuBnpnIzCzs7a7VfPjW8+kojKAxLVJq6gMcN7YYfTcd8/CBGeMMSavMrkGfjDwSOwWsiRWA52aF5JJZdgJ/enUrR1PPDiHhfOXElXlgL7dOO/yYQw6cr9Ch2eMMSZPMkngAqTotBtw7hNPleBNFhzYrzu33nNeocMwxhhTQJmcQl8CHJlsoYiUAUOxrlSNMcaYnMskgT8FDBKR65Ms/ynO9e8nmh2VMcYYY1LK5BT6H3EGKxnvtkgXABH5PXAMMAR4G7g/20EaY4wxZkdpJ3BVrRWR4cBfgDFsP3r/KaDAROAqVU3cz6cxxhhjsiajntjcHtYuEJEbgMOBPYCNwHxVXZOD+IwxxhiTQJP6QlfVtcALWY7FGGOMMWlK2YhNRC4SkYPzFYwxxhhj0tPYEfgjON2nfhCbISIXAxer6vG5C8uUiq2b65gx6S1em/IukXCEgcMO5NSLjqF9p90KHZoxxrRoTTmF3gsYluU4TAla9vEqfnrOPTSEwtTXhgD4bPFKnrl/Fj9/4DIOO75vI2swxhjTVBmNRmZMTLghwk3f/QubN9RuS94ADcEwwboQt499mHWrbVRZY4zJFUvgpknefOkDgnWhpMujkSjP//M/eYzIGGN2LZbA86RuSz1TH5zJNceO43sDfsbvLr2XT975rNBhNdnCuZ9QtzV5t/cNoTALZlmvusYYkyvpXANPPgC1Scvald/wwxG/YcuGrQTd082rPv+aedPeY9QPT+LCm84scISZ8/ka/+2XThljjDFNk8437DgRicQm4JcA3nlxUzi3IZeeW8/9EzVfbdyWvAE0qgTrQky+50XenflhAaNrmiNG9qOqVUXS5eWVAYaeNiCPERljzK4lnQQuGU522OXx2aLlfPnJaqKRxCOxBmtDTLz7+W2v62uDvPXiQuY+O5/Vn3+VrzAzNnDYgbTbsy1lSY6yAwEfJ45JOnidMcaYZkp5Cl1VLRk308cLGr/O/em7n6OqPHrrJCb9YQo+v7Pbw6EwfYbsz82PX0/7zu1yHWpGysrKuGPSdfzfqD9Rs3bTtuvhVa0q8Ad83P7ED2jTrlWBozTGmJarSV2pmvT5Az6kTFKW8fl93H/jozx/3ysEa3dsGPbfuR9z7ZCbeXDReKrbVOUy1Ix16LI7D8y5hXdmf8QbL75PuCHCgGMO4JhTB1BeGSh0eMYY06JZAs+xQSP6EQlHki4v85Ux8PiDmPK36TTU7zyQWyQcYeO6Tbz08EzOvv6UXIbaJD5fGYeP6MvhI6zTFmOMySc7RZ5je3Rpx9AzBic9Ig2U++nWqyNlkvwoPVgb4vn7X85ViMYYY0qQJfA8uOGvlzHoW/0orwpsu75d2aqCqtaV/GridQQCPkIJjr69tmyozUeoxhhjSoSdQs+D8ooAv3r8Or74eCVznn2b2k117Nt/L4456zAqqsrZvHYjla0qqNtSn3Qd3ffvmseIjTHGFDtL4HnU88BuXHhTt53mH3nGYcjY+5LWq2xVwXd+fFouQzPGGFNi7BR6ESivCHDLxBuoqC4n/lJ4RXUFQ04bzJBTBxUmOGOMMUXJEniROOzEAYx/7TYOP3kgZb4yRKDLPp24esIl3PTYdUiKRm7GGGN2PaJamK7OBw8erAsWLCjItoudqhKNRvH5fIUOxRhjTJ6JyDuqOrixcnYEXoRExJK3McaYlCyBG2OMMSXIErgxxhhTggp2DVxE1gJfFGTj+dcBWFfoIArM9oHtA7B9ALYPwPYBpN4HPVW1Y2MrKFgC35WIyIJ0GiS0ZLYPbB+A7QOwfQC2DyA7+8BOoRtjjDElyBK4McYYU4IsgefHA4UOoAjYPrB9ALYPwPYB2D6ALOwDuwZujDHGlCA7AjfGGGNKkCVwY4wxpgRZAs8iETlRRD4RkSUi8rMEyytE5El3+XwR6ZX/KHMrjX3wIxFZLCIfiMgMEelZiDhzqbF94Cl3joioiLS422nS2QciMtr9LHwoIo/nO8ZcSuP/YC8RmSUiC93/hZMLEWcuicjDIvK1iPw3yXIRkT+5++gDERmY7xhzLY19cL773heJyBsickhGG1BVm7IwAT5gKbAPUA68DxwUV+Zq4D73+RjgyULHXYB9cBxQ7T6/alfcB265NsAcYB4wuNBxF+Bz0BtYCLRzX+9Z6Ljz/P4fAK5ynx8ELCt03DnYD8cCA4H/Jll+MvAiIMAQYH6hYy7APjjK8z9wUqb7wI7As+dwYImqfqaqIWAicEZcmTOAf7rPJwMjpGWNE9roPlDVWapa676cB3TPc4y5ls7nAODXwB1AfT6Dy5N09sHlwF9VtQZAVb/Oc4y5lM77V6Ct+3w3YFUe48sLVZ0DfJOiyBnAo+qYB+wuIl3yE11+NLYPVPWN2P8ATfg+tASePd2ALz2vV7jzEpZR1TCwEdgjL9HlRzr7wOt7OL/AW5JG94F7qrCHqr6Qz8DyKJ3Pwf7A/iLyuojME5ET8xZd7qXz/scBF4jICmAacG1+QisqmX5ftHQZfx/6cxSIMSmJyAXAYGBYoWPJJxEpA8YDlxQ4lELz45xGH45z1DFHRPqr6oaCRpU/5wKPqOrdInIk8C8R6aeq0UIHZvJPRI7DSeBDM6lnR+DZsxLo4Xnd3Z2XsIyI+HFOna3PS3T5kc4+QES+BfwcOF1Vg3mKLV8a2wdtgH7AbBFZhnPtb0oLa8iWzudgBTBFVRtU9XPgU5yE3hKk8/6/BzwFoKpvApU4g1vsStL6vmjpRORg4CHgDFXNKB9YAs+et4HeIrK3iJTjNFKbEldmCnCx+3wUMFPd1gstRKP7QEQGAPfjJO+WdN0zJuU+UNWNqtpBVXupai+c616nq+qCwoSbE+n8L/wb5+gbEemAc0r9s3wGmUPpvP/lwAgAEemDk8DX5jXKwpsCXOS2Rh8CbFTV1YUOKp9EZC/gGeBCVf000/p2Cj1LVDUsItcA03FaoT6sqh+KyG3AAlWdAvwd51TZEpyGDWMKF3H2pbkP7gJaA5Pc9nvLVfX0ggWdZWnugxYtzX0wHfi2iCwGIsCNmR59FKs03/+PgQdF5AacBm2XtLAf84jIEzg/0jq41/p/BQQAVPU+nGv/JwNLgFrg0sJEmjtp7INf4rSD+pv7fRjWDEYos65UjTHGmBJkp9CNMcaYEmQJ3BhjjClBlsCNMcaYEmQJ3BhjjClBlsCNMcaYEmQJ3Jgi4o5OpiIyPG5+L8+yXgUJzhhTVCyBG2OaTER2F5Fx7rR7oeMpBBE51P1h9WXjpY3JHuvIxZjS0AB84nleLHbH6ZwC4BFgV+nL3Cs20thzBY3C7HIsgRtTAlR1JXBgoeMwCVkCNwVhp9CNMaaJ3L6sB+AMDTy7sNGYXY0lcFOURKSHiNwpIu+JyEYRqRORpSLynIhcJCKVCer4ROQyEZkpIutEJCgiK0VkUnyjsCTbHO6WXenWXSciM0TkUhHxJakzzr3+Odt9fY6IvCwiX4tIVETGxZVvJyJ3ue+lXkRWu9sc1EhsSRuxuXGriKj7ej8ReVhEvnTfxwoReVBEEo61LCJlIjJCRP4kztjcK0QkJCLrReQ1EblSRAIJ6s0GPvfM+twT47Z9ElenXESuFpFZ7v4Nicga9+96Uor3XyUiPxGRN0WkRkQaRGStiCwWkX+KyDmp9l+C9d3rxrghWaNAEbnKLRMWkWOTrCp29P2iqja4f99at97oRmL4tVvuM3E7wjYmI6pqk01FNQEXAnU4gzwoEATW4Vz7jc07NK7ObsAsz/IwUANEPfPuSrHN8Z5yUbdu2DNvBtAmQb1x7vLZwN2e+t+49cd5yvYClsW9r42e56d7lg2P204vz7JeccuGe5YdB2x2n2+K22crgW4J3oN33erW3xA3bw5QFVfvGZwRtGJl1gJrPNMzceV7Av+N28/x27k3QXxtgPcS/H28721Zhp+xKk8sbwD+uOX92P4ZvDXFel51y3zXM+8Rd96rKer5cIZUVeDmQv/P2VSaU8EDsMkm7wScwvakOxdngPsyd1m5+/oB4KC4epM9ifBaoNqd3xlnFLjYF/2VCbZ5jWf5/UBnd34r4IeeRDExQd1xnqSnwO+Bju6yCqCn+9yHM8yk4iT378SSBnCQmyBrPHEMj9tOL8+yXnHLhnuWfYNzLfZAzz4bjZPMFXg0wXvoDjwGnAa098xvDVyCk/gVGJ+gbtK44sq1Aj5yy80ChgEV7rLdgBs8+/D6uLq3uPPXA2d76pUBXXF+8D3QhM9aP5xRsBS43TPfm9znAr4k9Xd3PxtBoK1n/hFs/6GxT5K6p7llGmKfN5tsynQqeAA22RSbcBpVfuZ+sf0HKE+zXuwLU4GxScrEEvxaoNIzv8pNDAo8nqTutZ71D4pbNs6z7O4UMY72lBuRYHk1zrCKzU3gM3F/8CR5D7XEHW2msX8Hu3W3ePddY3HFlfsF289UBJKUOcvzN/J75k9z59+Ug8/cle66I8Bx7rz73Hk1wF4p6p7nlpueYNlCd9nvktSd6i5/Ohf/SzbtGpNdAzfF5Dhgb/f5DaoaSrPed93HFcBDScr8wn3sAIz0zB8JtHefj0tS92/Aavf5eUnKRIE7UsQYG/v9dVWdEb9QVWuBO1PUT9dvVTWaYH6shXQV0DuTFarqAuBrnKPoQ5sY1/fcx/Gqmuw2uH/jnCnoAHjbBMRuTevSxG0npc6YzM/gHM0/JiJjgSvcxZer6vIU1VO1Pr/Xfbwkvv2A2xYhdr3//iYFbgzWiM0UbAVo2gAABnNJREFUl6PcxzVu0kjXYPdxVpLkhap+hHMq2Fve+/xLVf00Sd0IzpFtfF2vJar6dRoxzkxRJtWydM1PMn+V53n7+IVu47Ir3QZ4q9zGb9sapAF7ukW7ZxqQm7B6ui//7jZa22nC+ZHU2i3X07OK593Ha0TkCRE5U0Q6ZBpHCt8HluOcjo8l1IdUdXKyCiJSjpOEFZiSoMjjOJcEOuOcLve6DOeSyufAK82K3OzSLIGbYtLZffwiw3qx5LIyZSnnCN1bvrl1vVIl73S3syLFsrSo6uYk88Oel/FHhHsCC3COGkfiHOlGcRoOfuVOsR9GrZoQVlfP8w5ApxRT7Dup2hP748A9OMlyDPAssFZE/icif22sBX9jVLUG+IFn1mfA9Y1UOx6ncd07qrrT301Vt+C0KwAYG5svImVsPxvxoKpqU+M2xhK4KSal/GUWKXQAzTAB6I/TFuAyoIuqVqlqR1XtrKqd2X4E35Tbnby34PVRVUljesS7AlX9IXAAcDPwIs5p9f2Aq4EFIvLHJsTldbnneTd33amk03lL7DT6SM+tat/GObsQBv6RWYjG7MgSuCkma9zHnilL7Sx29NvY6d3Ycu/RcnPqZiJWL+G92Gksywn3+uzZ7strVPUfqromrowP58i5qbzry/Rvu42qLlHV36nqycAewJE4180BrheR05uyXhG5BucWvgiwGOfugYkiUp2kvLjlIUUCV9VFOLeoeY+6Yz8Unovfz8ZkyhK4KSZvuI+dRSTZteZEYtfLj3NPUe5ERA5ke4J8O0Hd7iKyf5K6PpwGdvF1M7EtxhRljm/iupujIxDrFGdhkjJDPWXiedscJDw6V9VlbL90EH89uElUNaqq84BRONevYcfGiWkRkf7AXe7L24CTcY7u++CcmUjkMJzLAp+5STqV2FH4ZW5bgNj7fyDTWI2JZwncFJNZONcfASa4DYXSMdF97IbTICmR29zHdTidb8S8gnPqGJK3Qr+C7ddxn0gzpnhPuo9DJUGvcCJSBdzYxHU3R+z+cIBD4heKiB+4vZH6MalGI3vQffyeiAxIFZCItI97XZGsrNvAMHa3QsIGjCm2U4Xz2anEud/7dlX9gu3XrMcm6eEtk77PJ+F8vrriNGwLYI3XTJZYAjdFw/0yjnWqMhSYISJDY0fVbkvp4SLymIgc5Kn3FvC0+/LPInJN7PSniHQWkQdxOk4B+IWq1nvq1rE9cZ8rIveJSCe3brWIXAfErq8+qarvNPHtPQ28G3suTperPnc7fXCu63Zs4rqbzG1s9br7cryIHO/Z3/1w7sEeDGxNUn8D24+uL3UTfiJ3A4twkuUs92+0R2yhOMOSniQij+L0AeA1X5xuXoeLSCtPna4i8me2X6+elubbjpmA04nOBuB89/OHqk7C6fwH4EER6RFXL+0ErqpBnJ7ZAGLdsVrjNZMdhb4R3Sab4ifgIqCe7R2E1JNeV6qzPcsbcHola2pXqt/EbW8mjXSlmsb72gfndK/3fcW6Es1KV6qNbD/ZugfhdNLijSt2ZN6A09PZMvf1JQnWe0tc3eVu+Ylx5boCb8bt5xq2dycbm/4XV29Zgjpb4urs1EtcI/vibE/dUQmWV7O957g5uL2xAfu689aRpIe2BOvaz/M5tJ7XbMraZEfgpuio6qM4Q2f+EadRURinA5IvcBotXYjz5eqtsxEYgdNYaDbOPbitcRpQPY3Ty1bSU9Sq+iOca9BP49w21dpdxyycltkjNcktWhm8r89wOkIZj3MaVXAS3mTgKFVNdD9xzqlzVuFw4CmcxFSG896fcuP6VyOr+C3ObVcLcBJUd5zGap29hVR1Fc6ZlXNx7p1ejZMoy3GS9FScrmvjBw4ZgzPm+Ayc/VaOcyr6C5xLEyPcv19a3CPqWIc/f9cE93ur07HOuTg/rI7B+ZEC24++X1D3iL0xqroEpy93sMZrJotE1c7kGGNMOkRkDk5CP0dVn0mzTmfgS5yugk9Q1ZdzGKLZhdgRuDHGpMHt/e0onLMm0zOoeiVO8l6CNV4zWZSswYkxxpgdtQd+g9PtbsJGffHc2yF/7L4cr3bK02SRnUI3xpgsE5FlOB3CxNoBLASO0OQDuRiTMUvgxhiTZe4AMOA0onwJ+JmqflXAkEwLZAncGGOMKUHWiM0YY4wpQZbAjTHGmBJkCdwYY4wpQZbAjTHGmBJkCdwYY4wpQf8f5+iMYpa8N5cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x504 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% plot the distributions\n",
    "\n",
    "pl.close(10)\n",
    "pl.figure(10, (7, 7))\n",
    "\n",
    "pl.subplot(2, 1, 1)\n",
    "\n",
    "pl.scatter(ys, xs, c=phi, s=70)\n",
    "pl.ylabel('Feature value a', fontsize=20)\n",
    "pl.title('$\\mu=\\sum_i \\delta_{x_i,a_i}$', fontsize=25, usetex=True, y=1)\n",
    "pl.xticks(())\n",
    "pl.yticks(())\n",
    "pl.subplot(2, 1, 2)\n",
    "pl.scatter(yt, xt, c=phi2, s=70)\n",
    "pl.xlabel('coordinates x/y', fontsize=25)\n",
    "pl.ylabel('Feature value b', fontsize=20)\n",
    "pl.title('$\\\\nu=\\sum_j \\delta_{y_j,b_j}$', fontsize=25, usetex=True, y=1)\n",
    "pl.yticks(())\n",
    "pl.tight_layout()\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create structure matrices and across-feature distance matrix\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%% Structure matrices and across-features distance matrix\n",
    "C1 = ot.dist(xs)\n",
    "C2 = ot.dist(xt)\n",
    "M = ot.dist(ys, yt)\n",
    "w1 = ot.unif(C1.shape[0])\n",
    "w2 = ot.unif(C2.shape[0])\n",
    "Got = ot.emd([], [], M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot matrices\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAFgCAYAAAAcmXr5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuQHWd55/HvM2cu0uhmSyMLYUu+YeK7BWidIhjWbIDCbCo2bIXY7LJOQsVkC1I4IQlekkpcSW2VN8UtISkvBhxMikDYAoKLpQi2wWvMErBky7Is+Y4vkmVdLNm6zoxm5tk/TgvP6ffpUc/MOdPnnPl9qlQz5z3d/b490+dRz/v2+7zm7oiIyNzqqboBIiLzkYKviEgFFHxFRCqg4CsiUgEFXxGRCij4iohUQMFXRKQCCr4iIhVQ8BWZ58ysz8z+wMx+ZmYvm9lRM9uYlfVX3b5uZZrhJjJ/mdnJwJ3A2cBngR9nb10B/B7wfnf/ekXN62oKviLzlJkZ8APgXOCt7v5I7v31wIvu/vMq2tfteqtugIhU5lrgcuCqfOAFcPcNc96ieUR3viLzlJltBnrd/fyq2zIfacBNZB4ys9OBi4CvVN2W+UrBV2R+uij7umWqjcxsjZndZWbbzOxhM/vrrK9YZknBV2R+WpZ93XWC7caAj7n7ecDrgF8G3tPKhs0XGnATmZ92Z19fPdVG7r4T2Jl9P5r1E69pcdvmBd35isxPPwEOAL8dvWlmlwVlK4CrgH9tbdPmBz3tIDJPmdnvATcDtwP/COyhPtniN4Cl7v6mSdsOAN8DvuPun6yguV1HwVdkHjOzK4E/ot6fC/AscA9wq7v/LNumBvwz8Ky7/2ElDe1CCr4iMiUz+wJQA37HFTCaRsFXRAqZ2ZuAe6k/kjaeFd/q7n9bXau6g4KviEgF9LSDiEgFFHxFRCqgSRbSMczsncDfUB/8+YK73zTV9kNDK/yMtWtnXqFPRK2INix/zIngmNFs3Wi7orrGx8rVXdTMsdFyxyzqohwrV79Hxyw6z+CYPhHUHxSNHgnOh7j50W/z4Nh4UlYrmFF9cDxt/14m9rr7ynCHSRR8pSNkjzv9PfB2YDtwn5nd7u5bi/Y5Y+1aNtx794zr9JEjaWGtLy2bSD+s2RHSoqOH0rK+YLGI4cPxIYNg5Qf2BtsFdRe00/c+nxa+FBxzeDhu0/5g26ie/fvTwiPBzxiY2LsvKRuPgmpwns/evz085uhIEFR70z/+7919IClbXIs7Ce55+WhS9jkOPRNunKNuB+kUlwJPuPtT7j4KfA24suI2icyYgq90ilOB5ya93p6VNTCz68xsg5lt2LP3xTlrnMh0KfhKV3H3W9x9vbuvXzm0ourmiBRSn690ih00ZtM6LStrHQvuTaKBl8L0tkF5LfjIRfX0FHw0LRig6h1Iy6LBwqK+6f5g/4GF8baRhYvK1b8w6DMuGHDrGVyQFkY/52AUbdFg0C8P9Ab9trXe9Jgn9daSsgUFv+MVwbaUHP/Una90ivuAc8zszGw586upJ4QR6Ui685WO4O5jZvZh6ukMa9SnuD5ccbNEZkzBVzqGu38X+G7V7RBpBnU7iIhUQHe+IkV6ogG3qGwaM9yiY/YEgzZRWf2NtCgaMIumcxUNuA0MpmUL0skDhRaOlKrfhtNjFv7kRtJj9kQ/u6CewcE4rPUGg2vRJIulwcDcQPR7B1b1a8BNRKSjKPiKiFRAwVdEpAIKviIiFdCAm3Qvn0gzkxUMnEQDYdYbZBsruW9hk3qWpPUE+3vZugHrj2ajRekX4+EtW3xyWjgaDI4Vpa48kmYBi/jhg2ndUeY4gIMvpWUjwQy54JyWrl4dH3M0yIpWSwfM3vLAtqSsZyCeNXfB47uSso/cl2Zki+jOV0SkAgq+IiIVUPAVEamAgq+ISAU04CZdzNJlf4rSPxYNxDVbMMssHAYrmo0WDZqNHYs2LLcvxINrR9LBMY4VrI0WDbhFKSUPvZyWBbPeAHgpSIQfzHoL69lbsKxRNOAWDHYe3ZEud1Trj0Pl7t3TmAmYr3rGe4qIyIwp+IqIVEDBV0SkAurzlY5hZk8DB4FxYMzd11fbIpGZU/CVTvNWdy8YUcnzdOCqcMAtGIyaxsy10oJUkeEMt4mCuqPmB7O0QkUDbtFsugXBumy98Swviwa9gro8qr/gmESz6fqDGW6RZcvi8pIz3AZWLE7KehbE7Tz5QDAI+MyUrXvlmOU2ExGRZlLwlU7iwPfNbKOZXVd1Y0RmQ90O0kkuc/cdZnYKcIeZPeLu90zeIAvK1wGsPe20KtooUorufKVjuPuO7Otu4FvApcE2t7j7endfv3Jo+Vw3UaQ03flKRzCzRUCPux/Mvn8H8Jcn3jM/yFMw4NZuigYGo0Grom3L7Fu0f9kyKFjXbhb1QLyGXTSwODGNgdKoPDim9UV1x8es1WZ+PSn4SqdYBXzL6h/WXuCf3P171TZJZOYUfKUjuPtTwCVVt0OkWdTnKyJSAQVfEZEKqNtButfEBBw91FhWK7jko1lmwXprYarHaHCoQDSbLdyuqJ0BLzsTr2jALajL+gbS7SaCmWwAC4OfU2RJsFbcWJymkpPSNJUepbQMZtdZ0Yy/aP/gdzcQ/Zz64hluQ0PBuW98JK4/X3WprUREpKkUfEVEKqDgKyJSAfX5Svcyg75cxq6i5YLKZhuLqmlB9rMwA1j9jbQsygBWdl+I+12j/tGCpY38WMnlffL97wCjcaYyj5YcitoZ1bNvT3jMsn2+4TJEBX2+I8+nSw6VpTtfEZEKKPiKiFRAwVdEpAIKviIiFdCAm3SviQkYPtxY1lM0ySIdePFoeZ1g0Kl4yZ8041U0eSIaXLPCDGJpuVvQzqLBtZLHDCejFEyysOjnFOlfkJYVDRZGEzfGgoG9IKtZOBkDSi8jxPDRpMgKBtwG+qNz/0lcf47ufEVEKqDgKyJSAQVfEZEKKPiKiFRAA27SVszsVuDXgN3ufmFWthz4Z+AM4Gngve5eYmqRp4NEVpCZq+x9SLiMzzS2LX3MWS53FO0/nWWEopMqXPIn+NmF5xQtN1Twcw+XDArKot2Lssz1BoNm0TJCUd1FWeYKBuLK0J2vtJsvAe/Mld0A3OXu5wB3Za9FOpqCr7SVbCn4fbniK4Hbsu9vA66a00aJtICCr3SCVe6+M/v+BeqLaYbM7Doz22BmG/bsm3nSE5FWU/CVjuL1GQmFnanufou7r3f39SuXBysniLQJDbhJJ9hlZqvdfaeZrQZ2l9prfAw/kEsP2BssjwPQn5Zb/8J0u7FjaVnRsjXRbLRocCqY5RXOWitQOBuuRHvqdQWDRtFAWJS+EeLBqGjAbVqz5oI2FaS0TIwGM+EAH09/dxbNbDwSpL7MpyY9vn/UzpJ05yud4Hbg2uz7a4FvV9gWkaZQ8JW2YmZfpT45/pfMbLuZfQC4CXi7mT0OvC17LdLR1O0gbcXdryl461fntCEiLaY7XxGRCujOV7pbPuVg0aBROPurbFkLTGs22lyZRt2zbmfJWXfR76MwHWfJGXbR/oWz+2Z+nrrzFRGpgIKviEgFFHxFRCqg4CsiUgENuEn3ctJZUUWzpKLyaNCrbFkz9i97zFYMwk0nJWXZgbjomEUDVtFMwGhwLZxJF884tGjbKP1kNGutaCZbwcy3MnTnKyJSAQVfEZEKKPiKiFRAwVdEpAIacJPuNTaK732+sSxIHQnAwGBSZIuDfMCjR9Oy3oJBl2iAKUqrODZabt+C8jAl5DSEKSktSLVYtN5aODhWcruClJLhQFi4bTqIZkuH4mNGsxujn+fw4XS7ogG3WQx26s5XRKQCCr4iIhVQ8BURqYCCr8yYmWnMQGSG9OGR2ThsZg8DD2b/NgEPuvuMlw02s1uBXwN2u/uFWdmNwO8Ce7LNPu7u3z3hwcbH4KXcGm4DwbpsAAuCgbRgcM2PHAz2XRQfMxiMsb5gwO9YMOAWDczVjxAUlbyHKhzEK1iDLtm9/P4eziaL2l4w4BbOhivYNq9/QVweDdhFv6MlwUBrNAAIePS7K0nBV2bjWuBi4BLgj4HVgJvZDrJAnP3b5O5PlDzml4C/A76cK/+0u3+iGY0WaQcKvjJj7v414GvHX5vZEPVAvC77+uvAnwB9ZnbY3ZeUOOY9ZnZGSxos0kYUfKVp3H0vcFf2DwAz6wcuoH6HPBsfNrP/CmwAPlrUtWFm1wHXAawdCv58FGkTGnCTlnL3UXd/wN1vm8VhbgbOpn5HvRP45BT13eLu6919/cqlBX2xIm1Ad77S9tx91/HvzezzwHdK7gjDwzOvd3wsLYwGWKYz+yka9InSWRbN/ApTPUbblk8JGc1cKxxcKynaP659GuutlVUwOBYPVgZlUZrIgmNa0YzJEnTnK23PzFZPevluYEtVbRFpFt35Slsxs68ClwNDZrYd+AvgcjNbR/3m6Wngg5U1UKRJFHylrbj7NUHxF+e8ISItpm4HEZEK6M5XutfYGOzPzXBbWPAExMKRtOzIgaTIgzILB7yIB40Wpo86+7G0bitKUxmlZSycDZcctfwxS856m454EK4oTWVQ5iUHAYt+HuEMuWCGW38wC7IgbaYHv8+ydOcrIlIBBV8RkQoo+IqIVEDBV0SkAhpwk/mlaHCsYPZXqf2L9i07SayoTeG2QV1R2SxnqEUpIWc76616JWe4hedZfk29snTnKyJSAQVfEZEKKPiKiFRAwVdEpAIacJOu5eNj+P5czvWFcYpJGw7WazscrNd26OV0u8IBt2AwJlof7OihtKxoHbJo1lzZGW5Fg0PR7K1gvbWiIcnZDMQV7RvOfCtbTa0gxWc0sBnVH80uLEhxaUUzJkvQna+ISAUUfEVEKqDgKyJSAQVfaStmtsbMfmhmW83sYTP7SFa+3MzuMLPHs69aHVM6mgbcpN2MUV+d+H4zWwJsNLM7gN8C7nL3m8zsBuAG4GNTHmliAo4cScsC0WCSjRxJC4OBuWmt4TYWrAE3GgwCRuvHQTzwU7TeW14wiFa4v5VcF44p0kLm957GwFzZbcPBzoL0j7NLSVk0WDnzEKo7X2kr7r7T3e/Pvj8IbANOBa4Ejq+AfBtwVTUtFGkOBV9pW2Z2BvA64KfAKnffmb31ArCqYJ/rzGyDmW3YeyRIkC7SJhR8pS2Z2WLgG8D17t6wfITX/9YMHzt191vcfb27rx8anPmy3iKtpuArbcfM+qgH3q+4+zez4l3Hl5DPvu6uqn0izaABN2krVh9p+SKwzd0/Nemt24FrgZuyr98+4cHGxpjYu6+hqGewYObYSNBFcfCltOylF9OyosGxnmAdtJOCdeGCWXPRWm8A1NJjWjjgFwwQFQ1ERe0MUy1OY721QDQwN9s0leG6cNH5QPnUm9GAW9G5TycdaI6Cr7SbNwHvBx4ys01Z2cepB92vm9kHgGeA91bUPpGmUPCVtuLu91I8i/9X57ItIq2kPl8RkQoo+IqIVEDdDtK1fMIZP5KbUVYwwNMTDUaNBDPPooG5/jhNZTQ45seCGW7RrLexgmeUJ4LBpInxYMPoPAuSQkYz3Hqms65cyUGzOVsCrmi9tZKF0eBa0WBl0UBcCbrzFRGpgIKviEgFFHxFRCqgPl/pXg5M5Po5i5b8icrDsmn0hebrLto/Kov2hWncLpU8n6Jt50jREkyznXwxK3NUt+58RUQqoOArIlIBBV8RkQoo+IqIVEADbtK1Ro+M8uz92xvKFg3GS/4MDqYfhaWrV6cb7t2bli1bFjcgeDDfgokX7NuTFIWTMSDOQDYaTMgIs3XF2b5s6VBa2B9kfyvKFhYuuxNtF/zsCyYvxJnJyg2EWdGEiJK8L8oDXbSMUMHPpATd+YqIVEDBV0SkAgq+IiIVUPAVEamABtykrZjZGuDL1FcnduAWd/8bM7sR+F3g+OjUx939u1Mdyx1GRxozfvXW4vuN3t5gQGU0GPQqWwbxYFI0kBaVFR0zWDLIx4+l2wXZtqxwdl80wy6ayVc06FRy23B2X8ExwyV/4k2bL8p0VlT5zBul4CvtZgz4qLvfb2ZLgI1mdkf23qfd/RMVtk2kaRR8pa24+05gZ/b9QTPbBpxabatEmk99vtK2zOwM4HXAT7OiD5vZZjO71cxOLtjnOjPbYGYb9o9HScZF2oOCr7QlM1sMfAO43t0PADcDZwPrqN8ZfzLaz91vcff17r7+5IJJBSLtQN0O0nbMrI964P2Ku38TwN13TXr/88B3TngcoNbbeH9RiwbWgu3qhUHwjgbRioJ8uH/Jsmkc06L9w6VwCo4ZDSaVLau/Mcv9S7Zp7kbcUi1Ix6k7X2krVk/k+kVgm7t/alL55Lm+7wa2zHXbRJpJd77Sbt4EvB94yMw2ZWUfB64xs3XUbzWeBj5YTfNEmkPBV9qKu99L/PfllM/0inQadTuIiFRAd77StQ6OjXPv7gMNZSf1xoNOS4OZb295YFtSdnTH/qRsYMXi8JjWl9Y1EA3cRGkqh4/GxwwG3PzIoWDD4I+HYHYcgA8fTndfEjzJ19cft6l/Ycn6g/2L0lFG5dEgYlBPnBISWpKSsif+mZTadcZ7iojIjCn4iohUQMFXRKQCCr4iIhXQgJt0rZoZi3MDaQsKZlkNBIM5PQPpAFWtP/3I9CyIB7KI0lf2BdsGZRZtB/FAVDQQNo0Bt7C87Ew8KFiHLZrhVnK7om1LD4QVHTNKcznzGWr13TXDTUSkoyj4iohUQMFXRKQCCr4iIhXQgJt0rYPjE9zzcuNMsRUFM9xW9aflFzy+KynbvTudeXbygZHwmLVaOsAzNLQkKRt5Ppg1118wcyoanAsHzILBpYIZauEssWBdOeuPZ475wvScomPawkXpdj0FISha7y0csAsUDQyGA3HBgFkwa61oYM2mkyYzX82M9xQRkRlT8BURqYCCr4hIBRR8RUQqYLOZoSHSzsxsD/BM9nIICHI3Np3qaf+6Wl3P6e6+8kQbKfjKvGBmG9x9veppz3rmsq65PKepqNtBRKQCCr4iIhVQ8JX54hbV09b1zGVdc3lOhdTnKyJSAd35iohUQMFXRKQCCr7S1czsnWb2qJk9YWY3tLiup83sITPbZGYbmnjcW81st5ltmVS23MzuMLPHs6/BWu9NqedGM9uRndMmM3tXE+pZY2Y/NLOtZvawmX0kK2/qOU1RT9PPaUbtU5+vdCszqwGPAW8HtgP3Ade4+9YW1fc0sN7dm/oAv5m9BTgEfNndL8zK/hrY5+43Zf+pnOzuH2tBPTcCh9z9E7M5dq6e1cBqd7/fzJYAG4GrgN+iiec0RT3vpcnnNBO685VudinwhLs/5e6jwNeAKytu07S5+z3AvlzxlcBt2fe3UQ8qrain6dx9p7vfn31/ENgGnEqTz2mKetqCgq90s1OB5ya93k5rP3wOfN/MNprZdS2sB2CVu+/Mvn8BWNXCuj5sZpuzbolZd29MZmZnAK8DfkoLzylXD7TwnMpS8BVpnsvc/fXAFcCHsj/jW87rfYet6j+8GTgbWAfsBD7ZrAOb2WLgG8D17n5g8nvNPKegnpad03Qo+Eo32wGsmfT6tKysJdx9R/Z1N/At6t0erbIr69M83re5uxWVuPsudx939wng8zTpnMysj3pA/Iq7fzMrbvo5RfW06pymS8FXutl9wDlmdqaZ9QNXA7e3oiIzW5QN6mBmi4B3AFum3mtWbgeuzb6/Fvh2Kyo5Hgwz76YJ52T1tXe+CGxz909Nequp51RUTyvOaSb0tIN0tewxos8ANeBWd/8fLarnLOp3u1BfG/GfmlWXmX0VuJx6KsRdwF8A/wJ8HVhLPW3me919VoNlBfVcTv3PcweeBj44qV92pvVcBvwIeAg4vljbx6n3xzbtnKao5xqafE4zap+Cr4jI3FO3g4hIBRR8RUQqoOArIlIBBV8RkQoo+IqIVEDBV0SkAgq+IiIVUPAVEamAgq+ISAUUfEVEKqDgKyJSAQVfEZEKKPiKiFRAwVdEpAIKvk1kZn1m9gdm9jMze9nMjmbref1BlsxbpBJm9mkz88nLwk967yQz25e9/8dVtG8+6q26Ad0iW4TvTuprQ30W+PPsrSuAm6gvX/P1alonwkXAQeA1ZlZz9/FJ7/0JcPzmYPOct2yeUjL1JsiWK/kBcC7wVnd/JPf+euBFd/95Fe0TMbNdwB3AfwZe6+6PZ+WrgCepL+FzDfDqKlZ1mI/U7dAc11JfbuX38oEXwN03KPBKVbIAewrwHep3v+dOevtPgQepL9uzV4F37ij4NscfUl+kryWLGIrM0kXZ183AVuA8ADNbC3yQegC+mPpaZzJHFHxnycxOp35xf6XqtogUuAgYAR4DHiYLvsCNwD3ufne2jfp755CC7+wdv6s44fLTZnazme0wM3W0y1y6mPpfZmPUg++5ZnYu8H7gT83sJGANuTtfM/tQ9gTEWbnyq8zsRTPbZGZbzexBMzt/js6layj4zt6y7OuuEtt+FXh9C9siErmIVwLrw9T7fP8K+D/u/jPgwuy9X9z5mtlS4KPARurBe7J1wGfdfZ27nw/cR32wTqZBwXf2dmdfX32iDd39HncvE6RFmsLMeoDzeSX4bgFOAt4D/FlWdjEwQT0wH3cD8GXqT/FcRKN11IMy2V3zhcCPWtD8rqbgO3s/AQ4Avx29aWaXzW1zRBqcAywkC77ufvx587909+NdZRcBT7r7EQAzOw24GvhEtl9053uTmW0GXgC+4e7fb/WJdBtNspgldz9kZh8DbjazbwP/COyhPtniN4ClwJsqbKLMb8fvWn/Rn+vuvxlsM3mw7X8AN2XX9hbqT0MAv7jTXeLuZ2SvzwE2mNnn3P1AC9rftRR8m8Dd/5eZ7QT+CPhSVvwscA9wa1XtEqEeWPdnd7xFLgS+D2Bm64D3Apeb2ccBA15tZgvd/Sj1u96tx3d098ezAeRl1P8ClJI0w60CZubublW3QyTPzO4A/t7d/2VS2WPANe6+0cyuB85z9w9m770X+HN3vzA+ohRRn+8cMrMvmNn27PvtZvaFqtskcpyZXQEsmxx4M4/ySr/vOuDK7DGz+4D3Ab8+h83sGrrzFRGpgO58RUQq0FbB18zeaWaPmtkTZnZD1e2R9qLrQ7pJ23Q7mFmN+tzztwPbyWbNuPvWon2Ghlb4GWvXlq/EJ4Ky3Pn35P4/Gh/Lt7Tx5UTumGOjuffHyfPDh6Zu10Rjm3x4pPHt0cZjHhjJtxH6rLGdB8Ya6+jJncbzPr7X3VcmB2oTM7o+li72M1ataCycKLjeFyxMy4aPpGW14AGh6LoC6Am2HT+WlvUFefaPBnVD2E4/dDAps/7gmP0L4mPWamlZb7B/dO6QfgbqrUqLeoJ6in52+c9ZkekMW0e/+/znHYLPfMF2kMYPYOODD5X6LLXTo2aXAk+4+1MAZvY14EomPdaSd8batWy49+7CAyb/sYweTTfKfxgGFjW+PvhibofG37aPNH5IfPdzjZsffilt18afNBYMDze+PtrYztHHnm14fejZfQ2vf/DY3qSOVX19Da/vfOlww+vBWuN5/NnIS88kB2kv078+Vq3gp5/6742FoyPhtvaa/CQu8EfuTzdcHnymhoPrCmDpyWnZ/j1p2arT0rJHHoyP+ZoLkiL/tx8mZbY6OOZZ56ZlgC1dnpYtDyZsLiuIJyOH07IogC1cmpaNxb8PLAh2UVkUFItuKI8Np2X9g2lZ8Jkt/I8r+M+0Z+XppT5L7dTtcCowOXJtz8oamNl1ZrbBzDbs2ZsPjNLFpn99vHwo/7ZI22in4FuKu9/i7uvdff3KoRUn3kHmlYbrY9niqpsjUqiduh12UE9rd9xpWdnM5f/8CPpfk36gpP81v0+ukynfxzua+9Mm+pM0161wom6HsQONr48cbfyTbt9Y2m+2qNbY7hfHGl8fmei4/3enfX1M7NnL0c/9Q0PZ+KH4z9xFl69Lyg7duTEpG1ib/oc/9nLcPzuw+qSkbPiZ9K+1wfOSG3h23P1oeMxXv/HMpOz+bz+clJ2xZklStuKXz0rKAHzVqqSs96OfSrfLX+sZW1DuPzkPuiKsP+hrJ+gyrBdGWwYNiq/tqJ0+lnYb2JK0G6YV2ukTeB9wjpmdma30ezX1daVEQNeHdJm2ufN19zEz+zDwr0ANuNXd0//SZV7S9SHdpm2CL4C7fxf4btXtkPak60O6SVsF35lo6BvK9QlZ7jEUj/qXcn26lnu+0fP9RLnnZ6238ZGuRPAYDy/vb3w9ku/zbexDHBwYaHjd/3zj40pv3J/2Y64cajzXnlxfdX/uQd/P7O2+hFQ+4UwcbezTGxsJnrMFGEl/hsPDwfPTR9N+z3wdvygP9h8/ktYzcbRc3fVjpnUdDZ6zHR5JxzcmhuM+257RuDw9QNEzuSW1Yk5BdMxppayqbp5DO/X5iojMGwq+IiIVUPAVEalAZ/f5+kTjlOFc/22+j9eiuem5suT5whM9w9jX2B+b9AFHzwzm6vBjuT633DRIO2V1Y5V7Xmh4fX4t7eSqDTVObV364FONx+zLzbP/7s60nR2uZ/Egg2/OLT+Wf6Y6Y+vekJQNHQv6clemU2wHjhTkYViRPhO8ZO0LSZmtPT0pS5/mrautyy+nBm/YnfbXLzxndVJm5xWs7r4yfc7Xj6X90FYwxdZL5nY44fjI5G0t6LiNyqZx/xg9O5wf44H4eeQwLwVMkZvixHTnKyJSAQVfEZEKKPiKiFRAwVdEpAIdPuDmjfk080ly8klxipJBNxwz14GeS9KRHwhwch3x+QTateAh7oWNg3jWl0/OkxsEXLIs937jYEhtRZrAhZMay/pXNuZStf6CAYQucnjPAX5y852NZeNBciXgjReniWzueTAdhDxt4UBSdmA0nhDxqpPSAapn96eJll67Ok2C88On9yVlAJf96Mmk7LOPpTmC3/xvzyZl6097JDzmilPTPLtL3vU7SZkfKEjhumhZWhYkt/EoT+5AkE8XCpJgRYshBGXhwBzQl/4+/GiQdjQ6nygJPhQPxJWgO18RkQoo+IqIVKCtuh3M7GngIDAOjLn7+mpbJO1E14d0k7YKvpm3unu6KFmkp6dxzbVc/0+SJCdM0JwxTSF5AAASU0lEQVTbJ9eHk3/YPDlCPsH0cON6Vp5Prg74c481FpwgAbs/9/PG93fvbnj58o+3JXUseFVjv9Uzmxsf7u/t7dg/ekpfH4tfvZxf+dj7Gso8n8g+YxencfyKTT9Ltxs6Jd35SMFyRcF6b2fvej495pp0SsV/2bQhPKRdkE6y+MQP7wqOma7hZuem678BYTt9OD0nWxqvHBMlJGciSFK+KB2b8KhvF6BWbkJGOBmjQJjMfXHQpmitx2hBUYj7pkvq2E+giEgna7fg68D3zWyjmV0XbaAFNOe16V0fhwpWFRZpA+0WfC9z99cDVwAfMrO35DfQAprz2vSuj8Xx+mAi7aCt+nzdfUf2dbeZfQu4FLincIfxMTg46e43n1gnnxQnSpKTex4x38drfelznQ1y73vuWUKbCPqZzs31L+bbne9DW3NO4+uDjXf8y4IELnZyY9lrz8s9x9qb+9X/Vfr8aLuZ7vXhR44y/sCDDWUTR+LE4b19aR/j6M82J2X9p6Y/64lDcWKdnpXBYps7didlvQfSxDhH/l+8QtJg8Kzr9nufSspWnp0+U7tgLH4e2U4J+rEv/Y9JUVH/bJQwJxpfiRLw5MdYpitMllPQDxwl1gqTAgXPAxf2LXfDc75mtsjMlhz/HngHsKXaVkm70PUh3aad7nxXAd/K/ofpBf7J3b9XbZOkjej6kK7SNsHX3Z8CLqm6HdKedH1It2mb4DszRuNqebl+mRLPAKa5GuZAvp/qhK9z/VL5/ARRv9WJtinIcdBVeoyewVz/XcEijrYgHZyrDQbPdi5I+wN7CvpSGUjHC2oLg+dX+9N6ehYVPFca1L9gQdrvWBsMxiqC9tTrj5KkRz+naa1MGew+y/3DQ7agTeGinHE94dyBktqmz1dEZD5R8BURqYCCr4hIBTq7z3diAh+Z9IxlLs9C8vxh8Mxuko83n6vhRM/55uWf+4v6ivI5THPPT1pfrj922VDDS6811mGvOjWtY1njApocPti4T/453y40emCY5+5szHtx5Gicl/U142m/+bYfpM8+r35VmuP30KH4mCtOeS4pe+H5NGfCqa9JZ2pu3JjmgAB4/eH0OeW7nkhTXZyzI63n/JcOJ2UAA2vS55Fr70v7kf3owaQMgufpiftiGz6rx8uKPl8l8rAA8XhHT8E9ZZQvIsjjYEGO4XiRULCiukrQna+ISAUUfEVEKqDgKyJSgc7u+BsbxXdP6lcLcudOFs1BT9Zcy+fjDeZ5N+6f639N+oCCtaz6c8+U5vuy8jmG8+vILWjsk5pYm+YGIJenNOmBq3X/Gm5m0FNrPPOenqJ5/+nvKShKjjfdY0bbhnUX5hJIy3uDbcP9oxMqqD/esBXP+RYds+Rzxi14dniu6M5XRKQCCr4iIhWYMvha/u/dJjCzW81st5ltmVS23MzuMLPHs68nT3UMaQ+6PkRm7kR9vi+a2ZvcfauZfQDYDGxx99ksEfAl4O+AL08quwG4y91vMrMbstcfO+GRJsbh8KS8pbm1z1i6vPH1ktxrgFpj31J+zbUoH2/jBvk+pxLxKJ9XNOnzzfV39ebakN8+WIfKBpc27rI0F69mmUc109bXR/+553H6vXeXqjTKVXtRyZ9RugJasWVBWZQf4M3TyCXwm1FO2+D508I8BMG5b7v4dUnZOb9/Zbj70e/dm5SNvpQ+07vs134lKRvb+kR4zNqyNNeGLQyS4w+mz+QyMpKWAcNbn0nKFlyWnucLX7s7KVtx7qrwmCM707zJZZ0o+P4hcHw05zPAIDBhZk9R/6A9lH3dnGWdOiF3v8fMzsgVXwlcnn1/G3A3ZYKvVE3Xh8gMTRl83f0fJr1cCpwFXARcnH19H/DnQI+ZHXL3pelRSlnl7senDr1APXertDldHyIzV/pRM6//zfJk9u9fjpeb2QLgwuzfrLm7m1lhnrZs4cTrANaeojXc2kVbXh9r1jSjSpGWmPWAibsPu/sGd//SLA6zy8xWA2Rf04WuXqnvlQUST5rpjZTMlUqvDy2wKm2sXSZZ3A5cC9yUff12mZ388CF8409eKTiaG+d5eX9uh+CGaWFjUhB/7rGG18lil/lj5JJwJBMogkX7ooX8puL57XN19Kw+K90pN6HEBnPJT5r/oEIrzej6mJ7OeFg/TFhTMqF3UeJxD8oXDATXbX+cBKe2KC3vjRIYRYnoo+TygEWJ36OyIBF90cSLnqCd0f4DQXL6onaGyfFLmvNPoJl9FfgJ8Etmtj0bJb8JeLuZPQ68LXst85CuD5kv5vzO192vKXjrV+e0IdKWdH3IfNFRf3uKiHSLdunznRmfgOFJkyKGc4l1Rhpf+7E0EbX15cryyXnyD6Dn+9eS96dOkjMz+YVBe6Z+HZUlSd7nx/+7SRLsokQsx9KkTEn/ff2Awc4Fx4y2jSZujAf9o73xApo+lm4bJowK902v/3r96USiB3a+nJStefSRcPdnNqWJ3/cfTOt6w1np/nsfSBPOAyxaGiw+ujgqS/uRfThObr99WzpOe+bS9Hf84JP7k7LXHo0nW+3ZM/P5RPPjEygi0mYUfEVEKqDgKyJSgc7u853wxmd788/5Hs0l9wj69RjNJ1NvPEbSx5br800WuzxRkpyZyPdT5l9H/Yj5xCrzsc/XJ9IFUYsWPBwJ+u6iBRej/tmiBDxl+3KDRQC86JhR33T03HiQLKdwsYGgL/j7+9PEOG97OO6fffjFdGHO7SNpH+nZD+1IyrZtDxYCAJYFyf4HB9PzHAyesx0ZDc4d2Ppiuqjo0JbtSdmPD6TXwkTBx3j7aNy/XMY8+ASKiLQfBV8RkQoo+IqIVKCj+3x9eITRx579xeuxXF/NYG4uuJ2yOj3Gksb01v7czxs3WHNObodcn+6yocY68n2p+UToRH10U/fpJnPyc/t7f7DIZ5KwPddp1cELD5ZmPVj0swn44JJ09+D52ah/tTBnQtBnHCY5D/p3raDP1wcWlas/aKcvSPcFYDx9fvb9p6Rp3wfffEm4+5tfSvtIXz6Q9iOveGua2O7STU+Gx+xdliZJD5/zXZT+fidG4n7YFQ+l/bvL/v1FSdlVO9I+7DPPCRZiAPbtTvuR2bI33DZPd74iIhVQ8BURqUAVWc2iBRJvNLMdZrYp+/euuW6XtAddHzJfVNHn+yXSBRIBPu3un5jOgSZGxzn07L5fvD6Sm3/d//yehtd9e15IDzKaW2xvd27+98EXG1+PNz5D6LnnEW1Brq8qmt+fzxmQ5GqYuo83Eb2f7zMsyBXQhr5Ek66PaSm5WGZR/264bdEzxemGpY8Z9dUn+SuK6i6qJzj35UuC62VZtPwnDJ6yOCnr7U37TaP9+1fGCyLMZgHNnoIFNBcOBYtdnpQuPrtiedqP3LcyHRMAOGksfqa4jDm/83X3e4B9J9xQ5iVdHzJftFOf74fNbHP2Z+fJRRuZ2XVmtsHMNrw4doJl3aWbTPv62LP3xaLNRCrXLsH3ZuBsYB2wE/hk0YaT1+ha0dvRT8pJeTO6PrSGm7Sztohe7r7r+Pdm9nngO2X2OzAyxg8ee+WZun1jjX1fb9zf2Pdzfi3tL6utaOzzefnH2xpeL1uR+wDn+tfsVac2vr02N1d9cdqnlKy5doLcu8lzvLk+3lJrwvXF6291gpleHwXHCsuj52rL9qUWHrPkemvhdlFuhqJ2BrkZnHLPGAPhmMGac1em+59+drj7wvPSnA8DLx9M9z8j3b93vKDPdFHwPHOwBhwDQT9wkLcbIHri205P1z8cuuDRpKz3NWvDY/atCJ7pvWdzuG1eW9z5Hl+ZNvNuYEvRtjL/6PqQbjTnd77ZAomXA0Nmth34C+ByM1sHOPA08MG5bpe0B10fMl+0ywKaX5zrdkh70vUh80VbdDuIiMw3bTHgNlN9Zqzqe2VgYVGtsfN+5VBjZ3xtKHhCKfeQ9YJXNT4IbifnBtzyAwTLcsfMDbDZYPAQeT5hSzLglnudHxApOSFgstILSYpQMIgYDEBmG8+8bDrHLNq2rJILG4TnXrBv0WBrGbrzFRGpgIKviEgFFHxFRCrQ0X2+B8YmuPOlVxJ4vJhLctGTS1K+9MGnkmPkE3s8s7kx+c5rz8s9cJ3vdzrc+DB5vifVl6b9zDaYS0SS78PNv873K+WT5AQTKPJ9vKWTvHS5wsTnwUQFC5IRlZ0kAdOYpBFMqChMph4ds+QEmqKJG4yn0/QfuO/5pOxX/t394e57/+/WpCxKpn72Sen4x5EHf56UAfQGiXVqi4Jk6ovLJ1N/aXOaTH158Lt79EdPJ2VnPh8k5QH27UkXGi1Ln0gRkQoo+IqIVEDBV0SkAh3d59tjMDgpWc6Ricb/S/p7cgtR9gULFfY3lvX25hfAzP2Ics/5Wv792gn6byFInn6i1zN4JlfP8UpZwbXSG10++Wv7eHH+MwPUeoIDBFkILUh2BWC9QV214F4xaJPV4r7t5LNdsH9v0KYodgDUwh9UObrzFRGpgIKviEgFFHxFRCpgs5mbXDUz2wM8AwwBQVbjWWn2Mdu9jae7e5pBu4NNuj6gNT//iOpp/7paXU+pz1JHB9/jzGyDu69v52N2Qhu72Vz9rFRP+9fVLp8bdTuIiFRAwVdEpALdEnxv6YBjdkIbu9lc/axUT/vX1Rafm67o8xUR6TTdcucrItJRFHxFRCrQ8cHXzN5pZo+a2RNmdkMTjve0mT1kZpvMbMMMj3Grme02sy2Typab2R1m9nj2NVhQbtrHvNHMdmRt3WRm75pJe7tZs6+PE9Q162un4LhNv56mUU/TrzEzW2NmPzSzrWb2sJl9JCtv6jlNUU97fG7cvWP/ATXgSeAsoB94EDh/lsd8Ghia5THeArwe2DKp7K+BG7LvbwD+ZxOOeSPwR1X/Htr1Xyuuj1ZfO3N1PU2jnqZfY8Bq4PXZ90uAx4Dzm31OU9TTFp+bTr/zvRR4wt2fcvdR4GvAlRW3CXe/B9iXK74SuC37/jbgqiYcU6bWltfHdLXieppGPU3n7jvd/f7s+4PANuBUmnxOU9TTFjo9+J4KPDfp9XZm/8N14PtmttHMrpvlsSZb5e47s+9fAFY16bgfNrPN2Z+Ms/7Ts8u04vqYSquunUirrqdIy64xMzsDeB3wU1p4Trl6oA0+N50efFvhMnd/PXAF8CEze0uzK/D630HNeMbvZuBsYB2wE/hkE44pM9fyayfSxOsp0rJrzMwWA98Arnf3A5Pfa+Y5BfW0xeem04PvDmDNpNenZWUz5u47sq+7gW9R/9O1GXaZ2WqA7Ovu2R7Q3Xe5+7i7TwCfp3lt7RZNvz6m0sJrJ9L06ynSqmvMzPqoB8SvuPs3s+Kmn1NUT7t8bjo9+N4HnGNmZ5pZP3A1cPtMD2Zmi8xsyfHvgXcAW6beq7TbgWuz768Fvj3bAx6/UDPvpnlt7RZNvT6m0uJrJ9L06ynSimvM6ss9fxHY5u6fmvRWU8+pqJ62+dxUPeI323/Au6iPYj4J/Oksj3UW9RHxB4GHZ3o84KvU/5w5Rr2f8QPACuAu4HHgTmB5E475j8BDwGbqF+7qqn8f7favmdfHXFw7c3U9VXmNAZdR71LYDGzK/r2r2ec0RT1t8bnR9GIRkQp0ereDiEhHUvAVEamAgq+ISAUUfEVEKqDgKyJSAQVfEZEKKPi2iJldbWYj2QwbEZEGCr6tcwmw1d2PVd0Qkekys/dkNw+Lqm5Lt1LwbZ1LgAeqboTIDP0U+GV3P1x1Q7qVgm/rrKM+nREAM/sPZvaimX3GzGoVtkvkhNx9h7tvOvGWMlMKvi1gZiupZ9HflL3+feA71LP0X+/u41W2T+REzGyPmV1fdTu6WW/VDehSl1BP6LHVzD4PvAd4l7vfXWmrREows1cDQ0z6y02aT8G3NdZRX47lm8BK6n1nT1TbJJHSLsm+PlhpK7qcgm9rXAIY8GbgbQq80mEuAZ5z9/1VN6SbKfi2xiXA3wKvAf7BzN7g7nsqbpNIWReju96W04Bbk5nZAHAu9UTN1wF7gf9tZvqPTjrFJSj4tpyCb/OdD/QBD7n7UerLlFwAfLrSVomUkN08/BIKvi2n4Nt8lwCHqS9bg7s/A/wm8N/M7LcqbJdIGRcANRR8W07LCInIL5jZbwN/Byzx+uq+0iK68xWRyd4E/EyBt/UUfEUEMzvVzK4GrqY+G1NaTN0OIoKZfQ74T8DXgevdfbTiJnU9BV8RkQqo20FEpAIKviIiFVDwFRGpgIKviEgFFHxFRCqg4CsiUgEFXxGRCvx/XTXmOvJULLsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%%\n",
    "cmap = 'Reds'\n",
    "pl.close(10)\n",
    "pl.figure(10, (5, 5))\n",
    "fs = 15\n",
    "l_x = [0, 5, 10, 15]\n",
    "l_y = [0, 5, 10, 15, 20, 25]\n",
    "gs = pl.GridSpec(5, 5)\n",
    "\n",
    "ax1 = pl.subplot(gs[3:, :2])\n",
    "\n",
    "pl.imshow(C1, cmap=cmap, interpolation='nearest')\n",
    "pl.title(\"$C_1$\", fontsize=fs)\n",
    "pl.xlabel(\"$k$\", fontsize=fs)\n",
    "pl.ylabel(\"$i$\", fontsize=fs)\n",
    "pl.xticks(l_x)\n",
    "pl.yticks(l_x)\n",
    "\n",
    "ax2 = pl.subplot(gs[:3, 2:])\n",
    "\n",
    "pl.imshow(C2, cmap=cmap, interpolation='nearest')\n",
    "pl.title(\"$C_2$\", fontsize=fs)\n",
    "pl.ylabel(\"$l$\", fontsize=fs)\n",
    "#pl.ylabel(\"$l$\",fontsize=fs)\n",
    "pl.xticks(())\n",
    "pl.yticks(l_y)\n",
    "ax2.set_aspect('auto')\n",
    "\n",
    "ax3 = pl.subplot(gs[3:, 2:], sharex=ax2, sharey=ax1)\n",
    "pl.imshow(M, cmap=cmap, interpolation='nearest')\n",
    "pl.yticks(l_x)\n",
    "pl.xticks(l_y)\n",
    "pl.ylabel(\"$i$\", fontsize=fs)\n",
    "pl.title(\"$M_{AB}$\", fontsize=fs)\n",
    "pl.xlabel(\"$j$\", fontsize=fs)\n",
    "pl.tight_layout()\n",
    "ax3.set_aspect('auto')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute FGW/GW\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Relative loss|Absolute loss\n",
      "------------------------------------------------\n",
      "    0|4.734462e+01|0.000000e+00|0.000000e+00\n",
      "    1|2.508258e+01|8.875498e-01|2.226204e+01\n",
      "    2|2.189329e+01|1.456747e-01|3.189297e+00\n",
      "    3|2.189329e+01|0.000000e+00|0.000000e+00\n",
      "Elapsed time : 0.005539894104003906 s\n",
      "It.  |Loss        |Relative loss|Absolute loss\n",
      "------------------------------------------------\n",
      "    0|4.683978e+04|0.000000e+00|0.000000e+00\n",
      "    1|3.860061e+04|2.134468e-01|8.239175e+03\n",
      "    2|2.182948e+04|7.682787e-01|1.677113e+04\n",
      "    3|2.182948e+04|0.000000e+00|0.000000e+00\n"
     ]
    }
   ],
   "source": [
    "#%% Computing FGW and GW\n",
    "alpha = 1e-3\n",
    "\n",
    "ot.tic()\n",
    "Gwg, logw = fused_gromov_wasserstein(M, C1, C2, p, q, loss_fun='square_loss', alpha=alpha, verbose=True, log=True)\n",
    "ot.toc()\n",
    "\n",
    "#%reload_ext WGW\n",
    "Gg, log = gromov_wasserstein(C1, C2, p, q, loss_fun='square_loss', verbose=True, log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize transport matrices\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6oAAADxCAYAAADC3Uq3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu8XGV56PHfI+ClSAGbiFwSsK3lBKgEjdgW5OClFDgk2NZaaKtS6IkWaMsp1FsPSPFzrG1Jq1a80EKjFvFaaLDYSi8WtIoEBOQiBRFMAiWJyEUBLfKcP9baMox7Zs3sPXvNu2f/vp/P/mTmXe9+59krM8+8z6w1643MRJIkSZKkUjxp3AFIkiRJktTJQlWSJEmSVBQLVUmSJElSUSxUJUmSJElFsVCVJEmSJBXFQlWSJEmSVBQLVUmSJElSUSxUJUmSJA0tIv44Ik4ZdxzjEBFfioh9xx3HJLNQ1USLiBsj4tBZ/P5YE3BE3BERLxuwrwlTmqfGnWvaMGg+M5dJvdWvo4cj4tsdP7vV246JiCsj4jsRsbm+fWJERL39TRHx6a7xbu3RdswAsSwGXg28v6Nt54jIiLizq+8eEfFQRNw787/+CePtEBFvi4jbIuLBiPh6RLy7jqmt8c4GzprZX6BBWKguAKNOTOM0TOEGkJn7ZuZnZ/hYT0jAbSXfWTBhakFpmpTNF9NN9ur2EiZi42Auk/pbmZlP7/i5KyJOBd4J/BnwLGAX4HXAQcCT69+7HPi5iNgGICJ2BbYDDuhq+8m6b5PjgEsz8+GOtuXAJmDniNiho/1twEbgun4DRsSZEXFmQ5+dgCuA/wEckZk7AC+q/5Y9B4h7VOOtA14cEc8a9jE1GAvVhWHUiWnORMS243z8LsfxxAQ8q+TbAhOmFowBJ2Wd/UvKLd2Oo2uyV9BEbBzMZdIQImJHqg93TszMT2Tmg1n5cmb+emZ+t+56FdVrfnl9/0XAvwG3dLV9LTPvGuChjwD+vattOfBl4EZg3zq+5wE/B3y+3jZbfwHcC7wiM28FyMyNmfnazFzf1niZ+QhwNfALM3hMDcBCdWEYKjFFxBsj4mv1J+43RcQvdg4WEW+IiE319lsi4qUN7btFxCcjYkv9Cf7vdo13R/271wPfiYhtpxsrIj4ELAUuiepUl9cPOPbLOm6fFhHXR8T9EfHRiHhqn/3WnYBnlHwjYllEfDYi7ovqVORV08TYGFdE/EFEfLKr7V0R8U4wYWrhGHRS1iO39Hw91v3/oH4tficizouIXSLi03Uu+ueI2Lmj/7Rj1Y/5ia6Y3xkR7+rxJ0032StiItZtgP1nLpPa97PAU4C/79cpM78HXAkcUjcdQvUB1ue62gY9aPHTVHPJTgcA1wLXA/vVbWuANwH71NtmLCKWAK8C/jAzH5vNWCMa72Zg/9nGoelZqC4AM0hMX6MqXHcE/gj42/qIKxGxN3Ay8IL6E/lfAO7o0/4k4BKqo427Ay8FTomI7gnIscD/AnYCfmK6sTLzVcA3qE95oTo9bJCxO70SOBx4NvBcqiMZvXQn4KGTb0RsV8f4GeCZwO8AF9T7a9i4/hY4vD4yMnWE6Bjggx19TJhaCAaalNU6c0vQ/Hr8ZeDngZ8CVgKfBt4MLKZ6z/xdaHxtfwQ4MuozL6I6c+WVwId7xPiEXFPgRGxqnEHymblMmnsX1x8W3RcRFwOLgK2Z+ehUh4j4j3r7wxFxSMfv/juPz/1eRDUfvKKrrfuDs152Ah7salvO43OlfesPs54CXESV62Z7RPVlwJbM/MJ0GyNix6i+5/7tiNhvuj7DjFePeWBEfCEiLo+IC+tcOOVBqv2gOWChunAMnJgy8+OZeVdmPpaZHwVuBQ6s+32fKuHsExHbZeYdmfm1Pu0vABZn5lmZ+b3MvB34K6pJSad3ZeaG+tS3XmN1G3Ts7se5KzPvpZpwLe/TtzsBzyT5/gzwdODtdYz/CnyKavI8VFyZeTfVhwm/UjcdTvXGdHVHNxOmFoJhJmWduWWQ1+NfZuY9mbmJKkdeWR+pfYTqtX5A3a/nWJl5J3ANMHU2ykuAhzLziz3+nu5c0+pErGES1mmQ/Wcuk+beyzNzp/rn5cA3gUXR8RWHzPy5zNyp3tY5378cODginkE1h7oV+A+qr4g9g+qD+EGPqH4L+MFXoSLiKcAyHp8rLQf+BPh9qg//nkT1IdQTRMSnpgpv4I3AGzsK8U91dd+F6qBFLw9RfTj5iT59hhkPYAPwksw8BLgDOLpj2w7AfQM+loZkobpwDJyYIuLVEXFtR9LYj2piSGbeBpwCnAlsjoiPRMRuvdqpvvu0W0fCuY/q6MQuXfFtmLrRZ6xug47d6b86bj9ENenq5QcJeKbJF9gN2NB1FONOqiPAM4nrA8Bv1Ld/A/hQ13YTphaCYSZlGzpuD/J6vKfj9sPT3J96bTaN9WEeL+B+jd5HU6Frskf7E7F+k7BOg+w/c5nUvi8A36X3a7e7747A/6b62hKZ+QBwV912V2Z+fcDHvZ5qDjRlP6rX/e3AV4BDgevqD+kOAG7o/IBxSmYeNVV4A2+n+jBsqhA/qqv7N4Dd6zP2fkhm/ndmbhkw/sbx6jHvzsevIfA9oDMHLmP81yiZWBaqC8dAiSki9qQ6Knky8GN10riB6pQ56t/7cGYeTFUoJlXB1qt9A/D1joSzU2bukJlHdsWXT7jT4zG6+g069kx1JuAZJV+q/bukKwEupboo00xcDDy3PopyFHBB13YTphaCYSZlnTljlK/HprE+DhwaEXtQHVntV6h2T/ZanYg1TMI6jXL/mcukEcnM+6i+qvWeiHhFVFf4flJELAe27+r7MLCe6oP2Kzo2fa5uG+aimpcC/7Pj/gHA9Vm5D3gx1XwSHj8rbbamjrC+PSJ+FCAifqq+DsCifr8YEWsjYu1Mx6vnyIdRnS1C/R385wOXzeYPUm8WqgvEEIlpe6qJ3RaAiPhNHv8+JhGxd0S8pD7C+AjVEYbHerUDXwIejOriIk+LiG0iYr+IeEGvWPuMBdXRjR+vbw899pA6E/BMk++VVAXu6yNiu6jWdF1J9R22odWnH36CatL7pcz8wVESE6YWimEmZV1G+XrsO1ZdSH4W+BuqD9SmO+NiSvdkbywTse5J2DRGtv/MZdJoZeafUs3nXk81V7qHasmrN1CdQdfp36m+Z/65jrYr6rZhCtUPUn0f/2n1/SfMhzLzs5m5tb57ACMoVOuDLC+h+nDv1vpsuouA73Q8Vi9LqA/WDDtenTs/BByXmf9dN68EPpuDXSFZM2ChurA0JqbMvInqAkFfoEpyP80TX9RPoTotYyvVKV7PpLqY0LTtmfl9qk/LlwNfr7f/NdXR3V56PQbAHwP/t04k/2cGYw+jMwHPKPlmdSGrlVRX9dwKvAd4dWZ+dRZxfYDq/6X7VDkTphaMISdlU78zstfjgGN9mOr7of2OpkLXZG8cE7Eek7CZ/M3DMJdJQ8rMvTLzn3tsuyAzD8zMH8nMxZn5wsw8t37tdvZ7U2ZGZl7T0faxuu39Pzxyz1i2UuWv19b3T87M3+nR92WZ+e4BxjwzM89s6POfmfnyzNylPptu38x8c7/fiYgnU319Ye2w49VfM/kI8EeZ2XmRzdOAM5r+Js1cZGZzL2mBioi3AZsz8x3jjmVKRCwFvgo8q56ATrVfCZyQmTeMLThJMzKTXFMfOT176jVfT8SuA57bq9jsMc62VGuXrsnMfxkq8Fkwl0maCxFxKdUBhjuB92fm2lmO9yrgHVRf/QJ4b1YXG9Ucs1CV5pH6u2F/DvxoZh4/7ngkjccoJ2LjmISZyyRJTSxUpXkiIranOr3xTuDwzNzQ8CuSVBxzmSRpEBaqkiRJkqSieDElSZIkSVJRtm3uMv8tWrQo99xzr3GHoUJ9+eZv9N1+wLKlLUUyOTbc/0jf7Ut2fOqcx3DnnXewdevWaO45vzTls6Z9D7D1rs19t/ucl8qyUPOZpHa1NYe45pqrt2bm4qZ+RRWqEXE48E5gG+CvM/PtXdufQnUZ7OcD3wR+NTPvaBp3zz334vNXrh99wJoIO7/g5L7bP39l49XU1eW0S/otGQlnr1w25zEc9MIVc/4YTeYipzXls6Z9D3DeWef03e5zXirLQs1nktrV1hziadvFnYPEU8ypvxGxDXAO1fps+wDHRsQ+Xd1OAL6VmT8J/AXwJ+1GKUmDMadJmhTmM0njUEyhChwI3JaZt9cLE38EOLqrz9FUC4QDfAJ4aURM3KkwkiaCOU3SpDCfSWpdSYXq7kDnJeo31m3T9snMR4H7gR+bbrCIWB0R6yNi/ZatW+YgXEnqa2Q5zXwmaczMZ5JaV1KhOlKZeW5mrsjMFYsXNX5XV5KKZT6TNCnMZ5IGVVKhuglY0nF/j7pt2j4RsS2wI9UX9iWpNOY0SZPCfCapdSUVqlcBz4mIZ0fEk4FjgHVdfdYBr6lvvwL418zMFmOUpEGZ0yRNCvOZpNYVszxNZj4aEScD/0R16fPzM/PGiDgLWJ+Z64DzgA9FxG3AvVSJUpqVb11VxlIcp667qe/2Nau6L7BYrlEsP9O0bFAp/2+9zFVO+/LN3+i7b44//cTG2Erfd5LK4hxNWhgG+WypaQ7RNH8bRjGFKkBmXgpc2tV2RsftR4BfaTsuSZoJc5qkSWE+k9S2kk79lSRJkiTJQlWSJEmSVBYLVUmSJElSUSxUJUmSJElFsVCVJEmSJBXFQlWSJEmSVBQLVUmSJElSUYpaR1Uq0WmX3NzY5+yVy2b9OGtW7TPrMdowyELOTYtBD2IUY0yiA5Yt5fNXjn/fND0Pjj/9xMYx5stzXpKk0g0yP2t6bx7F+/Ig87enbXfOQGN5RFWSJEmSVBQLVUmSJElSUSxUJUmSJElFsVCVJEmSJBXFQlWSJEmSVJRiCtWIWBIR/xYRN0XEjRHxe9P0OTQi7o+Ia+ufM8YRqyQ1MadJmhTmM0njUNLyNI8Cp2bmNRGxA3B1RFyWmTd19bsiM48aQ3ySNAxzmqRJYT6T1LpiCtXMvBu4u779YETcDOwOdCdBqVWjWCO1FIOssXXCGSf13e76poOZzzmtrbVyR+HUdf13p2u1SrM3n/OZVIKm96rz3/qexjGa3ndLeV8epWJO/e0UEXsBBwBXTrP5ZyPiuoj4dETs22eM1RGxPiLWb9m6ZY4ilaRms81p5jNJpTCfSWpLcYVqRDwd+CRwSmY+0LX5GmDPzNwf+Evg4l7jZOa5mbkiM1csXrR47gKWpD5GkdPMZ5JKYD6T1KaiCtWI2I4qAV6QmX/XvT0zH8jMb9e3LwW2i4hFLYcpSQMxp0maFOYzSW0rplCNiADOA27OzD/v0edZdT8i4kCq+L/ZXpSSNBhzmqRJYT6TNA7FXEwJOAh4FfCViLi2bnszsBQgM98HvAL47Yh4FHgYOCYzcxzBSlIDc5qkSWE+k9S6YgrVzPwcEA193g1M3iWtJE0cc5qkSWE+kzQOxZz6K0mSJEkSWKhKkiRJkgpTzKm/0lzY+QUnN/aZxAWSeynlb21a+Bpgzap9WohE0ynleTLI6/eEM05qIRJJkqY3iveqUt53S+MRVUmSJElSUSxUJUmSJElFsVCVJEmSJBXFQlWSJEmSVBQLVUmSJElSUSxUJUmSJElFsVCVJEmSJBXFdVQ10UpZl2qhrefatE6qa6ROvtMuubmxz3lnndN3+yS9JiRJk6mt96qmueQgcTS9N5+9ctlQMc01j6hKkiRJkopSXKEaEXdExFci4tqIWD/N9oiId0XEbRFxfUQ8bxxxSlIT85mkSWE+k9S2Uk/9fXFmbu2x7QjgOfXPC4H31v9KUonMZ5ImhflMUmuKO6I6gKOBD2bli8BOEbHruIOSpBkwn0maFOYzSSNVYqGawGci4uqIWD3N9t2BDR33N9ZtTxARqyNifUSs37J1yxyFKkl9mc8kTQrzmaRWlVioHpyZz6M6heSkiDhkJoNk5rmZuSIzVyxetHi0EUrSYMxnkiaF+UxSq4orVDNzU/3vZuAi4MCuLpuAJR3396jbJKko5jNJk8J8JqltRRWqEbF9ROwwdRs4DLihq9s64NX11eV+Brg/M+9uOVRJ6st8JmlSmM8kjUNpV/3dBbgoIqCK7cOZ+Y8R8TqAzHwfcClwJHAb8BDwm2OKVRrYKBZhPu+sc0byOG1Ys2qfcYdQggWdzwZZNPzslWU8XyfptSfNkQWdzzR/7fyCk/tuP+GMkxrHGOT9rA2jeJ8p5W8ZVFGFambeDuw/Tfv7Om4n0PyskqQxMp9JmhTmM0njUNSpv5IkSZIkWahKkiRJkopioSpJkiRJKoqFqiRJkiSpKBaqkiRJkqSiWKhKkiRJkopS1PI00kLWtLbVIGtONq0XNsgaXKeuu6nvdtdI1XzS9JoAOP70E/tud41USZqfmvJ305wHRjO30sx4RFWSJEmSVBQLVUmSJElSUSxUJUmSJElFsVCVJEmSJBXFQlWSJEmSVJRiCtWI2Dsiru34eSAiTunqc2hE3N/R54xxxStJ/ZjTJE0K85mkcShmeZrMvAVYDhAR2wCbgIum6XpFZh7VZmySNCxzmqRJYT6TNA7FHFHt8lLga5l557gDkaQRMKdJmhTmM0mtKOaIapdjgAt7bPvZiLgOuAs4LTNvnK5TRKwGVgMsWbp0ToKUppSyGPQoHmfNqn1GEEl/p11yc2Ofs1cum/M4WjSrnGY+m7m2XntNOeD4009sHKON1540AuYzteLUdTf13X7+W9/TOEbTe8AgeXfNqnbeR/TDijuiGhFPBlYBH59m8zXAnpm5P/CXwMW9xsnMczNzRWauWLxo8dwEK0kNRpHTzGeSSmA+k9Sm4gpV4Ajgmsy8p3tDZj6Qmd+ub18KbBcRi9oOUJKGYE6TNCnMZ5JaU2Kheiw9TimJiGdFRNS3D6SK/5stxiZJwzKnSZoU5jNJrSnqO6oRsT3w88BrO9peB5CZ7wNeAfx2RDwKPAwck5k5jlglqYk5TdKkMJ9JaltRhWpmfgf4sa6293XcfjfgN5olzQvmNEmTwnwmqW0lnvorSZIkSVrALFQlSZIkSUUp6tRfqURN63hBe2s1TopRrJHatG7ld2/5xqwfQxpUUw4YZO3gUtZjlqS51pTvAE4446S+282Jk88jqpIkSZKkolioSpIkSZKKYqEqSZIkSSqKhaokSZIkqSgWqpIkSZKkolioSpIkSZKKYqEqSZIkSSpK30I1IixkJUmSCuMcTdKk27Zh+zcj4qDMvCkiTgCuB27IzIdbiG1kNtz/CKeuu6nn9jWr9mklDhdzn5/aen5oOE2vl4Ne+MWWIpGanb1y2QB9+j+nm95DAE4446RZx6F5YyLmaJp/+s2pp5z/1vf03e6cV4No+jTu94EH6tvvAL4IPBARt0TExyPijIh4eUT8+DAPGhHnR8TmiLiho+0ZEXFZRNxa/7tzj999Td3n1oh4zTCPK0mjZj6TNCYjn6OZzySVpG+hmpl/k5kb67s/CvwU8CvABXXbrwGfAG6LiAemGaKXtcDhXW1vBP4lM58D/Et9/wki4hnAW4AXAgcCb+mVMCWpJWsxn0lq2RzN0dZiPpNUiKZTf38gMxP4Wv1z8VR7RDwV2K/+GXSsyyNir67mo4FD69sfAD4LvKGrzy8Al2XmvfVjX0aVUC8c9LElaZTMZ5LGbVRzNPOZpJIMXKj2kpmPAOvrn9nYJTPvrm//F7DLNH12BzZ03N9Yt/2QiFgNrAZ4+qJdZxmaJA1lzvLZkqVLRximpEk2ojma+UzSWBR5xbj6k8Gc5RjnZuaKzFzxtB2fMaLIJGk4o85nixctHlFkkjQc85mkNpVUqN4TEbsC1P9unqbPJmBJx/096jZJKon5TNKkMJ9JGouSCtV1wNRV4l4D/P00ff4JOCwidq6/pH9Y3SZJJTGfSZoU5jNJYzHr76jORERcSPXF/EURsZHqSnFvBz5WrwV2J/DKuu8K4HWZ+VuZeW9EvBW4qh7qrKkv7vezZMenFrEW5iStGdW0htYo9vcgawZO0j7V/NR2PtPCNUi+a8rN5lX1Yz4TNOeJpvWawTyi0RhLoZqZx/bY9NJp+q4Hfqvj/vnA+XMUmiQNxXwmaVKYzySVpKRTfyVJkiRJslCVJEmSJJXFQlWSJEmSVBQLVUmSJElSUSxUJUmSJElFsVCVJEmSJBXFQlWSJEmSVJSxrKOqsjUtGA+wZtU+s36cpgWlXSxakobTlJvXrGrOq025+YQzTmoc4+yVyxr7SBq90y65ue/28846p3EM518qhUdUJUmSJElFsVCVJEmSJBXFQlWSJEmSVBQLVUmSJElSUVovVCPi/IjYHBE3dLT9WUR8NSKuj4iLImKnHr97R0R8JSKujYj17UUtSdMzp0maFOYzSSUZxxHVtcDhXW2XAftl5nOB/wTe1Of3X5yZyzNzxRzFJ0nDWIs5TdJkWIv5TFIhWi9UM/Ny4N6uts9k5qP13S8Ce7QdlyTNhDlN0qQwn0kqSYnrqB4PfLTHtgQ+ExEJvD8zz+01SESsBlYDLFm6dORBTrKImPUYTevwget0acGYdU4zn6lNTbl5kLW2XSd7YpnPxmiQudXxp5/Yd7uvPc0nRV1MKSL+EHgUuKBHl4Mz83nAEcBJEXFIr7Ey89zMXJGZKxYvWjwH0UpSf6PKaeYzSeNmPpPUtmIK1Yg4DjgK+PXMzOn6ZOam+t/NwEXAga0FKElDMKdJmhTmM0njUEShGhGHA68HVmXmQz36bB8RO0zdBg4DbpiurySNkzlN0qQwn0kal3EsT3Mh8AVg74jYGBEnAO8GdgAuqy9r/r66724RcWn9q7sAn4uI64AvAf+Qmf/YdvyS1MmcJmlSmM8klaT1iyll5rHTNJ/Xo+9dwJH17duB/ecwNEkamjlN0qQwn0kqSRGn/kqSJEmSNMVCVZIkSZJUFAtVSZIkSVJRWv+OaokGWbx8zap9Wohk9gZZDLppseezVy6bdRwuKC1Jk2mQ98M1q/q/BwzyXnXCGSf13T6K9yqpLaN4zju30kLjEVVJkiRJUlEsVCVJkiRJRbFQlSRJkiQVxUJVkiRJklQUC1VJkiRJUlEsVCVJkiRJRbFQlSRJkiQVxXVUmT9rpELzOlxtrbF12iU3993u+naSpF4Gea9qWuN8FOuGa/KV8jwZxXNeWmg8oipJkiRJKkrrhWpEnB8RmyPiho62MyNiU0RcW/8c2eN3D4+IWyLitoh4Y3tRS9L0zGmSJoX5TFJJxnFEdS1w+DTtf5GZy+ufS7s3RsQ2wDnAEcA+wLERMX/O2ZU0qdZiTpM0GdZiPpNUiNYL1cy8HLh3Br96IHBbZt6emd8DPgIcPdLgJGlI5jRJk8J8JqkkJX1H9eSIuL4+7WTnabbvDmzouL+xbptWRKyOiPURsX7L1i2jjlWSmowsp5nPJI2Z+UxS60opVN8L/ASwHLgbWDPbATPz3MxckZkrFi9aPNvhJGkYI81p5jNJY2Q+kzQWRRSqmXlPZn4/Mx8D/orqFJJum4AlHff3qNskqSjmNEmTwnwmaVyKKFQjYteOu78I3DBNt6uA50TEsyPiycAxwLo24pOkYZjTJE0K85mkcdm27QeMiAuBQ4FFEbEReAtwaEQsBxK4A3ht3Xc34K8z88jMfDQiTgb+CdgGOD8zb2wr7qZFmNesaufidqUsXn72ymXjDgFoXsi7lP2lyTVfc5pUuqb31TWrZp/fm95DAI4//cSGOCbn4raTmM8GmQc0PQ+angMwmufBJD2XpFFovVDNzGOnaT6vR9+7gCM77l8K/NBl0SVpXMxpkiaF+UxSSYo49VeSJEmSpCkWqpIkSZKkolioSpIkSZKKYqEqSZIkSSqKhaokSZIkqSgWqpIkSZKkorS+PM18FRHjDmFkJmnt0fkUqySpXaN4vzvtkptHFY4K1fQ8GOQ5MElzK6kUHlGVJEmSJBXFQlWSJEmSVBQLVUmSJElSUSxUJUmSJElFsVCVJEmSJBWl9av+RsT5wFHA5szcr277KLB33WUn4L7MXD7N794BPAh8H3g0M1e0ErQk9WBOkzQpzGeSSjKO5WnWAu8GPjjVkJm/OnU7ItYA9/f5/Rdn5tY5i06ShrMWc5qkybAW85mkQrReqGbm5RGx13Tbolqs9JXAS9qMSZJmypwmaVKYzySVZBxHVPt5EXBPZt7aY3sCn4mIBN6fmef2GigiVgOrAZYsXdr3QQdZyPnslcv6bm9a6BnKWey5lDikBWAkOW2YfCbpcaN4v2t6/z913U2NY6xZtU/f7U1ziO/e8o3Gx2jByPMZ2z29799+/OknNgbVtG9Hoek5UPVxbiWNWmmF6rHAhX22H5yZmyLimcBlEfHVzLx8uo51gjwX4PnPX5GjD1WSGo0kp5nPJBVg5PnsST/yTPOZpJ6KuepvRGwL/BLw0V59MnNT/e9m4CLgwHaik6ThmNMkTQrzmaRxKKZQBV4GfDUzN063MSK2j4gdpm4DhwE3tBifJA3DnCZpUpjPJLWu9UI1Ii4EvgDsHREbI+KEetMxdJ1SEhG7RcSl9d1dgM9FxHXAl4B/yMx/bCtuSZqOOU3SpDCfSSrJOK76e2yP9uOmabsLOLK+fTuw/5wGJ0lDMqdJmhTmM0klKenUX0mSJEmSLFQlSZIkSWUpbXmaObHh/kf6rpU6yPpYTUaxVtt8WotVkiRBRDT2aXp/b3pvP+iFXxwqpvnigGVL+fyVvf/2Qda5n+2+lVQuj6hKkiRJkopioSpJkiRJKoqFqiRJkiSpKBaqkiRJkqSiWKhKkiRJkopioSpJkiRJKoqFqiRJkiSpKBaqkiRJkqSiRGaOO4Y5FxFbgDvHHYekVu2ZmYvHHcSomc+kBWkh5bNFwNYxhDMs4xy9+RLrfIkTyo11oJy2IApVSZIklS8i1mfminHH0cQ4R2++xDpf4oT5Fet0PPVXkiRJklQUC1VJkiRJUlEsVCVJklSKc8cdwICMc/TmS6zzJU6YX7H+EL+jKkmSJEkqikdUJUmSJElFsVCVJEmSJBXFQlWSJEljFRGHR8QtEXFbRLxx3PH0ExF3RMRXIuLaiFg/7niYuI1vAAAEC0lEQVSmRMT5EbE5Im7oaHtGRFwWEbfW/+48zhin9Ij1zIjYVO/XayPiyHHGWMe0JCL+LSJuiogbI+L36vai9mufOIvbp8PwO6qSJEkam4jYBvhP4OeBjcBVwLGZedNYA+shIu4AVmTm1nHH0ikiDgG+DXwwM/er2/4UuDcz315/ALBzZr5hnHHWcU0X65nAtzPz7HHG1ikidgV2zcxrImIH4Grg5cBxFLRf+8T5Sgrbp8PwiKokSZLG6UDgtsy8PTO/B3wEOHrMMc07mXk5cG9X89HAB+rbH6AqXsauR6zFycy7M/Oa+vaDwM3A7hS2X/vEOa9ZqEqSJGmcdgc2dNzfSNmT7AQ+ExFXR8TqcQfTYJfMvLu+/V/ALuMMZgAnR8T19anBRZymPCUi9gIOAK6k4P3aFScUvE+bWKhKkiRJgzs4M58HHAGcVJ/GWrysvu9X8nf+3gv8BLAcuBtYM95wHhcRTwc+CZySmQ90bitpv04TZ7H7dBAWqpIkSRqnTcCSjvt71G1FysxN9b+bgYuoTl0u1T319xenvse4eczx9JSZ92Tm9zPzMeCvKGS/RsR2VMXfBZn5d3Vzcft1ujhL3aeDslCVJEnSOF0FPCcinh0RTwaOAdaNOaZpRcT29cVqiIjtgcOAG/r/1litA15T334N8PdjjKWvqcKv9osUsF8jIoDzgJsz8887NhW1X3vFWeI+HYZX/ZUkSdJY1ctmvAPYBjg/M//fmEOaVkT8ONVRVIBtgQ+XEmtEXAgcCiwC7gHeAlwMfAxYCtwJvDIzx34Rox6xHkp1imoCdwCv7fge6FhExMHAFcBXgMfq5jdTff+zmP3aJ85jKWyfDsNCVZIkSZJUFE/9lSRJkiQVxUJVkiRJklQUC1VJkiRJUlEsVCVJkiRJRbFQlSRJkiQVxUJVkiRJklQUC1VJkiRJUlEsVCVJkiS1KiJ+KSK+GxHbjzsWlSkyc9wxSJIkSVpAImJ3YHFmXjvuWFQmC1VJkiRJUlE89VeSJElSqyJiS0ScMu44VC4LVUmSJEmtiYjdgEWAp/2qJwtVSZIkSW3av/73urFGoaJZqEqSJElq0/7Ahsz81rgDUbksVCVJkiS16bl4NFUNLFQlSZIktWl/LFTVwEJVkiRJUisi4inA3lioqoGFqiRJkqS27Atsg4WqGlioSpIkSWrL/sBDwG3jDkRls1CVJEmS1JaDgC9l5mPjDkRls1CVJEmSNKciYveIOAY4BvjUuONR+SIzxx2DJEmSpAkWEe8Hfhn4GHBKZn5vzCGpcBaqkiRJkqSieOqvJEmSJKkoFqqSJEmSpKJYqEqSJEmSimKhKkmSJEkqioWqJEmSJKkoFqqSJEmSpKJYqEqSJEmSivL/ASQ5WP1N4VHtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 936x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% visu OT matrix\n",
    "cmap = 'Blues'\n",
    "fs = 15\n",
    "pl.figure(2, (13, 5))\n",
    "pl.clf()\n",
    "pl.subplot(1, 3, 1)\n",
    "pl.imshow(Got, cmap=cmap, interpolation='nearest')\n",
    "#pl.xlabel(\"$y$\",fontsize=fs)\n",
    "pl.ylabel(\"$i$\", fontsize=fs)\n",
    "pl.xticks(())\n",
    "\n",
    "pl.title('Wasserstein ($M$ only)')\n",
    "\n",
    "pl.subplot(1, 3, 2)\n",
    "pl.imshow(Gg, cmap=cmap, interpolation='nearest')\n",
    "pl.title('Gromov ($C_1,C_2$ only)')\n",
    "pl.xticks(())\n",
    "pl.subplot(1, 3, 3)\n",
    "pl.imshow(Gwg, cmap=cmap, interpolation='nearest')\n",
    "pl.title('FGW  ($M+C_1,C_2$)')\n",
    "\n",
    "pl.xlabel(\"$j$\", fontsize=fs)\n",
    "pl.ylabel(\"$i$\", fontsize=fs)\n",
    "\n",
    "pl.tight_layout()\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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}