summaryrefslogtreecommitdiff
path: root/notebooks/plot_gromov.ipynb
blob: 6a237e62f0d034d3879b5cbf373b972515a7abf0 (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
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Gromov-Wasserstein example\n",
    "\n",
    "\n",
    "This example is designed to show how to use the Gromov-Wassertsein distance\n",
    "computation in POT.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Erwan Vautier <erwan.vautier@gmail.com>\n",
    "#         Nicolas Courty <ncourty@irisa.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "from mpl_toolkits.mplot3d import Axes3D  # noqa\n",
    "import ot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample two Gaussian distributions (2D and 3D)\n",
    "---------------------------------------------\n",
    "\n",
    "The Gromov-Wasserstein distance allows to compute distances with samples that\n",
    "do not belong to the same metric space. For demonstration purpose, we sample\n",
    "two Gaussian distributions in 2- and 3-dimensional spaces.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_samples = 30  # nb samples\n",
    "\n",
    "mu_s = np.array([0, 0])\n",
    "cov_s = np.array([[1, 0], [0, 1]])\n",
    "\n",
    "mu_t = np.array([4, 4, 4])\n",
    "cov_t = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])\n",
    "\n",
    "\n",
    "xs = ot.datasets.get_2D_samples_gauss(n_samples, mu_s, cov_s)\n",
    "P = sp.linalg.sqrtm(cov_t)\n",
    "xt = np.random.randn(n_samples, 3).dot(P) + mu_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the distributions\n",
    "--------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsvXl0XWd57//Z+wyaj+bhaJY827Ila7ATCqGhGFLKcCHhUlpIm5KbNoVeEtILK73cRVIWiykJIYU2pbQ0lBZK+AENEEIphDgJxYNkWbImS9Y8WvN45v3+/pD3ydF49paOpCPr/aylhS3t8+5XJ+a7n/N9n0ERQiCRSCSSmwd1pzcgkUgkksgihV0ikUhuMqSwSyQSyU2GFHaJRCK5yZDCLpFIJDcZUtglEonkJkMKu0QikdxkSGGXSCSSmwwp7BKJRHKTYd2Jm2ZkZIji4uKduLVkD1BbWzsmhMjcodvLUm7JVqIYuWjTwq4oSgHwTSCbxX/UXxNCfHm91xQXF3Px4sXN3loiWRVFUXp2eg8SyU4SiYjdDzwkhKhTFCUJqFUU5edCiOYIrC2RSCQSk2zaYxdCDAkh6m78eRZoAfI2u65EIpFINkZED08VRSkGTgLnIrmuRCKRSIwTMWFXFCUR+P+AB4QQM6v8/D5FUS4qinJxdHQ0UreVSCQSyTIiIuyKothYFPV/FUJ8f7VrhBBfE0JUCyGqMzN3KmFBIpFIbn42LeyKoijAPwItQognNr8lSTgeeWSndyCRSKKZSETsvwV8EHiToij1N77eFoF1JWvw6KM7vQOJRBLNbDrdUQjxCgaT5iUSiUSy9ciWAruERx4BRVn8gtf+LG0ZiUSyHGUnhllXV1cLWXm6cRQF5AzytVEUpVYIUb1Dt5f/ZSRbiSF3REbsEolEcpMhhX0X8qlP7fQOJBJJNLMj3R0lm0P66pJoQtM0XC4XiqJgt9uxWCwoisyn2EmksEskkg0jhMDr9RIIBNA0Da/Xi6IoxMTEYLPZUFVVivwOIIVdIpFsCF3UhRAoihIUcSEELpcLl8uF1WrFbrcHRV6yPUhhl0gkphFC0NfXh8ViISMjY8nPFEXBarUihAjaNC6XC1VViYmJwW63yyh+i5GPUMkKosXDj5Z9SJaiR+rz8/N4vd41r9OjeIvFgqqq9PX10d3dzfT0NAsLC/j9fnYi3XovIIVdsoJoaVkQLfuQvIYu6pqmBW0XIyiKEvxSVRWv18vc3BwzMzO4XC4CgYAU+QgirRiJRGKI5aK+UTtFURQsFktwTY/Hg8fjQVVV7HY7drtd+vGbRL57EiB6WhZEyz4kS1lN1M1E7Guhi7yqqgghcLvdzMzMMDs7i8fjQdO0CP0GewvZUiAMjzyy90QlWloWbHQfsqVAZBFC4PP5CAQCSyL13t5eFEWhoKBgSXbMWvT19WG1WnE6nWHvpx+8KoqCzWbDbrdjtVrloatsKRAZpM8r2cusJerLMSK4RiN83Ye3Wq2oqorP5wv68fLQ1RjSY5esIFpaFkTLPvYqQgj8fv+aoh4JKyYcy/34yclJrl+/TmlpabAISv+55DVkxL4Ke93njZbfM1r2sRfRRd3v968ZqW+3LaLfz+fzoSgKLpdL+vFrICP2VQj11aPFb5ZItgsjoh567Xaj72m1Iiir1UpMTMye9+NlxC6RSIIIIZifn6ezszOsqG+HFROO5UVQgUCA+fn5PV8EJYU9DNLnlewl/H4/Xq+XiYmJqIx418u8WS7ye7kISloxYZA+r2SvoNsvek55OMxG7NsZ4e/1Iigp7BKJBL/fHzyUNBqpR2NEvxq6yOv58W63m4GBARwOBw6H46bsPCmFXSLZ4ywXdaMRO2z/4Wm4Iqj1CH1oTUxMLCl6utmKoKSwSyR7mLUi9a2wYqINvQhKL8Lyer2oqhoU+d08CUoKu0SyRwkEAquKejQLdqT2FRr5L/fjvV4vHo8Hi8USHBKy24qgpLBLJHuQ0DF2G60o3akHQCSi6LUsnVCRX54fv5smQUlhl0j2GJqm4fF41q0oNSrsZlAUJWqqQ4149aqqBs8bdlsRlBR2iWQPEU7UwVwkLoRYLM3eJsHeCismHKHvlRAiWAQFRO2hqxR2iWSPoHdJjIuLi0hFqQIkvfAC1l/9Covfj+/MGfx/8idg3VpZ2Uorxsi9Q0W+qamJnJwckpOTg03J9KHeO0n0m0USiWTTaJrGzMwMbW1tYUXHqCjZz58n7fvfh7Q0RHY21uefx/qDH0Riu1vOZtImdfTX69G62+1mdnaW2dlZ3G73jtpOUtglkpuc0IPSSGJvbkaz28FuB4sFkZKCpa4uovdYzk5YMeuhaVowQl9vEtR2HzJLK0YiuYnR0/eEEFgslohGkVp6OlafL/h3ZWGBQGbmmtdHKosmUg+orbB0lls1mqaxsLCAy+XCZrPhcrlwOBzY7fZN33s9ZMQukdykhIq6HlVGMnL0nDmDNy8P+vpQBgYQSUn4/+APIrb+VqKP3dvKdVabBPXQQw9x6dKlTd83HDJil+zJua43O6HDp/W864jnnScmMvDQQzh8PgIeD76DB1FSUiK3/ipEykKByEXsRvLadavG6/USExOz6fuGQ0bsEjnX9SZjNVGHxbzsSFoxiqIgYmIQNTUMlZQw5vNFTZ56OCL1gFv+HofD4/EQGxsbkXuvhxT2DSCjW0m0Eirq2zGjVAjB4OAgg4ODTE5OcuHCBTo6OoJ53tFKpCJ/s+tIYY9iboYId6/Pdb0ZWS7qWy3s+tzRvr4+Tpw4wYEDB6ipqSE5OZlr165RW1tLf38/vpAD1s2yU4IcqXU8Hs+2WDER8dgVRfkn4O3AdSFEWSTWlGwtcq7rzUU4UYfIC/vs7CxTU1O84Q1vCI6lU1WVzMxMMjMz8Xq9DA8PU19fT1xcHLGxsdhstojdfzNEOt3RKLvNY/9n4I4IrRWVyAhXEq3obWfXE3WI7GCM2dlZurq6SE9PXzN1z263U1hYSHV1NYWFhczOztLX18e1a9dYWFiI2F42yk4cwm6XFRORiF0IcVZRlOJIrBWt3MwRrpzrunsRQjA2NobVag3bKiBSuFwuGhoaOHjwIGNjY0v2stYnBYfDgdPpxO12ExcXR3t7O4FAgJycHLKysrAabEMQyayYnWC7InaZ7iiRnzp2KUII/H4/w8PDJCcnEx8fv+X39Hq9XLp0iWPHjpmatKSjqirZ2dlkZ2fj8XgYGhqirq6OxMREnE4nKSkpu1q4jbAdbX+3TdgVRbkPuA+gsLBwu267JcgIV7LT6KJuZvj0ZvH7/dTV1XHw4EFSUlKYnZ3d1H1jYmIoLi6mqKiI6elphoaGaG9vJzMzk5ycHOLi4iK4+73Ftgm7EOJrwNcAqqurd7WRISNcyU4SKuqhg5q3Ek3TqK+vp7CwkIyMjIiurSgKKSkppKSkEAgEuH79Oq2trQBBqyZ0wtFujui3q2eMtGIkkl3EclHXv7ayMEgIQUNDA5mZmeTm5ga/bybLRr16lbRnnkFzubDcdReB171u1essFgtOpxOn04nL5WJ4eJja2lqSkpLIzc2N2pF9RtjOvUfE7FEU5dvAfwOHFEXpVxTlQ5FYVyKRLGW5qEPkK0p1hBAIIWhpaSEhIYGioqINraN0dhLz8MPENjYS09GB/bOfxfLyy2FfFxcXR0lJCTU1NeTk5DAwMEBnZyeTk5N4PJ4N7SUa2I5PHJHKinl/JNaRSCRrs5qow9YIux6Nd3Z2IoRg//79a14TDsvZsxAIEEhPJxAIIPx+rD/6EYE3vMHwXlJTU0lNTWVgYIDJyUmamppQVRWn00lGRsauGTa9XTaStGIkkl2A3+/H5/OtWVG6FcLe29vL7OwsFRUVa6YxGmL54a4QiA0KscViweFwUFhYyMLCAkNDQ3R3d5OSkoLT6SQpKWnbxFP/RGPm+u1CCrtEEuWsJ+oQ4Yjd5UJpbSWmtZVRTaPqllvWFUojYhV405uw/uQnWEZGAFCsVvx33rmh7YUensbHx7Nv3z5KS0uZmJigt7cXl8tFdnY2OTk5W97z3GhnRx2fz7dtlbdS2CWSKCacqIP5VgFrZpaMjmJ94AH8fX0cmJsj8fRpOHkS1kg7NHpfkZ+P57HHcP37v+Obm8Ny551ox48b3m84FEUhPT2d9PR0fD4fIyMjNDY2YrPZcDqdpKenb0nuuNkMHbfbvS3FSSCFXSKJWgKBQFhRB3MRuy7Gq61n+Yd/INDXx3RcHEpcHJamJsT3v4/2h3+45lqr/Xk1REEBM3ffzcLCAknFxYb2uuo6yx8kQqAMD4PbjcjPx2azkZ+fT35+PnNzcwwPD9PZ2UlaWhpOp5PExMQN33s5G2nZu13CLrs7SrYUmfO/MbxeLzMzM2FFHUwI+8gIh774RexvfjOWBx6AycklP/Z3dTEjxGL1p6ouzjIdGFh3yZ1IPwy+H0Jg+9zniLvzTuI++EFiP/hBlNHR4HWJiYns37+fmpoaUlNT6erq4uLFi/T19UVk39Ha2RGksEu2mJuhxfF2o2kao6Oj9PT0GBIOQ5aI243tfe8j7fx5lMFB1Oeew/rBD0IgcOPHbrrT0ki+MRsVTQO3G8rWbta604VCll/8Att//AciPR2Rno7a24v9859fcZ2qqmRkZHD8+HFOnDgBwMLCAo2NjYyNjW34fMKsx75dfWJACrtEElVomobH4zFVTWokYldaW1GGhggkJCBiYiApCeXqVejrC/Z/SXngAZTbb0cZGsI+MYHvPe9Bu2P9pq3bHbGH3k/t7Lzxh0UZEw4Halvbuq+32+3k5+eTkJBAcXExExMTGx4OYnZuqozY9zBmrItotTlki+ONEwgEgm0CjEaShqyYmJgVKYeKpqFZrVy6dIl9+/aR6nQS+Ou/xvfcc7Q+9RSe++8PiuZqmD20VRQFZXYW67PPYv3a11DPndtQm9SY9nZsTz2FcuUKeL2Lny4AZWYG7cCBsK/XLZSkpCQOHjxITU0NDoeDjo4OamtrGRgYwO/3h13HrMe+nYenUtijDDPWRbTaHI88svj/V/3/s/qfpbCHR1XV4HR7swei6yEOHkScPo1lYQFlZgbm5/HfcQeXhofJz88nKytLXwySktBiYyMejSsuFzlf+hLW734X60svYf/iF7E8/7ypNWJaW8l67DEsr7yCOjgImoYyOIgyPo6Wm4v3E58wtpdlBV5ZWVmUl5dTVlaG3+/n0qVLNDU1MTExseb7YNZj304rRmbFSCRRiD6RyOi1YR8CFgv+f/xH+j/zGfJmZlArK7lcVkZaSgp5eXkrLjfysDAbsdubm7GPjCAOHVr8hseD7bvfJfC2t7328S4Mif/5nwibDZGdDYAmBIGaGvx33okoLAQDwrnenmNiYigqKqKwsJCZmRmGhobo6OggIyODnJycJa2R5eGpZF3MWBe7zeaQLY43hpmI3XDbXrudsXe/m/lPf5qWmhpi4uMpKSlZ9dKtEHZF05YKuKqCActjCUIstYdUFZGcjDhwwJCoLy4RXpAVRSE5OZnDhw9TVVVFfHw8V69epa6ujqGhIfx+f1SnO8qIPQowM51pt01yitYHTrSykcZeZloKKIoSHDBdFprx4vOBxwMJCXAjxTLSVoz30CECSUkow8OIuDiUyUn873634WgdYO7224lrbEQZH1/01oUg8OY3m9qH2UjbYrGQk5NDTk4Obrc7OBzEbrcHD7mNrCeFXSLZ45iN2I1e63K5cLvd1NTUBMVI+clPsPzTP4Hfjzh2jMAnPrElEbvmcDD04IOU/vd/o4yP43/nOwm8/e2GXw/gKitj7MEHcdbWLtpL73jHYrRugs30dI+NjaWkpITi4mL6+voYGhriwoULZGZm4nQ6151nKoV9D2PGupA2x83LVlgx169fx+VyUV5eHrQQlOZmLF/7GuTkgM2G0tKC+vd/j/K2txkSdjMoioIvIwPfRz9q6nXL8R4/jvctb9nw6yMxrENRFBISEsjIyKCoqIjR0VGam5tRFAWn00lmZuaKjpPbNcgapMceddwM6Y6SzWMmGjZixUxMTHDt2jUyMjKW+MJKd/fiH+z2RUskKwu1sdF4H5ho9wJXIVJTmHSP3Wq14nQ6qays5NChQywsLFBbW0trayvT09PB92ijWTGBQICTJ0/ydhOfbqSwS3YlN/tDzYzwhIvuZ2ZmaGlp4eTJkysKn0RaWtCrBmB6GpGba9iK2W4iIcqREvbV1omPj6e0tJSamhqysrLo7+/n4sWL9PT0MDU1tSFh//KXv8yRI0dMvUYKu2RXEq05/JtlI4KzngjrpfMVFRXExsaueAiIU6fQ3vhG6O9f7AtjtxO4//4tjdjV5mbsjzyC/ZFHUJuaTL9+s0Q6Yl8NRVFIS0vj2LFjnDx5Eo/Hw7/+67/yd3/3d/zyl780fI/+/n5+8pOfcO+995ram/TYJZJdzloi5fF4qK+vp6ysjISEhOC1S8RYVdEeegjxrneBy4UoLgaHA6W1dUsi9tjWVuIefngxAwew/eAHuL75TbTyctNrbZStjNhXw2q1cvDgQf7n//yfpKSkELdGG+TVeOCBB/jCF77A7Oysqb3JiF2ya9htOfw7ic/no66ujkOHDpGcnBz8/qqRuKIgDhxAnDgBDsfa10WA9GefRfh8iJSUxS+fD9vXv770ovl5LL/+NZZXXoGZmSU/iiYrxmweu9frZd++fdx6662Grv/xj39MVlYWVVVVpvcmI3bJrmG35fDvFIFAgPr6ekpKSkhPT1/yMzOCHfGWAoqC4vWuLFJyu1/7+9QUsQ8+iDI4uLiHjAw8Tz6JyMyM6F4igdnujmbTHV999VWee+45nn/+edxut97G+VtCiA+Ee62M2CWSmwghBA0NDcHxcMsxmkZpduKQ0QfG5NvehgKwsAALCyhC4H/f+4I/tz77LMrgIMLpRDidKOPj2L75zeDPd6KP+lpsJGI3I+yf/exn6e/vp7u7m+985zu86U1vwoiogxR2yS5lr+Twmx1519TUFBz2vBpbFbF7vV5cLlfY6+ZuuQX3E0+gHTqEdugQ7sceW1I5ql6/vph6qe8hLg5uzErViYQVE4lRedHcK0ZaMZJdyc3qqy/vOqhp2opCl7Vob29HVVVKS0vXXT+sYPf1kfWVr5DY2orV50OUl6Pdfz/i4MFVL/f7/dTV1QVTKXNzc1ct0NEJvPWtBN761tV/Vl2N5Ve/WmxxcKPNr3bq1Pr7NYrLhTI9HZm12N7ReL/927/Nb//2bxu+Xgq7RBKlmBF2j8eDy+XixIkT60aRYYuZJiex/tVfkdjcTEx3Nyogrl5FfeUVfN/9LixrGiaECPr5qampeDwehoaGqK2tJTU1ldzc3GBGjhECb3kLvsFBbM8+C0Lgv/POxX4yIffbSMSutrVhe/pp8HpJDgSYevvbITOTmM98BrW9ncCpU3j/8i/BRGXoRiL27ao8lcIukUQZelRt1A8fHBzE7/dz/PhxQ/NR14vYlbY2mJvDNjaGsNkW2wz4fIi5OdSf/hTtz/88eK0QArfbTW5uLtnZ2Xi9XuLi4igtLQ1OJ+ro6CAQCOB0OrFareE/LSgK/nvuwf9Hf6RvOOzvHxa3G9vTTyNiYyEzE21igqzvfIf4Rx9F7e8HTUOtrUVtasL9b/9muCnZRjx2KewSyR7HiLCPjo7S19dHfHx8ZOaj2u3BiUTBCSmh+aUhdN4YTVdcXLxiTX3OaEZGRrAj4tDQEFarlbm5ORITE9ffaCQE/QbKzMxi3vyNzBotLg77+DjK8PDifSwWFCGwnj2LMjqK0IeOhCGaPXZ5eLpLuFk9ZcnahBP2qakp2tvbg60CjET34awYUVaGOHwYf0oKqsezmEe+sIAyOgpzc8Hh14ODg0xNTREbGxtW3PSOiIcOHSIuLo7Ozk5qa2sZGhoyPEwkuL8NWDHC4Vjs1T43ByxOchKqaqpd8GpEcz92Key7hJu1hF6yNusJ+9zcHE1NTZw8eRK73W64w2PY6+x2Ao8+yvjHPsb87bdDUhKioADt9a9H/c//RP3hD5mYmKCnp4fy8nJTIqsoCnFxcZw4cYJjx47hdrupra3l6tWrzN0Q3S0hNhbfffehLCygDAygzMxw/d57Ebm5i59OfD4E4L/tNlP58tJjl9w0hBYJSbaWtYTd5XJx+fJlysvLg+XpRodtGMqKiY3FdfvtxF67RrzfD7o1YbXiP3uWlpwcqqqqsFo3Lh+hfc0nJibo7OzE5/PhdDrJzs42nAlkFO3IETyf+QzK5CTjPh9uYOH551cenpoQ6mj22GXEHsVEYwm9/OSw9aw3Rcnr9XLp0iWOHTu2xKc2etBqNI9dURT8ycmLlaI3EAsLDHm9nDhxImICpSgK6enpnDhxgrKyMrxeL7W1tbS1ta3aH2VTxUUJCYj8fLSEhMU10tLwPP44rueew/vJT5rKiNnIXgKBwKYehmaQEXsUI0vo9zbLxVrPFz9w4AApKSkrrjVqxRh9AMydOUNmYyP09yOAKU0j7v77SUpKMv276IQbJF1cXExRURETExN0d3fj9XpxOp1kZWVtmygaRdM0U8IeqYpXI8iIXRKWaPzksBcIFWFN06ivr6ewsJDMVXzgiFoxN64LOBz4/+ZvCDz8MFff8x5mv/Ql0k6cWLnPH/8Y++23k3j6NDEf+9hiu4A11jSCHsUfP348GMXX1dXR2tqK50ZXyM1g1kJZi0hVsG4F0fUIlKzJTpbQy08OO4Mu7EIIGhsbycjIIDc3d91rw2FG2IUQkJRES04Oam4uuatUniY3NmJ54gmU2VmEqmL77ndR+voW88FttvC/ZBhCo/jJyUna2tqYmZnB7XaTnZ29o1G8GWHf7klT0fm4kaxARsd7D1VVCQQCtLS0EBcXR3Fx8brXGhVsMw+A3t5ePB4PB9YYGJ1aW4syPY1ITETExyPi41GvXEHt6Ah7DzPogysyMzMpKSkJ2lKtra3MzMyY7qkTqSZgG5n7uh3IiF1iir3SfCsaUFWVkZERbDbbmsKqY1SwzTwA5ubmcLvdVFdXrylI/ri41wqaYPHPMTHBfPdII4TAZrORnZ1NYWEhk5OT9Pb24na7gxk14aL4nRL27RwlKCN2iSl26pPDXvrEogvA9PQ0CwsLHDt2zFCrgEhaMS6Xi/HxcSoqKta1G0bOnEHk5MDU1GKFp89H4ORJtH37wt5jo+jvhR7Fl5WVceLEiWAU39LSsm4UH8kJSmY89u20Y2TELtkVPPro3hL34eFhZmZmyM/Pj0yrAIPXKa+8gvjnf8YxNMTBt74Ve5iOgr60NHw/+xnqF74AXV0Eysvx/dmfwSrj37ZqKhOA3W6nqKiIwsJCpqam6Ovrw+VykZOTQ05OzpIofifa9gYCgW09aI2IsCuKcgfwZcACfF0I8blIrCvZGWQR0s4yNjZGV1cXxcXF+P1+Q68xM0BjreuUy5dRP/1pJjSNhIQEHM8+i1pSgvaOd6y/qNOJ//HH8fl8W243GJnDmpqaSmpqKl6vl+HhYS5dukRiYiK5ubk4HI6Iph0aXWc72wlABKwYRVEswFeB3wWOAu9XFOXoZteV7BzRUoS0V9MsvV4vlZWV2Gw2Q2IN5rzzNa979VVmXC7is7OxJiXhT0pCefFFM1vfFoyKqd1up7CwkOrqanJycujv7+fixYtMTk4afl8jxXYLeyQi9lNAhxCiE0BRlO8A7wKaI7C2ZA+zV9Msc3NzCQQChht7webz2IUQDE5NkaWqxMbF4fZ4UL1eMNFLfV2EwH7xIlkvvYS1pITA7/4uwumMzNphCI3ifT4fzc3N9PT0MD09TV5eHg6HY8s/aey6iB3IA/pC/t5/43tLUBTlPkVRLiqKcnF0dDQCt5VEkr0aHUcj67UUWIvNHp52dnYy98Y3EpOXh9Lbi2VgAAFoH/xg2DWNfFKwnD1L4j/9EzG9vVjOn8f+2c+ijI+Hfd1q99qMCNtsNpKTk9m3bx+5ubkMDAxw8eJF+vr68Pl8G143HLtR2A0hhPiaEKJaCFG9WuWcZOtYT5z1nz3yyGvtt+G1P0eLsO/FNMutEPbVrtNb8O6/9VYCTz1F4MMfxvVHf0T///2/a47DW431BN7ywgtoGRkEUlIWh1TPzqI2NBheO5Loh6cpKSkcPXqUiooKAOrr62lubmZqairih7y7UdgHgIKQv+ff+J4kSljPM48WPz0c0fKA2U7MCPtGs2JCW/Cqqgrp6WjvfCeed7wDb3Z2xO4b/CgYGjnsYDl+aNRvs9koKCigurqavLw8BgcHIx7Fezwe7CFDureaSHjsF4ADiqKUsCjovw/8QQTWlewQezE6jkb0ytNIXhsqxHNzc7S0tKzagtd25QoZP/0p6v79i1kxBj5lr2eR+N/+dixf+QpWVV0cdJGaSmBZ3xllcDA4yDrwutchVinKikRGy1oPIkVRSE5OJjk5GZ/Px8jICPX19cTHx5Obm0tKSsqG772dLXshAsIuhPArivIR4Gcspjv+kxCiadM7k2yKRx5ZGo3r/x510V7rZzLVMXowa8UYSY3U1/R4PDQ0NKzaglc5exbHxz9OjNeLxW5H/cEP8H/jG5CRsWK90Pms66HdeisLQuD61a9I3LePwJvfDKmpr60zNIT9M59ZrFhVVawvv4z3wQfRjh0z9PubwcjDwWazkZ+fT15eHjMzMwwODtLR0UF2djY5OTmmo2+3273rsmIQQjwPPB+JtSSRIVxGyV7MNtkthB6eGvV6zWTF6J0iDx48uGoLXsvTTxOIi8OfmIhwOFAGB1F/8Qu0971v1fWMWjG+sjLGMjJIP3RoxY/V3/wGvF5EwQ1Xd3wcy89+tkLYI+F9m4n6Q6N4v9/P8PAwDQ0NxMbG4vf7Da/l9Xp3nccukUi2gK04PAWYn58nLy+PjFUicADcbljebyUC7XJXZWEBy9mzWM6fX2xJoKOqS3vQhBAJK2Yja1itVvLz86mqqiI/Px+/38+FCxfo7e3FGzKQZDW2O2KXwr4HWM8zl3569GJW2I1Es+3t7aiqSn5+/prXaP/jf6BMTaEsLMDkJMTEIF73ulWvNdMmYMW1Hg+2Z57B8l//BZqG2tWF2tAQGX9eAAAgAElEQVSAMjaGMj1N4Hd+x9C6ZtmsT68oComJiSQkJFBZWYnFYqGhoYErV64wMTGx6vux3RG77BWzBzCS7iiJPswcnhqxYnp6egyl3Wkf+AA+vx/f978PBQX4/+zPEPv3r3nfjaL09qIMDCCKigDwJyai1tWhHThA4Lbb0MrLV7xmu62YcGtYrVby8vLIzc1ldnZ2iRfvdDqDXrzZQdZut5vbbrsNj8eD3+/nrrvu4lETKWwyYo8QUiDlexBpImnFXL9+nZGREY4fPx5e1FQV3+//Ph2f/jT+v/s7xCoCG8pGxVYRYsnwaJGRgVZRge8jH1lV1IOvW2//LheEOUSOhLAvn8KkKAoOh4PDhw9TWVmJ1WqloaGBxsZGJiYmTEfsMTEx/PKXv+Ty5cvU19fzwgsv8Jvf/Mbw66WwbxJdzHZLPvhWIt+DyBDaltbMa9YS2OnpaTo6OsK24N0IK/aoaVj+8z+xfutbqFevrvtaLT8fkZKCMjgI09MofX0Ebr11idgbZnYW+yc/Sdx730vcXXdh+dGP1rw0khH7auhRfHV1NUVFRfT39/PlL3+Zl156ieHhYUPr63YPgM/nM91gTQr7JpFiJokG1orYXS4XV65coaKiwlSK3oaaimnaorD+0R8R+3/+D/FveAPW555b+8Xx8fg+9CECFRWIrCz873znYhrkOqwlqPann8ZSX4/IzUWkpmJ/+mnUxkZTa5jB6NxUh8PBiRMn+IM/+ANSUlJ4br33YxmBQICKigqysrI4c+YMp0+fNvxaKeybIDRlUP/fvdJfJbQVgewxs/OsJuw+n4/6+nqOHTtGfHy86TWNWiz6ddaf/QzLr3+NMj+/+OVyEXv//cF82lU/VaSkEHjXu/DffTfaLbdsuBpVra9HZGQs/uO78QBTOzvX3O9OZNacOXOG++67z/D1FouF+vp6+vv7OX/+PFeuXDH8WinsG0AXs+XR+qc+FV39VbYS/XeP9h4ze4XloqnnqpeUlJCSkrLp9da7LsjQ0MqiiLm5sJ53JBA5OSizszf+sviPUKSlbdn9zI7F20xWTEpKCrfffjsvvPCC4ddIYd8Aq4mZ/n2JZCcIjdiFEDQ1NZGRkUFOTs6G1jOTxiiEoKuri7qYGLSQLB5hsaAdOQI224b2sNa9VhNU70c+grDZUEZGUAYHCbzudYt+vYk1zO7DzHmF2ayY0dFRpqamgEU77ec//zmHDx82/Hop7BFiL+SDh7NdjL4H8gEYeUI98c7OTiwWC8XFxRtez0zEPjY2xvj4OGXvfz9TX/gCgZgYhKriLSlh7t//fcN7MIMoKcH99NN4PvlJPF/4At6HH15ZZKVfu0PCbiZiHxoa4vbbb+fEiRPU1NRw5swZ3v72txt+vcxj3yS6mO0FsTLapiAce21+qVmWi44RIdLz2PUWvCdPntx0EY4RYff7/fT09FBTUwOA5Z57WPjjP8YzM8PgxATXh4dJ9/nIy1sxoiHypKSgVVeHvSxSh6dm1jAr7CdOnODSpUsb2RoghX3TSIGSbCW6wIYTEVVV8Xg89Pb2Ul1dvem0RiPC7vV6mZ6e5sSJE9jt9tfK6hWFmORkSpKTKSoqYmxsjNbW1mAmyWaEdacOPje7huwVI9kVmLWeZPaMOcxOUXK5XCwsLFBRUbGiBe9qGBkKvd41mqZx+fJlkpKSSFhnfJ6qqmRlZXHy5EkKCgrweDzB/ipbObFoPbYz3VFnNw7akOxBzAqyzJ7ZGEaE3ePx0NjYSFxcnKEDOiPReLhrWlpayMjIIDY21vAha3x8PA6Hg8rKSlRVpb6+npaWFmZCm39tAzsRsZs9PN0sUtglkigmnLAHAgHq6+s5dOiQ4QjSqLCvRW9vL4FAgOLi4iVrGW1/oHdJrK6uJicnh56eHmpraxkaGlq7N44QKN3dxHV0oMzPG7rPWuxUxL6rBm1IJGbZCxlEkWI9YRdC0NDQEGzBezVMCb+OmVTG5YyPjzM0NER1dfUScdQ0LfilqioWi8XQPlJTU0lNTcXtdgdH0mVkZJCbm0tcXJz+i2J78kmsP/4x+9xuYp55Bt+XvoQoKdnQ77BTEbu0YiQ3NdJ+Mc56wt7W1kZiYuK6LXjNrrkeCwsLtLa2UlFRERRuPRtH0zRsNhsWiwVN0/D5fAQCAcP3iY2NpbS0lJqaGhITE2lpaaGhoYHx8XHUc+ewPvccIiMDX0oKysICMZ//vOn96+xExC7b9kokkiUj51YTx56eHrxeL8ePH9/w2mbw+/3U19dz/PjxFQIVCARQFCUYqVsslqCo6z/THwTh7quqKtnZ2WRnZzM7O8vAwAAzr77Kfq8X6w0xFg4Hak+Pqf2HIj12iUSyo6wm7HoL3rKysg0JlFlhF0Jw+fJlSktLcTgcS75vt9vp6upieno6uKaqqthsNmJiYrDb7SiKQiAQCI6SMxrFJyUlcfjwYYpuuw2A6fHxxTXGxtBWGa+3nUR7xC6FXSKJYpYLu96C9+TJkxvOVTczSxXg6tWrOByOJe0JdIEuKSmhqKiI3t5eLly4wODg4JIDUIvFEhT48fFx4uPjgyJvuNd8TQ3K/feTqmnETk8zm5zMxd/7PUZGRjZkKcH2j9cTQhg6d4gU0orZJYRWfUr2DqHCrrfgPXnyJLZN9F8xOvgaYGBgIJgfH0roQWnoAajeiTAzM5P8/Pyg/TA2Nsbc3Bzl5eUIIQgEAsEvi8WCqqooo6MonZ1gsaAdOADJycH7+T/wAfzveAetFy6w/7d+i/3A4OAg3d3dZGZmkpubu61Wh9mIfbuRwr5LkGX4exNd2DfbgjcUM+0Cent7qampWRKd6tG6oihLvh8bG8v+/fspLS1lZGSExsZG7HY7mZmZ9PT0UF1dHYxarVYrmqbh9/sX/fiBAWJeeAHsdhRNQ21pwf+ud0GI9UNyMp70dLBYiLPb2bdvHyUlJVy/fp2mpibsdjv5+fmkpKRsOiIPRyR8+q1ECrtEEsWoqho8uCwtLd1QC97lGBF2l8uF2+2mpqZmsZLV5ULp7UUEAvizs1GSk9cUNlVVcTqdOJ1OxsfHaWhoICYmhpGREZxOZ1DcVVXFbrcvPiRaWtDi4hA3onTL9euonZ1oyz4p6PsPvVdOTg45OTnMzMwwMDBAR0cHTqeTnJwcQ1W4G8FMxB6JOa1mid7PEhJZhr+HCR2P19fXR2ZmJtnZ2WFfY0REwqU7+v1+Ll++TGxs7GIu+fw8yg9+gPLii3D2LNbvfx9lYiLsfYQQ9PX1ceTIEaqrq/F6vZw/f56rV6/icrmW7McCWGw2LFbrolUE+L1eUx66w+HgyJEjVFRUoGkadXV1tLW1MTc3Z3gNo2wkYt/OCF9G7FFMuG6KkpufiYkJFEWhqKgo7LW6dx7ukG69B4Dw+7n2ox+xLyEBff6Q0tEB8/NoeXmLkerkJOrly2hvetO69+ns7CQhISF46FpaWkpxcTGjo6M0NTVhsVgoLCwkLS0NcewYlp/+FAFYAgGExUKgtDR40GqxWAw/uGw2G4WFhRQUFDA+Ps61a9fQNC1YyBUJpMe+S5CHk5JoY3BwEI/HQ35+vqFoT4/ENyzsfj/Tf/VXFDQ340hJwTI3B0eOgNe7ODRDb1Vrt4PHs+49rl+/zvT0NCdPnlyxx9A89b6+Ptrb28nNzSXvLW/B1tEBVitaWRm2zEwsN3LhA4EAQojglxEURSEjI4OMjAxcLhcDAwN0dXXh8/k2XQlqJmKXVswOEu1DqWUZ/t5ifHyc3t5eCgsLDQvDhgZQhzD5059iqa8n6ehRRH4++P2o//zP+PPyEG53cI4pk5OIAwfWXH9+fp7Ozk6OHz++rvglJSVx9OhRqqqq0DSN88PDNBcUMHvqFGRmBn8nm81GbGwsXq83mEoZ+mcjxMXFsX///mBL48bGRq5cucLU1NSGhNdMxO7z+TaVxbQRZMS+S5CfJvYWHo+HiooKJicnDbe3NZrGuJrHPjMzw0hbG0dTU4MDpQPx8SgjI2iZmfC7v4t66RIEAmi3347Yt2/Vtf1+P42NjRw7dsywmNlsNoqLiykqKmJ0dJS2tjYACgoKyMjIQFEU/H4/LS0twcpXPYrXP0WoqmpIaC0WCzabjaqqquBhq/6JITs72/Bhq5mIfbv7xMAeF/ZHHlkaqev/nfbSVCRJdJKXlxf0lg0X8hjsAbM8Ytfb/lbefjvqhQsIjwfsduyTk/jf9KZFASssRCssXHddIQRXrlyhuLiYpKQkQ3tevq+srCyysrKYm5ujv78/mOEyNTVFUVERiYmJwd9Vf290kff7/VitVkMCrygKycnJJCcn4/V6GRwcpK6ujpSUFPLy8tbtMQ/mInYp7NvMeoeTMvtEEg2YKSbaiBWjaRr19fUcPnyYuPR0Av/rf2H55jfB52Pm6FEWzpwh0+Belx+WbobExEQOHz6Mz+ejubmZqakp7HY7SUlJS8RdVVWsVuuSgqdAILBY8HQjkg+H3W4PfmIYGxujvb0dIQT5+fmkp6evuoaZiN3tdkthl0gkr2GmE6MZK0Y/hGxqaiInJ4f09HQAxJvfjO/229G8XhxeL319fXSeO0dubi65ublrWhWjo6OrHpZulvn5eTweD294wxuYnJykvb0dTdPIz88nMzMzKLp68zFt2WGrLv5Go/jMzEwyMzNZWFgIHrZmZWWRm5uL3W4PXmtmmPV294kBeXga5FOf2tt543vhd9yNmBF2M1aMpml0d3ejqiqFyywWDdCsVhISEjhy5AhVVVUIIbhw4QKtra0sLCwsuX5+fp5r166FPSw1i9frpaWlhQqLBfszz5D9s59RmZ3NkSNHmJqa4ty5c3R1db02a5XXDlvtdntQiAOBAD6fz1ROfHx8PAcOHKCqqgqbzUZDQwNNTU3BZmdmhllLK2YH0YVtr+aNy5YF0YlZYTdqxUxPTzM7O0tVVVXYdgE2m42ioiIKCwsZHR2lpaUl+EBwOBymD0uNoPv1RzweEr71LURSEorfj3rpEvF/+ZccOnQIv9/P8PAwly5dIjExkYKCgmD3ydBIPrR1gR7JG/XILRYLeXl55ObmMj09TX9/Py6XC5/PZzibRgq7RCIBzA+z1l9j5Fqv18vw8DCve93rloibEAK/37+iB0zo+vrh5uzsLL29vTQ0NJCdnb3p/jXL6ezsxOFwkHH2LCI1FVJSEIDS17dYHHWjXUB+fj55eXlMTk4Go/f8/Hyys7ODv1to64KRkRHi4uJWNCALh6IopKSkkJKSgsfj4eLFi9TX15OWlkZeXt66v7+0YqKIvZA3vpetp91CpLNivF4vQ0ND5OXlrfCM9bxwIxZDUlISsbGxOJ1OYmNjOX/+PO3t7bjdbkN7XY+xsTGmpqbYt2/fyo/OQgTTMXUURSEtLY3y8nLKysqYm5vj3LlzXLt2DU9IIZXX66W7u5sjR45gs9lQVTVo05jJiY+JiSEmJoaqqiqSk5Npa2vj8uXLjI2NrRrFy8PTKGIviJtsWRD9RNKK0TSNy5cvk5WVtURodGvCzIFg6GGpoigUFxczMjJCQ0MDsbGxFBYWkrxOo7C1cLlctLe3By2iwO/8Dtannwafb/ErNhatvHzN18fFxXHgwAFKS0sZHh7m8uXLxMXFkZ+fT0dHB4cPHw7+7qE2jX7oqrcuCPc+6P3V9U8w8/Pz9Pf309nZSXZ2Nrm5uUFratdZMYqivBd4BDgCnBJCXIzEpiQSySJ6VGmEcFZMa2sr6enp2O32JWvqomZU1PXD0lB/PrSj49TUFL29vbjdbgoKCpbYIuuhadqir37kSPDThCgrw/8Xf4FaWws2G4E3vAGyssKutdwbb25uxu/343K5cDgcq3aYXG2c33r7Dn1oJSQkLPH96+vrSUhIID8/37Sw9/X1cffddzMyMoKiKNx333189KMfNfx62HzEfgV4D/D3m1xHssPsBetpNxKpiL23txe/309JSQlDQ0P4/X6AoJjped/hMFJZqnvRbrebvr4+uru7yc7OJj8/f4n9s5yrV6+SlZW1ojWxOHSIwAZH4elj+WJiYqisrGRgYIDz58+TkZFBQUFBcDhHaEqkXuwU2oDM6PSjUN9/amqKrq4unnrqKQoKCnC73YaGgVitVh5//HEqKyuDB9xnzpzh6NGjhn/vTXnsQogWIUTbZtaQRAd7wXraTUTy8HR8fJyhoSGOHTsWPBjVrZfVBmashdnK0tjYWA4cOMCpU6ew2+1cunSJpqYmZmdnV1w7PDyMx+NZkXq5WTweD1evXqWsrIzY2Fj27dvH6dOnSUxMpLGxkcuXLzMxMbHkgRg6zk8fCOLz+UyN81MUhdTUVMrLy3nPe97D5OQkn/3sZw291ul0UllZCSyeZRw5coSBgQFTv7f02CWSKMbM4OnVHgILCwu0trZSVVUVjDr1NXXLwagPvtHKUovFsiR7paOjA03TKCgoIDMzk/n5ebq7u6muro5oHrxegHXgwIElVkiobTQ9PR3sMJmXl7diEMharQvMNA6LjY3lHe94h2k7BaC7u5tLly5x+vRpU68LK+yKovwXsNp/yf8rhPgPozdSFOU+4D4g4k9lycaR7YqjGzNCt1zY9YEZerQaumZoZaYRIlFZqmevpKWlsbCwQF9fH9euXcPn83H8+PGITzvq6ekhKSlp3R7sof1idJsmPT2d/Pz8YApjaOsCTdOYm5tDURR8Pp+h1gUej2dDvXPm5ua48847efLJJ4P5+UYJ+04KId5sekerr/M14GsA1dXVMv8iSpCFSTcPoVaMEILLly9TXFxMcshQaCFEcEydfsAYrrBotcPSdfF6sXznO1hefRXhcOC/5x7EMn84Pj6egwcP0tDQgKqqtLW1kZKSQkFBQdgGXEaYmppidHSUqqoqQ9fb7XZKSkqCHSabm5uXDAIJnWjV3t7OgQMHsFgswQfkeoetXq/X9KBtn8/HnXfeyR/+4R/ynve8x9RrQVoxEslNg6qqwRa/V69exeFw4HQ6l1yjaRoJCQnU1NQwNDTExYsXSUtLo6CgYNUim4204bV861tYfvQjRE4OjI5i++u/xvfEE4jc3CXX9ff3Y7PZOHr0KEIIxsbGaG1tDVa1hgqqGXw+H62trZSXl5uecrTeIJDc3Fz6+/txOBzB3joWiyVoa61V9GQ2K0YIwYc+9CGOHDnCxz72MVP7D/4eG3rVDRRFebeiKP3ArcBPFEX52WbWk2wPsjDp5kTPihkcHGRhYYH9+/cv+XnoYak+Pu6WW24hNTWV5uZmLl++vGTwxEbb8FpeemlRxGNjITUVfD6U5uYl10xPTzM0NMShG9kuegOuqqoqDhw4wMjICOfOnaOvr89U8ZAQgubmZkpKShbntW6C5YNAzp07R09Pz5IzBt2D13vT6Omp+iAQTdNMV56++uqr/Mu//Au//OUvqaiooKKigueff97U3jcVsQshfgD8YDNrSLYfWZgU/WwkUlUUhYWFBUZGRqipqVnRA2a1dgGhbQKmp6fp7e2lvb2dgoIC5ufniY+PN31YKhITF0fn6amNQkCIyHq9XpqbmykvL181jTAxMZGjR48u8b0zMjLIz88PK9b9/f3Y7fawg7/NYLPZKCgoYHh4mMLCQtrb24Glg0BgZU68/kAaHx83dX7w+te/ftPj9GRLAYnkJsHv93P9+nUqKiqWCInRdgHJyckcP36c48ePMzIyQk9PDzabzfAEJ53AH/8xytQUSm8vSnc3Yt8+tBtet56psn///rD9ZXTf+/Tp0yQlJXHlyhUuX77M5OTkqsI3OzvL4OAgBw8eNLVfI+hZMwUFBVRWVnLo0CHGx8f5zW9+Q09Pz5L3KLTDJMAvf/nLYN3AdiE99j2OLEyKXkJTHcMNdggEAnR1dZGenr4kqt1Iu4BAIIDL5eLWW29ldHSU2tpaUlJSKCwsNNTsS6upwff5z6O2tCASEtBuuWXRlgG6urpITEwkM9Po+I5FoczJySEnJ2fJp4r8/HxycnKC9kdTUxPHjx83XExklNHRUVwuV9A2gsVK08OHD+P3+xkcHKS2thaHw0FBQUHQtlJVlccff5y77757Qwegm0EK+x5H+urRj57GuJZgCSFobGwkMzNzRR672XYBoYelcXFxFBYWUlBQEMwU0b35lJSUdR80Yt8+Asvmoo6PjzM5ObmplEn9U4Xb7aa/v59z584Fe7UUFhZGJKMmFK/XS0dHx5oZQVarNfgejY+PB3P0U1NTuX79Or/61a/41a9+FdE9GUFaMRJJlBOu+rSzs5OYmBiys7OXXGe2XYBukyw/LNV9+OrqaoqLi+nv7+fChQsMDQ0ZrsR0u93BClCzmSqrERsby/79+zl16hQej4eJiQkmJiaYmZnZ9No6+kHs/v37122FAIvvUUZGBidPnuTIkSP8/Oc/533vex+nTp3C5XJFbE9GkRG7RBLlrCfsIyMjTE5OUllZydzc3BLrxky7AFi0SeLi4tY9LA2NmPU+ME6nc918eE3TaGxsXNJZMVK43W5mZmb4rd/6LWZnZ+ns7MTv9werWjfzEBkYGCAmJsaUbQSLHSa7urp46KGHSE1NNTW5KVJEnbDLSkiJZClrCfvMzAzXrl2jpqYmGJXrfrqZ3uqw6CNPTU1RUVFh6Hq9D0xJSUnQY17Lh29vbyczM5PU1FRDaxtF0zSampo4evQoNpstWNXqcrno7e2ls7Mz7ENnLfQ2vDU1Nab39fLLL9PU1MSXvvSliPv9Rok6K+bRR3d6BxJJdLBeIzCPx0NjYyPl5eVB0dKvM9suQK8s3YhNonvMp0+fJi0tLZgPr2eujIyM4Ha7KSoqMrWuEa5evUpOTs6Kcvu4uDgOHToUfODV1tbS0tLC/Py8oXVDHxhmhXlmZoZPfOIT/OM//uOOiTpEYcQukUiWslzY9YEZhw4dWnJYGBqxG43UQw9Lw/nI6xGaDz8zM0NPTw9tbW34fD5OnToV0eZeANevX8ftdi/JVFlO6MHm2NgYbW2LjWgLCwtJT09fc0+dnZ1kZmaa7s8ihODhhx/mf//v/01xcbGp10aaqIjYZSWkRLI2ocKuH3BmZ2evaG5lsViYn5+ns7MTr9cbdt21Dks3i8PhCLYJSEtLo66uju7ubtP58Gvhdru5du0aR48eNfTA0KtaKysrOXjwIKOjo8Gq1uX55ZOTk0xNTW1ImF944QXGx8e55557TL820iibrXDaCNXV1eLixdWHLclKyJ3jZjnfUBSlVghRvUO3j9i/Xp/Ph6ZpNDc343Q6SU1Npbu7m7m5uWBv9eBNb1SW6vZHX18fSUlJFBUVrZkCqB80RrqgR39gpKamkpeXh9/vZ2hoiIGBAVP58KuhaRp1dXXs27dvU569z+djYGCAoaEh0tPTKSgowGq1UltbS3l5uel2BGNjY/ze7/0eP//5z01X6prE0EcfKeySIDfLe3+zCXtra2swR72rq4vq6uolXvhq7QKEEIyPj9PT04PFYqGoqGhJ7vno6Ch9fX1UVFREJP0wlP7+fqanpzl27NiS7wshgvddbU9G6OjoQFVVSktLI7JXTdOCe1pYWCA3N5d9+/aZ2pMQgrvvvpvf//3f573vfW9E9rUOhjYWdR67rISUSJaiqioLCwv09vYGDwR11sqA0fOqMzIygp53R0dHsIjn2rVrVFZWRlzUZ2ZmGBwcXLVd7nIfPrQvjZG5qOPj40xPTwenC0UCvZuj3jzN7XZz4cIF8vPzyc7ONnQA+t3vfpfY2FjuuuuuiO1rs0SFxx7KzWAF7Cbk+Ub0Ehp9d3Z2cuLEiSUHnEbbBTgcDo4fP05ZWRmTk5OcO3eO9PT0iGdt+Hw+mpubKSsrC7u2w+GgrKyMEydOMDc3x7lz5+jq6lrThw8dcRfpg1i3201XV1fwPSovL8flcnH+/Hk6OjrweDxrvnZwcJAvfelL/M3f/E3E97UZos6Kkewc0oqJCBF7B/1+Pz6fj1deeYWsrCwOHz685Od6ZamZmaUNDQ1kZGTg8/kYGhoiMzOTgoKCTRcOCSGor68nLy+PrKws068PBAIMDg4GffjQgRtCCC5dukRhYeG605A2uu+6ujpKSkpIS0tb8jNN0xgeHqa/v5/4+HgKCgqWDC3RNI0777yTBx98kDvuuCOi+1qH3WnFSCSS12hrayMhIWFF1opZUYfXKkvz8vKAxbS/4eFh6uvrwx60hqO7u5uEhIQNiTosZvQUFBSQn58fHLih+/BTU1NhR9xtlJ6eHhwOxwpRh0WbJjc3F6fTydTUFN3d3Xi9XgoKCsjKyuIb3/gGpaWlvPWtb434vjaLFHZJEHm+EV309vbi9XrJyMhYkse+kXYBq1WWhgrX+Ph4cHpRcXGxqUPNiYkJxsfHI+J966mJmZmZzMzM0NHRwdTUFIcPHw72vYkUs7OzXL9+nerq9T/cKYpCamoqqampuFwu+vr6+NCHPkRfXx8vvvhiVFkwOlHnsUt2DumrRxdJSUlBvzo0j91su4BwlaX6QWtVVRX79u2jv7+fixcvMjIyEnbgg9vtpq2tjePHj0f8IDYuLg6v10tlZSULCwthfXgzBAIBmpubOXr0qKl9x8XFUVpaiqZpvOtd7+Lzn//8pveyFciIXSKJUtLS0vD7/cFZphvpre73+7ly5YrhylL9oDW030peXh55eXkrDkQ1TePKlSscOnQo4s29QkfcpaSkkJKSsqIvzWYGX+tzTBMTE02/9itf+Qq33HJL1Io6SGGXSKIefZCEpmlBC8YIeqFQYWGh6cpSvd+KXshz/vz5FQetHR0dpKenr+pPb5bVRtyt5cMXFhaSmppq+H0ZGxtbMTjDKE1NTfzwhz/k5ZdfNv3a7UQKu0QS5YQK+0YOS51O54bvbbPZKC4uXnHQmpCQwMLCAuXl5Rteey30EXdred+hPvzs7GwwR99IPrzX66W9vZ3KykrT3rjX69lRJYkAABYJSURBVOUjH/kIf//3fx/xTyiRRnrsEkmUo6oqs7OzwepSI+iHpfv374/YHnJzczl16hTJyclcu3YNTdOYmpra9ODlUPQRd0Zy4eG1c4gTJ04wPz+/rg8fOjhjI8L8uc99jne+852bmgAVytTUFHfddReHDx/myJEj/Pd//3dE1gUZsUskUY0QgqSkJNLT06mtrSU9PZ3CwkJib8wQXQ39sHQrKks1TaO/v5/q6moURVlS0brZwRYAra2tGxpxp09UCvXhk5OTl6w1MDCA3W43PTgD4MKFC7z66qu8+OKLpl+7Fh/96Ee54447+N73vofX62VhYSFia8sCJclNx81SoBQIBHC73cHDUk3Tgg2+EhISKC4uXiGAfr+f2tpajh49GtGOjTpNTU0kJyeTn58f/J7b7aa3t5fx8fE1D1qNMDQ0xNjYWESqS4UQjI2N0dvbG2wb0NPTw6lTp0zvbWFhgbe85S18+9vf3pAvvxrT09NUVFTQ2dlp9nc1dLG0YiSSKOWJJ57g2WefDaY6qqqK0+mkpqaG7OxsWltbuXz5MtPT08DmDkuNMDAwgBAiWOCkExsby8GDB6murkbTNEOl+MuZn5+np6eHI0eORCQvXPfhq6qqKC0t5erVq8HOl2ZG1Qkh+NSnPsU999wTMVGHxfOPzMxM7rnnHk6ePMm9995reBCIEaSwSyRRyvve9z7q6+u57bbb+Id/+IfgUOTQvPPi4mK6u7upra3lypUrxMbGbuqwdC1mZmbo7+9fV3j1g9bTp08THx9PfX09TU1NzM3Nrbt26MQiqzXy7vDo6ChFRUVUVVWZzod/6aWXaG9v58Mf/nBE9+T3+6mrq+P+++/n0qVLJCQk8LnPfS5i60thl0iilMLCQp588kn+67/+i8nJSd74xjfy2GOPMTU1FbwmOTmZ8vJysrKygkMihoaGIjpA2UxzL1h60JqTk8PVq1e5dOkSExMTqx60rjXiLhJMTU0FB2fExMSwf/9+Tp06hc1mCzsyb3p6mocffpivf/3rET+ryM/PJz8/n9OnTwNw1113UVdXF7H1pbBLJFFORkYGjzzyCOfOnSM5OZk77riD//f//h/Dw8PAoo0xMDDALbfcQkVFRbBbYm9vb7BKdaPo9k5JSYnpA01FUUhPT6eyspL9+/czODjIhQsXGB4eDj549BF3BQUFm9rnavj9flpbW1cMJbFYLEFRzcjIoLW1dcWDRwjBxz/+cT72sY9RWFgY8b3l5ORQUFAQHNf3i1/8gqNHj0ZsfXl4KrnpuFkOT9fC5/Px7W9/m6eeeoqysjJGRkb4xje+sSTiDZ0QlJWVRUFBwYZmmnZ3d+PxeCLmL4cetGZlZTEyMkJ1dfWm5q2uRVNTE2lpaYasqdnZWXp7e5mfnycuLo5r167xve99j+9973sRj9Z16uvruffee/F6vZSWlvKNb3zDyFSo3TlBSSLZLDe7sOsEAgFe//rX4/P5KCoq4sEHH+TkyZNLolNN0xgcHKS/v5/U1FQKCwsNj32bmJigs7NzS9ImvV4v58+fB16LXiNZ9DMyMsLIyAjHjx83dRjr8Xj46le/yhNPPMGf/umf8slPfnLDbQu2CNm2VyK5mfH5fPz5n/85H/jAB3j55Zf53Oc+h9vt5sEHH+SNb3wjqqqiqir5+fnk5eVx/fp1GhsbiY+Pp7i4eN0+KR6Ph7a2Nk6ePLklEWtvby+5ubkUFxcHK1oTExMpKiraUP+WUNxuN52dncFcezPYbDYuXrzI3/7t3zI3NxfxYSTbhYzYJTcdeyViX3FjIWhsbOSLX/wi7e3t/MVf/AXvfOc7l4iTEILJyUm6u7tRFGXVFr2apnHp0iWKi4tJT0+P+D7Hx8fp7u5eUtYvhGBiYoKenh4URaGoqMhU/xed9QZnGOHf/u3fOHv2LM8880xUtuNFWjGSvcpeFfZQurq6ePzxx3n11Ve59957ef/737+iWnV2dpbu7m7cbjdFRUVkZmaiKArt7e1YLJaIDYwOxePxUFdXR2Vl5ZrWi97/ZWFhgcLCQrKysgx/aujp6cHr9XLgwAHTe+vv7+e9730vL730EikpKaZfv01IYZfsTaSwv8bo6ChPPfUUP/zhD3n/+9/Pn/zJn6xIK3S5XPT09DA1NUVKSgoLCwsrvPpIYHbE3fKK1tzc3HXz3GdnZ2lpaaG6utq0faRpGu9+97v5+Mc/zpkzZ0y9djWKi4tJSkrCYrFgtVqJoN7JylOJZK+TmZnJpz/9aX79618TExPDmTNn+NSnPsXIyEjwmri4uGAjqpGREdxuNz09Pfj9/ojupaenx9SIu+UVrRcuXKC9vX3VitaNDs7Q+frXv86hQ4d485vfbPq1a/Hiiy9SX18fSVE3jBR2iWQPkJSUxEMPPURtbS1Hjhzhrrvu4oEHHqCzsxNYFMa2tjYqKio4deoUiqJw4cIF060B1mJqaorR0VH27dtn+rWhFa2JiYmrVrR2dHTgdDo3dPDa3t7ON7/5Tb7whS9Eq69uGmnFSG46pBUTnkAgwHPPPcdjjz1GXl4eVquVD3/4w0ta0mqaxvDwMH19fTgcDoqKioiPjzd9L5/PR21tLeXl5YZTLddj+UFramoqExMTG7KP/H4/b3vb23jssce45ZZbNr03nZKSkuDh75/+6Z9y3333RWrprbdiFEX5oqIorYqiNCiK8gNFUaL2xEEikbyGxWLh3e9+Ny+//DKlpaX85je/4dFHH+Xs2bNLmo7prQEyMjJobm6moaGBmZkZw/cJHXEXCVGHpRWtRUVFdHV14fV6TTf4AnjyySe57bbbIirqAK+88gp1dXX89Kc/5atf/Spnz56N6Prh2KwV83OgTAhxArgKPLz5LUkkku1C7xjZ0NDA5z//eZ555hne+ta38txzzy0Zmp2ZmUl1dTWFhYV0dnZSV1fH+Ph42CEbq424ixRCCPr+//buPaaqe0vg+HcVEKtoMfQyKNyqafFJEQcVUxXBjh183XuZOk0xOG2NtbYZxcRWOpqp3jZqhptMNU0gpffeNrT2Ecc0GuU69Ip1WlsOYkVSBtuIilDpkVGQRzW8fvMHcIqPW177nH04rE9yknMOm73WJmRl57d/6/errGT69OnExMRQX1/vWkqhN88HSkpKyM3NZfv27Zbn1rUCZmhoKMnJya5mLE8ZUGE3xuQZY7r+ggVAxC8dr5TyPmlpaQQFBTFz5kw++ugjcnJyOHHiBAsXLiQnJ+e2Mfbg4GBiYmKYNGkS1dXVnDp1CqfTec8C37XF3aRJk9yS95UrVwgICCA0NLRPD1qhY9rlhg0byM7Otnybu6amJhoaGlzv8/LyiIqKsjRGT6x8eLoG+IuF51NK2eDhhx8mKyuL3NxcKioqiI+P56233nIVK4CgoCDXlnR1dXU4HA6qqqpcd/l93eKur3766ScqKyvvWsOmNw9aAXbt2sXKlSuJjo62PDen08n8+fOZMWMGc+bMYdmyZSQlJVke55f0+PBURP4KhN3jR9uMMQc7j9kGzAL+yfyNE4rIOmAdwEMPPRRbUVExkLyV+pv04am1bty4wdtvv83777/P8uXLWb9+/V3by7W0tFBZWYnT6SQsLIzGxkZCQkIYN26c5fm0t7dz+vRpJk2axAMPPPCLx975oDU0NJSKigreeOMNjh075pb1393MMw1KIvIs8ALwuDGmV5v26awY5U5a2N3j1q1b5OTkkJWVxdy5c0lLS7trSdu2tjbKysqoqakhPDy8x/1Z+6O8vBwR6XNnbENDA3v27OHdd98lPT2dDRs2+GxhH+ismCRgC/Cb3hZ1pdTgNHz4cNatW0dRURGLFi1izZo1rF27ltLSUtcY+61bt2hsbGTevHmMHj2akpISSktLLdv2ra6ujtraWiZOnNjn3w0KCuL69ets3LgRp9Pp2lLQFw3ojl1EzgOBwLXOrwqMMet7+j29Y1fupHfsntHe3k5+fj4ZGRmuefAOh4OXXnrJtWxB11DIpUuX8PPzcy061h+tra0UFRX1ez78sWPHyMzM5MiRI25bY90DdK0YNTRpYfcsYwynT59m9erVBAcHs3nzZpKSku4qnjdu3ODSpUu0trYyfvx4QkJC+tRQVFpaypgxY/o1bl9bW8uSJUvIzc0lImJQT97TtWKUUu4nIkRERDB79mzee+898vLySEhI4MMPP7xtw+iu/VmnTJnC1atXKSws7PX+rE6nk7a2tn5t1G2M4ZVXXiE9Pd3Sot7W1sbMmTNZvny5Zee0ihZ2pdSAhYWFkZOTw+TJk8nOzubQoUN8//33LFiwgMzMzNvG2EeOHMm0adN6vT9r18YZU6dO7ddaLgcPHqSlpYVVq1b1+/ruZe/evUydOtXSc1pFC7tSynLjxo0jIyODEydO0NzczKJFi9i1axfXrl1zHRMYGEhkZKSrqaiwsJDy8nKam5tdx3QtSTB58mQCAgL6nIfT6WT37t1kZmZausBXVVUVR44cYe3atZad00pa2JXlduywOwPfcuvWLebMmcOMGTOYPn26W1rg3WXMmDFs3bqVwsJCwsPDWbFiBenp6VRVVbmO6d5UFBgYyDfffMO5c+e4efMmly9fJigoqF+7IbW3t5OWlsbOnTvvmnc/UJs2bSIjI8NrH8J6Z1ZqUPv97+3OwLcEBgaSn5/P2bNnKS4u5ujRoxQUFNidVp/cf//9vPjiixQVFfHYY4+xevVqXnjhBcrKylxTJbv2Z42Li2PMmDGcPXuWixcvEhZ2r/7Inu3bt48HH3yQFStWWHkpHD58mNDQUGJjYy09r5W0sCvl5UTEtc54S0sLLS0tg3bdcH9/f1JSUvj6669JTU1ly5YtpKSk4HA4XAVeRFybcURGRnL+/HnOnDlDbW1tj4uOdbl8+TKZmZm8+eablv+tTp48yaFDh5gwYQJPP/00+fn5pKamWhpjwIwxHn/FxsYa5Vu2bzcG7n5t3+75XIAiY8P/defLLVpbW82MGTPMyJEjzZYtW9wVxuPa29tNQUGBSU5ONgsXLjSffvqpaWhoMCdPnjRlZWWmqanJNDU1mR9//NEUFBSYzz//3Fy8eNE0Nja6fnbnq76+3iQmJpr8/Hy353/8+HGzbNkyt8fpplf/h3rHriyxY8fP5Rx+fq/j7dbw8/OjuLiYqqoqCgsL+fbbb+1OyRIiQlxcHAcOHCArK4uDBw8SHx9Penr6bUMwo0aN4tFHHyUqKorr16/jcDj44Ycf7jlVMjs7m+joaBISEjx4Jd5FG5SU5UR+LvD2xPftBqXXX3+dESNG8PLLL7s7lMddv36dBQsWMG/ePBwOB8899xypqal37dzU3NzM5cuXqampYezYsURERODv7893333H888/zxdffGHZxh5eRhuUlD0G0aSNQaGmpoa6ujoAbt68yWeffcaUKVNszso9Ro8ezf79+8nOziY/P5/6+noSExPJyMigtrbWddywYcN45JFHmD17tmt/1j179rBx40aysrJ8taj3mhZ2ZTkdfrFWdXU1iYmJREdHM3v2bBYvXuyV3Y5W8Pf3Z9q0aQCEhITw2muv4XA4CAkJYenSpWzdupUrV67cdvz48eOJi4ujrKyMCxcusH//frvS9xqDbs1KpYaa6Ohozpw5Y8m52tramDVrFuHh4Rw+fNiSc7rbiBEj2LBhA+vXr+eTTz4hJSWFqKgo0tLSiIyMREQoKSmhvLyc8vJyLl68aHfKttM7dqWGEG9ug+9JQEAAqampOBwOnnzySTZt2kRqaipfffUVGzdu5J133mH48OGD9vqspIVdqSHC29vge+u+++5j+fLlHD9+nM2bN/Pqq68SHR3N9OnTB3zuwdzl250OxSjL7Nih4+verKsNvvvepYOZiDB//nyKiop63bjUk64u36CgIFpaWpg/fz5Llixh7ty5lpzfU/SOXVlGlxLwXoOhDX4grOou9ZUuXy3sSg0Bg6IN3ku0tbURExNDaGgoixcvJi4uzu6U+kwLuxqQHTs6GpK6bmq63uuQjHfZvXs3VVVVXLp0iY8//phFixbxwQcf2J2WV/KFLl8t7GpAdCkB5auCg4NJTEzk6NGjdqfSZ/rwVKkhJiEhYUDrqEyYMIFRo0bh5+eHv78/vrQ8SE1NDQEBAQQHB7u6fNPT0+1Oq8+0sCvLDNKZYaofjh8/7lpa15dUV1fzzDPP0NbWRnt7O0899dSg7PLVwq4so8MvarCzssvXTjrGrpTqExHhiSeeIDY2luzsbLvTUfegd+xKqT758ssvCQ8P5+rVqyxevJgpU6YQHx9vd1qqG71jV0r1SXh4OAChoaEkJydTWFhoc0bqTlrYlVK91tTU5FqSoKmpiby8PKKiomzOSt1Jh2KUUr3mdDpJTk4GoLW1lVWrVpGUlGRzVupOtmyNJyI1QMUvHPIg8H8eSscb4toZ2xevebwx5lduOK8axETk10AO8Hd0bGGYbYzZa29W7mFLYe+JiBTZsWelXXHtjD0Ur1l5BxEJBv4IRNFRaNcYY752Y7yxwFhjzDciMgo4DfzOGPO/7oppFx2KUUrZZS9w1BizUkSGASN6+oWBMMZUA9Wd7xtEpAwIB7SwK6XUQInIA0A88CyAMaYZaPZg/AnATMDhqZie5K2zYuzqerCz20KvWQ0lE4Ea4F0ROSMifxSRkZ4ILCJBwAFgkzGm3hMxPc0rx9iVUr5NRGYBBcA8Y4xDRPYC9caYf3dz3ADgMPDfxpj/dGcsO3nrHbtSyrdVAVXGmK6hkP8C/t6dAaVjK6Q/AWW+XNTBSwu7iLwhIiUiUiwieSIyzoOx/yAi5zrjf9r55N4Tcf9ZREpFpL3zbsYTMZNE5DsROS8ir3oiZmfcP4vIVREZfDsYKEsYY34EKkVkcudXj+P+h5jzgNXAos7aUiwiS90c0xZeORQjIqO7xr5EZCMwzRiz3kOxnwDyjTGtIvIfAMYYty/ILCJTgXbgbeBlY4xbF7kWET/ge2AxHXdPp4AUT0z9EpF4oBHIMcZo2+IQJSIxdEx3HAZcAJ4zxtTam5Vv8MpZMXc80BhJxxxXT8XO6/axAFjpobhlYN2mvL0wBzhvjLnQGfdj4Ld4YOqXMeZ/OmclqCHMGFMMaB+DG3hlYQcQkZ3AvwA3gESb0lgDfGJTbHcLByq7fa4CBt+uvUqpu9hW2EXkr0DYPX60zRhz0BizDdgmIv8G/Ctg2f48PcXuPGYb0Ars82RcpZQaKNsKuzHmH3p56D4gFwsLe0+xReRZYDnwuLHwIUQfrtkTfgB+3e1zROd3SqlBzltnxUR2+/hb4JwHYycBW4DfGGN+8lRcG5wCIkVkYmc799PAIZtzUkpZwFtnxRwAJtMxS6QCWG+M8cjdpIicBwKBa51fFXhiRo6IJANvAb8C6oBiY8w/ujnmUmAP4Af82Riz053xusX9CEigY3VHJ7DdGPMnT8RWaijwysKulFKq/7xyKEYppVT/aWFXSikfo4VdKaV8jBZ2pZTyMVrYlVLKx2hhV0opH6OFXSmlfIwWdqWU8jH/D/n7sbPc7qYcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5c871e15f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = pl.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax1.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples')\n",
    "ax2 = fig.add_subplot(122, projection='3d')\n",
    "ax2.scatter(xt[:, 0], xt[:, 1], xt[:, 2], color='r')\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute distance kernels, normalize them and then display\n",
    "---------------------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAC7CAYAAAB1qmWGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXt0XXW1779zv7J3svNukqZN0jelBVoKKW+lCIeDeBBxqIAeBPGK9171yBWPIB6ur8EdoCLXK+i5RTjleERBEUHhIhQRUORRSoG+aNPSRx7NO81Odnb263f/aDijK99ZmrYhTRbzM0ZHsmfXWr/f+q25fntlfX9zTnHOwTAMw5j6BI52BwzDMIzxwSZ0wzAMn2ATumEYhk+wCd0wDMMn2IRuGIbhE2xCNwzD8Ak2oRuGYfgEm9ANwzB8whFN6CJygYi8KSJNInLDeHXKMI425tvGVEQON1JURIIAtgD4OwDNAF4GcLlzbuP4dc8wJh7zbWOqEjqCfU8B0OSc2w4AIvIrABcDOKDTR6TARVHksWVqimg7yfO++TDbnPL3RUGPsjMABIVM2Sjbgmn+gnPC2+Ui3EQkkSNbPsKdzBTx8SIJblcyfLx9/8H7I50mU7YsxrsqwxNMcTsuoIxXEZ9LXvGgUJLPJR9W+gwglOQOSTpLtlwRD3im3LtvtrMPucSg3tChMS6+7UoLeUPl+Skb4y6HE5of6w9fuVjwQN3yEFB8O6fcAxrBYaVtZVfJcL+HK7h/coDnSIkq91CO/S68d2z3pEYwo9gG2efyUe53LsLtavcUAISUY6ZL+YYJDfFgjJ4jMn09yCUP7ttHMqHPBLB7v8/NAE59px2iKMKpcq7H1nbFGbRdcIj3TdZqzsjbzX9gQG07G+dvhO7FfICSncqFVS5ifz1f7BnP9JEtWRcnW0cjD/vMZ1Jki7T2kw0AXJj3l12tZOv68HF8zAH2vuKmBNlyMR6vzpP5yzc1jftXtY7HcLBGn3SmvcrXK9TSTbb+U+rJ1vpR75dY6013qm0cBuPi28NnLaftAln2484lPBPNfIavCbL6zNGzpISNyq1fvIu/9HsWFfCuOe6jdl9obUTbkmRr+lQx2QIZfW4qOHYv2QZ6+Ytx5qPsT/2zFB9Tmom38DiWre0kW3J+BdkS9dqEzG0AQOWLfMzdF1Xzdhv4G6Z9uff+27nyh3ojo/sypq2OABG5BsA1ABCF8sRiGFMU821jsnEkomgLgP0fm+pGbB6ccyudc43OucYw+GnAMCYh5tvGlORIntBfBrBAROZgn7NfBuCT77RDpqaIXrHU3vY8bbd11clk+6flT5GtPtxDth//+VK17V2X8ru5QAf/6VW5gf+07LiW/6Ya6OMnsm21pWTLlHK7pZu4f5kivhT9p1fxhgCGy/jvyHyE/zwcbOC2w738Hd5+Cvcbyl/301/i45Xs5D/P25VXSrN/r78+evNz/J4/vnU22cKD3M55Czd7Pv8uyq+tDpND9m1XWkivWAoee5m2a7meXzEGlHe6bWcpryn4jQkA/RVJcQtfq1CSG4p18au1nsWK7qO8stT6XZHj61n2Jvtr73H6S/T0RvZFmcknnlBeeZZv4Xs32sH37vaP8WvQ8GAl2xJ8vKqX+JVSV6PyygtA55l8/ybr+MaqXsu29DHefrvoAV7Uj+KwJ3TnXFZEvgjgjwCCAO5xzm043OMZxmTBfNuYqhzRO3Tn3GMAHhunvhjGpMF825iKWKSoYRiGT7AJ3TAMwye868sW90fyvMZcE0AXXPUK2X70b+eSTYIsrMQXKhFIAErWsL3m4l1k2xqdSbai53m9+sI71pFt278s5YZD3MezruDze/7uk3jXA6xv1dbkRztZeCpr4HXx0+9kcSU5p4xsu89j0Skxk21D1dyXufezWL39Jj3qo+RvfF0yrFlh6CQejM3fOcHzOdX6J7WNCcHxGnNNAJ15Ky8CaPr5MrLN+wlfp/AuXtcMAInldWTrOIlv7ZLtHEfQsYKVTRngfcs3kwnlL7crx5tOtvxF7A+S0Jd5XnbOS2T7Y8si3r+DAyCaL+NzyScUkXa9sqhACXwbnME+2/VBXs0UTOlr6uc8yOftAuVk23MaHzPQPOq+So/t2due0A3DMHyCTeiGYRg+wSZ0wzAMn2ATumEYhk+YUFE0H2ZBT4sA1QTQYz7DQmKwnAWGbf/KSZwAoOoBFke2bKsl27TZvWQreYij9pp+dgzZCl9icaSwjUW/oRPYlr2ABczeTkUdBBDq4cuWZb0L6Tc4enTz14fJFtvE/SmiQHfg3KteINsZxU1k+2rdx8m2YtYWPiCATY9wArHhCkW0yvCzR+3XvW1v2MrnNlFkY0JJtrRISk0AnX/Fq2TbdttpZMuHdd+OzhgkW6pTER2F/Wb6amUKUDS+du4Oei9h/6x4kEXyxHMchVm9W498/H3x8WSbUcJRxps/xPdk7W9ZxNxzOp/MME8bSE1jwT9dxn0s6ObjTVuvJC4DsPXTvNjgjDM5Pq35mzyXtJ3mPZcDZXQcjT2hG4Zh+ASb0A3DMHyCTeiGYRg+wSZ0wzAMnzChoqgLcJUhLQWuFgGqCaC5XhYwC6M1atvZGItEwTirVhUxTo+ZirG4EVDqWGU1HUoRM+pi3O/NUa5kMhhVSjIByCml4LTv5nyBMo4h7lBeSeXtFM8oD/PYVIVYsAqEuY3aAq5EAwBvKOXXsoVKJGzRAXLHThLCiTxVGdJS4GoRoJoAOu86FqCDC+aqbbf8A0dnKpdALSPXu4jHP9bOtoV38X2armYlPsldQfh9XIGqrVsX/JdXdJHt5fXzyDb3AU4P/NZHlAMG+Jynv8DieWQv+1f7cr5+ibk8sF0n6NPowtt3kO3N13kRQEq5BvnIqH6PsbCiPaEbhmH4BJvQDcMwfIJN6IZhGD7BJnTDMAyfcESiqIjsAJAAkAOQdc41vtP2BT15zH9gwGPTaoBqKXC1CFBNAK2+WMnzCaDjC5zKNPoaq5iBx1kI6bqQhymyhkWd4l0smAxV8Xfmg/edTbaGVRxxWRbRhcDUAj7v4DCLRG1n8PmV/FVJ1ZnlfTOFrMKsvv59ZPtD5QqyzVDqfz5/36lkAwBwsC7mPMyi1ZAiwK39mNcnkmk9Re/hcKi+DTgg673+Wg1QLQWuFgGqCaC5rdvVlotaWVBPNLDfFe3h69I/n4+XquRrP9TA90WsOUG24QquCTqQVFT3A6SDLQrxoEUrOXVyulRZgcDBnggoqW2TNUp91CyPTaxbuZ9ruN/BA5WyzfP+OWUoAkqgaSA7RhV0FOOxyuUc5xxL04Yx9THfNqYU9srFMAzDJxzphO4APCEir4jINePRIcOYJJhvG1OOI33lcpZzrkVEqgE8KSKbnXPP7r/ByM1wDQBEI/x+zTAmKebbxpTjiJ7QnXMtIz87ADwE4BRlm5XOuUbnXGM4pOR4NYxJyKH7tl4j0zAmksN+QheRIgAB51xi5PfzAXznHXcKCrJxr8K861JeYaEVdNbymWvh/NpqFgCovpOL8+7+Bm+78yJOMTD7/27lfT+zgGzxXazGRxK88mLXZXzOzZ/i5QaF7XoS5FBKCY3vZls4wbahSv4Or1o7QLZsnPv91keV1T7VnIs7/gSvAOo5kUwAgGNv5PzQqdM4P3T/bF7CkM14bc4d3sqA0RyOb+diQfQsKfEeJ6fkBlcKOmv5zLVwfm01CwAU389pAnq/xb5d2KEsu5nBPhZdwMs2WqN8X8x4jle+DNbwdcr0cAqLcIW+NCSnXMP3z+IVYH9ZxHnlXYzPL9LKftyuFMYeXsd9rHmB01UMVfJfYin9sqDz/DlkS8zi7UJ8+ZGe751LXMHYEqIfySuXGgAPicjbx7nPOff4ERzPMCYL5tvGlOSwJ3Tn3HYAS8exL4YxKTDfNqYqtmzRMAzDJ9iEbhiG4RMmNB96NiroXuwVHwId/LK/5uJdZNMKOmv5zLVwfkAXQOtvZqG096rTybbtn1iwXPT+bWTrbppNtkgfx/U2/IqFo9b3sYCWVfOeA9Vr+Lz7Z3NMce3jXOm56VbO7Z6YwyJm3dPc7+h0Vm9Ki1gI7vwAn1/ZsyxqA0Dz57gocLyZfWLGn1iguv1Ld3k+fyLOYfUTyig9r7iFxe+Ok/iW0wo6a/nMtXB+QBdAG77Fvt3x33m7iifY7/qO4WtVyusC0LeAFy/ktOwLET6Z+TV6AO6637A/BJXa39N2sn/2DXLj8TZue7hS8c9tfE91NpaQLc+njNR0vUh0UQu3U7iHRd+Msvjv5NneObDvAGlARmNP6IZhGD7BJnTDMAyfYBO6YRiGT7AJ3TAMwydMqCgaTDuUjBIzKjewoLA1OpNs02ZzYWWtoLOWzxzQI0A1AbR81d94u1t5u9c2cshXjVLweLCahZrCThZqGp5g5ScX0b9vQ0kW25T05dhzPo9jaB1vGFYi1VIViqCjjO1QIQtHMSXn8/TVbWwE0Hwxi93RXvYJF+axuLn1Qs/nPZnfqW1MBIG0Q/Eur3AVSrLQVrJdUcCEb0OtoLOWzxzQI0A1AbT6JyyUNt/I2+Vi7J/J6Tz+JTt5u3gL+3HnMIusTZ1KyCSAwKmcY726lCOZe1az3wR4uJHXcqSnFWEyzuc3/WHOP59eMINshZ2Kw0PPsd5+KtvqVvP9vGamdyFGMqkXjKc2x7SVYRiGMemxCd0wDMMn2IRuGIbhE2xCNwzD8AkTKoo6EeQjXkGi41qONCx6ngWAkodYkEvFOOpRK+gM6ClwtQhQTQCdez0LpbtvUtKTdrI4FU6yKtP9KVYh3XoWF6MdZBo5ppI+t4+FlVyUxZ+gEnBWsZmNuQL+rm+/gMXKmiqO4Ez+kYtYb76JRWkAWPQtjmbta+TUsXtO5eu/fcdsz+eB4fErEn2o5KKCnkVecSzWxWGFHUrq1umr2Wd7F/G10wo6A1BT4GoRoJoAWve/WCjt/+RpZEs0cLN7zmafq1zD926qivuSKed9ASC+hq/znjKOZA4q2vKs37Mvdizn+0q7B1Kl7O+dX+RC3dFuRVA9QCp8p0xFUssLOYbLeMwK67zicCCij9do7AndMAzDJ9iEbhiG4RNsQjcMw/AJB53QReQeEekQkfX72SpE5EkR2TryU39BahiTGPNtw2+MRRRdBeAOAP++n+0GAE85524RkRtGPl9/sAPlIkB/vVckHOhjRWHhHevI1vQzrjMZCLLAFFnDAgqg1wDVUuBqEaCaAFr/XRaTtn2fBdV8DUfOXTJnI9n+9CcWorTaoQDQvUQRRTv5UpadvYdspdexUJc4hsXllnO53dIXWbzprWbbnKc5qnfHmbqr7byUo1mzhcp5H8sRhHV3e0XIDkWwOgirME6+DXAN0Z7F/LwkA5pSxqZYOxtTlfr5aTVAtRS4WgSoJoCW3Mc1SrNXs28PzuL+qJGZi/naBXfo92nRCl4JMLi9kvfv5IZazuF6n4PH8/3XcD9fl0Q9X5d0NS8CSE9n34y9pYvxaq3QnHJdy7k/6S1eMdellIFVOOgTunPuWQA9o8wXA7h35Pd7AXxkTK0ZxiTCfNvwG4f7Dr3GOfd2co492FdU1zD8gPm2MWU5YlHUOecA6O8GAIjINSKyRkTW5IaUv0EMY5JyKL6dNd82JgGHO6G3i0gtAIz8PEAIDOCcW+mca3TONQZjSjSAYUwuDsu3Q+bbxiTgcCNFHwFwJYBbRn4+PJadIokcZjzT57Ftq2UhY9u/LCVb4UssJmSVCK3iXUoxRgDxXRyRqtUA1VLgahGgmgA67585ojQ0h0XWxy9jISqWVoTOXv1cFv5UqZ3plGi851ns3P0hFssqNrL4s/D6DWTb+u0TyJarZGG69QO8MCTMQwMAmLl69CtsILCXn3aTx/Kbj9DX2j2f5b8q+VMPncPy7eAwp4bOxFmALt/M+7azO2DhXTwuQw16aujWKI+3VgNUS4GrRYBqAmjFPXwBw4Pc8SFFuK0pZVG0Z/7YIh8BINbKgmDDoyy8b7mS5xIkNBGa76uqtdzHVCVHmQ7Vcb+HK/X7dN6vec4Z3MmLCPqUCOBMhdeXXOiAfyh6GMuyxV8C+BuAhSLSLCKfxT5n/zsR2QrgvJHPhjGlMN82/MZBn9Cdc5cf4L+UhW2GMXUw3zb8hkWKGoZh+ASb0A3DMHzChKbPzUcCSNZ5I8QypYo4oggAhW0sMCnaBoaq9O+oSIKjuSJ9LAZqNUC1FLhaBKgmgGbf2km2QI6jI4eqlKg7pY4mABS28IqKXBGPT6S1n9up5vqHgz18fvFijuTLV7M4XF7GAmaiRhnDhB7lmKlkZTsQV2o0Cu+/pMybene9Ejk8YQgo4lOrcVn+cjvZei/hsU5X8zWONbNwBwAznmOxtG8B+4NWA1RLgatFgGoCaPH9HFE69AWOqu5LshAfCuqi6NzSbrK9spR9ZGgdj5mrZP8MtrMvdS5lf696nftS8Sb3saOI78lcVBcsAymeXwJsQl4LNB09B44xCNqe0A3DMHyCTeiGYRg+wSZ0wzAMn2ATumEYhk+YUFE0UyToaPQ2WbqJtzvrilfINnQCizx1MY4We/C+s9W2d13GAkfDr1gcKexk4UirAaqlwNUiQDUBdMb3OPXujps5Oi/DgWr7jplhQUgTVYvaWBCqWNRFtu4aRZQrnccNO07T2j/AglfxIr4uZT/R06Vuu5RdsHgLR9NVv8pt31a71vP5uTDXa5woJJNHtM3bfkWOx6ZjBddLrXiQRbUkb4bhCiUSEsBgDftxThHa4i0s5Gs1QLUUuFoEqCaAVt/Jvr0nwtsFBnUhccOF3HhJEV/7rhPYn6Y9zb40XMb9Dik1eZ0iumtZfCrXs7F4u57HZ3AO97H5fN5/8a1cV7f5Yu+80ZEcmypqT+iGYRg+wSZ0wzAMn2ATumEYhk+wCd0wDMMniFPSrr5blJTUucblX/DYMkUsZCTqWBjJXtBHtuIoizylVw6obTd/inNUpqr43Bue4GO2nM3CUVTJkh1UUuBqYmU2ztvN/oaSendWPTcCYO/yGbztIIu5bUodz9q/cqhaLsrf64NKqtWqV1n8GarhsQkNcV9S5br+PlTJ7ZTsViJ4FeGv+zRvKOae7/wYwzuaD7mw6HhQUF/vZv6Paz22sjeV6N+LOC1u5jmumRl+H0dMDiSVCFoAmR6+BojwNSh9jZVS7R7QaoBqKXC1CND8nyvINv12Fkp3f4OFUgBIL+aUs99ufIRsq5p5//5768iWU4as/LJmsjU1sQpduIMXYlS/yuG/Wj1SANh7Dp/Lf1vyLNnu/97fky31Ue981/SVn2GoqfWgvm1P6IZhGD7BJnTDMAyfYBO6YRiGTxhLxaJ7RKRDRNbvZ/uWiLSIyLqRfxe+u900jPHHfNvwG2OJFF0F4A4A/z7Kfrtz7geH0phkcpTStf/0Ku4Uawno7eSoq8Eoi0FlEU6hCQCF7SwSZbVUmBG2aQJoKDW2GqBaClwtAlQTQLM7d/OGACILqtnWy9F04QGOLAylOGI20sdjlouw4JWq4vFOVrFYWb6Zj5edwQITABR1cH8KulmYTpUrBWQzo8b20PX9VRgv33ZAIOPVrHqP4w5Jgs+jerciaHcrkbVp/fkrXMHXfn4NRwQ3dXJ650w5j39wB7et1QDVUuBqEaCaAFp/MwulAHDNlu1kO7GglWy/DjWSbcdiPl54vlLPdAMvKkARn0uygcX53ZXs73N+p0xYALpP5em1IsSLNpxyWQsf8N67gV4lfFfhoE/ozrlnAbA0bxhTHPNtw28cyTv0L4rI6yN/tnLZccOYuphvG1OSw53QfwpgHoATAbQBuO1AG4rINSKyRkTWpLNHL3mSYYyRw/Lt3KCeoMkwJpLDmtCdc+3OuZxzLg/gLgCnvMO2K51zjc65xkhIeQ9qGJOIw/XtYBGXjDOMieawJnQRqd3v4yUA1h9oW8OYSphvG1OZg65yEZFfAlgBYJqINAP4JoAVInIi9q0r2AHg82NqTQQu7G1Sy1ecrGWlPNTDXc0pq1RSC2rUprVVKdVrOIw3lGS1O6zkT+5ewraFP+0km1bQWctnroXza6tZACC8mvPFB5UC1dNf4Ljn9kZeqVKyk8+57BUuZrzj8lqypct4hUYmzn+JFbUpFb0BhAfYHtrCodmVe6eRrepq74qpvtihFYkeT9+WaA4Fx+712NIbeZXRZee8RLbfFx9PtuUVvEqlKKSv4Mo5vofW/YaPGTiVV3zE13CB6aIVyrIuBa2gs5bPPD3IPqetZgGAlcfMJdu22z5JtnADv+Iq6FOKW/+Zl5R9/HPPkO3n607l45Xwaqv0MN/PVbdyIXgAyA/x2P6hcwnZes7nVTJNK1Z5Pp/yBs8tGged0J1zlyvmu8d0dMOYxJhvG37DIkUNwzB8gk3ohmEYPsEmdMMwDJ8woUWikU5DdnnDePMRzp8c7VRyiKurwvj7KDjMAh8ARLtZxOyfzaJhoZJxONrHx4x2KkOn5JbPFXHIu5YjPd7C4qAWzg/oAmj2LRZmhho5x3Osk/tY1KK0U8C5s4NKhHNAGbDCDj6XdImeyrlkGx/UDbP4l6pncau9yzu2w9mxhUe/G+RzAQz0esVgmcnn8ceWRWSbUdJPtpfXc5HuaKUeYv7+WU1kC7Keh+pSDjvfU6ak1NjO+dljrTy2ryxl8Vsr6Pw1JZ+5Fs4P6ALovOteIFvrVzmdwNAyHp9hYX//jzd4JWq4meeCypO5BkPJUhaCF8R1EflvG7kGQ6SD541py1jwXDfsvYBJpy8qGI09oRuGYfgEm9ANwzB8gk3ohmEYPsEmdMMwDJ8woaJotiyGrg8f57ENNihRig0sRqTfYPE0X8CCR9sZer6YcIK3rX28hWx7zp9JtlyUBb2ys/eQLfN8GdlG538HgKI2FmC0gs5aPnNAjwDVBND4r18k29ZVJ5NtoJ4V54bH+bpklUjD6mKO2Gup4mtVvIlFVgDoOZ6j6WLTWTgs6GWB8fsn/sbz+SuxXrWNiSC8VzDzUa9wmKhnIVE6OOJ184d4DOY+wOOfLtV9+y+LlpFt2k7O5d2zmiN9g8pig2An97vhUR7boXUsqHadwLZVMRYwtXzmgB4BqgmgM37A+dS33snRnqEEP7MWt/D9XL2G223v4nuqbRaLk5uKlfzqAOb9Qsk1n2TReAe4JsTV2U97txkaW7ybPaEbhmH4BJvQDcMwfIJN6IZhGD7BJnTDMAyfIE6Jbny3iFfUuyXnfdlj61zK3ynz7uPUoZu/zsJRMMQCRfVDLBgCwFAltzN4DgshIUXoCSpZS+se54ix3R/iCLuhau5jxSLet/AOFlS1gs4A0H4ypyPVIkC7zuGOL7iKU+/mzzqRbC3nsABXuYH7M1DLAlrNyjVka7rlJLIBQOkWPYJ0NHuP4fMbve+bD96OZMfusR1wnIlNr3fz/vErHlv5FhYmmy/jFL+1D7Fg3Hq20sgBAmFdjK9LyWuKCK2MTO1ze8nWcg6L8ckZ7Meukv1r2tN8/4kS5NijFHQG9BS4qRM4AjTfx+e34Au8CGDgE6eRrXUFd6h4Gy9KyPJthrJtvG+ySn8uHqxnn61ZwmmpQz/meaPzam91tx3/vBJDTa0H9W17QjcMw/AJNqEbhmH4BJvQDcMwfMJBJ3QRqReRp0Vko4hsEJEvj9grRORJEdk68rP83e+uYYwf5tuG3xhLpGgWwHXOubUiUgzgFRF5EsBVAJ5yzt0iIjcAuAHA9e90oGAqh+Imb7Rh+ymKADOHBcLYJk5Dm1f0z0BWFxKr1nLq0MQcFkDDrJOiYjOLP4ljuI8VG1kEG+xhJau7htstiPJ3a6RPryGp1QDVUuBqEaCaABr4yzqylc9gMam/gc9lqIaFH3fysdxGWtdztOuSLmXBK13CF3ugwdt2Xg9GfSfGzbcBkOgY7VDEvESMbHtOV8YmwOMaSOljGGnlE48rNVzziqjasZzTEg8er+TeTfBUEWzna6LVCA4P8rmE53PUMaDXANVS4GoRoJoAGn+AU+/GL2dFVrbw/Vz6liKebucJomsJ388AEFBK3LZ18nxXVsdju6zWG8XeER5bvdyDPqE759qcc2tHfk8A2ARgJoCLAdw7stm9AD4yphYNY5Jgvm34jUN6hy4iswEsA/AigBrnXNvIf+0BUDOuPTOMCcR82/ADY57QRSQO4EEA1zrnPBmn3L7F7OqCdhG5RkTWiMiadDapbWIYR5Xx8O1cUnlXZxgTzJgmdBEJY5/D/8I599sRc7uI1I78fy0AtQ6Tc26lc67ROdcYCenZ4gzjaDFevh0sVGskGsaEclBRVEQEwN0ANjnnfrjffz0C4EoAt4z8fPhgx3IBQS42StxUosh2n8fqTRFnuoVTep/RioICyMZZOKp7mkXMVAW3nSvg772Wc7mNhddvIFu8mAWTdCnXixycrtQjjbCABgBlr3C0mVYDVEuBu+uDLDppAqgmJrX9iLeTChbQms/hc44pXQYAF+TrFd3JqVprsixalX13l+dz172KmPcOjKdvBzNcF3b7x3gcytbz+Q4ra2imv8DnkqzhhQEA0L6CBbPhSvZjTZjWoqAb7lee85Rwz86l3EYoyX5cflkz2Xo26ClnP/65Z8im1QDVUuC2rmB/1wTQGZdsJNvef2Tf7uBm0fkPfM6RTfqcU/4mj1lHkMNP06W8/xv3e/ud7HlcbWM0Y1nlciaAKwC8ISJvL4e4Efuc/QER+SyAnQA+MaYWDWPyYL5t+IqDTujOub9AzQIBAFCeUw1jamC+bfgNixQ1DMPwCTahG4Zh+ISJrSlaFEDnyd7VANNfYiEjMZOFh3OvYpGuPMzLIFdf/z617bc+yqcanc5LzQof5zS97ReweFr6IosbW799Atny1Yrq5Diqc8FPFYG2SsnfCWDH5UptSA5KVGuATvslj7cWAaoJoAu+zNcgNGcWH28Z12JM1Om5X5su5ZVPsXYWE4fLWGzb/ppXXE4O6amTJ4LgYBZlazs9tvAgp0XNh5X0sNMUoW0v+00gq6e6Hl7HflK2jYXSTJyf31KlbEvU871StZZ9qep17osTPr+mJvYHFOkR3T9fx3VBw81Rd8S3AAAMHklEQVR8XbUaoPkIrzTSIkA1AbT0P9i3g2nebmAG++uBIpRLN7C4n45zvd3gsCKenjEqClqfCgh7QjcMw/AJNqEbhmH4BJvQDcMwfIJN6IZhGD5hQkXRfAhITfPaSnay0DNUzbYzipvIVhXqJ9sfKleobUeqWUQpLWIlcaiQIylrqrjuYm81qxS5Shaiysu43f4BjgAdquHjJat0ITFdxiJKQImQrS7mtgdqWfTVUuBqEaCaAJp9aye3ceFMsuUP5GlV3M5wWhnbuHLOZaOEw6ASdjxB5KNBJOd7Ba9wQkmnPENJDaxcz/blfJ1i3fr51bzA/tnZyH48/eHtvN0X53J/qhWBvpKPV/GmImwqum3hDr74yQZuAwAKStgfKk/uI1t7FwutWg1QLQWuFgGqCaBatHTB+Y1k612gq6J9x3MIcD8HiaOcA1e5TqySTlnDntANwzB8gk3ohmEYPsEmdMMwDJ9gE7phGIZPmFBRNJR0qFrnFUPaG7kLc+/vIdtX6z5OtkCYBY8ZSv1CAIg/wdGHnR9g0TGmBBsm/8gFa+Y8zVFgrR9gESRRw4JJ8SLeNzTE41Cu1DIFgEyco9UKO3gsWqo4Km3hyjVk02qAailwtQhQTQCtvvN5sm3/3ulkA4Dyp8cWArf3GH72iG/0isvBgaP3fJKLCEVYVr3EkcxdH2QHK+hmQTsxl6/nUI1+fkOVXKcyr2TaTS/glLVRpe20ksp5qI4F0I4i7k/let63+lVeLLBbSe8LAOlhjvYsWdpNtrZZPD41rGGqNUC1FLhaBKgmgIaf4PsnOEv37T7FZ4MLeSFH5CUlWrrEG00uYxT87QndMAzDJ9iEbhiG4RNsQjcMw/AJB53QRaReRJ4WkY0iskFEvjxi/5aItIjIupF/F7773TWM8cN82/AbYxFFswCuc86tFZFiAK+IyJMj/3e7c+4HY20sHxYM1ngFidm/Z5Fg+00sJK6YtYVstQUcIff8fZx+EwB6TmRb2bMcsTl9dRvZNt/EYueOM3nown/jNsIJFp3KfsKCY6qcj5edodeQLGpjgSRdwu0Ub+JxbLrlJLJptSa1GqBaClwtAlQTQOd+TRmcA2yrRc7VP8kC8Yd/9JTn8+1Psz8chHHzbckDoVGBx11KtGYwxWM9bT1HTXadwAMb5KzLAIBUtWKbzscs7GRBNqPUbY+9xX4zXMk+l4uyAKqJkN1L2d/n/E7J9wyg6laOPF4Q5xrdm4pZ4E1Wcb+7lnDbWg1QLQWuFgGqCaCVd+m+nbr+DLJVlXIa4oJmHseiR72+E+zTReTRjKUEXRuAtpHfEyKyCQAvbTCMKYb5tuE3DukduojMBrAMwIsjpi+KyOsico+IKLXLDWNqYL5t+IExT+giEgfwIIBrnXP9AH4KYB6AE7HvKee2A+x3jYisEZE12SH+c8wwjjbj4tsp823j6DOmCV1Ewtjn8L9wzv0WAJxz7c65nHMuD+AuAEoOM8A5t9I51+icawzFOGjAMI4m4+bbUfNt4+gzllUuAuBuAJuccz/cz75/YctLAKwf/+4ZxruH+bbhN8S5d86zKyJnAXgOwBsA3pa6bwRwOfb9SeoA7ADw+RGR6YCUxGe6U4//vMe25bMc+l2yiVd3xDpZZc/GWK3W8jEDQPWvN5Ct+XPHH6CnXup/20K2nZeydlb/GKcsyFTyMoJtl7IWXfY624o69EK64QEei4IuXjXQczzn1M4pqQ2q1g6QzQWVYr9KQWctn7kWzr93Ie8K6Ktf3BlLyda1hNvuf7/3nFtu/AmGt7coTqEznr5dGqt1p8++ymPrPLOKtpu2htM+bP00FzJeePsubiSvh393nj+HN1VWbRTs5ZujZxE/04WUt0e1f2VjIKXke5/Dq0r2fIx9JJvU12PMnc0rWra/xak35v2C742dF7JzBzjrAKa9xuOgFXTW8plr4fwBPUMHZt7KKTD6ruBVMlnltkqMuqTN/+d2pJp3H9S3x7LK5S8AtAM9drB9DWMyY75t+A2LFDUMw/AJNqEbhmH4BJvQDcMwfMKE5kOXdBahFm9u4/jW2bRdhnUVDFfwq85sIYsbcx5mAQYAUqcdQ7Z4M4tM0V4WevoalYK0StuBvYpwFGehpniLIgTv5nYLuvVzCW1pJpsbZmUmNn0R2QZmcghxupQVtOhOFoli7XxhtILOGlo4P6ALoPL8a2SrDC4jW+gi73i3H8Ui0enSEHZf5I3BT9Zxf1yAhbYzzmTB/s3XjyObJmgDQIJrd6NwD98v7acqxcBrOWd7Osf7Du7k6xxQ6jw3n89tfGnJs2SrCLEQDwB/6FxCtkiHkgYhyXkQapZwMem2Ts4V3xHkc0nHuXaAVtBZy2euhfMDQF8rC6BlP+dFAG/9ku+BM+d4C3o/+m96qoTR2BO6YRiGT7AJ3TAMwyfYhG4YhuETbEI3DMPwCRMqiuaKIug/pd5jCytFnYdOYgEgn+HvnmgRC4FD1XpOjf7ZLAbO+BPnz3ZhbmfPqRxxiWNZCEkeyxFtEBaYql9lQadvHiteqXIlhAxA5d5pvG09594u6OXxaTlHEYRKuO2aLEcvDpfxtcrFWfjTCjpr+cwBPQJUE0ADz71Ktk/+q1cEu73g6CXICg05VG7whiVWr+Wx2XMaj3XzN1mwTy1iv9FESECP7Mwot0Hdao6uHC5jf0iV8/Xrm8/H06JRF9/KUdX3P//3ZHMHeJTsOZ/v/WnLOsm2AxyFW/tjPpeyOp7i0qU8tsFhvlaakK8VdNbymQNAP9eYVgXQOZfzIoAXv+HNpT6w90naRsOe0A3DMHyCTeiGYRg+wSZ0wzAMn2ATumEYhk+YUFE0U55H60e94th5CzfTdpu/cwLZar/eNKY21n6sXrVnMyyK3v6lu8h2cysXeN++YzbZ6u5mcSv0Na6svKSMRaLbateSbc5j/4VsUIRgAKi6mqPV2rs45fD3T/wN2f7n9z9DtoEGFnXKvsvpW7e/xqFzgTIWO+Mbufj26ILOb3PHayvINjoCFGABFAD+33Fe4Xavm1B39pApErQv916D9DEs8AUUAa3tNFYX8xElEjmrZ09Nz+d2Tp7N12/NTFY2C+tY3E9vYYE9U6EosiHuY/PFSknWFRx1XPgAR3ACQNOKVWRbN8wR01dnP022zmP5+i+r5fvvjfsXk63jDD4XF2MROVbCCxpGF3R+m9EpcAGOAAVYAAWA+pu9qXdb3dgEf3tCNwzD8Ak2oRuGYfgEm9ANwzB8gk3ohmEYPuGgNUXHtTGRTgA7AUwD0DVhDb+72LlMHmY55ziEcAIw3570TPVzGZNvT+iE/p+NiqxxzimBsVMPOxdjf/w0hnYuUw975WIYhuETbEI3DMPwCUdrQl95lNp9N7BzMfbHT2No5zLFOCrv0A3DMIzxx165GIZh+IQJn9BF5AIReVNEmkTkholu/0gQkXtEpENE1u9nqxCRJ0Vk68hPLus+yRCRehF5WkQ2isgGEfnyiH3Knctkwnz76PNe9+0JndBFJAjgTgAfBLAYwOUiwplyJi+rAFwwynYDgKeccwsAPDXyebKTBXCdc24xgNMAfGHkOkzFc5kUmG9PGt7Tvj3RT+inAGhyzm13zqUB/ArAxRPch8PGOfcsgJ5R5osB3Dvy+70APjKhnToMnHNtzrm1I78nAGwCMBNT8FwmEebbk4D3um9P9IQ+E8Du/T43j9imMjXOubaR3/cAUAqLTl5EZDaAZQBexBQ/l6OM+fYk473o2yaKjiNu35KhKbNsSETiAB4EcK1zzpNkfaqdi/HuMtX84b3q2xM9obcA2L8CRd2IbSrTLiK1ADDys+Mo92dMiEgY+xz+F865346Yp+S5TBLMtycJ72XfnugJ/WUAC0RkjohEAFwG4JEJ7sN48wiAK0d+vxLAw0exL2NCRATA3QA2Oed+uN9/TblzmUSYb08C3uu+PeGBRSJyIYD/DSAI4B7n3M0T2oEjQER+CWAF9mVuawfwTQC/A/AAgAbsy7b3CefcaHFpUiEiZwF4DsAbAPIj5hux713jlDqXyYT59tHnve7bFilqGIbhE0wUNQzD8Ak2oRuGYfgEm9ANwzB8gk3ohmEYPsEmdMMwDJ9gE7phGIZPsAndMAzDJ9iEbhiG4RP+P6nvFvs/vjKQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5c9220df60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "C1 = sp.spatial.distance.cdist(xs, xs)\n",
    "C2 = sp.spatial.distance.cdist(xt, xt)\n",
    "\n",
    "C1 /= C1.max()\n",
    "C2 /= C2.max()\n",
    "\n",
    "pl.figure()\n",
    "pl.subplot(121)\n",
    "pl.imshow(C1)\n",
    "pl.subplot(122)\n",
    "pl.imshow(C2)\n",
    "pl.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute Gromov-Wasserstein plans and distance\n",
    "---------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "    0|3.731009e-02|0.000000e+00\n",
      "    1|1.846414e-02|-1.020678e+00\n",
      "    2|1.752056e-02|-5.385587e-02\n",
      "    3|1.470479e-02|-1.914863e-01\n",
      "    4|1.371582e-02|-7.210441e-02\n",
      "    5|1.263823e-02|-8.526431e-02\n",
      "    6|1.190590e-02|-6.151022e-02\n",
      "    7|1.014664e-02|-1.733837e-01\n",
      "    8|1.011396e-02|-3.230516e-03\n",
      "    9|1.011363e-02|-3.245532e-05\n",
      "   10|1.011363e-02|-3.245682e-07\n",
      "   11|1.011363e-02|-3.245683e-09\n",
      "   12|1.011363e-02|-3.245598e-11\n",
      "It.  |Err         \n",
      "-------------------\n",
      "    0|7.847531e-02|\n",
      "   10|2.424218e-03|\n",
      "   20|4.250670e-02|\n",
      "   30|5.679949e-05|\n",
      "   40|1.219762e-08|\n",
      "   50|1.121975e-11|\n",
      "Gromov-Wasserstein distances: 0.01011363084946804\n",
      "Entropic Gromov-Wasserstein distances: 0.0063386519916289385\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEtCAYAAAAsgeXEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHstJREFUeJzt3XuYZHV95/H312EUhvtFxxkERhFdCa6oI7qIigZdYkjQjYugIqyGway4EnHVYLI0Wd1FI7i6GJ5AJICKSlQUXTUaFfGCyMCjIowGxEEuw4zcBASVy3f/OL+WM013V03/6tr9fj1PPV11zqlT3z5d59ufOnXqV5GZSJIkaW4eMewCJEmSxplhSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYJhSpqDiHh1RHxl2HVIoy4ido2IuyNi0bBr0XBExJURsf+w6+gnw9QQRMShEXFJRPw6IjaU6/81ImLYtXUrIjYrDfLZrWmvjoicZtpPhlNl9yJiIiI+2u3ymfmxzHxJP2uSuhURayPi3rJPTl5O7fK+F0bEn/ertsz8RWZulZkPbOp9I2JZRJwRETeV3+naiDgrIv5dP2rtl4j4q4j40pRpV88w7dDBVrfpSp9/YrfLZ+YfZOaFfSxp6AxTAxYRxwEfAP4OeCywFHgD8FzgkTPcZ+Re0WXm/cDFwPNbk58P/GSaaRcNsLRpRcRmw65B6rM/KaFl8nJML1Y6rH0nInYEvgssAZ4HbA08A/gm8OIZ7jOq+/lFwL6TvTwilgGLgadPmfZE7JfjKTO9DOgCbAv8GvizDsudBZwGfLEsf0C57znAL4HrgL8GHlGWPxL4DvB+4A7gWmDfMv16YANwxJQ6HrYu4FHl/nu1ln00cC/wmGnq/Bvg863bV5XHnDrtNeX6PjQB7A5gHXAq8MgyL0r9G4A7gSsm6wBeWtZzF3Aj8NbW+g8CflDW+V3g37fmrQXeDvwI+C2wWbl9Y1nXT4E/BA4EfgfcB9wN/LC1nT5car0ReBewqLXNv916rKQJxVeXWj4ExLCfc14WxqU81w+YYd6RwLeB9wG3Az8H/qjMezfwAPCb8tw/tUxP4I3l+fzzMm1f4FLgV+Xnvq3HuBD438D3y/77OWCHMm9FWd9m5fYOwD8BN5V6PjtD3e8CfkjpczMsM7nu1wO/AC4q0/8UuLLsixcCT5myrf576Qu/Lvv4UuBLpS/8K7B9a/lp11V6yaem1PMB4IPT1PlI4B7gmeX2IWUbfHPKtGumrOv6sj0vA57XmrcPsLrMWw+cUqZvDnwUuLXUeymwtMzr1M8m/4fcWuY9sdT3K+AW4JNl2YvKNv81zXPmlWV6p158QLk+AZxH8z/orrJtVw57H6reB4ddwEK60PzTvn+yqcyy3FnlCfxcmpCzeXnifY7m1dkK4N+A15fljyzr/S/AorIj/ILmH/qjgJeUJ+1WZfnZ1nUm8O5WLW8EvjxDnS8Abis17kQTzJaUnXtyWgK7luWfCTyHJtSsANYAx5Z5/7E0jO1ogtVTgGVl3jpKIwG2B55Rrj+dJnw9u/zeR5Sd9lFl/tqyc+8CbAE8maY5LS/zVwC7l+sTwEen/H7nA/8AbAk8huYfxdGtbT41TH2h1L8rTVA9cNjPOS8L40LnMHUfcFTZT/6CJshEmX8h8OdT7pPAV2mCzxbl5+3A4WX/Pazc3rG1jhuBvcr+8unJ/YmHh6n/B3yy7MuLgRfMUPf3gIkOv/fkus8pj7sF8CSaf/QvLut/G3AND71wW1vWvRTYufSQy0s/2Rz4OnBCWXbGdQG70QSkrcuyi2h61XNmqPUbwF+W66cCr6MJs+1pZ7aWfw2wY9nexwE3A5uXeRcDh5frW00+JnA08HmaPryIpuduU+Z16mf3A28qj7cF8HHgnTz0P2i/Kc+PJ7Zud9OL22HqNzQvkhfRhPDvDXsfqt4Hh13AQrqUnePmKdO+S5Pk7wWeX6adBZzTWmYRzZGTPVvTjgYuLNePBK5uzXtqebIvbU27Fdi7i3UdAPysNe87wGtn+H02LzvF04CXAx8r07/XmvbzWbbHscD55fqLaELdc5jySpQmGB492RRa008D/ueUaT+lNOeyA7+uNe+JZYc/AFg85X4TtMIUTaP9LbBFa9phwDda23xqmGo3m/OAdwz7OedlYVzKc/3u0ksmL0eVeUey8RGPJeX5+thy+0KmD1Mvat0+HPj+lGUuBo5sreOk1rw9S59ZRCtMAcuAB2kd+Znld7oGeEPr9p+W3+su4Ctl2uS6n9Ba7m+A81q3H0ET9PZvbatXt+Z/GjitdftNlKNlXazr25T+SBO4fjbL7zPBQ/3uh8AeNC+w29OOmOX+twNPK9cvAk4EdpqyzOuYclSoTO+mn/1iyn3OAU4HHjdNLVPDVDe9uB2m/nXKc+XeYe9DtRfPmRqsW4Gd2u9HZ+a+mbldmdf+e1zfur4Tzaui61rTrqN5VTVpfev6vWXdU6dt1cW6vgEsiYhnR8QKmgB2/nS/TGb+hubVzfPL5Vtl1rdb037//n9EPCkivhARN0fEncD/KvWQmV+neWX2IWBDRJweEduUu/4ZzauY6yLimxHxH8r03YDjIuKOyQvNUajlrTJ/vx0z8xqaADdRHuMTEdFetm23sp3Wtdb9DzSv6GZyc+v6PTTbWxqUl2Xmdq3LGa15v39uZuY95Wqn52e7By1n454BD+9B10+Zt5iyf7fsAtyWmbd3eGxoeuKyyRuZeUHplX/Jw88vnbHWzHywzJ+tX07XK7tZ17k0oQTgVeX2TC4C9ouIHYBHZ+bVNMFn3zJtLzbul2+NiDUR8avSf7bloe35epqjZj+JiEsj4qAy/SPAvwCfKCftvzciFtNdP2tvQ2iOwgXw/fJpvNfN8rt104vbpvbKzcf9PC3D1GBdTPPq4OAuls3W9VtoDtPv1pq2K80rpE0167qy+cTNeTQN4jDgC5l51yzru4gmND2Ph8LUt1rT2idTnkZzgvoembkNcDzNzkp57A9m5jNpXqk8iea8BjLz0sw8mGbH/2ypD5qd/91T/oEsycyPtx6zvR3JzHMzc7/y+yfwnumWK+v+Lc0rv8l1b5OZfzDLtpDG0dTn/nTTb2LjngEP70G7TJl3H02/abse2CEituuirq8BL4uIbv5PzVhr+ZT0LsytX3Za1z8D+0fE42iOxM8Wpi6mCURH0RzxJzPvLI9xFHBTZv68PM7zaMLMITRH8bajOfUjyv2uzszDaHrie4BPRcSWmXlfZp6YmXvSnON2EPBauutnU3vlzZl5VGYup3ln4O9n+QRfN714XjNMDVBm3kFzaPbvI+IVEbF1RDwiIvameR97pvtNBpx3l/vsBryF5kTDTa2hm3WdC7wSeDWzNwdowtILaRrMVWXad4D9aY5qtcPU1jQnTN5dPtr8F5MzIuJZ5WjYYppzFH4DPBgRjyzDK2ybmfeV+z9Y7nYG8IZyv4iILSPijyNi6+kKjYgnR8SLIuJRZf33tta1Hlgx2bgzcx3wFeDkiNim/J12j4gXdNge0rhZDzyhwzJfBJ4UEa8qw6K8kuZFzxday7wmIvaMiCXA39KcnL3RcAhlv/oSTQ/cPiIWR0T7079tp9CcV/WRsu9F2bf37lDrecAfR8Qfln5yHE2Q+G6H+23yujLzlzRvcf4TzSkNa2ZaUWbeS3PS+Ft46IUnNEfy38LDe+X9NOdebhYR/wOYPFJPRLwmIh5djpTdUSY/GBEvjIinlk8I3kkTaB+cSz+LiP9cQiI0bzEmG/fL9nNmk3rxfGSYGrDMfC/NjvM2mifkeprDrW9n9p39TTQh41qane9cmpPF52LWdWXmJWX+cprGN5vv0rzauiTLG+CZeQtNE9hQDmVPeivNofC7aHa+T7bmbVOm3U5zWP1WmuEjoDlfY215a/ANNCGPzFxN84ru1HK/a2je+5/Jo4CTaF4t30zzqu6vyrx/Lj9vjYjLy/XX0rydcFVZ/6dove0gjZjPx8bjTE379vw0PgC8IiJuj4gPTrdAZt5Kc5TjOJp9823AQWVfn/QRmvM9b6Y5n/K/zfB4h9P8k/8JzTmMx87wmLfQnEP5G5o+dRfNB0q2pvVCbJr7/ZTm/NT/S7Ov/wnNsBG/m+k+les6l+Y8zE4vPKH5dNxjyu8z6VtlWjtM/QvwZZrzSK+j2Qbtt+EOBK6MiLtp/n6HlrD2WJo+dSfNB3y+SfN3gU3vZ88CLimPcQHw5sy8tsybAM4ub+kdModePO9MfppDkqQ5iYgLaT7A8Y/DrkUaBo9MSZIkVTBMSZIkVfBtPkmSpAoemZIkSapgmJIkSapQNeJoRBxI87HMRcA/ZuZJsy+/JJuvLtM4Wca6jsusc8QAzWjdLZn56GFXMZ1N6WH2r409c/fZ+8JlP7MnzFu7zjSwefGLmwZTx0B017/mHKbKoGAfovk+ohuASyPigsy8auZ7bQesmutDakhWcWLHZU7076oZnTj1a0hGwqb3MPtX2+q/m70vxH9yW81bx0/MPv8NHeaPle76V83bfPvQfHnmtWUAs0/Q3dekSNIosIdJ6omaMLUzG4/IegMbf5GkJI0ye5iknuj7tzRHxCp+f2x8234/nCT1jP1LUjdqjkzdyMbfEv44pvlW7sw8PTNXZuZKWFLxcJLUUx17mP1LUjdqwtSlwB4R8fiIeCRwKM2XIUrSOLCHSeqJOb/Nl5n3R8QxNN9uvQg4MzOv7FllktRH9jBJvTLQr5OJWJ7z5aPFJ3QYLuBEThhQJdKoO/Gy5m2y8Taf+pc2jf1+IeuufzkCuiRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUoW+f9HxfOUgbePJwfekAfnyRMdFcseYdX48azT2R/uCOvHIlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXHmdKC4ngxvdd57C4tSAdOdFwk3B9VYZTGDfTIlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgXDlCRJUgUH7dTYGKUB2vSQztvdYTsl+1fvjdI288iUJElSBcOUJElSBcOUJElSBcOUJElSBcOUJElSBcOUJElSBcOUJElSBceZGiLHHdk0bg9Jm26LLpa5t+9V2L8G7+Quxrg7rkd/l6owFRFrgbuAB4D7M3NlL4qSpEGwh0nqhV4cmXphZt7Sg/VI0jDYwyRV8ZwpSZKkCrVhKoGvRMRlEbGqFwVJ0gDZwyRVq32bb7/MvDEiHgN8NSJ+kpkXtRcoDao0qW0rH06SemrWHmb/ktSNqiNTmXlj+bkBOB/YZ5plTs/Mlc2JnUtqHk6SeqpTD7N/SerGnMNURGwZEVtPXgdeAvy4V4VJUj/ZwyT1Ss3bfEuB8yNicj3nZuaXe1KVJPWfPUxST8w5TGXmtcDTeljLgtNpELdOg3p2sw5J07OHjb4vdtEDX9qxB/Z/QM7eeEoXy6zpexWjMshpL/RqQM5uODSCJElSBcOUJElSBcOUJElSBcOUJElSBcOUJElSBcOUJElSBcOUJElShdrv5lMfOYaU+q3TWGY+BzVMnceQmk8GMYaU+3y/eGRKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgoN2ykHcFjD/tprvTurQ3wDesYD2g4W1z2/RxTL39uSRPDIlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwXGmxlwvxohaWOOOSBofT+limTWzzh2XMaQ+0cV4WIeOye8yOnozhlQ3PDIlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUoeOgnRFxJnAQsCEz9yrTdgA+CawA1gKHZObt/Suze50GsYT5NUjlfPpdpH4Ytx42X+QVnXtx3JizL3DgRG+KGQPdDci5uMP8+3pRiuagmyNTZwEHTpn2DuBrmbkH8LVyW5JG0VnYwyT1UccwlZkXAbdNmXwwcHa5fjbwsh7XJUk9YQ+T1G9zPWdqaWauK9dvBpb2qB5JGgR7mKSeqT4BPTMTmPGN74hYFRGrI2I13FP7cJLUU7P1MPuXpG7MNUytj4hlAOXnhpkWzMzTM3NlZq6EJXN8OEnqqa56mP1LUjfmGqYuAI4o148APtebciRpIOxhknqmY5iKiI8DFwNPjogbIuL1wEnAiyPiauCAcluSRo49TFK/RXO6wIAeLJYnrBrY40kL2eiMuXbiZc3bZONtZPrX2nfOPn/FaV2sZOqHG6XRsiFP7rjMY+K4AVTSXf9yBHRJkqQKhilJkqQKhilJkqQKhilJkqQKhilJkqQKhilJkqQKhilJkqQKhilJkqQKmw27APXX6AzcqEHz7zqKFneYf1/nVaw4pcMC93ZbzLxwcoced5z7wVgazICcveORKUmSpAqGKUmSpAqGKUmSpAqGKUmSpAqGKUmSpAqGKUmSpAqGKUmSpAqOMzXPjcpYQ53GuxqVOqX+6mIcqY4W1jhSnTiOlGp0/t/UHY9MSZIkVTBMSZIkVTBMSZIkVTBMSZIkVTBMSZIkVTBMSZIkVTBMSZIkVTBMSZIkVXDQTg2Eg3JuzEFMJc1XnfobjE6P61xHd8N2emRKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgmFKkiSpgmFKkiSpQsdxpiLiTOAgYENm7lWmTQBHAb8six2fmV+sLWacxqaQavg8HpxB9TD7l9RYiM/zbo5MnQUcOM3092fm3uVSHaQkqU/Owh4mqY86hqnMvAi4bQC1SFLP2cMk9VvNOVPHRMSPIuLMiNi+ZxVJ0mDYwyT1xFzD1GnA7sDewDrg5JkWjIhVEbE6IlbDPXN8OEnqqa56mP1LUjfmFKYyc31mPpCZDwJnAPvMsuzpmbkyM1fCkrnWKUk9020Ps39J6sacwlRELGvdfDnw496UI0n9Zw+T1EvdDI3wcWB/YKeIuAE4Adg/IvYGElgLHN3HGiVpzuxhkvotMnNwDxbLE1YN7PHmyvFipF468bLmbbLxFvH4ZNb9fu2gSplVvqhz/4qv279UY4cO83vx4dmJHi1Tq7v+5QjokiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFQY6aOfyiJxtyE4HwpTmo3kyaOf2K5P9V8+8wGcnBlbLSDhmYvb5p3aYrxG1uOMS5+als85/Vezdq2JGgIN2SpIk9Z1hSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYJhSpIkqcJAx5mKWJ4w20hT0tydwIkdl3Ess2GYJ+NM9aB/vavDc/SvfX7OQadxke4bSBWarxxnSpIkqe8MU5IkSRUMU5IkSRUMU5IkSRUMU5IkSRUMU5IkSRUMU5IkSRUMU5IkSRU2G3YBUq8stAE5Ow1SutC2xzhYWINydhpME3ozoObCGZQzz5t9n4+buhiE+9iJ3hSjjXhkSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYJhSpIkqYLjTEljahDjSHUay2pQdWiwevN3XzjjPw1KHNJpm08MogxNo+ORqYjYJSK+ERFXRcSVEfHmMn2HiPhqRFxdfm7f/3IlqXv2L0mD0M3bfPcDx2XmnsBzgDdGxJ7AO4CvZeYewNfKbUkaJfYvSX3XMUxl5rrMvLxcvwtYA+wMHAycXRY7G3hZv4qUpLmwf0kahE06ZyoiVgBPBy4BlmbmujLrZmDpDPdZBaxqbm07tyolqZL9S1K/dP1pvojYCvg0cGxm3tmel5kJTPsNi5l5emauzMyVsKSqWEmaC/uXpH7qKkxFxGKaRvSxzPxMmbw+IpaV+cuADf0pUZLmzv4lqd+6+TRfAB8G1mTmKa1ZFwBHlOtHAJ/rfXmSNHf2L0mD0M05U88FDgeuiIgflGnHAycB50XE64HrgEP6U6IkzZn9S1LfRXO6wIAeLJbn78/lHKJOA9I5CKHUSyde1pxzNN5GpX/dmu+Zdf6O29/beSV3TPSmmJGwRYf5XWwPjZ5XTHRe5lNdLFOtu/7l18lIkiRVMExJkiRVMExJkiRVMExJkiRVMExJkiRVMExJkiRVMExJkiRV2KQvOp4vHEdKM+k0Bhn4/FH/dPP82zE6Pf8melLL+HAcqXlpIGNI9Y5HpiRJkioYpiRJkioYpiRJkioYpiRJkioYpiRJkioYpiRJkioYpiRJkioYpiRJkiosyEE71XudBhscl4Eux6VOzU8+/4ZjvvQvDY9HpiRJkioYpiRJkioYpiRJkioYpiRJkioYpiRJkioYpiRJkioYpiRJkio4zpR6wnFYeqvTuDfgNpd6ZWz2pfdNzDo7l0fHVcSrxuR3HTMemZIkSapgmJIkSapgmJIkSapgmJIkSapgmJIkSapgmJIkSapgmJIkSapgmJIkSarQcdDOiNgFOAdYCiRwemZ+ICImgKOAX5ZFj8/ML/ar0HHjoIuq4XOjN+xfmlfeOjHr7BijvtHpf+S49cBuRkC/HzguMy+PiK2ByyLiq2Xe+zPzff0rT5Kq2L8k9V3HMJWZ64B15fpdEbEG2LnfhUlSLfuXpEHYpHOmImIF8HTgkjLpmIj4UUScGRHb97g2SeoZ+5ekfuk6TEXEVsCngWMz807gNGB3YG+aV34nz3C/VRGxOiJWwz09KFmSNo39S1I/dRWmImIxTSP6WGZ+BiAz12fmA5n5IHAGsM90983M0zNzZWauhCW9qluSumL/ktRvHcNURATwYWBNZp7Smr6stdjLgR/3vjxJmjv7l6RB6ObTfM8FDgeuiIgflGnHA4dFxN40HzdeCxzdlwolae7sX5L6LjJzcA8WyxNWDezxpH6Yb+Oj9N+JlzVvk4230elfizvMv28gVQzGNl0sc2ffq1hIPtDFGIlvXlA9rrv+5QjokiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFQxTkiRJFRy0U+qxToN6wkIb2HN+DNq5cufI1W+ceX68cyH9TQGe0mH+moFUoSF43MTs82/oML8rT+1imSt68DidOGinJElS3xmmJEmSKhimJEmSKhimJEmSKhimJEmSKhimJEmSKhimJEmSKozdOFOdxvBZWOP3SONgfowzNZ/GybOPSt1ynClJkqS+M0xJkiRVMExJkiRVMExJkiRVMExJkiRVMExJkiRVMExJkiRVMExJkiRV2GzYBWwqB5ObnzoNIgj+7aVecV/qLfuXPDIlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUwTAlSZJUITJzcA8W8UvgutaknYBbBlbA3I1LnTA+tVpn741qrbtl5qOHXUStafoXjO42n8o6e2tc6oTxqXVU6+yqfw00TD3swSNWZ+bKoRXQpXGpE8anVuvsvXGqdb4Yl21unb01LnXC+NQ6LnXOxLf5JEmSKhimJEmSKgw7TJ0+5Mfv1rjUCeNTq3X23jjVOl+Myza3zt4alzphfGodlzqnNdRzpiRJksbdsI9MSZIkjbWhhamIODAifhoR10TEO4ZVRycRsTYiroiIH0TE6mHX0xYRZ0bEhoj4cWvaDhHx1Yi4uvzcfpg1lpqmq3MiIm4s2/UHEfHSYdZYatolIr4REVdFxJUR8eYyfaS26Sx1jtw2na/sX/XsX71l/xquobzNFxGLgH8DXgzcAFwKHJaZVw28mA4iYi2wMjNHbvyLiHg+cDdwTmbuVaa9F7gtM08qTX77zHz7CNY5Adydme8bZm1tEbEMWJaZl0fE1sBlwMuAIxmhbTpLnYcwYtt0PrJ/9Yb9q7fsX8M1rCNT+wDXZOa1mfk74BPAwUOqZWxl5kXAbVMmHwycXa6fTfMkHaoZ6hw5mbkuMy8v1+8C1gA7M2LbdJY6NRj2rx6wf/WW/Wu4hhWmdgaub92+gdHdmAl8JSIui4hVwy6mC0szc125fjOwdJjFdHBMRPyoHEYf+uH8tohYATwduIQR3qZT6oQR3qbziP2rf0Z2X5vGyO5r9q/B8wT0zvbLzGcAfwS8sRzyHQvZvIc7qh/XPA3YHdgbWAecPNxyHhIRWwGfBo7NzDvb80Zpm05T58huUw2N/as/RnZfs38Nx7DC1I3ALq3bjyvTRk5m3lh+bgDOpznEP8rWl/ekJ9+b3jDkeqaVmesz84HMfBA4gxHZrhGxmGYH/1hmfqZMHrltOl2do7pN5yH7V/+M3L42nVHd1+xfwzOsMHUpsEdEPD4iHgkcClwwpFpmFBFblhPkiIgtgZcAP579XkN3AXBEuX4E8Lkh1jKjyZ27eDkjsF0jIoAPA2sy85TWrJHapjPVOYrbdJ6yf/XPSO1rMxnFfc3+NVxDG7SzfOzx/wCLgDMz891DKWQWEfEEmldzAJsB545SnRHxcWB/mm/bXg+cAHwWOA/YleYb7g/JzKGePDlDnfvTHM5NYC1wdOt9/aGIiP2AbwFXAA+WycfTvJ8/Mtt0ljoPY8S26Xxl/6pn/+ot+9dwOQK6JElSBU9AlyRJqmCYkiRJqmCYkiRJqmCYkiRJqmCYkiRJqmCYkiRJqmCYkiRJqmCYkiRJqvD/AZXy1uABQvUUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5c850d1550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = ot.unif(n_samples)\n",
    "q = ot.unif(n_samples)\n",
    "\n",
    "gw0, log0 = ot.gromov.gromov_wasserstein(\n",
    "    C1, C2, p, q, 'square_loss', verbose=True, log=True)\n",
    "\n",
    "gw, log = ot.gromov.entropic_gromov_wasserstein(\n",
    "    C1, C2, p, q, 'square_loss', epsilon=5e-4, log=True, verbose=True)\n",
    "\n",
    "\n",
    "print('Gromov-Wasserstein distances: ' + str(log0['gw_dist']))\n",
    "print('Entropic Gromov-Wasserstein distances: ' + str(log['gw_dist']))\n",
    "\n",
    "\n",
    "pl.figure(1, (10, 5))\n",
    "\n",
    "pl.subplot(1, 2, 1)\n",
    "pl.imshow(gw0, cmap='jet')\n",
    "pl.title('Gromov Wasserstein')\n",
    "\n",
    "pl.subplot(1, 2, 2)\n",
    "pl.imshow(gw, cmap='jet')\n",
    "pl.title('Entropic Gromov Wasserstein')\n",
    "\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}