1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
|
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Gromov-Wasserstein Barycenter example\n",
"\n",
"\n",
"This example is designed to show how to use the Gromov-Wasserstein 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",
"\n",
"import numpy as np\n",
"import scipy as sp\n",
"\n",
"import scipy.ndimage as spi\n",
"import matplotlib.pylab as pl\n",
"from sklearn import manifold\n",
"from sklearn.decomposition import PCA\n",
"\n",
"import ot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Smacof MDS\n",
"----------\n",
"\n",
"This function allows to find an embedding of points given a dissimilarity matrix\n",
"that will be given by the output of the algorithm\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def smacof_mds(C, dim, max_iter=3000, eps=1e-9):\n",
" \"\"\"\n",
" Returns an interpolated point cloud following the dissimilarity matrix C\n",
" using SMACOF multidimensional scaling (MDS) in specific dimensionned\n",
" target space\n",
"\n",
" Parameters\n",
" ----------\n",
" C : ndarray, shape (ns, ns)\n",
" dissimilarity matrix\n",
" dim : int\n",
" dimension of the targeted space\n",
" max_iter : int\n",
" Maximum number of iterations of the SMACOF algorithm for a single run\n",
" eps : float\n",
" relative tolerance w.r.t stress to declare converge\n",
"\n",
" Returns\n",
" -------\n",
" npos : ndarray, shape (R, dim)\n",
" Embedded coordinates of the interpolated point cloud (defined with\n",
" one isometry)\n",
" \"\"\"\n",
"\n",
" rng = np.random.RandomState(seed=3)\n",
"\n",
" mds = manifold.MDS(\n",
" dim,\n",
" max_iter=max_iter,\n",
" eps=1e-9,\n",
" dissimilarity='precomputed',\n",
" n_init=1)\n",
" pos = mds.fit(C).embedding_\n",
"\n",
" nmds = manifold.MDS(\n",
" 2,\n",
" max_iter=max_iter,\n",
" eps=1e-9,\n",
" dissimilarity=\"precomputed\",\n",
" random_state=rng,\n",
" n_init=1)\n",
" npos = nmds.fit_transform(C, init=pos)\n",
"\n",
" return npos"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Data preparation\n",
"----------------\n",
"\n",
"The four distributions are constructed from 4 simple images\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/rflamary/.local/lib/python3.5/site-packages/ipykernel_launcher.py:6: DeprecationWarning: `imread` is deprecated!\n",
"`imread` is deprecated in SciPy 1.0.0.\n",
"Use ``matplotlib.pyplot.imread`` instead.\n",
" \n",
"/home/rflamary/.local/lib/python3.5/site-packages/ipykernel_launcher.py:7: DeprecationWarning: `imread` is deprecated!\n",
"`imread` is deprecated in SciPy 1.0.0.\n",
"Use ``matplotlib.pyplot.imread`` instead.\n",
" import sys\n",
"/home/rflamary/.local/lib/python3.5/site-packages/ipykernel_launcher.py:8: DeprecationWarning: `imread` is deprecated!\n",
"`imread` is deprecated in SciPy 1.0.0.\n",
"Use ``matplotlib.pyplot.imread`` instead.\n",
" \n",
"/home/rflamary/.local/lib/python3.5/site-packages/ipykernel_launcher.py:9: DeprecationWarning: `imread` is deprecated!\n",
"`imread` is deprecated in SciPy 1.0.0.\n",
"Use ``matplotlib.pyplot.imread`` instead.\n",
" if __name__ == '__main__':\n"
]
}
],
"source": [
"def im2mat(I):\n",
" \"\"\"Converts and image to matrix (one pixel per line)\"\"\"\n",
" return I.reshape((I.shape[0] * I.shape[1], I.shape[2]))\n",
"\n",
"\n",
"square = spi.imread('../data/square.png').astype(np.float64)[:, :, 2] / 256\n",
"cross = spi.imread('../data/cross.png').astype(np.float64)[:, :, 2] / 256\n",
"triangle = spi.imread('../data/triangle.png').astype(np.float64)[:, :, 2] / 256\n",
"star = spi.imread('../data/star.png').astype(np.float64)[:, :, 2] / 256\n",
"\n",
"shapes = [square, cross, triangle, star]\n",
"\n",
"S = 4\n",
"xs = [[] for i in range(S)]\n",
"\n",
"\n",
"for nb in range(4):\n",
" for i in range(8):\n",
" for j in range(8):\n",
" if shapes[nb][i, j] < 0.95:\n",
" xs[nb].append([j, 8 - i])\n",
"\n",
"xs = np.array([np.array(xs[0]), np.array(xs[1]),\n",
" np.array(xs[2]), np.array(xs[3])])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Barycenter computation\n",
"----------------------\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ns = [len(xs[s]) for s in range(S)]\n",
"n_samples = 30\n",
"\n",
"\"\"\"Compute all distances matrices for the four shapes\"\"\"\n",
"Cs = [sp.spatial.distance.cdist(xs[s], xs[s]) for s in range(S)]\n",
"Cs = [cs / cs.max() for cs in Cs]\n",
"\n",
"ps = [ot.unif(ns[s]) for s in range(S)]\n",
"p = ot.unif(n_samples)\n",
"\n",
"\n",
"lambdast = [[float(i) / 3, float(3 - i) / 3] for i in [1, 2]]\n",
"\n",
"Ct01 = [0 for i in range(2)]\n",
"for i in range(2):\n",
" Ct01[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[0], Cs[1]],\n",
" [ps[0], ps[1]\n",
" ], p, lambdast[i], 'square_loss', # 5e-4,\n",
" max_iter=100, tol=1e-3)\n",
"\n",
"Ct02 = [0 for i in range(2)]\n",
"for i in range(2):\n",
" Ct02[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[0], Cs[2]],\n",
" [ps[0], ps[2]\n",
" ], p, lambdast[i], 'square_loss', # 5e-4,\n",
" max_iter=100, tol=1e-3)\n",
"\n",
"Ct13 = [0 for i in range(2)]\n",
"for i in range(2):\n",
" Ct13[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[1], Cs[3]],\n",
" [ps[1], ps[3]\n",
" ], p, lambdast[i], 'square_loss', # 5e-4,\n",
" max_iter=100, tol=1e-3)\n",
"\n",
"Ct23 = [0 for i in range(2)]\n",
"for i in range(2):\n",
" Ct23[i] = ot.gromov.gromov_barycenters(n_samples, [Cs[2], Cs[3]],\n",
" [ps[2], ps[3]\n",
" ], p, lambdast[i], 'square_loss', # 5e-4,\n",
" max_iter=100, tol=1e-3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Visualization\n",
"-------------\n",
"\n",
"The PCA helps in getting consistency between the rotations\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.PathCollection at 0x7fcec1ef00b8>"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAJDCAYAAABHZBNLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3W+sJVWd//vPt7tp+B3ujUJDEIU+Db/LiJiZOyMnRsdk4h0xwX4AOjoJpjMDE80ZW41Pb5tO7iQkZJx5cpOJZMgJMfT86CBK4rX9TU8IiMQHE5BDAvLvIg0B7B7QFm9MTE9whO99sOs0+5yz965VVatqrap6v5KVs//U2atO72+v9a21VlWZuwsAAABp7Ei9AwAAAGNGMgYAAJAQyRgAAEBCJGMAAAAJkYwBAAAkRDIGAACQUJRkzMy+ZWa/NLOn57xvZvZPZnbCzH5qZh+KUS+GgxhCDMQRmiKGkEKskbG7JF2/4P1PSbqqKKuS/jlSvRiOu0QMobm7RByhmbtEDKFjUZIxd/+xpF8v2ORGSf/iE49IereZXRqjbgwDMYQYiCM0RQwhha7WjL1P0s+nnp8sXgNCEUOIgThCU8QQotuVegemmdmqJsO+Ov/886+9+uqrE+8R2vb444//yt0vjvmZxNG4tBFDEnE0NrRFaKpJDHWVjJ2SdPnU88uK1zZx9zVJa5K0srLi6+vr3ewdkjGzVwI3DYohiTgamwoxJBFHmIO2CE1VbIs26Wqa8pikvy7OQvmIpN+4+2sd1Y1hIIYQA3GEpoghRBdlZMzM7pH0cUkXmdlJSX8n6RxJcvc7JB2XtF/SCUlnJP1NjHoxHMQQYiCO0BQxhBSiJGPu/vmS913SV2LUhWEihhADcYSmiCGkwBX4AQAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICESMYAAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICESMYAAAASipKMmdn1Zva8mZ0ws0Mz3r/FzE6b2RNF+WKMejEsxBGaIoYQA3GErjVOxsxsp6TbJX1K0jWSPm9m18zY9F53/+Oi3Nm0XgwLcVTP0aPSvn3Sjh2Tn0ePpt6jdIghxEAcRUQDFSzGyNiHJZ1w95fc/XeSvi3pxgifi3Ehjio6elRaXZVeeUVyn/xcXR11e0cMIQbiKAYaqEpiJGPvk/Tzqecni9e2+qyZ/dTM7jOzyyPUi2Ehjio6fFg6c2bza2fOTF4fKWIIMRBHMdBAVdLVAv4fSNrn7n8k6QFJR2ZtZGarZrZuZuunT5/uaNfQI8TRlFdfrfY6JAXGkDSeOEIttEVlqjZQI5/SjJGMnZI0fVRwWfHaWe7+hru/WTy9U9K1sz7I3dfcfcXdVy6++OIIu4YeIY4q2rs3/PWRtHPRYqjYdhRxhG1oi2Ko2kCNfEozRjL2mKSrzOwKM9st6SZJx6Y3MLNLp57eIOm5CPViWIijim67TVpa2vza0tLk9WkjaueIIcRAHMUQ2kBJTGkqQjLm7r+X9FVJ92sSkN9x92fM7FYzu6HY7Gtm9oyZPSnpa5JuaVpvciMZaujKaOOogQMHpLU1aXlZMpv8XFubvD5tLO0cMYQYiKNIQhsoiTUXkuTuWZZrr73Ws3X33e5LS+6TgYZJWVqavD69zfKyu9nk5/R7OEvSuo81jjpitjlUN4pZ6j2Lo+0YcuJoFGiLOjCvX1xent1ILS+n29camsQQV+Cvo2yoIXReiNE1dKDK0g2gDpoylFrUL1aZ0hwokrE6yoZUQ+aFSNhQUUgozNqGdg5tGtGaRDSxqF+sMqU5VHWH1NouWQ/plg2phswLhQzLjmA6VEwNBAkNhXnb9DxMFmo7hnxAcdSGgcww0Ra1bejrJbxZDDEyVkfZUEPIvFDIgkWmQ1EIGWwtO/B8+WXp7bcnP8d0wDlUufy3Zu01NpkXmKyXWIhkrI6yIdWQeaEYCVvM6VBkLaTDo1Mcj5z+W9PH4izWhdVGMlbXoqGGkPnvGAlbjNE19EJIh1e1U8xlZAXV5fTfmj4WZ7EurDaSsbaUzQvFSNhiTYcieyEdXpVOMaeRFVQX+791k8ScPhZnlQXm1n5R4ohwQ93FZm2XrBc7drkaelFdIau6M19dKxbNBgsJu9DQzDwsKmk7hjzDOIr5/YU0I2NAWxRBlcAcYOA1iaFWG7AmJdvADQ2gsl4xVkIXUk/GAU8DmMaQTmwaejI26794zP/WXSfmuZ7ZO/i2qIt/+CqBOaQjwgLJWJdiXJIiVkIXKtfWz0fQAAbq+isaUjs45GSsi8uVdJmY53xsOOi2qMt/+NDAHNIRYYFkrEsxriHGNcbOGnQDGKiNdrLnA6aVDDkZa5I05zhlnfNBwKDbopT/8AO/BdI0krEuhQRQWcIWI6HrenStJYNuAAPFbpMGEhrBhpyM1R08qJJsd5mY5zwYMui2qOo/fMxZmUVDu0M5IiyQjHUpxqL5GAldrNG1xAbdAAaq00EtaisHeMC50JCTsXnf5Z49i/vKqjHQtO/NcRSuqkG3RakW1ocMKgzhiLBAMta1pnNAMRK6WLdcSmzQDWCgOh3novDJefShDUNOxmZ917t3u59zzuLmI9d1YDkfHw66LUq1sH5kjRHJWI6ank1Z9p8nxuhaBgbdABZir9+KMfDa1t+SwpCTMfft/+Z79pR/vzmvA9v69xw8mEdMDb4tanNh/YjWhS1CMtalLnujRXXFGF3LwNAbwDbWb5W1lW2NPuQ6qjH0ZGyrkL4yxnfVxUlxOcXU0NuibWIlUCNbF7YIyVhXclsZPYBT5obeALZxNlzocsFFeXyd8Mw1tx9bMjbve9j4LqZXQ9SNgTZmtWbVmVNMDb0t2iRmAjWydWGLkIx1Jdai+VxG1zIw9AawjbPhmuTYTX4311nvsSVjs77DKt9n7EH10CZv1jbz/oYUMTX0tmiTmAlUqjM1M0Qy1pWcLkkxkIAeegNY98i/rYPNJqNqOY1iTBtbMua+eFSp7DuJtdx0Ok727JmUqmd37tyZT0wNvS3aJGYCFTtz7zGSsa7EaMVijK7lNl3awNAbwLptT9NRqHlffZP1Zrm2o2NMxjZUTZoWJXBtHFOW7efG7+UQU4Nsi2IcWcXqj6rW20PJkzFJ10t6XtIJSYdmvH+upHuL9x+VtK/sM7Ns/GKM7+d0Bf8MTAfvUONoVntYlic3XWs276tvGlo55vdtx5BnEkez1Ema5jVBVZqPqvG5aPtcYmpwbVGsI6umi1Snxb6oYmaSJmOSdkp6UdKVknZLelLSNVu2+bKkO4rHN0m6t+xzc238Gq98jTG6NpBrjLm/E7xjiqMma2xC2qGyjm9o1yhrO4Y80zhyr9/cbP2eqw6s15nlyv3YcHBtUawjq5SXuuhD4ExJnYx9VNL9U8+/LunrW7a5X9JHi8e7JP1Kki363FwbvyBlpzHlcAX/TEw1gKOJo9D2qO4BYchUZIzlH7loO4Y80zjaUCdp2vhOF8VWjDipsq4stcG1RbHWhcVMoKomVz1rkFInY5+TdOfU87+S9M0t2zwt6bKp5y9KumjR52bZ+MUaLm06uhb7dKiEphrA0cRRrDy5jYX2PTsQdff2Y8gzjaMQdQdHQpqg3bs3v797d7zR3RQG1xbFWhcWO4Gq0o/27EzNwSRjklYlrUta37t3b2v/YLV0vWg+JGFrktBloo0GMHUctbkebLqOthbaVwnhHJZztNWJpo6jGA4enB1rBw82W1t4993bb8l0zjm9PB48a3BtUcx1YakSqFgJZUdSJ2PDGNItE7qIkWuMVaKBTQ20vR5sQw4L7TNo+9y9/RjyHNujQIviZNF7MU4K78lKibOG1ha5e3vrwrpKoGImlB1InYztkvSSpCv0zmLHD27Z5ivavNjxO2Wfm13jl9M1xgZkqgEcRBy1vR5sQ1uXvqgig7bP3duPIc+xPQq0KE4WvVf23Q7oHKKzhtYWVZJzAtXmiQaRJU3GJvVrv6SfFUO1h4vXbpV0Q/H4PEnf1eQ04J9IurLsM7ML3BiHgn0cXWvZdPAOIY5itwe5rAubtR8ZtH3u3n4MeY7tUaC6I2NlMRKrKcvJ0NqiSvqUQLXRKEaSPBlro2QXuDEWzTO6tk2T4A0puY6MhWhrXVjVfZxX1549yds+d28/hjzH9ihQkxha1IQMrBly9+G1RZX1IYFqc7FsBCRjXWm6aD6n0bVMDK0BjNlJtbUurGpbOm8/9uzJI8xIxsoTpzrvNamzj4bWFrUmZQKVw2LZBUjGctL0cHJkizWG2AA2zdk3tLUurGp4LNqPHDrksSdjsY69QuI29XfdpiG2Ra1ImUDlsjZiDpKxPslhdC0jY2wAQ5OhttaFxW5LUxt7Mhbj+4m5JruvxtgW1ZYqgcq8MSIZ60JI8MU6dGx7dC0jY2wAQ9umNteFVQnV3DvisSdjMfq6snip0qT0dQRtjG1RVDHXhTVdrJgIyVjbQgKgyxWtsebBMjDGBjBWx7bovZhrbMveS23syViMY6+yeOniACK1MbZF0cQaig/tazNtjEjG2hbS2sVaWN/F6FpGxtgAxuiwYsxmx96nVMaejMX47mKNjPVoUH6bMbZF0cQaiu9zAHmzGGq1AWtSsgrckMNCLltRy1gbwKZfYaxQCv28nI09GXNvHk+x1oz1aLnqNmNti6KI9cXHHtLvGMlY22KNjHFR2G1oAOsJabMGdJLSQmNKxto84zH0syX3nTvfaZp6ulx1G9qiEosCJNYX3/MhfZKxtsVaM1YWaF2OrmWCBrCe2J0enWj+cZTDGY857ENbaIsW6OqL7/mQPslYF2KcTRljoU+M0bWM0ADO1vU5GnSi+cdRzDMe29oH914NzG9CW7RAl198j4f0Scb6pOlpuzFG1zJCA7hdqqWDdKJ5x1GsMx7b3IcyOccYbdECufYrmQ0+kIx1oatWpIvRtYzQAG7Xs68wubEkY30ZGZsn99FX2qIF6nzxXfSZmQUVyVjbcls033R0LSM0gNtVPQjNebShC2NJxnJYrzWrDjP3gwfLfzf3gwzaogWqBleX/VBGDSDJWNu6uiRFijn3xGgAt6vSaY3waijbjCUZc++uCVnk4MHtBwwh/WyuM10baItKVAmu3DPvlpCMta2LS1KMtFelAdyuykFlrKuh9NmYkrFZum4S6vazuffPtEUR5Z55t4RkrG1dXJIiVq/as2SNBnC20K8xRuj13ZiTsRSJdt1+NveDAtqiirq47ljPkIy1rYtLUozwGmPuNIBNxQi9vhtzMpaiz2u6iD/XY0XaogpyWMCYoSYxtEMod+CAtLYmLS9LZpOfa2uT1yXpttukpaXNv7O0NHl9w969sz974/Wy9yXp1Vdnb7Px+uHD0pkzm987c2byOgbh6FFp3z5px47Jz/37m4ce+qusSWhDSHM3z4ED0ssvS2+/Pfm50YSiZ8r6mrI+E9vVzeLaLr07imh6SYoRXmPMnaPRKuaFyMGD3V4gNjdtx5BnHEepZoNyHuGqi7aogh72NV1oEkOtNmBNyqACd0PTU6EGdo0xdxrAKoY6PdTUmJOxoSfaXaItmmFew9HDvqYLJGNjMqBrjLnTAFbBwehsY07G3IedaHeJtmiLRf1JD/uaLjSJoUZrxszsQjN7wMxeKH5eMGe7t8zsiaIca1Ln6C1adNHTeXriKEzo2q+t68qOHm17z9IbcwyxDiueMcfRNovWhfW0r8lZ0wX8hyT90N2vkvTD4vks/+nuf1yUGxrWiUX62TITRwFCFk4fPSqtrkqvvDI5XH3llcnzESRkxBBiII42lJ0dMquvGeORYCRNk7EbJR0pHh+R9OmGn4dxIo4ChByMjviEWmIIMRBHG6qehj3iI8EYmiZjl7j7a8Xj1yVdMme788xs3cweMbPxBjfmIY4ClQ18prjUQSaIIcRAHG2oeg2TER8JxrCrbAMze1DSe2a8telf2N3dzHzOxyy7+ykzu1LSQ2b2lLu/OKOuVUmrkrSXiyANynXXXafXX3991lvvnn5CHDWzd+/kgHTW633XZQxJ446jIaMtCrRxpHf48ORobu/eSSI2b+nLiI8EY7DJCQA1f9nseUkfd/fXzOxSSQ+7+/tLfucuSf/T3e9btN3Kyoqvr6/X3jf0g5k9Lul/FXEUxcZMwfQB6tLSsNfWth1D0vjiaIxoixrat2/2keDy8mQYfwTM7HF3X6nzu02nKY9Jurl4fLOk72/dwMwuMLNzi8cXSfqYpGcb1othIY4iGfFJTsQQYiCO6mpyawY0Tsa+IemTZvaCpOuK5zKzFTO7s9jmA5LWzexJST+S9A13J3AxjTiKqJ8n1DZGDCEG4qiuER8JxtBomrJNgx/ShaRmw7ohiKPhazuGJOJoDGiL0FTKaUoAAAA0QDIGAACQEMkYAABAQiRjAAAACZGMAQAAJEQyBgAAkBDJGAAAQEIkYwAAAAmRjAEAACREMgYAAJAQyRgAAEBCJGMAAAAJkYwBAAAkRDIGAACQEMkYAABAQiRjAAAACZGMAQAAJEQyBgAAkBDJGAAAQEIkYwAAAAk1SsbM7C/N7Bkze9vMVhZsd72ZPW9mJ8zsUJM6MTzEEZoihhADcYRUmo6MPS3pLyT9eN4GZrZT0u2SPiXpGkmfN7NrGtaLYSGO0BQxhBiIIySxq8kvu/tzkmRmizb7sKQT7v5Sse23Jd0o6dkmdWM4iCM0RQwhBuIIqXSxZux9kn4+9fxk8RpQBXGEpoghxEAcIbrSkTEze1DSe2a8ddjdvx9zZ8xsVdJq8fRNM3s65udXcJGkX1FvVH8g6ZwZr38wdkWZxFGq7zJl3W3X21kMSaOPoyHHL23R8OtOVe/76/5iaTLm7tfV/fDCKUmXTz2/rHhtVl1rktYkyczW3X3uAso2pap7bPVu1B24aa/iKPW/6Zj+5jZiSBp3HI01fgM3pS3KvO4exNA2XUxTPibpKjO7wsx2S7pJ0rEO6sWwEEdoihhCDMQRomt6aYvPmNlJSR+V9K9mdn/x+nvN7LgkufvvJX1V0v2SnpP0HXd/ptluY0iIIzRFDCEG4gipND2b8nuSvjfj9f+QtH/q+XFJxyt+/FqTfWsoVd1jq1eS1gYaR8Rvh/W2HEPSCP9NE9Wbsm7aouHU3bt6zd1j7ggAAAAq4HZIAAAACWWTjKW8DYWZXWhmD5jZC8XPC+Zs95aZPVGU2gs2y/4GMzvXzO4t3n/UzPbVrativbeY2empv/GLker9lpn9ct5p3TbxT8V+/dTMPtSgriRxNJYYCqybOKpf7yjiiBjatF2vY6j4LOJo8/vV48jdsyiSPqDJNToelrQyZ5udkl6UdKWk3ZKelHRNhLr/UdKh4vEhSf8wZ7vfRqir9G+Q9GVJdxSPb5J0b0f13iLpmy18t38m6UOSnp7z/n5J/ybJJH1E0qN9i6MxxBBxRBzRFhFDxFE7cZTNyJi7P+fuz5dsdvY2FO7+O0kbt6Fo6kZJR4rHRyR9OsJnzhPyN0zvz32SPmG2+P4ckepthbv/WNKvF2xyo6R/8YlHJL3bzC6tWVeqOBpDDIXW3QriKDraou2IoeqIo+0qx1E2yVigtm5DcYm7v1Y8fl3SJXO2O8/M1s3sETOrG+Ahf8PZbXxyGvVvJO2pWV+VeiXps8Ww6n1mdvmM99vQ9e1F2qhvDDEUWrdEHNU1hjgihtqtr8sYkoijWSp/r40ubVGVdXhrpSp1Tz9xdzezeaeYLrv7KTO7UtJDZvaUu78Ye18T+oGke9z9TTP7W02OZP488T5tkyqOiKFgxFHNeqefjDyOiKGa9U4/GXkMST2JI6njZMw7vLVSlbrN7Bdmdqm7v1YMJf5yzmecKn6+ZGYPS/oTTeasqwj5Gza2OWlmuyS9S9IbFeupXK+7T9dxpyZrD7pQ9TY1SeKIGAqrmzhajDgihurWF1JvxzEkEUezVP5e+zZN2dZtKI5Jurl4fLOkbUc0ZnaBmZ1bPL5I0sckPVujrpC/YXp/PifpIS9WBTZQWu+WOe0bNLm6dBeOSfrr4gyUj0j6zdQwexvaiKMxxFBQ3cRRI2OII2LoHX2PIYk4mqV6HHnkswzqFkmf0WRe9U1Jv5B0f/H6eyUdn9puv6SfaZLBH45U9x5JP5T0gqQHJV1YvL4i6c7i8Z9KekqTMzaekvSFBvVt+xsk3SrphuLxeZK+K+mEpJ9IujLS31lW799Leqb4G38k6epI9d4j6TVJ/1V8x1+Q9CVJXyreN0m3F/v1lOaceZRzHI0lhogj4ogYIoaIo/hxxBX4AQAAEurbNCUAAMCgkIwBAAAkRDIGAACQEMkYAABAQlGSMevw5qsAALSF/gwpxBoZu0vS9Qve/5Skq4qyKumfI9ULAEBMd4n+DB2Lkox5hzdfBQCgLfRnSKGrNWNd33wVAIA20J8huk7vTVnGzFY1GfbV+eeff+3VV1+deI/Qtscff/xX7n5x6v0AgNjo08alSX/WVTIWdNNMd1+TtCZJKysrvr6+3s3eIRkzeyX1PgBABcE3gaZPG5cm/VlX05Rd33wVAIA20J8huigjY2Z2j6SPS7rIzE5K+jtJ50iSu98h6bgmN/Q8IemMpL+JUS8AADHRnyGFKMmYu3++5H2X9JUYdQEA0Bb6M6TAFfgBAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICESMYAAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABKKkoyZ2fVm9ryZnTCzQzPev8XMTpvZE0X5Yox6AQCIjT4NXdvV9APMbKek2yV9UtJJSY+Z2TF3f3bLpve6+1eb1gcAQFvo05BCjJGxD0s64e4vufvvJH1b0o0RPhcAgK7Rp6FzMZKx90n6+dTzk8VrW33WzH5qZveZ2eUR6gUAIDb6NHSuqwX8P5C0z93/SNIDko7M2sjMVs1s3czWT58+3dGuAQBQCX1aREePSvv2STt2TH4ePZp6j7oXIxk7JWn6qOCy4rWz3P0Nd3+zeHqnpGtnfZC7r7n7iruvXHzxxRF2DQCASujTOnT0qLS6Kr3yiuQ++bm6Or6ELEYy9pikq8zsCjPbLekmScemNzCzS6ee3iDpuQj1AgAQG31ahw4fls6c2fzamTOT18ek8dmU7v57M/uqpPsl7ZT0LXd/xsxulbTu7sckfc3MbpD0e0m/lnRL03oBAIiNPq1br75a7fWhirJmzN2Pu/sfuPt/d/fbitf+ryJo5e5fd/cPuvv/7u7/h7v/vzHqHauy+XXm3wGgPvq07uzdW+31oeIK/JkITaDK5teZfwcAVJXqIP6226Slpc2vLS1NXh8TkrEMVEmgyubXmX8HAFSR8iD+wAFpbU1aXpbMJj/X1iavj4m5e+p9mGllZcXX19dT70Yn9u2bBP9Wy8vSyy9vfm3Hjsl/lq3MpLffLn8/N2b2uLuvpN4PAGhTzn1alT4I8zXpzxgZi6TJEG+VBYxl8+vMvwMAquhiEf2iPpJ1ziRjUTQd4q2SQJXNrzP/DgCoou2D+EV9JOucJ0jGImi6TqtKAlU2v878OwCgirYP4hf1kaxzniAZiyB0iHfeUGzVBOrAgck8/ttvT35u3a7sfQAANrR9EL+oj2wyRTqk6U2SsQhChnjLhmKrJFBDCkAAQPvK+o2QPqhu37Ooj6w7RTq46U13z7Jce+21nsrdd7svL7ubTX7efXf59ktL7pOQmJSlpc2/t7y8+f2Nsrxcrd6QuvpEkytaJ483CoVCabOk7tMW9Rtt9z2Lfrfu54b0qV1r0p8lD9B5JVXg1g2MsmA2mx04ZtXqzTEAmyAZo1AoYygpk7FF/UZXfc+iPrLqAIh7eZ+aAslYRG0lO2WfG1pvaADWCe4USMYoFMoYSspkbFG/EbvvmdZmP5TjwEST/ow1Y1u0db2VsrNVQuuNsT4NADAei/qNmH3PtLb7oaFdxolkbIu2rrdSdrZKaL0hAcipwgCADYv6jZh9z7S2+6HBXcap7pBa22Voa8Zi1tt0fVpOxDQlhUIZQUk5Tek+v9+I2fdMa9IP9WWZzVZN+rPkATqvpD7zJPbZlG3UO0+Oc+nzkIxRKJQxlNTJ2CIbfY/kvnOnb1rcX1fZSQOLFvP39YoBJGMdWBQ8uSU/fQpmkjEKhTKGkluftlWMy1+EfN7Bg4vrya0/rYJkrGVlQVplODb0ei5NR8j6MsxLMkahUMZQcurTZolx+Yut/c7Bg9v7obJkq0/LbLYiGWtZrMtShAR0n0a1YiAZo1AoYyg59WmzNL38RWjfVZZsjXVkjLMpA5Sd+ht6lknI2SWhZ6BwSyQAQCxNL38R2neVnb05tEtWhCIZC1AWPKGn2IYEdMg2XEcMABBT08tfhF6vrCzZGtwlK0LVHVKbLpKul/S8pBOSDs14/1xJ9xbvPyppX9ln5jSkG2vqMGT4NdY2G/ud+7oxMU1JoVAyK0Pv0+ZpcvmLKtOLfeib6mjSn8UI2p2SXpR0paTdkp6UdM2Wbb4s6Y7i8U2S7i373NwCN9ai+hhrxkIWOPZl7RnJGIVCyamMpU+rqqwP7Euf06bUydhHJd0/9fzrkr6+ZZv7JX20eLxL0q8k2aLP7WPgxjpTsmybmKNnqZGMUSiUnAp9Wn2hgxaMjG0vMdaMvU/Sz6eenyxem7mNu/9e0m8k7YlQdzZC13EdOCC9/LL09tuTn7Pmwcu2CVng2NY9NgFg4OjTSsw7gSykf2PN82xZLeA3s1UzWzez9dOnT6fenU3Kzl7s8n6QIQsc27rHJgAgTM59Wl0hydSi/pJ7J88WIxk7JenyqeeXFa/N3MbMdkl6l6Q3tn6Qu6+5+4q7r1x88cURdi2OkOCLNRIVesmKGKNnAIBtBt+nNVGWTJX1l1X6ylFdwqnu/OZG0WS+/CVJV+idxY4f3LLNV7R5seN3yj43p/n1WGvMkqF6AAAgAElEQVS0ul4A2Yd5ebFmjEKhZFTG0Kc10fSirTEvkp6bJv1ZrODdL+lnmpyBcrh47VZJNxSPz5P0XU1OA/6JpCvLPjOnwI1x9mLsU4On9SHpmodkjEKh5FaG3qc10fR2RqFJVl9OQpuWPBlro+QUuDGu6xXyGXXuydXHo4dpJGMUCmUMJac+rYmyPifGLJF7P+9R2aQ/y2oBf65C118tWscVMk9eZ9E9iyEBAF0pO4EspL8MOetybCehkYwFiHF7hpDAqrPonktYAABiKls4vyiZqtJfLqpnbCehkYwFCsnkFwk9WpgXxPOCdmxHDwCA9jS9dEWseqoOgvT+zMu685ttl6HMr0+ru9B+0Rz9vPcOHuzHon6xZoxCoYyg9KVPK1vzFeNktZB6tlrUf+aydrpJf2aT38/PysqKr6+vp96NLOzbNzlq2Gp5eTJKd/ToZI3Yq69ORsT275eOHNm8lmxpqfrUahfM7HF3X0m9HwDQpr70aTt2TNKZrcwmM0Nl/VHZ+6H1TNsYRZvXp4XW2bYm/RnTlA2UDYvGGjYtWxe2dQr1+HEW9QMAqitb+lLWH4WuY66yxKbsRLUhrJ0mGaupbL479P5bIQlb1XVhQwhMAED3ytY3l/VHof1VlQX6ZX3aINZO153fbLvkPr8e4yrDofPcVefD+3SxPLFmjEKhjKDk3qdNa7I+q0p/FbqOuuk6tq406c+SB+i8knvgll2QLuSCdVWSpiqL/3MJzBAkYxQKZQwl9z6tipBb+8U8gSykT8vhTjQkYwnEGBlr8wrDOQRmCJIxCoUyhpJTn5ZL/1B1kCGHfV6kSX/GmrES89Z0lc13h8yHtznP3fS6aACA4Ym5nrnJSWqh+7Fh8H1a3Syu7ZLDUUTI3HiTodqmnz8EYmSMQqGMoOTQp7nHW8/cdDlMn9Y2h2rSnyUP0Hklh8DtIljmJVwx1n31IZkjGaNQKGMoOfRp7vHWMzftH/t4I/AyTfozpikX6OISEfOGXhddVyV0+LjKEDAAYPhClseE9H1N+8dBXI4iIpKxBVIGy7yA3kiqypKssovkAQDGJ9Z65qb949huBF6GZGyBNoIldMHjvIDeuTMsyeLCrwCArUJuwB3S9zXtH6veCHzw6s5vtl1ymV+Pue6q6sXwZm07a4591jx7XxZHijVjFAplBCWXPi1USN/X1rrkPqx3nqVJf5Y8QOeVvgVuiBh3qQ/9jL5c+JVkjEKhjKH0pU9LnQj1pe+apUl/xjRlh6pOHc5a3B86NMwQMACgihxO/BrremeSsQ7FOCGgSpI1+IvkAQCiySERGut6Z5KxDsU6IYAkCwAQWw6J0FgvedEoGTOzC83sATN7ofh5wZzt3jKzJ4pyrEmdfcbUIQDka+x9Wg6J0FgvedF0ZOyQpB+6+1WSflg8n+U/3f2Pi3JDwzp7jVEtAMjWqPu0HBKhsQ5aNE3GbpR0pHh8RNKnG34eAACpjLpPyyURGuOgxa6Gv3+Ju79WPH5d0iVztjvPzNYl/V7SN9z9/2lYLwAAsY2+TztwYBzJT25KkzEze1DSe2a8ten8Cnd3M/M5H7Ps7qfM7EpJD5nZU+7+4oy6ViWtStLeoa/WAwB0jj4NOSpNxtz9unnvmdkvzOxSd3/NzC6V9Ms5n3Gq+PmSmT0s6U8kbQtcd1+TtCZJKysr8/4TAABQC30actR0zdgxSTcXj2+W9P2tG5jZBWZ2bvH4Ikkfk/Rsw3oBAIiNPg1JNE3GviHpk2b2gqTriucysxUzu7PY5gOS1s3sSUk/0mR+ncAFAOSGPg1JNFrA7+5vSPrEjNfXJX2xePzvkv6wST0AALSNPg2pcAV+AACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICESMYAAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACChRsmYmf2lmT1jZm+b2cqC7a43s+fN7ISZHWpSJwAAbaBPQypNR8aelvQXkn48bwMz2ynpdkmfknSNpM+b2TUN6wUAIDb6NCSxq8kvu/tzkmRmizb7sKQT7v5Sse23Jd0o6dkmdQMAEBN9GlLpYs3Y+yT9fOr5yeI1AAD6hj4N0ZWOjJnZg5LeM+Otw+7+/Zg7Y2arklaLp2+a2dMxP7+CiyT9ino78f5E9QIYoRH2aSnb97H1abX7s9JkzN2vq/vhhVOSLp96flnx2qy61iStSZKZrbv73AWUbUpV99jq3ag7Rb0AxmlsfVrq9n1Mf3OT/qyLacrHJF1lZleY2W5JN0k61kG9AADERp+G6Jpe2uIzZnZS0kcl/auZ3V+8/l4zOy5J7v57SV+VdL+k5yR9x92fabbbAADERZ+GVJqeTfk9Sd+b8fp/SNo/9fy4pOMVP36tyb41lKrusdWbum4AOGugfdoY2/fe1WvuHnNHAAAAUAG3QwIAAEgom2Qs5W0ozOxCM3vAzF4ofl4wZ7u3zOyJotResFn2N5jZuWZ2b/H+o2a2r25dFeu9xcxOT/2NX4xU77fM7JfzTuu2iX8q9uunZvahGPUCQCqp+rSu+7Pis+jTNr9fvU9z9yyKpA9oco2OhyWtzNlmp6QXJV0pabekJyVdE6Huf5R0qHh8SNI/zNnutxHqKv0bJH1Z0h3F45sk3dtRvbdI+mYL3+2fSfqQpKfnvL9f0r9JMkkfkfRo6nikUCiUJiVVn9Zlfxb6N9Cnlfdp2YyMuftz7v58yWZnb0Ph7r+TtHEbiqZulHSkeHxE0qcjfOY8IX/D9P7cJ+kTVnJ/jkj1tsLdfyzp1ws2uVHSv/jEI5LebWaXdrFvANCGhH1al/2ZRJ82S+U+LZtkLFBbt6G4xN1fKx6/LumSOdudZ2brZvaImdUN8JC/4ew2PjmN+jeS9tSsr0q9kvTZYlj1PjO7fMb7beD2IgDGqI22r8v+TKJPm6Xy99ro0hZVWYe3oahS9/QTd3czm3eK6bK7nzKzKyU9ZGZPufuLsfc1oR9Iusfd3zSzv9XkSObPE+8TAGQpVZ9GfxasN31ap8mYd3gbiip1m9kvzOxSd3+tGEr85ZzPOFX8fMnMHpb0J5rMWVcR8jdsbHPSzHZJepekNyrWU7led5+u405N1h50ofb3CgCppOrTMurPJPq0WSp/r32bpmzrNhTHJN1cPL5Z0rYjGjO7wMzOLR5fJOljkp6tUVfI3zC9P5+T9JAXqwIbKK13y5z2DZpcXboLxyT9dXEGykck/WZqmB0AhqqNPq3L/kyiT5ulep8W+yyDBmcnfEaTedU3Jf1C0v3F6++VdHzLWQo/0ySDPxyp7j2SfijpBUkPSrqweH1F0p3F4z+V9JQmZ2w8JekLDerb9jdIulXSDcXj8yR9V9IJST+RdGWkv7Os3r+X9EzxN/5I0tWR6r1H0muS/qv4jr8g6UuSvlS8b5JuL/brKc0584hCoVD6UlL1aV33Z/P+Bvq0an0aV+AHAABIqG/TlAAAAINCMgYAAJAQyRgAAEBCJGMAAAAJRUnGWrlpJgAAHaM/QwqxRsbuknT9gvc/JemqoqxK+udI9QIAENNdoj9Dx6IkY86NoAEAA0B/hhS6WjPGjaABAENAf4boOr03ZRkzW9Vk2Ffnn3/+tVdffXXiPULbHn/88V+5+8Wp9wMAYqNPG5cm/VlXyVjQTTPdfU3SmiStrKz4+vp6N3uHZMzsldT7AAAVBN8Emj5tXJr0Z11NU3IjaADAENCfIbooI2Nmdo+kj0u6yMxOSvo7SedIkrvfIem4Jjf0PCHpjKS/iVEvAAAx0Z8hhSjJmLt/vuR9l/SVGHUBANAW+jOkwBX4AQAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICESMYAAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICEoiRjZna9mT1vZifM7NCM928xs9Nm9kRRvhijXgAAYqNPQ9d2Nf0AM9sp6XZJn5R0UtJjZnbM3Z/dsum97v7VpvUBANAW+jSkEGNk7MOSTrj7S+7+O0nflnRjhM8FAKBr9GnoXIxk7H2Sfj71/GTx2lafNbOfmtl9ZnZ5hHoBAIiNPg2d62oB/w8k7XP3P5L0gKQjszYys1UzWzez9dOnT3e0awAAVEKfhqhiJGOnJE0fFVxWvHaWu7/h7m8WT++UdO2sD3L3NXdfcfeViy++OMKuAQBQCX1ah44elfbtk3bsmPw8ejT1HqURIxl7TNJVZnaFme2WdJOkY9MbmNmlU09vkPRchHoBAIiNPq0jR49Kq6vSK69I7pOfq6vjTMgaJ2Pu/ntJX5V0vyYB+R13f8bMbjWzG4rNvmZmz5jZk5K+JumWpvUCABAbfVp3Dh+WzpzZ/NqZM5PXx8bcPfU+zLSysuLr6+updwMtM7PH3X0l9X4AQJvo07bbsWMyIraVmfT2293vT1NN+jOuwJ+R0LnzkO2YhwcA5Gzv3mqvDxnJWCZC585DtmMeHgBQV1cH87fdJi0tbX5taWny+tgwTZmJffsmSdNWy8vSyy9X2y70s3LANCWAMehLn7ZxMD+9lmtpSVpbkw4caKe+w4elV1+djIjddls79XShSX9GMpaJ0LnzkO36NA9PMgZgDPrSp/XpYD43rBkbgNC585DtmIcHAMxSNgX56quzf2/e603qCt1mDEjGMhE6dx6yHfPwAICtQtYTxzqYZ31zRe6eZbn22mu9T+6+23152d1s8vPuu9v7jJDtYuxPFyStewbxRqFQKG2WHPq05WX3SdqzuSwvv7PN3Xe7Ly1tfn9pqXofElJXyDZ90qQ/Y81YBKELHoe0UDEW1owBGIMc+rTQ9cQx+qqhrW8OwZqxxEKuIsxwLAAgpdApyAMHJov133578rPOoAHrm6shGYsgZMFj6G0fWPAIAGjDrPXEu3dLv/1ttf4kpA9ifXNFdec32y45zK+HCpn3Npu9jdk724TM1YfO57NmjEKhUPIpufRp033Dnj3u55zjldaHVVlTNqT1zSGa9GfJA3ReySVwQ4QEZ6zFjF0uwOwCyRiFQhlDybFPq7OAvsmi+yElXrM06c+Ypiw0mfo7cGCyWH95ebLwcHl5++L9kOHYkOnOmFOiAIDxqnNNsbrXIWPd9GIkY4oTJGULHkMStlgLHmNetA8AMEx1FtDXXXTf9iBB39dSk4ypu5GksoQt1oJHzlABAJSZ1Z+cc87iBf11F923OUgwhFE3kjHFCZIYWXnI6FmsKVEAwLht7U/27Jn8fOON+UlNSB80qz9sc5BgEEtz6i42a7t0udixbEFi2aLDHBfM92WhpFjAT6FQRlByXMC/VYwr4s/rDw8ebK+fDLlaQRea9GeMjGnxSFLI8GfMa4jFEuOifQCA8QiZJSrrx+b1h8ePl4+ohdax1SCW5tTN4touXR9FzBtJ6vIaYov2o+o2fSFGxigUygjKEEbGQvqxsv6wjZmmXGanmvRnyQN0XsklcEMSrS6vD8ZFXykUCqV/JZc+bZGy/qVpXxfrmpzz9j11n0cy1qJYSVRuSV0uSMYoFMoYSi59WplFSU3TWaBYM025atKfRVkzZmbXm9nzZnbCzA7NeP9cM7u3eP9RM9sXo94uhJyZGOsaYlz0FQDSG3KfVmbReuOQfmxRfxjSfw1i/VcNjZMxM9sp6XZJn5J0jaTPm9k1Wzb7gqT/z93/N0n/t6R/aFpvV0ISrY3tml5DjIu+AkBaQ+/Tmgi9bNK8/jCk/xrrpZlijIx9WNIJd3/J3X8n6duSbtyyzY2SjhSP75P0CTOzCHV3oizRCjnzI9b1wbjoKwC0avB9Wl2hgxPzxJppGqS685sbRdLnJN059fyvJH1zyzZPS7ps6vmLki5a9Ll9ml+PuT4rxtmUrBmjUCiUemXsfVqIkD5o3vtNF9rnsFB/nib9WVaBK2lV0rqk9b1797b2DxZTjIvktSHngJ1GMkahUHIqY+/TypQd7Lc5GJD7QEOT/izGNOUpSZdPPb+seG3mNma2S9K7JL2x9YPcfc3dV9x95eKLL46wa/HMm4oMXZ8VMpUZeqG70GlRLvoKAJWNok9bZFEfU3aCWJUTyKpe3HXQJ6fVzeI2iqRdkl6SdIWk3ZKelPTBLdt8RdIdxeObJH2n7HNzGtJteqpu7OuH5XxkUJUYGaNQKBmVMfRpi5T1MWWXngi9NEWdviz3y1406c9iBe9+ST/TZKj2cPHarZJuKB6fJ+m7kk5I+omkK8s+M6fA7eIidqHTnblOi9ZFMkahUHIrQ+/TFinrY5q+H1pPnX1LLXky1kbJKXCb3t4hJJsPzfirHHWwZoxCoVDyKDn1aYuE9Hcx1ozVGeXKfWaoSX/GjcIDlF0qomx9Vqzrh4VuF3JzcwAAtgrp7xZdeiL00hRVLsG0sbbsr/5K+m//TdqzZ/ZnV12DlpW6WVzbJaejiKbZeNdrxnIfyp0mRsYoFMoISk592iJdjT7FXiedw6hZk/4seYDOK7kFbhfXRgmtI8a0aC5IxigUyhhKbn3aIm0uc5n+7D17JmVRPX1aT00yFklf1lmVySEoQ5GMUSiUMZQ+JWNbxRooaPMMyhwGIZr0Z6wZKwxpndVY7+0FAIgrtG8M2a7OdcJirqfOGclYoc2LycVaVBj6OaO9txcAIKrQvjFku9CLpE8LHVzo/SBE3SG1tkvXQ7pNhjjL7sMVuvhwKPebrEJMU1IolBGUnKYpqyzJiTlNWHcJTaxp0rY16c+SB+i80nXgNgmSRUlSrCv0h+5f6mCsimSMQqGMoeSSjFU9sI+5gH6ogwobSMYiqBskZQEY62gh5HP6GOgkYxQKZQwll2Ss6sBDG5eg6NOAQRUkY5HUCZKyJClWohXzlko5IRmjUChjKLkkY6EH9tN94cGD5ctoNvqfnTv9bL8zpEQrRJP+jAX8U8qupD9L2RkcIYsKQ84CCfmcOosjAQDjUdbfzDor8siRSV8zq2+c3l6S3nrrnb5pXh/a6yvlt6VuFtd2yeUookzo1fVjLM4v+xxGxigUCiXPkkufFmOd87S2pj37qEl/ljxA55VcAjdEjDnwWJ/RtyAnGaNQKGMoOfVpi/qbqlcWqLp9HwcNQjXpz2zy+/lZWVnx9fX11LvRO0ePTq7r8sor0s6dkyHj5eXFQ8Ypmdnj7r6Sej8AoE196dP27XtnynHa8vJkirLp9jt2TNKvrcwm06B91qQ/Y81YYqFz51Uu+Lqxvuyttyav9fluAgCA9mztW/bvr3bx1KoXW+37lfJbU3dIre2S05BuW9q6G32fhoHFNCWFQhlBybFPm9e3lJ09OetzQrfv43KaUE36M6YpEwod3h3yMDDTlADGIMc+rWrfEsvGcppXX52MiOW6jKaqJv3Zrtg7g3Chl6KoesmKvXtn/wcb/TAwAOCsVJdDOnBgGMlXTKwZS6itu9H3/oapAIDW5bB+i2uOTZCMJdTW3egPHJDW1iZDzWaTn2trHIkAAN6R+sB91gVmx3qyWaNkzMwuNLMHzOyF4ucFc7Z7y8yeKMqxJnUOSWjSVCe5qnM3AQAYs7H1aakP3A8fls6c2fzamTOT18em0QJ+M/tHSb9292+Y2SFJF7j7/zlju9+6+/9S5bNzXOyI+FjADyAX9Gnd6tPJZiFSXmfsRklHisdHJH264ecBAJAKfVqHclizloumydgl7v5a8fh1SZfM2e48M1s3s0fMjOAGAOSIPq1Dqdes5aT00hZm9qCk98x4a9Osrru7mc2b81x291NmdqWkh8zsKXd/cUZdq5JWJWnvGFNjAECr6NPysbE2bYjXHKuq6Zqx5yV93N1fM7NLJT3s7u8v+Z27JP1Pd79v0XbMr48Da8YA5II+DU2kXDN2TNLNxeObJX1/6wZmdoGZnVs8vkjSxyQ927BeAABio09DEk2TsW9I+qSZvSDpuuK5zGzFzO4stvmApHUze1LSjyR9w90JXABAbujTkESj2yG5+xuSPjHj9XVJXywe/7ukP2xSDwAAbaNPQypcgR8AACAhkjEAAICESMYAAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICESMYAAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAIKFGyZiZ/aWZPWNmb5vZyoLtrjez583shJkdalInAABtoE9DKk1Hxp6W9BeSfjxvAzPbKel2SZ+SdI2kz5vZNQ3rBQAgNvo0JLGryS+7+3OSZGaLNvuwpBPu/lKx7bcl3Sjp2SZ1AwAQE30aUulizdj7JP186vnJ4jUAAPqGPg3RlY6MmdmDkt4z463D7v79mDtjZquSVounb5rZ0zE/v4KLJP2Kejvx/kT1AhihEfZpKdv3sfVptfuz0mTM3a+r++GFU5Iun3p+WfHarLrWJK1Jkpmtu/vcBZRtSlX32OrdqDtFvQDGaWx9Wur2fUx/c5P+rItpysckXWVmV5jZbkk3STrWQb0AAMRGn4boml7a4jNmdlLSRyX9q5ndX7z+XjM7Lknu/ntJX5V0v6TnJH3H3Z9pttsAAMRFn4ZUmp5N+T1J35vx+n9I2j/1/Lik4xU/fq3JvjWUqu6x1Zu6bgA4a6B92hjb997Va+4ec0cAAABQAbdDAgAASCibZCzlbSjM7EIze8DMXih+XjBnu7fM7Imi1F6wWfY3mNm5ZnZv8f6jZravbl0V673FzE5P/Y1fjFTvt8zsl/NO67aJfyr266dm9qEY9QJAKqn6tK77s+Kz6NM2v1+9T3P3LIqkD2hyjY6HJa3M2WanpBclXSlpt6QnJV0Toe5/lHSoeHxI0j/M2e63Eeoq/RskfVnSHcXjmyTd21G9t0j6Zgvf7Z9J+pCkp+e8v1/Sv0kySR+R9GjqeKRQKJQmJVWf1mV/Fvo30KeV92nZjIy5+3Pu/nzJZmdvQ+Huv5O0cRuKpm6UdKR4fETSpyN85jwhf8P0/twn6RNWcn+OSPW2wt1/LOnXCza5UdK/+MQjkt5tZpd2sW8A0IaEfVqX/ZlEnzZL5T4tm2QsUFu3objE3V8rHr8u6ZI5251nZutm9oiZ1Q3wkL/h7DY+OY36N5L21KyvSr2S9NliWPU+M7t8xvtt4PYiAMaojbavy/5Mok+bpfL32ujSFlVZh7ehqFL39BN3dzObd4rpsrufMrMrJT1kZk+5+4ux9zWhH0i6x93fNLO/1eRI5s8T7xMAZClVn0Z/Fqw3fVqnyZh3eBuKKnWb2S/M7FJ3f60YSvzlnM84Vfx8ycwelvQnmsxZVxHyN2xsc9LMdkl6l6Q3KtZTuV53n67jTk3WHnSh9vcKAKmk6tMy6s8k+rRZKn+vfZumbOs2FMck3Vw8vlnStiMaM7vAzM4tHl8k6WOSnq1RV8jfML0/n5P0kBerAhsorXfLnPYNmlxdugvHJP11cQbKRyT9ZmqYHQCGqo0+rcv+TKJPm6V6nxb7LIMGZyd8RpN51Tcl/ULS/cXr75V0fMtZCj/TJIM/HKnuPZJ+KOkFSQ9KurB4fUXSncXjP5X0lCZnbDwl6QsN6tv2N0i6VdINxePzJH1X0glJP5F0ZaS/s6zev5f0TPE3/kjS1ZHqvUfSa5L+q/iOvyDpS5K+VLxvkm4v9uspzTnziEKhUPpSUvVpXfdn8/4G+rRqfRpX4AcAAEiob9OUAAAAg0IyBgAAkBDJGAAAQEIkYwAAAAlFScZauWkmRoUYQgzEEZoihpBCrJGxuyRdv+D9T0m6qiirkv45Ur0YjrtEDKG5u0QcoZm7RAyhY1GSMedG0GiIGEIMxBGaIoaQQldrxrgRNJoihhADcYSmiCFE1+m9KcuY2aomw746//zzr7366qsT7xHa9vjjj//K3S+O+ZnE0bi0EUMScTQ2tEVoqkkMdZWMBd00093XJK1J0srKiq+vr3ezd0jGzF4J3DT4xqvE0bhUiCGJOMIctEVoqmJbtElX05TcCBpNEUOIgThCU8QQoosyMmZm90j6uKSLzOykpL+TdI4kufsdko5rckPPE5LOSPqbGPViOIghxEAcoSliCClEScbc/fMl77ukr8SoC8NEDCEG4ghNEUNIgSvwAwAAJEQyBgAAkBDJGAAAQEIkYwAAAAmRjAEAACREMgYAAJAQyRgAAEBCJGMAAAAJkYwBAAAkRDIGAACQEMkYAABAQiRjAAAACZGMAQAAJEQyBgAAkBDJGAAAQEIkYwAAAAmRjAEAACREMgYAAJAQyRgAAEBCJGMAAAAJRUnGzOx6M3vezE6Y2aEZ799iZqfN7ImifDFGvRgW4ghNEUOIgThC13Y1/QAz2ynpdkmflHRS0mNmdszdn92y6b3u/tWm9WGYiCM0RQwhBuIIKcQYGfuwpBPu/pK7/07StyXdGOFzMS7EEZoihhADcYTOxUjG3ifp51PPTxavbfVZM/upmd1nZpdHqBfDQhyhKWIIMRBHXTl6VNq3T9qxY/Lz6NHUe5RMVwv4fyBpn7v/kaQHJB2ZtZGZrZrZupmtnz59uqNdQ48QR2gqKIak8cQR/WEttEVNHT0qra5Kr7wiuU9+rq6ONgBjJGOnJE0fFVxWvHaWu7/h7m8WT++UdO2sD3L3NXdfcfeViy++OMKuoUeIIzQVLYaKbQcfR/SHM9EWxTQv2z98WDpzZvO2Z85MXh+hGMnYY5KuMrMrzGy3pJskHZvewMwunXp6g6TnItSLYSGO0BQxVBH94UzEUSyLsv1XX539O/NeH7jGyZi7/17SVyXdr0lAfpkzV3cAABQxSURBVMfdnzGzW83shmKzr5nZM2b2pKSvSbqlab0YFuKoubFPNxFD1dEfbkccRbQo29+7d/bvzHt94MzdU+/DTCsrK76+vp56N9AyM3vc3Vfa+vyxxNHGAeh0u7e0JK2tSQcOpNuvLrQdQ9Jw42jfvslgxVbLy9LLL3e9N2nRFrVgx47JiNhWZtL/+B+Da7SaxBBX4Ad6oGzUi+km1HHbbZP+b9rS0uR1oLFFo18HDkwSr+XlSXK2vNzrRKwpkjGgY1WnE0MWWTPdhDrm9YfSuKe8EUlZtn/gwGQI9u23J68dPjzaoCMZ65OxLwoagDpnr4WMerH8AnVN94cbU5OcYYkoQke/OK2XZKw3CNZBqDOdGDLqxXQTYmHKG1FtzfZnTUMSdCRjrYo5khUSrIycZa/OdGLIqBfLLxALU94j1FXfMa8ego5krDWxR7LKgpWRs16oM50YOuoVcgCKfuuiz2TKe2S66jsW1UPQkYy1Jvawa1mwMszbC3WmExn1ghS3z1yU1DHlPTIx+45FgbWoHoJOcvcsy7XXXuudu/tu9+Vld7PJz7vvrv9ZZu6TNnNzMatX3913uy8tbf6spaV3fiekvth/YwSS1n1ocVSi668gs688urZjyDOJo+Xl2f/Fl5erfU5ZU7KxzZBjZpbBt0XzvtTQviPk85v0UQMIuiYx1GoD1qR0HrghLVQVZS1nnfoWBWtISx37b4xg8A3gHF21Oxl+5dGNJRlb1JdViafQpG4AfWMlg26LFjUEsbL8ss+JVU/GSMZiiB0oZb1g1/W1UWcEg24A5whNkEI7wxg5ep873bEkY/O+yz17qiXcoYP2bSXxucbboNuiRQ1B1S+77gjbCI4MScZiiD2tWLZ9G9OKZdvGGo6OaNAN4BwxBzFjzAz0vX0cSzI277vas6c8nqaFxF9bx205x9ug26JYU4RNR9hyzcQjIRmLoY1pxSb1paqzY4NuAOcIyYlDv6qmMwMZhkRlY0nG3Gf3ZVWPsUKalbaO23KOt0G3RV1MRVbprwaalJGMxTCGacUMD0sH3QDOEfK1hnaGTUe+MhwsrWxMydgsdZqJsr6wraQp53gbdFsUK1GKMcKWYT8UC8lYLH2cVow5dZrAoBvAOWLm4U1nBnIeqQg19mQs9hrEss9s0oTkHG+Db4tiJEoxvsCcg6AhkrEu5DitOIAjjME3gHOUtYux1oyF7EfPQ2j0yZh7vHgq+8whx9uo2qJ5AROzz6m70L/HSMaaiLFwcUPsBYxdT50mMKoGcEroQWpoaDYZ7MxssLQykrFyXSwZCpVrvI2mLVrUr8SajWm60H/rZ+UYMDOQjNVV9TCt6bRirMPT0PpC9zuh0TSAU3IeHeijsSVjdf47xxqMaPo5GTdF42mLFiVDVRKluusfqo6u9aixJBmrK/bIUtenrqWYOo1sNA3glFjtXddijubFNKZkrO5/51iD9iGzWPM+I/OmaDxtUdnVg7u6pk5IQ9Gz2R+Ssbr6fsuimCvBExlNAzilShjEWpjddJuQfUnV2Y4pGav737nsu4nRB8daVZHqAGQ0bVGTjDr0M2IdcfZsfRnJWF1tLJCvO3Q7/ftdTp0mNpoGcEpoOxUrXGJsE7IvqfL+MSVj8/47b/yXrru2MEbfWfYZqa/8X2Y0bVHVacJZX3asJTkDWxdNMlZX14GQYiQr82AeTQM4JbSdinVx2BjbhOxLqrx/TMnYvO+pafIS47sr+4ycE3r3kbVFXSzCjzXv3aMLySZPxiRdL+l5SSckHZrx/rmS7i3ef1TSvrLPzOJsyiFcWyzzhRrTwdvrOKooRjvlHi9J6nNH2nYMeUZxNOu/c4x/8y5GxkKaopQD+WNti+Za9IXGmsHp4uzNDiVNxiTtlPSipCsl7Zb0pKRrtmzzZUl3FI9vknRv2edmEbhtTCs2rTP21GliG8E76DiqKdZAaoxt+rBmrK0Y8sziaPq/87xkbFbyEmNxfZM1Y2X74J7HyBhtUSHWIvymI2whMpkBSp2MfVTS/VPPvy7p61u2uV/SR4vHuyT9SpIt+twsAjfFtOLA5tDLTDWAw42jBsrau67WjIXsS+g2sbUdQ55xHIU2B7G+35CkfVHCl/MAB23RFjGGS8s+J9Y0ZCZro1MnY5+TdOfU87+S9M0t2zwt6bKp5y9KumjR52YTuF1PK5ZtP4Bri02bagCHHUczxPqKYiVJPQmZbdqOIc84jkL7sljHcHX7vD4s/RlzWzRTrEX4MUbYejJIMZhkTNKqpHVJ63v37m3tHyyqNqYVm9TXRp0taqMB7EMc9egrCja0TrQPceQe9u8ea+Cgbp+XSV+50FjbooViDJfmsiCxA6mTsXEP6XadsaeYOm2RRjo1EJpT5zLiFWO6tC1tx5BnHEehYjUJdb/nTGaRFhprWxSs7iL8GCNsbcxAtSB1MrZL0kuSrtA7ix0/uGWbr2jzYsfvlH1urwK362nFGFOnmZhqAEcVRzHar1jblMk9/287hjzjOAo17zs8eLB6/1Wnz+vD8eFY26IgTRfhNx1h60MA+TsxVKc0TsYm9Wu/pJ8VQ7WHi9dulXRD8fg8Sd/V5DTgn0i6suwzex2401JMK/YkcN03B++Y4ijGyH6MbWLMQrhndUmC6DHkGcfRhjojpAcPxlnUH7p/GcwiLTTWtihI6kX4fQgg3xxDVUuUZKyNkn3ghrZSMYcVYtaZiSbBG1JyjaOma15jbBMaJrEuPtuWtmPIM44j9/r/3btehpPBLNJCY22Lzmoyw9PFIvzcA8ibxVCrDViTknXgVm2lYkwrxq4zE2NuABd9RV2MjIUmULmfMzL2ZKxuIhzjQr/uvWlqSo25Lepk7XNPFuE3QTLWtdjDALF63h4adQO4QBdrxqosZwxdf5uiUx57MlZ3irisSWnjGDFno26LYiZKdRf6L/rdniAZ61rsRfkhgT7QQB51A1ii7bMpq+T3OYfW2JOxusdpMQZDhnSMOOq2qIvbEg0pWOYgGetaG/M2ZYE+0CHeUTeAifU0ZLYZezLW5Htc1OzEOkbsi1G3RbGOzGIt9O8pkrGupTjXv4s5/QRG3QBmIOcRr1BjT8bcZ3+PMb7bpseIfTLqtqjKWoS2r7bfYyRjKcRalF8lMJuc7ZKpsTaAbU9BjgnJ2HYxlvjErid3Y22LzgoJhLLseyhrH2oiGctR19OKPT1EHWMD2MXi/OltBtbebUMytl2VMyFjXBR4CDE2xrZorrqL8GONsPUUyViOup5W7Glwj7EB7OKyFe75nwUZC8nYdqED5aHLX/scH6HG2BbN1HQRfowRtp4iGctVjGnFKi1hD1vNMTaAXVzQ1b2d80xyRDK2XWhfF+u2XD1rdmYaY1s0U8xF+HVH2HqKZKyPxtJTlhhjA9jVyFjuV86PhWRsu9Cmo+kSoCE1UaNri5pMRYZk3yO8zAXJWB+lOCMzQ6NrAL27NWMh4TOEA9QxJ2Nll6YIOQGkyclxQ2qiRtUWdZEojfAyFyRjfZXijMzMjKoBnNLF2ZRjyffHmozF6s8WxVCM0deyOnIxqraoi0RphJe5IBkbqhHMEYyqAUygrL0bQAiNNhmrcsZkW5euGNJqi1G1RV0kSkM40quIZKxPqi7IH+CFXqeNqgGcEpIkdXXQ2PcD1LEmY13dO3Iso6+jaou6+FKqBF/fG6ECyVhf1GkZY5yRmbFRNYCFsjDoy0hCLsaajMU6EaSpGKstcjCqtqirRibGwsUeIRnri9gtY18OORcYVQNYaHoGGzYbazIW0oflkAj1JZ5H1xblMvzelwAJ0CSGdgjdefXVsNePHpX27ZN27Jj8PHp09u/ddpu0tLT5taWlyevIVlkYhIbJLKGhg/47cEBaW5OWlyWzyc+1tcnrG/bunf27815vA81Upg4ckF5+WXr77cnP6cCJ6ehRaXVVeuWVSZr1yiuT5xuNU5MGb0BIxroU0jKWBe60kNYY2SkLg7odaJXQwTCU9af798/+vXmvt4FmauQOH5bOnNn82pkzk9elPI4YMkAy1qWQQ8SywN2qq6MbRFMWBnVHEqqGDobv+PFqr7eFZmrEyka+GDqVRDLWrZBDxJAhW+aieq0sDOqOJDDaj62ICSRXNvLF0KmkhsmYmV1oZg+Y2QvFzwvmbPeWmT1RlGNN6uy9skPEssAd4FzUGOOoLAzqjCSMebR/jDEUYswxUQdx1IKQka/pBu+22ybD+SMbbGg6MnZI0g/d/SpJPyyez/Kf7v7HRbmhYZ3DVha4w5yLIo4iGPloPzE0w8hjog7iKLYqI18DHGwIVvc0zMlZnHpe0qXF40slPT9nu99W/ezsTgPu0sCvLTZN0jpxFM9Arp1YSdsx5D2PozHGRB20RRno+WUulPDSFpe4+2vF49clXTJnu/PMbN3MHjGzTzesc/gWzVENc96BOIpkxAuliaE5RhwTdRBHbVu05nnEixx3lW1gZg9Kes+MtzbNi7m7m5nP+Zhldz9lZldKesjMnnL3F2fUtSppVZL29ju5aM9tt02GbaenKnsw73Ddddfp9ddfn/XWu6efEEeYp8sYkoijoaItSmhjGnKj/9qYhpQmRwl7905e22oM/3Z1h9S8wjTllt+5S9LnyrZjSHeBAc07qMLUgBNHmKHtGHLiaBRoizpQNg3Z81sjKeE05TFJNxePb5b0/a0bmNkFZnZu8fgiSR+T9GzDesdtePMOxBGaIoYQA3HUprJpyBFf5qJpMvYNSZ80sxckXVc8l5mtmNmdxTYfkLRuZk9K+pGkb7g7gYtpxBGaIoYQA3HUppA1z8MbbAhSumZsEXd/Q9InZry+LumLxeN/l/SHTerBsBFHaIoYQgzEUct6uua5C1yBHwAAtG/E05BlGo2MAQAABDtwgORrBkbGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACAhkjEAAICESMYAAAASIhkDAABIiGQMAAAgIZIxAACAhEjGAAAAEiIZAwAASIhkDAAAICGSMQAAgIRIxgAAABIiGQMAAEiIZAwAACChRsmYmf2lmT1jZm+b2cqC7a43s+fN7ISZHWpSJ4aHOEJTxBBiII6QStORsacl/YWkH8/bwMx2Srpd0qckXSPp82Z2TcN6MSzEEZoihhADcYQkdjX5ZXd/TpLMbNFmH5Z0wt1fKrb9tqQbJT3bpG4MB3GEpoghxEAcIZUu1oy9T9LPp56fLF4DqiCO0BQxhBiII0RXOjJmZg9Kes+Mtw67+/dj7oyZrUpaLZ6+aWZPx/z8Ci6S9CvqjeoPJJ0z4/UPxq4okzhK9V2mrLvtejuLIWn0cTTk+KUtGn7dqep9f91fLE3G3P26uh9eOCXp8qnnlxWvzaprTdKaJJnZurvPXUDZplR1j63ejboDN+1VHKX+Nx3T39xGDEnjjqOxxm/gprRFmdfdgxjapotpysckXWVmV5jZbkk3STrWQb0YFuIITRFDiIE4QnRNL23xGTM7Kemjkv7VzO4vXn+vmR2XJHf/vaSvSrpf0nOSvuPuzzTbbQwJcYSmiCHEQBwhlaZnU35P0vdmvP4fkvZPPT8u6XjFj19rsm8Npap7bPVK0tpA44j47bDelmNIGuG/aaJ6U9ZNWzScuntXr7l7zB0BAABABdwOCQAAIKFskrGUt6EwswvN7AEze6H4ecGc7d4ysyeKUnvBZtnfYGbnmtm9xfuPmtm+unVVrPcWMzs99Td+MVK93zKzX847rdsm/qnYr5+a2Yca1JUkjsYSQ4F1E0f16x1FHBFDm7brdQwVn0UcbX6/ehy5exZF0gc0uUbHw5JW5myzU9KLkq6UtFvSk5KuiVD3P0o6VDw+JOkf5mz32wh1lf4Nkr4s6Y7i8U2S7u2o3lskfbOF7/bPJH1I0tNz3t8v6d8kmaSPSHq0b3E0hhgijogj2iJiiDhqJ46yGRlz9+fc/fmSzc7ehsLdfydp4zYUTd0o6Ujx+IikT0f4zHlC/obp/blP0ifMFt+fI1K9rXD3H0v69YJNbpT0Lz7xiKR3m9mlNetKFUdjiKHQultBHEVHW7QdMVQdcbRd5TjKJhkL1NZtKC5x99eKx69LumTOdueZ2bqZPWJmdQM85G84u41PTqP+jaQ9NeurUq8kfbYYVr3PzC6f8X4bur69SBv1jSGGQuuWiKO6xhBHxFC79XUZQxJxNEvl77XRpS2qsg5vrVSl7ukn7u5mNu8U02V3P2VmV0p6yMyecvcXY+9rQj+QdI+7v2lmf6vJkcyfJ96nbVLFETEUjDiqWe/0k5HHETFUs97pJyOPIakncSR1nIx5h7dWqlK3mf3CzC5199eKocRfzvmMU8XPl8zsYUl/osmcdRUhf8PGNifNbJekd0l6o2I9let19+k67tRk7UEXqt6mJkkcEUNhdRNHixFHxFDd+kLq7TiGJOJolsrfa9+mKdu6DcUxSTcXj2+WtO2IxswuMLNzi8cXSfqYpGdr1BXyN0zvz+ckPeTFqsAGSuvdMqd9gyZXl+7CMUl/XZyB8hFJv5kaZm9DG3E0hhgKqps4amQMcUQMvaPvMSQRR7NUjyOPfJZB3SLpM5rMq74p6ReS7i9ef6+k41Pb7Zf0M00y+MOR6t4j6YeSXpD0oPT/t3OHOA0EYRiG3zo0R6jiACgO0UvU9BwYToDDIzgEvnaDw3MIDIJxhCC6yUB5HjXZFf9O9hNfspntcly/rh7G+qZa+jyxsVT7E+Z92UN1W+3G+qJ6ql6rY7VdaZ8/zb2rXsYen6urleY+Vm/V+3jH++pQHcb9TXU/nmvpm5NHvzlH/yVDciRHMiRDcrR+jvyBHwBgor/2mRIA4KwoYwAAEyljAAATKWMAABMpYwAAEyljAAATKWMAABMpYwAAE30Awp+alJ2J0isAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fcec20be470>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"clf = PCA(n_components=2)\n",
"npos = [0, 0, 0, 0]\n",
"npos = [smacof_mds(Cs[s], 2) for s in range(S)]\n",
"\n",
"npost01 = [0, 0]\n",
"npost01 = [smacof_mds(Ct01[s], 2) for s in range(2)]\n",
"npost01 = [clf.fit_transform(npost01[s]) for s in range(2)]\n",
"\n",
"npost02 = [0, 0]\n",
"npost02 = [smacof_mds(Ct02[s], 2) for s in range(2)]\n",
"npost02 = [clf.fit_transform(npost02[s]) for s in range(2)]\n",
"\n",
"npost13 = [0, 0]\n",
"npost13 = [smacof_mds(Ct13[s], 2) for s in range(2)]\n",
"npost13 = [clf.fit_transform(npost13[s]) for s in range(2)]\n",
"\n",
"npost23 = [0, 0]\n",
"npost23 = [smacof_mds(Ct23[s], 2) for s in range(2)]\n",
"npost23 = [clf.fit_transform(npost23[s]) for s in range(2)]\n",
"\n",
"\n",
"fig = pl.figure(figsize=(10, 10))\n",
"\n",
"ax1 = pl.subplot2grid((4, 4), (0, 0))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax1.scatter(npos[0][:, 0], npos[0][:, 1], color='r')\n",
"\n",
"ax2 = pl.subplot2grid((4, 4), (0, 1))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax2.scatter(npost01[1][:, 0], npost01[1][:, 1], color='b')\n",
"\n",
"ax3 = pl.subplot2grid((4, 4), (0, 2))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax3.scatter(npost01[0][:, 0], npost01[0][:, 1], color='b')\n",
"\n",
"ax4 = pl.subplot2grid((4, 4), (0, 3))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax4.scatter(npos[1][:, 0], npos[1][:, 1], color='r')\n",
"\n",
"ax5 = pl.subplot2grid((4, 4), (1, 0))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax5.scatter(npost02[1][:, 0], npost02[1][:, 1], color='b')\n",
"\n",
"ax6 = pl.subplot2grid((4, 4), (1, 3))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax6.scatter(npost13[1][:, 0], npost13[1][:, 1], color='b')\n",
"\n",
"ax7 = pl.subplot2grid((4, 4), (2, 0))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax7.scatter(npost02[0][:, 0], npost02[0][:, 1], color='b')\n",
"\n",
"ax8 = pl.subplot2grid((4, 4), (2, 3))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax8.scatter(npost13[0][:, 0], npost13[0][:, 1], color='b')\n",
"\n",
"ax9 = pl.subplot2grid((4, 4), (3, 0))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax9.scatter(npos[2][:, 0], npos[2][:, 1], color='r')\n",
"\n",
"ax10 = pl.subplot2grid((4, 4), (3, 1))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax10.scatter(npost23[1][:, 0], npost23[1][:, 1], color='b')\n",
"\n",
"ax11 = pl.subplot2grid((4, 4), (3, 2))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax11.scatter(npost23[0][:, 0], npost23[0][:, 1], color='b')\n",
"\n",
"ax12 = pl.subplot2grid((4, 4), (3, 3))\n",
"pl.xlim((-1, 1))\n",
"pl.ylim((-1, 1))\n",
"ax12.scatter(npos[3][:, 0], npos[3][:, 1], color='r')"
]
}
],
"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
}
|
