summaryrefslogtreecommitdiff
path: root/notebooks/plot_gromov.ipynb
blob: f565bfbb0e815a7fd61b1f5ac36df66bcdec25c4 (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
{
 "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.make_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": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXl0pFd54P27tWpt7SWptEvd7W61pFa31G724yQ4QIwTDMZjvhDMdhwIzudhGQeGOWDPwBCTkAQCST4HmGSSCQ5wQuKYLRCHzR7TLanVai2tpVtq7fuuUq3v/f6ofqtLUpVqUUkqSfd3Th9LVfd931uS/Dz32YWUEoVCoVAcPQz7vQGFQqFQ7A9KASgUCsURRSkAhUKhOKIoBaBQKBRHFKUAFAqF4oiiFIBCoVAcUZQCUCgUiiOKUgAKhUJxRFEKQKFQKI4opv3ewHbk5+fLysrK/d6G4pDS2to6K6Us2IdHq/J7xW4iol2Y1AqgsrKSlpaW/d6G4pAihLi133tQKPYT5QJSKBSKI4pSAAqFQnFEUQpAoVAojihKASgUCsURRSkAhUKhOKIoBaBIGE8+ud87UCgUsaAUgCJhPPXUfu9AoVDEglIACoVCcURRCkCxI558EoTw/4M7Xyt3kEKR/IhkHgrf3NwsVSXwwUEISOI/py0IIVqllM378OgD9FNSHECibgWhLACFQqE4oigFoEgYn/70fu9AoVDEglIAioSh/P6K7fB6vaysrOByudA0bb+3oyDJu4EqFIrDgdfrxePx4PP5WF9fZ319HaPRiNVqxWw2I0TUbmtFAlEKQKFQ7Cq68AcwGAwIIZBSomkaDocDALPZjMViwWQyKWWwhygFoFAodg2fz4fH49ki1IUQgdeklHg8nsA6XRkYjUalDHaZhMQAhBBfF0JMCyE6w7x/jxBiSQjRfvvfpxLxXIVCkbxomobb7d4g7EMhhMBoNAYEvtvtpqOjg6WlJZxOJz6fbw93fbRIlAXwN8CXgf+9zZqfSynfnKDnKRSKJEbTNFwuV0ThvxldGSwtLSGlxOl0sr6+jslkwmKxYDabMRhU7kqiSIgCkFL+TAhRmYh7KRSKg028wn8zBoMBg8GwIV4ghMBkMmG1WlW8IAHsZQzglUKIq8A48DEpZdcePluhUOwB27l9RkZGGBwcJD8/n6KiIjIyMrYV4HqXAv1eujLwer14vV4VL0gAe6UA2oAKKeWqEOI3gH8GToRaKIR4FHgUoLy8fI+2p1Aodoou/IEtwnh8fJypqSmam5tZWlri1q1brK+vU1BQQFFRESkpKRvWhxPmuosI/ArC7XYHFI6eUqpnGikisycKQEq5HPT194QQfyGEyJdSzoZY+wzwDPh7Ae3F/hQKxc7wer0sLS2Rlpa2xUc/NTXFyMgI58+fR9M08vPzyc/Px+PxMDMzQ3d3N0IIioqKKCgowGSKTixtVgbr6+s4nU4MBsMGZaAIz54oACFEETAlpZRCiLvxZx/N7cWzFQrF7qIL366uLi5evLjhvdnZWW7evElzczMmkylgIYA/999ut2O321lfX2dqaoq2tjbS09PxeDzE0qhSjw1IKQP72Rw8VlbBVhKiAIQQ3wDuAfKFEKPApwEzgJTyr4AHgQ8KIbzAOvCwTOY2pAqFIiqklLhcrpDvLSws0NfXR3NzM2anEyYnEVYrFBVtWZuamkplZSUVFRWsrKzQ0dFBa2srubm5UcULdDbXF/h8vkCxma4IVPD4DonKAnpHhPe/jD9NVKFQHBJ0H7yUEqPRuOHEvrS0RE9PD+fPn8cyMoLxj/4InE6sXi/eBx7Ad//9Ie8phODYsWOkpqbS0NDA4uIiQ0NDOJ1ObDYbRUVFWK3WqPa3WRkExwtU8NiPcpDFiGp4plDcEaiapm1o7wCwsrJCZ2cnjY2NpFitGL/8ZTCZoLQUWVyM+TvfQQwNRXyGwWCgoKCA+vp6GhsbMZlMdHV1ceXKFSYmJvB6vVHvN1Sx2crKCoODg0e62EwpgBhRc28VR53Nwh/uZO04HA46Ojo4e/YsaWlp4PUi5uYgK8t/sckEBgNiYSGq5+iYzWZKSko4f/48p06dwul00tbWRldXF3NzczHHC4xGIyaTKaAAVlZWjmSnUqUAkhxlcSiSiWDhv9l14vP5uHLlCvX19WRkZPhfNJuR5eUwPe3/3ukEKZEh4gDBbOeWSU1NpaqqigsXLlBaWsrs7CyXLl2iv7+f1dXVmD6PrgwMBkOg2Gx5eZnV1dWYA9EHEaUAomA/594qi0ORLGwW/sFC2u1243A4qK2t5dixYxuu8z32GOTkwMgIYmEB9/vfjywu3vF+hBBkZWVx1113ceHCBbKzsxkcHOTy5csMDw+HDU6Hu5fBYMBkMmEwGPB6vaytrbG8vIzD4cDr9R5KZaC6gUbBk0/eEfYHbe6tQpEI9I6doYS/x+Ohra2NlJQUcnJytl5cVITvc59DLizgtlgQm4q+EoEeLygoKMDj8TA9PU1nZydGozFQX6DXDETiKBWbKQsgCdlPi0Oh2Iwu/H0+3xbh7/V6aWtro6amZvsCLqMR8vIgygyenaDHC5qamjh58iTr6+u0trbS3d3N/Px8XPECPXi8vr7OysoKq6urhyJeoCyAGNmLubfK4lAkC3rvnVDCX/f5l5eXY7PZGBgY2MedhiYtLY2qqioqKytZXl5mcnKS/v5+8vLyKCwsjOlewcVmmqYdimIzpQBiRJ3CFUeFzY3XgoWbpmm0t7dTXFxMcQL8+buNHi/IyspC0zTm5uYYHBxkbW2N4eFhCgsL464vOMjFZkoBJDl7YXEoFJuRUrK0tMTq6ir5+flbhH9HRwd5eXmUlpbu6h52Az1ekJeXR0tLCwaDYUfxgnDFZroySOZiM6UAkhxlcSj2Gv3kv7KywuLiIgUFBRve6+rqIjMzk8rKyl3bw14ITL2CubS0lNLSUhwOB5OTk7S2tpKRkUFRURE5OTlR72Vz8NjlcuFyuTAYDBuUQTKhFIBCodiA7vbR++/rSCnp6enBYrFQXV29jzvcHdLS0qiurqaqqoqlpaUN8QK9H1G0BCsDTdOYmprC4XBQWlqaVJPNlAJQbAg6K442m33+wQqgv78fgJMnTyatSyMWpJQhP4cQguzsbLKzs/H5fMzNzXHz5k3cbjeFhYUUFhZisViifo7BYMDn8wVqKJJpstn+qyDFvqOKzRTgF/4ej2fDBC49zfHGjRs4nU5Onz59KIQ/hFcAwRiNRmw2Gw0NDTQ0NCCE4Nq1a1y9epWpqamoewhJKQMjLpOp2ExZAAqFYovwBwIWwK1bt1heXubs2bOHRvjHg8ViCcQL1tbWmJqa4tatW2RmZlJUVER2dnbYn09w3ySIXGy2tLQUc5pqPCgL4Iiiis0UOqGEP/iF1OrqKrOzs5w9ezYpfNaJJBoLIBzp6elUV1dz4cIFiouLmZyc5PLly9y4cYO1tbWYnhWq2OyNb3xjXPuKFWUBHFFUsZkC/MVcoYQ/wNzcHA6Hg9e97nWHTvjDzhSATqh4wcDAAB6PZ0O8YLMFsN39Qv0udgulABSKI4o+xD2UwJmZmWFycpKcnJx9SV3cHIA+COjxApvNhtvtZmpqio6OjsDwmczMzKju43K5oi5K2ylKAShUsdkRRNM0XC5X2JP/wMAAtbW1DEUxuOWgkggLIBwWi4WysjLKyspYW1ujr6+P+fl5VlZWIsYLXC5XTFlGO+Hw2XW7yGH1jx/Wz6UIzXbCf3Fxkd7eXv8oR4sl4afwZAoi76YCCCY9PZ2cnBxqamooKiqKGC/YSwtAKYAYUOmSioPOdm6f5eVlurq6OHfuHFar9UC6YZIVveo4JyeH06dP09TUREZGBgMDA7S0tDA6Oorb7QbiVwCLi4s8+OCDCCGuCyF6hBCvjHSNcgEpFEcETdOYmJggLy9vi/BfXV3l2rVrNDY2kpqaCsTnh9+rU3Ui2Mu9bn6W0WgMBImD4wUWi4Xe3t64XECPP/44b3zjG/n2t799SghhAdIiXZMQC0AI8XUhxLQQojPM+0II8SUhxIAQokMIcT4Rz90LVLqk4jCgn/z7+vq2CD2Hw8HVq1dpaGggPT098PrmVhCRiKQwfD4fIyMjId0e+8FeWjfbZQHp8YLm5mZqamp46aWX+OlPf8qHP/zhqO+/tLTEz372M973vvcBIKV0SykXI12XKBfQ3wDbJa6+CThx+9+jwF8m6Lm7zpNP+lMk9b8V/WulABQHBb3QKJTAczqdtLe3U1dXtyVLRQiRsIEnUkquXbuG0+lkYGCA1tZWxsbG8Hq9Idfv1cl8vyyAcKSnp/POd76T+++/n0cffTTq+w8ODlJQUMB73vMehBBXhBBfFUKkR7ouIQpASvkzYH6bJb8F/G/p52UgWwiR/E3EFYoDjt6VUm9FEIzb7aatrY1Tp06RlZW15dpYXUDh1utN5DIyMjhx4gRnz56lrq4uME2sq6uLhYWFLdfu9gl9L11A0dYBgF8pp6SkcPr06ajvr/8sP/jBDyKlPAesAR+PdN1exQBKgJGg70dvvzaxR89PCCpdUnGQCD75bxY++hzfkydPkpubG/L6RCmAGzduIKWkpqYmUHRmtVqpqKigvLyc5eVlxsfH6e/vp6CggKKiotg+aJzspQsoFmXjdrtjDgLrLSouXryov/RtolAASZcFJIR4VAjRIoRomZmZ2e/tbEC5fRQHBV34hzp56qfFqqoq8vPzw94jEQpgZGSElZUVamtrw3bezMrKCmTGpKSk0NPTw/LyMjMzM1E3W4uXZLUAYlUARUVFlJWV0dvbq7/0a0B3pOv2SgGMAWVB35fefm0LUspnpJTNUsrm4EEUCoUiOrYT/lJK2tvbKSsri9hsbKcKYGpqiomJiUAXzUgYjUaKi4s5f/486enpOBwOWltb6e3tZXl5OeEn9v3MAtoOt9tNSkpKzM/48z//c377t38bIUQH0Aj8z0jX7JUL6DngMSHEs8BFYElKeaDcPwrFQSBY+G8WOPog87KyMux2e8R7xSocg9fPz89z8+ZNmpub42olYTQaqaio4MSJE8zPzzM8PMz6+jqFhYUUFRUlrFI2GS2AeOsAGhsbaWlpAWiI9pqEKAAhxDeAe4B8IcQo8GnADCCl/Cvge8BvAAOAA3hPIp6rUCjusFn4Bws3PQvHaDRSVla2zV3iR7cAVlZW6OnpoampCbPZvON75uXlkZeXh8fjCeTLW61WiouLyc3NjbtR3WGKAcRLQhSAlPIdEd6XwIcS8SyFQrEVKSUejyes8O/q6iItLQ2Hw7Grrg+Hw0F3dzeNjY1xuTF0Qu3PbDYHgp2rq6tMTExw48YN8vLyKC4u3lDDEA3JnAV0oBSAQqHYP3Th7/P5Qgr/3t5ejEYjx48fZ35+ftdOvrqiqa+vDymM9X1FK3i326eeUqppGrOzswwMDOD1eikqKqKwsBCTKbJoS1YLwOVyRd05dKcoBaBQHGCklHi93pDCHwgIxjNnzgTe1zQt4S2evV4vS0tLYWsKdBLdX8hgMARaMLtcLiYnJ2lrayMjI4Pi4uJtu27q+9kLYrEA4g0Cx4NSAArFAUUX/sFD3IMZHBzE4XBsyMIJnvObKDRNo729nbS0NLKzsxN671jYrraguLh4i1BN5iwg5QJSKBRhkVKytrbGyMgIVVVVW4TL8PAwCwsLNDY2bngv1v4+0ezj2rVr5Ofns7gYsfXMnqDXFmRlZeHz+Zienqa7uxuDwUBxcTH5+fkYjcakdQGpGIBCoQiLfvL3eDwsLCxQXV294f3x8XGmpqY4f/78FrdDovv7XL9+nZSUFCorK7l69WrStY/WawuKi4txOBxMTk4yNDREdnY2aWkRm2UmFGUBKBSKhODz+TAajVuE+dTUFCMjIzQ1NYX08yfSBTQ4OIjX66Wurg6AzKtXSfmnf8Jgs6E99BDssJAz0cokLS2N6upqqqqqmJ+fZ2hoiPX1daxWa2B2bzKgRkIqFIqwCCEwGo1bFMDs7CyDg4M0NTWFzYJJVBB2dHSUhYUFzp07hxACw7/+KxV/+IeYrVaMUmL4l3/B++yzEKbPUCR20zev1xZIKVlYWEAIkbDagkTgcrn2LAicdL2AFApFdASf5hcWFujv7+f8+fPbFl8lwgKYnp5mbGyMxsbGgKA0fOUraGlpyNxcKCxETE9jeOGFHT1nt5FSYjKZKC0tpbm5mcrKSubn57l8+TIDAwP7NrdAuYAUCkVE9NP80tISPT09gTm+2xFPEDg4gLmwsMCNGze2uJiE2w2bT823RxwmK5sDs5mZmWRmZm6pLSguLsZms0VVW5AI1EzgQ47qKqpIBEIIfD4fnZ2dUVfexhoEDnYZra6u0t3dzblz57YoGt/b3oZxZQWxtgbz85CWhvbqV8f2gZIEvbZAn1ugt87u7u4OObcg0exlHYBSAPvAQR0urxRXcrG2tsb6+jpnz56NOqMl3jGP6+vrXL16lbNnz4YUTtoHPsD0O9+Jz25HNjXh/eu/hoqKqJ+zH0STmqnXFly4cAG73c7k5CSXL19mcHAQp9O5K/tSaaCKpOSpp5QSSBb0UY6pqalkZGREfV2sMQCDwYDb7ebKlSvU1taGf5bRyMJb34rlgx8MO2Am2YglN18IQXZ2NtnZ2Xi9XmZmZkLWFiQCZQEcEGIRhmq4vCKRDA4OcubMmZiFTjx1AFevXuX48ePk5OREvHeii8x2k3jvbzKZAnMLTp48yerq6rZzC2J9jsoCOiDE4so5qMPlleJKTmpra+NquxCLC0jTNFZWVigqKsJms0Vcn8jUzYMyFD4tLY2amhouXLhAfn4+w8PDtLa2MjIygvt2EDyWPkCgsoAUScSTT94R9kLcUWCKg0m0FoDe2dNiscQ0ozfZKoG3I5G9gLabWxCN8gxmN5r1hUNZADGSiBOxGi6vSCSxCN1oYwB9fX2YTCbS09Ojvn+iXUC7zW7tVZ9bEFxbsLq6uq+1BeFQCiBGEuHKOajuE6W4ko9Yg7rRCOmhoSGcTienTp2KyWV00BQA7L6rKTMzk+rqarKysjh27BgDAwO0trYyPj6O1+vd1WdHg3IBKaLmoCquw8jm9s7RugwMBgM+ny/s++Pj48zOznL+/PkN8wOi3dNBUgBSyj1p+aD/foLnFkxMTIScWxCvW6qyspLMzEyMRiNXr15tkVI2R3OdUgA7QJ2IFftNPGmd4U6es7OzDA8P09zcHBCMsQj1SGt9Ph/d3d2kpaVht9uTpvnabrNZqFutViorK6moqGBpaYmJiQn6+/ux2WwU3G6gF48S+I//+A/y8/MBohL+oFxAOyKZTsTJtBfF3hGqI+h2hFMYS0tL9PX1cf78+Q0tD2JVAOHQ5wYcO3YMs9lMR0cH165dY25uLuz99yINdC+yjcJlAem1BadPn+b8+fNYrVY+/elPMz09zT/+4z/i8Xh2fW8JUQBCiDcKIXqFEANCiI+HeP/dQogZIUT77X/vT8RzFXc4qNXFip2ht4OIZf1mBbC2thZoJ7H5VB6rWyfc2r6+PlJTUykrK8NutwcCpLOzs1y+fJmhoSFcLteG58aFy4UYHERMTERMWdsrBRDNc/Tags9//vNkZWXR3d0dc2zn13/912lqakII8Wi01+3YBSSEMAJfAe4FRoHLQojnpJTdm5b+o5TysZ0+T6FQ3CHW1g6b1+sVxQ0NDSHbScTiYgqnLIaHhwMtK4LdT5mZmdx11114vV6mp6fp7OzEbDZTUlIS1+lfTE1h+dznEHNzoGl4X/96vO9+952UvX0iljoAr9dLeno6T8V4ovvFL35BSUkJ09PTFBYWfkgIcV1K+bNI1yXCArgbGJBS3pRSuoFngd9KwH0VEVBFWkeXeGf8BlsAHo+HK1eucPr0aTIzM8Ou30kMYGZmhsnJSerr68Oegk0mE3a7naamJqqrq5mbm2N+fp7R0dGY+u2Yv/51xOIisqQEWVKC6d/+DcOVK2HXJ5MFoBNvJ9CSkhIAvebgO/jlckQSoQBKgJGg70dvv7aZtwkhOoQQ3xZClCXguUeeg1pdrEgc8cYAfD4fV65cobq6etvePTtRAMvLy/T399PY2Bh1llJGRgYnT54kNzcXq9VKV1cXHR0dzM7ORvycYngYqberMBhACL81EIb9jgGEwuVyxRwcX1tbY2VlJfA18OtAZzTX7lUW0L8C35BSuoQQvwv8LfCroRbe9l89ClBeXr5H21MoDiaR0jo3o1sAHR0dFBcXU1hYGHF9PEHg9fV1rl27FrJ1dLT3KigooLKyktXVVcbHx7lx4wYFBQXY7fbQHUmPH8d49SrSbgevF6REblPFvFcpq7s9EH5qaooHHngAQHexfVdK+YNork2EAhgDgk/0pbdfCyClDFbDXwU+H+5mUspngGcAmpubD05S8T6jUlKPJvG4gObn5ykuLqasLLIhHmuMQUqJx+Ohvb2dM2fOJGTwum4V+Hy+QBdOo9GI3W4nLy8vcLr2vOc9iD/6IwyjowB4H3wQ7fa84nAkmwUQTx+g6upqrl69GvzSZ6O9NhEK4DJwQghRhV/wPwz8P8ELhBDFUsqJ29/+JtCTgOcqglBun6NJrApgbMx/NqupqYlqfawWgKZptLe3U1VVFVezumA2P9doNFJUVERRURFra2uMj49z8+ZN8vPzsdvtpObm4v7MZxDz80iLBbKyIt4/GWMAe9UJFBKgAKSUXiHEY8APASPwdSlllxDivwMtUsrngP9XCPGbgBeYB9690+cqFEeZeILAejZOZmZmTH3wY1Ew4+PjASG9EyLtLz09nRMnTgSsguvXryOEwG63k5+fv69D3TcTawxgrzqBQoJiAFLK7wHf2/Tap4K+/gTwiUQ8S6FQ3CFaBTA5OcnU1BSnT5+mv78/6vvHYgHMzs4ihKCysjLq+++UYKvA4XAwPj7O4OAg+fn5FBcXb+uCSkYLwO1272mFdPKoSYViE8qtFZloFMDc3ByDg4OcO3cOk8kU10jISIyPj7O+vo7NZotfqC4vw+JiIK3NMDWF6R/+AfPXvrZtOqdOWloax48f58KFC2RmZtLb28uVK1eYmpoK+TNKxiygvRwHCaoXkCKJUSMoI2MwGLZtGbC8vExvby9NTU2YTCZ8Pl/Cu4fOzc0xPDxMeXn5xr3MzSHGxpAlJZCXF/4GPh/mL3wB8/f8TgTfq16F+d57Sf/udzGaTEiLBfPVq3g8HrS7I6e360PdbTYbDoeDiYkJhoaGyMvLw263JyQwHQuxNJ07kC4ghUKxP2xnATgcDq5du0ZjY2NAqOy0cngzq6urXL9+nebmZubn5wNTsAzPP4/xYx8LTBHyfeELaPfdF/Iepn/5F8zPP4+02cBgwPiLX1A0NYVISUGeOAGAtFgw/uQnUSmAYPSJXVVVVczOztLX14eUkuLiYjRN2zMLILi/0nbsdRBYuYAUSYWqbo6OSEFgl8tFe3s7dXV1pKenb7gu3srhUM+4evUqDQ0NWK1WhNvNse98B+P/+B+YHnvMX4xltfqF+kc/CvPzIe9j6OoCiwWMRhACmZFB6u1UzgCa5r9fnOhWQWNjI6dOncLhcDA1NcXo6OiuD2mJNQZwoLKAFIpEokZQxkYoBeD1erly5Qp33XUXWZtSIXfSOiLUM06dOuVvI+F2U/jII5iuX8fodoPTibDZ/OmYFgvC6fQ3aAvRckIrLwe32//LFgIcDtbq6zGlp8PoKMJiAYcD72/+ZtT73o7U1FSqq6vxeDxYrVYGBgbw+XwUFxdjs9kSPo5xtyuBd4JSAIeQYCGqONxsFuiapnHlyhUqKirIC+F3j7W7Z6j1Uko6OjooKysLPMPw4x9j6O9HMxohIwOcTpiehuxscLmQQvgrdEPgffvbMb78Mobr1/0WQFkZU7/92xiLi8nq6ACHA3n2LNqpU1HvO1pycnKorKzE6XQyPj5OS0sLOTk52O12MjIyEvKMvegFFC9KARxCDkvwVFU3b48QYoMC0AWzzWajuLg47DWxECoGcP36dY4dOxZoQAbA7V40COF31WRmwvIyYn0daTDg/eIXIScHQgWs09JwfelLGPr6wOdDO3kS3+AgMjcX3/33x7TfWAgWzCkpKVRXVwdm+N68eROPx4Pdbt+xVRCrBRDsstttlAJQJC2HQYntNroCkFLS09NDeno6FRUVCbv/ZhfQ0NAQXq+XU5tO49orXoHBYEC4XAEloL3+9fieftp/8r9dFRxWAZnNaGfOJGzf8WIwGMjPzyc/Px+n08nExAStra1kZWVht9vDdk3djlgtgNtTvfYEFQQ+JKjg6dFEVwA3b95E0zSOHz+e0PsHu4CmpqaYmZnhzJkzWwVaRQVzX/0q7pISsFrR7r0X79//PbK2NiD89fslC5EEc0pKClVVVVy4cIH8/HyGhobiGui+272AdoKyAA4JKnh6NDEYDKyuruL1emlsbEy4gNUVwOLiIjdv3twwL3gz3uZm+r/xjS3WQTDJNDQ+2pO5EIK8vDzy8vI2DHTX3WCRrAJVCaw4kCjrIfmZn58PTNvajf43QgjcbjddXV00NjZiNpu3XZtMAn430Ae6X7hwAZvNxtDQEC0tLYyNjYW1CmKtBFZ1AIodkajgqZoznNwsLCwwPDxMZmZmwlMXdTRNY2pqirq6OlJTU7ddm2gFkMxD4YUQ5ObmUl9fT0NDA16vl7a2Nnp6elheXt6w91gtAOUCUuwIdXI//EgpGRkZoaGhgZ6e3emurmkaAwMDZGVlbakn2G2i7lY6M4OhpQXhcqEdP45WWxt1wViiFIzFYqGiooLy8vKAUnY6nYEmdbFmASkLQLFvqGDywUAIwdmzZ0lPT49pIli0SCm5du0aeXl5EU/+wXvaUxfQ4iKm55+HxUWkz4fxxRf9VcUxkMiYiW4V1NXV0dDQgKZptLW1sby8zOrqalQ/m722AJQCUGxAzRk+WMQrdCNd09/fj9VqpaioaEdD4eN5drQYJieRmubPMkpLQ7PZ/LUEUbKb3UAtFgvl5eVcuHABi8XC5OQkLS0tjIyMbNu8TxWCKRSKqIlHgOnFXeGuHRkZYW1tjcbGxi3+7Eh7ibR2eHgYn89zeHu8AAAgAElEQVSH3W7fNqAcDdJkgmDrx+OBGNwne2GtCCEwmUyBzKjJyUna29tJT0/HbreTlZW14fdw4CaCKQ4vqhI3udlJADOcX3pmZobx8XGam5sDlcaJUgCjo6MsLi6SlZUVEIKlpaUcO3Ysrs8hy8ogPx8xMgJmM2gavte+NqZ77OU8AKPRSFlZGaWlpSwtLTE+Pk5/f38gVmA2m+OyAHw+H83NzZSUlPD888/HdK1SAIqwKLfP4SRcQ7jl5WX6+/tpbm4OZBXF0j10O2G6sLDA6Ogo586dQ0pJaWkpi4uLgYCp3W6nsLBww3MjYrXive8+xK1b/iCw3b793IFN7NdEMCEE2dnZZGdn4/F4AlbBwsJCoEFdLHzxi1/k9OnTLC8vx7w3FQM4BChBrYiFUKd6p9PJtWvXOHv27IZCpFhjDKHWOhwOuru7aWxs3CDgc3JyqKuro76+HpfLRWtrK/39/TgcjrD32oLVijx5Eq2+Pibhv5dsp2jMZjNlZWU0Nzdjs9kYGRnhjW98I52dnVHde3R0lO9+97u8//3vj2tvSgEccJ58UuXrK2Jj86leb+1cW1u7pRFZLAog1Fqv18vVq1epq6sL69u2Wq1UVVXR3NxMdnY2fX19zM7OsrCwEFPr6ljZKwsAIls0QggaGxvJysrihz/8ITU1NVHd9z//5//M5z//+biLAJUCOOAo4a+IlWALQNM02tvbqaysJCcnZ9u1kdisAKSUXL16lcrKyqjqCAwGAwUFBTQ2NpKdnc3q6iotLS0MDg7icrmi/HTRs5cKIFrcbjf5+flRpd4+//zz2Gw2mpqa4n5eQhSAEOKNQoheIcSAEOLjId63CiH+8fb7vxRCVCbiuUcZPV9fR+XrH21i7fGvdxDt7u4mNzd32/bR8SqA3t5esrKywt57O0wmE2VlZTQ1NZGSkkJnZyednZ0sLCwkLHsnWdtWRHuaf/HFF3nuueeorKzk4Ycf5oUXXuCd73xnbM+KZ4PBCCGMwFeANwG1wDuEELWblr0PWJBSHgf+FHh6p889yoRz+3z600oBHCUijYUMh75+cHAQIQRVVVXbPiOWILAuVEdHR3E6nVG7MsJhNBopLi6mqamJ8vLyQD796OhoTB05t9vzQeVzn/sco6OjDA0N8eyzz/Krv/qr/P3f/31M90iEBXA3MCClvCmldAPPAr+1ac1vAX97++tvA78mDvJPfp/ZXKwFqljrKBOPApiZmWFhYYHTp09vKwTjKTSbn59ndHSU+vr6hArYY8eOcfr0aRobGwNVttevX2dFH0YTI8nmAtoPiyQRaaAlwEjQ96PAxXBrpJReIcQSkAfMJuD5Rx6Vr3+0iVUBuFwuFhYWeMXtIS7bEasLyOv10tPTQ1NTU/QN6tbXMX3rWxgGB9Fqa/E+8MC2gtlsNlNeXk5ZWRkLCwsMDg7i8XgoKSnBZrNF7UJJVhdQPErpnnvu4Z577on5uqSrAxBCPAo8ClBeXr7Pu0l+lNtHEYsCWF1dZWFhgdraWkymyP/7xxIE9vl8LC8vc/fdd4fN+Nki3Hw+rP/lv2C8cgUsFvi3f/P383nHOyI+V++9k5ubG5jpe/nSJSpv3qToxg1M+fl43/pWZGFh+HssLyPcbmReHkTx8zhsJMIFNAaUBX1fevu1kGuEECYgC5gLdTMp5TNSymYpZXNBQUECtnewiSTclfBXRKsAXC5XYGZwtG0YorUANE2jq6uL1NTU0Bk/c3OIvj6YmNi49xs3MF67hiwsRObmIm02TP/+75gWF6Pan44+0/eVU1PYv/51nD//Oevf+hb8/u8j50KIGinJf+EF0h5/HOsTT2D91Kcg1LpDTiIUwGXghBCiSghhAR4Gntu05jngkdtfPwi8IJPV/koyVJqnIhyxBIF9Ph/t7e2cPHmSlJSUuAK726Fn/IRSLKKrC+MXvoDhf/0vTF/8Isaf/Sx4Y2y4u57OFmeHU8u3v425pISs6mqs1dVoMzPc+MY3uHXrFm63O7DO0N1N4U9/ilZUhCwtRUxOYn7mGYwvvojxhRcQk5NxPX8zsYg5TdP2PCaxY5vntk//MeCHgBH4upSySwjx34EWKeVzwNeAvxNCDADz+JWEQqFIAJEUgJSSjo4OSkpKyM/PZ2lpKaH+75GREdxuN2fOnKG1tXXjmx4Phm9+E5mfD6mp4HZj/tGP8NTXI/Py0GpqkJWViJs3ITUVsbaG78IFvAmo6jWbTFjS0zleU8OYycTVq1dJS0ujtLSUnMlJJCBuKyyZmYn5n/4J7do1vwKyWnH91/+K3OGM5WSeBwwJqgOQUn5PSnlSSlkjpfzs7dc+dVv4I6V0SinfLqU8LqW8W0p5MxHPPawclJ78ybafo0okBdDb20tGRgalpaVAbKmdkZibm2N8fJy6urrQ1oLTCW63X/gDmExoAGtr/u8tFpxf/CLe++9Hq67G+5rXYPrhDzl/8SJF99yDobc3pv143/pWxMwMzM8jJieR2dlw992UlJQEGqaNjo7SPTOD9PnQbqeSGvr6QEpkRQWyvByMRszPPruzHw6xZRo5nc49nQcMSRgEVhycAe9PPaWUQDKwnQLQXR/19fVRrY+a+XncP/kJM0NDnP2d38FoNAaKyzaQno602WB6GgoK0BYXkWYz7sxMjPrpODsbzxNPwOIiGXV1iNtNzUyDgxjvu4+1ri6I8mTsfctbkJmZGF96CZmdjfftbw/0CApuwuY+fpzR69cRLS2YUlJIN5kwBtVDyNRURJzppcHEOg1sry0ApQAUigNOOIE+PT3N9PQ0TU1NG06hsWT2hGR0FNPb345vaoozFguG55/H861vIXJzt97XYEB717sw/J//A7duQUYG2nvfi+HYMXw+Hz6fL9Aq2bhpmpeQEhwODENDaHfdFd3ehMD3a7+GduKEP75gs4VcZrFamX3Tm6j+4AdZnJzk1sAAZV//OpbpaVIyMzFMT+N5eOee6lgsAKUAFFtIthz/zVXI+t+2Skfde7YLAi8tLTEwMMCFCxe2nEB36gIy/PEf456cxJSbi8FshrExjH/1V/g++cnQF+Tl4fnQh9CcToTVikEILPhPx16vF03T/NZDdrZ/qEswbjcyRI+isLhcpPzu72J86SUQAu30adb/7u8gM3PrWiGgsJDswkKyz57FbbPh/tu/ZXVsDO8995D+67/OTsVxsscAlAJIcpJNqB4U99RRYrMCcDgcdHZ2cu7cuZBZOQaDYduxhJFY6+8nxWrFpN/baNyS3hmMlNKf4WK1brFELBYLmqb5rYG77sL1wANY//mfkR4PmEx4PvjBsKf4UJj/+q8x/uIXgZiDoasLy9NP4/7MZyJea3nta+G1r8Xs8zE9Pc1gTw8Wi4WSkhJycnLiytBRFoBCodhVghWAx+Ohvb2duro60tLSQq7fyfD24eFhzI2NlA8M+FM1pfS7Wl7zmpDrpZT4fD6EEGEFocFgCJySPX/5l3jvv5+Zl14i5cIFUu+/nyjrif33unbN/0VQBoVRfy1K9P5DxcXFrKysMDo6yo0bNzZM7ooWFQNQHFqSzT11VNEVgN7auaamZtv2y/EGgefm5piYmKD5U59CEwLDt7/td7O8731oDz20Zb0u/KWUUQtBo8nExIULTFdUcPr06cDn0sdTRrqPduoU/Pu/3zFNNQ3t9OmYP6tOZmYmp0+f3jC5S8+oygzlVtqEsgAUh5Zkc08dVXSXTmdnJzabjcJtWh/o62O1ANbW1rh+/bp/XKTViu8zn8H31FP+k3YIoay7fWJtuLa8vMzw8DBNTU2YTKaAAtADxl6vF5PJFFYReD7wAYwvv4yxvd2vnE6cwPXxLR3qY0af3FVaWsrCwgJDQ0O43W7sdjs2my1s3yNlASgUil0hOAg8OztLRkYGFRUVUV0XiwWgD3apr6/fKKC2afYWfHKPVgG43W66u7tpaGgI9CnST/26MvB6vRuyh3TLIEBqKs5/+AcM/f3g86GdPOkfGB+CeHz6wf2HXC4XY2NjtLa2kpubi91u3+J2UxaA4kATHPRVJCcLCws4nU6am5ujWh+LC0jTNBwOBw0NDRw7diyqawJB3xiEv6ZpXLt2jePHj2+NXayuYvzWtzD39GDJycHz0EN4y8vx3W4X4fP5MBqNdxSB0eh3Be0yVquV6upqKisrmZ2dpbe3FyEEpaWl5OXlBRRtLFlA4Zro7RZqJKRiW1QvouRGn51rs9miFraxuICuX7+OyWTCFkMmjtfrjUn4AwwMDJCbm0t+fv6W94zPPouhuxtpt4PXi+WrX8W8uorFYsFsNmMwGPD5fHg8nl2dIRwOg8GAzWbj3LlznDhxgvn5eS5fvhxwE8VSCXwgW0Eojh7KKth/vF4vAwMDHD9+PK6RkJG4desWXq8Xq9Ua1f31oO/8/HzUewGYmJhgfX2dysrKrW96vRiuX0eWlvpjDVlZ4PMhJiYCBWQWiwWLxRJQBG63G5/PF/Ez7kY/yvT0dE6ePElTUxMWi4XBwUFmZ2dZXFyM+LwD2wtIcbiIpheRsgz2H5PJxMWLF7FarTFPBIskjGZnZ5mcnKSurg6j0Rhxve72qaurY2pqikuXLjE6Ohpw04RjZWWF4eFhzpw5E/qkbDIh09LA4fB/r2ng9d7pLRT0mXRFoMcLdmQVuN0Yf/ADzH/2Z5j++q8Rg4MxXW40GrHb7VRXV5OVlcX4+DgtLS2MjY2FHWWpLADFrhLtqX3zyEn9a3XqTz70IGksQi6SBbC6ukpvby+NjY2BQGskBaAHfTMzMzlz5gznzp3D4/Fw6dIl+vv7WV9f33KN2+2mq6uLurq6bYfT+B5+GDE/jxgeRgwPo128iAxlLeD/eZjNZlJSUjCbzQghAoogkjIKxvjTn/rbQo+PYxgcxPTNbyLimBcgpSQtLY3a2lrOnj2L1+ulra2N3t5eVldXN6zdjxiACgIfIXbavE21gUhO4h0KHwq3201HRwcNDQ2B02gkhaG7W4L9/haLhaqqKioqKpienqazsxOLxUJ5eTnZ2dkAdHZ2UlNTQ3p6+rb7lbW1eD72McTEBGRkIGtq7vzxbYPRaAw0qdMzh4JrCra99vvfx9jdjbRa/RZHSgq+++7zTw6LgeAsIIvFQkVFBeXl5czPz3Pjxg18Ph8lJSUUFBTErACcTieve93rcLlceL1eHnzwQZ6K0TRXCkCxLcHFXqoNRHISjwUQ6kSvaRpXr17l+PHjG4qctnMZRcr4MRgMgQrapaUlRkZG6O/vx2QycezYMaKe+ldYuO1ox+3QraRgZeD1egN736IMpMTQ14fMyICMDADE0BBMTkJdXUzPDnV/IQR5eXnk5eWxvr7O+Pg4PT09vPzyy+TFoGCsVisvvPACGRkZeDweXvOa1/CmN72JV7ziFVHfQ7mADjjRjIzcyWwBdbJPbvTT7E4tACklPT095Ofnb8n4CacwomnzEExWVhZ1dXUUFxezvr7OzMwMAwMDOJ3OqPe+E/R6AqvVytTUFDk5OaGDxlIGMo5YWoLFRX+jujgUUKQ6gNTUVGpqanjFK15Beno6f/EXf8HXvva1qO4thCDjtoLyeDx4PJ6YaxuUAjjgRLL4dsufr9pAJA+JUADDw8NomhYyEyeUCyi4zUMsQmdlZYWxsTEuXrzIxYsXSU9Pp6Ojg46OjqgyZRLB0tIS09PT3HXXXRuCxl6v19+dFPC9/vXIwkK08nJkURFafT1a0LyAaIm2DiAlJYWioiL++I//mEceeSTieh2fz0djYyM2m417772XixcvxrQ/5QJSxIWyDJKHeBRAsKCdnZ1lamqK5ubmkMJ8swUQb5sHj8dDV1cX9fX1gaCv3nRtcXGRkZER+vr6KC0tpaioKOoCqljweDz09PRw9uzZwP1195CUMlBpvP7mN2PNyPDHAfLy8L75zQF3UCzE8jPS00C3C4hvxmg00t7ezuLiIg888ACdnZ3UxeCmUhbAASRet446tR9OdpIFtDnjJ9z6YAUQT5sHKSXXrl2juro6ZNA3Ozub+vp6Ghoa8Ny8yeDTTzP1d3+HS0//TABSSrq6uqiuriY1RBqp0WjEarX6awqsVlz33ovj8cdxvutd+OKcURxrL6B4s4Cys7P5lV/5FX7wgx/EdJ1SAAeQeN066tR+OInXAgjO+NluFm2wAoinzQP4K30Lx8Yo/rd/w/jd7/pnBYcg7coVTv2n/8SZL3+Z8ieegDe8gc4rV1haWtq6r4EBLA8/jPU1r8H8e78Hs7Pb7mF4eJjU1NSIVc3BNQXBgeN4agp2sxfQzMwMi4uLAKyvr/OjH/2IUzG2wNiRC0gIkQv8I1AJDAEPSSkXQqzzAXpT7mEp5W/u5LkKhcJPrIJYv0bP+Dlx4kTEtsa6wtBdJLE+c2pqCvP3v0/1t7/tf0HT/GMkv/hF2KR4LB/4AKyvI4xGhJRkDwxwvLWVPqMRl8tFaWkphYWFGJaXsbz3vbC6CunpGF98EfF7v4f72WdDdiddXFwMjMeMluDsIT3mof/b0HtoG3azG+jExASPPPJIIID90EMP8eY3vznq62HnMYCPA/8upfxDIcTHb3//ByHWrUspG3f4LEUIlFtHEStSSpxOZyD/PBK6wtALqWIR/qurqwzevMnrnn/eP9oxJcWfZtndjeHyZbRXv/rOYrcbMTkZ/GBwuchYXqahoYH1+XlW/uEfmLp+nUzAsrSEuO2akbm5GPr6YG4ONn0m3e+/nZtrO/RrdGtAjxN4vV6MRuO2dQWxxgBicQE1NDRw5cqVqNeHYqcuoN8C/vb2138LvGWH9ztU7IXLZb/dOvv9fEXs3Lp1CyCq1tFAoJo2lsEuQGBGQd2ZMwiX685pXw9aBVcHr65i/vjHkRaLfy6w0+mfOGa1ol24AEDGd7+LfWqKkpMnMTud+KanWV9Z8Ssmn89/zxDtmLu6uqipqdni94+HzS0ndKsonHso2ecB7FQBFEop9WGgk0C4RNkUIUSLEOJlIcSRURKJ6JeT7AJW9QQ6WMzMzDA9PU1KSkrUJ1MhBPPz8zGlaEop6ezspLq6moxjx/Ddc4//dO9ywcICpKSgBWWrGJ9/HjEwgPb61yNzc/Wb4Plv/w3tta8Fl8vf47+8HENqKul3342pqgrr8jK+iQlc09PMv+MdaJuEfLR+/1jRW07oQWNdSeo1BcE/h92yABJBRBeQEOLHQFGItz4Z/I2UUgohwv2FVEgpx4QQ1cALQohrUsobYZ73KPAoQHl5eaTtHXp22r5BodBZXV2lv7+f5uZmWlpaorrG5/NRXl7OyMgIly5doqioiJKSkm2DxgA3btwgb2GBovV15OIi3ieegPR0DC+/jKyqwvPhD0PRHbEixsf97qG0NLQ3vxlmZ5HV1Xg/9CH/ApMJaTT6FYjVClKi/cqv4KutxeB04iwtZcpqZeYnP6GgqoqSkhIcDkfMfv942K7lRCwWgNPpTD4FIKV8fbj3hBBTQohiKeWEEKIYmA5zj7Hb/70phPgJcA4IqQCklM8AzwA0NzcfuGYDR6FfzlH4jIeNaDN+gtEzfqxWKydOnKC6upqJiQna2to4duwYFRUVIVM6p6enMf7gBxx/+WX/H4em4Xvb2/xKIAxaXR3GF15Aer1gMCDW1vCdO3dngdGI7777MH7nOwiDAXw+fE1NaG99K6K7m4z3vY/Tq6uclpKJxx+ntb4ep9NJfX39rtQThCJUy4ngKuNI+ziI7aCfA/SytUeAf9m8QAiRI4Sw3v46H3g10L3D5yYtiai83Wn7ht1GdQtNHoLdC+FcNPqw+BMnTgRaB0QiVJsHo9FIaWkpFy9exGaz0dvby5UrV5ibmws8e3V1lVudnZz45S+RxcXIsjJkSQnG73wHZmbCPk97/evxPvAAYnISMTaG79578T3wwMY9NTfj+8AH8L3lLXjf9S60Bx4AKbH87u8iVlYgPR0sFoq/9CWOzc0FLJfW1lamp6f3pMoY7rScWF5exmg0kp6evqEjabhUUk3Tws4W3i12mgX0h8A3hRDvA24BDwEIIZqBD0gp3w+cBv4/IYSGX+H8oZTy0CqARKCariliRc/V3+xv1oOgNpst6sZrwW0eQp1ahRDk5+eTn5/vF/i3bjEwMEBxcTFjY2OcrazEIARSn8VrMiENBoTDQdg/ZYMB3/veh++d7/T/wYdxhciyMmRZ2Z0XFhb8bZr1cZUWCx6nk5zZWew1NQA4HA5GRka4ceMGxcXFlJSUYA4zJzhRuN1uBgYGAoNhgt1D+u8pVCppPHOKd8KOFICUcg74tRCvtwDvv/31S0D9Tp5zUDkKKZpH4TMeBPRisM0CZWhoCCFE1Bk/sbZ5yMjI4MyZM7hcLlpaWtA0jQmPh4zcXExTU8iCAn9qZlYWMppAbKwukKws/8l/fR1SU/G53WhuN8V33x1YkpaWxl133YXX62V8fJzW1laysrIoKyuL2iKKBV3pnjhxIuBuS0RNwW6gKoF3kUS4RJJdwCq3T3IQqhp4enqa2dlZamtrtwjz7VpCx1PpOzo6is1m49WvfjWpmZm0vOENTKem4h0eBpsN7x/8wZYpXgnBYMD9pS+BEMj1dbzLy/h+//cRtbVblppMJsrLy7l48SIFBQX09fXR1tbGzMxMQt1DY2NjpKamhpxvHG6Mpd7Nc69RzeCSHCVgFdGwWQGsrKwwMDBAc3NzWDfOZp9zvG0epqenWVpa4ty5cwghsNvtFBcXM3/33bTduoUmJeUWC/kxNo+LFu1Vr2L9xz9m4Ic/JK+2ltzz57ddH+zCWltbC7iH7HY7drs9pmZsm1lbW2NsbIzm5uaIa/WaAk3TaGlpYXl5Oe7nxotSAArFASY4QKsrAD3j5+zZs2EzfnSFoSuAmNo8uN1gNoMQrK2tcfPmTZqamjZcFzz0ZG1tjeHhYW7cuEFJSQl2uz3hwc5ba2vIu+8m9+TJmK5LT0/n1KlTeDwexsfHuXz5Mjk5OZSVlUWcVLYZTdPo6uri9OnTMX0+p9PJRz/6UX70ox8drBiAQqFIDnRXgp7xc9ddd23r3w5uCa37pCFCEHJ5GfMnPoHx5ZfBbMb50Y/SWV7OmTNntgZVp6YwjI0hs7NJr6ri9OnTeDweRkdHuXTpEvn5+ZSVlSUk731xcZGZmZkd5fubzebAuMbZ2VmuX7+OwWCgrKyMvLy8qATzzZs3sdlsHNMD0lEgpeSpp57id37nd2Ju5JYIlAJQKA4BugLo6uqisLAwpP85GN0FFBz0jRSINH/2sxhfftkf2HW7kU89xYk/+ZMtzeQMLS2Y/vRP/dk8Ph++Bx7A9/DDmM3mwJzgqakpOjo6SE1NpaKiIiahGYzb7d5Rn5/NCCEoKCigoKCA1dVVhoeHGRgYoKSkhOLi4rDuoYWFBZaWljgfwf20mZ///Of09PTwZ3/2ZzveezwoBaBQHAIMBgPj4+MYDIaoKuh1F5CenhiN8DT88pf+hm5CsK5pGIGCqSl8wYt8Pkx//ufIrCx/Xx6fD+N3voP2qlchb+/LYDBQXFxMUVERi4uLDA4O4vF4KC8vp6CgIKYZA11dXRw/fjwhfX42k5GRQW1tLR6Ph7GxMS5fvkxubi5lZWWkBfUc8ng89Pb2cvbs2ZhcOMvLy/zBH/wBzz333J7n/+soBaBQHAKcTidOp5NXvOIVUQkhIQRerxez2YzBYIjqGllUhBgexpOaitfjIc1qxbN5UIrD4W/XoNccGI3+f4uLsEkxCSHIyckhJycnkKt/8+bNqIOxt27dIj09PfrB8nFiNpuprKykvLycmZkZuru7MZlMlJWVkZubS29vLxUVFTEpISklH//4x3n88cejTtHdDY5UGqjKqFEcNoQQrKyssLi4SFVVVdRuEN1lFEvGj+eTn0QDPNPTZLrdaM3N+O69d+OijAxkSQliasr//coKGI3IkpJt763n6jc1NSGl5PLly/T19bEe3DE0CN3vf/z48aj2nggMBgOFhYU0NzdTU1PD5OQkL774ImtrazE3m/v+97/PwsIC7373u3dns1Ei9qo8Oh6am5tltE2rokFV1SqCEUK0Sikj5+slnoT9FTqdTl566SVycnLIy8ujsDBcQ96gh0vJ0NAQk5OTlJeXY7PZolIcXq+Xaz/6EbVeL6n5+f42zaFO6ZOTmL/wBcTICGRk4Hn8cWR9bLWgmqYxMzPD8PAwVquV8vJysrOzAb/fv7W1lcbGxl1x/USL0+mkra0Nm83GzMwMeXl5lJWVRdzT7Ows9913Hz/+8Y+j+n3FQdR+KKUAFIea4LYamzkMCsDlcrG4uMj8/DypqakUFxdv/+CgNg9ut5vh4WHm5uYC6Znh3C5SSjo6OrDZbBGfcfsCf0//lJQ7Ta3iZGlpiVu3bgUmgk1MTFBWVrbrrp/tkFLS1tZGVVUVubm5AYU1MjKC2WymrKyMnJycLdaVpmk88sgjvOMd7+DBBx/cre1F/QM/9C6gZG+spthdDvu8AovFQnZ2dlRzgTe3eUhJSeHkyZM0NzejaRqXL1+mv78fl8u15dqhoaGoFEwAIfyVvwnIa8/KyqKhoYG6ujrGxsZYXl5mbW1tXypndW7dusWxY8fIvT27INg9VFVVxfj4OJcuXWJsbGzDfIBvfvObpKWl8ba3vW2/tr4BZQEoDjXb/c4PgwUQfJIXQlAW3ChtE8F96kP5/TVNY2pqipGREdLT06moqCAjI4PZ2Vlu3brFuXPn9q1nDfhTLQcGBjh79iyTk5OMj4+TnZ1NeXn5hqyc3WZlZYWenp6wVdY6LpeLsbExpqamyMrKYmVlhccff5yf/vSn5OTk7OYWlQWgOLocJatPF+SRLIBo2jzo6ZkXLlygqKiIvr4+WlpauH79OnV1dfsq/N1uN9evX6e+vh6LxRLo6ZObm0tPTw/t7e0xTy2LB5/PR3d3N7W1tRF/Hlarlerqai5evMjKyoXeC64AABndSURBVArvec97yMnJYTJ47vE+c6TSQJO9sZoiMRzFdtoGgyGsSySmNg/caeOQlZXFpUuXSE9P5+rVq5SVlVFYWLjniiA43z+4clgIgc1mw2azsby8HCja2s19DgwMYLfbY+oiajAY6Ojo4Dd+4zd417vetW85/6E4UgrgMJ4AFQoIbwFE3eYhxHVdXV1UVVVRXFyM0+lkZGSEoaEh7HY7JSUlO2qaFgtDQ0NkZGRsG/Q9duwYdXV1OJ1ORkdH+eUvfxn1+MpomZubw+FwcDLGfkMDAwP8zd/8Db/4xS9i7i+02ygXkOJQcxSsPiFESAUQS5uHzQwNDZGSkhII+qakpHDixAkuXLiAECIQMHY6nQn7HKFYWFhgdnaWmtvDXSKRkpLC8ePHufvuuzGbzbS1tdHd3c3a2tqO9uF2u+nr6wvZWns7vF4vjz32GF/+8peTTvjDEbMAFEePo2L1hVIAsbR5CGZ2dpb5+XnOBc/kvY3eU7+0tJTp6Wk6OjpIS0ujoqJiS08gACYnsXzkIxh6etDq63H/yZ9AhD5FOrrfP57gsz6+sqSkhLm5OXp7exFCUF5eTm5ubszWUE9PDzU1NTHP7P3Sl77Eq171Kl796lfHdN1eoRSAQnEI2KwAggeRxyLsHA4HAwMDnD9/fluhazAYKCoqorCwMJCdI6WkoqLijoB1uUh5wxsQo6P+nkAjI6Rcv47zpZdCF5AFEc7vHyvhxleWlpZSVFQUlT9+fHwcs9kcc7VvZ2cnzz33HD//+c/j3f6uoxSAQnEICFYA8Q528fl8XLt2jdra2qj95kIIcnNzyc3N3SBgy8rKsI+NIWZmQNPAYACfDzEygrhxA3nXXdveNxq/f6zo4yvdbjcjIyNcunQJm81GaWlp2JO93qPowoULMT3L5XLx2GOP8cwzz8RsNewlKgZwwDkqLg7F9ugKYKdB37KysrhbM+sCtrGxEYfDwbW+PjSvd+Minw8iKJeFhQXm5uai9vvHisVioaamhosXL5Kamkp7ezv9L7zA+osvwupqYF28A14APve5z/GWt7yFxsbGRG8/oSgFcMA57JWuisjoQWDd7RNP0PfWrVtYLBbsdvuO92O1Wjl+/Di1Dz2Es74en8GA9HqRZjO+e+5BVlaGvVb3++9F3YHBYMBut/PqH/6Q+ocf5tiDD2I6dYql//gPpJQMDg4G0mFj4dKlS7z88ss88cQTu7TzxLGjn7AQ4u1CiC4hhCaECFtRKYR4oxCiVwgxIIT4+E6eqVAothLc3z/WsYJzc3PMzc3FnN4YCaPZjPje9/B98pM47ruPm+9+N22f+ATLKysh10sp6ezs5MSJEwmZFBYNhv/7fzH/5V8iNA2T14t5bY3cRx/lpZdeCvQcioW1tTU+8pGP8NWvfnXP0mR3wk5VbCfwVuBn4RYIIYzAV4A3AbXAO4QQtTt87pHmKFW6KqLDaDQG5vPG0iNnfX2dvr6+3TtxW634PvIRDN/4BsVPP01JZSU3b96ktbWV2dnZDZW7Q0NDZGZmRpxmlkhEX9/tL27/z2Q0Ypmdxahp2Gw2Wlpaok53lVLyqU99ive+970JV6a7xY5UlJSyByL6Gu8GBqSUN2+vfRb4LaB7J88+yhzFSldFePRsn4sXLzI5OUlbWxvZ2dkRh5QEB333IlAZPABmbW2NW7ducePGDcrKyrBarczNzcU8UnGnyBMnbn8h/f8z+Xy48/Mpr6nBbrcH+iNFM77yJz/5CTdu3OArX/nKHn6CnbEXMYASYCTo+9HbrykUih3icDj40Ic+RHd3NwaDgZKSkkCPnM7OTjo7O1kJ4XKRUtLd3U1JSUnMPu5EkJ6eTm1tLefOnWN1dZX29naysrI2dM7cC7RXvQrPBz/on1pmtaJlZtL72c8GCuCC+yOVlpYyODhIS0sL09PTG6yXxcVFPvGJT/C1r31tX3smxUpEC0AI8WOgKMRbn5RS/kuiNySEeBR4FIhqtulR5yhUuirCk5qaykMPPcQnP/lJLBYLH/7wh3nlK1+JzWajoKCAhYUF+vv7EUJQWVlJdnY2QgiGh4cxm82URJjUtduYzWZWV1epr6/H5XLR2tpKTk4O5eXlezbsxfupT+F773vxjI3RtrZG42tes8Wrsd34yoKCAp544gk+9rGPxRwz2G8S0g5aCPET4GNSyi29m4UQrwSelFK+4fb3nwCQUn4u0n0T3Q5aoQjmMLSDDtxQSlpaWnj66aeZnp7m8ccf5w1veEPgNLqyssLQ0BBOp5O8vDzm5uZoamr6/9u78+Cqy3OB49/HhARMwiLmIITkBJQtJTSQIG2vjWRAXMq0gFi3jrm4UCvK4sVeuXdu1Tp0SsQSpIhaaua6F2/Fa1VAHdyghhCBiCBQQk5ONuBAEiH7Wd77R5ZyJZCc5Jyck5znM8NMlpPf+/yc8X1/7/s+v/cJ+NNqUVERLpeLMS1LMcYYHA4HxcXFREZGYrVae2SGYoxh7969WK1Whn63zvEFOJ1OSkpK2gq7bN26lcSLZDj1oKA6Dno3MEZERolIBHAb8E4PtKtUyBARpk6dyptvvsmf/vQntm7dSkZGBq+//jpOp5OYmBiSk5O56qqrsNvtuFwuKioqOiwi40+VlZXn5fu3nvA5depUEhISsNlsfPnllzgcDr8e9dxaA6GznT80z16io6Pp168fy5cv529/+5vf4vOXbs0ARGQusA6IBaqBfcaY60VkBLDRGHNTy+duArKBMOBFY8zKzlxfZwDKn/rSDKA95eXlrFmzhm3btpGZmcm8efP44osvyMjIYMCAAdjtdhwOR4+f7gn/rOs7efLkDlM+6+rqKC4u5ttvvyU+Pr7TRzh0Vk1NDQcOHCAtLc2r63o8Hu68804WLFjAnDlzfBaPD2hNYKU60tcHgFZVVVWsX7+e559/nmuuuYannnqqrZShy+WirKyM8vJyYmNjSUhI8NnxyRfSutySkJDgVcpnU1MTpaWlnDhxgmHDhjFy5Mhux+rxeMjPz2fChAntH2Z3Ea+++io7d+4kJyfH63cv2pOYmEhMTAxhYWGEh4fTjb6v08EE/5sKSqluGTJkCGlpadhsNqZMmcJPfvITrr32Wh566CHi4uKwWq3Ex8e3pZAOGjQIq9XqtzKLRUVFDBw40Ot8/4iICEaPHo3VaqWioqIt3bU7JSELCwsZNmyY151/SUkJf/zjH/n000990vm3+vjjj3v2PQidAahQFSozAPhnYZjw8HBcLhebNm1i7dq1jBs3jiVLljB+/HhEBGMMp06dwmazERkZSWJiYpfPBmpPZWUlx44d6/C00c5o3TC22+1ERER4vWFcWVmJzWZj8uTJXnXiHo+HOXPm8OijjzJz5syuhN6uxMRE8vPzfTEA6BKQUh0JpQGgPR6Ph23btvHUU08RHR3Nww8/3FbwBZqXjmw2G8YYEhMTGTJkSLeedhsbG9mzZ0+n1v299e2331JcXExTUxMJCQnExsZeNFan00l+fn6XYnn++ecpKipi7dq1Pn36HzVqVNt/41/+8pcsXLiwq5fSAUCpjoT6ANDKGENubi5ZWVlUV1ezdOlSZsyY0faEXlNTg81mo66uDqvVisVi8brj60qaZVfU1dVht9uprq5m5MiRDB8+vN2N3f3792OxWBg2bJhX1z9y5Ah33303O3bs8PkSWVlZGXFxcZw8eZLrrruOdevWkZ6e3pVLBVUaqFIqiIkIP/zhD3nrrbdYv349b7/9NhkZGWzatAmn00l0dDQTJ04kOTmZ6upqdu3aRWlpqVdv7bau+/uz8we49NJLGT9+PKmpqTidTvLy8igsLKSpqantMxUVFYiI152/y+Vi0aJFbNiwwS/7I60v5VksFubOnUteXp7P2/guHQCUUkDzQJCUlEROTg5vvfUWBQUFpKen88ILL1BfX8+AAQMYN24cU6ZMobGxkby8PIqKijo8fK6yspLKykpGjx7dQ3fSnKM/atQorr76avr378+ePXv45ptvqKyspLi4mPHjx3t9zTVr1pCRkcG0adN8Hm9tbW3bkR21tbV88MEHTJw40eftfJdmASmlzhMfH8+aNWs4ffo069evZ/r06cyfP5/77ruPwYMHc+WVV5KYmEhZWVnbxmVCQsJ5h8o1NjZy+PDhLtX19YWwsDDi4uIYMWIEp06doqCggOjoaGpqahg8eHCnr1NQUMCWLVv47LMLHnzcLSdOnGDu3LlA80zjjjvu4IYbbvBLW+fSPQAVsnQPoPPq6urYuHEjGzduZMaMGTz44INtB6a1nphpt9uJiYnBarUSFRXVY+v+nVVUVITb7SY2Npbi4mIaGxuxWq0dbhg3NDQwa9YscnJySE5O7sGIu0w3gZXqiA4A3nM6nbzxxhs888wzJCcns3jxYsaMGdOWQnr69GmKi4sJDw+nX79+REREcNVVVwU6bM6cOcOhQ4dIS0trm4nU19djt9upqqpqmyW0t2H8m9/8BovF0isqfLXQAUCpjugA0HUej4f33nuP1atXM3ToUJYtW8aUKVPanqRLSkooLCwkJiaGxMRELrvsMp+mTHrD7Xaze/dukpOTiYqKOu/3TqeT0tJSjh8/jsViIT4+vu0N47///e88+eSTbN++3afHT/iZDgBKdUQHgO4zxrBjxw6ysrKoq6tj2bJljB07liNHjvCjH/0It9uNzWajtraWhIQELBZLj+8FHDp0iOjoaEaOHHnRz3k8Ho4fP05JSQnR0dHU1dWxdOlS3nzzzaCYxXhB00C7S8srKtUxEeHHP/4x77zzDtnZ2bz++uvMnDmTvLw8wsPDiYqK4nvf+x6TJk3izJkz7Nq1i5KSkh4r/OJwOGhoaOhU3YPWIvGtmUP3338/Ho+HqqqqHog0MHQAuIAnngh0BEr1HiJCcnIyP/jBD5g7dy5VVVWkp6fz4osv0tDQQP/+/Rk7dixpaWm4XC7y8vK8rl/sraamJo4ePUpSUpJXy08iQlFREaNGjSInJ4f6+nq/xRhougR0AVprt+/TJSDfq62tpX///oSFheFwOFi3bh2bN2/m1ltv5Z577mk7q8ftdlNeXk5ZWRlDhgzBarX69HgIYwwFBQXExcURGxvr1d9WVVVx4403smXLFp9VTHO73aSlpREXF8e7777rk2tehC4BdcXjjzd3/K0PC61f63KQUp0TFRXVtlkaGxvLb3/7W7744guioqK4/vrreeyxxzh+/DhhYWHEx8czbdo0Bg0axFdffcWBAweoqanxSRxlZWVERkZ63fkbY1i+fDmPPvqoT8tlrl27lgkTJvjser6iA8A5Hn+8+am/9cm/9WsdAJTquujoaJYtW0Z+fj5JSUnccsstLF68mMLCQkSEK664gqlTp3LFFVdw+PBh9u3bR3V1dZfbq62tpbS0lLFjx3r9t2+//TZut5vbb7+9y+1/V2lpKe+99x733nuvz67pKzoAKKV6REREBAsWLCAvL4/Zs2fzq1/9iszMTAoKCgAYOnQoqampjB49GrvdTn5+vtelID0eDwcPHiQpKcnrtM3jx4+zatUqnn32WZ+mrC5dupSsrKyA119uT/BFFCQeeyzQESjVN4WFhTFnzhw+//xzFi1axJNPPsm8efP47LPP8Hg8DBw4kEmTJpGUlITD4SAvL4/y8vJO1S8+duwYl19+udc1DDweD4sXL2blypU+Lcjy7rvvYrFYSE1N9dk1fUk3gVXI0k3g4NC6YZuVlYXNZuOhhx5i9uzZbU/wjY2N2O12Tp06RVxcHHFxce0+3VdXV3P06FFSU1O9foJ/6aWX2L17Nxs3bvTp0/+KFSt4+eWXCQ8Pp6GhgTNnzjBv3jxeeeUVn7XRDn0RTKmO9MUBoKSkhLvuuosTJ04gIixcuJAlS5b4qzmfKyws5OmnnyY3N5eFCxdy6623th0w53K5KC0tpaKi4rw3dl0uF/n5+Xz/+99nwIABXrVpt9u57bbb+Pzzz70uDemNTz75hNWrV/edLCARuUVEDoiIR0Qu+D+SiNhEZL+I7BMR7dGV8pPw8HCefvppDh48SG5uLuvXr+fgwYOBDqvTrrzySp599lm2bNmC3W4nPT2dZ555hrNnzxIeHk5iYiLTpk0jMjKSPXv2cOjQIerr6zl8+DBWq9Xrzt/tdvPAAw+wZs0av3b+waq7ewBfA/OAzpyRmmGMSQnQE5dSIWH48OFMmTIFgJiYGCZMmEBZWVmAo/LesGHD+N3vfsfOnTsJCwtj5syZPPHEE5w8eZJLLrmEkSNHMm3aNIYMGcKePXuorKwkOjra63aee+45UlJSmD59uu9v4jumT5/eE0//XunWAGCM+cYYc9hXwSilfMdms7F3716/FDDpKQMHDuSRRx4hPz+fMWPGMHfuXB5++GGKi4sREQYNGsQll1zCuHHj+Mc//sHevXupqqrqVObQoUOH+Mtf/sLKlSsDdlBdoPVUQRgDfCAiBnjeGPNCD7WrVEiqqanh5ptvJjs72+uMmGAUGRnJvffey4IFC9i8eTP33HMP8fHxnD17luzsbCwWCxaLhbNnz2Kz2Th69OhFz/p3Op08+OCDPPfcc14vG/UlHc4AROQjEfm6nX8/86Kda4wxU4AbgUUicsFKxyKyUETyRSTf4XB40UTfoi+fqa5yOp3cfPPN3HnnncybNy/Q4fhUWFgY8+fPZ8eOHQwdOpSioiKWLFnCjh07MMYQExNDcnIyEydOpLKykl27dlFWVnZeCunq1auZNWsWaWmhvSLtkywgEfkEWG6M6XCDV0QeB2qMMas7+mwoZwHpWUT+1xezgIwxZGZmctlll5Gdne2vZoJCVlYWixcv5sCBA6xatYry8nKWLFnCDTfc0JYm2tTUhN1ux+FwMGLECOLi4ti/fz+PPPIIn376Kf369QvwXfhF8JwFJCJRIhLT+jUwi+bNY9UNOkNQ7dm5cycvv/wy27dvJyUlhZSUFN5///1Ah+UXv/71r+nfvz+pqals2rSJnJwcPvroIzIyMnjllVdoampqq0g2depUAO6//37uu+8+fv/73/fVzt8r3ZoBiMhcYB0QC1QD+4wx14vICGCjMeYmERkNbG75k3DgNWPMys5cP9RmAI8/3v4x1I89dn6HrzOE7uuLMwAFFRUVZGdns2XLFu666y4yMzPbahSvWLGCU6dOUVZWxrZt23x6AmkQ0RfBeruOOngdALpPB4C+rbq6mg0bNvDaa6/x05/+lJSUFDZs2MCHH37Y7fKODQ0NpKen09jYiMvlYv78+TwRPEVEgmcJSPmOHletVOcNHjyYFStWsHv3buLj43nggQfYsGGDT2r7RkZGsn37dgoKCti3bx9bt24lNzfXB1H3rF45AIRCh9feYXR6XLVS3mst73jy5EnGjRvnk2uKSNuLZ06nE6fT2SvfJeiVA0DwzLT8Rzt1pXzL1x202+0mJSUFi8XCdddd1ytfuOuVA4DS46pVcHC73UyePJnZs2cHOpQeFxYWxr59+ygtLSUvL4+vv+59yY29ZgDQ9e//L1TvWwWXYC112JMGDx5MRkYGW7duDXQoXutVA4CufysVPIK51KG/ORyOtrKV9fX1fPjhh4wfPz7AUXmvp84CUkr1Ma2lDs+ePRvoUHpcRUUFmZmZuN1uPB4PP//5z3vlMlivHAB0/VupwDq31OEnn3wS6HB63KRJk9i7d2+gw+i2XrMEdC5d9lEqsHbu3Mk777xDYmIit912G9u3b+cXv/hFoMNSXtI3gVXI0jeBfaMHSx2qztE3gZVSSl1cr9wDUEoFj+nTp/dISUXlezoDUEqpEKUDgFJKhaig3gQWEQdQ7OdmLgdO+bmNQOnL9wbdvz+rMSbWV8Eo1dsE9QDQE0QkP0CZIH7Xl+8N+v79qX8SkcHARmAizVlUdxtjvvBDO/HAS8CwlnZeMMas9XU7wUI3gZVSvcFaYKsxZr6IRACX+qkdF/Bvxpg9LaVsvxSRD40xB/3UXkDpHoBSKqiJyCAgHfgzgDGmyRhT7Y+2jDEVxpg9LV+fBb4B4vzRVjDQAQBeCHQAftSX7w36/v2pZqMAB5AjIntFZKOIRPm7URFJBCYDu/zdVqCE/B6AUiq4iUgakAv8izFml4isBc4YY/7Lj21GA58CK40xb/mrnUDTGYBSKtiVAqXGmNYn8f8BpvirMRHpB/wVeLUvd/6gAwAAIvKUiBwSka9EZHNLxkGvJiI3iMhhETkqIo8GOh5fEZF4EflYRA6KyAERWRLomJR/GWOOAyUi0lrQdwbgl01Zaa4b+WfgG2PMH/zRRjDRJSBARGYB240xLhFZBWCM+fcAh9VlIhIGHAGuo/npaTdwe1/IZBCR4cDwc7M0gDl94d7UhYlICs1poBHAMWCBMabKD+1cA3wO7Ac8LT/+D2PM+75uKxhoGihgjPngnG9zgfmBisVHrgaOGmOOAYjIG8DP8NNTU08yxlQAFS1fnxWR1iyNXn9v6sKMMfsAv7/zYYzZgRenafZ2ugR0vruBLYEOopvigJJzvi+lD6ayhUKWhlL+FDIzABH5CLiinV/9pzHmf1s+8580vwjyak/GprzXkqXxV2CpMeZMoONRqjcKmQHAGDPzYr8XkX8FZgMzTO/fGCkD4s/5fmTLz/qEUMrSUMqfdBOY5owZ4A/AtcYYR6Dj6S4RCad5E3gGzR3/buAOY8yBgAbmAy1ZGv8NVBpjlgY6HqV6Mx0AABE5CkQCp1t+lGuMuT+AIXWbiNwEZANhwIvGmJUBDsknQi1LQyl/0gFAKaVClGYBKaVUiNIBQCmlQpQOAEopFaJ0AFBKqRClA4BSSoUoHQCUUipE6QCglFIhSgcApZQKUf8Hlcqa7NPKroIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXlwXPWV77+n1WpJrV2WJWuxvAqDDcY2BjtgdgwkZGJnSCYQKoHAi/OyvEnGvJchVGV5k8pkqQQymUmlwg4PAiFAgACBuBKHBLDB2Bjvqyxb+2LJLbXUUm+/94fFlFvfY9y2hZbL+VS5JB3fe3+/+7vn/vr2/f7OOeKcg2EYhjHx8Y11BwzDMIyRwSZ0wzAMj2ATumEYhkewCd0wDMMj2IRuGIbhEWxCNwzD8Ag2oRuGYXgEm9ANwzA8wmlN6CJyrYjsFpF9InLHSHXKMMYa821jIiKnGikqIhkA9gBYDqARwAYANzrndoxc9wxj9DHfNiYq/tPY9wIA+5xzdQAgIk8AWAHguE4fkCyXjdwUW6Ikl7aL5/GHTCDEx0tkCdty9A+oQBfbooW8vy/O2/n7+ZixUqWPrcp2BRl8wKRiCirH6+btAMDXF2WjKGORm6nszNtlHInwvoU5ZMsu5+36unm7QDf3L5Eb4L4A8MX5vLXrmlQ8NbM3dSAHBo8gGuvjnU+eEfFtVxik7Zwy/rE8Pp5PucSZfYrjAIhn8xftwKRBsiWb2B9iebyvP6Jck0xlWBWTL8Y27dq547wbCPTwDZjI4QNk9A6QLVqazV1U7men9Ee7x5PKOWcM8DWIB/WT8SXYJgluxxflY8byUueNaG8X4pET+/bpTOhVABqO+bsRwJL32yEbuVgiV6bYQtctpe3alvEJTvsDD8SRmeygoQXKnQBg5mO8/6FreZIJdPOYlW3iY7bcxjfM9B9zvxuuLiSbn+dF9C7g49U8qTtKcEM92SSLz6V3cRXZtJu/8MXtZAtdM49s81ZvJdv6p84lW82TDWTrXsp9AYBgG49taGYW2SKT+bpUre1N+fvNbb9W2zgFRsS3By69gLaLZ/N5tFzCx8s7wA8C5RsVxwHQNYcnsqqbDpAt8p0KsjVfxB/IkzfzrByu4qlCm+TzWngWi5Swz8Xy9Lmp6uV2soXmTyJb4dp9ZGu4eQ7Zsrr4vo+Uafc4n3N/GZ9z0Z4+srWfp3wiA8g5zPNBoIfHJ6cpTLbWZcUpf+/97V1qG8M5nQk9LURkFYBVAJANfmIxjImK+bYx3jgdUbQJwNRj/q4esqXgnLvHObfYObc4E/zkZRjjEPNtY0JyOk/oGwDUisgMHHX2GwB89v12SJTk0iuWwkfX03aR0gvJFnxrL9kO/sMMsuXUHeddbYy/rk56l7+OlWzqJJvL5K+/ea8Uk00S/KJ/ylv8KsVl8Fe+YBu/Pgoe1F+i9y3h8x4s5M/m7C7+elf8Ds1L6L6OX6+EK/l4r66dT7by/dxG6PxKshVt4q/SALD3i+VkS5Txa5j8d3jC7Bj2VTdeN2KrcE/at11hkF6xZP/hLdoudBO/YkQ+f93vm8o+0p7k1yMAUPFaD9l6P83j1VvLr2Zy2vkeOHw2+2K0gLcrfZevfcsy7nfuIbb1V+haV/8svq+Sft6/59LZZMvg1+romcU2SXLbPTU8FQ4Wc7slW/ilfGiO8rIcQMUvFJ/vZDGv+fN8/4XmpraTeD69xSunPKE75+Ii8jUArwDIAPCAc45fxhrGBMN825ionNY7dOfcSwBeGqG+GMa4wXzbmIhYpKhhGIZHsAndMAzDI3zgyxaPJZ7naI25JoBO+fkbZGv437ydJFiMmH4Pi6cA0Ho9iygVL/J66barp5LNKbFB8SALJvUri8g2UMmCl2Rxv/0tSsBJbgk3DKDrHBZIpITF18mvsDDWf1k12TqWsdCTv4s/6+NTWKzMPcBK1IFPFZCt7XwWPwFg0tkdZPvKzL+S7cHHV5Kt+ZJU8S6p6+GjgvMJrTHXBNDCx3gRQMf5vF1mD49/4jgLaXpn8jro/n6+puUbWTwdqOTllnmsm8MX4zXVPTN4wLM62Y8zlaCdQEhfhx6awYJs5bP1ZOtYPo1sOR3cx6K9fK8lA9y2FvDVzWsAEKngQMj8/coEAaD7AvZ5/0AZt60MRfG7qcfsiKQXL2dP6IZhGB7BJnTDMAyPYBO6YRiGR7AJ3TAMwyOMqigaCHGSLS0CVBNAK3/KQmndTz5CtgNfrlXbnv4cR122X8UCaMkujijNbOLoroZ/ZHFx5n0HyRZewEmpRFFZe2rIhPI1jWwEULKTkxVFi1hMygyzUBqoYxEyt42TNvXUsJDlItzvvhksyFX/mcXTrMYjZAOA+k+zcPS9/deTrXAWt128I1UEa9ZzV40KsTwlyZYSAaoJoLO/wUJp4rJFZGtexpGeABD9PPtn7oMccVm/ksdQy46oZkxUEndmdyqJr+awz5Vv4H1LtilpEAHsu4GjYQdum862GdxOxSvcyeaLeYor4LxlanbLmy9+lWyPH76MbMW79SyY0Rv5urQ28sKJjH7e3x9OfdbWxl/DntANwzA8gk3ohmEYHsEmdMMwDI9gE7phGIZHGFVRNJElVGVIS4GrRYBqAujMb64j276fK+lJAbRexCJRaB6303MtCxSZm1gAve4GFmmfL2Axt/byOrJF4qxwrCjdT7Yniy4jG6CXtoor9RVyG/ny9lzDQrAWtTe8GhAAvPqde8l2QfNqsuU18HPCgVu4chMATF6jCGsz+AQz+/mYwbZU9c4XO7X6uCOBL8pVhrQUuGoEqCKAZvx1E9kmlehFk9rOZrHUdy63k9vI4xNmd0DRXt5uQEklq0VX5m7ncFZ/P1f58e+s54YB5NdzKlktklL2cTvRXO53sFXrN29XcIhF1ucOnkO27MPcl8I9XHEIAPZ084KBzCMsTBdw8SVK8ZvB3VOxJ3TDMAyPYBO6YRiGR7AJ3TAMwyPYhG4YhuERTksUFZF6AL0AEgDizrnF77d9IschtCA1JEurAaqlwNUiQDUBVIu6A4C2f2bBcs69LPw1XcGRXD4lqO35Z/l48RwWW3Y1cyRkIIsP+NCeZWQr1LUWhObx/oFiTmPrq+dUn3mcMRhHlvK+rdF8sn1kw61kq/5LP9kOfEKpfdmrh7r5b2oj292z/kS2//zNZ8jWcFWq78R2pJdiNB1O1rcz+5Io35gaqqrVANVS4GoRoJoAGnzmTbXtxDK+D/IWcJRiYCcvDEgol8oXZz9O5PCznyjifPFe9s3wVG7EV3km7wxgyjp2+tYlLC5WrGeVsHUJD65WZ3SglG2F+9l3om9wRHZRs1JDt5b7BwCBAzxmOUqZ0fK1zWRruSY1etul+eg9EqtcLnfOcWVlw5j4mG8bEwp75WIYhuERTndCdwD+JCIbRWTVSHTIMMYJ5tvGhON0X7ksc841iUgZgDUisss597djNxi6GVYBQMYkfj9tGOOUk/LtrCw9cMowRpPTekJ3zjUN/WwH8HsAFyjb3OOcW+ycW5yRxyKdYYxHTta3A5nm28bYc8pP6CKSC8DnnOsd+v1qAP/2fvsEuoCZj6Uq6L4YJ7HWCjpr+cy1cH5tNQsAlP+CQ/WP3MirA/IbWMUu/Auvutl3+xyyzX6U+9i8XCn0rKRPDpTxyoLKJ/SC12W1lWQbKOOVEnmvKzHFpfwtqWg/P136Yrw8oP1S7mMyk58Jah/m+OjwGfq3s97tU8i2+qybyFaTw4M2/YXUFTbtR/S81CfLqfh2PNuHrjmp16DiNS7KrBV01vKZa+H82moWAJi9mld2Va7nVUrvTOZVG8U7+Jp2LOAVH8lM3i6vgbdrWKGkbWjhcPe4UjwdAAo28vjkHObr2nWmslxIyfzQX6EUVFdW52grZOL56e3bvlKPyz/zW1x7IF7Bc9aeVVyPIFmdev8l/pieb5/OK5dyAL8XkfeO8xvn3MuncTzDGC+YbxsTklOe0J1zdQDOHcG+GMa4wHzbmKjYskXDMAyPYBO6YRiGRxjVfOjRQsGha1PDtSe9y8JDxYscn64VdNbymWvh/IAugBY8zmLSwMdpMQNCV3DagaI93EbzlSyARhSxMxlgW04rf7b2XDKTGwHQqeS6jpawaFINzjXvi3LbjVewaJXbyG2ED3H4/tRDHLpf/xkWbbUc7gCQXBIi22dnbSHbay/x9Wu8MnVlSbRu7J5PApMGUXVTavXh3k+z0Nbfz6HxWkFnLZ+5Fs4P6AJo81K+DxK/50LdfbtZEM+YyeH37CHAYJjb/eJ5r5Ht8UevJJu/TxE1AUxeyff+/oYysi2fu4Nsbz66kGzBi1igD21hcXigXKnB8Klfk212zS1kS4Y4fQkATHuGg4zXt/BqqHOKePFCW3/q2Hb4j3MDDcOe0A3DMDyCTeiGYRgewSZ0wzAMj2ATumEYhkcYVVHUFwcC3anRZSWbWDhou5oF0JJdHFGqFXTW8pkDegSoJoBmv/AW2aLXns/bdSk5yXt5OIv3cR9jQSVPcidHZiYz9PzeOW1K3ucIy1YxJbqysIHFspx2FsYGS1g8zd/PbfSew4KVn1Oko3SrHk13cAELres6Wcw9PI+3c8O7M3Lp0E+aZFMmIt9JjfjrreVoz/KNHD1av1ITpXn8tXzmgB4BqgmglZ9kIbHrVi6+XvYr3jdZzAJo+xLu4ws/uJxsBY79MJqrX6zwfVVkm7ODhfNdU7mAc5bis8XfZr/JmcZzQbCRnba2/8tkq/kLR7jGjnMu67Zy8e+sEI9FZ5jvv8xhixekLb2p2p7QDcMwPIJN6IZhGB7BJnTDMAyPYBO6YRiGRxhVUdTf71C2KbVItMtkQYjELgCZTRwll7mpmmxaQWdAT4GrRYBqAmjg5Q1ka7yT0/RO+10r2dov4yLRWprPI7NZQKv+sV4UeHKEBaFIGQulBZu5P4iw+Fq6JUi2RDYLPU2fYkGobDM3UfEGC9ixfD2aLriR266bxdtWHVDSsvalCkzNYWVgR4lYng/NF6UWQ85p5/4MVPL5+pRMsmFeF6AWdAb0FLhaBKgmgJY8sI5sjV9n385QIowL67jjhz7H16nsJSVidoouJGp+l9fI90bLhTxJzPwdC87NV/A4lOzifnfNY9G3YiHfPwPv8P3cfZb+XJydZjXapsuVAtyx1HGIbklP8bcndMMwDI9gE7phGIZHsAndMAzDI5xwQheRB0SkXUS2HWMrEZE1IrJ36Kce8WAY4xjzbcNriHPvLySJyCUAwgAecc6dPWT7CYAu59yPROQOAMXOuX89UWPZs6pczY+/lGLLe4VrCMbyWQBwykfPRz/HdUKff1avKZrI4vPUUuDmN0TJptUbrP53bvvAv7PolMfZQBEt4POb+jLXI234qD6X9M3mPgby2ebqOFVnrlIHMryMo+QCWSwcTfkli1Od5/LY9J7LUaF5W/R0qZ/4/N/JtiyPL8xXXrmFbLmVqVGvdbffi8i+5rTjRUfSt/OLqt3CZf+cYjt8NkcpTtrK49o5n7crqOeIQl9cv1c7Fiii2hxOgTvr/3AEaONKVl+n/Af7dsZZvICg6erJZCuq41UJycz07mcAGMzn/yg4xL4dKeUx0yI2J69X6tvWcjR5bj1HUO++rYBsFeyuyG3khQYA0LCc778svs1RpIjL2c2p9+T6Xfegp+/Evn3CJ3Tn3N8ADF9isgLAw0O/Pwxg5YmOYxjjDfNtw2uc6jv0cudcy9DvrThaVNcwvID5tjFhOW1R1B19Z3Pc9zYiskpE3haRtxM9fafbnGGMGifj27Go+bYx9pzqhN4mIhUAMPSz/XgbOufucc4tds4tzijgd0qGMc44Jd/ODJhvG2PPqUaKPg/gZgA/Gvr5XDo7BVodpv84VeyRBKfGrF/JosXM+w5yJwpYAI3n6A9Usx9lNUKrAaqlwNUiQPcpAuiMOznqrvcGroUZbGfB6+DHWQCd/juu1wkAEmPhKVHKAk4iyGKNP6SksX2dx6zrXL4GDav4KbT4RY7qrLxLeVr1cfQoALzUv4xsT1VcTLbCDt638KXUSfRQ54iswj0l305kCsJVqb4TLeBx9cX42idZ38NAMetfiRz9/JKZ3I5WA1RLgatFgGoCaGInR1rHV7Io2jeFW9YiYZPHmXn6y/m8C+uVSNgKHotAj3LvN/M9lF3M0boDFbw4wz+Z75+BEt430KNcQADxoHL943x+oWm8f3b78PMboUhREXkcwDoAc0SkUURuw1FnXy4iewFcNfS3YUwozLcNr3HCJ3Tn3I3H+S8u5W0YEwjzbcNrWKSoYRiGR7AJ3TAMwyOMavrcWEEGGq5OTWc55S0W6QYqWUUJL+Bag7WX15FtV7O+bLh5OQugkTKlBqFSA1RLgatFgGoCaP4T68mWUco1IP2RmdzuxVyvEwASWUrknXIlo6yTouAAR2xGyvhzXYv2DG5m4ejIGTyG4ak81jUvKiFyAHouYbFUk39C01hsC50xTIR8dwyLispRYfRYSt/lVLI9M1hEzu7kMXQ+PhfhwwEA8pTo38FwejVAtRS4WgSoJoBW/5AjSvf8mtNPz/5/3EZolp4LONiqpByexKKhJiSXbOPo2LrVc8mm3bvBdmVwRZkfdvN90bqUI6gBIO8Q2wrreSw6lLq6ez+bKvgP/iy9Z297QjcMw/AINqEbhmF4BJvQDcMwPIJN6IZhGB5hVEVRJAH/MA3MZSjiTxYLFKIUGo3EWUwIZB2nqChrnUgGWPSIBZXPOCUATUuBq0WAagJoopNTeg4WcXTe8fAPKCKaogdG1TTEbMuIKMdL8HZa5KOmYGqRgRLTFb1Ev+KCyfR8whcdtt3YlRSFLwbktaT2sWUZn0dWJ9sic1hoy93O4nXxXt23G1bw2HzxvNfI9sIPLiebVgO0+rfcjhYBqgmgZ3yJ6+8e+q5So1QJWAZ0gbi7ltsu2cX93nMrC61lrNui62yl4SS3EXiHFwG0XcC7akIuAHSex/NBuIb9feazLOZGylOF1o7e9JzbntANwzA8gk3ohmEYHsEmdMMwDI9gE7phGIZHGFVRNBl06F2QqoYE21jY9LewcNRTw8dbUbqfbA/t4XSsABBQokJzWvnzLKeTU2Yemc2RYDUvceSjlgJXiwDVBNCC33BEafM39fqofVUstrgcFon83Xx5Y0ra7u5FLIIVlbJQU/oIi06Nl7ItXKvUzTyfxWEAuHr+ZrJNyx5eFQ64752LyJY1I7UOpE8RTkeLpB+IlKT6U+4h9uPMfvbDctYR4e/nFMThqXp0ZWYLC3qPP8r5xQoc+03ZSyy+aul4NaFbiwDVBNCa/6vUKJ03hw8IoGMJRxmX7FYibpXI4TMe4LqgHYtY2Jyyjschfw/XW637Lkf1liv3gBynLnNAuf+0KNXO+XxTDo8GTwZGKH2uYRiGMTGwCd0wDMMj2IRuGIbhEdKpWPSAiLSLyLZjbN8TkSYR2Tz072MfbDcNY+Qx3za8Rjqi6EMA/gvAI8PsdzvnfnoyjQW6gZonUz9DggdZXIzlsjBSvqaRbE8WXUa2QtbyAACVT3BNxJ5LWLBMKpGr1T9+k2wN31xCNq0G6PFS4A5HE0Arf6KEuQFIXryQbJEyFnAK32khm+th4aiggcchkVVItv23sgg26z4l/W1cqZuZrQuWbzzF5/KnmSzSVq3h6+IfSBXGGk++puhDGCHfdj4glpfax/4KFssCIT6Pkm18vv6d9WTzVZ6pth1X0k37+1jsjOZy2/1T2Fa8l6+VVgNUS4GrRYBqAmhi+27eEEDsSr4PEll8XeOKPtw7g8VFLRI9mq8cr5AXPgSzeYFEuJJF1tLN+qQzWMHHlAQvAsltYj/xDxfPlUh3jRPeAc65vwHgZQeGMcEx3za8xum8Q/+aiGwZ+trK6/UMY+Jivm1MSE51Qv8VgFkAFgBoAfCz420oIqtE5G0ReTsW5bW1hjHOOCXfTkTMt42x55QmdOdcm3Mu4ZxLArgXgJKD7L+3vcc5t9g5tzgzoES1GMY44lR9OyPHfNsYe05pQheRimP+/CSAbcfb1jAmEubbxkTmhKtcRORxAJcBKBWRRgDfBXCZiCzA0QzU9QC+lE5jvr4oghvqU2x9S2bQdl3nsOpbspNDx33KwonQPD1ndFltJdk6z1VC/9t4dcDkyDlk65sdJZvEuG2toLOWz1wL59dWswCA7+/vkK2wis+v71wurA3hYr/B1/eQreu6s8hWW91OtraFnJOh6nlekdQ/R1/tU6SsqPD3sVuGeCEOyjekru6QNFcC/Pf2I+jbgZ44ql5OHZ/+Wfz6PTSDVznsu4GXbOTXzyPblHX6aoqCjbzyYvJKjjEP38f+kMhWCkwrq0D6y5X8/0oecC2fuRbOr61mAYDyX/DKrs5VHyHbtOc6yXZwRSnZnPLIGp7KtkAvr0jxPx0kW9HefrI1Xc4FuQGg/FV2yJxOXpGUGeK5pOGa1G98SV7EpnLCCd05d6Nivj+9wxvG+MV82/AaFilqGIbhEWxCNwzD8Ag2oRuGYXiE0S0SLQLJSn27P1jInylSwvHD0SIWk+KsWSBQzOG6ADBQxqJHtEQp6hxRilGXsVAayOeQ90RpAdmcMsJaQWctn7kWzg/oAmi8qZn3v3Ka0hCblGHEYBFfl0nZvNb6EOtd6jgksvVnB19MEda6+bocruZBi+WnDq4mgI0WiRw/QvNThfukn/tc+Ww92QZum042zUdal7D4CQA5h3m89jewCD1nR4hseY18XyQDPJCF9XydBibxPakWdFbymWvh/IAugJbes45sHbfyduUbeN5IZnI7/kh6efMbP8rnUljPx8s+rOdD19KIaIskjiziJa/B5tRj+lg3VbEndMMwDI9gE7phGIZHsAndMAzDI9iEbhiG4RFGVRRN5Gaid3FqtFp2FwsUk19hETIzzIJHbiN331ev59TIe30f2arBUaqxHBaYCja3kq3tAhYmE0EWZKOsDyKar0SPKgVltXzmgB4BqgmgRY+wmOQ7lyNAQ8s5z3bZRo5K3HUNC23T/sBCm8tkMWmwgG0AEJnMY9F7NitAtfdyhN3wXNeS0MWp0SCjdwCFa1N9rOfS2bRdx3K+TgMz2LdlH98DFeuVZOMAus7kbZfP3UG2XVM54rnlQr4upe8qkcwV/OyXZE0UJbvSK+is5TMH9AhQTQAteYB9u/77vF1Oh5IDvoLvtcnv8DlX/5FtCUVkHSzRCzjnNinF3BWhdMrrXKC68aqiYfupTRD2hG4YhuERbEI3DMPwCDahG4ZheASb0A3DMDzC6EaK+gTxYRGDxe800Wb9l1WTLVDXQbaeazgPZh5nDT1KaRGZfFEWPQobuIgyIix25jYowmaIRauCAyxYOZ8S9ahouVpBZwBqClwtAlQTQJPv7uRdz1xKtiNncIeOHGIlqzRLEYILOcK1oF6P4D0yhyMVfZksrIWncdvJYUJRYvvYPZ9ES7PRcHNqMeQM5ZRzOlgoq3iF1cVoLl/Q1iXsSwDUa//mo5x6OauEN5z5ux6yHV7ISn6gR0lpvY2F8z238nU64wH2Y62gM6CnwNUiQDUBdPq3WSht+1+cpnfay3xhYkGeCmU1p4uO/noK2QIhXYzvmcn+ODwCFADar+eC7BnD1gWkGwVtT+iGYRgewSZ0wzAMj2ATumEYhkc44YQuIlNFZK2I7BCR7SLy9SF7iYisEZG9Qz+53pZhjGPMtw2vkY4oGgdwu3Nuk4jkA9goImsA3ALgz865H4nIHQDuAPCv73egjCMRFL64PcXWfR3XTuxYxrU5c9sqyBYIKakol+riW9F+Fh4ar+Dwq5x23q50CyeYDS/j2oJ4nQWPSBl/ZmZEeLvuRXzOBQ1KIU3oNUC1FLhaBKgmgOb/dj335xYWnZYu5Ha3NLLwWvN7FpPCZyl5dgFUvsbn3bePhdIjs/laV7yRKpZpqXhPwIj5tsSBrK7U9ntm8XZaDdXmi/k2DLby+WoiKwD0V/B5By86TLbib7P42nwF+3v1HzlaE81tZKpbPZdsZVwSFB2LOO2vFjEJ6OKflgJXiwDVBNDy/+QODX70fG4jwMdramOfLSjXUupyFDMA9Mzi69pXxe3kNfL1K1vXlfJ3y5H0Uv6e8AndOdfinNs09HsvgJ0AqgCsAPDw0GYPA1iZVouGMU4w3za8xkm9QxeR6QAWAngTQLlz7r1kI60Ayke0Z4YxiphvG14g7QldRPIAPA3gG865lMWrzjkHdTUsICKrRORtEXk76o7zndEwxpCR8O34AFdzMozRJq0JXUQycdThH3POPTNkbhORiqH/rwDAL04BOOfucc4tds4tDgi/GzWMsWSkfNufrQfKGMZockJRVEQEwP0Adjrn7jrmv54HcDOAHw39fO5Ex0oU5iB0TaoIGq7kz5T8XWzrqeGHpKq1HIHWGs1X2/bFlGjPRi0VppIyM5uFjEAWCyFd53I0au+5HOXmEny8olKOuktksWAFAF3XsRCp1QDVUuBqEaCaAFr8EEfdnb2a66gerJtDtmglRxp2LNBdLV+pVdlxIQultQ/zePu7Up+KfbH0hKP3GEnfdn4gUjYsnW+Sz00T3woOKMfzKTU8OYjyaDvKaYe2TCJbzjTesGQXj2u4lv04u5hldy0qu+tstk1Zx9Gx0Xz9WTLMwd9qDVAtBa4WAaoJoFl/3EC2gX+4gGxnVXPa7ObkdLLFgvq5uAy+hr44X/8wB8bDvyB1YVW8Mb38uemscrkIwOcAbBWRzUO2O3HU2Z8UkdsAHATwT2m1aBjjB/Ntw1OccEJ3zr0GQF9jBFw5st0xjNHDfNvwGhYpahiG4RFsQjcMw/AIo5o+N7s8gnmrt6bYXl07n7aLT+Gaki7CosCr37mXbB/ZcKvadvulLFCED3HkXP5+bqfpUywcTf8lr9hpWMVL14KbOUouWsB9KX2E047uv1WPQKut5kUXk7K5ba0GqJYCV4sA1QTQV+fzvi0P87VaPOsg2erfrCUbADz1bz8lW42fx+wfz1xOtkVFqarc1s9wn0cLf79D2abU69VTw7eXljrZx0OIgkMsphfu198OaWl1B8pZSAw2cnRz1zxeRDBpM9e4HKhP4jxnAAAMcElEQVTgaxJsV9TY4TmNAeTv4ePFC/UVb4He9FbCaTVAtRS4mgitCaDZf3iLbDuuWUK2WTuUpdes+QIActr4fslt5o2zu3kc4zmpz9pynDaGY0/ohmEYHsEmdMMwDI9gE7phGIZHsAndMAzDI8jRVBWjQ86UqW7mzatTbEX7WRDIPcARjn0zWJQ5PI8FmOq/KGltoafgzDrURbbec1hIzDrC4mT7IhY8sjt5LI+cwX1xihSdEWHxpnqtLvJ1LOS2B5XstNP+ECJbMosbb76Eo0cL61iFabmO1bvamzeRresLHHkaz9UFvSNzOSp03lwOQYzfzifYOaz25a7f343+DqXY6yiQWzrVzf34v6TY+iq4K/1T2d9vvvjvZHvu4Dlki77B0Z8AEC1iv9vz+V+RrfaRL5OtYiFHQzbu5nvAP1kRA4XbDbzD92n8PI7oDmYrSjAA/9N8nTsuYh+p/iPfz/6v8LkcUlLgahGgO7bVkK32a2+SreHbnKI3YxGLvgCANzjiNsHldpF9WFmwcWnqPNZ4568wsL/phL5tT+iGYRgewSZ0wzAMj2ATumEYhkewCd0wDMMjjKooWphV7i6svCnFFjq/krbrWMifM9V/ZhGlt4YVhm7OLAsAqH2YayweXMH5SP2KplrxBos6e77Kbc+5i6M1D32cRRmfEgAarmVj7YO6cJTZ0k22RCmnrHV+JT1wMUcVBuv4eFoK3JLvcwTogftZ9S15kFPv9t7AtUwBINDDImHnORzBm9POfjr576ni1rpDjyA0oBTjHAUK8qrckvn/M8XmG2AxL1LBAnTrUj7fbHZX5DXr6YG19Lktn+ZI05qHWBAfKGablvZ3oIR9qXg3t9F2AUd6luzgcQhX6ulgJ23jG9ApCxoSii1apNTwVGqAauc3SYkAbbqEFx9M/T7XKG39OgulAFCyi+/fnDpeiBFaqER0z049v/r770Kk5cSCvz2hG4ZheASb0A3DMDyCTeiGYRge4YQTuohMFZG1IrJDRLaLyNeH7N8TkSYR2Tz072MffHcNY+Qw3za8Rjrpc+MAbnfObRKRfAAbRWTN0P/d7Zzj/KfHIZEbQPfSqhRb0SZOBdt2fjnZsho5GuvALUrNzV4WmAAgfAZHbfkUMal0Kws9sXwWQPO2sLgIH0d21rzIgqMotS87z+cowGS2LoL1z2ERJZGtCKAFLAgV1LP4Ez5Lic5TaoBqKXDzlQhQTQDNf2I92QDgwA+VqNIgR6nGCridaFFF6t+P6df+fRgx344HfWg/LzVKMjSHr5+Wnrl4N59v4R6Olg7VchQmALSvZJ9NhthnY8q16j6L/WbqKyzuB3p4bFuXsgAabGXBUZSFF6Wb+fwAoOlyTuerRVIOlih1fkO8XWE9LzZQa4Aq6Wm1CFBNAJ3yHyyUAsCeB88jm7+N57aSHbxv9V9Sx6epN738uemUoGsB0DL0e6+I7ARQ9f57Gcb4x3zb8Bon9Q5dRKYDWAjgvSQHXxORLSLygIgUH3dHwxjnmG8bXiDtCV1E8gA8DeAbzrkeAL8CMAvAAhx9yvnZcfZbJSJvi8jbsUH+KmcYY81I+HY8Yr5tjD1pTegikomjDv+Yc+4ZAHDOtTnnEs65JIB7AXBdp6Pb3eOcW+ycW5yZxUEVhjGWjJRv+3PMt42xJ51VLgLgfgA7nXN3HWM/VpH6JIBtI989w/jgMN82vEY6q1wuAvA5AFtFZPOQ7U4AN4rIAgAOQD2AL53oQL64Q7AtNRx27xdZ9Z10dgfZ6j/N201ew6q2/6Y2te3e7VPIllzC+cIPLmA1P7gxSLZPfI5zWL/Uv4xsPZfwypdEPw/71fM3k+2NpxaSDQCK9vLqCV+MxyIymVcCHJnDKxMqX+PQ7Px6Pp5W0PnSF1aTrfAF3ldbzQIAM77FaQL6rufivG1L+FwmDVuRVN9/0mksRs63E0DO4dSVCBW/4BVc3RewH0dv5HDwPd28oiVwQH/+OvNbfL9Me6aTbOu2LiJbNm+GhuX8bSMe5LHNO8T7dp7HqzEC3ezvgxV6MejyV3n/ZAZf+9wm3q5nJo9PzyylUHcGn4tW0BlvsE0L59dWswDAGV/YSLaMeXPIdnAFr3DruCJ1Hhr8bnoZLdJZ5fIaAO1oL6XVgmGMU8y3Da9hkaKGYRgewSZ0wzAMj2ATumEYhkdIRxQdMRJZgtDM1JD5RBmLDF+Z+VeyfW//9WSLzGBx8O5Zf1LbXn3WTWT77KwtZFvXOYNsdbM4jHpZ3h6yPVVxMdlUKSPJ1mnZLIz9aSaLlQDg7+PLlt3NIlHv2Ty2vkwes759LFB1XMht1/hZqNMKOtcf5DHUwvkBXQDNfZqL8876KgdwHgxPS/k7tm1MUqEDACThOLd7J19T/wCnbWht5LQUmUc4RUAOa6wAgHgFxz2tb2FhMyuUXvi4E37O88V5bLWw+nAN+2YeuwgkoadpyOnkYyayuG2nCKXBZhY7+6p4O+1ccpsVkXUaj4OWz1wL5wd0ATSxfTfZApdxOoFIaNj4JNLzbXtCNwzD8Ag2oRuGYXgEm9ANwzA8gk3ohmEYHmFUi0TnVEx1M25JjSzM4PTcajTWEUWYzFQiA/MP6YWVEznp5UA+PI/FmsIDLCQ2X8H7Fu5kISu0mHNVQ9Ib88rndeEoNJPbGSzhY05/ngvuhqdx9NvwgrQAUL2W9837YRPZ+v+FRb7D87nAdM9sMgEAksopzlrMIYjuCm57eN71ba/8HOHDJy6k+0FQGKx0S+f8jxRby6UsVjqld3017IhFO3jD8rUtatt7VlWQ7ZwL95Gt82csVjddruRDX8P+HprGFyrKWi7l8QaAzvks0Mpx9NmS7ex3HYt4/ymvc67yuuu5PkJeI7cRrmZbxXpeBNAzlQXenC6+z+LZusuFq7Wc7bxd2S85n3rXF1Ijq3c9dzf6Oq1ItGEYxocGm9ANwzA8gk3ohmEYHsEmdMMwDI8wqpGimb1JVK3tTbF1nMfRh82XsABTvINVlGAbR5U1XMXiKQBMf4HFlsYrWWxxrDcis4/bzq3kCjWFL/HxQmfwEPuirG1kzeglm39A6QyA8g1K4dt8LU0ot5NUDlnxBgu3/i4+v0VFHPL37EJWO8v+2kq24QWd32N4ClyAI0ABoPiGqWQbXnja58aualAsLwOty1JF0NBcFtqK3+UL4A/zc5UkWXxruUYfw2Q1ryxo6+diy5lRpYBzjH0ku5nvlex27uPez7K/R8o56liL9PQfJ9VxwzV8TC0CtPEqVmQzlPUQZeuUyM4FLFbHlUUT4Ut5HAa3cv80IRjgFLiAEgEKwP8FTi1d8mBqWumMNH3bntANwzA8gk3ohmEYHsEmdMMwDI9gE7phGIZHGNVIURHpAHAQQCkApZrhhMTOZfwwzTk3eSwaNt8e90z0c0nLt0d1Qv/vRkXeds4tHvWGPwDsXIxj8dIY2rlMPOyVi2EYhkewCd0wDMMjjNWEfs8YtftBYOdiHIuXxtDOZYIxJu/QDcMwjJHHXrkYhmF4hFGf0EXkWhHZLSL7ROSO0W7/dBCRB0SkXUS2HWMrEZE1IrJ36CcnihhniMhUEVkrIjtEZLuIfH3IPuHOZTxhvj32fNh9e1QndBHJAPBLAB8FMBfAjSIydzT7cJo8BODaYbY7APzZOVcL4M9Df4934gBud87NBbAUwFeHrsNEPJdxgfn2uOFD7duj/YR+AYB9zrk651wUwBMAVoxyH04Z59zfAAxP37YCwMNDvz8MYOWoduoUcM61OOc2Df3eC2AngCpMwHMZR5hvjwM+7L492hN6FYBjc7A2DtkmMuXOufeKPbYCKB/LzpwsIjIdwEIAb2KCn8sYY749zvgw+raJoiOIO7pkaMIsGxKRPABPA/iGc67n2P+baOdifLBMNH/4sPr2aE/oTQCOrVRQPWSbyLSJSAUADP1sH+P+pIWIZOKowz/mnHtmyDwhz2WcYL49Tvgw+/ZoT+gbANSKyAwRCQC4AcDzo9yHkeZ5ADcP/X4zgOfGsC9pISIC4H4AO51zdx3zXxPuXMYR5tvjgA+7b496YJGIfAzAzwFkAHjAOfeDUe3AaSAijwO4DEczt7UB+C6AZwE8CaAGR7Pt/ZNzjutejSNEZBmAvwPYCuC9+np34ui7xgl1LuMJ8+2x58Pu2xYpahiG4RFMFDUMw/AINqEbhmF4BJvQDcMwPIJN6IZhGB7BJnTDMAyPYBO6YRiGR7AJ3TAMwyPYhG4YhuER/j/HhnTjJsVj2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "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.135088e-02|0.000000e+00\n",
      "    1|1.414971e-02|-1.215656e+00\n",
      "    2|9.934682e-03|-4.242736e-01\n",
      "    3|7.486162e-03|-3.270729e-01\n",
      "    4|7.415422e-03|-9.539599e-03\n",
      "    5|6.433930e-03|-1.525492e-01\n",
      "    6|6.020392e-03|-6.868966e-02\n",
      "    7|6.012738e-03|-1.272962e-03\n",
      "    8|6.012661e-03|-1.278831e-05\n",
      "    9|6.012660e-03|-1.278889e-07\n",
      "   10|6.012660e-03|-1.278889e-09\n",
      "   11|6.012660e-03|-1.278958e-11\n",
      "It.  |Err         \n",
      "-------------------\n",
      "    0|7.283759e-02|\n",
      "   10|4.751585e-03|\n",
      "   20|4.981526e-09|\n",
      "   30|3.401818e-14|\n",
      "Gromov-Wasserstein distances: 0.0060126599835825445\n",
      "Entropic Gromov-Wasserstein distances: 0.00591822665714488\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEtCAYAAAAsgeXEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHyJJREFUeJzt3Xu0ZHV14PHvtmmV5t0+2gbR9oEmhIyoHXR8okFDEhJ0mUGJIkSlMaOORhwlJBlvZySDjo/RQVlplQAqKhERdNRojIhPpGH5AjQQaERsuuWhgKLy2PPH+V05fbn3nur7q+e9389atW7V+Z06Z9e5dXbtOnVqV2QmkiRJWph7jToASZKkSWYxJUmSVMFiSpIkqYLFlCRJUgWLKUmSpAoWU5IkSRUspqQFiIgXRsTnRh2HNO4i4iERcWtELBt1LBqNiLgkIg4cdRyDZDE1AhHxgoi4ICJ+HhFby/X/GhEx6th6FRE7lAT5hNa0F0ZEzjLt+6OJsncRMRURH+x1/sz8UGY+e5AxSb2KiE0RcVvZJ6cvJ/V43/Mi4mWDii0zf5iZO2fmndt734hYHRHvjYgfl8d0ZUScGhG/NYhYByUi/joiPjNj2uVzTHvBcKPbfiXPP7LX+TPzdzLzvAGGNHIWU0MWEccC7wT+N/AgYBXwcuDJwL3nuM/YvaPLzDuArwNPa01+GvD9WaadP8TQZhURO4w6BmnA/qQULdOXV/ZjoaPadyLifsDXgBXAU4FdgMcBXwKeNcd9xnU/Px940nQuj4jVwHLgsTOmPRLz5WTKTC9DugC7AT8Hntcx36nAycCny/wHlfueDvwEuBr4W+BeZf6jgK8C7wB+ClwJPKlMvwbYChw5I457LAu4T7n/fq15HwDcBjxwljj/Dvhk6/alZZ0zp72oXD+ApgD7KbAZOAm4dxmLEv9W4Gbgu9NxAH9UlnMLcC3wutbyDwG+VZb5NeA/tcY2AW8AvgP8Ctih3L62LOsHwO8DBwO/Bm4HbgW+3dpO7y+xXgu8CVjW2uZfaa0raYriy0ss7wZi1M85L0vjUp7rB80xdhTwFeCtwE3AVcAflrETgDuBX5bn/kllegKvKM/nq8q0JwEXAj8rf5/UWsd5wP8Cvln233OAlWVsTVneDuX2SuCfgB+XeD4xR9xvAr5NyXNzzDO97JcCPwTOL9P/FLik7IvnAb89Y1v995IXfl728VXAZ0pe+Fdgj9b8sy6r5JKPzYjnncC7Zonz3sAvgMeX24eVbfClGdOumLGsa8r2vAh4amvsAGBjGdsCvL1Mvy/wQeCGEu+FwKoy1pXPpl9Dbihjjyzx/Qy4Hvhomff8ss1/TvOceX6Z3pWLDyrXp4AzaV6Dbinbdu2o96HqfXDUASylC82L9h3TSWWe+U4tT+An0xQ59y1PvHNo3p2tAf4deGmZ/6iy3L8AlpUd4Yc0L+j3AZ5dnrQ7l/nnW9YpwAmtWF4BfHaOOJ8O3FhivD9NYbai7NzT0xJ4SJn/8cATaYqaNcBlwGvK2B+UhLE7TWH128DqMraZkkiAPYDHleuPpSm+nlAe95Flp71PGd9Udu69gR2BR9Mkpz3L+BrgEeX6FPDBGY/vbOAfgZ2AB9K8UBzT2uYzi6lPlfgfQlOoHjzq55yXpXGhu5i6HTi67Cd/SVPIRBk/D3jZjPsk8HmawmfH8vcm4Iiy/x5ebt+vtYxrgf3K/nLW9P7EPYup/wd8tOzLy4GnzxH3N4Cpjsc9vezTy3p3BB5F80L/rLL81wNXcPcbt01l2auAvUoOubjkk/sC/wa8scw757KAh9IUSLuUeZfR5KonzhHrF4G/KtdPAl5CU8y2p53Smv9FwP3K9j4WuA64bxn7OnBEub7z9DqBY4BP0uThZTQ5d9cy1pXP7gBeVda3I/Bh4G+4+zXoKTOeH49s3e4lF7eLqV/SvEleRlOEf2PU+1D1PjjqAJbSpewc182Y9jWaSv424Gll2qnA6a15ltEcOdm3Ne0Y4Lxy/Sjg8tbY75Yn+6rWtBuA/XtY1kHAf7TGvgq8eI7Hc9+yUzwGeC7woTL9G61pV82zPV4DnF2uP5OmqHsiM96J0hSGx0wnhdb0k4H/OWPaDyjJuezAL2mNPbLs8AcBy2fcb4pWMUWTaH8F7NiadjjwxdY2n1lMtZPNmcBxo37OeVkal/Jcv7XkkunL0WXsKLY94rGiPF8fVG6fx+zF1DNbt48Avjljnq8DR7WWcWJrbN+SZ5bRKqaA1cBdtI78zPOYrgBe3rr9p+Vx3QJ8rkybXvbDW/P9HXBm6/a9aAq9A1vb6oWt8bOAk1u3X0U5WtbDsr5CyY80Bdd/zPN4prg7330b2IfmDXZ72pHz3P8m4DHl+vnAeuD+M+Z5CTOOCpXpveSzH864z+nABuDBs8Qys5jqJRe3i6l/nfFcuW3U+1DtxXOmhusG4P7tz6Mz80mZuXsZa/8/rmldvz/Nu6KrW9OupnlXNW1L6/ptZdkzp+3cw7K+CKyIiCdExBqaAuzs2R5MZv6S5t3N08rly2XoK61pv/n8PyIeFRGfiojrIuJm4B9KPGTmv9G8M3s3sDUiNkTEruWuz6N5F3N1RHwpIv5zmf5Q4NiI+On0heYo1J6tMH+zHTPzCpoCbqqs4yMR0Z637aFlO21uLfsfad7RzeW61vVf0GxvaViek5m7ty7vbY395rmZmb8oV7uen+0ctCfb5gy4Zw66ZsbYcsr+3bI3cGNm3tSxbmhy4urpG5l5bsmVf8U9zy+dM9bMvKuMz5cvZ8uVvSzrDJqiBODPy+25nA88JSJWAg/IzMtpCp8nlWn7sW2+fF1EXBYRPyv5Zzfu3p4vpTlq9v2IuDAiDinTPwD8C/CRctL+WyJiOb3ls/Y2hOYoXADfLN/Ge8k8j62XXNw2M1fed9LP07KYGq6v07w7OLSHebN1/Xqaw/QPbU17CM07pO0177Ky+cbNmTQJ4nDgU5l5yzzLO5+maHoqdxdTX25Na59MeTLNCer7ZOauwPE0Oytl3e/KzMfTvFN5FM15DWTmhZl5KM2O/4kSHzQ7/wkzXkBWZOaHW+tsb0cy84zMfEp5/Am8ebb5yrJ/RfPOb3rZu2bm78yzLaRJNPO5P9v0H7NtzoB75qC9Z4zdTpNv2q4BVkbE7j3E9QXgORHRy+vUnLGWb0nvzcLyZdey/hk4MCIeTHMkfr5i6us0BdHRNEf8ycybyzqOBn6cmVeV9TyVppg5jOYo3u40p35Eud/lmXk4TU58M/CxiNgpM2/PzPWZuS/NOW6HAC+mt3w2M1del5lHZ+aeNJ8MvGeeb/D1kosXNYupIcrMn9Icmn1PRPxZROwSEfeKiP1pPsee637TBc4J5T4PBV5Lc6Lh9sbQy7LOAJ4PvJD5kwM0xdIzaBLMpWXaV4EDaY5qtYupXWhOmLy1fLX5L6cHIuL3ytGw5TTnKPwSuCsi7l3aK+yWmbeX+99V7vZe4OXlfhERO0XEH0fELrMFGhGPjohnRsR9yvJvay1rC7BmOnFn5mbgc8DbImLX8n96REQ8vWN7SJNmC/Dwjnk+DTwqIv68tEV5Ps2bnk+15nlRROwbESuAv6c5OXubdghlv/oMTQ7cIyKWR0T7279tb6c5r+oDZd+Lsm/v3xHrmcAfR8Tvl3xyLE0h8bWO+233sjLzJzQfcf4TzSkNl821oMy8jeak8ddy9xtPaI7kv5Z75so7aM693CEi/gcwfaSeiHhRRDygHCn7aZl8V0Q8IyJ+t3xD8GaagvauheSziPgvpUiE5iPGZNt82X7ObFcuXowspoYsM99Cs+O8nuYJuYXmcOsbmH9nfxVNkXElzc53Bs3J4gsx77Iy84IyvidN4pvP12jebV2Q5QPwzLyeJglsLYeyp72O5lD4LTQ730dbY7uWaTfRHFa/gaZ9BDTna2wqHw2+nKbIIzM30ryjO6nc7wqaz/7nch/gRJp3y9fRvKv76zL2z+XvDRFxcbn+YpqPEy4ty/8YrY8dpDHzydi2z9SsH8/P4p3An0XETRHxrtlmyMwbaI5yHEuzb74eOKTs69M+QHO+53U051P+tznWdwTNi/z3ac5hfM0c67ye5hzKX9LkqVtovlCyC603YrPc7wc056f+X5p9/U9o2kb8eq77VC7rDJrzMLveeELz7bgHlscz7ctlWruY+hfgszTnkV5Nsw3aH8MdDFwSEbfS/P9eUIq1B9HkqZtpvuDzJZr/C2x/Pvs94IKyjnOBV2fmlWVsCjitfKR32AJy8aIz/W0OSZIWJCLOo/kCx/tGHYs0Ch6ZkiRJqmAxJUmSVMGP+SRJkip4ZEqSJKmCxZQkSVKFqo6jEXEwzdcylwHvy8wT559/RTY/XSZpEqxmc+c8mzu7RWy+PjMf0J+I+mt7cljEyoQHzzXMPftTasnYba5G3y0/+/G8w4/fp3tfu+hyO7MMX2/5a8HFVGkK9m6a3yP6EXBhRJybmZfOfa/dgXULXaWkIVvH+s551nfu0+tn/gzJWNj+HPZgmt6Vc7ErwJL11KnueT41/zwbT+re1+IPfP0cvt7yV83HfAfQ/HjmlaWB2Ufo7WdSJGkcmMMk9UVNMbUX23Zk/RHb/pCkJI0zc5ikvhj4rzRHxDp+89neboNenST1zbb5yzpL0uxqjkxdy7a/Ev5gZvlV7szckJlrM3MtrKhYnST1VWcO2zZ/rRxqcJImR00xdSGwT0Q8LCLuDbyA5scQJWkSmMMk9cWCP+bLzDsi4pU0v269DDglMy/pW2SSNEDmMEn9MtSfk4nYM2tbI7yx46va63lj1fIl9dv6i5qPySZbP/IXz5maf/wTHeMakDUd45uq15Bvnv+1K97ga9d46i1/2QFdkiSpgsWUJElSBYspSZKkChZTkiRJFSymJEmSKlhMSZIkVbCYkiRJqmAxJUmSVGHimnZKo2bj2O1l087+mqoc1/bKfToabl7uPr942bRTkiRp4CymJEmSKlhMSZIkVbCYkiRJqmAxJUmSVMFiSpIkqYLFlCRJUoUdRh2ANGmG0Ueqq5fVsOJQf+XxHf2K/qGX/+lUX2JR7/rSR+pFU/OPf7BjXGPNI1OSJEkVLKYkSZIqWExJkiRVsJiSJEmqYDElSZJUwWJKkiSpgsWUJElSBYspSZKkCpGZw1tZ7JmwbmjrkzR4XQ1G18NFmbl2SOEMjPlrMtkAV3XW95S/PDIlSZJUwWJKkiSpgsWUJElSBYspSZKkChZTkiRJFSymJEmSKlhMSZIkVdhh1AH0mz1FpOHq3p+690n1Lh82//aMq8xvbf3J9z7H+29lx/iNQ4miX6qKqYjYBNwC3AncsRga80laOsxhkvqhH0emnpGZ1/dhOZI0CuYwSVU8Z0qSJKlCbTGVwOci4qKI8EerJE0ac5ikarUf8z0lM6+NiAcCn4+I72fm+e0ZSoIqSWq3ytVJUl/Nm8PMX5J6UXVkKjOvLX+3AmcDB8wyz4bMXNuc2LmiZnWS1FddOcz8JakXCy6mImKniNhl+jrwbOB7/QpMkgbJHCapX2o+5lsFnB0R08s5IzM/25eoJGnwzGGS+mLBxVRmXgk8po+x9MVSa8jZ1aR0qW0PqVfjmsO6jE1TzrVT849v7BgH8mUdDUjf954eAtnSwzy1bMrZf5PVlLOLrREkSZIqWExJkiRVsJiSJEmqYDElSZJUwWJKkiSpgsWUJElSBYspSZKkCrW/zacRs4/U9rEvl9QnPfSR6hLv69rfhtFDqhcrO8Ynp2fSOzty4KvNgQvikSlJkqQKFlOSJEkVLKYkSZIqWExJkiRVsJiSJEmqYDElSZJUwWJKkiSpgsWUJElSBZt2akk1slxMj0XSsIxHU85+5Gqbcg6GR6YkSZIqWExJkiRVsJiSJEmqYDElSZJUwWJKkiSpgsWUJElSBYspSZKkCvaZkr2X+qyrFwy4zTX5zspvzjv+vDigh6X8Tcf4W3pYxu09zLM4DCVvfGOqe54n9jDPEuORKUmSpAoWU5IkSRUspiRJkipYTEmSJFWwmJIkSapgMSVJklTBYkqSJKmCxZQkSVKFzqadEXEKcAiwNTP3K9NWAh8F1gCbgMMy86bBhbk02fxxMvk/GS+TlMMmaZ/vrSlnlxP6sAz1lQ05F6SXI1OnAgfPmHYc8IXM3Af4QrktSePoVMxhkgaos5jKzPOBG2dMPhQ4rVw/DXhOn+OSpL4wh0katIWeM7UqMzeX69cBq/oUjyQNgzlMUt9Un4CemQnkXOMRsS4iNkbERvhF7eokqa/my2HmL0m9WGgxtSUiVgOUv1vnmjEzN2Tm2sxcCysWuDpJ6quecpj5S1IvFlpMnQscWa4fCZzTn3AkaSjMYZL6prOYiogPA18HHh0RP4qIlwInAs+KiMuBg8ptSRo75jBJgxbN6QLDsWdErptnfFz6p0jqp/UXNR+TTbaIPRPmy2CSYLL6pXXrLX/ZAV2SJKmCxZQkSVIFiylJkqQKFlOSJEkVLKYkSZIqWExJkiRVsJiSJEmqYDElSZJUYYdhrmwzq1lv0ztpKBZX4zxpdl3Pc5/jw7cUt7lHpiRJkipYTEmSJFWwmJIkSapgMSVJklTBYkqSJKmCxZQkSVIFiylJkqQKQ+0z1cW+OFL/uK9oLoupN9MkxarFyyNTkiRJFSymJEmSKlhMSZIkVbCYkiRJqmAxJUmSVMFiSpIkqYLFlCRJUgWLKUmSpApj1bRzkpqvLaamd5KWlmHkJ5swaynxyJQkSVIFiylJkqQKFlOSJEkVLKYkSZIqWExJkiRVsJiSJEmqYDElSZJUobPPVEScAhwCbM3M/cq0KeBo4CdltuMz89ODCnIc9aM/ir2qpMEbrxy2vGP89sGH0C9vnZp3eP3relnIQR3jX+phGZOyzSbjf3/XDd39we51P1+bZurlyNSpwMGzTH9HZu5fLkuqkJI0UU7FHCZpgDqLqcw8H7hxCLFIUt+ZwyQNWs05U6+MiO9ExCkRsUffIpKk4TCHSeqLhRZTJwOPAPYHNgNvm2vGiFgXERsjYiP8YoGrk6S+6imHmb8k9WJBxVRmbsnMOzPzLuC9wAHzzLshM9dm5lpYsdA4Jalves1h5i9JvVhQMRURq1s3nwt8rz/hSNLgmcMk9VMvrRE+DBwI3D8ifgS8ETgwIvYHEtgEHDPAGCVpwcxhkgYtMnN4K4s9E9YNbX0L1dX/CewBJfVu/UXNx2STbe3OkRv3m3s8LjAnLFk7THXPc0fXPCt7WJFfSh2+3vKXHdAlSZIqWExJkiRVsJiSJEmqYDElSZJUwWJKkiSpgsWUJElSBYspSZKkChZTkiRJFTo7oC9FNuScTDZb1SBd9PPVxAVzNx0+sYfn33E+/xanzoacvbAh57aW9zDP7QOPolcemZIkSapgMSVJklTBYkqSJKmCxZQkSVIFiylJkqQKFlOSJEkVLKYkSZIq2GdqkVtKvZcWy+PQZLKH1EKsmXf0HP6icwmHut1bVvYwz6T0sxqfHlK98MiUJElSBYspSZKkChZTkiRJFSymJEmSKlhMSZIkVbCYkiRJqmAxJUmSVMFiSpIkqYJNOxe5xdTIsqsB6WJ6rNLSsGne0WE15PxAR245oi9x7Ngxflsf1tGPhpyrephnSx/Ws7h4ZEqSJKmCxZQkSVIFiylJkqQKFlOSJEkVLKYkSZIqWExJkiRVsJiSJEmqYJ8pjYWuHlJgHymN2s7Ak+cZv7iHZfSjl5D6rbaPVF7bnb9ir451nDTVvaJX9jBPNXtILUTnkamI2DsivhgRl0bEJRHx6jJ9ZUR8PiIuL3/3GHy4ktQ785ekYejlY747gGMzc1/gicArImJf4DjgC5m5D/CFcluSxon5S9LAdRZTmbk5My8u128BLgP2Ag4FTiuznQY8Z1BBStJCmL8kDcN2nTMVEWuAxwIXAKsyc3MZuo45ftAnItYB65pbuy0sSkmqVJ+/HjjoECVNqJ6/zRcROwNnAa/JzJvbY5mZQM52v8zckJlrM3MtrKgKVpIWoj/5yzeDkmbXUzEVEctpEtGHMvPjZfKWiFhdxlcDWwcToiQtnPlL0qD18m2+AN4PXJaZb28NnQscWa4fCZzT//AkaeHMX5KGoZdzpp4MHAF8NyK+VaYdD5wInBkRLwWuBg4bTIiStGDmL0kDF83pAkNaWeyZvzmXc4S6GkTaHFLqp/UXNeccTbZxyV/dduxhHpuHLkYbemh+vM7Xtxm69pfjespf/pyMJElSBYspSZKkChZTkiRJFSymJEmSKlhMSZIkVbCYkiRJqmAxJUmSVGG7fuh4sRhGH6muXlbDikPSUlPfQ8r8NZnsIbUQ/em55pEpSZKkChZTkiRJFSymJEmSKlhMSZIkVbCYkiRJqmAxJUmSVMFiSpIkqYLFlCRJUoUl2bRzGJZaQ7uuJn9LbXtIk2x89tdde5jn5oFHsbjs2DHenyaWXTZ0vGZMWgNSj0xJkiRVsJiSJEmqYDElSZJUwWJKkiSpgsWUJElSBYspSZKkChZTkiRJFewzpb4Yn740dbr6ZcHieazaXsuYv+/RePQ7yt/qfg7H9yflOTwe23RxGU4fqS7D6SM11ad5unlkSpIkqYLFlCRJUgWLKUmSpAoWU5IkSRUspiRJkipYTEmSJFWwmJIkSapgMSVJklShs2lnROwNnA6sAhLYkJnvjIgp4GjgJ2XW4zPz04MKVBqGpdaQs6tJ6aRvj/7mrzuZhCaSk9OQU6rzpo789bdDigN664B+B3BsZl4cEbsAF0XE58vYOzLzrYMLT5KqmL8kDVxnMZWZm4HN5fotEXEZsNegA5OkWuYvScOwXedMRcQa4LHABWXSKyPiOxFxSkTs0efYJKlvzF+SBqXnYioidgbOAl6TmTcDJwOPAPaneef3tjnuty4iNkbERvhFH0KWpO1j/pI0SD0VUxGxnCYRfSgzPw6QmVsy887MvAt4L3DAbPfNzA2ZuTYz18KKfsUtST0xf0katM5iKiICeD9wWWa+vTV9dWu25wLf6394krRw5i9Jw9DLt/meDBwBfDcivlWmHQ8cHhH703zdeBNwzEAilKSFM39JGrhevs33FSBmGbKn1BLR1YsIJr8f0VK12P9v5i9pMvXyuvO3Y5S/7IAuSZJUwWJKkiSpgsWUJElSBYspSZKkChZTkiRJFSymJEmSKlhMSZIkVbCYkiRJqtBLB3QtcYu9seModDWkc5tLGlfDyF+TlgM9MiVJklTBYkqSJKmCxZQkSVIFiylJkqQKFlOSJEkVLKYkSZIqWExJkiRVsM+UhsK+Sttaao9X0jAs72Ge26vX0p/81RVrfZy9OLHjtem4HpfjkSlJkqQKFlOSJEkVLKYkSZIqWExJkiRVsJiSJEmqYDElSZJUwWJKkiSpgsWUJElSBZt2aihsUtl/NkLV7Kb6NI8mz3AaXfbH4GPNh82fIwHiqq482b0M8MiUJElSFYspSZKkChZTkiRJFSymJEmSKlhMSZIkVbCYkiRJqmAxJUmSVCEyc3gri/gJcHVr0v2B64cWwMJNSpwwObEaZ/+Na6wPzcwHjDqIWrPkLxjfbT6TcfbXpMQJkxPruMbZU/4aajF1j5VHbMzMtSMLoEeTEidMTqzG2X+TFOtiMSnb3Dj7a1LihMmJdVLinIsf80mSJFWwmJIkSaow6mJqw4jX36tJiRMmJ1bj7L9JinWxmJRtbpz9NSlxwuTEOilxzmqk50xJkiRNulEfmZIkSZpoIyumIuLgiPhBRFwREceNKo4uEbEpIr4bEd+KiI2jjqctIk6JiK0R8b3WtJUR8fmIuLz83WOUMZaYZotzKiKuLdv1WxHxR6OMscS0d0R8MSIujYhLIuLVZfpYbdN54hy7bbpYmb/qmb/6y/w1WiP5mC8ilgH/DjwL+BFwIXB4Zl469GA6RMQmYG1mjl3/i4h4GnArcHpm7lemvQW4MTNPLEl+j8x8wxjGOQXcmplvHWVsbRGxGlidmRdHxC7ARcBzgKMYo206T5yHMWbbdDEyf/WH+au/zF+jNaojUwcAV2TmlZn5a+AjwKEjimViZeb5wI0zJh8KnFaun0bzJB2pOeIcO5m5OTMvLtdvAS4D9mLMtuk8cWo4zF99YP7qL/PXaI2qmNoLuKZ1+0eM78ZM4HMRcVFErBt1MD1YlZmby/XrgFWjDKbDKyPiO+Uw+sgP57dFxBrgscAFjPE2nREnjPE2XUTMX4MztvvaLMZ2XzN/DZ8noHd7SmY+DvhD4BXlkO9EyOYz3HH9uubJwCOA/YHNwNtGG87dImJn4CzgNZl5c3tsnLbpLHGO7TbVyJi/BmNs9zXz12iMqpi6Fti7dfvBZdrYycxry9+twNk0h/jH2ZbymfT0Z9NbRxzPrDJzS2bemZl3Ae9lTLZrRCyn2cE/lJkfL5PHbpvOFue4btNFyPw1OGO3r81mXPc189fojKqYuhDYJyIeFhH3Bl4AnDuiWOYUETuVE+SIiJ2AZwPfm/9eI3cucGS5fiRwzghjmdP0zl08lzHYrhERwPuByzLz7a2hsdqmc8U5jtt0kTJ/Dc5Y7WtzGcd9zfw1WiNr2lm+9vh/gGXAKZl5wkgCmUdEPJzm3RzADsAZ4xRnRHwYOJDm17a3AG8EPgGcCTyE5hfuD8vMkZ48OUecB9Iczk1gE3BM63P9kYiIpwBfBr4L3FUmH0/zef7YbNN54jycMdumi5X5q575q7/MX6NlB3RJkqQKnoAuSZJUwWJKkiSpgsWUJElSBYspSZKkChZTkiRJFSymJEmSKlhMSZIkVbCYkiRJqvD/AY+VzKLTcnK+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}