summaryrefslogtreecommitdiff
path: root/notebooks/plot_free_support_barycenter.ipynb
blob: b8df5898f3ae5d38a1dc256f215c879931b6283f (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
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 2D free support Wasserstein barycenters of distributions\n",
    "\n",
    "\n",
    "Illustration of 2D Wasserstein barycenters if discributions that are weighted\n",
    "sum of diracs.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Vivien Seguy <vivien.seguy@iip.ist.i.kyoto-u.ac.jp>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data\n",
    " -------------\n",
    "%% parameters and data generation\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 3\n",
    "d = 2\n",
    "measures_locations = []\n",
    "measures_weights = []\n",
    "\n",
    "for i in range(N):\n",
    "\n",
    "    n_i = np.random.randint(low=1, high=20)  # nb samples\n",
    "\n",
    "    mu_i = np.random.normal(0., 4., (d,))  # Gaussian mean\n",
    "\n",
    "    A_i = np.random.rand(d, d)\n",
    "    cov_i = np.dot(A_i, A_i.transpose())  # Gaussian covariance matrix\n",
    "\n",
    "    x_i = ot.datasets.make_2D_samples_gauss(n_i, mu_i, cov_i)  # Dirac locations\n",
    "    b_i = np.random.uniform(0., 1., (n_i,))\n",
    "    b_i = b_i / np.sum(b_i)  # Dirac weights\n",
    "\n",
    "    measures_locations.append(x_i)\n",
    "    measures_weights.append(b_i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute free support barycenter\n",
    "-------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "k = 10  # number of Diracs of the barycenter\n",
    "X_init = np.random.normal(0., 1., (k, d))  # initial Dirac locations\n",
    "b = np.ones((k,)) / k  # weights of the barycenter (it will not be optimized, only the locations are optimized)\n",
    "\n",
    "X = ot.lp.free_support_barycenter(measures_locations, measures_weights, X_init, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot data\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYU2X2wPHvSTKNAQap0hEUBIYmMAooRVRcRMG2/hSkWJC1N1RUXF1Xd91VsS82FGzAioiriIgFEVAERKoC4lBGytCZnvL+/kgyTkmbSTKTDOfzPPNMkvvm3pObzJk3733vuWKMQSmlVM1hqe4AlFJKRZYmdqWUqmE0sSulVA2jiV0ppWoYTexKKVXDaGJXSqkaRhO7UiEQkbEi8m202pd57kgRWViB9m1ExIiIrTLbUzWPJvYoEpFMEckXkWMiclhElonIBBEJab/rH2x8iPT7ZIx5xxhzXiTWFS/C+UeoytPEHn0XGmPqAK2BfwL3Aq9Xb0ixS/+JBRbt/ROv+z9e444WTexVxBhzxBjzEXAFMEZE0gFE5AIR+VFEjorIThF5uMTTvvH8PiwiOSLSR0TaiciXInJARPaLyDsiUs/fdj09yRtFZIvnm8OjnnUs82xztogklmg/TETWlPiG0bXEsvtE5FfPejaKyMUllp0sIotF5Ignrlmex8v1ZkXkaxG5znN7rIgsFZEpInIAeNjz+DUisklEDonIZyLS2vO4eNru88S/zrsvfbz2cZ51HBORbSJyQ4llA0Vkl4jc5VnXbhEZV2J5AxH5yLONFUA7f/vY1/tUYj1Pel7DbyLypxKPp4nI657tZonI30XEWmKffFuirRGRm0RkC7AlQBzXiMjvnnXeXeL5GSKy3POe7haRF8q856XWLyIvishTZfblRyJyh+d2SxH5QESyPZ/DF0q08/m+ldjOBM9n8bBnOyIiHYGpQB/P/jvsaZ/k2X87RGSviEwVkZQy79+9IrIHeCPAfjn+GGP0J0o/QCZwjo/HdwB/8dweCHTB/U+2K7AXGOFZ1gYwgK3Ec08GzgWSgEa4k8ozAWIwwDygLtAZKAS+ANoCacBGYIynbQ9gH3A6YAXGeF5Dkmf55UAzT6xXALlAU8+y94AHPMuSgTMDvIavges8t8cCDuAWwAakAMOBrUBHz2MPAss87YcAq4B6gHjaNPXz2i/AnZAFGADkAaeV2O8O4G9AAjDUs/wEz/KZwGwgFUgHsoBv/WzH12scC9iB6z378i/A74B4ls8FXvasvzGwArihxHO/LfMefg7UB1ICbP89z/q6ANl4PntAT+AMz75sA2wCbve3fiDDE6vFs7yhZ9808byWn4Apnm2VfK/9vm8ltvOx571r5YnxfF+v2fPYFOAjT1x1gP8B/yjz/j2B+2+h3H45nn+qPYCa/IP/xP4d8ICf5zwDTPHcLpcwfLQfAfwYYLkB+pW4vwq4t8T9p/D8YwD+Azxa5vm/AAP8rHsNMNxzewbwCtCiTJtyr4HyiX1Hmed8Clxb4r7Fk1haA2cDmz2JylLB9+ND4DbP7YFAfpm49nnWa8WdlE8tsezxsoknyGscC2wtcb+Wp82JuBNkYclkBFwJfFXiuWUT+9kBXpd3+yXj/Rfwup/2twNzA60fd/I/13P7ZmC+53Yf3Am53Gcy0PtWYjtnllg+G7jPz2sW3B2HdiUe6wP8VuL9KwKSK/p3eTz86FBM9WgOHAQQkdNF5CvP19ojwATcPSSfRKSJiMz0fH0/CrwdqL3H3hK3833cr+253Rq4y/M1+bDnK3FL3L10RGS0/DFMcxh3T9a77Xtw/zGuEJENInJN0L3wh51l7rcGni2xnYOedTc3xnwJvAC8COwTkVdEpK6vlYrIn0TkOxE56FnPUErvqwPGGEeJ+3mefdEId4+zZFzbK/B6vPZ4bxhj8jw3a3teXwKwu8RrfBl3z92fsvsoWJvt/PG+tReRj0Vkj+cz8zjlPzNl1z8dGOW5PQp4y3O7JbC9zH7z8vu+lWizp8Rt7/72pRHuf4arSqxvgedxr2xjTIGf5x/XNLFXMRHpjfuD7h1DfRf3182Wxpg03GON4lnmq/Tm457Huxhj6uL+oxMf7SpjJ/CYMaZeiZ9axpj3PGOlr+LuvTUwxtQD1nu3bYzZY4y53hjTDLgBeElETsbd6wL3H6nXiWW2W/Z17sQ9LFEyjhRjzDLPtp4zxvQEOgHtgYllX4iIJAFzgCeBJp545xPavsrG/TW/ZYnHWgVoX9ESqTtx99gblnh9dY0xncPcRtl4f/fc/g/wM3CK5zNzP+X3Q9n1vw0MF5FuuIdWPiwReyvxfbAy4PsWRNnt78fd6ehcYl1pxpjaAZ6jPDSxVxERqSsiw3CP3b5tjFnnWVQHOGiMKRCRDOCqEk/LBly4x8Mp0T4HOCIizfGR1MLwKjDB8y1CRCRV3Ad36+AeTzWemPAcaCw+aCkil4tIC8/dQ562LmNMNu7x6VEiYvX05AMdiAT3P7dJItLZs+40Ebncc7u3J74E3P80CnDvo7IScY+9ZgMOz4HLkKYQGmOcwAfAwyJSS0Q64T7e4I+v9ynQ+ncDC4GnPJ8Li7gPaA8I5fkBTPbE2xkYB8zyPF4HOArkiMipuMf7g8W4C/gBd099jjEm37NoBbAb+Kfn85EsIv08y/y+byHYC7TwHtQ1xrhwfx6niEhjz/qai8iQENd3XNPEHn3/E5FjuHszDwBP4/6j87oR+JunzUO4xx2B4q/vjwFLPV9HzwAeAU4DjgCf4E5AEWGMWYn7YN8LuJPzVtxjnxhjNuIej1+O+4+wC7C0xNN7A9+LSA7ubyC3GWO2eZZdj/sf0AHcB3AD9uCMMXNxHxSb6Rk6WA94Z5TUxf0Hfwj3cMMB4N8+1nEMuBX3/jyE+x/mRyHtCLebcQ8T7AHeJMCsCz/vUzCjcf/z2eiJ732gaQXi82Ux7vfsC+BJY4z3JKe7cb/+Y7j33SzfTy9nOu732TsM4/2ndyHug/g7gF24D6QHe9+C+RLYAOwRkf2ex+71vJ7vPOtbBHQIcX3HNe8ReqWUKkVE+uMekmltNFHEFe2xK6XK8Qx13Qa8pkk9/mhiV0qV4jlh6DDuoaFnqjkcVQk6FKOUUjVMRHrsInK+iPwiIltF5L5IrFMppVTlhN1jF3d9i824T3P3TpG60jOLwqeGDRuaNm3ahLVdpZQ63qxatWq/MaZRsHaRqIiWgfvU6W0AIjITd80Iv4m9TZs2rFy5MgKbVkqp44eIhHQGdCSGYppT+nTkXZQ+hdgb0HgRWSkiK7OzsyOwWaWUUr5U2awYY8wrxphexphejRoF/SahlFKqkiKR2LMoXaOihecxpZRS1SASif0H4BQROclT5+H/qNip20oppSIo7IOnxhiHiNwMfIa7jvU0Y8yGsCNTSql4ZwxkLoFfvwJHITTuCOmXQGJqVDcbqYvvzsddElUppRTA7p9g5kjIPwhFnurVCanw6UQY9CD0uQkkUhW3S9MLwCqlVKTt2wRv/OmPhO5l99z/6jFwFsJZd0Vl81orRimlIu2Tu8on9ZLsefD1E5B7ICqb18SulFKRdHgHZK0K3k4EVk+PSgia2JVSKpL2bgBrYvB2jgLYuSIqIWhiV0qpSBILIV+OVaKTgjWxK6VUJDXrAY6i4O0SakG7QVEJQRO7UkpFUu3G0O5sEGvgdsZAt/+LSgia2JVSKtIueBKS0/wnd1sKXPgMJNWJyuY1sSulVKSltYAbFkPLDLAlu4ddbMmQWBtqN4FLXo5abx30BCWllIqOeq3gmgVw4Ff47Rtw2qFRe2jTHyzR7VNrYldKqWhq0M79U4V0KEYppWoY7bErpVSkZa2GZc/B5s/cNWFSG0PGeOg5FmrVj/rmtceulFKRtPQ5eHMobJznrgnjcsKx3bD4X/BCL9i/JeohaGJXSqlI2fwZfP042PPBuEovc+RD3kF48wJ3bfYo0sSulFKR8uVj7qTul3FXfdwY3YvMaWJXSqlIOPo77P85eLuiHFj5elRD0cSulFKRkLs/tKqOALnZUQ0lrMQuIpeLyAYRcYlIr0gFpZRScadWffdJSCG1bRjVUMLtsa8HLgG+iUAsSikVv9JaQIOTg7dLrA29xkU1lLASuzFmkzHml0gFo5RScW3Q/ZCQEqCBQEIydBoR1TCqbIxdRMaLyEoRWZmdHd3xJaWUqhanXuC+QHVCLUBKL7MmQUo9GPOxO7lHUdDELiKLRGS9j5/hFdmQMeYVY0wvY0yvRo0aVT5ipZSKZf0nwtVzof35YE1wP1arAZx5O9z0AzTuGPUQgpYUMMacE/UolFKqJml1Blw1033b5Yp6NceydLqjUkpFUxUndQh/uuPFIrIL6AN8IiKfRSYspZRSlRVWdUdjzFxgboRiUUopFQE6FHOcMsZUdwhKqSjReuzHke+3HWDq4l/5dut+7E5Dw9qJjOnbhlGnt+aE1BBPhVZKxTxN7MeJpxb+wmtLfqPA7sTbV9+fU8QLX27ljaWZvD+hD20b1a7WGJWKO3vWw/IXYMvn4LJDvdbQ9xboMBQw7rNMRYKuJtKkOr6S9+rVy6xcubLKt1udnC7DN5uz+SHzIC5jSG+exnmdTiTRFv3RsI/X/s7E/64l3+70uVwEGtdJ4tt7zybBqqNzSoVkydOw+Al3fRhT4m9LLO5a7GIDWxKcdjX0vRXSmoe9SRFZZYwJWpdLe+xVYNmv+7n1vR/JL3KSW+T+AKQmWZn0wToev7gLF3ZrFtXtP71ws9+kDmAM5BQ6WLhhLxd0bRrVWJSqEdZ/AN/8CxwF5Zd5L7BhHGB3wA+vw5p3Yewn0LRrlYSn3bMo+37bAa598wf25xQVJ3WA3EInxwocTHz/J/73U1bUtp+5P5ffjwQq/P9HPO+t2BG1OJSqMYyBLx4JckGNElx2KDwKM4aH/pwwaWKPImMM98xZS77d5bdNgd3FpA/WU+Tw3yYcB/OKQh5e2Z8T3ct1KVUj7N0AOfsq/jxnIWyomtnhmtij6Medh8k+FjxZGgyfbdgTlRgapCZid4b2T6NRnaSoxKBUjZKzByyVGMUuyoXV0yMfjw+a2KNoQ9YRXK7gB6dzC538tOtwVGJo3SCVFifUCtouNcnKyNNbRSUGpWqU5HrlL1QdqvxDkY3FD03sUSQi5Sp3+mOJ4pSou89rT0qCNcC2oW5yAud0bBK1GJSqMZr1AFsly+7WqZrJCZrYo+i0VieElNdTk6z0blM/anGcn96UGwe2IyXBWm5KbbLNQv3URGaOPwObTnVUKjiL1V2CNyH4N+FSEmtD7+ujE1MZ+pccRZ2a1aV1g9Sg7RKtVs4+tXFUY7ll8Cm8c/3pDOnUhOQEC1YRTqybzB3ntueLuwaGFKdSyuOMm6DDnyqQ3MVdk739+VENq3hreoJSdG38/SiXTV1GXpHveeTJCRZeHd2Ls04pffGRAruTvCIndZJtetKQUrHIGPjpXfjwJiBYHhW4+QdoeEpYm9QTlGJEp2Z1eX9CX26f9SM7D+ZjjMEAVotQPzWRf13Wlb7t/rhi+bdb9vPCV1v4IfMQVs8Y/QVdmnLToHac3LhO9b0QpVRpIoC4h2ZcjsBtE1Jgz9qwE3uoNLFXgU7N6rLwjgGszzrC6h2HcLkMnZun0av1Ce4DrB7Pf7GFl77+tfgsUW9Vl4/W/M6C9Xt4dXQvzjyloc9tKKWq0L5NsHY2bJwXPKmD+8SkQ9ujH5eHJvYqlN48jfTmaT6XfbM5u1RSL8lpDPl2J+PfWsniiYN0vrlS1SXvIMwaBVmrwVlUukZMQAbE/8y0SNPB2xjx3BdbAtZzAXchMT3tX6lqYs+HaefDrh/AkV+BpO7R+NToxOWDJvYYkFPoYM3O4CcoFTpcvL9qVxVEpJQqZ827cGSHu6deYQL1qu4EwHCvefpvEflZRNaKyFwRqRepwI4nOQUObNbQTlDKLQxhPE8pFXnLnqt8ES+LDeq3i2w8gTYX5vM/B9KNMV2BzcCk8EM6/tSrlYAzhNIDAA1q65WOlKpyxsDhMA5+tu4Ltqr72w0rsRtjFhpjvF3I74AW4Yd0/ElOsDL41MZBL7RSK9HKmL5tqiQmpVRZlS37IfCnJyIaSTCRHGO/BvjU30IRGS8iK0VkZXZ2dgQ3WzPcMvgUkgJcTUmAlEQrI7qHfxUWpVQFiUDT7pV5IpzUHxp3jHhIgQRN7CKySETW+/gZXqLNA4ADeMffeowxrxhjehljejVq1Mhfs+NW52ZpPH/laaQkWEkoM96ekmClfmois2/oQ2qSzlBVqlqceTskVKD0hjUJ6reFy9+MWkj+BM0SxphzAi0XkbHAMGCwqY76BDXIuZ2a8MVdA5ixfDsfrskir9BBw9pJjOnbhktOa06d5ITqDlGp49epF0LbWfDrF74viedlTYLEVDj9BuhzEyRV/RnjYdWKEZHzgaeBAcaYkMdXjqdaMUqpGuTAVnipj/8pj9Yk6DQcLp7qLjUQYaHWigl3jP0FoA7wuYisEZGpYa5PKaVi15ePBS4h4CyETR/BgV+rLiYfwhqwNcacHKlAlFIqpuUfgl/mB796kssB3/0HLpxSNXH5oGeeKqVUKA5sA2sIc9FdDvh9dfTjCUATu1JKBeNyuWvE2HNDa5+bDYXHohtTAJrYlVIqkJxsmNoXvvwbuEIs/HVsL0w9C3L3Rzc2PzSxK6WUP047vDkU9m+FohB76wDGAUd2wXv/F73YAtDErpRS/vwyH45mgcte8ee67LB3PexZF/m4gtDErpRS/ix/qWI99bIcRbDu/cjFEyJN7Eop5c/hMC9sY5yQsy8ysVSAJnallPInITm851sSIK3qC/fFTUWpAruTRZv2svNgPikJFgad2pjWDSpQkEcppSqq8yWw7Hn3GaWVYbFCtysjG1MIYj6xG2N4dck2nl20BcSd4G0WC//49Gd6tKrHc1f2oHGdMP+rKqWUL72vheUvVu651iRocxY0qLorJ3nF/FDMUws3M+XzLeQWOcktdOJ0ua/9WehwsTLzEBc+/y0HcytzDUKllAqibjMY8R8qnCoTUqFJZ7j8jaiEFUxMJ/btB3J5dck28u2+TwpwuAwHc4p4+vNfqjgypdRxI/1iaNIp9PYntIOLnoNrF1ZLyV6I8cT+5rJMXEHKCttdhjmrssgvCvGMMKVUUNOnT2fWrFnVHUbsaDcotDox1iS4dgF0uQys1Xf9hJhO7N9tO4DdGbxevNUC2/bnVEFEStV8OTk53HLLLUyYMIGCggAXlDie9L4OJFi6FGg7EGo3roKAAovpxB76NUAqe5FZpVRZzz//PE6nE7vdzssvv1zd4cSGE9pAz3GQUMt/m8RUGPJYlYUUSEwn9t5t6mOzBE/aDpeLkxrq1EelwpWTk8M//vEP8vLyyM3N5eGHH9Zeu9eQxyFjvHu4xZbyx+OJtaFucxj3KTQ8pfriKyGmE/u4fm2wBknsNoswontzaiXG/MxNpWKSw/HHFYG8vXUv7bWXYLHAuY/AXT/D4IfgtNFw+l/girfhjg3QtGt1R1gsrGueVlZFrnn66Mcbeff7HT5nxlgtwgm1Eph/21k6l12pSnj66ad5//33Wbp0Kbm5uTRr1oxjx0rXEa9Xrx67d+8mOVn/xqpblVzzVEQeFZG1nuudLhSRZuGsz5cHL+jIhAFtSUmwUivRfXFYm0VITrCQ3qwu/7vlTE3qSlVCbm4ujzzyCKtXr+brr78u11v30l57/Amrxy4idY0xRz23bwU6GWMmBHteRXrsXrmFDuav283Og3mkJNoY3LEx7ZtUzxxRpWqCJ554gr/97W/k5eXRrVs3tm3bVq637qW99tgQao893ItZHy1xNxWI2rhOapKNy3u1jNbqlTqu5Obm8vjjj5OXlwfAxo0bsVj8f4H39tpvu+22qgpRhSHsg6ci8piI7ARGAg+FH5JSKtpeeOGFUgdN7XY7hYX+C13pDJn4EnQoRkQWASf6WPSAMWZeiXaTgGRjzF/9rGc8MB6gVatWPbdv317poJVSlec9SHr06NFyy1JTU7HZfH+RP3bsGG+++SZXX311tENUfkRsKMYYc06I23wHmA/4TOzGmFeAV8A9xh7iOpVSEVa2t15S8+bNee211xDxPc24e/fu0QxNRUhYY+wicooxZovn7nDg5/BDUkpF0v79+3niiSf417/+RV5eXqmx9bKysrJwOBwMGjSoiqNUkRTuWT3/FJEOgAvYDgSdEaOUqlqPP/44U6ZMYcCAAWzYsMFvbx3cwzR33303K1eu9NtrP+447bDpI1j6HBzYAmKF1n2h763u3zG4n2L+BCWlVOVlZ2fTunVr8vPzOeWUU9izZ4/fKY1eycnJzJ8/X3vtAIXHYPqFkL0Z7CUvai2QkAJd/wzDnqmy5F4l0x2VUrHtH//4B97O265du7BYLNSuXTvo85YsWaKJHWD2WNi70cel8QzY82DtbDjhJDjz9uqIzi9N7ErVUNnZ2UydOrV4imJ+fj4dOnRg06ZNOswSiv1bYPu3ga93as+Db5+GPjdVa/31smK6CJhSqvJK9ta9srKy+OSTT6opojiz5l1w+T8eUczlgm2Lox9PBWhiV6oG+OWXX1i2bFnx/bK9da+cnBzuvvvucglf+XBkV2iJ3bggd1/046kATexK1QCjRo1ixIgRFBW5L+zuq7fupb32ENVuREgp0mKB5HpRD6ciNLErFee+/vprNm3aRF5eHm+++abf3rqX9tpD1OXPkJAUvJ3L6b4kXgzRxK5UnJs4cSK5ubnk5uby4IMP8ve//z1o0tZeewiadYdGp4IlwBwTWwr0vhYSA1wyrxroPHal4tj//vc/rrjiCvLz8wF3rZeCggKSk5P91nwBKCgooFu3bnz//fdVFWp8OrYXXj8HcvaBo8w3oIRa7hOUrpxZZTNidB67UseBsWPHFid1cJ85mpaWxvTp00lICJxsmjWL+HVxap46TWDCUvjhdfjuJcg7ABho2ME9d73L5WCxVneU5WhiVypOzZw5k4MHD5Z73OFwsHfvXsaPH18NUdVAyXXhrDvcP/Z8d0kBW2J1RxWQDsUoFacaNGjgM7EDNGrUiF27dpGYGNsJSFVMlVzzVClVPfz11r28M2TU8alG9tgdThdZh/NxugzN6qWQnBB7Y2BKhSNQb91Le+01z3F58DSvyMHUxb8yY9l2ipwuBHAZuOS05tw2+BQa19UL8ar4N2fOnIBJ3Xvt0oMHD/L2229zzTXXVFVoKkbUmMSeU+jg0peWkXkgl0KHq9SyWT/s5NP1e5h3Uz9a1o+t+aZKeb399tts27aNhx4KfOngd999F4vFgsvlKrcsJSWF++67r3iqY48ePaISq4ptNSax3//BOn7bn0ORs/zQksNlOJxXxLg3f+DzO/prZTsVcwoKCrj11lvJy8tj/PjxnHiir8sMw/bt2/n00099JnVw99YbNGjATTfdFM1wVYyrEQdPD+UWsWDDHp9J3ctlIOtQPj/uPFyFkSkVmldeeYWioiKMMTz66KN+202ePBm73e53eW5uLg899BCFhQFKzaoar0Yk9m+2ZJNgCd4LL7A7+XTd7iqISKnQFRQU8PDDD5Obm0tRURHTpk1jz5495dpt376dd955J+Cl7QCOHj3Ka6+9Fq1wVRyoEUMx+UVOXCFM7jHAsYIQynAqVYW8vXUvl8vFo48+yosvvliqnTGGq666CqfTGXSdelbp8S0i0x1F5C7gSaCRMWZ/sPaRnu74zeZs/vLOKnILA3/gk2wWbj+nPX8Z2C5i21YqHAUFBTRr1oxDhw6Vejw5OZnffvvN71i7Oj5V2QlKItISOA/YEe66KqtvuwYkWIO/FANcelrz6AekjhvGGD799NNKl8At21v38vbalaqMSIyxTwHuwZ03q4XNamHieR1ICXAiUnKChYt7NNe57CqiPvvsM4YOHcrChQsr/NySY+tlBRprV1XIGHAUun/HkbASu4gMB7KMMT+F0Ha8iKwUkZXZ2dnhbNankWe0ZsKAtiTZLNhKHEi1CKQkWBl0amP+PiI94ttVxy9jDHfddRdApS5c4a+37qW99mq0axXMvAoebQiPnQiPN4OP74SDv1V3ZCEJOsYuIosAXwN9DwD3A+cZY46ISCbQqzrG2Ev6bX8ubyz9jW+37MdpDOnN07j+rLZ0a5Gm89dVRC1YsIDLLruM3NxcUlNTmTNnDkOGDAnpuf7G1svSsfZqsOJV+Hwy2AsoNRBhSQBrIoz8L7TpVy2hhTrGXumDpyLSBfgCyPM81AL4HcgwxgT8/qjVHVW8M8aQnp7Oxo0bix9LT09n7dq1IXUg5syZw+WXX05qamrAdvn5+Tz++OPcc889YcesQrB9Obx9sbs8rz+JteHWNZ5rolatqNeKMcasAxqX2GAmIfbYlYp3n332Gdu3by/12G+//cbChQtD6rWPGDGCzZs3h7Stli1bVipGVQmLnwic1AFcDlg9HfrfXTUxVULEqjvGylCMUtHmq7fuVZFeu4ox9nz4Rwt34g6mXmu4fW30YyqjyuuxG2PaaG9dHQ989da9vL12FYcKjwW+cHVJBUeiG0uYakRJAaWqincmjK8piuCu1VKZGTIqBiTVBRP8rF4AajWIbixh0sSuVAUE6q17aa89TiUkwynnAUGG0RJqQe9rqySkytLErlQF3H333X576165ublMnDixiiJSEdV/ojvBB2JNhO4jqyaeSqoRRcCUqipXXHEFu3e7K4QeOnSI//73v6WKclmtVoYPH64XuIhXzXrA8JfgwxvdB1FdJUok25LBlgRjPoKUetUXYwg0sStVAZMnTy6+ffHFF5db7nK5yM7O5sEHH6zKsFQkpV8CTbvB8hdhwwdgz4OU+tDrWug1DlIbVneEQdXIi1krFW0///wzPXr0oKCgoNyyWrVq8eWXX3L66adXQ2SqJqvy6Y5KHU8mTZrk90pG+fn5OsauqpUmdqUq6Oeff2bBggV+L3hhjGHVqlV8//33VRyZUm6a2JWqoEmTJgWsygiQl5envXZVbfTgqVIVJCJ07NgxaLukpKQqiEap8jSxK1VBH3zwQXWHoFRAOhSjlFIHf7g5AAAeoElEQVQ1jCZ2pcLgcrmqOwSlytHErlQluVwuevbsyezZs6s7FKVK0cSuVCXNmzePDRs2cOedd/qd+qhUddDErlQluFwuJk6ciN1u58iRI8ycObO6Q1KqmCZ2pSph3rx57N27F4CcnBzuvfde7bWrmKGJXakK8vbWc3Jyih/TXruKJWEldhF5WESyRGSN52dopAJTKlaV7K17aa9dxZJI9NinGGO6e37mR2B9SsUsX711L+21q1ihQzFK+fHuu++W65n76q17aa9dxYpIJPabRWStiEwTkRP8NRKR8SKyUkRWZmdnR2CzSkXPtm3buPrqq5k0aVLxY4F6617aa1exIGhiF5FFIrLex89w4D9AO6A7sBt4yt96jDGvGGN6GWN6NWrUKGIvQKlomDx5MhaLhffee4+dO3cC7t6697J4/mivXcWCoEXAjDHnhLIiEXkV+DjsiJSqZtu2beODDz7A4XAgIvz1r39l2rRpLFy4EKfTGbRq4+HDh/n1119p3759FUWsVGlhXRpPRJoaY3Z7bt8BnG6M+b9gz9NL46lYNnLkSGbPno3D4QAgOTmZzZs307Jly2qOTB3vqurSeP8SkXUishYYBNwR5vqUqlYle+teTqeTv/71r9UYlVIVE1ZiN8ZcbYzpYozpaoy5yNt7VypeTZ48uVRSB7Db7aXG2pWKdTrdUSkPX711L+21q3iiiV0pD1+9dS/ttat4ooldKQL31r20167ihSZ2pYDHHnsMh8NBUlKS3x8RYcaMGX7PPFUqVujFrJUCxowZQ48ePYK2s9ls1KlTpwoiUqryNLErBfTv35/+/ftXdxhKRUTMJHa73c6uXbsoKCio7lCUIjk5mRYtWpCQkFDdoShVYTGT2Hft2kWdOnVo06YNIlLd4ajjmDGGAwcOsGvXLk466aTqDkepCouZg6cFBQU0aNBAk7qqdiJCgwYN9Nujilsxk9gBTeoqZuhnUcWzmBmKqYjtB3J5dck2Pvzxd3ILHaQm2RjRoxnXn9WW1g1Sqzs8pZSqVjHVYw/FV7/s4/xnljBzxU5yCh0YIKfQwcwVOzn/mSV89cu+Sq+7b9++kQvUIzMzk3fffTfi61VKKX/iKrFvP5DLjW+vJt/uxOEqXW7Y4TLk253c+PZqth/IrdT6ly1bFokwS4m3xB7ozEulVHyIq8T+6pJt2J2ugG3sThevLfmtUuuvXbs2AF9//TUDBw7ksssu49RTT2XkyJF469a3adOGe+65hy5dupCRkcHWrVsBGDt2LO+//365dd13330sWbKE7t27M2XKlFLb+/rrrxkwYADDhw+nbdu23HfffbzzzjtkZGTQpUsXfv31VwCys7O59NJL6d27N71792bp0qUArFixgj59+tCjRw/69u3LL7/8AsCGDRvIyMige/fudO3alS1btpCZmUl6enrxtp988kkefvhhAAYOHMjtt99Or169ePbZZ/1uTykVH+JqjP3DH38v11Mvy+EyzP0xi0dHpAdsF8yPP/7Ihg0baNasGf369WPp0qWceeaZAKSlpbFu3TpmzJjB7bffzscf+79w1D//+U+efPJJv21++uknNm3aRP369Wnbti3XXXcdK1as4Nlnn+X555/nmWee4bbbbuOOO+7gzDPPZMeOHQwZMoRNmzZx6qmnsmTJEmw2G4sWLeL+++9nzpw5TJ06ldtuu42RI0dSVFSE0+kMehp8UVER3oufXHXVVT63p5SKD3GV2HMLQxsmyC0KfzghIyODFi1aANC9e3cyMzOLE/uVV15Z/PuOO8K7tkjv3r1p2rQpAO3ateO8884DoEuXLnz11VcALFq0iI0bNxY/5+jRo+Tk5HDkyBHGjBnDli1bEBHsdjsAffr04bHHHmPXrl1ccsklnHLKKUHjuOKKK4pv+9ue91uIUiq2xVViT02ykRNCck9NDP9llbyupdVqLTX2XHIqnPe2zWbD5XIPE7lcLoqKiiq8HYvFUnzfYrEUb9PlcvHdd9+RnJxc6rk333wzgwYNYu7cuWRmZjJw4EDA3eM+/fTT+eSTTxg6dCgvv/wy7du3L44PKDdHOzX1j9lE/ranlIoPcTXGPqJHM2yWwPOLbRbh4h7NoxrHrFmzin/36dMHcI+9r1q1CoCPPvqouPdcp04djh07Ftb2zjvvPJ5//vni+2vWrAHgyJEjNG/ufq1vvvlm8fJt27bRtm1bbr31VoYPH87atWtp0qQJ+/bt48CBAxQWFgYcPvK3PaVUfAg7sYvILSLys4hsEJF/RSIof64/qy0J1sAhJ1gtXHdWdE8DP3ToEF27duXZZ58tPiB6/fXXs3jxYrp168by5cuLe8Bdu3bFarXSrVu3cgdPQ/Xcc8+xcuVKunbtSqdOnZg6dSoA99xzD5MmTaJHjx6lvlHMnj2b9PR0unfvzvr16xk9ejQJCQk89NBDZGRkcO6553LqqadWeHtKqfgg3tkelXqyyCDgAeACY0yhiDQ2xgSdSN6rVy/jPVDntWnTJjp27Bh0m1/9so8b316N3ekqdSDVZhESrBZeGnUagzo0rvBrCVWbNm1YuXIlDRs2jNo2VGwI9TOpVFURkVXGmF7B2oXbY/8L8E9jTCFAKEk9XIM6NGbB7WdxZUYraifZEIHaSTauzGjFgtvPimpSV0qpeBDuUcb2wFki8hhQANxtjPnBV0MRGQ+MB2jVqlVYG23dIJVHR6SHPaWxMjIzM6t8m0opVRFBE7uILAJO9LHoAc/z6wNnAL2B2SLS1vgY3zHGvAK8Au6hmHCCVkop5V/QxG6MOcffMhH5C/CBJ5GvEBEX0BDIjlyISimlKiLcMfYPgUEAItIeSAT2hxuUUkqpygt3jH0aME1E1gNFwBhfwzARd3AbLHsB1s6GohxIrA1d/wx9b4b6baO+eaWUimVh9diNMUXGmFHGmHRjzGnGmC8jFZhfWz6H//SD1TOg6Bhg3L9Xz3A/vuXzSq9ay/YqpWqCuDrzlIPbYPZosOeBy156mcvufnz2aHe7StCyvVq2V6maIL4S+7IXwGkP3MZph+UvVmr1WrZXy/YqVRPEVREw1s4u31Mvy2WHtbPggqfC2pSW7dWyvUrFq/hK7EU5kW0XgJbt1bK9SsWr+ErsibU9B0xDaBcmLdurZXuVilfxNcbe9c9gSQjcxpIAXa8I3CZMWrZXKRXL4iux970ZrEESuzUB+twU1TC0bK9SKpaFVba3ssIp28uWz91TGp320gdSLQnupP7nGXDKuRGO+A9atvf4oWV7a5bfjvzGB1s+YOexnaQlpTH0pKFknJhRamg11oVatje+xtjBnbT/stQ9pXHtrBJnnl7h7qnrmadKqRLyHflMXDyR73Z/h9PlxGHc324X/LaABikNmHrOVFrVDa/ibKyJv8QO7uR9wVNhT2msDC3bq1T8cBkX4xaM4+eDP+M0zlLL8hx5FBwrYNT8UXww/AMaptScb+HxNcaulFIh2p2zm+EfDmfDgQ3lkrqXCxfH7Md4bd1rVRxddGliV0rVOPvy9nHFx1eQeTQzaFuHy8HcLXOxBzurPY7E51CMUiouGWNYvW810zdMZ93+dQB0a9SNsZ3H0q1Rt4gdyHx65dMcLToacnuXcXGg4AAnpvq6plD8icvEvvPoTqZvnM7H2z4mz55HrYRaDGs7jDGdxtCybsvqDk8p5YPD5WDSkkks3rWYAkcBBveMvC93fMmy35dxTqtz+PuZf8ciFnYe3ck7P7/Dkl1LcBon7U9oz+hOo+nZpGfQ5H+s6BiLdizyO/zii8u4sFniMh36FHevZMmuJdy5+E4cTkfx0e1cey5zNs9h3q/zeHrA05zV4qxKrbtv374Rr/CYmZnJsmXLuOqqqyK6XqVinTGGn7J/4q2Nb7H50GYOFR4ipzAHJ6UTrsGQ78jn8+2f06hWI9IS03jpp5dKzWD5Ped3vtv9HT0b9+SZs58hyZrka5MAbD28lQRLAoXOwpBjrZ9SnwbJDSr3QmNQXI2x7zy6kzsX30mBo6D4DfdyGAcFjgLuXHwnO4/urNT6tWyvlu1VkVHkLOLmL29m/Ofj+Xz752QezeRI4ZFySb2kAmcBMzbM4D8//YdCZ2Gpv3Fv8v9h7w/cv+T+iMaabE1mXOdxcTWfPZi4SuzTN07H4QyceBxOBzM2zqjU+rVsr5btVZVT5Czik22fMG7BOEZ8OIJz/nsOy39fTr4jv3jIJRQO46DAWeB3eaGzkK93fc3OY/47b+3qtcMerAqshxUr6Q3TubzD5SHHGA/iaijm420fl+upl+UwDj7e9jEPnPFAWNvSsr1atleF5tfDv3LtZ9eS78gnz5EX9e25jIu5W+Zy62m3+lxeN7Eug1sN5rPMz4KOsw8/eTgPnvEgCcFqUMWZsBK7iMwCOnju1gMOG2O6hx2VH3n20D40ufbcsLelZXu1bK8K7kD+AcZ8OoajRUcr1DMPh8PlYMexHQHb3NnzTpb/vpwjRUdwGVe55YmWRMZ0HuP3n0O8CyuxG2OKs4GIPAUcCTuiAGol1AopaacmpAZtE4yW7dWyvSq4d39+t8LDLZFQJ7FOwOVNUpvw3rD3uOvru9h6eCtO48TpcpJiSwHgxu43MrrT6KoItVpEZIxd3Nntz8B7kVifP8PaDsMmgf8X2cTGsLbDohmGlu1VymPWL7MocoXWiQnGJjYsEjwl1bLV4vw25wdt17x2c2YOm8l7F7zHrT1u5cbuN/JQn4dYfMVixnQeU6MOlpYVqYOnZwF7jTFb/DUQkfEislJEVmZnZ1dqI2M6jcFmDZLYrbao/yfWsr0qHhhjcLpCn8tdUQ6Xg6OFoZ8EFEiSNYlWdVtxUbuLAk5lFIQTkk8g48SMkNd9ygmnMC59HBO6TeCCtheQbKv530SDlu0VkUWAr9OxHjDGzPO0+Q+w1RgTUlWucMr2+prHDu7/9jarLax57KHQsr3Hj3gs21voKOTVda8yf9t8duXswmBoXKsxV3e8mss7XB6RYUovYww93upRoROByrKIhXpJ9RjVcRQjO47EIhbGLRjHlsNbys1Dt4qV1IRU3hn6Dm3S2oQZfXyKWNleY8w5QTZkAy4BeoYeXuWd1eIsPrjwA2ZsnMHH2z4m155LakIqw9oOY3Sn0XrmqTouOVwOXlzzIm+sf6Ncot2Xt4/nf3yemb/M5O2hb0esiqGIcHrT01n2e+XO/7CJjbeGvkWnBp1KDcFM/9N03tr4FjM2ziDPnodFLLiMi2Fth3FDtxtqzGn/0RT2hTZE5HxgkjFmQKjPCetCG0pVkXj5TDpdTm796laWZS0LOB3YJjban9CeWRfOiti2f9jzAzcuujHg3HNfEiwJ9GvWj+cHP++3jcu42J27G6fLSeNajY+LIZRgqvJCG/9HlA+aKnU8ybPn8eHWD5mxcQZ7cvdgtVjp3aQ316RfQ0bT8mPL83+bzw97fgjpHI/fjv7G+v3rSW+YHrBtqHqf2JurO13NWxvfCjm5p9hSaFevHU/0fyJgO4tYaF67eSTCPO6EndiNMWMjEIdSCsjOy+bqT6/mQP6B4kTpdDpZ+vtSVu9bzSUnX8K9GfeWmtHx+rrXyXfkh7T+Qmch87bOi1hiB7j1tFvpUL8DL615iaycLBIsCdhddtrUbcOAlgNYsXsFGw9sxGA4Ke0kxqWPY0ibITXupKBYEldnnipVkxljGP/5ePbm7vXZ+8535DNn6xza1mvLnzv8ufixUGqOe7mMi/35+yMVcrEhbYYwpM0Qdh3bxZHCI9RPrk/T2u4T7+gR8c2pIOI6se/evZszzzyTpUuXcuKJekBFVa1fDv7CjI0zWLXXff5Ct0bdGN15NJ0bdC5us/nQZnbn7CbFlkK3xt0CTuVbuXclWTlZAYdUChwFTP1pKpe1vwyLWHC4HO7ee4iHyixYaJLaJLTGldCiTgta1GkRtfWr0MRVEbCyHn30UTIzM3n00UfDXtfOnTsZNGgQnTp1onPnzjz77LM+291xxx0888wzxfeHDBnCddddV3z/rrvu4umnnw47nnAdPnyYl156KaS2ffv2rdC6Bw4cSNmD37FkzZo1zJ8/P2rrN8bw+PePM2r+KD7Z9glZOVlk5WSxIHMBYz8dy4PfPsjXO7/mwrkXMmr+KO5bch+3fXUb/Wf256mVT1Hk9H1Cz9wtcylwBB+nzrXnsmH/BgBqJ9Smlq1WyLEnWhMZcfKIkNur+BS3iX337t288cYbuFwu3njjDfbs2RPW+mw2G0899RQbN27ku+++48UXXyxVL8WrX79+xeV9XS4X+/fvZ8OGDcXLly1bVuFEGQ5/ZXYrktijUa7Yn6ooC1yZxF6RuF5e+7I7CTsLSk0tdBkXBc4CPtn2Cbd/dTuZRzPJd+STY88hx55DniOP935+j2s/u9bnZdj25e0L6dR8i1g4VHgIcE85vPLUK0m0JAZ9nlWsdGzQkVPr+z85TdUMcZvYH3300eLaJ06nM+xee9OmTTnttNMAdxmAjh07kpWVVa5d3759Wb58OeAuj5uenk6dOnU4dOgQhYWFbNq0idNOO42cnBwGDx7MaaedRpcuXZg3bx4Aubm5XHDBBXTr1o309PTi8gT33XcfnTp1omvXrtx9992A/3K9Dz/8MFdffTX9+vXj6quv9lmm97777uPXX3+le/fuTJw4EYB///vf9O7dm65du/LXv/61+DWFUq64rLfeeovu3buTnp7OihUrAP9lhN98800uuugizj77bAYPHszo0aP58MMPi9c1cuRI5s2bh9Pp5O677yY9PZ2uXbsWlzVYtWoVAwYMoGfPngwZMoTdu3cD7m8O9957LxkZGbRv354lS5ZQVFTEQw89xKxZs+jevTuzZs0iNzeXa665hoyMDHr06FH8XpSNKxT5jnymrZ8WcAaIwzj8nrRT6Cxk08FNvLnxzXLLQp1f7jIu6iXVK74/quMo6ibVDXg6viC0rtOa58/2P71Q1RxxOcbu7a17C20VFRXxxhtvMHny5IiMtWdmZvLjjz9y+umnl1vWrFkzbDYbO3bsYNmyZfTp04esrCyWL19OWloaXbp0ITExEYvFwty5c6lbty779+/njDPO4KKLLmLBggU0a9aMTz75BHDXezlw4ABz587l559/RkQ4fPgwgN9yvQAbN27k22+/JSUlhVtuuaVcmd5//vOfrF+/vrjOy8KFC9myZQsrVqzAGMNFF13EN998Q//+/Uu9vkDlikvKy8tjzZo1fPPNN1xzzTWsX7/ebxlhgNWrV7N27Vrq16/P4sWLmTJlCiNGjODIkSMsW7aM6dOn88orr5CZmcmaNWuw2WwcPHgQu93OLbfcwrx582jUqBGzZs3igQceYNq0aYC7p71ixQrmz5/PI488wqJFi/jb3/7GypUreeGFFwC4//77Ofvss5k2bRqHDx8mIyODc845p1xcofhqx1cI4dUYKXQW8vbGt7mm8zVYLdbixy8+5WK+2vlV0NK3KbaUUrNa6iXX492h73LDohvYk7un3AyZZGsyd/a8k0vaXxJwjF/VHHGZ2Ev21r28vfYXX3wxrHXn5ORw6aWX8swzz1C3bl2fbbyX0Fu2bBl33nknWVlZLFu2jLS0NPr16we4x2Hvv/9+vvnmGywWC1lZWezdu5cuXbpw1113ce+99zJs2DDOOussHA4HycnJXHvttQwbNoxhw9xFzPyVzwW46KKLSElxV6oLpUzvwoULWbhwIT169Ch+nVu2bCmX2AOVKy7JW7q4f//+HD16lMOHD3Ps2DGfZYQBzj333OLkOWDAAG688Uays7OZM2cOl156afE/gwkTJmCzuT+W9evXZ/369axfv55zzz0XcL/P3jLHAJdccgkAPXv2JDMz0+f7tXDhQj766COefPJJwF3ZcseOHeXiCsWevD1+x8grIt+Rz66cXbSu27r4sYwTM2hSqwk7ju3w2+NPtiYzvuv4cr3zprWbMm/4PFbvW83/fv0fBwsO0qx2My455RLan9A+7HhVfIm7xF62t+4ViV673W7n0ksvZeTIkcUJY+fOnVx44YUATJgwgQkTJhSPs69bt4709HRatmzJU089Rd26dRk3bhwA77zzDtnZ2axatYqEhATatGlDQUEB7du3Z/Xq1cyfP58HH3yQwYMH89BDD7FixQq++OIL3n//fV544QW+/PLLgOVzS5bZ9VWmt23btqXaG2OYNGkSN9xwQ8B9EKhccUllK+OJCJMnT/ZZRrhsvACjR4/m7bffZubMmbzxxht+4zHG0Llz5+LhL3/xBorVGMOcOXPo0KFDqce///77cnEFk2pLxWqxBr2SVzAWLOXG2UWEV857hVHzR3Gk8Ei54Z4UWwrD2g7jylOv9LlOEaFnk570bFIl1T1UDIu7MXZfvXWvcMbajTFce+21dOzYkTvvvLP48ZYtW7JmzRrWrFnDhAkTAHeP/eOPP6Z+/fpYrVbq16/P4cOHWb58efGB0yNHjtC4cWMSEhL46quv2L59OwC///47tWrVYtSoUUycOJHVq1cXXzRj6NChTJkyhZ9++gkIvXyurzK9ZcsFDxkyhGnTphX3+LOysti3b1+l9hX8Ubr422+/JS0tjbS0NL9lhH0ZO3Zs8eyiTp06Ae7e88svv1ycoA8ePEiHDh3Izs4uTux2u73UwWpffL32559/vvh4wY8//ljBV/uH/i36+z3uUBF2Y/dZ8+TE1BOZO3wuN3W/iSa1miAIVrHS+8TeTBk4hclnTK7R5WZVZMRVYvfXW/fy9torM0Nm6dKlvPXWW3z55Zd0796d7t27+51Z0aVLl+Jx85KPpaWlFVd9HDlyJCtXrqRLly7MmDGjuEzuunXrig90PvLIIzz44IMcO3aMYcOG0bVrV84888zi6ZKhls/1Vaa3QYMG9OvXj/T0dCZOnMh5553HVVddRZ8+fejSpQuXXXZZWHXik5OT6dGjBxMmTOD1118H/JcR9qVJkyZ07Nix+BsOwHXXXUerVq3o2rUr3bp149133yUxMZH333+fe++9l27dutG9e/egs3gGDRrExo0biw+eTp48GbvdTteuXencuTOTJ0+u9OtuWrspGU0zwjpr0iIWBrcaTO1E31ekqpNYh7HpY1l0+SJ+Gv0Ta0avYdqQafRr3k+TugpJ2EXAKqOyRcBuvPFGXn/99YBXJ0pMTOS6664Le6xdRVdeXh5dunRh9erVpKWlVXc4Pvn7TB4pPMKVn1zJ3ty95S4yUfJiEf4uQJFqS2X2hbNpVbdV5INWNVqoRcDiqsf+0UcfBb3kXFFRUfF0NhWbFi1aRMeOHbnllltiNqkHkpaUxuxhsxnTeQx1EuuQYkshxZZCakIqozqNYv4l8+nTrA9J1iSs8sesl1q2WtRPrs+086dpUldRFVc9dqWqUiifSYfLwZ5c99Bfk1pNSLD+MUSz9dBWZm+eTeaRTGon1mboSUMZ2HIgNkvczVlQMaIqy/ZGjDFGxxBVTAi1w2Oz2PzWRjn5hJO5//T7IxmWUiGJmaGY5ORkDhw4EJEZB0qFwxjDgQMHfE4zVSoexEyPvUWLFuzatYvKXuhaqUhKTk4uPlFLqXgTM4k9ISGBk046qbrDUEqpuBczQzFKKaUiQxO7UkrVMJrYlVKqhqmWeewikg1sr/INR1ZDIPIXj6xZdB+FRvdTcLqP3FobYxoFa1Qtib0mEJGVoZwocDzTfRQa3U/B6T6qGB2KUUqpGkYTu1JK1TCa2CvvleoOIA7oPgqN7qfgdB9VgI6xK6VUDaM9dqWUqmE0sSulVA2jiT1MInKXiBgRaVjdscQiEfm3iPwsImtFZK6I1KvumGKFiJwvIr+IyFYRua+644lFItJSRL4SkY0iskFEbqvumOKBJvYwiEhL4DxgR3XHEsM+B9KNMV2BzcCkao4nJoiIFXgR+BPQCbhSRDpVb1QxyQHcZYzpBJwB3KT7KThN7OGZAtwD6BFoP4wxC40x3itbfwdoLVy3DGCrMWabMaYImAkMr+aYYo4xZrcxZrXn9jFgE9C8eqOKfZrYK0lEhgNZxpifqjuWOHIN8Gl1BxEjmgM7S9zfhSasgESkDdAD+L56I4l9MVOPPRaJyCLgRB+LHgDuxz0Mc9wLtJ+MMfM8bR7A/bX6naqMTdUMIlIbmAPcbow5Wt3xxDpN7AEYY87x9biIdAFOAn7yXKO1BbBaRDKMMXuqMMSY4G8/eYnIWGAYMNjoiRNeWUDLEvdbeB5TZYhIAu6k/o4x5oPqjice6AlKESAimUAvY4xWnytDRM4HngYGGGP0uoceImLDfTB5MO6E/gNwlTFmQ7UGFmPE3XOaDhw0xtxe3fHECx1jV9H2AlAH+FxE1ojI1OoOKBZ4DijfDHyG+4DgbE3qPvUDrgbO9nx+1ojI0OoOKtZpj10ppWoY7bErpVQNo4ldKaVqGE3sSilVw2hiV0qpGkYTu1JK1TCa2JVSqobRxK6UUjXM/wOLEmgPgpS5HgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.figure(1)\n",
    "for (x_i, b_i) in zip(measures_locations, measures_weights):\n",
    "    color = np.random.randint(low=1, high=10 * N)\n",
    "    pl.scatter(x_i[:, 0], x_i[:, 1], s=b * 1000, label='input measure')\n",
    "pl.scatter(X[:, 0], X[:, 1], s=b * 1000, c='black', marker='^', label='2-Wasserstein barycenter')\n",
    "pl.title('Data measures and their barycenter')\n",
    "pl.legend(loc=0)\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
}