plot_free_support_barycenter.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 2D free support Wasserstein barycenters of distributions\n",
    "\n",
    "Illustration of 2D Wasserstein barycenters if distributions are weighted\n",
    "sum of diracs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": 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\n",
    "from src.main import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate data\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": 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": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# len(measures_locations[0])\n",
    "# b.shape\n",
    "# X_init.shape\n",
    "\n",
    "# type(measures_weights)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute free support barycenter\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> \u001b[0;32m/mnt/data1/jjs/anaconda3/envs/tutorial/lib/python3.7/site-packages/ot/lp/__init__.py\u001b[0m(509)\u001b[0;36mfree_support_barycenter\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m    507 \u001b[0;31m            \u001b[0mT_i\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0memd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmeasure_weights_i\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mM_i\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m    508 \u001b[0;31m            \u001b[0;32mimport\u001b[0m \u001b[0mpdb\u001b[0m\u001b[0;34m;\u001b[0m \u001b[0mpdb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_trace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m--> 509 \u001b[0;31m            \u001b[0mT_sum\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mT_sum\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mweight_i\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1.\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mT_i\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmeasure_locations_i\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m    510 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m    511 \u001b[0;31m        \u001b[0mdisplacement_square_norm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msquare\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mT_sum\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  np.reshape(1. / b, (-1, 1)).shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10, 1)\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>   weight_i.shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** AttributeError: 'float' object has no attribute 'shape'\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  weight_i\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3333333333333333\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  T_i.shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10, 19)\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  measure_locations_i\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[-1.74405746, -2.33088241],\n",
      "       [-1.53027203, -0.38910966],\n",
      "       [-2.46831313, -3.67876652],\n",
      "       [-1.71605865, -2.60562761],\n",
      "       [-1.55449653, -1.8179163 ],\n",
      "       [-1.73519639, -2.8968618 ],\n",
      "       [-1.5627412 , -3.85464208],\n",
      "       [-2.14509887, -3.76208661],\n",
      "       [-1.92377825, -3.79063364],\n",
      "       [-2.29606906, -4.42458486],\n",
      "       [-1.83793195, -2.12813137],\n",
      "       [-1.68909875, -1.97687479],\n",
      "       [-1.42504216, -1.57986962],\n",
      "       [-2.18904942, -3.58416849],\n",
      "       [-2.03579307, -2.42157027],\n",
      "       [-1.72013441, -1.73526201],\n",
      "       [-2.31027248, -5.12986654],\n",
      "       [-2.0682932 , -2.95264218],\n",
      "       [-1.87229562, -2.82886523]])\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  measure_locations_i.shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(19, 2)\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  np.matmul(T_i, measure_locations_i).shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10, 2)\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  np.reshape(1. / b, (-1, 1)) * np.matmul(T_i, measure_locations_i).shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.],\n",
      "       [100.,  20.]])\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  np.matmul(T_i, measure_locations_i)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[-0.18446482, -0.22869986],\n",
      "       [-0.20671176, -0.35879348],\n",
      "       [-0.1453548 , -0.16322875],\n",
      "       [-0.24105302, -0.40314513],\n",
      "       [-0.16715926, -0.13910934],\n",
      "       [-0.22755264, -0.37047855],\n",
      "       [-0.17090573, -0.24158213],\n",
      "       [-0.20082868, -0.28151532],\n",
      "       [-0.17619252, -0.28836052],\n",
      "       [-0.15627412, -0.38546421]])\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  np.reshape(1. / b, (-1, 1)) * np.matmul(T_i, measure_locations_i)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[-1.84464825, -2.28699856],\n",
      "       [-2.0671176 , -3.58793481],\n",
      "       [-1.45354801, -1.6322875 ],\n",
      "       [-2.41053023, -4.03145135],\n",
      "       [-1.6715926 , -1.39109337],\n",
      "       [-2.27552645, -3.70478553],\n",
      "       [-1.70905731, -2.41582132],\n",
      "       [-2.00828682, -2.8151532 ],\n",
      "       [-1.76192518, -2.88360523],\n",
      "       [-1.5627412 , -3.85464208]])\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  np.reshape(1. / b, (-1, 1)).shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10, 1)\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  np.matmul(T_i, measure_locations_i).shape\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10, 2)\n"
     ]
    },
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "ipdb>  type(np.matmul(T_i, measure_locations_i))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ndarray'>\n"
     ]
    }
   ],
   "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": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_bary = (measures_locations, measures_weights, X_init, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot data\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABBe0lEQVR4nO3deXhTZfbA8e9J0zYtm4CA7BVk3woCsikgKA4woqLjhg4q+sNdR3FBVFxnHHQUt3EDHLcRB1QQnBEQEGQRAUE2kZ2yl0ILbWmzvb8/kpTSJm3apE3Tns/z5CHJvbn35DacvHnve88rxhiUUkpFL0ukA1BKKRUaTeRKKRXlNJErpVSU00SulFJRThO5UkpFOU3kSikV5TSRKxUEEflQRJ4vq/ULvPYdEXmyBOuPFpEfS7MvVTloIi9DIrJbRE6JyEkRSReR5SIyVkSCOu4ikiQiRkSsZR2rKr1wJ1JjzFhjzHPh2l40COWLT2kiLw9/NMbUAJoDfwMeBaZENqSKS7+0ilbWxydaj3+0xh0umsjLiTEmwxgzG7gW+LOIdAQQkWEi8ouInBCRFBGZmO9lS7z/potIpoj0FpGWIrJQRNJE5KiIfCoiZwXar7dFf5eIbPP+MnjOu43l3n1+ISJx+dYfLiLr8v2C6Jxv2WMissO7nc0icmW+ZeeJyA8ikuGNa7r3+UK/KkRksYiM8d4fLSLLRORVEUkDJopIvIi8LCJ7ReSwt6shwbv+2SIyxxvfMRFZGugXjohM9h7TEyKyRkQuzLdsove9f+R9P5tEpHu+5V1FZK132XTAFmAf7YB3gN7ev1F6vsW1RWSudxs/iUjLfK9rKyLzve9hq4j8Kd+yvNapiAwQkX0i8qiIHAKmBf5Ty5ve4/+biAzKt+AWEdnijWOniPxfvmWFti8iG0Xkj/nWifX+Tbt6H/fzfjbSvcd3tPf5ov5uvv08JCJHROSgiNziXXYHcCPwiPcYfuN9vpGIzBSRVBHZJSL3Ffj7zRCRT0TkBDA6wHGpGowxeiujG7AbGOzn+b3And77A4BOeL5UOwOHgSu8y5IAA1jzvfY84BIgHqiHJ9m/VkQMBpgF1AQ6ALnA90ALoBawGfizd92uwBHgAiAG+LP3PcR7l18DNPLGei2QBTT0Lvs38IR3mQ3oV8R7WAyM8d4fDTiBewErkAC8CswG6gA1gG+Av3rX/yuexBnrvV0ISID3Pgqo693uQ8AhwOZdNhHIAYZ63+tfgZXeZXHAHuBB7z6uBhzA8wH2Mxr4scBzHwJpQE/v/j8FPvcuqwakALd4l3UFjgLt8732+XyfDyfwkvdvnhBg/8588V4LZAB1vMuHAS0BAfoD2UC3QNsHHgGm59v+CGCD935z4CRwvXdfdYFk77Ki/m6+/Tzrfd1Qbxy1C75n72MLsAZ4yvv3aAHsBIbk+/s5gCu86xY6LlXpFvEAKvONwIl8JfBEgNe8BrzqvZ9EgSToZ/0rgF+KWG6AvvkerwEezff4FbxfBMA/gecKvH4r0D/AttcBI7z3PwLeA5oUWKfQe6BwIt+bb5ng+YJome+53sAu7/1n8XwxnVeKv8dxoIv3/kRgQb5l7YFT3vsXAQfI9wUBLKfkifyDfI+HAr95718LLC2w/rvA0/lemz+R2/F+ARWx/4LxrgJuCrD+18D9gbaP58v6JFDT+3gG8Ij3/uPAV362WdzfbQBwqsDn4AjQq+B79j6+IP/nIt++p+X7+y0p6Wegst60ayUyGgPHAETkAhFZ5P35mAGMBc4O9EIRaSAin4vIfu9Pyk+KWt/rcL77p/w8ru693xx4yPuTOd3bTdAUz39sRORmOd3tkg50zLfvR/D8Z17l7aa4tZiY8kvJd78ekAisybef/3mfB5gEbAfmebsJHgu0URF52NulkOHdTi3OPFaH8t3PBmzi6QJqBOw33ozhtacE7yfQ9vMf5wsKHOcbgXMCbCfVGJNTzL78xev7u/1BRFZ6u3HS8Xyp5D8OZ2zfGHMAWAaMFE+33R/w/KIAz+dhh5/9F/d3A0gzxjjzPc5/TApqDjQqcIzGAw3yrZPi95VVUJU+QRAJItIDTyL3jXL4DHgT+IMxJkdEXuP0fzJ/pSlf9D7fyRhzTESu8L4+HFKAF4wxL/iJuznwPjAIWGGMcYnIOjzJG2PMIeB277r9gAUisgTPT3zw/Cc/4b1fMGHlf59H8Xy5dDDG7C8YhzHmJJ5ukofEc55hoYj8bIz5vkC8F+L5chkEbDLGuEXkuC/eYhwEGouI5EuOzfCfwArGH4wU4AdjzCVBrh/M9v3FO1tE4oGZwM3ALGOMQ0S+5szj4G/7/wLG4MkRK/L9LVLwdBcVVOTfLQgFY0jB05pvVYLXVFnaIi8nIlJTRIYDnwOfGGM2eBfVAI55k3hP4IZ8L0sF3Hj6B8m3fiaQISKNgXFhDPN9YKz3V4KISDXxnIytgadf13hjwnuiqmO+93eNiDTxPjzuXddtjEkF9gOjRCTG21JvSQDGGLc3jldFpL53241FZIj3/nDxnFgVPF8SLjzHqKAaePpkUwGriDyF5zxBMFZ4X3uf90TfVfhPXj6HgSaS76RxMeYArUXkJu/2Y0Wkh3hOnJZW/XzxXgO0A77F078cj+c4OEXkD8ClQWzva6AbcD+ebjOfT4HBIvInEbGKSF0RSS7u7xaEw5z5OV8FnBTPSdgE72eno7chpArQRF72vhGRk3haGE8A/8BzksvnLuBZ7zpPAV/4FhhjsoEXgGXen5e9gGfw/AfLAOYCX4YrUGPMajyt6jfxJOPteEcDGGM24+lPX4HnP10nPD+/fXoAP4lIJp4TXvcbY3Z6l92O5wsnDc8J1+XFhPKod98rvd1HC4A23mWtvI8zvbG8bYxZ5Gcb3+H5af87nm6GHIL8KW6MsQNX4Xnvx/D0aRd1nBcCm4BDInI0iO2fxJNMr8PTt32I0ycbS+snPMfmKJ7PzNXGmDTvvu7D87k6jqehMDuIGE/hacmfS773bozZi6dr5iE8x2Yd0MW7uKi/W3GmAO29n/OvjTEuYDiQDOzyvq8P8HSPqQLkzG41pZTy8P6KaW2MGRXpWFTRtI9cKVWIiNQBbgNuinQsqnjataKUOoOI3I6nG+q/xpglxa2vIk+7VpRSKsppi1wppaJcRPrIzz77bJOUlBSJXSulVNRas2bNUWNMvYLPRySRJyUlsXr16kjsWimlopaI+L3CWLtWlFIqymkiV0qFxBhDliMLh9sR6VCqLB1HrpQqla3HtvKvTf/iu93f4TZu3MbNebXP47aOt3FJ80uIjYmNdIhVRkSGH3bv3t1oH7lS0Wvqhqn8c/0/cbgduIzrjGUJ1gSa1WjGlCFTqBWvV9SHk4isMcZ0L/i8dq0opUrky21f8s/1/yTHlVMoiQOccp5iZ8ZOxswbg9Pt9LMFFW6ayJVSQXO4Hbz888vkuIouj+5wO9h7Yi8/7PuhnCKr2jSRK6WCtmjvItzGX9XgwrKd2UzdMLWMI1KgiVwpVQI/HfyJLGdW0OtvSttUhtEoH03kSqmgFdelUpDbuNF6TmVPhx9WQpsPnODNRdvIcbgZ1asZF7dtUPyLlApCk+pNiLXEBj1mvGZcTTyTOamypC3ySmb7kZNc/c5yvt1wiIW/HeGuT9fyv40HIx2WqiQuP+/yoBNzrCWWka1GlnFECjSRVzoz1uznlP30kLAch5t3l+ws4hVKBa9x9cacX/98rFL8j/kYieH6dteXQ1RKE3klYxEo2GCyBDVxvFLB+dtFf6NuQt0ik7ktxsbEPhM5p9o55RhZ1aWJvJK5tkdTEuKsecncFmvhnovPi2xQqlKpY6vDF3/8gr6N+xIXE0ecJQ4AQUi0JnJOtXN4ZcArDGsxLMKRVh16iX4ltCM1k/eX7CTH4eJPPZrSp+XZkQ5JVVJHso/wv13/43D2YRKtifRs2JPuDbrrCc4yEugSfU3kSikVJbTWilJKVVKayJVSKsppIldKqSinV3aWEWMMy3ek8d6SnWw6kIHDZahps3Jlt8aMuqA59WvaIh2iUsrrWM4xvt72NZuPbcbhctCkRhOuOO8KWtVuFenQgqInO8vAtsMnGT3tZ45n28m2n1mvOd5qwQDXdm/KxMs7EGPRs/tKRUqmPZNnVjzDwpSFCEKuKxfwXMxktVhpUasFz/d7nta1W0c4Ug892VlOth46yRVvL+NA+qlCSRwg1+nG7nQzY00Kd326BrdbCwopFQkn7Ce4bu51LNy7ELvLnpfEAVzGRa4rly3HtnDTtzfxa+qvEYy0eCEnchGxicgqEVkvIptE5JlwBBaNnC43o6b8RFaui+LS8ymHmyW/H+VfK3aXR2hKqQIeWvwQBzIPYHfbi1wv25nN2PljOWE/UU6RlVw4WuS5wMXGmC5AMnCZiPQKw3ajzoItR8jODX5qq1MOF28v3qGtcqXK2a6MXfxy5Jegqzg63U5mbZ9VxlGVXsiJ3Hhkeh/Gem9VMjO9+8MOsvx0pxQlO9fJj9uPllFESil/Pt3yKS538P9XT7lO8eGmDytsbfWw9JGLSIyIrAOOAPONMT/5WecOEVktIqtTU1PDsdsKZ/PBkv/0yra7WJ+SHv5glFIB/XzoZ5ymZBNDH8s5RqYjs/gVIyAsww+NMS4gWUTOAr4SkY7GmI0F1nkPeA88o1bCsd+KxuEKbi7D/AyQWYLuGKVU6OyuovvF/YmRGHJdudSgRolel2nPZPaO2czcNpOM3AxqxNVgRMsRXNnqSmrF1ypxHP6EddSKMSYdWARcFs7tRouE2JgSv8ZqEWpXiyuDaJRSgdS21S7xa5xuJzXiSpbEfzr4E4P+M4jX1r7G78d/53D2Ybanb+etdW8x+D+DWZyyuMRx+BOOUSv1vC1xRCQBuAT4LdTtRqOL29anpMPCrRZhUNv6ZROQUsqvK8+7kgRrQole071Bd+Jj4oNef1PaJu75/h6yndmccp46Y1mOK4ccVw7jfhjH6kOhX1MTjhZ5Q2CRiPwK/Iynj3xOGLYbdW6/qAXx1pK1yls1qE6rBiX7lldKhWZYi2ElOnGZYE3g1k63lmgfL//8crGTVee4cvjrqr+WaLv+hGPUyq/GmK7GmM7GmI7GmGdDjirCjmXZmb3+AF+sTuHn3ceC/oN3bnIWnZvUIs4a3GG1xVp49LJ2oYSqlCqFxNhE7j//fmwxxZfKiLPE0fHsjvRqGPyo6oOZB9lwdENQ6+49sZdtx7cFvW1/tNZKPtl2J49/uYH/bjxErEUwxjNtWs2EWF68qhMD2xTfBfLBn7sz8p/L2ZOWTa4z8MlPW6yFCcPa069VZCd9MMawdu9x9qfn0LR2AslNz9JJAVSVMKrdKE7mnmTqxqnkunIxfkZNJ1gTaFenHW9e/CYWCb7du/X4VmItsWdcLRqIRSz8duy3kOq6aCL3ynG4+NO7K9h2OBO7003+c9pZdhd3frKG16/ryqUdip6DsIYtlq/v7ssLc7cwc+0+LCJ5l+pbBOKtMTSoGc9Tf2zPxW0blOE7Kt76lHTu/HQNGdkOEDAG6laP472butOuYc2IxqZUebgz+U6qx1Vn5u8z2XtyL/Ex8QiCw+2gcfXG3NrpVoaeOxSrpWKnyoodXTn6z+oUdhzJCtiKznG4efg/61nTtj6xMUV/MyfGWXnhyk6MH9qOr9ftZ11KOtl2Fw1q2Phjl4YVotWbciybG95fWegCpuxjp/jTuyv4/i/9tUKjqvTm7Z7H62tfx23cJFgTmNR/ErGWWOol1COpVlKpt9uqdqugrxp1G3fIRbk0kXu9u2QnpxxFX+nlchvmbz7M0E4Ng9pmtXgrN17QnBsvaB6OEMPqgx93YQ/wpWV3uvnXit2MG9K2nKNSqnytOLAi74Sk1W2ljq0ObeuE/rlvXL0xHep2YO2RtcWu27RGU9rUaRPS/rT6IZ4Evf/4qWLXy7K72Lg/oxwiKnvfbzmMI0CNl1ynm+82Hi7niJQqf0POHYItxkY1azVq22pzbq1zw7btcT3GFXsy1RZj47Gej4W8L22RA4LnpGZxg1MEiCumWyVaFNexo+c7VVXQq2EvPh/+ObszdnNBwwtKNE68OB3P7sjkgZN5YPEDGGPOGIro64t/sd+L9GzYM+R9aSIHLBahc5OzWFdMzZOEuBh6t6xbPkGVscs6nsOHy3fjcBX+9oq3WhjeuVEEolKq/LU8qyUtz2pZJtvu07gP86+ez1fbvmLmtpmcsJ/Iu0R/ZOuR1LHVCct+NJF73TWgJQ9MX+d3MgjwtGDPrh5Pz3PDc+DLwym7i2PZdpwuN7USYqmVEJt3kvW2fi34/OcUHK4z67yIeEoNjOrVLBIhK1Xp1IqvxeiOoxndcXSZ7UMTudcl7RtweZdGzFp3oNBJT4t4Tlx+8OfuER9tUhxjDD/tOsZ7S3aydFsqVosFEU9Br8a1E7izf0su79KYc2rZmDG2D/d8tpZ9x09hjREcLjdJdavx9o3dqFs9fD8xVfRJz0ln78m9xFpiOe+s84iNiY10SKoIOmdnPsYYZqzZx+sLt3HkRC5Wi+B0Gy7rcA4PD2lD0zqJkQ6xSCnHsrl56ioOn8jhlN3/LEWJcZ4SAi+N7MQfuzQG4LdDJziYnkOT2glaLqCK23NiD6+sfoVl+5cRFxOHMQYR4bq21zG2y9iw9iGrkgs0Z6cmcj+MMRzMyCHH4aJBTRvV4iv+D5ftRzK56u1lZOY6CWbCoXir8PjQdozuE76z9Cq6bT22lT//789kO7ILXeUYHxNPq9qt+PCyDzWZR5BOvlwCIkKjsxJoUa96VCTx9Gw71723gpNBJnGAXKfhr99uYdFvR8o2OBUVjDHct/A+shxZfi9Vz3Xlsu34Nt5Z/04EolPF0UReCXyycg+ZOc5ih08WlOs0PD17Y4WdvkqVn58P/Ux6bnqR6+S6cvn8t89xuIK7YlGVH03kUc7lNkxdtpucIgp0FeXwiVx+0anmqrxlB5aR7cwudj1jDNvTt5dDRKokKn6/gQJgR2omk/73G0u2HSXGIvyxcyP+cmlrNuzLINdZsgmf88t1uvlgyU7eHnV+GKNV0SaYKn3g6XYMtoaIKj+ayKPAjtRMLn/zR7Ltrrzuky9Wp7Bw6xGuOb8JpwKMfQ+WtshVm9ptSLAmFJrJpiC7y07TGk3LKSoVLO1aiQIv/e+3M5I4gNNtOJ5lZ+WOtKBPcAaSZdfJn6u6IUlDcLmLbhAIQt/GfUs136UqW5rIo8CP2476PZGZ63SzOy2rxPOEFmQr4fR0qvJJjE3ksnOLnjM9MTaRv5z/l3KKSJWEJvIoEFNEpk6Mt5IQF1oibnOOXgSkYHSH0VilcG9rgjWBhtUa8q/L/hVSjW5VdkJO5CLSVEQWichmEdkkIveHIzB12rBODbH6SeYJsRZu7ZuEJYSyAVYLjO1fNgWDVPGOHz9Oenp6pMMAPJMhPNLjEWItscTFxBFDDAObDuTNi9/ku5HfhVwzW5WdcJzsdAIPGWPWikgNYI2IzDfGbA7DthXw0KVtWLT1COnZjrwZjBJiLbQ5pwbX9WzGoYxcpvy4C7ur5EMQaybE0qeSVHSMRldccQWxsbEsWLAg0qEAcH2767k06VIOZx+mcfXG1IqvFemQVBBCTuTGmIPAQe/9kyKyBWgMaCIPk3o14vnugYv4eOUe5v56kDirhWu7N+Xq7k2It8Ywum8Sn/y0p8SJ3GoRJgxrX+ELgVVWq1atYvXq1Z4JsNeupVu3bpEOCYC6CXWpm6Bf7tEkrH3kIpIEdAV+Cud2FZyVGMe9F7fifw9cxOx7+nFjr+bEe09SNqhp41+39iQhNvi+cqtFuK3fuVzVrUlZhVyp2e12Ro4cyaFDh0q9jXHjxnHq1ClycnJ45JFHwhidqmrClshFpDowE3jAGHPCz/I7RGS1iKxOTU0N126VV7dmtfnP2N7UToylWhEnP2ME4mKEcUPa8PjQduUYYeXy4Ycf8tVXX/H888+X6vX5W+PGGJYvX87atcXP76iUP2GpfigiscAc4DtjzD+KW7+iVz+MZrlOF//beIh3ftjBztQs4qye72q3MVgtFv7cpzmjLmhO/ZpFzyWoArPb7TRp0oTU1FQSEhLYuXMn55xzTom20b9/f5YuXZpX50ZEuPjiiytMX7mqmAJVPwy5j1w8HaxTgC3BJHFVtuKtMYxIbsyI5MakHMvmyMlcHC43ZyXGcl696lgryZyjkfThhx+Sne2pS+J2u3n++ed58803g359/ta4T/5WeUXpK1fRI+QWuYj0A5YCGwDf2bbxxphvA71GW+QqWuVvjfuUtFVesDXuo61yVZwyq0dujPnRGCPGmM7GmGTvLWASVyqa5W+N+/ha5cHw1xr30b5yVVr6O1upINntdiZMmEBWVtYZz+fm5jJ16tSgRrD4RqoEoiNYVGloIlcqSP5a4z7BtMpXrVrFsmXLipzIwxjD4sWLWbNmTUixqqpFy9gqFYRArXEfX6t8woQJAfvKT5w4QY8ePYqdkUlEOHGi0AjeiEvNTuXH/T+SkZtBXEwcrWq3onuD7npBWQWgiVypIBTVGvcpbgTL4MGDGTx4cFmEV6Y2Hd3E2+veZuWhlcRIDA6XA4tYsFqs1IirwegOo7m27bXEWmIjHWqVFZZx5CWlo1ZUtGnXrh27d+8mLi4u4DpOpxOXy0VWVhYxMZWjNPCcHXN4ZsUz5Lpy/U7KDGCLsdG+bnv+OfifJMYmlnOEVUuZjSNXqir4/vvvOXbsWLHrVatWrdIk8R/3/8gzK54hx5VT5Ho5rhw2Ht3IfYvu471L3sMieuqtvGkiVyoIjRo1olGjRpEOo9wYY3h6+dPFJnEfu9vOr6m/smz/Mi5scmEZR6cK0q9OpSoJpzN8U/b9dOgnMu2ZJXrNKecpnlnxDA6XTs5c3qIikW87fJLxX25g2OtLGfLqEm6Z9jPfbzmMK9TJKpWqJLZv3079+vVJSUkJy/Y+/+1zsp1Fn9z153D2YW76703YXfawxKGCU6G7Vo5m5nLnJ2vYsD8Dh8uNr9z21sMnWbU7DZs1htev70rf886ObKBKRdiECRPIyMhg4sSJTJkyJeTtpZws/RfC78d/59kVz/J8v9JVhlQlV2Fb5GmZuQx/40d+2ZtOjuN0EvfJynWRlmXntn/9zKLfjkQmSKUqgO3btzNr1izcbjefffZZWFrloYxmc7gd/HfXfzmeczzkOFRwKmwiv+/fv5B2MhdnMd0nOQ43d326luNZ+lNOVS7ffPMNe/fuLXa9CRMm5PWPu1wuJk6cGPK+G1ZvGNLrRYQvt31Z7HrpOelM/206r699nX9t+hcHMw+GtN+qqkIm8r1p2azecxxHkH3gBsP0n4v/wCsVLY4ePco111zDvffeW+R6vta4L5E7HI6wtMqvbXMtidbSjwnPdeWy4sCKgMtdbhd//emvDJ4xmJdXv8z7G97n9bWvM/zr4dy/8H6yHSXvn6/KKmQi/3jlbtwl+GmX43AzZdnusgtIqXL20ksvATB//nw2bw48/W3+1rhPOFrlfRv1xWYNbfKRooYuTlg2gS+3fUmuKzdvPbvbjt1l58cDP3LbvNtwuHX0S7AqZCLfeOAEDlfJ+uiOZnomUFAq2h09epS3336b3Nxc7HY7jz/+uN/1CrbGfcLRKo+xxDD+gvHYYkqfzOvY6vh9fkvaFhbsWRAw0dtddnak72DBHq3LHqwKmchLM6xQSvk6pWbNmhXSJMrh9tJLL+F2exolLpcrYKvcX2vcJxyt8iFJQ3jw/AdLlcwTrYmMaDnC77JPt3yK3V30Oa1TzlNM2zitxPutqipkIk+qm4ilhAXVEmJjsJVgFnmlAA4ePMjVV1/NuHHjIh0KcLo1npNzurXqr1UeqDXuE66+8hva3cAbg96gWY1mJXpdrCWW/k37+1225dgW3Kb4X897Tuwp0T6rsgqZyG+4oDnx1uCTcmyMcE33pmUYkaqsnnvuOUSEGTNmsGvXrkiHc0Zr3Mdfq3zChAk4HEX3ITudzrCMYOnVsBdfX/E1LWq1CKqOii3GxvgLxmO1+L9MJdgqiYFerwqrkIm8S5NaNDwr+J9zFhFu6ZtUdgGpSungwYNMmzYNh8OB0+nkqaeeimg8/lrjPgVb5evWrSMhIYHq1asHvNlsNsJVZTTWEstHf/iIVme1CtjVIgi2GBsPnv8gQ1sMDbitwc0HEx8TX+T+LGKhT6M+IcVclVTIRC4ivHVDNxLjim+VJ8TGcO+g82het1o5RKYqk+eeey6v9et0OiPeKvfXGvcp2Cr/7bffOHnyZLG39evXhy2+WvG1+HTYp4zrMY7G1RuTYE2gemx1qsdWJ84Sx8XNLmbaZdO4od0NRW5nZKuRCEX3ncZZ4hjdYXTYYq/swlKPXESmAsOBI8aYjsWtH2w98l/3pXPzlFXYXW6y7a4zlsVZLQjwwOBWjO3fUmcpUSVy8OBBWrRocUbr12q1ct111/Hxxx+XezxHjx6lefPmRU5eERMTw7Bhw5g1a1Y5RuafMYYd6TtIy0kjLiaOpJpJ1LbVDvr183fPZ/yP4/2OXLHF2Lij8x3c3vn2cIZcKQSqRx6uRH4RkAl8FM5EDpDjcDH314O8u2QHu49m43IbzkqM5doeTbmpd3Ma1koIMXpVFd11111MmTIFu/3M0RM2m43Nmzdz7rnnlms8TzzxBH//+9+xWAL/SDbG4HA42LZtG+edd145Rlc21h1Zx+trX2d96npiY2JxuB0k1Uzi7uS7ubjZxZEOr0Iq00Tu3UESMCfciVypcPPXGveJVKt8xYoVrFy5stj1LBYLN910E3Xq+B+jHY2O5Rzj2Klj1IirQYNqDSIdToUW8UQuIncAdwA0a9bs/D17dGiRioxArXGfSLXKlSpOoERebic7jTHvGWO6G2O616tXr7x2q9QZfCNVAiVxICwjWNxuN7/88ktI21AqWBVy1IpSZeWll17C4XBgtVoD3owxfPbZZ+zbt6/U+/n888/p3r0727dvD2P0SvmnI+5VlXLFFVfQuHHjYteLiYmhRo0apdqHy+Xi0UcfBTwnMadPn16q7SgVrLAkchH5NzAAOFtE9gFPG2NCn6ZEqTAbMGAAAwYMKNN9TJ8+nfT0dNxuN7Nnz2b79u2VYpSJqrjC0rVijLneGNPQGBNrjGmiSVxVVb7WeGamZ+Jip9PJE088EeGoVGWnfeRKhZGvNe7jdDrzWuVKlRVN5EqFScHWuI+2ylVZ00SuVJgUbI37aKtclTVN5EqFQaDWuI+2ylVZ0kSuVBgEao37aKtclSVN5EqFqLjWuI/dbtdWuSoTekGQUiXgcrmIiYkp9NyAAQPIyMgo9vVt27Ytq9BUFaaJXKkgpaSk0LVrV37++eczCmrFxcVFpIa5Uj7ataJUkCZOnMjx48d58sknIx1K5XLkN5h1N3w0An6eAs7ABc2Uf2ErY1sSWo9cRZuUlBRat25NTk6OlrkNp/1r4MPh4MwF44LYRGjaE276GnTWr0IiXsZWqWg2ceJEXC7PdINOp7NQq3zXrl3MmTMnEqFFt/kTwZHtSeLguZ+yCvavjWhY0UYTuVLFSElJ4bPPPsPhcACeRD5z5swzJmq+8847ueaaazh+/HikwoxOab8Xfk4E0raVfyxRTBO5UsXI3xr3yd8qX79+PUuWLAFg0qRJ5R5fVGvSE6RAGnK7oPH5kYknSmkiV6oIBVvjPvlb5Y8++ii5ubnk5OQwefJkbZWXxCXPgq0WWL2TqMcmQtdRcHaryMYVZXT4oVJF8Nca93E6ndxzzz0sWbIEt9sNeKZ4mzRpEi+++GJ5hhm96pwL9/0CG2bAiQNw3iBo3jfSUUUdHbWiKjVjDPv376dJkyYlfm3+kSqBiAgikpfIARITE9m3bx+1a9cuVcxKBaKjVlSV9O2333LuueeSkpJS4tcW1Rr3McackcThdKtcqfKiiVxVWsYYHn74YdxuNxMnTizx63/66SdiY2NJTEz0e7NY/P/30b5yVd40katK69tvv2Xfvn243W4+++yzErfKN27cSFZWlt/b8uXLiY+PD/habZWr8qSJXFVKvta4ryKhy+UqVas8EN9IlUC0VR7AyUMwZQhMOg8WPAsROEdXGYUlkYvIZSKyVUS2i8hj4dimUqHwtcZ9HA5HqVrl/mzZsoXvvvuOmJgY4uLiAt7sdjuTJ08OeX+Vyqy7Yd9qyEqFn96BbfMjHVGlEPLwQxGJAd4CLgH2AT+LyGxjzOZQt61UaRRsjfv4WuVTpkwJafv169dn8uTJhU5y+tOnT5+Q9lXpZKSAcXruGzecPBDZeCqJkIcfikhvYKIxZoj38eMAxpi/BnqNDj9UZWnu3Llcd911fid6sNls/P777zRt2jQCkSnW/RvmPggSA1Yb3LUCqtePdFRRoyyHHzYG8v9e3ed9rmAAd4jIahFZnZqaGobdqqrMGIO/Rkig1rhPuPvKVQklXw+3zYcr/gn3/KxJPEzK7WSnMeY9Y0x3Y0z3evXqldduVSX12GOPMWrUqELPF+wbLyicfeWqlM7pBO0vh8Q6kY6k0ghHIt8P5P+d2sT7nFJlIjU1lTfeeIMvv/yS338/XT2vuNa4j8Ph0Fa5qlTCUWvlZ6CViJyLJ4FfB9wQhu0q5deLL76I2+3G5XIxfvx4ZsyYAXgSdKNGjUhISCh2GwXn3VQqmoWl1oqIDAVeA2KAqcaYF4paX092qtJKTU2lefPmnDp1CvCcvFy/fj2tW7eOcGRKlb0yrbVijPnWGNPaGNOyuCSuVCh8rXEfp9PJ+PHjIxiRUpGnV3aqqJGamsq77757xhWVTqeTuXPnntFXrlRVo4lcRY2CrXEfbZWrqk4TuYoK/lrjPtoqV1WdJnIVFQK1xn20Va6qMk3kqsJLTU3lnXfeKbLaoNPpZM6cOdoqV1WSztmpKrzly5cDFFn/22fx4sU6FFFVOZrIVYU3YsSIvHHjSqnCtGtFKaWinCZypZSKcprIlVIqymkiV0qpKKeJXCmlopwmcqWUinKayJVSKsppIldKqSiniVwppaKcXtmplIITB+HX6XB8N9hqQbs/QuPzQSTSkakgaCJXqipz5sLse2HT157HrlxAYNV7UKspXP9vqNsykhGqIGjXilJVldsN/74eNs/2JHCXr7qkAUc2pG2DDwZBekpEw1TFCymRi8g1IrJJRNwiUmhCUKVUBbZ9AaSsBGeAgmTGDTknYMHEcg1LlVyoLfKNwFXAkjDEopQqL24XLHkZ7FlFr2dc8NscOHW8fOJSpRJSH7kxZguA6AkRpaJDxj5Y+gqs/zc4giwNHBMHR7ZA8z5lG5sqtXI72SkidwB3ADRr1qy8dquU8jmyBaYO8bTC3c6SvdaYsolJhUWxiVxEFgDn+Fn0hDFmVrA7Msa8B7wH0L17d/1UKFWe3G749GrIySj5a125UK9t+GNSYVNsIjfGDC6PQJRSZWjnIshOK/nrxAKthkC1uuGPSYWNDj9Uqir4/bvg+8TzCMRVh0ueKZOQVPiEOvzwShHZB/QG5orId+EJSykVVsd3l2z9uGpQqwncNg/qtCiTkFT4hDpq5SvgqzDFopQqK87c4tcBqH4OtPkDdLgSzr1IL9GPEnqJvlJVQZ1zYdfiYlYSGPm+J4GrqKJ95EpVBS0vBmt80euc1QySLiyfeFRYaSJXqrJLWQX714KliB/gFqunQJZ2pUSlCtO14nA42LdvHzk5OZEORVVxNpuNJk2aEBsbG+lQQrfxS5h1FzhyPEMJAaw2cHr/n4kVYqzwx9ehQYfIxalCUmES+b59+6hRowZJSUl6yb+KGGMMaWlp7Nu3j3PPPTfS4YTum/tPDzs0LkCg9rkQa/OcAD33QrhgrI5MiXIVJpHn5OSUKokbYzTxq7AREerWrUtqamqkQwlN5hHY9BXkniiwwHj6yu9YHImoVBmpMIkcgiu+ZYxh7d503l+yk8Vbj5DjdGOzWhjQpj539G9B16ZnaWJXIYnqz4/bDd89Dqunne5Kyc9qg+a9yz8uVaYqVCIvjt3p5sEv1rFwyxFynS7c3ootOU438zYf4offU7m4XX1e/VMycVY9j6uqoPkTYO1H+SaJKOCsZjDo6fKNSZW5qMl2xhge/GId3285zCnH6STu4zZwyuHi+y2HefCLdZhSVGvr0yf8ZTp3797NZ599FvbtKsX+NTDnQfjoCpjcFZ5vCCve8szuE4g9y9MqV5VK1CTyX1LSWbjlCDkOd5Hr5TjcLNxyhHUp6SXex/Lly0sZXWDRlsidzhKWN1XlL2M//LMvfDjc04WycxEc3wnOIhK4T046HNpQ5iGq8hU1ify9H3aS43QFtW6u08V7S3aWeB/Vq1cHYPHixQwYMICrr76atm3bcuONN+a18JOSknjkkUfo1KkTPXv2ZPv27QCMHj2aGTNmFNrWY489xtKlS0lOTubVV189Y3+LFy+mf//+jBgxghYtWvDYY4/x6aef0rNnTzp16sSOHTsASE1NZeTIkfTo0YMePXqwbNkyAFatWkXv3r3p2rUrffr0YevWrQBs2rSJnj17kpycTOfOndm2bRu7d++mY8eOeft++eWXmThxIgADBgzggQceoHv37kyePJk1a9bQv39/zj//fIYMGcLBgwdLfCxVGclMhfcHemqLO7KBEv7ylBhPMleVStT0kS/eeiTo2vZuA4u2Hglpf7/88gubNm2iUaNG9O3bl2XLltGvXz8AatWqxYYNG/joo4944IEHmDNnTsDt/O1vf+Pll18OuM769evZsmULderUoUWLFowZM4ZVq1YxefJk3njjDV577TXuv/9+HnzwQfr168fevXsZMmQIW7ZsoW3btixduhSr1cqCBQsYP348M2fO5J133uH+++/nxhtvxG6343K5OHz4cJHv1263s3r1ahwOB/3792fWrFnUq1eP6dOn88QTTzB16tTSH0wVPt8/4ylHa4Jr1BTisnv6yVWlEhWJ3BhDjrPoLpWCchzukIYm9uzZkyZNmgCQnJzM7t278xL59ddfn/fvgw8+WKrt+/To0YOGDRsC0LJlSy699FIAOnXqxKJFiwBYsGABmzdvznvNiRMnyMzMJCMjgz//+c9s27YNEcHhcADQu3dvXnjhBfbt28dVV11Fq1atio3j2muvBWDr1q1s3LiRSy65BACXy5UXn4qw3JOwYUbJZ/fJr0FHqJ0UtpBUxRAViVxEsFktJUrmtlhLSMPI4uNP16WIiYk5o+84/3Z9961WK263Jz63243dbi/xfiwWS95ji8WSt0+3283KlSux2c48SXXPPfcwcOBAvvrqK3bv3s2AAQMAuOGGG7jggguYO3cuQ4cO5d1336V169Z58QGFrqCtVq0a4PnS7NChAytWrAgqflWODqyDmNjAs94XSTylaS9/PdxRqQogavrIB7SpjyXIvGwRGNimfpnFMn369Lx/e/f2jMlNSkpizZo1AMyePTuvdVyjRg1OnjwZ0v4uvfRS3njjjbzH69atAyAjI4PGjRsD8OGHH+Yt37lzJy1atOC+++5jxIgR/PrrrzRo0IAjR46QlpZGbm5uwK6eNm3akJqampfIHQ4HmzZtCil+FSbBlqL1sVghvqbnAqBmveG2+XoZfiUVNYn8jv4tiLfGBLVuvDWGOy4qu0uOjx8/TufOnZk8eXLeCczbb7+dH374gS5durBixYq8Fm7nzp2JiYmhS5cuhU52Buv1119n9erVdO7cmfbt2/POO+8A8Mgjj/D444/TtWvXM34xfPHFF3Ts2JHk5GQ2btzIzTffTGxsLE899RQ9e/bkkksuoW1b/3MwxsXFMWPGDB599FG6dOlCcnJymYzmUaVwVjNwO4JbN646DJ4Io76E+9bBrf+FBu3LMjoVQVKa8dah6t69u1m9evUZz23ZsoV27doFfI0xhnv+/Qvfbzlc5BBEW6yFQe0a8Ob1XcvkCr2kpCRWr17N2WefHfZtq4qjuM9jxLzdyzNipTixiTBuB8Qlln1MqtyIyBpjTPeCz0dNi1xEePVPyQxq14CE2JhC3SwWgYTYGAa1a8Crf0qO7suslQrk4ichNqHodWIToef/aRKvQqLiZKdPnNXCm9d35ZcUT62VRVs9FwjZYi0MbFOfOy5qQXIZ11rZvXt3mW1bqWK1HQYXPQpLXvI/mXJsIrS6BAY9Wf6xqYiJqkQOnpZ5t2a1+eeo8wGtfqiqoAsfhKY9YOkrsGupZySLywH120K/B6HDVTpBRBUTUiIXkUnAHwE7sAO4xRiTHoa4ShJDee5OqYohqZ/nlpsJp45DfHVIqB3pqFSEhNpHPh/oaIzpDPwOPB56SEqpM5w4CP99BCadB39rDp9d55m6DTwJ/KymmsSruJBa5MaYefkergSuDi0cpdQZju2E9y/2tLx9Qw9//x/sWgxXvQ/t/hjR8FTFEM5RK7cC/w20UETuEJHVIrI65NlX7FnwwySY1AqeqeP594dJYA+i+lsRtIytqnC+eQByMgqMHzeeE51f/Z9nLk5V5RWbyEVkgYhs9HMbkW+dJwAn8Gmg7Rhj3jPGdDfGdK9Xr17pI7ZnwZRLYenLkHXEUzwo64jn8ZRLQkrmWsZWy9hWKFlpsHcFmEDXTQhs+65cQ1IVU7GJ3Bgz2BjT0c9tFoCIjAaGAzea8ri6aMVbkLb99CzgPs4cz/Mr3ir1prWMrZaxrVCy0yAmLvByt9MzN6eq8kIdtXIZ8AjQ3xgTWr9GsFa9XziJ+zhzYNV70H9cyLvRMrZaxjbiajYqutKhxQpnF1/ZUlV+oY4jfxOIB+Z7hwGuNMaMDTmqomSnhbY8SFrGVsvYRlx8deg40lO6ttAcnAK2WpB0UURCUxVLqKNWzgtXIEFLrOvpEy9qeRhoGVstY1sh/OHvcHijp9vQnuV5LjbRU9Fw1EywRE2VDVWGou9T0POOwJPHWm2e5WVMy9iqchNfHcYshCvf8ww1bDkYBj0N96+H+hWwqJeKiOhL5L3vhrrnFU7mVpvn+d53l3kIWsZWlasYK7QbDtd+AjfNhF5jPd0qSnlFTRnbM9izPaNTVr3n6RNPrOtpife+u8wrvmkZ26qhwpaxVVVaoDK2UVc0C/Ak6/7jwjI6RSmlol10JvII0jK2SqmKJvr6yJVSSp1BE7lS5c3tOj2UUKkwiOquFWMMOa4cbDE2rUuuosPKd2D+BE/9lLqtPGPBazWJdFQqykVlIt96bCsfbPiABXsX5M0QNLjZYMZ0GkObOm0iHZ5S/u38ARY87ZnNByD1d/jsWrhzWWTjUlEv6rpW5u6cy6hvRzFvzzycbicu48LpdjJvzzxGfTuKuTvnlnrbWsZWhV3WUU/C3r4QZt9boE6Q23PV5tHtEQtPVQ5Rlci3HtvKxOUTyXHl4C5Q2tNt3OS4cpi4fCJbj20t1fa1jK2WsQ2brf/1TAjxj3bw7oXwyZWQvsf/uu8PhCNbyjc+ValEVSL/YMMH2N1F1zCxu+1M2TilVNvXMrZaxjYsvhsPM26F/WvAZQ9crdMn9yR8chW4A9UdV6poUdNHboxhwd4FhVriBbmNmwV7Tvedl5aWsdUytqWy7t+weho4SlLV2UDOCdjxPbS6pMxCU5VX1CTyHFcOwZYT8HWzJFgTSr0/LWOrZWxLzBhY9GIJk7iXPRPWfqyJXJVK1CTyvCGGQeRyi1iwxQSokBgkLWOrZWxL7OC60OrhF1WeWakiRE0fuW+IoUWKDtkiFgY3H1ym48q1jK3y68SB0OqDV6sfvlhUlRI1iRxgTKcxxFmKmMMQiLPEcVvH28o0Di1jq/yy2oBSNiDiqkO3m8Iajqo6oq6M7dydc5m4fCJ2t/2ME58WsRBniWNin4kMazEs7DH7aBnbqqFUZWxzTsDLrcF5qoR7E6jZEB7YpDP+qCIFKmMbdZ+aYS2G8cnQTxiSNIRYSywxEkOsJZYhSUP4ZOgnZZrElSqSraZnjk1LbPCvEQvE14BRX2oSV6UW0slOEXkOGAG4gSPAaGPMgXAEVpQ2ddrw94v+HpFaK1rGVhXpkmc9wwgzU8EUcXFVXHVwO6HNUBj8NNROKrcQVeUT6qiVScaYJwFE5D7gKWBsyFEFSURCGmKoVNhVqwt3/ABf3OwZxeJ2em5YINYGNRvDhQ95CmU16ACJdSIdsaoEQkrkxpgT+R5WI6jBgUpVcjUawG3fQepW2PAfOHkQEs+GDldAo66Rjk5VQiGPIxeRF4CbgQxgYMgRKVVZ1GsDF0+IdBSqCij27IqILBCRjX5uIwCMMU8YY5oCnwL3FLGdO0RktYisTk1NDUvwBw8epGXLlhw6dCgs21NKqWhUbCI3xgw2xnT0c5tVYNVPgZFFbOc9Y0x3Y0z3evXqhRo3AM899xy7d+/mueeeC3lbKSkpDBw4kPbt29OhQwcmT55caB1jDGeffTbHjx8HPF8kIsKPP/6Yt069evVISwvh6r4wSE9P5+233w5q3ZKW7h0wYAAFh45WJOvWrePbb7+NdBhKlauQxjuJSP4iHiOA30ILJ3gHDx5k2rRpuN1upk2bFnKr3Gq18sorr7B582ZWrlzJW2+9dUZ9E/CcXO3Vq1feVY/Lly+na9eueRfMbN26lbp161K3bt2QYglWoJKzJUnk5XmxT3mUyC1NItfSvSrahTpw9W/ebpZfgUuB+8MQU1Cee+65vNohLpcr5FZ5w4YN6datG+C5rL5du3bs37+/0Hp9+vTJS37Lly/nwQcfPCOx9+3bl8zMTAYNGkS3bt3o1KkTs2Z5frxkZWUxbNgwunTpQseOHfMu9X/sscdo3749nTt35uGHHwYCl66dOHEiN910E3379uWmm27yW7L2scceY8eOHSQnJzNu3DgAJk2aRI8ePejcuTNPP/103vsJpnRvQR9//DHJycl07NiRVatWAYFL6n744YdcfvnlXHzxxQwaNIibb76Zr7/+Om9bN954I7NmzcLlcvHwww/TsWNHOnfunFeSIFBJ3QEDBvDoo4/Ss2dPWrduzdKlS7Hb7Tz11FNMnz6d5ORkpk+fTlZWFrfeeis9e/aka9eueX+LgnGVyOFN8PXd8Eo7mHQefDAYNn19euYfpcqbMabcb+eff74paPPmzYWeC+TAgQPGZrMZPKNkDGASEhLMwYMHg95GUXbt2mWaNm1qMjIyCi1bvHixGThwoDHGmH79+pmTJ08a3/sZM2aM+eCDD4zD4ch7bWpqqmnZsqVxu91mxowZZsyYMXnbSk9PN0ePHjWtW7c2brfbGGPM8ePHjTHGXH/99Wbp0qXGGGP27Nlj2rZta4wx5umnnzbdunUz2dnZxhhj7rnnHvPJJ58YY4zJzc012dnZZteuXaZDhw55+/nuu+/M7bffbtxut3G5XGbYsGHmhx9+MMYYU61aNWOMMYsWLTI1a9Y0KSkpxuVymV69euXtP7/+/fvnvYcffvghbz8ZGRnG4XAYY4yZP3++ueqqq4wxxkybNs00btzYpKWl5R2/ESNG5L3/pKQk43A4zNtvv21GjhyZt420tDRjt9tN7969zZEjR4wxxnz++efmlltuyYvjL3/5izHGmLlz55pBgwbl7e/uu+/Oi/fxxx83H3/8cd6xbdWqlcnMzCwUV0F+P49utzHfPmLMc/WNmVjbmKdrnr690MiY17oYk7Hf7/aUCgdgtfGTU6Om+mF++VvjPr5W+VtvvRXStjMzMxk5ciSvvfYaNWvWLLS8R48e/PLLL2RlZeFwOKhevTotWrRg+/btLF++nIceeghjDOPHj2fJkiVYLBb279/P4cOH6dSpEw899BCPPvoow4cP58ILL8TpdGKz2bjtttsYPnw4w4cPBwKXrgW4/PLLSUjwjJ8PpmTtvHnzmDdvHl27ds17j9u2beOiiy46Y72iSvfm5yvje9FFF3HixAnS09M5efKk35K6AJdccgl16njGS/fv35+77rqL1NRUZs6cyciRI/PqqY8dOxar1fORrFOnDhs3biyypO5VV10FwPnnnx/wQq158+Yxe/ZsXn75ZcBT9XHv3r2F4grKohdh7Uf+J4qwZ4LjFEy9DO5aCXGJwW9XqRBF3TXBvr7xgmVi7XZ7yH3lDoeDkSNHcuONN+YliZSUFJKTk0lOTuadd94hMTGRVq1aMXXq1LyumF69evHtt99y5MgR2rRpw6effkpqaipr1qxh3bp1NGjQgJycHFq3bs3atWvp1KkTEyZM4Nlnn8VqtbJq1Squvvpq5syZw2WXXQacLl27bt061q1bx/79+/O6QXwFucBTsnb27NkkJCQwdOhQFi5cWOh9GWN4/PHH87a1fft2brutcGGxokr35lfwKloR4cknn2TgwIFs3LiRb7755owyufnjBbj55pv55JNPmDZtGrfeemvAv4fxltT1xb1hwwbmzZtXKN6iYjXGMHPmzLxt7N27N6+GSsG4ipSTActfL7rWuHF55ujcOCPwOkqVgahL5P5a4z6h9JUbY7jtttto164df/nLX/Keb9q0aV4SGDvWc9Fqnz59eO211/JK2Pbu3ZvJkyfTq1cvRISMjAzq169PbGwsixYtYs8ez1yNBw4cIDExkVGjRjFu3DjWrl2bN0HE0KFDefXVV1m/fj0QuHRtQf5K1hYsnTtkyBCmTp2a16Lfv38/R46Uvva1r2//xx9/pFatWtSqVStgSV1/Ro8ezWuvvQZA+/btAU/r+N13381LyMeOHStVSV1/7/2NN97I6+//5ZdfSvZmfX79wlMXpTiOLFhWeMSTUmUpqhJ5oNa4Tyit8mXLlvHxxx+zcOHCvBZ4oNEPffv2ZefOnXmJvFu3buzbty9vKN+NN97I6tWr6dSpEx999FFeydgNGzbknZh85plnmDBhAidPnmT48OF07tyZfv368Y9//AMIXLq2IH8la+vWrUvfvn3p2LEj48aN49JLL+WGG26gd+/edOrUiauvvjqkGuk2m42uXbsyduxYpkzxzI8aqKSuPw0aNKBdu3bccsstec+NGTOGZs2a0blzZ7p06cJnn31WqpK6AwcOZPPmzXknO5988kkcDgedO3emQ4cOPPnkk6V70wfXBz/zz/EAkywrVUaiqoztXXfdxZQpU4qcfScuLo4xY8aE3Feuyk52djadOnVi7dq11KpVK9Lh+FXo8/jN/bDmw+BeHBMHT4bnojel8qsUZWxnz55d7BRqdrs9b4iZqngWLFhAu3btuPfeeytsEvereT9PxcJgnNOpbGNRqoCoGrWyb9++SIegQjR48OC8cwZRpf3lMCeIibZjq0G/0CbkVqqkKlSLPBLdPEoV5PdzaI2HYS9DbBFlk2PiPa3x1n8ou+CU8qPCJHKbzUZaWpomcxVRxhjS0tKw2WyFF3a5Dv4wCawJEJtvnLhYPfN1tugPN30JMVH1Q1dVAhXmE9ekSRP27dtHuCojKlVaNpst78KoQrrdBO1HwK/TYet/wZUL9dpBjzFQ3/+E1kqVtQozakUppVTRKsWoFaWUUoVpIldKqSiniVwppaJcRPrIRSQV2AOcDRwt9wAqLj0eZ9LjcZoeizNV1ePR3BhTaIq1iCTyvJ2LrPbXcV9V6fE4kx6P0/RYnEmPx5m0a0UppaKcJnKllIpykU7k70V4/xWNHo8z6fE4TY/FmfR45BPRPnKllFKhi3SLXCmlVIg0kSulVJSrEIlcRO4Vkd9EZJOI/D3S8USKiEwUkf0iss57GxrpmCoCEXlIRIyInB3pWCJJRJ4TkV+9n415ItIo0jFFkohM8uaNX0XkKxE5K9IxRUrEE7mIDARGAF2MMR2AlyMcUqS9aoxJ9t78TxpahYhIU+BSYG+kY6kAJhljOhtjkoE5wFMRjifS5gMdjTGdgd+BxyMcT8REPJEDdwJ/M8bkAhhjSj+9u6qMXgUeAar8WXljzIl8D6tRxY+JMWaeMcY30/dKIEDt4cqvIiTy1sCFIvKTiPwgIj0iHVCE3eP9qThVRGpHOphIEpERwH5jzPpIx1JRiMgLIpIC3Ii2yPO7FfhvpIOIlHIZfigiC4Bz/Cx6AngBWATcB/QApgMtTCUdF1nMsViJp36EAZ4DGhpjbi3H8MpdMcdjPHCpMSZDRHYD3Y0xlbq+RlHHwxgzK996jwM2Y8zT5RZcBARzPETkCaA7cFVlzRvFifg4chH5H/CSMWaR9/EOoJcxpkpPFSQiScAcY0zHSMcSCSLSCfgeyPY+1QQ4APQ0xhyKWGAVhIg0A76tqp8PHxEZDfwfMMgYk13M6pVWReha+RoYCCAirYE4qmZVM0SkYb6HVwIbIxVLpBljNhhj6htjkowxScA+oFtVTuIi0irfwxHAb5GKpSIQkcvwnD+5vConcagYLfI4YCqQDNiBh40xCyMaVISIyMd4joMBdgP/Z4w5GMmYKoqq0rVSFBGZCbQB3HjKQI81xuyPbFSRIyLbgXggzfvUSmPM2AiGFDERT+RKKaVCUxG6VpRSSoVAE7lSSkU5TeRKKRXlNJErpVSU00SulFJRThO5UkpFOU3kSikV5f4fUJ2/pQEsPT8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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_i * 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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}