wasserstein-julia/wass-bary-julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )\n",
    "\n",
    "# Implementing Wasserstein Barycenter problem using JuMP in Julia."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this notebook is to implement a model to calculate the Wasserstein barycenter in Julia with JuMP and then solve it using MosekTools. For additional info about the data used, theoretical explanation of the calculation of barycenters, references and for more insight in construction of the model, please consult the corresponding [Python notebook](https://nbviewer.jupyter.org/github/MOSEK/Tutorials/blob/master/wasserstein/wasserstein-bary.ipynb). Data files can be found in http://yann.lecun.com/exdb/mnist/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Plots\n",
    "pyplot()\n",
    "\n",
    "using JuMP\n",
    "\n",
    "using MosekTools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzkAAAGwCAYAAACKMqsgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X18VOW97/3v5IExYIiGmCcTYkSwllDqI4hsDEqjUVHqE6h3b9jbQ20FrAVqT+zxNHa7SQ89t9pbWmpbNw8q4D5uQdxYNCiEUsQKiEB0Y9CgoRBTLSQh4JCQdf7oJu1wDZCZrJU1ufJ5v17zatd3Vtb1m/h7JfyyJlcCjuM4AgAAAABLJPhdAAAAAAC4iSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYJcnvAk7U3t6uffv2KTU1VYFAwO9yECccx1Fzc7Nyc3OVkODtbE4P4kT0H/xGD8JP9B/8FGv/xd2Qs2/fPuXn5/tdBuJUXV2d8vLyPF2DHsTJ0H/wGz0IP9F/8FO0/Rd3Q05qaqokabRuUJKSfa4G8aJNrdqgVzv6w0v0IE5E/8Fv9CD8RP/BT7H2X9wNOcdvTSYpWUkBmhv/xfnr/3THrWt6EAb6D36jB+En+g9+irH/2HgAAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYJaohp6KiQpdffrlSU1OVmZmpCRMmaNeuXWHnFBcXKxAIhD0mTZrkatEAAAAAcDJRDTlVVVWaNm2aNm3apMrKSrW1tamkpEQtLS1h502dOlX79+/veDz99NOuFg0AAAAAJ5MUzcmrV68OO16wYIEyMzO1ZcsWjRkzpiPv27evsrOz3akQAAAAAKLQpd/JaWxslCSlp6eH5c8//7wyMjI0dOhQzZ49W83NzV1ZBgAAAAA6Lao7OX/PcRzNnDlTo0ePVlFRUUd+zz33qLCwUNnZ2dq5c6fKysr03nvvqbKyMuJ1QqGQQqFQx3FTU1OsJQExoQfhJ/oPfqMH4Sf6D16J+U7O9OnTtX37di1dujQsnzp1qsaNG6eioiJNmjRJL774otasWaOtW7dGvE5FRYXS0tI6Hvn5+bGWBMSEHoSf6D/4jR6En+g/eCWmIWfGjBlauXKl1q5dq7y8vFOee8kllyg5OVk1NTURny8rK1NjY2PHo66uLpaSgJjRg/AT/Qe/0YPwE/0Hr0T1djXHcTRjxgwtX75c69atU2Fh4Wk/prq6Wq2trcrJyYn4fDAYVDAYjKYMwFX0IPxE/8Fv9CD8RP/BK1ENOdOmTdOSJUv08ssvKzU1VfX19ZKktLQ0paSk6KOPPtLzzz+vG264QRkZGXr//fc1a9YsXXzxxbrqqqs8eQEAAAAA8Peierva/Pnz1djYqOLiYuXk5HQ8XnjhBUlSnz599MYbb+i6667ThRdeqAceeEAlJSVas2aNEhMTPXkBAAAAAPD3on672qnk5+erqqqqSwUBAAAAQFfEvIU0AHu1X32Jkf3zwt8Y2QM/mW5kZy94y5OaAAAAOqtLfwwUAAAAAOINQw4AAAAAqzDkAAAAALAKQw4AAAAAqzDkAAAAALAKu6t1o2PXXmpkP/z1YiN7/M6JRta+ZacnNQGJ6Wcb2UP/avblFUHzZyJ/KTKvZ14N6LqkzHOMLPuVI0a26ZWvGVneYxs9qQk4UdKAAUZWsLrFyH557ttGVvHFECNb89mFRub8LMvIkl97p7MlAkrKO9fIMl9sNrLCvp8b2b8tKzayc+fE59dY7uQAAAAAsApDDgAAAACrMOQAAAAAsApDDgAAAACrMOQAAAAAsApDDgAAAACrsIV0d2o3o2tTjhrZkl98amT7RnpRECAFgkEjG5vS6kMlwMm9/1ihkb2S/7SRDf26eR7QXRommFs+v3LuL4wswj8H9Mz2UUY2eOou88RjDZ26HnAyDU/3M7JXBv6Hka09kmxkG+eaH+u4U5bruJMDAAAAwCoMOQAAAACswpADAAAAwCoMOQAAAACswpADAAAAwCrsrtaNgh+bO6L84uB5RjbxnLeN7OfDbzOy9vfed6UuAIh3FcUvGtnrR84wsvO/u9/I2jypCDAdm/CXmD928NyQkbUfPtyVcgAlDcwzslvydxjZR61mr/3k+983sjPa/uhOYd2AOzkAAAAArMKQAwAAAMAqUQ05FRUVuvzyy5WamqrMzExNmDBBu3aF/6GqUCikGTNmKCMjQ/369dPNN9+svXv3ulo0AAAAAJxMVENOVVWVpk2bpk2bNqmyslJtbW0qKSlRS0tLxzkPPvigli9frmXLlmnDhg06dOiQbrrpJh07dsz14gEAAADgRFFtPLB69eqw4wULFigzM1NbtmzRmDFj1NjYqGeeeUbPPvusxo0bJ0l67rnnlJ+frzVr1ui6665zr3IAAAAAiKBLu6s1NjZKktLT0yVJW7ZsUWtrq0pKSjrOyc3NVVFRkTZu3BhxyAmFQgqF/rajSFNTU1dKAqJGD8JP9B/8Rg/CT/QfvBLzkOM4jmbOnKnRo0erqKhIklRfX68+ffro7LPPDjs3KytL9fX1Ea9TUVGhRx99NNYyepS2T+qMbHHtCCN7++IXjOz+b6ca2eBp7tTV2/WmHkT8of9Mjd+60si+2W+zkQ1f8ICRFfx5oyc12YwejA8Jez8zsnYf6uhu9J+7kgryw46HvGRuq/9whvknSMZ/eIeRHTzfHBOyu1Bbd4t5d7Xp06dr+/btWrp06WnPdRxHgUAg4nNlZWVqbGzseNTVmYMA4CV6EH6i/+A3ehB+ov/glZju5MyYMUMrV67U+vXrlZf3tz8ylJ2draNHj+rAgQNhd3MaGho0atSoiNcKBoMKBoOxlAG4gh6En+g/+I0ehJ/oP3glqjs5juNo+vTpeumll/Tmm2+qsLAw7PlLL71UycnJqqys7Mj279+vnTt3nnTIAQAAAAA3RXUnZ9q0aVqyZIlefvllpaamdvyeTVpamlJSUpSWlqZ7771Xs2bN0oABA5Senq7Zs2dr2LBhHbutAQAAAICXohpy5s+fL0kqLi4OyxcsWKApU6ZIkp544gklJSXpzjvv1JEjR3Tttddq4cKFSkxMdKVgAAAAADiVqIYcx3FOe84ZZ5yhp556Sk899VTMRfUmZ85LM8NnzOjp6//VyJ742u1Gdmz7B26UBQC+cSL8TCxB5uY1gbZuKAY4iYS+fY1sQN8WIzvmmHukfbXqvxnZBc3V7hSGXu2Tu8N3V3slZ6VxzvPNmUbmfPdMI2ua2WpkvWJ3NQAAAACIRww5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKwS1RbScF/S4WNG1nDssJFdm2J+bEWOud1f8nZXygIA35z14REj+6gtZGSBoubuKAeIqO2KrxjZaxeZfwPiQLvZu+ff/a6RmRtNA6e27wejjOzN+392QmJudT73mTuMLPeDjUY2ZGrMpcUF7uQAAAAAsApDDgAAAACrMOQAAAAAsApDDgAAAACrMOQAAAAAsAq7q/ksoWqrkX3vk1uMbOn5ld1RDgD47uAQczvJQUlBI3N2pnZHOUBEn1x/ht8loJf7+q3vG1lGYvhuaqPfM3dSG/jMLiNrc6+suMGdHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBW2kAYAADiFljtGGln1t56KcKb5s+P799wc4bzPu14UepX6B0cZ2W/OnWtkH7Qmhh2f8cRZxjltX+x2r7A4xp0cAAAAAFaJeshZv369xo8fr9zcXAUCAa1YsSLs+SlTpigQCIQ9Ro40fwICAAAAAF6IeshpaWnR8OHDNW/evJOec/3112v//v0dj1dffbVLRQIAAABAZ0X9OzmlpaUqLS095TnBYFDZ2dkxFwUAAAAAsfLkd3LWrVunzMxMDRkyRFOnTlVDQ4MXywAAAACAwfXd1UpLS3XHHXeooKBAtbW1euSRR3TNNddoy5YtCgaDxvmhUEihUKjjuKmpye2SepwPVgwxw5mVRrRnomNkg1/zoiK70YPwE/1naioM+F1Cr0IPnt6fLzZ/JpwYMLP9x1qMrO4XFxhZKrurdaD/TEdvvMLIXplp7qSWl9TPyB7+7OKw4+TXN7tXWA/j+p2ciRMn6sYbb1RRUZHGjx+v3/3ud/rwww+1atWqiOdXVFQoLS2t45Gfn+92ScAp0YPwE/0Hv9GD8BP9B694voV0Tk6OCgoKVFNTE/H5srIyNTY2djzq6uq8LgkIQw/CT/Qf/EYPwk/0H7zi+R8D/eKLL1RXV6ecnJyIzweDwYhvYwO6Cz0IP9F/8Bs9CD/Rf/BK1EPOoUOHtHv33/5Sam1trbZt26b09HSlp6ervLxct912m3JycrRnzx49/PDDysjI0De/+U1XCwcAAACASKIecjZv3qyxY8d2HM+cOVOSNHnyZM2fP187duzQ4sWLdfDgQeXk5Gjs2LF64YUXlJqa6l7VAAAAAHASUQ85xcXFchxzV6/jXnuN7b0AAAAA+Mfz38lB9M7afaxT5w0+r97jSoBTu2DpISM7+Y9AAFNS3rlG9oOJL/lQCfBXSYMKjex/3vZCpz72obrxRpa6dFOXa0Lv0nCx+c/zgUnmO6La1W5kq5aOCjvO1Ub3CuthPN9dDQAAAAC6E0MOAAAAAKsw5AAAAACwCkMOAAAAAKsw5AAAAACwCrur9WA1n2YZ2WDt9aES9FZNQ840stTNPhSCHqv1vHOM7J/6f2ZkNa1fmh8c8KIi9HYZz/3FyO5K/bxTH/vB8xcZ2Tm9eHcrnF7SublGNuOelUYWaSe1ob+ZZmTnP1cbdtzWhdp6Ou7kAAAAALAKQw4AAAAAqzDkAAAAALAKQw4AAAAAqzDkAAAAALAKu6vFoZQVb+vIhBGnPe+8pcyo6Lq2AnOXvs46lGtub5XalWIASa2OuR/QH78caJ7odEMxsF5i+tlhx3++USrbvO60H7e8Jc3IsteZu7Adi7ky9AZN/3qGkX07rc7IatsOG1nu748aWdu+/e4UZgH+lRyHOjPgAACArjlxwJE6N+AAiH8MOQAAAACswpADAAAAwCoMOQAAAACswpADAAAAwCoMOQAAAACswhbSPdgXFyUbWfZrPhSCHu2jO/vG/LH5z+42MnPzX+Dkkvc3GtnqI/2N7M4zG4xsyW/2GBn9h2h8/P2LjOyqM940slbH3Aj6hyvvMbJB1W+5UxisFBp/hZEtvOiJCGea35dv33avkWWu2exGWdbiTg4AAAAAqzDkAAAAALBK1EPO+vXrNX78eOXm5ioQCGjFihVhzzuOo/LycuXm5iolJUXFxcWqrq52rWAAAAAAOJWoh5yWlhYNHz5c8+bNi/j83Llz9fjjj2vevHl65513lJ2drW984xtqbm7ucrEAAAAAcDpRbzxQWlqq0tLSiM85jqMnn3xSP/rRj3TrrbdKkhYtWqSsrCwtWbJE9913X9eqBQAAAIDTcHV3tdraWtXX16ukpKQjCwaDuvrqq7Vx48aIQ04oFFIoFOo4bmpqcrMkq2Xf/KkZPtn9dfR0vb0HL3yyzgwnde5ja++7wMjyf2LugoWT6+3915qTZmTXp5ifg4b2I0bW9qd9ntTU2/SmHjx864iw43f/KdI3TXPn0rs/vsHIBs1iJzU39Kb+Ozb9cyMrTDJ3Unvs86FGlnu/+Y4odpM8NVc3Hqivr5ckZWVlheVZWVkdz52ooqJCaWlpHY/8/Hw3SwJOix6En+g/+K239OCJAw7iQ2/pP3Q/T3ZXCwQCYceO4xjZcWVlZWpsbOx41NVF+Kky4CF6EH6i/+A3ehB+ov/gFVffrpadnS3pr3d0cnJyOvKGhgbj7s5xwWBQwWDQzTKAqNCD8BP9B7/Rg/AT/QevuHonp7CwUNnZ2aqsrOzIjh49qqqqKo0aNcrNpQAAAAAgoqjv5Bw6dEi7d+/uOK6trdW2bduUnp6ugQMH6sEHH9ScOXM0ePBgDR48WHPmzFHfvn119913u1o4AAAAAEQS9ZCzefNmjR07tuN45syZkqTJkydr4cKFeuihh3TkyBHdf//9OnDggEaMGKHXX39dqamp7lUNAAAAACcR9ZBTXFwsx3FO+nwgEFB5ebnKy8u7UheAbrLvloFGlhDhnaztajeywqd3GxlbWsILY597yMgKtdGHStBT7Ss2N0AKBsztohc3ZRvZ4ekZEa7Idvk4udD4K4xs4UVPRDjT3EJ6+Z6vGVnm3v90o6xexZPd1QAAAADALww5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKlHvrgbALs2F5m6JkXZSA/x0wZK/GNkxH+pAz3G09PKw45W3RNrZ6gwjWfLfSo0s8N42t8pCL/Hnr5n/xC5MMndS+/zYYSM752cpntTU23AnBwAAWOXEAQdA78OQAwAAAMAqDDkAAAAArMKQAwAAAMAqDDkAAAAArMKQAwAAAMAqbCHdQyREmEcTApG2/gWis2jCLzt13t62FjNsp+PQPT67Kt3IMnb4UAjiUtKFF4Qff3xAj6x+ISy7KNncLvr1I2aWvOfPRtbWxfpgt8QzzzSyn/7jok597Lif/8DIcn6/scs1gTs5AADAMicOOAB6H4YcAAAAAFZhyAEAAABgFYYcAAAAAFZhyAEAAABgFXZXi0OfXZFoZO0R9k2r2TrQyAbpT57UBHv1CRyLkJo9eOOW+4ws98/VHlQEmJouMHeTzPChDsSnD6aHd8Pd/zFNu2+bf9qPe3zK3UYW2Puua3Whd2j+90wju7HvISM73H7UyHL+EGHnUriCOzkAAMAqnRlwANiNIQcAAACAVVwfcsrLyxUIBMIe2dnZbi8DAAAAABF58js5Q4cO1Zo1azqOExPN9/cDAAAAgBc8GXKSkpK4ewMAAADAF54MOTU1NcrNzVUwGNSIESM0Z84cnX/++RHPDYVCCoVCHcdNTU1elAScFD0IP9F/8Bs9CD/Rf/CK60POiBEjtHjxYg0ZMkSfffaZHnvsMY0aNUrV1dUaMGCAcX5FRYUeffRRt8vo0TLeM7eLjmTg620eV9I70IOd01qd5ncJVqL/OmfwogNGFmnzc0Svp/Vg7f+60sgev/7ZsOOVLam6uV9zWLa3zdyqN3HzfxpZ574Dwy09rf+Oll5uZM9f9ISRPfrnkUa2+ol/MLKz3nrLncJgcH3jgdLSUt12220aNmyYxo0bp1WrVkmSFi1aFPH8srIyNTY2djzq6urcLgk4JXoQfqL/4Dcbe/DEAQfxy8b+Q3zw/I+B9uvXT8OGDVNNTU3E54PBoILBoNdlACdFD8JP9B/8Rg/CT/QfvOL538kJhUL64IMPlJOT4/VSAAAAAOD+kDN79mxVVVWptrZWb7/9tm6//XY1NTVp8uTJbi8FAAAAAAbX3662d+9e3XXXXfr88891zjnnaOTIkdq0aZMKCgrcXgoAAAAADK4POcuWLXP7kr1OWrW5i9DKllQfKkFv8D8KL+vUeQXa6HElwF899vlwv0tAPGsPGFGkjQYOtB8JO5583/eNc/ocece9utArNFySbGR5Sf2MbMXCMUaWvZDvo93J89/JAQAA6E4nDjgAeh+GHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWcX0LaXTdsZ27jGz+4AuMLFlsfQmgZwts2GZkm4ZH+tb0n94Xgx6hsMzchre07OLTflwfvmfCBXn/EqH//sXsv2z+7ILvuJMDAAAAwCoMOQAAAACswpADAAAAwCoMOQAAAACswpADAAAAwCoMOQAAAACswpADAAAAwCoMOQAAAACswpADAAAAwCoMOQAAAACswpADAAAAwCpJfhdwIsdxJEltapUcn4tB3GhTq6S/9YeX6EGciP6D3+hB+In+g59i7b+4G3Kam5slSRv0qs+VIB41NzcrLS3N8zUkehAm+g9+owfhJ/oPfoq2/wJOd4zlUWhvb9e+ffuUmpqqQCDQ5es1NTUpPz9fdXV16t+/vwsVsoYf13ccR83NzcrNzVVCgrfvsuxpPWhDb8T7GvSfv2vY8Bq6ukZP7cF4/7z2pjV6Y/9J8f15ZY3OibX/4u5OTkJCgvLy8ly/bv/+/T37j8Ya3XN9r396dFxP7UEbeiOe16D//F/DhtfQlTV6cg/G8+e1t63RG/tPit/PK2t0Tiz9x8YDAAAAAKzCkAMAAADAKonl5eXlfhfhtcTERBUXFyspybt357FGfFw/XtnweWWNnsuGz6sNr6G71og3tnxebVijN/afZMfnlTWiF3cbDwAAAABAV/B2NQAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWSfK7gBO1t7dr3759Sk1NVSAQ8LscxAnHcdTc3Kzc3FwlJHg7m9ODOBH9B7/Rg/AT/Qc/xdp/cTfk7Nu3T/n5+X6XgThVV1envLw8T9egB3Ey9B/8Rg/CT/Qf/BRt/8XdkJOamipJGq0blKRkn6tBvGhTqzbo1Y7+8BI9iBPRf/AbPQg/0X/wU6z9F3dDzvFbk0lKVlKA5sZ/cf76P91x65oehIH+g9/oQfiJ/oOfYuw/Nh4AAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWiWrIqaio0OWXX67U1FRlZmZqwoQJ2rVrV9g5xcXFCgQCYY9Jkya5WjQAAAAAnExUQ05VVZWmTZumTZs2qbKyUm1tbSopKVFLS0vYeVOnTtX+/fs7Hk8//bSrRQMAAADAySRFc/Lq1avDjhcsWKDMzExt2bJFY8aM6cj79u2r7OxsdyoEAAAAgCh06XdyGhsbJUnp6elh+fPPP6+MjAwNHTpUs2fPVnNzc1eWAQAAAIBOi+pOzt9zHEczZ87U6NGjVVRU1JHfc889KiwsVHZ2tnbu3KmysjK99957qqysjHidUCikUCjUcdzU1BRrSZBU+9NRRnbpP+wyssYbjxnZsYMHPakp3tGDpsSLhhjZ15Z8aGQ/zXrPyAYv+U7Y8fmz33KvMAvRf/AbPQg/0X/wSsx3cqZPn67t27dr6dKlYfnUqVM1btw4FRUVadKkSXrxxRe1Zs0abd26NeJ1KioqlJaW1vHIz8+PtSQgJvQg/ET/wW/0IPxE/8ErMQ05M2bM0MqVK7V27Vrl5eWd8txLLrlEycnJqqmpifh8WVmZGhsbOx51dXWxlATEjB6En+g/+I0ehJ/oP3glqrerOY6jGTNmaPny5Vq3bp0KCwtP+zHV1dVqbW1VTk5OxOeDwaCCwWA0ZQCuogfhJ/oPfqMH4Sf6D16JasiZNm2alixZopdfflmpqamqr6+XJKWlpSklJUUfffSRnn/+ed1www3KyMjQ+++/r1mzZuniiy/WVVdd5ckLAAAAAIC/F9Xb1ebPn6/GxkYVFxcrJyen4/HCCy9Ikvr06aM33nhD1113nS688EI98MADKikp0Zo1a5SYmOjJCwAAAACAvxf129VOJT8/X1VVVV0qCAAAAAC6IuYtpOG/j392pZH9ftLPjCwzsa+R3Vhwj3nBXrqFNEwNozOM7J8zlxpZa4Sfe/xh4v8OOx5X/wPjnJz/vTH24tArJQ77ipHd9eIbRjYqZY+R3fePD5jXe2OLK3XBPu1jLjayf170WyO7Z8NUIzv/N+YXxYT177pTGICodOmPgQIAAABAvGHIAQAAAGAVhhwAAAAAVmHIAQAAAGAVhhwAAAAAVmF3tR6ibc1AI3v/onlGlhgwd1IDutPZCWeEHR8a3OZTJbBJ49CzjOye1AYjW3skzcj6/PFDIzvmTlmwUO23zeziPuauae9f82sj23yV+TcBf3K+uVsbEG9qfnWFkT017lkj++XesUbWVrzPk5q6ijs5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKmwh7bM//fdRRnbDxE1GNjfrFSP7+twZRnZ0dLOR7Ry1OMbqACA+/OUrnfuZ3Ihgi5Ht+6dhRpb1841drgl2+mDsb42svZMf287PjhFvAmZP/qlspJFtu+lxI3uh+Xwjq192npFliC2kAQAAAMBzDDkAAAAArMKQAwAAAMAqDDkAAAAArMKQAwAAAMAq7K7mkqSC/LDjIxdmGefsudWcKX9/41wjy0nsZ2QXLLvfzJ5628hqLh1+yjoBIN4l5Z1rZD+c9O+d+ti+CX2M7OhZXS4J6JT7FnzXyPLFTn7wz96HzZ3Udkz7hZH94uCFRrb6lsuMLKOm5/Qzd3IAAAAAWIUhBwAAAIBVohpyKioqdPnllys1NVWZmZmaMGGCdu3aFXZOKBTSjBkzlJGRoX79+unmm2/W3r17XS0aAAAAAE4mqiGnqqpK06ZN06ZNm1RZWam2tjaVlJSopeVvf2H6wQcf1PLly7Vs2TJt2LBBhw4d0k033aRjx465XjwAAAAAnCiqjQdWr14ddrxgwQJlZmZqy5YtGjNmjBobG/XMM8/o2Wef1bhx4yRJzz33nPLz87VmzRpdd9117lUOAAAAABF0aXe1xsZGSVJ6erokacuWLWptbVVJSUnHObm5uSoqKtLGjRsjDjmhUEihUKjjuKmpqSslAVGjB+En+g9+owfhJ/oPXol5yHEcRzNnztTo0aNVVFQkSaqvr1efPn109tkUj+2KAAAY8klEQVRnh52blZWl+vr6iNepqKjQo48+GmsZnks880wjS119hpFNz/2PsOOrzmjv1PVXHTa3mq545FtGNmjpW0bmdGoF6Y0j5paqCX8+YGSdq9g+8d6DsBv9Z6r9pwIj+3/7/0eEM03bjx41soKVB42st369i4QejM2CpvOMLP+fe872uvGC/uuCgPlbJ3X/I3zL6E33PW6cE3G76FuvMLK2mpouFOe/mHdXmz59urZv366lS5ee9lzHcRQIBCI+V1ZWpsbGxo5HXV1drCUBMaEH4Sf6D36jB+En+g9eielOzowZM7Ry5UqtX79eeXl5HXl2draOHj2qAwcOhN3NaWho0KhRoyJeKxgMKhgMxlIG4Ap6EH6i/+A3ehB+ov/glaju5DiOo+nTp+ull17Sm2++qcLCwrDnL730UiUnJ6uysrIj279/v3bu3HnSIQcAAAAA3BTVnZxp06ZpyZIlevnll5WamtrxezZpaWlKSUlRWlqa7r33Xs2aNUsDBgxQenq6Zs+erWHDhnXstgYAAAAAXopqyJk/f74kqbi4OCxfsGCBpkyZIkl64oknlJSUpDvvvFNHjhzRtddeq4ULFyoxMdGVggEAAADgVKIachzn9Pt5nXHGGXrqqaf01FNPxVxUXEkyP0VDzvzMyFITvgw7vqf2BuOcPb8cYmTplbvNazVsiqbC08pObDbDM3j/K4D49PPJv4n5Y+uP9Tey9nff70o5sNhH/9+VEdIt3V4HEIuEr33FyHZ+9xdhx786ONg4Z/UtlxlZT99JLZKYd1cDAAAAgHjEkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKlFtId0bHTt40Mj++HXzb/5sHXhb2HHbp3uNc/rrLSNr60JtnfXsgVHmuh/v6YaVAQCIDwkpKUYWPC/Cn1jopG3NAyOkh2O+HnAqSRecb2Tf+j+rjexA+5Gw4xXfvtY4J1Czzb3C4hh3cgAAAABYhSEHAAAAgFUYcgAAAABYhSEHAAAAgFUYcgAAAABYhd3VXBJpNzUAQPcq+/m9RpapjT5Ugrhzkbk71daRCyOcGOjU5d78eLCRFeq9KIsCOuejOWca2Z1nfmFkX9s0Nez43A29Yye1SLiTAwAAAMAqDDkAAAAArMKQAwAAAMAqDDkAAAAArMKQAwAAAMAqDDkAAAAArMIW0r3A/9l8qZEN0Ts+VIKe4sVH5kZIU7q9DiBagWN+V4CeJCHCdtHJAfOfRq1OW3eUA0iSDk0caWQ7rvqFkW0JOUaWf8/HYcft7pXV43AnBwAAAIBVoh5y1q9fr/Hjxys3N1eBQEArVqwIe37KlCkKBAJhj5EjzYkUAAAAALwQ9ZDT0tKi4cOHa968eSc95/rrr9f+/fs7Hq+++mqXigQAAACAzor6d3JKS0tVWlp6ynOCwaCys7NjLgoAAAAAYuXJ7+SsW7dOmZmZGjJkiKZOnaqGhgYvlgEAAAAAg+u7q5WWluqOO+5QQUGBamtr9cgjj+iaa67Rli1bFAwGjfNDoZBCoVDHcVNTk9sl9XrB/cl+lxDXelMPBpLMXvj0R5cbWVbiH7ujHKh39V8kDTNGGdmVwbcjnNnH+2J6qd7cg+0yd6eKtJNapPMyXmbHSTf05v6TpKSsTCO799EVRlbdavblD78z3ciSj7B77nGu38mZOHGibrzxRhUVFWn8+PH63e9+pw8//FCrVq2KeH5FRYXS0tI6Hvn5+W6XBJwSPQg/0X/wGz0IP9F/8IrnW0jn5OSooKBANTU1EZ8vKytTY2Njx6Ours7rkoAw9CD8RP/Bb/Qg/ET/wSue/zHQL774QnV1dcrJyYn4fDAYjPg2NqC70IPwE/0Hv9GD8BP9B69EPeQcOnRIu3fv7jiura3Vtm3blJ6ervT0dJWXl+u2225TTk6O9uzZo4cfflgZGRn65je/6WrhAAAAABBJ1EPO5s2bNXbs2I7jmTNnSpImT56s+fPna8eOHVq8eLEOHjyonJwcjR07Vi+88IJSU1PdqxoAAAAATiLqIae4uFiOY+4yctxrr73WpYIAAAAAoCs8/50cAPGj5ZZLjey9bz8V4cxE74sBJH05wMz6JrBdNIDe4YOfnGdkU/qbNwyGLHvQyAa99pYXJVnD893VAAAAAKA7MeQAAAAAsApDDgAAAACrMOQAAAAAsApDDgAAAACrsLsaAFdtOxq+xfyFvzxknNPeXcXAGgkRfibXTicB6EECiebOpUMvrDOyXzfmG9kFD202spP/QRdI3MkBAAAAYBmGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBV2V7PM2We3GNmh2jN8qATxaO/1xzxfIzXhaNjxEyt/q29/7/vhJxWOUMqKtz2vBfZgJzV0VfvWaiVcMjQsmzB+sla+stg4NzkQ/s+jVqfNvGAg4Gp9sF/NguFmNuQZI7tgxXeMbHAb3zOjxZ0cAJ4yBhyJAQdAtztxwJHUqQEHQM/EkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKmwhYpkFRYuM7I6tM32oBH5LOq/AyP5t3PwIZ7q7DeqgpGDY8d5Sxzhn8ApXlwSAmLTL/PoUabvoSOfJiZAB/yWQmGhkIy+oNbJfN+Yb2VfK/tPIvP8DEPbhTg4AAAAAqzDkAAAAALBK1EPO+vXrNX78eOXm5ioQCGjFivD3nTiOo/LycuXm5iolJUXFxcWqrq52rWAAAAAAOJWoh5yWlhYNHz5c8+bNi/j83Llz9fjjj2vevHl65513lJ2drW984xtqbm7ucrEAAAAAcDpRbzxQWlqq0tLSiM85jqMnn3xSP/rRj3TrrbdKkhYtWqSsrCwtWbJE9913X9eqBQAAAIDTcHV3tdraWtXX16ukpKQjCwaDuvrqq7Vx48aIQ04oFFIoFOo4bmpqcrOkXmdP2wC/S+hxbO3B9x/OMrKv93F3J7VI9rR9GXb82vVP6B9nzfJ83Z7K1v5Dz9FberB9a7USzzorLLvu/duN8yq/+u+dut6Rc8w3w6TGVlqvZmv/7Z57hZG9et4vjeyrz0wzsoLGjZ7U1Nu4uvFAfX29JCkrK/wfV1lZWR3PnaiiokJpaWkdj/x8cys9wEv0oLcYcE6N/oPfeksPnjjgID70lv5D9/Nkd7VAIPynxY7jGNlxZWVlamxs7HjU1dV5URJwUvQg/ET/wW/0IPxE/8Errr5dLTs7W9Jf7+jk5OR05A0NDcbdneOCwaCCwWDE54DuQA/CT/Qf/EYPwk/0H7zi6p2cwsJCZWdnq7KysiM7evSoqqqqNGrUKDeXAgAAAICIor6Tc+jQIe3evbvjuLa2Vtu2bVN6eroGDhyoBx98UHPmzNHgwYM1ePBgzZkzR3379tXdd9/tauEAAAAAEEnUQ87mzZs1duzYjuOZM2dKkiZPnqyFCxfqoYce0pEjR3T//ffrwIEDGjFihF5//XWlprLnCAAAAADvRT3kFBcXy3Gckz4fCARUXl6u8vLyrtQV15LOKzDDlpaww7Y/f95N1YR78K2JRlZ09cdG9qWRoEcb9XUj+t11T0Y40d33Pf8hlGxkP/no1rDjfi9ucnVN4FR+OfspI7uvzwwjy36CLVp7m2MHDxrZp3svME/8aueud993VxrZ8v8/I9qyYKnffvNpI/v82GEjG7Rwv5G1eVJR7+PJ7moAAAAA4BeGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYJWod1ezRcLF5vYptQ+bn46rBpo7k/1L7hIj2/RllpG98hdzx6sTrf1j0WnPkaTMdwJG1pxvZnNHPmdkXw+aO3fcMOcHRnbBb/cZWdvHezpVH/x1LJhoZIOSvP8L0t/+47eMrHDSe56vC5zMFUHzZ3et/X0oBD3C4H/cbGTJ+8x/C7Q65n5X9/b/xMgWvnqlkaXdUBNjdegpks4/z8jOSXzLyEa++n0jG/zRH70oCeJOjmc6M+AAAAAAcB9DDgAAAACrMOQAAAAAsApDDgAAAACrMOQAAAAAsApDDgAAAACr9NotpPf8yNxy9/1Rizr50X2N5OZ+zWHH41LeMM5JDphrJuf/vnNL3ta50yIz6/1gyi+MbOltGUb24meXGVnD4TONLPWmT43MaWvtbIHoAb6/f5SRXXDfHiM75n0psEjwgJkdbj9qZH0T+nRDNYB0wb/dZ2TVd8zr1MdWDV9qZDePvNc8cRNb7fcEbePMfwPVjzC/FoXOaTeyocnmv71y1nJvoTvx2QYAAABgFYYcAAAAAFZhyAEAAABgFYYcAAAAAFZhyAEAAABglV67u9pZK8wdwq6onGZkCUedmK5/9q4jRhYaYO7I0ZJl/ic4nGVer9++2Oo4mUP5ASNr62eu8W8TnzSyr/UxX8c/3P5dIztz2aYYq4PfHqofYWS1/89AIzvWWNMd5cBiWT/faGRvPWB+fb42xdxxLZLqoyEjy3+tJfrC0Gv1+Yu7P//96A5zl61BfHvsEZLWbDayvDXmeaHK84zsD1+aO66dtf4TI2uLqTJ0BndyAAAAAFiFIQcAAACAVVwfcsrLyxUIBMIe2dnZbi8DAAAAABF58js5Q4cO1Zo1f3vTYmJiohfLAAAAAIDBkyEnKSmJuzcAAAAAfOHJkFNTU6Pc3FwFg0GNGDFCc+bM0fnnnx/x3FAopFDob7vhNDU1eVEScFL0IPxE/8Fv9CD8RP/BK64POSNGjNDixYs1ZMgQffbZZ3rsscc0atQoVVdXa8CAAcb5FRUVevTRR90u47T6P/dWt68Z7GSW7nUhks7u5Hk//KG5lXAkZ6rn7ofpVw+6KXHtFiO76dxLu3DFSJtasl20F2zoP7d9/5mpRvbze39tZGNTWo3sm38wt7MftOlddwqzFD0Y7tx15p+A2DzFfNv9ZcFjnbreoFnd/++NnqSn9V9S3rlGtvKrLxrZz74YZmRt+/Z7UhMic33jgdLSUt12220aNmyYxo0bp1WrVkmSFi1aFPH8srIyNTY2djzq6urcLgk4JXoQfqL/4Dd6EH6i/+AVz/8YaL9+/TRs2DDV1ET+KXAwGFQwGOl+BtA96EH4if6D3+hB+In+g1c8/zs5oVBIH3zwgXJycrxeCgAAAADcH3Jmz56tqqoq1dbW6u2339btt9+upqYmTZ482e2lAAAAAMDg+tvV9u7dq7vuukuff/65zjnnHI0cOVKbNm1SQUGB20sBAAAAgMH1IWfZsmVuXxIA0IucO2ejkc2dU2RmET52kNhJDV2TsN7soZ+cf7EPlSAuJQSMKDlg7r63sLLYyAaJnfa6k+e/kwMAAAAA3YkhBwAAAIBVGHIAAAAAWIUhBwAAAIBVGHIAAAAAWIUhBwAAAIBVXN9CGgAAALBR26d7jaw019xinO2i/cedHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWYcgBAAAAYBWGHAAAAABWSfK7gBM5jiNJalOr5PhcDOJGm1ol/a0/vEQP4kT0H/xGD8JP9B/8FGv/xd2Q09zcLEnaoFd9rgTxqLm5WWlpaZ6vIdGDMNF/8Bs9CD/Rf/BTtP0XcLpjLI9Ce3u79u3bp9TUVAUCgS5fr6mpSfn5+aqrq1P//v1dqJA1/Li+4zhqbm5Wbm6uEhK8fZdlT+tBG3oj3teg//xdw4bX0NU1emoPxvvntTet0Rv7T4rvzytrdE6s/Rd3d3ISEhKUl5fn+nX79+/v2X801uie63v906PjemoP2tAb8bwG/ef/Gja8hq6s0ZN7MJ4/r71tjd7Yf1L8fl5Zo3Ni6T82HgAAAABgFYYcAAAAAFZJLC8vL/e7CK8lJiaquLhYSUnevTuPNeLj+vHKhs8ra/RcNnxebXgN3bVGvLHl82rDGr2x/yQ7Pq+sEb2423gAAAAAALqCt6sBAAAAsApDDgAAAACrMOQAAAAAsApDDgAAAACrWDvklJeXKxAIhD2ys7O7dM3169dr/Pjxys3NVSAQ0IoVK8KedxxH5eXlys3NVUpKioqLi1VdXe3a9adMmWK8ppEjR0b1GioqKnT55ZcrNTVVmZmZmjBhgnbt2hV2TigU0owZM5SRkaF+/frp5ptv1t69e127fnFxsfE6Jk2aFNXriHc9sf86s0ZXe9Dr/uvsGvRg9Og/99ag/2LD92B3rk//xYavge6t0V09aO2QI0lDhw7V/v37Ox47duzo0vVaWlo0fPhwzZs3L+Lzc+fO1eOPP6558+bpnXfeUXZ2tr7xjW+oubnZletL0vXXXx/2ml599dWoXkNVVZWmTZumTZs2qbKyUm1tbSopKVFLS0vHOQ8++KCWL1+uZcuWacOGDTp06JBuuukmHTt2zJXrS9LUqVPDXsfTTz8d1evoCXpa/3VmDalrPeh1/3V2DYkejBb9R/9Fq6d9DeR7sF16Wv91Zg2Jr4FRcSz14x//2Bk+fLhn15fkLF++vOO4vb3dyc7Odn760592ZF9++aWTlpbm/OpXv+ry9R3HcSZPnuzccsstsRcdQUNDgyPJqaqqchzHcQ4ePOgkJyc7y5Yt6zjnT3/6k5OQkOCsXr26y9d3HMe5+uqrne9973tdLz6O9fT+i7SG47jfg173X6Q1HIce7Cr6L/Y1HIf+cwPfg2O7vuPQf27ga2DsazhO9/Wg1XdyampqlJubq8LCQk2aNEkff/yxZ2vV1taqvr5eJSUlHVkwGNTVV1+tjRs3urbOunXrlJmZqSFDhmjq1KlqaGjo0vUaGxslSenp6ZKkLVu2qLW1Nex15ObmqqioKKbXceL1j3v++eeVkZGhoUOHavbs2VH9pKOnsLH/JHd70Ov+i7TGcfSge+i/zq9xHP3nLr4Hd+76x9F/7uJrYOfXOK47etDaP3k7YsQILV68WEOGDNFnn32mxx57TKNGjVJ1dbUGDBjg+nr19fWSpKysrLA8KytLn3zyiStrlJaW6o477lBBQYFqa2v1yCOP6JprrtGWLVsUDAajvp7jOJo5c6ZGjx6toqIiSX99HX369NHZZ59tvI7jr7Er15eke+65R4WFhcrOztbOnTtVVlam9957T5WVlVG/hnhlY/9J7vag1/13sjUketDtHqT/Or+GRP/11K+BfA/uOWzsP4mvgdGydsgpLS3t+P/Dhg3TlVdeqUGDBmnRokWaOXOmZ+sGAoGwY8dxjCxWEydO7Pj/RUVFuuyyy1RQUKBVq1bp1ltvjfp606dP1/bt27Vhw4bTnhvL6zjZ9adOndrx/4uKijR48GBddtll2rp1qy655JKo1ohXNvaf5G4Pet1/p1qDHvSmB+m/zq1B//XMr4F8D+45bOw/ia+B0bL67Wp/r1+/fho2bJhqamo8uf7xXTtOnHQbGhqMyd4tOTk5KigoiOk1zZgxQytXrtTatWuVl5fXkWdnZ+vo0aM6cOBA2PnRvo6TXT+SSy65RMnJyZ79t4kHNvafFHsPet1/p1ojEnqwa+i/zq8RCf3XdXwP7tz1I6H/uo6vgZ1fIxKverDXDDmhUEgffPCBcnJyPLn+8dtuf3+r7ejRo6qqqtKoUaM8WfOLL75QXV1dVK/JcRxNnz5dL730kt58800VFhaGPX/ppZcqOTk57HXs379fO3fu7NTrON31I6murlZra6tn/23igY39J0Xfg173X2fWiIQe7Br6r/NrREL/dR3fgzt3/Ujov67ja2Dn14jEsx70fGsDn8yaNctZt26d8/HHHzubNm1ybrrpJic1NdXZs2dPzNdsbm523n33Xefdd991JDmPP/648+677zqffPKJ4ziO89Of/tRJS0tzXnrpJWfHjh3OXXfd5eTk5DhNTU1dvn5zc7Mza9YsZ+PGjU5tba2zdu1a58orr3TOPffcTl/fcRznu9/9rpOWluasW7fO2b9/f8fj8OHDHed85zvfcfLy8pw1a9Y4W7duda655hpn+PDhTltbW5evv3v3bufRRx913nnnHae2ttZZtWqV85WvfMW5+OKLO3X9nqIn9t/p1nCjB73uv86sQQ/G1oP0H/0XjZ74NZDvwfTfqfA1sOd9DbR2yJk4caKTk5PjJCcnO7m5uc6tt97qVFdXd+maa9eudSQZj8mTJzuO89ctBH/84x872dnZTjAYdMaMGePs2LHDlesfPnzYKSkpcc455xwnOTnZGThwoDN58mTn008/jeo1RLq+JGfBggUd5xw5csSZPn26k56e7qSkpDg33XRTp9c53fU//fRTZ8yYMU56errTp08fZ9CgQc4DDzzgfPHFF1G9jnjXE/vvdGu40YNe919n1qAHY0P/ubMG/Rc7vgd3/fr0X+z4GujOGt3Zg4H/KggAAAAArNBrficHAAAAQO/AkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKgw5AAAAAKzCkAMAAADAKv8XNIYLnjOzaH4AAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f323cef0198>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Define the number of images for the barycenter calculation.\n",
    "n = 20\n",
    "\n",
    "#Read the images from the file.\n",
    "function read_idx(filename)\n",
    "    f = open(filename,\"r\")\n",
    "    data_layout = zeros(UInt8,4)\n",
    "    readbytes!(f,data_layout,4)\n",
    "    data_zero = reinterpret(UInt16,data_layout[1:2])\n",
    "    data_type,data_dimensions = reinterpret(UInt8,data_layout[3:4])\n",
    "    data_shape = Int32[]\n",
    "    for i = 1:data_layout[4]\n",
    "        s = zeros(UInt8,4)\n",
    "        readbytes!(f,s,4)\n",
    "        s = map(hton,reinterpret(Int32,s))\n",
    "        push!(data_shape,s[1])\n",
    "    end\n",
    "    idx_data = zeros(UInt8,cumprod(data_shape)[length(data_shape)])\n",
    "    read!(f,idx_data)\n",
    "    idx_data = reshape(idx_data,Tuple(reverse(data_shape)))\n",
    "    return(idx_data)\n",
    "end\n",
    "\n",
    "data = read_idx(\"train-images-idx3-ubyte\")\n",
    "labels = read_idx(\"train-labels-idx1-ubyte\")\n",
    "\n",
    "#Select the images.\n",
    "mask = labels .== 1\n",
    "train_ones = data[:,:,mask]\n",
    "train = train_ones[:,:,1:n]\n",
    "\n",
    "x = [i for i=1:28]\n",
    "y = reverse(x)\n",
    "f,ax = PyPlot.plt.subplots(2,5,sharey=true,sharex=true,figsize=(10,5))\n",
    "PyPlot.plt.xticks([5,10,15,20,25])\n",
    "\n",
    "for i = 1:10\n",
    "    rand_pick = rand(1:size(train_ones)[3])\n",
    "    ax[i].pcolormesh(x,y,transpose(train_ones[:,:,rand_pick]))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Barycenters using JuMP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "wasserstein_barycenter (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function single_pmf(data)\n",
    "    #Takes a list of images and extracts the probability mass function.\n",
    "    v = vec(data[:,:,1])\n",
    "    v = v./cumsum(v)[length(v)]\n",
    "    for im_k in 2:size(data)[3]\n",
    "        image = data[:,:,im_k]\n",
    "        arr = vec(image)\n",
    "        v_size = size(arr)[1]\n",
    "        v = hcat(v, arr./cumsum(arr)[length(arr)])\n",
    "    end\n",
    "    return v,size(v)[1]\n",
    "end\n",
    "\n",
    "function ms_distance(m,n)\n",
    "    #Squared Euclidean distance calculation between the pixels.\n",
    "    d = ones(m,m)\n",
    "    coor_I = []\n",
    "    for c_i in 1:n\n",
    "        append!(coor_I,ones(Int,n).*c_i)\n",
    "    end    \n",
    "    coor_J = repeat(1:n,n)\n",
    "    coor = hcat(coor_I,coor_J)\n",
    "    for i in 1:m\n",
    "        for j in 1:m\n",
    "            d[i,j] = norm(coor[i,:]-coor[j,:]).^2\n",
    "        end\n",
    "    end\n",
    "    return d\n",
    "end\n",
    "\n",
    "function wasserstein_barycenter(data)\n",
    "    M= direct_model(Mosek.Optimizer())\n",
    "\n",
    "    if length(size(data))==3\n",
    "        K = size(data)[3]\n",
    "    else\n",
    "        K = 1\n",
    "    end\n",
    "    v,N = single_pmf(data)\n",
    "    d = ms_distance(N,size(data)[2])\n",
    "    \n",
    "    #Define indices\n",
    "    M_i = 1:N\n",
    "    M_j = 1:N\n",
    "    M_k = 1:K\n",
    "\n",
    "    #Adding variables\n",
    "    M_pi = @variable(M, M_pi[i = M_i, j = M_j, k = M_k] >= 0.0)\n",
    "    M_mu = @variable(M, M_mu[i = M_i] >= 0.0)\n",
    "    \n",
    "    #Adding constraints\n",
    "    @constraint(M, c3_expr[k = M_k, j = M_j], sum(M_pi[:,j,k]) == v[j,k])\n",
    "    @constraint(M, c2_expr[k = M_k, i = M_i], sum(M_pi[i,:,k]) == M_mu[i])\n",
    "    \n",
    "    #Objective\n",
    "    W_obj = @objective(M, Min, sum(d[i,j]*M_pi[i,j,k] for i=M_i,j=M_j,k=M_k)/K)\n",
    "                    \n",
    "    return M,M_mu\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "show_barycenter (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function run_model(data)\n",
    "    @time begin\n",
    "        M,M_mu = wasserstein_barycenter(data)\n",
    "        optimize!(M)\n",
    "    end\n",
    "    println(\"Solution status = \",termination_status(M))\n",
    "    println(\"Primal objective value = \",objective_value(M))\n",
    "    mu_level = value.(M_mu)\n",
    "    return mu_level\n",
    "end\n",
    "\n",
    "function show_barycenter(bary_center)\n",
    "    bary_center = reshape(bary_center,(28,28))\n",
    "    x = [i for i=1:28]\n",
    "    y = reverse(x)\n",
    "    PyPlot.plt.pcolormesh(x,y,transpose(bary_center))\n",
    "    PyPlot.plt.title(\"Non-regularized Wasserstein Barycenter\")\n",
    "    PyPlot.plt.show()\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem\n",
      "  Name                   :                 \n",
      "  Objective sense        : min             \n",
      "  Type                   : LO (linear optimization problem)\n",
      "  Constraints            : 31360           \n",
      "  Cones                  : 0               \n",
      "  Scalar variables       : 12293904        \n",
      "  Matrix variables       : 0               \n",
      "  Integer variables      : 0               \n",
      "\n",
      "Optimizer started.\n",
      "Presolve started.\n",
      "Linear dependency checker started.\n",
      "Linear dependency checker terminated.\n",
      "Eliminator started.\n",
      "Freed constraints in eliminator : 0\n",
      "Eliminator terminated.\n",
      "Eliminator - tries                  : 1                 time                   : 0.00            \n",
      "Lin. dep.  - tries                  : 1                 time                   : 1.00            \n",
      "Lin. dep.  - number                 : 19              \n",
      "Presolve terminated. Time: 24.73   \n",
      "GP based matrix reordering started.\n",
      "GP based matrix reordering terminated.\n",
      "Problem\n",
      "  Name                   :                 \n",
      "  Objective sense        : min             \n",
      "  Type                   : LO (linear optimization problem)\n",
      "  Constraints            : 31360           \n",
      "  Cones                  : 0               \n",
      "  Scalar variables       : 12293904        \n",
      "  Matrix variables       : 0               \n",
      "  Integer variables      : 0               \n",
      "\n",
      "Optimizer  - threads                : 20              \n",
      "Optimizer  - solved problem         : the primal      \n",
      "Optimizer  - Constraints            : 17440\n",
      "Optimizer  - Cones                  : 0\n",
      "Optimizer  - Scalar variables       : 1395520           conic                  : 0               \n",
      "Optimizer  - Semi-definite variables: 0                 scalarized             : 0               \n",
      "Factor     - setup time             : 25.48             dense det. time        : 1.58            \n",
      "Factor     - ML order time          : 0.06              GP order time          : 20.69           \n",
      "Factor     - nonzeros before factor : 1.55e+06          after factor           : 1.44e+07        \n",
      "Factor     - dense dim.             : 0                 flops                  : 1.64e+10        \n",
      "Factor     - GP saved nzs           : 1.75e+06          GP saved flops         : 3.15e+09        \n",
      "ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  \n",
      "0   5.5e+03  3.2e+02  8.3e+07  0.00e+00   1.182613040e+07   0.000000000e+00   4.9e+01  53.40 \n",
      "1   6.8e-01  4.0e-02  1.0e+04  -1.00e+00  1.120231124e+07   -1.638527155e+05  6.1e-03  54.77 \n",
      "2   5.8e-02  3.4e-03  8.8e+02  2.61e+01   1.558624793e+04   -2.985004056e+03  5.2e-04  56.57 \n",
      "3   4.5e-02  2.6e-03  6.8e+02  1.08e+01   3.256531610e+03   -6.820367822e+02  4.0e-04  57.82 \n",
      "4   4.2e-02  2.5e-03  6.4e+02  5.25e+00   2.407863337e+03   -4.995626509e+02  3.8e-04  58.96 \n",
      "5   3.8e-02  2.2e-03  5.7e+02  4.40e+00   1.571023934e+03   -3.171235831e+02  3.4e-04  60.16 \n",
      "6   3.3e-02  1.9e-03  5.0e+02  3.48e+00   1.104556063e+03   -2.125212676e+02  3.0e-04  61.35 \n",
      "7   2.9e-02  1.7e-03  4.3e+02  2.91e+00   7.928892890e+02   -1.415543156e+02  2.6e-04  62.37 \n",
      "8   1.2e-02  7.1e-04  1.8e+02  2.47e+00   2.252273636e+02   -1.976506734e+01  1.1e-04  63.67 \n",
      "9   5.2e-03  3.1e-04  7.9e+01  1.44e+00   9.137869757e+01   -2.802281268e+00  4.7e-05  64.92 \n",
      "10  1.2e-03  6.9e-05  1.8e+01  1.16e+00   2.223918686e+01   2.047776185e+00   1.1e-05  67.16 \n",
      "11  7.7e-04  4.5e-05  1.2e+01  1.04e+00   1.552926524e+01   2.419699999e+00   6.9e-06  68.42 \n",
      "12  3.9e-04  2.3e-05  5.9e+00  1.02e+00   9.413697273e+00   2.768959100e+00   3.5e-06  69.85 \n",
      "13  1.3e-04  8.1e-06  2.0e+00  1.01e+00   5.274785616e+00   3.000446303e+00   1.2e-06  71.81 \n",
      "14  5.5e-05  3.3e-06  8.2e-01  1.00e+00   3.991063103e+00   3.070294704e+00   4.9e-07  73.18 \n",
      "15  2.5e-05  1.5e-06  3.8e-01  1.00e+00   3.515580593e+00   3.095346906e+00   2.2e-07  74.61 \n",
      "16  1.2e-05  7.0e-07  1.8e-01  1.00e+00   3.304330090e+00   3.106260422e+00   1.1e-07  76.16 \n",
      "17  7.1e-06  4.3e-07  1.1e-01  1.00e+00   3.229749859e+00   3.109863797e+00   6.4e-08  77.29 \n",
      "18  3.1e-06  2.4e-07  4.8e-02  1.00e+00   3.165710662e+00   3.112133239e+00   2.8e-08  78.38 \n",
      "19  1.2e-06  9.4e-08  1.9e-02  1.00e+00   3.135044854e+00   3.113925074e+00   1.1e-08  79.69 \n",
      "20  2.5e-07  1.9e-08  3.9e-03  1.00e+00   3.119045127e+00   3.114732420e+00   2.3e-09  80.96 \n",
      "21  1.2e-09  9.1e-11  1.8e-05  1.00e+00   3.114878903e+00   3.114858640e+00   1.1e-11  82.84 \n",
      "22  3.7e-12  2.5e-13  4.7e-08  1.00e+00   3.114858824e+00   3.114858771e+00   2.8e-14  83.82 \n",
      "23  1.8e-13  5.7e-14  4.2e-13  1.00e+00   3.114858772e+00   3.114858772e+00   2.9e-18  84.82 \n",
      "Basis identification started.\n",
      "Primal basis identification phase started.\n",
      "Primal basis identification phase terminated. Time: 0.13\n",
      "Dual basis identification phase started.\n",
      "Dual basis identification phase terminated. Time: 0.11\n",
      "Basis identification terminated. Time: 1.09\n",
      "Optimizer terminated. Time: 101.45  \n",
      "\n",
      "246.621270 seconds (373.68 M allocations: 35.441 GiB, 28.30% gc time)\n",
      "Solution status = OPTIMAL\n",
      "Primal objective value = 3.1148587717952356\n",
      "******\n"
     ]
    }
   ],
   "source": [
    "bary_center = run_model(train)\n",
    "println(\"******\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGxCAYAAABfrt1aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XlclXX+///nUeEACiggHBBF3HPJcUlNTcCdJC3TsbJGG3Mqtxz1UzlakplMTlNNY6XVjNZoWU1Zlma5m6blkpZa/rRwS3HLwBUB398//HHGI7igl74FH/fb7dzqXNd1XteL61zC87yv6zqXyxhjBAAAYEEp2w0AAIDrF0EEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BpASaOnWqXC6XAgICtH379gLzExMTVb9+fQudXdtcLpdSU1Mdrblt2za5XC5NnTrV0boXkr8PbNu27ZzLPPfcc3K5XPr66699pp86dUphYWFyuVzavHmzz7yTJ08qKChI3bt3vxJtXzOWLVum1NRUZWVlXXKN0aNHq0yZMo71lJubK5fL5fMoW7asbrjhBj399NM6duyYY+sqLiZOnKi33nrLdhu4TASREiw7O1ujR4+23cZ1LTo6WitWrFCXLl1st1JAUlKSJGnRokU+09evX69Dhw6pbNmyBeZ9/fXXOn78uPe1JdWyZcv01FNPXVYQeeihh7R8+XIHuzqtV69eWrFihVasWKGPP/5Y3bt3V2pqqu6//37H13WtI4iUDASREqxz5856++23tX79etutnFdOTo5yc3Ntt+GovLw8ZWdny+12q0WLFqpYsaLtlgpo1KiRypcvr8WLF/tMX7x4sWJiYtS1a9cCQSR/2eIYRPLfk6slNjZWzZs3d7yux+NRixYt1KJFC7Vv317PPPOMevXqpQ8++EA5OTmOrON6HF3JZ4zRiRMnbLdxXSGIlGCPPvqowsPD9dhjj11w2RMnTmjkyJGKj4+Xv7+/KlWqpIEDB+q3337zWa5q1apKSUnR3Llz1bhxYwUGBqpOnTr697//fVE9LV68WC6XS//5z380fPhwVapUSW63W1u3bpUkZWRk6MEHH1RsbKz8/f0VHx+vp556qkBQ2bVrl3r06KHg4GCVL19evXv31qpVqwocBklMTFRiYmKBPvr27auqVauet9f9+/drwIABqlu3rsqVK6fIyEi1bdtWX375pc9y+YdfJkyYoHHjxik+Pl5ut1uLFi0q9NDM2cPrZz7OPJSyevVqde3aVWFhYQoICFCjRo303nvvFehz5cqVatWqlQICAhQTE6ORI0de1B+kUqVKqU2bNlq+fLnP9l28eLESExOVkJBQaEipWLGi6tWr55325JNPqlmzZgoLC1NISIiaNGmiqVOn6uz7ac6fP18JCQkKCwtTYGCg4uLi1KNHD59f+i+//LIaNGigcuXKKTg4WHXq1NETTzzhU2f37t3q37+/KlWqJH9/f1WvXl3jxo1TXl6ed5mtW7fK5XLp73//u8aOHauqVavK7Xbryy+/VF5ensaOHavatWsrMDBQ5cuX14033qiJEydKOn1IZeTIkZKkypUre9+bZcuWeeu/8847atGihYKCghQcHKzOnTsXCPyFHZqJjY3V7bffrjlz5qhRo0YKDAzUDTfcoDfffPOC79f5hIaGqlSpUipV6n+/0j///HN17dpVsbGxCgwMVM2aNTVgwAAdPHiwQJ8ul0vr169X9+7dVb58edWuXVtTpkyRy+XSqlWrCqzvySeflL+/v/bu3eudNmfOHCUlJSk0NFRBQUGqW7euJkyY4PO6b775RikpKapQoYICAgLUuHFjffDBBz7LvPHGG3K5XFq6dKkefPBBhYeHKzw8XD169FBGRobPtty8ebMWLFjgfY9q1KjhnZ+Zmanhw4d7f6fFxsZq2LBhPiEr/3DX0KFD9corr6hOnTry9/fXtGnTivgO4HI4dwAT15zg4GCNHj1ajzzyiBYuXKi2bdsWupwxRrfffrsWLFigkSNH6pZbbtF3332nMWPGeIeA3W63d/n169dr+PDhevzxxxUVFaU33nhD/fr1U40aNdSmTZuL6m3kyJG6+eabNWnSJJUqVUqRkZHKyMhQs2bNVKpUKT355JOqXr26VqxYoXHjxmnbtm2aMmWKJOno0aNKSkrSr7/+qmeffVY1atTQ3Llz1atXr8vfaGf49ddfJUljxoyRx+PRkSNHNHPmTCUmJmrBggUFAs5LL72kWrVq6bnnnlNISIhq1qxZaN0VK1b4PD9+/Ljuu+8+5eXlKSwsTNLpwyWdO3dW8+bNNWnSJIWGhmrGjBnq1auXjh07pr59+0qSNm3apHbt2qlq1aqaOnWqgoKC9Morr+jtt9++qJ8xKSlJs2bN0qpVq3TzzTfr1KlTWrp0qZ599lm1adNG+/bt06ZNm1S3bl2dPHlSK1asUEpKilwul7fG9u3b9fDDD6ty5coyxmjlypV6+OGHtXv3bv3lL3+RJP30009KSUlRUlKSpk6dqpCQEP3yyy+aO3eucnJyFBAQoGnTpmnQoEF65JFH1KVLF7lcLm3dutXnPJXdu3erWbNm8vf3V2pqqqpVq6bly5fr6aef1vbt2/X666/7/HwvvPCC6tSpo+eff17BwcGqVauW0tLS9PTTT+uJJ55Q69atdfLkSf344486dOiQpNOHVA4dOqRXXnlFs2bN8o5m5YevsWPHKjU1VQ888ICeeOIJZWdna8KECWrdurVWr16t2rVrn3ebr127Vo8++qgef/xxRUZGavLkyerbt69q1qypli1bXvA9M8Z4g+ORI0e0aNEiTZs2Tffcc49Kly7tXW7r1q1q1aqV+vfvr9DQUKWnp+vvf/+72rRpo/Xr1xcISd26dVPv3r01YMAAHTt2TJ06ddJjjz2ml19+2SdI5+Tk6PXXX1ePHj0UFRUlSZo8ebIeeughtW3bVpMmTVJkZKQ2b96sH374wfu6+fPnq0uXLmrZsqVee+01BQcH65133lGPHj30n//8R/fee69PP3/84x9122236Z133tH27dv16KOP6g9/+IO++OILSdInn3yiO+64Q5GRkXrppZckSQEBAd7tcssttygjI0OjRo1S/fr19f3332vMmDHasGGDPv/8c599+L///a8iIyOVmpqqqKgo78+Fq8SgxJkyZYqRZFatWmWys7NNtWrVTNOmTc2pU6eMMcYkJCSYevXqeZefO3eukWQmTJjgU+fdd981ksxrr73mnRYXF2cCAgLM9u3bvdOOHz9uwsLCzIMPPnjB3hYtWmQkmTZt2hSY9+CDD5py5cr51DbGmOeee85IMhs3bjTGGPPyyy8bSeazzz4r8HpJZsqUKd5pCQkJJiEhocC6+vTpY+Li4nymSTJjxow5Z++5ubkmJyfHtGvXztxxxx3e6enp6UaSqV69ujl58qTPa/LnndnT2TW7detmypUrZ9asWeOdXqdOHdOoUSOTk5Pjs3xKSoqJjo42eXl5xhhjevXqZQIDA01GRoZPzTp16hhJJj09/Zw/jzHGrFu3zkgy48ePN8YYs2bNGiPJ/Pjjj8YYY6KioszEiRONMcYsWbLESDKvvPLKOevl5eWZnJwc8+STT5rIyEjv9BkzZhhJZsOGDed87UMPPWQiIiLO22+/fv1MSEiI2blzp8/0v/71r8blcpnNmzcbY4zZsmWLkWRq1apVYBt27tzZNG3a9LzrSUtLM5IKrCc9Pd2ULl3a/PnPf/aZnpWVZSIjI80999zjnTZq1ChTunRpn+UqVapkgoKCzK5du7zTjh07ZkJDQ83AgQPP21NOTo6RVOgjJSXFHD169JyvPXXqlMnJyTE//fSTkWRmz57t06ckM3bs2AKvGzVqlAkICDD79+/3Tps+fbqRZJYvX26MMSYzM9MEBwebxMRE7++YwtSoUcPcdNNNJjc312d6586dTaVKlbyvff31140kM2TIEJ/lxo8fbySZffv2eafVrl3btGvXrsC6nn76aVO6dGmzdu1an+n5++EXX3xhjPnfNq1QoYL57bffztk7riwOzZRw/v7+GjdunFavXl3osL4kLVy4UJK8n7Lz9ezZU2XLltWCBQt8pv/ud79TlSpVvM8DAgJUq1Ytnyt0cnNzfR7mrGH6O++8s0Afn376qZKSkhQTE+Pz2uTkZEnSkiVLvP/NHw4/0913332+TXFJJk2apMaNGysgIEBlypSRn5+fFixY4PNJL1/Xrl3l5+dXpPqDBg3S7Nmz9f7776tx48aSTn+S/fHHH9W7d29Jvtvy1ltv1Z49e7yjBIsWLVK7du18PsGVLl36okeHbrzxRoWHh3sPwSxevFgej8f7qb5Nmzbe80TOdX7I/Pnz1a5dO4WGhqp06dLy8/PT2LFjtW/fPu9hgEaNGsnPz08PPPCA3nrrLaWnpxfopVmzZjpw4IB69+6tWbNmFTiEIJ3eR9q1ayePx1NgHzHGePeRfN26dSvwyb9Zs2Zas2aNBg0apC+++KJIJ6TOnTtXeXl5+sMf/uCz/sDAQN1yyy0FDmUVpnHjxqpUqZL3ef5hk8KucCvM3XffrVWrVmnVqlVaunSp/vGPf2jFihVKTk72OSS3d+9e/elPf1JsbKx3361evbokFbr/FvZvcsCAAcrLy9O//vUv77SJEyeqUaNG3tGbZcuW6fDhwxowYIDPKMOZfvzxR23dulW9e/f2juicuU//8ssv3sOz+bp27erz/MYbb5Qk7dix44Lb6NNPP1XDhg3VoEEDn3Xl/844+31q3769QkNDL1gXVwZB5Dpw1113qXHjxho1alSh5w4cPHhQZcqUKXBCpcvlksfjKfAHITw8vEANt9ut48ePe5/7+fn5PM4+Bh4dHV2gxt69e/XJJ58UeG3+kPiBAwe8/RY2dOr0cOrzzz+vhx9+WM2bN9cHH3yglStXatWqVercubPPz3q+n+l8xo0bp0mTJmny5Mk+oSr/uPuIESMKbIsBAwZI8t0WHo+nQO3CphXG5XIpISFBy5cvV05OjhYtWqSEhATv/ISEBC1ZskTGGC1atEgej0d16tTxzl+xYoU6d+6s0qVL64033tBXX32lVatW6fHHH5ck73aqVauW5s2bp/DwcD388MOqVq2aatSo4T0vQzodhN944w39/PPP6t69uyIjI9WiRQufILxv3z7NnDmzwHZp2LChz3bJV9h7Mnr0aE2YMEHLli1T586dFR4ervbt22vt2rUX3F75701+sDrz8cEHHxRYf2Eu5t/P+URGRqpp06Zq2rSpbrnlFg0ZMkQvvviili5d6r2CJC8vT+3bt9esWbP0+OOPa8GCBfrmm2+857lc7P4bExOjHj166NVXX9WpU6f07bffasWKFRo0aJB3mf3790s6fc7GueRvt6FDhxbYbkOGDJFU8L07ezvlHx6+mO20d+9erV27tsC6ypcvX+i6ivpvF87iHJHrgMvl0rPPPqsOHTrotddeKzA/PDxcubm52r9/v08YMcYoIyNDN910U5HXefYJbvHx8QV6OltERIRuvPFGPfPMM4XWjImJ8fb7zTffFJh/5ols+QICApSZmVlg+sX8wZg2bZoSExP16quv+kw/fPhwocuf69NgYaZOnaonnnhCqamp+uMf/+gzLyIiQtLp82jO9X0d+SMW4eHhhf7chU07l6SkJH344Yf6+uuv9eWXXyotLc07LyEhQQcOHNCaNWu0cuVK3XHHHT6vfeedd+R2u/Xpp5/K39/fO/2///1vgfUkJCQoISFBubm5Wr16tf7xj39o8ODB8ng86tGjh1wul/r166d+/frpyJEjWrJkicaMGaOUlBRt2bJFsbGxCg8PV7NmzfTUU08V+rOcOdIgFf6e+Pn5acSIERoxYoQOHTqk+fPna+TIkerUqZN27tzpPc+gMPnvzUcffVRgXeda39WQP1qQf8Ls+vXrtWHDBk2bNs07siadHpk4l3P1PmTIEL3zzjv69NNP9dFHHyksLMxn9DH/d8auXbvOWTt/uz3xxBMFRjrynRlwL1dERITKly9f4JyhfIV96II9BJHrRPv27dWhQweNHTtWlStX9pnXrl07TZgwQdOmTdOf//xn7/QPPvhAR48eVbt27Yq8vqZNmxb5NSkpKZozZ46qV6+uChUqnHO5hIQEvffee/rss8+8h20kacaMGQWWrVq1qt5//33vpbTS6VGEr776SiEhIeftx+Vy+ZykK0nfffedVqxYUWAbFsXcuXPVv39//fGPf9SYMWMKzK9du7Zq1qyp9evXa/z48eetlX+y6d69e70jQnl5eXr33Xcvup/8Qy0vvPCCMjMzfU7CrVevnsLDw5WWlqYTJ04UOCzjcrnk5+fnc7XGsWPHznvVQZkyZdSiRQvVrFlTM2bM0Nq1a9WjRw+fZcqVK6cuXbroxIkT6tGjhzZt2qTY2FilpKRo/vz5qlmzpiND6RUqVFDPnj21Y8cOjRgxQjt27FCtWrXO+ek7f/Tnp59+Urdu3S57/U5Zt26dpNOjJdL//rCevf9Onjy5yLVbtGihZs2aKS0tTevXr9fAgQMVGBjond+6dWsFBwdr0qRJ6tmzZ6E16tatq/j4eK1bt05jx44tcg/ncq6RpJSUFD333HOqWLGi4uLiHFsfrgyCyHXk2WefVZMmTbRv3z6fyy87dOjgPUM+KytLrVq18l4106hRI913331Xpb+xY8dq3rx5atmypYYMGaLatWvrxIkT2rZtm+bMmaNJkyYpNjZWffr00QsvvKB7771X48aNU40aNfTZZ5/p888/lySfP4r33XefJk+erHvvvVf9+/fXwYMHNWHChAuGEOn0L7Onn35aY8aMUUJCgjZv3qyxY8cqPj7+kr/3JD09XT179lS1atV0//33a+XKlT7zGzVqJLfbrcmTJys5OVmdOnVS3759ValSJf3666/64YcftHbtWr3//vuSTh9mmDVrltq2basnn3xSQUFBevnll3X06NGL7qlevXqKjIzUzJkzVbFiRd1www3eeS6XS23atNHMmTMlFTw/pEuXLnrppZd077336oEHHtCBAwc0YcIEBQUF+Sz38ssv68svv1RycrKqVKmi48ePe887aN++vSTp/vvvV0hIiFq1aiWPx6M9e/Zo/PjxqlChgpo0aSLp9OGsBQsWqGXLlho8eLBq166t48ePKz09XbNnz9a///3vCx6WuvXWW/W73/1OTZo0UcWKFZWenq5//vOfqlatmqpVqyZJatCggSTpxRdf1L333is/Pz/VqVNH1atX15gxY/T4449r69at6tSpk8qXL6+MjAx98803Cg0N1ZNPPnnR2/5SZGRkePebEydO6Ntvv9W4ceMUFhamPn36SDr9nlatWlWPPvqocnNzVb58ec2aNavA+V4X65FHHlHv3r1VqlQp7+HBfCEhIfrb3/6mhx56SB06dNADDzygyMhIbdmyRRs2bNBLL70kl8ul1157TV26dFFycrL+8Ic/KCYmRocOHdKmTZu0fv36IoXnfA0aNNAHH3yg9957T/Hx8QoMDFT9+vU1bNgwzZw5U23atNHQoUPVoEED5eXlaceOHfriiy/02GOPXdKHJVwhVk+VxRVx5lUzZ7vnnnuMJJ+rZow5feXLY489ZuLi4oyfn5+Jjo42Dz/8sDl06JDPcnFxcaZLly4F6p7r6pSz5V818/777xc6f//+/WbIkCEmPj7e+Pn5mbCwMNOkSRMzatQoc+TIEe9yO3bsMN27dzflypUzwcHB5s477zRz5swxkszHH3/sU/PNN980N9xwgwkICDB169Y177777kVdNZOdnW1GjBhhKlWqZAICAkzjxo3NRx99VOC1+VfG/O1vfyvw85x91Uz+z3+ux5lXuaxfv978/ve/N5GRkcbPz894PB7Ttm1bM2nSJJ91LF++3LRo0cK43W7j8XjM//3f/5nXXnvtoq6ayff73//eSDI9evQoMO/FF180kkylSpUKfe3rr79uatWqZdxut6lWrZp59tlnzeTJk32uOlm+fLm5/fbbTZUqVYzb7Tbh4eEmKSnJ5+qNf//73yYpKclERUUZf39/ExMTY+66664CV9rs3bvXDB482FStWtW7jzRt2tSMHj3aHDt2zBjzv6tmXnjhhQL9TpgwwbRs2dJEREQYf39/U6VKFdO/f3+zY8cO7zKnTp0yjz76qImOjjalSpUyksyXX37pnf/hhx+ahIQEExISYtxut6latarp2bOnWbhwoXeZc101061btwI9tWrVqtCrP85U2FUz/v7+plq1aqZfv37m559/9ll+w4YNpn379iY4ONhUqFDB9OrVy2zbts1IMk8//bRPn5IK/Fs/0/Hjx42fn59JSUk55zKffPKJueWWW0zZsmVNUFCQqVu3rnnuued8lvn2229Njx49TMWKFb2/Z9q1a2def/117zL5V818++23Pq+dN29egffh559/Nh06dDDlypXzXrmW7/Dhw2bUqFGmdu3axt/f34SGhpobb7zRDBs2zHvlTf42feSRR875c+HKcxlz1uUMQDE1fvx4jR49Wjt27DjviXMAimbmzJnq3r27Pv/8c3Xs2NF2OyhhODSDYin/aos6deooJydHCxcu9B4iIIQAzti0aZO2bdumESNGqGnTpoQQXBEEERRLQUFBeuGFF7Rt2zZlZ2erSpUqeuyxx7jJH+CgP/3pT/rmm2/UpEmTy/4aeuBcODQDAACs4QvNAACANQQRAABgDUEEAABYc82drHrq1Cnt3r1bwcHBfO0uAADFhDFGhw8fVkxMjM8XS17INRdEdu/efVlfnw0AAOzZuXNnkb5G4ZoLIsHBwZKk1rpVZVS0W6oDAAA7cpWjZZrj/Tt+sa65IJJ/OKaM/FTGRRABAKBY+P+/DKSop1VwsioAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsKVIQSUtL00033aTg4GBFRkbq9ttv1+bNm32WSUxMlMvl8nncddddjjYNAABKhiIFkSVLlmjgwIFauXKl5s2bp9zcXHXs2FFHjx71Wa5///7as2eP9zF58mRHmwYAACVDmaIsPHfuXJ/nU6ZMUWRkpNasWaM2bdp4pwcFBcnj8TjTIQAAKLEu6xyRzMxMSVJYWJjP9OnTpysiIkL16tXTiBEjdPjw4XPWyM7OVlZWls8DAABcH4o0InImY4yGDRum1q1bq379+t7pvXv3Vnx8vDwejzZs2KCRI0dq/fr1mjdvXqF10tLS9NRTT11qGwAAoBhzGWPMpbxw4MCBmj17tpYtW6bY2NhzLrdmzRo1bdpUa9asUePGjQvMz87OVnZ2tvd5VlaWKleurER1UxmX36W0BgAArrJck6PF+liZmZkKCQm56Ndd0ojI4MGDNWvWLC1duvS8IUSSGjduLD8/P23ZsqXQIOJ2u+V2uy+lDQAAUMwVKYgYYzR48GDNnDlTixcvVnx8/AVfs3HjRuXk5Cg6OvqSmwQAACVTkYLIwIED9fbbb+vjjz9WcHCwMjIyJEmhoaEKDAzUTz/9pOnTp+vWW29VRESENm3apOHDh6tRo0Zq1arVFfkBAABA8VWkq2ZeffVVZWZmKjExUdHR0d7Hu+++K0ny9/fXggUL1KlTJ9WuXVtDhgxRx44dNX/+fJUuXfqK/AAAAKD4KvKhmfOpXLmylixZclkNAQCA6wf3mgEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1hBEAACANZd80zsAKA52jW7pSJ2c4Eu6LVcBNacccKSOJOX+uMWxWoAtjIgAAABrCCIAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABruPsugGvPzQ0dK5XQba0jdT5b38CROj+OCnakjiTVuM+xUoA1jIgAAABrCCIAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACs4aZ3ABzjKuPnSJ2tPYIcqSNJO+Y0dqRO6O8OOVLn2HF/R+oAJQUjIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBruvgvAMaVqxztTqLRxpo6kgEbO3DW3uWeHI3XmL2/oSB2gpGBEBAAAWEMQAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhTpCCSlpamm266ScHBwYqMjNTtt9+uzZs3+yyTnZ2twYMHKyIiQmXLllXXrl21a9cuR5sGAAAlQ5GCyJIlSzRw4ECtXLlS8+bNU25urjp27KijR496lxk6dKhmzpypGTNmaNmyZTpy5IhSUlKUl5fnePMAAKB4K9K9ZubOnevzfMqUKYqMjNSaNWvUpk0bZWZm6l//+pf+85//qH379pKkadOmqXLlypo/f746derkXOcAAKDYu6xzRDIzMyVJYWFhkqQ1a9YoJydHHTt29C4TExOj+vXr66uvviq0RnZ2trKysnweAADg+nDJd981xmjYsGFq3bq16tevL0nKyMiQv7+/KlSo4LNsVFSUMjIyCq2Tlpamp5566lLbAHAt2bPfkTIh1Z07j75mhQOO1NmfXc6ROu4DLkfqACXFJf9rHzRokL777ju98847F1zWGCOXq/B/fCNHjlRmZqb3sXPnzkttCQAAFDOXFEQGDx6sWbNmadGiRYqNjfVO93g8OnnypA4dOuSz/L59+xQVFVVoLbfbrZCQEJ8HAAC4PhQpiBhjNGjQIH2FbyxGAAATuklEQVT44YdauHCh4uPjfeY3adJEfn5+mjdvnnfanj17tGHDBrVs2dKZjgEAQIlRpHNEBg4cqLffflsff/yxgoODved9hIaGKjAwUKGhoerXr5+GDx+u8PBwhYWFacSIEWrQoIH3KhoAAIB8RQoir776qiQpMTHRZ/qUKVPUt29fSdILL7ygMmXK6Pe//72OHz+udu3aaerUqSpdurQjDQMAgJKjSEHEGHPBZQICAvTPf/5T//znPy+5KQAAcH3gXjMAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsOaSb3oHAGfLrRvnSJ1bojc5UkeSsk8582tu+S/xF17oIpzyc6QMUGIwIgIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALCGIAIAAKzh7rsAHLO/UZAjdYZVWO9IHUl6d39zR+qULnXKkTpRX2Y7UgcoKRgRAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1nD3XeA6V6ZaVcdqHa/oTJ3UzV2dKSTp7rjVjtRZ+nVdR+p4jhxzpI4kGccqAfYwIgIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALCGm94B1zlTNsCxWjlVTzhSp3LIb47UkaSfHboTX/g6lyN1XN9vdaSOxE3vUDIwIgIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrihxEli5dqttuu00xMTFyuVz66KOPfOb37dtXLpfL59GiRQvHGgYAACVHkYPI0aNH1bBhQ02cOPGcy3Tu3Fl79uzxPubMmXNZTQIAgJKpyDe9S05OVnJy8nmXcbvd8ng8l9wUAAC4PlyRu+8uXrxYkZGRKl++vBISEvTMM88oMjKy0GWzs7OVnZ3tfZ6VlXUlWgJwDieiyzlWyy/wpCN1fst27o7AO7PKO1Lnt9rO3H037IzfdwCuwMmqycnJmj59uhYuXKi///3vWrVqldq2besTNs6Ulpam0NBQ76Ny5cpOtwQAAK5Rjo+I9OrVy/v/9evXV9OmTRUXF6fZs2ere/fuBZYfOXKkhg0b5n2elZVFGAEA4DpxRQ7NnCk6OlpxcXHasmVLofPdbrfcbveVbgMAAFyDrvj3iBw8eFA7d+5UdHT0lV4VAAAoZoo8InLkyBFt3brV+zw9PV3r1q1TWFiYwsLClJqaqjvvvFPR0dHatm2b/vKXvygiIkJ33HGHo40DAIDir8hBZPXq1UpKSvI+zz+/o0+fPnr11Vf1/fff66233tJvv/2m6OhoJSUl6d1331VwcLBzXQMAgBKhyEEkMTFRxphzzv/8888vqyEAAHD94F4zAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALDmin/FO4Brm3+mM3fMlaSeNb91pM5XB+IdqSNJ+w868x1GFdIdKaNSDn6nUt5vvzlWC7CFEREAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANd98FrnMm7VfHar3/yS2O1Ln/jnmO1JGk99/r4Egd1ylHynDHXOAsjIgAAABrCCIAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACs4aZ3QDFVpmZ1R+qkL4pypI4kxSXtcKTO5DVtHKkjSa3u3+RInW1/q+1IHQC+GBEBAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWcPddoLg6esyRMieqnXSkjiRlngh0pE5ohaOO1JGkFStvcKRO9T3ObG8AvhgRAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYUOYgsXbpUt912m2JiYuRyufTRRx/5zDfGKDU1VTExMQoMDFRiYqI2btzoWMMAAKDkKHIQOXr0qBo2bKiJEycWOn/ChAl6/vnnNXHiRK1atUoej0cdOnTQ4cOHL7tZAABQshT5XjPJyclKTk4udJ4xRi+++KJGjRql7t27S5LefPNNRUVF6e2339aDDz54ed0CAIASxdFzRNLT05WRkaGOHTt6p7ndbiUkJOirr74q9DXZ2dnKysryeQAAgOuDo3ffzcjIkCRFRUX5TI+KitL27dsLfU1aWpqeeuopJ9sArgtHbqriSJ34ynscqSNJP//scaROtWoZjtSRpJy9YY7UKbPfmcPLuY5UAUqOK3LVjMvl8nlujCkwLd/IkSOVmZnpfezcufNKtAQAAK5Bjo6IeDynPw1lZGQoOjraO33fvn0FRknyud1uud1uJ9sAAADFhKMjIvHx8fJ4PJo3b5532smTJ7VkyRK1bNnSyVUBAIASoMgjIkeOHNHWrVu9z9PT07Vu3TqFhYWpSpUqGjp0qMaPH6+aNWuqZs2aGj9+vIKCgnTPPfc42jgAACj+ihxEVq9eraSkJO/zYcOGSZL69OmjqVOn6tFHH9Xx48c1YMAAHTp0SM2bN9cXX3yh4OBg57oGAAAlQpGDSGJioowx55zvcrmUmpqq1NTUy+kLAABcB7jXDAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACscfQr3gFcPSfKO/M5Iil8hyN1JGn3b6GO1Anxz3akjiRl7Xeo0OEjDhUCcCZGRAAAgDUEEQAAYA1BBAAAWEMQAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYQRAAAgDXcfRe4yspERTpSZ3/bk47U+WJXHUfqSNJDNyxzpM6LK9o7UkeSqm13Zjvl7jvgSB0AvhgRAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1nD3XeBqc/s7UuaGuD2O1AkonetIHUla+mtNR+oEpfs5UkeS3Ks2OlInz5xypA4AX4yIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArOGmd8BVllMlwpE6/1d5uiN1hnx7lyN1JKmUyzhSJ2S7M3UkKS8z07FaAJzHiAgAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACscTyIpKamyuVy+Tw8Ho/TqwEAACXAFbnXTL169TR//nzv89KlS1+J1QAAgGLuigSRMmXKMAoCAAAu6IoEkS1btigmJkZut1vNmzfX+PHjVa1atUKXzc7OVnZ2tvd5VlbWlWgJuGaUOXTMkTpv7WvlSJ0fWk5zpI4k1Z38sCN1Knx3yJE6kpTnWCUAV4Lj54g0b95cb731lj7//HO9/vrrysjIUMuWLXXw4MFCl09LS1NoaKj3UblyZadbAgAA1yjHg0hycrLuvPNONWjQQO3bt9fs2bMlSW+++Wahy48cOVKZmZnex86dO51uCQAAXKOuyKGZM5UtW1YNGjTQli1bCp3vdrvldruvdBsAAOAadMW/RyQ7O1s//PCDoqOjr/SqAABAMeN4EBkxYoSWLFmi9PR0ff311+rRo4eysrLUp08fp1cFAACKOccPzezatUt33323Dhw4oIoVK6pFixZauXKl4uLinF4VAAAo5hwPIjNmzHC6JAAAKKG41wwAALCGIAIAAKwhiAAAAGsIIgAAwBqCCAAAsIYgAgAArLniX/EOwFfexv/PkTr7bnakjDqpoTOFJFXWV47U4Y65wPWDEREAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANYQRAAAgDUEEQAAYA1BBAAAWEMQAQAA1hBEAACANQQRAABgDUEEAABYQxABAADWEEQAAIA1BBEAAGANQQQAAFhDEAEAANaUsd3A2YwxkqRc5UjGcjMAAOCi5CpH0v/+jl+say6IHD58WJK0THMsdwIAAIrq8OHDCg0NvejlXaao0eUKO3XqlHbv3q3g4GC5XC7b7VyzsrKyVLlyZe3cuVMhISG22ynx2N5XF9v76mJ7Xz0leVsbY3T48GHFxMSoVKmLP/PjmhsRKVWqlGJjY223UWyEhISUuJ35Wsb2vrrY3lcX2/vqKanbuigjIfk4WRUAAFhDEAEAANaUTk1NTbXdBC5N6dKllZiYqDJlrrkjbCUS2/vqYntfXWzvq4dt7euaO1kVAABcPzg0AwAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgUM6mpqXK5XD4Pj8dju60SY+nSpbrtttsUExMjl8uljz76yGe+MUapqamKiYlRYGCgEhMTtXHjRkvdFm8X2tZ9+/YtsK+3aNHCUrfFX1pamm666SYFBwcrMjJSt99+uzZv3uyzTHZ2tgYPHqyIiAiVLVtWXbt21a5duyx1XLxdzPZOTEwssI/fddddljq2hyBSDNWrV0979uzxPr7//nvbLZUYR48eVcOGDTVx4sRC50+YMEHPP/+8Jk6cqFWrVsnj8ahDhw7emzXi4l1oW0tS586dffb1OXO4GealWrJkiQYOHKiVK1dq3rx5ys3NVceOHXX06FHvMkOHDtXMmTM1Y8YMLVu2TEeOHFFKSory8vIsdl48Xcz2lqT+/fv77OOTJ0+21LFFBsXKmDFjTMOGDW23cV2QZGbOnOl9furUKePxeMxf//pX77QTJ06Y0NBQM2nSJBstlhhnb2tjjOnTp4/p1q2bpY5Kvn379hlJZsmSJcYYY3777Tfj5+dnZsyY4V3ml19+MaVKlTJz58611WaJcfb2NsaYhIQE88gjj1js6trAiEgxtGXLFsXExCg+Pl533XWXfv75Z9stXRfS09OVkZGhjh07eqe53W4lJCToq6++sthZybV48WJFRkaqVq1a6t+/v/bt22e7pRIjMzNTkhQWFiZJWrNmjXJycnz275iYGNWvX5/92wFnb+9806dPV0REhOrVq6cRI0Zcl6OrfL9sMdO8eXO99dZbqlWrlvbu3atx48apZcuW2rhxo8LDw223V6JlZGRIkqKionymR0VFafv27TZaKtGSk5PVs2dPxcXFKT09XU888YTatm2rNWvWyO12226vWDPGaNiwYWrdurXq168v6fT+7e/vrwoVKvgsGxUV5d33cWkK296S1Lt3b8XHx8vj8WjDhg0aOXKk1q9fr3nz5lns9uojiBQzycnJ3v9v0KCBbr75ZlWvXl1vvvmmhg0bZrGz64fL5fJ5bowpMA2Xr1evXt7/r1+/vpo2baq4uDjNnj1b3bt3t9hZ8Tdo0CB99913WrZs2QWXZf++fOfa3v379/f+f/369VWzZk01bdpUa9euVePGja92m9ZwaKaYK1u2rBo0aKAtW7bYbqXEy7866exPh/v27SswSgLnRUdHKy4ujn39Mg0ePFizZs3SokWLFBsb653u8Xh08uRJHTp0yGd59u/Lc67tXZjGjRvLz8/vutvHCSLFXHZ2tn744QdFR0fbbqXEyx9CPXPY9OTJk1qyZIlatmxpsbPrw8GDB7Vz50729UtkjNGgQYP04YcfauHChYqPj/eZ36RJE/n5+fns33v27NGGDRvYvy/BhbZ3YTZu3KicnJzrbh/n0EwxM2LECN12222qUqWK9u3bp3HjxikrK0t9+vSx3VqJcOTIEW3dutX7PD09XevWrVNYWJiqVKmioUOHavz48apZs6Zq1qyp8ePHKygoSPfcc4/Froun823rsLAwpaam6s4771R0dLS2bdumv/zlL4qIiNAdd9xhsevia+DAgXr77bf18ccfKzg42DuyFxoaqsDAQIWGhqpfv34aPny4wsPDFRYWphEjRqhBgwZq37695e6Lnwtt759++knTp0/XrbfeqoiICG3atEnDhw9Xo0aN1KpVK8vdX2U2L9lB0fXq1ctER0cbPz8/ExMTY7p37242btxou60SY9GiRUZSgUefPn2MMacv4R0zZozxeDzG7XabNm3amO+//95u08XU+bb1sWPHTMeOHU3FihWNn5+fqVKliunTp4/ZsWOH7baLrcK2tSQzZcoU7zLHjx83gwYNMmFhYSYwMNCkpKSwzS/Rhbb3jh07TJs2bUxYWJjx9/c31atXN0OGDDEHDx6027gFLmOMuZrBBwAAIB/niAAAAGsIIgAAwBqCCAAAsIYgAgAArCGIAAAAawgiAADAGoIIAACwhiACAACsIYgAAABrCCIAAMAagggAALDm/wFhm2Tz5D5ooAAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f323fb1ed68>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_barycenter(bary_center)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href=\"http://mosek.com\">Mosek ApS</a>. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.3",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}