wasserstein-julia/wass-bary-reg-julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )\n",
    "\n",
    "# Implementing Regularised Wasserstein Barycenter problem using JuMP in Julia.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this notebook is to implement a model to calculate the Wasserstein barycenter by solving an entropy regularised minimization problem 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-reg.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 Mosek\n",
    "using MosekTools\n",
    "\n",
    "using Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzkAAAGwCAYAAACKMqsgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3Xl8VPX1//EzWRgCDYEQskGIAQMqQVxAFikExWhUFBcUlxaqploWRUT7TdU2tv5Iv9StgkVcyqIsWhWkItSgEoqISpB9Cxg0CDGCmIQAQya5vz/6NXVyLjIzmcmd3Lyej8c8dN65cz+fJIeZnMzkjMMwDEMAAAAAwCbCrN4AAAAAAAQSTQ4AAAAAW6HJAQAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgKxFWb6Churo6OXDggERHR4vD4bB6OwgRhmFIVVWVJCcnS1hYcHtzahANUX+wGjUIK1F/sJK/9RdyTc6BAwckJSXF6m0gRJWWlkqXLl2CugY1iFOh/mA1ahBWov5gJV/rL+SanOjoaBERGSxXSoREWrwbhAq31Mgaebe+PoKJGkRD1B+sRg3CStQfrORv/YVck/PDU5MREikRDoob/8f4z3+a4qlrahAK9QerUYOwEvUHK/lZfwweAAAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArdDkAAAAALAVn5qc/Px86devn0RHR0t8fLyMHDlSdu3a5XFMZmamOBwOj8vo0aMDumkAAAAAOBWfmpzCwkIZP368rFu3TgoKCsTtdktWVpZUV1d7HJeTkyMHDx6sv8yaNSugmwYAAACAU4nw5eAVK1Z4XJ89e7bEx8dLUVGRDBkypD5v06aNJCYmBmaHAAAAAOCDRv1NTkVFhYiIxMbGeuTz58+XuLg46dWrl0yZMkWqqqoaswwAAAAAeM2nZ3J+zDAMmTx5sgwePFgyMjLq89tuu03S0tIkMTFRtm7dKrm5ubJp0yYpKCgwPY/L5RKXy1V/vbKy0t8tAX6hBmEl6g9WowZhJeoPweL3MzkTJkyQzZs3y8KFCz3ynJwcGT58uGRkZMjo0aPljTfekJUrV8qGDRtMz5Ofny8xMTH1l5SUFH+3BPiFGoSVqD9YjRqElag/BItfTc7EiRNl6dKl8uGHH0qXLl1+8tgLLrhAIiMjpbi42PTjubm5UlFRUX8pLS31Z0uA36hBWIn6g9WoQViJ+kOw+PRyNcMwZOLEibJ48WJZtWqVpKWlnfY227Ztk5qaGklKSjL9uNPpFKfT6cs2gICiBmEl6g9WowZhJeoPweJTkzN+/HhZsGCBvP322xIdHS1lZWUiIhITEyNRUVGyd+9emT9/vlx55ZUSFxcn27dvlwceeEDOP/98ufjii4PyCQAAAADAj/n0crWZM2dKRUWFZGZmSlJSUv3ltddeExGRVq1ayfvvvy+XX3659OzZU+69917JysqSlStXSnh4eFA+AQAAAAD4MZ9frvZTUlJSpLCwsFEbAgAAAIDG8HuENDyVTR7kcd19sXcjEGPanFBZQlv9vkLFh+JUVrslRmXdntMDHtzl33q1F7RMXzwxUGV/ufYVlQ2POqyyNmGtAraP8tpjKrtz6G0qc3+xL2BrAgDQLAzoo6LK7m1Oe7OhUz5W2Z8TNqmsxnB7tY3Ltt+gMtcc/Xf37ebrdZtao94MFAAAAABCDU0OAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArTBd7UfCz+6hstIReqpZdUqdyuZeNdPj+sWt9TGNcqZJNkBHH94eqbJpv7hdH7h2Y+P3hGbHrMYNhz7uid/pmpn5tZ4EeDzev3ep/vZ8/b5Z23KeU9nu3+iJLd0e3OfXmgier97orbJtg+Z5ddswL3/XVif6PtXsttO/76qypz+8QmU9/16tF9nxhV73mJ76h+AontlfZeFV+nvcbbH+nnw9rK3Kjp/leZ/10bDp6piOYa292ts9+4eqbNUufX9qJu4DfT8Z+7qebkWt4QfFc/qqbOwFa1X2YEf/fparMfRjsFtqvbrt8nNeV1mvrN+orN183/cVaDyTAwAAbKVhgwOg5aHJAQAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWaHAAAAAC20mJHSO9+qZ/KVmY9o7L2YXq+7glDjzK9dtOdnsHSjuqY6K9qfNihp9Lb3Sp7bfALKhsWpW/70l9KVVZx+c9UVnv0qH+bQ7NRu2O3yro/4P/5TMrNK9GxA706ru1+k/nWCDmd5rdR2WXtr1XZ62ctUlk7h39jyEVEak3GSt8To8dA3zPyb/rGI3X08/+ZoLKYVz72a2/w3SuXP6+yvk6Tsba36ShM9H1FnRgNEl1r+hhzz3cp1Lftssqr28qlOup99kSVpf2WWrO7iIR4lR18KVZluy6YpTKzEc9VdfrnyreO9vS4flnbXeqYruF6hHRjXNdbj7LeGtAV/MMzOQAAwFa8bV4A2BdNDgAAAABbockBAAAAYCs0OQAAAABshSYHAAAAgK20iOlqxuDzVLbm8qdVNqHkBpWdzIlWmXvXHpV1FD29IpC6r9DZIxf+SmXL/vmKyuanrVTZiHZX6hMyXQ1NZPi9a1VWXntMZUlrKlWm52nBalGLP9HhYh398qwxKjMiAvu7tiPn6UlFF00qUtkTSXqS1fz/94TKJrw2TGV1J0/6uTv84PuxesLiBc5PTY70rz7CxBHSwwfcidRQS1Ry95kqK7pAT/YV8W762YBl96ss/R7Pf0cz/ucadcznE6d7dX5vffwXPbE4WtYFdA1/8EwOAACwlVBucAA0DZocAAAAALbiU5OTn58v/fr1k+joaImPj5eRI0fKrl2eL9NyuVwyceJEiYuLk7Zt28o111wj+/fvD+imAQAAAOBUfGpyCgsLZfz48bJu3TopKCgQt9stWVlZUl1dXX/MpEmTZPHixbJo0SJZs2aNHD16VK6++mqprTV512IAAAAACDCfBg+sWOH51++zZ8+W+Ph4KSoqkiFDhkhFRYW8/PLL8sorr8jw4cNFROTVV1+VlJQUWblypVx++eWB2zkAAAAAmGjUdLWKigoREYmN/c80m6KiIqmpqZGsrKz6Y5KTkyUjI0PWrl1r2uS4XC5xuVz11ysr9TQlIJioQViJ+oPVqEFYifpDsPjd5BiGIZMnT5bBgwdLRkaGiIiUlZVJq1atpEOHDh7HJiQkSFlZmel58vPz5bHHHvN3G15xt9Gf5gnDobLKP6WoLGLX+qDsKRC+vLKd1VuwhaaowZYkIr6Tx/VDs9urY/7QaaHKBk59UGWdivSoabtpSfXn3lkc9DXabdXZe2cO0mGOHiHdNby1Ps5h//k8VtRg3MffquxLt1tlaRGt/Dp/qI+Q7px4xOothIyWdB94xpuHVHZp8b36QIf+GVUMXc9nvbtDZQ3/OOR4Ust94wW/770nTJggmzdvloUL9Q8rDRmGIQ6zb5iI5ObmSkVFRf2ltLTU3y0BfqEGYSXqD1azYw2GcoMDT3asP4QGv57JmThxoixdulRWr14tXbp0qc8TExPl5MmTcuTIEY9nc8rLy2XQIJPfpImI0+kUp9PpzzaAgKAGYSXqD1ajBmEl6g/B4tMzOYZhyIQJE+Stt96SDz74QNLS0jw+fuGFF0pkZKQUFBTUZwcPHpStW7eesskBAAAAgEDy6Zmc8ePHy4IFC+Ttt9+W6Ojo+r+ziYmJkaioKImJiZE777xTHnjgAenYsaPExsbKlClTpHfv3vXT1gAAAAAgmHxqcmbOnCkiIpmZmR757NmzZezYsSIi8vTTT0tERITcdNNNcvz4cbn00ktlzpw5Eh4eHpANAwAAAMBP8anJMUwmOzTUunVrmT59ukyfPt3vTQVa5Ht6QtqEETkqi9gUupPUvs7VL/dbd/dTJkfq17X2+OAOlaUf2haIbaGFi0hKVFmPdzynx7yW8I46ZuDU+1XW6Tn7T1JD8FXePlBlG3KeMTlS/+Kt54p7VNbDvSEQ20ID7l17VDZqxhSVrZ/0rFfn21ajJ7NdXzjO43q7z/XjY8K6apWZKf6N/nFp+sULVDY8qsqr812ZrB+DV0mUV7dF81W7bbfKohvz41h7Pb1076vne1zfOew5dUykw+THf5Mf8futv11lCdfqiW7Rsu4nNmkd+8/GBAAALUrDBgdAy0OTAwAAAMBWaHIAAAAA2ApNDgAAAABbockBAAAAYCs0OQAAAABsxacR0nZSt2m71VsQEZHwmBiV7fndOSrb+Qs9AtBsXPQLFSkq6/mH71XmPnnSuw0C/8dsXPSZ//xOZYOiPUfDjsiZoI7ptJxx0Wi8I3focdH3PLRYZREm46LnVyWr7Ow/lavMXVvr5+7gq+S/6PuFa/7SV2XfjtdvqRC36ZjH9TOlVhxrNgZsbynxF6ksa/hRkyMdXp3vxQ2DVZYuRb5uCzZVPEfXvdmI56h2J1S2ZcAsj+t6uLrI9CPpKnt54RUqO+OFYpWZnS9U8UwOAACwlUA2OACaJ5ocAAAAALZCkwMAAADAVmhyAAAAANgKTQ4AAAAAW2mx09WCLaKzntyzY2qSyq7stVVl7yT/zas19tYcU9nLT1yjsti9TLKCb46OHqCyix5cr7I9VZ1UtuueER7XW63/LHAbQ4t14Ld6otbDdy5S2Q1tv1XZnEo9dXLxjXq6Ve2+3X7uDk2p03PBfUwLO19POM19eq7K6szGXcH2wtu3V5nRrbNXty25Tk/UfWjUWyr7RfQslbklcJMeC8rPVlmXx/W/q+Y0Sc0Mz+QAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgK4yQDpLDl6aq7JWf69HQA1v7v0b3yDYqW/PHZ1U2/8EuKpuVf53K4lbsVZn7m3I/d4dAiEiIV1nVoDS/z/fdWeEq63WVHpu7JE3XUd+/TVJZ6l83q8w4etDP3QH/EZ7RU2XLxk9TWVJ4lMpeP6r/zSwe9XOV1W7b5efuYCdmtXbn6++o7NLW+i0b6hqx7r8v+avKsh55SGVmY31hrdqz9M93A2YVqex/4jY0YhX9WB1IV8RvV9ms32erLO35YpW5y/WY/lDFMzkAAAAAbMXnJmf16tUyYsQISU5OFofDIUuWLPH4+NixY8XhcHhcBgzQbywIAAAAAMHgc5NTXV0tffr0kRkzZpzymCuuuEIOHjxYf3n33XcbtUkAAAAA8JbPf5OTnZ0t2dn6dXs/5nQ6JTEx0e9NAQAAAIC/gvI3OatWrZL4+Hjp0aOH5OTkSHk5f7wOAAAAoGkEfLpadna2jBo1SlJTU6WkpEQeffRRueSSS6SoqEicTqc63uVyicvlqr9eWVkZ6C1ZImbexyqbuiJLZZU/7+7V+Y4m6X7UlVmlsr9fOEdlY9vpaVdj8/Wkt6WPRKvsfzZer7K023aqrO7kSZU1F6Fcgzv+eIbKdl2tv3fhDl0ftYb/s39WnWirMoehj3O0j9Hh0aN+r9sShXL9BdqJa/urrOaewyr7ddoqlZlNUhu6+WaVtf9tpMpqt+7wcoctU0uqwYZ2/qa9yq5u853JkY6ArhsXpkervpnzhMpu+X6KypJW6X8ztVub77TA5lZ/7mh9H3N21AELduK/ce33qOzXdz+jsuvfHqtv3JKnq918881y1VVXSUZGhowYMUKWL18uu3fvlmXLlpken5+fLzExMfWXlJSUQG8J+EnUIKxE/cFq1CCsRP0hWII+QjopKUlSU1OluFjP2hYRyc3NlYqKivpLaWlpsLcEeKAGYSXqD1ajBmEl6g/BEvQ3Az18+LCUlpZKUlKS6cedTqfpy9iApkINwkrUH6xGDcJK1B+Cxecm5+jRo7Jnz39fy1dSUiIbN26U2NhYiY2Nlby8PLnhhhskKSlJ9u3bJ7/73e8kLi5OrrvuuoBuHAAAAADM+NzkrF+/XoYNG1Z/ffLkySIiMmbMGJk5c6Zs2bJF5s2bJ99//70kJSXJsGHD5LXXXpPoaP1H7QAAAAAQaD43OZmZmWIYJmOW/s+//vWvRm0IAAAAABoj6H+Tg/9ym4zda/Omd6P42piFM3T0B7lQZcev02Nb9490q+yDYX9V2c6LX1HZul36to+OzVFZWOEGvUH45MxX9GjugRsmqmzUxAKVvfnUcJUd7mPyC4pYvcZd532ksjfu1uNNw+/W57vxmQc9ric9+4k6xqit1fuA7UVU6+/7E2ctUtl5rbwb11tTF/TZObC5MFfo1FD3CP13Ketyn1XZO/fFquylG69WWd2m7YHZGDxEFKxX2eNzRqts0D3TVLb2RFeVzf16oF7Dod8Conz2GSqruKpaZWmdPEeMP5L6jjrm4tYmdW/y48F9b72psqdv05+rrNuksxAQOv+6AQAAACAAaHIAAAAA2ApNDgAAAABbockBAAAAYCs0OQAAAABshelqLUDUYj3dKn2xPm7c2XeorPj3eq7b1iEvqeyXL/xTZQsvOkdltRUVp9omTERs3KOy73+TrrKlj1+qsg6LPtaZl+uukiiTbJDKyibrbP2U6R7Xe3Yfp47pcZ+eTsPENfuLWKm/7w9N0PXx7Z16YtC/L9L3O2vPW6iyu2ZdorLyoa1UVndSTxVEy9N9sr6fHPHWXSr7OlM/FmZd96nK0qO+UdklbXfpdU0mqUU69I9kNYbJNNNXbldZyqa1KkPT6TxVf/2vqX5IZTEl+vvZeqmuI7N7p/ZyUGdz9XEN57I9/qGevLe0x1KVuUU/Brd21KjsRLyu3dZ6GyGBZ3IAAAAA2ApNDgAAAABbockBAAAAYCs0OQAAAABshSYHAAAAgK0wXQ31anfsVlm3W/Rxw5aPVtn/9nxTZQ6nnmgE3xwa1Vtl/c7Q36fvFh1uiu0oyTM/V9nZKRM8roeJyCejnvQ8aKTIXV9cp257bIieTAR7a7VMTxbqvEwfd7vJndFZb5Sq7KWuH6hs4Fu3qqzj1XriFSAi4lizUWVd1ujjNhb1U9nu9SUqe7u7nvhXnKMfH3dd8bzK6sRQWdqcr1SmZ3bBagl/Df7Eu7AoPQlVzu7mef0BkbjnvvaI7vjqUvl71/c9sggJVxPW+jtPyNxKz/OZTYMLVTyTAyCoVIMjNDgAADSFhg2OiKgGR8R8hHTDBqe5ockBAAAAYCs0OQAAAABshSYHAAAAgK3Q5AAAAACwFZocAAAAALbCCGn8pIgz9WSNDlHHVPZ9XRuVucu/DcqeWpK0O4tV9un6dJWdKdaMkK47flxl3e//2OP6xV3uUcdsGTRbZT2fG6ey9PGfNGJ3sBP3F/tUtnNkF5Xd9+bFKnu7z8squ+5XD6qsw+yPVQacSqvln6nMdJTzt4dUFHXpoMBvCC1Oye/OV9nndzzjxS3DvTr/tMIrVZYujJAGAAAAAEvQ5AAAAACwFZ+bnNWrV8uIESMkOTlZHA6HLFmyxOPjhmFIXl6eJCcnS1RUlGRmZsq2bdsCtmEAAAAA+Ck+NznV1dXSp08fmTFjhunHp02bJk899ZTMmDFDPvvsM0lMTJTLLrtMqqqqGr1ZAAAAADgdnwcPZGdnS3Z2tunHDMOQZ555Rh5++GG5/vrrRURk7ty5kpCQIAsWLJC77767cbsFAAAAgNMI6HS1kpISKSsrk6ysrPrM6XTK0KFDZe3ataZNjsvlEpfLVX+9srIykFtCI315Y6LKNvd4TmX/e/ispthOUIRyDd4Yv15lO0t6WLAT/6WO2ixhrVp5ZP3+cbs6ruhaPRHmlrm/Vpnx6ZbAbS4EhHL9hTr3V/tVtmaRnloVP/kjlS187AmVTdiQo7K6Tdv93F3zESo1GBYVpbKyu/T0qFZX6Mmdxz/opLKkJ9cGZmNBcDyxzuothIxQqb9Q4r6sr8pa/e6gypZ1/4vJrSNPf36plQgvJqz1nHlUZc2pcgM6eKCsrExERBISEjzyhISE+o81lJ+fLzExMfWXlJSUQG4JOC1qMLgaNjjwRP3BanaswVBucODJjvUX6rxpcOwgKNPVHA6Hx3XDMFT2g9zcXKmoqKi/lJaWBmNLwClRg7AS9QerUYOwEvWHYAnoy9USE//z0qaysjJJSkqqz8vLy9WzOz9wOp3idDoDuQ3AJ9QgrET9wWrUIKxE/SFYAvpMTlpamiQmJkpBQUF9dvLkSSksLJRBg3h3XwAAAADB5/MzOUePHpU9e/bUXy8pKZGNGzdKbGysdO3aVSZNmiRTp06V9PR0SU9Pl6lTp0qbNm3k1ltvDejGAQAAAMCMz03O+vXrZdiwYfXXJ0+eLCIiY8aMkTlz5shDDz0kx48fl3HjxsmRI0ekf//+8t5770l0dHTgdg0AAAAAp+Bzk5OZmSmGYZzy4w6HQ/Ly8iQvL68x+0KQmU28Kn1AjyxcfI8eT3i0Tk/lWHn/z1UWIXr8MXzz6CL9DGjNhScs2In/6k6eVFnCtTtU1uFAG5Xtvru1ytI/Dcy+YE+dp+v7nVuvu0JlC9JWqKyuzelHryJ4vnxIj4v+/NfPenXbNWc1uK8YJ/LUFdd6RO7ivX7vrTEi0rurbNaVL1uwE/uL6NpFZcdf1j+zPNRN//t/7Pd3+L1u1NgDXq0x4WP9mL5j2Esqi3RsUlmN4TZZWT9GeuP64qtVVj77DJW13/SxX+cPFUGZrgYAAGCVhg0OgJaHJgcAAACArdDkAAAAALAVmhwAAAAAtkKTAwAAAMBWfJ6uZhcRHTuq7ItZySpLmBulMuc/m/94p10vnKuyPVnPmRypJ15duP4WlcWvZJJaMKQt/l5loxZ+oLK/TrpeZYnPrA3KngIhoktnldXJ5/rAY3oqDvBTzKb59W2/r+k3Ap+lvViispI79Peze4RTZcOiajyvF74h/T+/2SOLvaqRG/SC2eTSiJeqVZbZWk/JjHToH8n6fT5KZbFf7fZzd/ZnNklt2Tn/8Oq2H0z7q1fHRYhewy21Xt12y7BZJrc1YTLE2Ns1vHEy86DK2ovOmjueyQEAALbSsMEB0PLQ5AAAAACwFZocAAAAALZCkwMAAADAVmhyAAAAANgKTQ4AAAAAW2mxI6TrzkhS2bZB81S2sHecyp5M0VNbOn3eYETkx5v83ltEaorKjp2T6NVtv7xK961rr31KZR3DilR2tK5GZee9fZ/Ket6vR/3WebU7+Kru8+0qe2qeHhc9896ZKpto/EZlnf95QGXuL/b5tzlvDeijonGvvq6yj07o2j1r1hGVBW6IZstVN+R8lX2R41BZwrt6HG70wnVB2VOg7H1qoMqWdpiusofKBqgsfMtelVFvTcd9QI+wHbN1jMpWn7dQ39jwfBRae95CWefy/BHnzwmX6TW/Kfdxl/8VcVa6ynY82F5nZ+r7Z7PHzH8cjVFZ7J/021jg1G7p8pnKdtXoecw9I/X9XSj5uva4ylZUn6Wyfx/RNfjdb/TPt5r+2cKOeCYHAADYSsMGB0DLQ5MDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgKy33L/O26yk64/ZfrLK/dflIZVf87gmVnWgw2eXefXoC1ufFXVX2yKB3VNat1WaVDWnt3Yyf/e5qk1RPEXniu3NU9sbfLlVZ+sy1KmOSmrU6T9Xfk3ur9CS1Vx7QU/Xkfh1N2nOTyko/66KyxE91Dda01b8nOXKtZw1+OnCWOmZ7TaTKfv+LO1Xm2LpRZWi8sgF6YtPOS/QEslf6dVbZ40OuVdnZTx1SmbtY38c2RnivHio78YxLZTvOmaGy3W49OXLnr/WkorqjW/3cHYKlwx91rb7zaqzKrmnjOYlxgNMtdeI5VStv3bvqdsurzlXZvPeHqqx//50qe7TzXJWlReiJhGYO1Z1QWd7cW1XWZZ2+v8epvXF2vMpm5YxU2UsPP62yQE9c216jzzfq/XFe3bZNsX6M7Pxns1rQE0jNs5aJZ3IAAICtNGxwALQ8NDkAAAAAbCXgTU5eXp44HA6PS2Kid29kCQAAAACNFZS/yenVq5esXLmy/np4eHgwlgEAAAAAJShNTkREBM/eAAAAALBEUJqc4uJiSU5OFqfTKf3795epU6dKt27dTI91uVzicv13Qk5lZWUwtgScEjUIK1F/sBo1CCtRfwgWh2EYAR1Bsnz5cjl27Jj06NFDvvnmG3n88cdl586dsm3bNunYsaM6Pi8vTx577DGVZ8q1EuHQI/SCKSK9u8p2je+ksk9u1KN5O4TpMZeB1PvjX6isbks7laUuP6qPa6VfLuj49+eB2VgTcRs1skreloqKCmnXTn/ejRFKNRhIEakpKtvxULLKOp/5rcou6LhfZU8nf6KyWuP0A8V7vqPHW5/9+30qc39TftpzWcVu9ecweQnxnnl6lO6OzJe8Ot9etx7lfN9ePZr8yKu6Jr8d6FZZj+4HVfaPHm+oLMrk69Pr33eoLP4Nff/c9o11KgtldqvBxjhyx0CV/fHhv6tsWGuzt1TwT5jJWzF4O8FtR40ev3/LHD3Pv2te6I6Ltlv9nbjmIpWZvSWCKYfJqGmTH6Ujq/XjY+uln3q3Bjz4W38BHzyQnZ0tN9xwg/Tu3VuGDx8uy5YtExGRuXP1PHkRkdzcXKmoqKi/lJaWBnpLwE+iBmEl6g9Ws2MNBrLBQXDZsf4QGoL+ZqBt27aV3r17S3FxsenHnU6nOJ3OYG8DOCVqEFai/mA1ahBWov4QLEF/nxyXyyU7duyQpKSkYC8FAAAAAIFvcqZMmSKFhYVSUlIin3zyidx4441SWVkpY8aMCfRSAAAAAKAE/OVq+/fvl1tuuUUOHToknTp1kgEDBsi6deskNTU10EsBAAAAgBLw6WqNVVlZKTExMSE71QXWCOZkl4aoQTTUEuovrE0blVWM7KOysqwav9dYe+lfVbazJlpld6zWE9LMJKzUX5/2b21SWd2xY16dL5S1hBpslAG6Vu999XWP68Ojqvw+vbfT1c5dc6fK4l/X0/3avMl0v1NplvWHoAqZ6WoAAAD4kQt8AAAgAElEQVRWatjgAGh5aHIAAAAA2ApNDgAAAABbockBAAAAYCs0OQAAAABshSYHAAAAgK0E/H1yAADNj9mY5egFH5tk/q8xRi726rh0We/3GnV+3xLN2jrP0eHPntlTHfJsE2zjDNEjzAFYg2dyAAAAANgKTQ4AAAAAW6HJAQAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgKxFWb6AhwzBERMQtNSKGxZtByHBLjYj8tz6CiRpEQ9QfrEYNwkrUH6zkb/2FXJNTVVUlIiJr5F2Ld4JQVFVVJTExMUFfQ4QahEb9wWrUIKxE/cFKvtafw2iKttwHdXV1cuDAAYmOjhaHw9Ho81VWVkpKSoqUlpZKu3btArBD1rDi/IZhSFVVlSQnJ0tYWHBfZdncatAOtRHqa1B/1q5hh8+hsWs01xoM9a9rS1qjJdafSGh/XVnDO/7WX8g9kxMWFiZdunQJ+HnbtWsXtG8aazTN+YP926MfNNcatENthPIa1J/1a9jhc2jMGs25BkP569rS1miJ9ScSul9X1vCOP/XH4AEAAAAAtkKTAwAAAMBWwvPy8vKs3kSwhYeHS2ZmpkREBO/VeawRGucPVXb4urJG82WHr6sdPoemWiPU2OXraoc1WmL9idjj68oavgu5wQMAAAAA0Bi8XA0AAACArdDkAAAAALAVmhwAAAAAtkKTAwAAAMBWaHIAAAAA2ApNDgAAAABbockBAAAAYCs0OQAAAABshSYHAAAAgK3Q5AAAAACwFZocAAAAALZCkwMAAADAVmhyAAAAANgKTQ4AAAAAW6HJAQAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArdDkAAAAALAVmhwAAAAAthJh9QYaqqurkwMHDkh0dLQ4HA6rt4MQYRiGVFVVSXJysoSFBbc3pwbREPUHq1GDsBL1Byv5W38h1+QcOHBAUlJSrN4GQlRpaal06dIlqGtQgzgV6g9WowZhJeoPVvK1/kKuyYmOjhYRkcFypURIpMW7QahwS42skXfr6yOYqEE0RP3BatQgrET9wUr+1l/INTk/PDUZIZES4aC48X+M//ynKZ66pgahUH+wGjUIK1F/sJKf9cfgAQAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArdDkAAAAALAVmhwAAAAAtkKTAwAAAMBWaHIAAAAA2ApNDgAAAABb8anJyc/Pl379+kl0dLTEx8fLyJEjZdeuXR7HZGZmisPh8LiMHj06oJsGAAAAgFPxqckpLCyU8ePHy7p166SgoEDcbrdkZWVJdXW1x3E5OTly8ODB+susWbMCumkAAAAAOJUIXw5esWKFx/XZs2dLfHy8FBUVyZAhQ+rzNm3aSGJiYmB2CAAAAAA+aNTf5FRUVIiISGxsrEc+f/58iYuLk169esmUKVOkqqqqMcsAAAAAgNd8eibnxwzDkMmTJ8vgwYMlIyOjPr/tttskLS1NEhMTZevWrZKbmyubNm2SgoIC0/O4XC5xuVz11ysrK/3dEuAXahBWov5gNWowcMLatNFZ184q23VPR5XV/qzWqzV6/L1Gh2s3enXbUET9IVj8fiZnwoQJsnnzZlm4cKFHnpOTI8OHD5eMjAwZPXq0vPHGG7Jy5UrZsGGD6Xny8/MlJiam/pKSkuLvlgC/UIOwEvUHq1GDsBL1h2Dxq8mZOHGiLF26VD788EPp0qXLTx57wQUXSGRkpBQXF5t+PDc3VyoqKuovpaWl/mwJ8Bs1CCtRf7AaNQgrUX8IFp9ermYYhkycOFEWL14sq1atkrS0tNPeZtu2bVJTUyNJSUmmH3c6neJ0On3ZBhBQ1CCsRP3BatQgrET9IVh8anLGjx8vCxYskLfffluio6OlrKxMRERiYmIkKipK9u7dK/Pnz5crr7xS4uLiZPv27fLAAw/I+eefLxdffHFQPgEAAAAA+DGfXq42c+ZMqaiokMzMTElKSqq/vPbaayIi0qpVK3n//ffl8ssvl549e8q9994rWVlZsnLlSgkPDw/KJwAAAAAAP+bzy9V+SkpKihQWFjZqQwAAAADQGH6PkMZPi+gUp7J9d/dQ2XlX7lDZ/DNWqaxO6rxaN/39O1XW+c1IlUV/fkBl7i/5Yz8AABo68NAglbXN/FZlqe2OqGxe2kKVmQkTh8rqRP9yefnQ9iqbOeZGfcJmPFYaCIRGvRkoAAAAAIQamhwAAAAAtkKTAwAAAMBWaHIAAAAA2ApNDgAAAABbYbqaH8wmp+185EyP689f/bI65tKoAq/O790cNXPFl+p15VIdPfZtb5Wtv7a7ytz7vmzEbhAMER07quzrX/RU2YaHnvN7jTAvf/9x/8ELVbZsZT+P6+mzDqpj3F/s82tfsF7YBb1UtvSf81QW6dAPLzWG26s1fr5ptMrOjv1GZZ+8m+HV+UzpQVaSNrNYZe5yPUEL9lZ94wCVvT/xLyrrENZaZWbT0AItu833Krt/bCuV9Vgb9K20SGGt9Ne6bkW8yt7p+bZX5/P2vnJOZYrKpn50lcrO+e0XHtfdhw97tQ874pkcAAAAALZCkwMAAADAVmhyAAAAANgKTQ4AAAAAW6HJAQAAAGArNDkAAAAAbIUR0qczsI+KRs95V2W3RZ9+PPT0I91U9tLuQV5to8YdrjLn2p+pLP7qUpX96+ylKvtDpy0q65nfV2XdbmGEdKgpv06Pi/7swekqq23EJNNaLweZT0v8RGe3e2ajB2WrY6qHRarMcNd4uTtYqfQRnZmNzTUbgerteN3CPgtVFmYy87nu7g+8Op8Zs/P16T1WZSk3MkK6pTnRUf/+NybMacFOvDe4926VlVuwj5bAERWlsu7Rh1Tm7f2dt/eVv2z3lc6yZ6ps4rlDPa5/fV2yOsb99QGv9tbc8UwOAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArdDkAAAAALAVpqv9SETXLip79R/Pe3Xb86c94HE9+W/r1TFGnZ6Wkeze5uXuvPSM7lv7Lr1FZesv1NOL/tL3DZXN6nONyuo2bfdzc2hKZtNZnvu+h8o+PtJdZbsPdVLZqr4ve7XuzxytPK4v6rZcHXPOK3eqrNstG706P6yVGnvE6i0EzawLXlVZ3nBdqxEr9f07/uP4yP4qi1qipzCGso6z1qqsb8y9Kut02dcq++qbWJXFvddaZTHzPlbZkV8NVNlHj8845T5/7O9d31fZ1XKhV7eFb2orKlRWck+Gyq7o1FtlB+9wqSzp795N7nPdp+97Pzx3kcqe67za4/rAqyaqYzq+wHQ1AAAAAGh2aHIAAAAA2IpPTU5+fr7069dPoqOjJT4+XkaOHCm7du3yOMblcsnEiRMlLi5O2rZtK9dcc43s378/oJsGAAAAgFPxqckpLCyU8ePHy7p166SgoEDcbrdkZWVJdXV1/TGTJk2SxYsXy6JFi2TNmjVy9OhRufrqq6W2tjbgmwcAAACAhnwaPLBixQqP67Nnz5b4+HgpKiqSIUOGSEVFhbz88svyyiuvyPDhw0VE5NVXX5WUlBRZuXKlXH755YHbOQAAAACYaNR0tYr/mzARG/ufaSJFRUVSU1MjWVlZ9cckJydLRkaGrF271rTJcblc4nL9d9pEZWVlY7YE+IwahJWoP1iNGoSVqD8Ei99NjmEYMnnyZBk8eLBkZPxndF5ZWZm0atVKOnTo4HFsQkKClJWVmZ4nPz9fHnvsMX+3EVAnu8WrbOjTD6gsoeiEyn6WWOdx3TivpzrG+HRLI3bnnYiUZJX96ey3vbrt8KjDKps87mcqS7/b932FslCqQW90+uQ7lT1Upse2vvNhP5V1m6LHlop8q5JEk2y06PGm4dHRKjt3dZXH9T/FF6ljnun3msqeFf1vpiVobvVX+s8zdDi5ybcRFBc5a1RWepfO0lY2xW6aTiBrsG1ptcrqTI5rbpKn6bHSMk1H3eVLv9c4luhQmdlbAZj58ERbv9e1WnO7DzRTV7RVZa1MjktdYRJ6qSp1kA7P1dFXtZ4/o3bYedz/RZs5v6erTZgwQTZv3iwLF+r3W2nIMAxxOPQ/XhGR3NxcqaioqL+Ulpb6uyXAL9QgrET9wWrUIKxE/SFY/HomZ+LEibJ06VJZvXq1dOny3zfQTExMlJMnT8qRI0c8ns0pLy+XQYNMOlARcTqd4nR690ZIQDBQg7AS9QerUYOwEvWHYPHpmRzDMGTChAny1ltvyQcffCBpaWkeH7/wwgslMjJSCgoK6rODBw/K1q1bT9nkAAAAAEAg+fRMzvjx42XBggXy9ttvS3R0dP3f2cTExEhUVJTExMTInXfeKQ888IB07NhRYmNjZcqUKdK7d+/6aWsAAAAAEEw+NTkzZ84UEZHMzEyPfPbs2TJ27FgREXn66aclIiJCbrrpJjl+/LhceumlMmfOHAkPDw/IhgEAAADgp/jU5BjG6ad8tG7dWqZPny7Tp0/3e1NWCVu1QWVJq7y7bcMZZN7NQwkCk9e1xodXmRyom85zV45TWfrdnwZiVwig2i07Vbb9Qn1cNzGbpBbgvVTp2nqzwHMK259u09PV0HwlPamnTF3zZF+VHc7RL1Hu+KLJhKoA+zpXr3vD6NUqeyRuk1fnm9znfZUtljjfN9ZCmE2ZghZ2YYbKkob7/wf3v3l/jMp6CI/fzZUx+DyV/e9DL6gsTPRQr8v/PcHjevfVnwduY82M39PVAAAAACAU0eQAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFvxaYQ0QktEUqLKHn3vdZVd6NTjon/1ZabKzn5gn8rcfu0MLYZD/56kU+9vT3uzY3V61DnspSnGRZvp+s/DKrvx7vUmR0YGfzPAKRx4tE5l689aojJ9FOwmrE0blbX7836VDW59QmXvn2irsu63tdyR0Q3xTA4AAAAAW6HJAQAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWmq4WgsFatdJbSWWWJC/QUq4ucum/9yl2tsgO/7abXOMxEjpbIGHyeziK8+/3H1xNqVLbp3Dmnvd3Db96qsjSxZhoXmi/3ZX1VlvjHYpX1iPT/oe6vr1+jsq7UKsT8sdqVea7KHp45R2VDWhepLNxkWqUYer7a+K+HqOyc3L0qYzpq89B3bZXKHum0WmW3l2SrrPruDiZn3B2IbdkCz+QAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgK4yQbkJhbdqobN9v9fjeAZdvVdnLXd/ye93PTnRRWellUSo741hvlRnrt/i9LpqO2SjTL3+rx+sOGbFRZdOSn1dZlCMyMBs7BUc3PdYc+CnHR/ZXWdj4b1T2Qtf3/V6j1z8mqKzH64dUVuv3CmiuIs7Ub7uwa0K8yraNmuHV+fRgaDEdF73PfUJlX953pr7t4U1erYumE3ZhhsoOPKq/x7/v9IrKzOrjnOiDKnv9pnSVdV4d7XE9/H09rryl4JkcAAAAALbic5OzevVqGTFihCQnJ4vD4ZAlS5Z4fHzs2LHicDg8LgMGDAjYhgEAAADgp/jc5FRXV0ufPn1kxoxTPyV7xRVXyMGDB+sv7777bqM2CQAAAADe8vlvcrKzsyU7O/snj3E6nZKYmOj3pgAAAADAX0H5m5xVq1ZJfHy89OjRQ3JycqS8vDwYywAAAACAEvDpatnZ2TJq1ChJTU2VkpISefTRR+WSSy6RoqIicTqd6niXyyUul6v+emVlZaC3FDIO5uhJattyngv6ujf87IjO7tTrHhp7TGWD3piisu73fxyYjYWI5laDjov0FDzX/9N73nTOdC/PGNxJamY2Dn5RZRnP3quy9Ml6KozhrgnKnqzS3Oov0I78aqDKjqY4VLb9npkqqzHcJmfUtzXzxHfnqKznE6Uqc+//2qvzNWctuQYjOierbPekVJVNvXahyka2Pawy06lpXlrn0j+SPX7bOH3gx/aapGbX+tt1r/6Zd1e/F1QW7jB5vsFk0t4jcfr7/shdOgvP8Txfz8JfqWPOvHO3yuqO6Z8Bm7uAP5Nz8803y1VXXSUZGRkyYsQIWb58uezevVuWLVtmenx+fr7ExMTUX1JSUgK9JeAnUYOwEvUHq1GDsBL1h2AJ+gjppKQkSU1NleLiYtOP5+bmSkVFRf2ltFT/Jg0IJmoQVqL+YDVqEFai/hAsQX8z0MOHD0tpaakkJSWZftzpdJq+jA1oKtQgrET9wWrUIKxE/SFYfG5yjh49Knv27Km/XlJSIhs3bpTY2FiJjY2VvLw8ueGGGyQpKUn27dsnv/vd7yQuLk6uu+66gG4cAAAAAMz43OSsX79ehg0bVn998uTJIiIyZswYmTlzpmzZskXmzZsn33//vSQlJcmwYcPktddek+jo6MDtGgAAAABOwecmJzMzUwzDOOXH//WvfzVqQwAAAADQGEH/mxz8V9LM9So7J268ytxtdROZ9rZLZY2xf1hrlb3xqydVtuOmGSrrXT1RZamPrA3MxuAhPCZGZQce1mNzPzvnrabYjnLP/qEq2zD/XI/rRwfqsZTbh7yssp036LHmVy+8Qy9qs/GpLUlE1y4qy7z3E5X9Kf5TldUYejR0nZz6F26nM3/BpSrrvJ/7MTuLOLObyv7nvcUqu8jp7Zh678aVm+m1KkdlPX//vT5wL/d3zVWbnfrvjO4752KVffriBX6vkXL7XpW93v1dj+tbh7ykjhl8i37LhtiX7Xf/F/TpagAAAADQlGhyAAAAANgKTQ4AAAAAW6HJAQAAAGArNDkAAAAAbIXpak2o7uRJlaU+as00i66FOpvw2X0qe/q56Sp7+VY9BStv5Z0qC1u1wb/Nod6uP56tsp199dc/0Dae1FOrbn1DT9U78+EilcWf9Kzp5CV6otaRj06orEOYnviH5uvIrwaqzNtJak0hdf5XKtNzC2EnB65MVNkAp9l33bupaXvdeurpi4d/rrKdo7qqrPvez1VG/dlL53z9893efH1cR/H/58BjL+hszo4Uj+u/bKfv6+T6QzrTQ0+bPZ7JAQAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWaHAAAAAC2wnQ11Gu1TE85eumxISqbEv++yiK/Paqy2sBsq0Uzfub/V3GfW08we/Kby1T2YcF5Kuv+5A6VdfvuY5XVebEP91f7VfaLlIvl76VrPLLyuuMSHxblkRX/Qk9cS9fbQAg6mqInVFk1Sc3MwHf2quzjK7upzL3/66bYDppA/LNr5cvHB3lkvWdPlE2/etav891//a9VVle01eTIEr/OD/hj6kdXeV4XkZ3ZMz2yNectkMEbb/XIvlvWQ2Kv2h3s7TUpnskB0OQaNjgiohocAAikhg2OiPjd4ADNRcMGR0RUgyMitmtwRGhyAAAAANgMTQ4AAAAAW6HJAQAAAGArNDkAAAAAbIUmBwAAAICtMEIa9b7K05NnZsVPU9mDX12nstpt9pvKEQrOzt2nsgGfTfTqtgkffaey2i07VXaGrNXHebWCdyLOSFVZpHzk1W3bbwsP4E7QlLbfoyf61Bh6rLS3Ih364arGcPt9vkfi9L+F7I7n6wMZIe2T8JgYlZWP7uXVbTvsOK6ysNWfN3pPP+j+6rc6/FXATg+EpDDx/363ueOZHAAAAAC2QpMDAAAAwFZ8bnJWr14tI0aMkOTkZHE4HLJkyRKPjxuGIXl5eZKcnCxRUVGSmZkp27ZtC9iGAQAAAOCn+NzkVFdXS58+fWTGjBmmH582bZo89dRTMmPGDPnss88kMTFRLrvsMqmqqmr0ZgEAAADgdHwePJCdnS3Z2dmmHzMMQ5555hl5+OGH5frrrxcRkblz50pCQoIsWLBA7r777sbtFgAAAABOI6DT1UpKSqSsrEyysrLqM6fTKUOHDpW1a9eaNjkul0tcLlf99crKykBuCadQ8cuBKiu88y8qiwlrrbJdi9NVliQmU2uaiVCuQXe5/rrGPe/d1zqQE9Iaw73vSzn4gOfkvg4mdWWm04bqYGwppIRy/TXGU0fSVDa+vXdTGKcfOUtl7992kcq+Ht5eZZ9NftarNUwnsz1ToaKIMV1U5v5qv1drNBeBrEFHUrzKPvq9d9+Tr9wule2uiVPZfetGqyxyb5TKuj27yzP49jupa/Cqkqt36mmh75y1+HRbFRGRA4/WqSxxpFc3xY/Y9T4wVPRYfo/szPacdrnmvAUyeOOtFu2o6QR08EBZWZmIiCQkJHjkCQkJ9R9rKD8/X2JiYuovKSkpgdwScFrUYHA1bHDgifqD1exYg6rBEVENDkKDHesvlDRscESkRTQ4IkGaruZweM7kNgxDZT/Izc2VioqK+ktpaWkwtgScEjUIK1F/sBo1CCtRfwiWgL5cLTExUUT+84xOUlJSfV5eXq6e3fmB0+kUp9MZyG0APqEGYSXqD1ajBmEl6g/BEtBnctLS0iQxMVEKCgrqs5MnT0phYaEMGsRLVgAAAAAEn8/P5Bw9elT27NlTf72kpEQ2btwosbGx0rVrV5k0aZJMnTpV0tPTJT09XaZOnSpt2rSRW29tGa//AwAAAGAtn5uc9evXy7Bhw+qvT548WURExowZI3PmzJGHHnpIjh8/LuPGjZMjR45I//795b333pPo6OjA7RoAAAAATsHnJiczM1MMwzjlxx0Oh+Tl5UleXl5j9hV0YU49wtborUcjhxV/pY87o7PK6jZtD8zGGimsTRuV7X30PJWtuV2Pi44L17c9s+AulaU/udbP3aFFGORZb9PGvawOcZsMuL50s362t92nW1V26nsfhJKVN/VV2Ys3XeHVbdNeP6Syum36PjbliB7vPPLqa1W2pMfbXq1rdtywiyeoLNpmI6QDafsUPdY7TMwHDzV0RoR+XE6PPKGyrGEv6RsPa3BdP3RJuEO/Qr/W0GOgxcv9Vn8R49VxsL+T2f28Oq7V8s+CvBORM9PMpxm3REGZrgYAAAAAVqHJAQAAAGArNDkAAAAAbIUmBwAAAICt0OQAAAAAsBWfp6vZxVdTLlDZlvHPqex/D5+lsjlLe6nsjE2B2ZeIiJhMgAlr1UplFdfrqWkDH9STO5Yl6s+rxtDvLtzrhfEqO3v6LpW5VYKWquqWASp7+E9zPa5f2vqYOmbKwcEqi87eozImqTVftdt2q6zrH3Rmelsv13CbTDkLe0DfP8s/vTwhGq3L8nCV1WX7/y+5xtCPOHX+3jOYTFLz+1wi8vQ181Q2Y7Ke0gr7Ozbhe5XVmUzpi10e/L28e9aSBolD1blzTmzwNxICeCYHAAAAsKHGNPLNHU0OAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArdDkAAAAALCVFjtC+tExi7w6rvxktMralejjqkfpUbreOJCpp170PEePRV3W8x2TW3/s1RrLjv1MZQ/PHKOyrk+uVRnjopuHMGdrlZVO1mPS2+/RI1SrOuvfdVSn6rocl/Uvlf2m/bMqixDPEbIvVpyhjtlzU7LKRPaZZEDjhZmMco106Ic/s3HF4tC3xam1W/eVyn5ZcoXK5qWtaIrtBNWrZQNN0kNNvg+EprXnLVTZRRMnqix+uv7Zy0x4TIzKahfrn1HDHZ97XH/qux7qmJgN36jMjj/v8UwOAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArdDkAAAAALCVFjtd7Ym/jFZZSu5zKnsy6TN94z+aZBb46ITuUXPm/0Zl3f/2hcqSDno3zQOhJ7x9e5U5lrRR2efp01VWJ3pqmtnkKR92o5JPXZ51+eLfRqhj4r+g/pqK+7K+Ktt/V43KYpfqGurwznaV1VZUBGZjAXB8ZH+VRU8qVZlZ3ZtNUjM7TgyTDKfk/vqAyiqzdG1dn6Ifgw8N7KSy8otNZj4Zft5nOcy+v96dq/1m/eNS8sJd/u0DtnNob0eV7cs4obK7xy1V2Wt7sr1a4+R936nsg56vqexI3UmP63e13yQj773fI2uz9xOv1mzueCYHAAAAsKGGDU5LQpMDAAAAwFYC3uTk5eWJw+HwuCQmJgZ6GQAAAAAwFZS/yenVq5esXLmy/np4uH7dPgAAAAAEQ1CanIiICJ69AQAAAGCJoDQ5xcXFkpycLE6nU/r37y9Tp06Vbt26mR7rcrnE5XLVX6+srAzGloBTogZhJeoPVqMGYSXqD8ES8Canf//+Mm/ePOnRo4d888038vjjj8ugQYNk27Zt0rGjHrGXn58vjz32WKC3cVodX9QjbKcWjFTZ3mkx+sYmEyefucBzjF9WlB4d+NR36Sp7accglbmOtFZZ12V60Z8V7VdZ6n79eZkM4MSPWFWD/vr2xnNUttZkXLSZxo2L1s56c7zOni33uB6/h3HRPyXQ9ddwZPQfZ72ojunrrNU3HKyjuyYMV9nuWQNVFrd8r97HN+UqMxt/bnTrrBc2EfGUHp86q9uTKksK1/ef3tpf61JZZHWd3+drLoJ9H1h37JjOdu1RWXuzbE4wdhQYPLYGRnN7DDZz5r3rVHbrpgdVduW9q1X23kvPq8x0nL2XMv86xeN68lst9zE44IMHsrOz5YYbbpDevXvL8OHDZdmyZSIiMnfuXNPjc3NzpaKiov5SWqrf4wAIJmoQVqL+YDVqEFai/hAsQX8z0LZt20rv3r2luLjY9ONOp1OcTmewtwGcEjUIK1F/sBo1CCtRfwiWoL9Pjsvlkh07dkhSUlKwlwIAAACAwDc5U6ZMkcLCQikpKZFPPvlEbrzxRqmsrJQxY8YEeikAAAAAUAL+crX9+/fLLbfcIocOHZJOnTrJgAEDZN26deSv2IsAAApnSURBVJKamhropQAAAABACXiTs2jRokCfssm4932pstSbvLvt03J2g+veSZXNXh6pMdmlZWq/+7jKehbcrbL7+r2vsnHt9fSiXv++Q2VGaRuVnblQj/VM31ykMre7RmVoOvvv8vz6m05S89JLXVfq8P/p7K679RS2j3b1VVnXzodU9q9zzIfSNGQ2GbBO/J+kZub6Z/Q0pMSlLXcyEQD/xb6s7zvWL0lQWY8/36Oy3w5errLnd/9cZY6CDipLnsF91g+C/jc5AAAAANCUaHIAAAAA2ApNDgAAAABbockBAAAAYCs0OQAAAABshSYHAAAAgK0EfIQ0gOAKW/25ytJX6+PelRiT7EKVpckmr9at8+ooWO2Mmz2/n9eIHuUceN+rJF3We3XLptmfdxKF0asAgsd9+LDKeuTobLHEqSxBdgRlT3bGMzkAAAAAbIUmBwAAAICt0OQAAAAAsBWaHAAAAAC2QpMDAAAAwFZocgAAAADYCk0OAAAAAFuhyQEAAABgKzQ5AAAAAGyFJgcAAACArdDkAAAAALCVCKs30JBhGCIi4pYaEcPizSBkuKVGRP5bH8FEDaIh6g9WowZhJeoPVvK3/kKuyamqqhIRkTXyrsU7QSiqqqqSmJiYoK8hQg1Co/5gNWoQVqL+YCVf689hNEVb7oO6ujo5cOCAREdHi8PhaPT5KisrJSUlRUpLS6Vdu3YB2CFrWHF+wzCkqqpKkpOTJSwsuK+ybG41aIfaCPU1qD9r17DD59DYNZprDYb617UlrdES608ktL+urOEdf+sv5J7JCQsLky5dugT8vO3atQvaN401mub8wf7t0Q+aaw3aoTZCeQ3qz/o17PA5NGaN5lyDofx1bWlrtMT6EwndrytreMef+mPwAAAAAABbockBAAAAYCvheXl5eVZvItjCw8MlMzNTIiKC9+o81giN84cqO3xdWaP5ssPX1Q6fQ1OtEWrs8nW1wxotsf5E7PF1ZQ3fhdzgAQAAAABoDF6uBgAAAMBWaHIAAAAA2ApNDgAAAABbockBAAAAYCu2bXLy8vLE4XB4XBITExt1ztWrV8uIESMkOTlZHA6HLFmyxOPjhmFIXl6eJCcnS1RUlGRmZsq2bdsCdv6xY8eqz2nAgAE+fQ75+fnSr18/iY6Olvj4eBk5cqTs2rXL4xiXyyUTJ06UuLg4adu2rVxzzTWyf//+gJ0/MzNTfR6jR4/26fMIdc2x/rxZo7E1GOz683YNatB31F/g1qD+/MNjcGDOT/35h/vAwK3RVDVo2yZHRKRXr15y8ODB+suWLVsadb7q6mrp06ePzJgxw/Tj06ZNk6eeekpmzJghn332mSQmJspll10mVVVVATm/iMgVV1zh8Tm9++67Pn0OhYWFMn78eFm3bp0UFBSI2+2WrKwsqa6urj9m0qRJsnjxYlm0aJGsWbNGjh49KldffbXU1tYG5PwiIjk5OR6fx6xZs3z6PJqD5lZ/3qwh0rgaDHb9ebuGCDXoK+qP+vNVc7sP5DHYXppb/Xmzhgj3gT4xbOoPf/iD0adPn6CdX0SMxYsX11+vq6szEhMTjT//+c/12YkTJ4yYmBjj+eefb/T5DcMwxowZY1x77bX+b9pEeXm5ISJGYWGhYRiG8f333xuRkZHGokWL6o/5+uuvjbCwMGPFihWNPr9hGMbQoUON++67r/GbD2HNvf7M1jCMwNdgsOvPbA3DoAYbi/rzfw3DoP4Cgcdg/85vGNRfIHAf6P8ahtF0NWjrZ3KKi4slOTlZ0tLSZPTo0fLFF18Eba2SkhIpKyuTrKys+szpdMrQoUNl7dq1AVtn1apVEh8fLz169JCcnBwpLy9v1PkqKipERCQ2NlZERIqKiqSmpsbj80hOTpaMjAy/Po+G5//B/PnzJS4uTnr16iVTpkzx6TcdzYUd608ksDUY7PozW+MH1GDgUH/er/ED6i+weAz27vw/oP4Ci/tA79f4QVPUoG3f8rZ///4yb9486dGjh3zzzTfy+OOPy6BBg2Tbtm3SsWPHgK9XVlYmIiIJCQkeeUJCgnz55ZcBWSM7O1tGjRolqampUlJSIo8++qhccsklUlRUJE6n0+fzGYYhkydPlsGDB0tGRoaI/OfzaNWqlXTo0EF9Hj98jo05v4jIbbfdJmlpaZKYmChbt26V3Nxc2bRpkxQUFPj8OYQqO9afSGBrMNj1d6o1RKjBQNcg9ef9GiLUX3O9D+QxuPmwY/2JcB/oK9s2OdnZ2fX/37t3bxk4cKB0795d5s6dK5MnTw7aug6Hw+O6YRgq89fNN99c//8ZGRnSt29fSU1NlWXLlsn111/v8/kmTJggmzdvljVr1pz2WH8+j1OdPycnp/7/MzIyJD09Xfr27SsbNmyQCy64wKc1QpUd608ksDUY7Pr7qTWoweDUIPXn3RrUX/O8D+QxuPmwY/2J/P927pilkTUK4/h7i1FEJBAQMhoUsbERUSstUlhYpbIRq7QK0SYfQOys/Ap+gxSCnTixEUFwihgsFKOxCAhBMBAxCs8We1fWu3PXSSaza8b/D6ZJhnPyMg/vcEgy7IGtivTP1X7W399vJicnzeXlZSj1fzy147+T7v39/S+TfafYtm1GR0fbWtP6+rrZ29szjuOYZDL59noikTDNZtM8PDy8O7/VdfxffS8zMzPGsqzQrs1nEMX8GdN+BsPO3+96eCGDwZA//z28kL/guAf7q++F/AXHHui/h5ewMvhlhpzn52dzcXFhbNsOpf6Pr91+/qqt2Wyao6MjMz8/H0rPWq1m7u7uWlqTJJPNZk0+nzeHh4dmbGzs3fuzs7PGsqx366hWq+b8/NzXOj6q76VUKpmXl5fQrs1nEMX8GdN6BsPOn58eXshgMOTPfw8v5C847sH+6nshf8GxB/rv4SW0DIb+aIO/JJfLqVAo6Pr6WicnJ0qn0xoYGNDNzU3bNev1ulzXleu6MsZoZ2dHruvq9vZWkrS9va1YLKZ8Pq9isaiVlRXZtq3Hx8fA9ev1unK5nI6Pj1Uul+U4jubm5jQ8POy7viStra0pFoupUCioWq2+HY1G4+2c1dVVJZNJHRwc6OzsTAsLC5qamtLr62vg+ldXV9ra2tLp6anK5bL29/c1MTGh6elpX/W7RTfm76Menchg2Pnz04MMtpdB8kf+WtGNeyD3YPL3O+yB3bcHRnbIWV5elm3bsixLQ0NDWlpaUqlUClTTcRwZY345MpmMpO+PENzc3FQikVBvb69SqZSKxWJH6jcaDS0uLmpwcFCWZWlkZESZTEaVSqWlNXjVN8Zod3f37Zynpydls1nF43H19fUpnU777vNR/UqlolQqpXg8rp6eHo2Pj2tjY0O1Wq2ldXx23Zi/j3p0IoNh589PDzLYHvLXmR7kr33cg4PXJ3/tYw/sTI8/mcF//v1AAAAAABAJX+Y/OQAAAAC+BoYcAAAAAJHCkAMAAAAgUhhyAAAAAEQKQw4AAACASGHIAQAAABApDDkAAAAAIoUhBwAAAECkMOQAAAAAiBSGHAAAAACRwpADAAAAIFIYcgAAAABEyjcOE5eMr5LZVgAAAABJRU5ErkJggg==",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f7ac7859208>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Define the number of images for the barycenter calculation.\n",
    "n = 2\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 .== 3\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": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "reg_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 reg_wasserstein_barycenter(data,lambda,relgap)\n",
    "    #Calculation of wasserstein barycenter by solving an entropy regularised minimization problem.\n",
    "    #Direct mode model\n",
    "    #M = direct_model(Mosek.Optimizer(MSK_DPAR_INTPNT_CO_TOL_REL_GAP=relgap))\n",
    "    \n",
    "    #Automatic mode model\n",
    "    M = Model(with_optimizer(Mosek.Optimizer,MSK_DPAR_INTPNT_CO_TOL_REL_GAP=relgap))\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",
    "    if lambda==nothing\n",
    "        lambda = 60/median(vec(d))\n",
    "    end\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",
    "    #Auxiliary variable for the conic constraint\n",
    "    M_aux = @variable(M,M_aux[i = M_i, j = M_j, k = M_k])\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",
    "    #Adding conic constraint\n",
    "    @constraint(M,cExp_cone[i=M_i, j=M_j, k=M_k],[M_aux[i,j,k],M_pi[i,j,k],1] in MOI.ExponentialCone())\n",
    "    \n",
    "    #Non-linear objective in the case of Regularized barycenters.\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) - \n",
    "            sum(M_aux[i,j,k] for i=M_i,j=M_j,k=M_k)/lambda)/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_regularised_model(data,lambda=nothing,relgap=1e-7)\n",
    "    @time begin\n",
    "        M,M_mu = reg_wasserstein_barycenter(data,lambda,relgap)\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(\"Regularized 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                   : CONIC (conic optimization problem)\n",
      "  Constraints            : 3691072         \n",
      "  Cones                  : 1229312         \n",
      "  Scalar variables       : 6147344         \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.02            \n",
      "Lin. dep.  - number                 : 1               \n",
      "Presolve terminated. Time: 7.65    \n",
      "Problem\n",
      "  Name                   :                 \n",
      "  Objective sense        : min             \n",
      "  Type                   : CONIC (conic optimization problem)\n",
      "  Constraints            : 3691072         \n",
      "  Cones                  : 1229312         \n",
      "  Scalar variables       : 6147344         \n",
      "  Matrix variables       : 0               \n",
      "  Integer variables      : 0               \n",
      "\n",
      "Optimizer  - threads                : 20              \n",
      "Optimizer  - solved problem         : the primal      \n",
      "Optimizer  - Constraints            : 1138\n",
      "Optimizer  - Cones                  : 1229312\n",
      "Optimizer  - Scalar variables       : 3687936           conic                  : 3687936         \n",
      "Optimizer  - Semi-definite variables: 0                 scalarized             : 0               \n",
      "Factor     - setup time             : 2.62              dense det. time        : 0.00            \n",
      "Factor     - ML order time          : 0.01              GP order time          : 0.00            \n",
      "Factor     - nonzeros before factor : 2.79e+05          after factor           : 3.41e+05        \n",
      "Factor     - dense dim.             : 0                 flops                  : 1.17e+08        \n",
      "ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  \n",
      "0   6.3e+02  5.1e+02  2.4e+07  0.00e+00   2.213290906e+07   -1.586952925e+06  1.0e+00  16.82 \n",
      "1   2.2e+02  1.8e+02  1.4e+07  -9.87e-01  1.997448319e+07   -3.251325127e+06  3.5e-01  21.07 \n",
      "2   1.4e+02  1.1e+02  1.0e+07  -8.50e-01  1.706335688e+07   -4.132584961e+06  2.2e-01  26.16 \n",
      "3   8.6e+01  6.9e+01  6.1e+06  -5.35e-01  1.255536087e+07   -4.350508428e+06  1.4e-01  33.89 \n",
      "4   3.9e+01  3.1e+01  2.3e+06  -5.91e-02  6.634369363e+06   -3.290261645e+06  6.2e-02  38.48 \n",
      "5   1.6e+01  1.3e+01  6.8e+05  4.56e-01   3.047700570e+06   -1.738022157e+06  2.5e-02  42.26 \n",
      "6   9.2e+00  7.4e+00  3.2e+05  6.93e-01   1.871313131e+06   -1.107859520e+06  1.5e-02  45.87 \n",
      "7   3.8e+00  3.0e+00  9.0e+04  7.82e-01   8.078383522e+05   -5.102774510e+05  5.9e-03  50.28 \n",
      "8   1.0e+00  8.1e-01  1.4e+04  8.51e-01   2.027574293e+05   -1.783130565e+05  1.6e-03  54.31 \n",
      "9   4.2e-01  3.4e-01  4.0e+03  8.87e-01   8.245554128e+04   -8.500164850e+04  6.6e-04  58.08 \n",
      "10  1.2e-01  9.9e-02  6.9e+02  9.00e-01   2.254804899e+04   -2.962342498e+04  2.0e-04  61.72 \n",
      "11  3.6e-02  2.9e-02  1.2e+02  9.02e-01   5.918580033e+03   -1.024328330e+04  5.7e-05  65.07 \n",
      "12  1.1e-02  8.5e-03  2.0e+01  9.00e-01   1.600621502e+03   -3.421588743e+03  1.7e-05  68.39 \n",
      "13  3.3e-03  2.6e-03  3.7e+00  9.03e-01   4.338656334e+02   -1.203318166e+03  5.2e-06  71.91 \n",
      "14  9.2e-04  7.4e-04  5.9e-01  9.15e-01   9.099381179e+01   -3.977761008e+02  1.5e-06  75.28 \n",
      "15  2.5e-04  2.0e-04  8.5e-02  9.22e-01   1.264148878e+00   -1.351942993e+02  3.9e-07  78.65 \n",
      "16  7.1e-05  5.7e-05  1.4e-02  9.31e-01   -1.815177735e+01  -5.892407032e+01  1.1e-07  81.93 \n",
      "17  1.6e-05  1.3e-05  1.6e-03  9.35e-01   -2.349891579e+01  -3.326506346e+01  2.6e-08  85.27 \n",
      "18  3.2e-06  2.6e-06  1.5e-04  9.38e-01   -2.441821029e+01  -2.642392235e+01  5.1e-09  88.70 \n",
      "19  7.5e-07  6.0e-07  1.7e-05  9.44e-01   -2.452958618e+01  -2.501499156e+01  1.2e-09  92.04 \n",
      "20  1.8e-07  1.4e-07  2.0e-06  9.49e-01   -2.454255521e+01  -2.466071465e+01  2.8e-10  95.32 \n",
      "21  3.6e-08  2.9e-08  1.9e-07  9.48e-01   -2.454079676e+01  -2.456580776e+01  5.7e-11  99.05 \n",
      "22  9.3e-09  7.5e-09  2.6e-08  9.50e-01   -2.453852209e+01  -2.454516580e+01  1.5e-11  102.62\n",
      "23  2.3e-09  1.8e-09  3.3e-09  9.51e-01   -2.453774799e+01  -2.453943385e+01  3.6e-12  106.06\n",
      "24  8.2e-10  5.0e-10  4.7e-10  9.51e-01   -2.453754484e+01  -2.453801397e+01  9.8e-13  109.44\n",
      "25  8.0e-10  4.1e-10  3.6e-10  9.54e-01   -2.453752652e+01  -2.453792006e+01  8.1e-13  120.02\n",
      "26  1.3e-09  3.7e-10  3.1e-10  9.54e-01   -2.453751783e+01  -2.453787528e+01  7.4e-13  135.61\n",
      "27  1.2e-09  3.6e-10  2.9e-10  9.54e-01   -2.453751378e+01  -2.453785435e+01  7.0e-13  149.61\n",
      "28  1.4e-09  3.4e-10  2.7e-10  9.54e-01   -2.453750992e+01  -2.453783439e+01  6.7e-13  161.23\n",
      "29  1.4e-09  3.1e-10  2.3e-10  9.54e-01   -2.453750254e+01  -2.453779623e+01  6.0e-13  173.52\n",
      "30  1.4e-09  2.9e-10  2.1e-10  9.54e-01   -2.453749919e+01  -2.453777893e+01  5.7e-13  186.35\n",
      "31  1.4e-09  2.8e-10  2.1e-10  9.54e-01   -2.453749759e+01  -2.453777069e+01  5.6e-13  199.37\n",
      "32  1.4e-09  2.8e-10  2.1e-10  9.54e-01   -2.453749749e+01  -2.453777019e+01  5.6e-13  215.66\n",
      "33  1.4e-09  2.8e-10  2.1e-10  9.54e-01   -2.453749746e+01  -2.453777006e+01  5.6e-13  233.18\n",
      "34  1.4e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749669e+01  -2.453776604e+01  5.5e-13  246.21\n",
      "35  1.4e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749649e+01  -2.453776505e+01  5.5e-13  261.34\n",
      "36  1.5e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749644e+01  -2.453776480e+01  5.5e-13  279.47\n",
      "37  1.5e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749643e+01  -2.453776474e+01  5.5e-13  298.25\n",
      "38  1.5e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749641e+01  -2.453776462e+01  5.5e-13  315.76\n",
      "39  1.5e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749640e+01  -2.453776460e+01  5.5e-13  335.01\n",
      "40  1.5e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749640e+01  -2.453776457e+01  5.5e-13  354.66\n",
      "41  1.5e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749640e+01  -2.453776457e+01  5.5e-13  373.38\n",
      "42  1.5e-09  2.8e-10  2.0e-10  9.54e-01   -2.453749640e+01  -2.453776457e+01  5.5e-13  392.86\n",
      "Optimizer terminated. Time: 419.08  \n",
      "\n",
      "649.126111 seconds (508.56 M allocations: 33.802 GiB, 22.71% gc time)\n",
      "Solution status = SLOW_PROGRESS\n",
      "Primal objective value = -24.5374963989306\n",
      "******\n"
     ]
    }
   ],
   "source": [
    "bary_center = run_regularised_model(train)\n",
    "println(\"******\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGxCAYAAABfrt1aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xt01PWd//HXJJlMLoQhISSTCASkIAjRelsQkQTFQFhERQp4ORvcymK5FATWI1YFPAqrbbVnyyrYrYAWCq4XsIuLxQJRpKyAKwKiBeUqBJRLbpDrfH5/+MvUkSAEPuGTy/NxzpyT+c533vOZT76ZvOb7/c68PcYYIwAAAAciXA8AAAA0XwQRAADgDEEEAAA4QxABAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEETQrCxYskMfjCV2ioqKUlpamkSNHaufOna6HF2bGjBnyeDzW644aNUodOnSwXvdsOnTooFGjRp11ve/+fjwej+Lj49WtWzfNnDlTpaWl9T/QBub555/XggULXA8DqDdRrgcAuDB//nx17dpVZWVl+uCDD/TUU09pzZo1+uyzz5SYmOh6ePXqscce08SJE10P4wcNGzZMU6ZMkSSVlJQoPz9fTzzxhD755BO9/vrrjkd3cT3//PNKTk4+pxAHNEYEETRLPXr00LXXXitJys7OVnV1taZPn65ly5bpvvvuczy6+nHy5EnFxcWpU6dOrodyVqmpqerVq1foev/+/bV3714tWrRIZWVliomJueDHqJmP5sgYo7KyMsXGxroeCsChGUBSKJQcPnz4tNuWLl2q66+/XvHx8WrRooUGDBig//u//zttvd/97nfq0qWLfD6fLr/8ci1evPi0wyBr166Vx+PR2rVrw+67Z88eeTyes+6CX7p0qXJycpSWlqbY2Fh169ZNDz/88GmHLEaNGqUWLVpo69atysnJUUJCgm6++ebQbd8dU80hoNou330XXlFRoSeffFJdu3aVz+dTmzZtdN999+nrr78Oe+zKyko99NBDCgQCiouLU58+ffThhx/+4PM6F36/Xx6PR5GRkaFlq1at0m233aa2bdsqJiZGP/rRjzRmzBh98803YfeteY4fffSRhg0bpsTERHXq1EmvvPKKPB6P/vrXv572eE888YS8Xq8OHjwYWrZy5UrdfPPN8vv9iouLU7du3TR79uyw+23atElDhgxRUlKSYmJidNVVV+nVV18NW6fmEOGaNWv0s5/9TMnJyWrdurWGDh0a9ngdOnTQ9u3blZ+fH/qdfPd3V1RUpKlTp6pjx46Kjo7WJZdcokmTJp22PXg8Ho0fP15z585Vt27d5PP5tHDhwnOffKAesUcEkLR7925JUpcuXcKWz5o1S48++qjuu+8+Pfroo6qoqNAvf/lL3Xjjjfrwww91+eWXS5JefPFFjRkzRnfeeaeee+45FRYWaubMmSovL7c6zp07d2rQoEGaNGmS4uPj9dlnn+npp5/Whx9+qNWrV4etW1FRoSFDhmjMmDF6+OGHVVVVVWvN+++/XwMHDgxb9sYbb+iXv/ylunfvLkkKBoO67bbb9P777+uhhx5S7969tXfvXk2fPl3Z2dnatGlT6N316NGj9fLLL2vq1Km65ZZbtG3bNg0dOlTFxcXn/DyNMaHx1hyaWbhwoUaOHCmv1xta74svvtD111+v+++/X36/X3v27NGzzz6rPn36aOvWrWHrStLQoUM1cuRIPfDAAyotLVVubq4eeugh/cd//Ieuv/760HpVVVWaN2+e7rjjDqWnp0uSfv/732v06NHKysrS3LlzlZKSor/97W/atm1b6H5r1qzRwIED1bNnT82dO1d+v19LlizRiBEjdPLkydMOr9x///36x3/8Ry1evFj79+/Xv/7rv+ree+8N/S7ffPNNDRs2TH6/X88//7wkyefzSfp2j05WVpYOHDigRx55RFdccYW2b9+uxx9/XFu3btW7774bdo7RsmXL9P777+vxxx9XIBBQSkrKOf8+gHplgGZk/vz5RpLZsGGDqaysNMXFxWblypUmEAiYvn37msrKytC6+/btM1FRUWbChAlhNYqLi00gEDDDhw83xhhTXV1tAoGA6dmzZ9h6e/fuNV6v12RkZISWrVmzxkgya9asCVt39+7dRpKZP39+aNn06dPND/2JBoNBU1lZafLz840ks2XLltBteXl5RpJ56aWXTrtfXl5e2Ji+7/333zcxMTHmnnvuMcFg0BhjzB//+Ecjybz++uth627cuNFIMs8//7wxxpgdO3YYSebBBx8MW2/RokVGksnLyzvj49aQVOslNzfXlJSUnPF+NfOxd+9eI8ksX748dFvNXD7++OOn3W/69OkmOjraHD58OLRs6dKlRpLJz883xnz7O2/ZsqXp06dPaE5q07VrV3PVVVeFbUfGGDN48GCTlpZmqqurjTF/3w7Hjh0btt4zzzxjJJlDhw6FlnXv3t1kZWWd9lizZ882ERERZuPGjWHLX3vtNSPJvP3226Flkozf7zfHjh0749gBVzg0g2apV69e8nq9SkhI0MCBA5WYmKjly5crKurvOwnfeecdVVVV6Z/+6Z9UVVUVusTExCgrKyt0eOXzzz9XQUGBhg8fHvYY7du31w033GB13F9++aXuvvtuBQIBRUZGyuv1KisrS5K0Y8eO09a/884761R/x44dGjJkiHr37q2XXnop9I76v//7v9WqVSvdeuutYXPx4x//WIFAIDQXa9askSTdc889YXWHDx8eNrdnM3z4cG3cuFEbN27Ue++9p3//93/Xpk2bNHDgwLC9TEeOHNEDDzygdu3aKSoqSl6vVxkZGXWaj5/97GeSvj20VmPOnDnKzMxU3759JUnr169XUVGRxo4de8ZPMu3atUufffZZ6Ll/d54GDRqkQ4cO6fPPPw+7z5AhQ8KuX3HFFZKkvXv3/vAE6dvfSY8ePfTjH/847LEGDBhQ6+G/m266qcmfiI3GiUMzaJZefvlldevWTcXFxVq6dKnmzZunu+66S//zP/8TWqfmfJHrrruu1hoREd/m+KNHj0r69gTL70tNTQ0d9rlQJSUluvHGGxUTE6Mnn3xSXbp0UVxcnPbv36+hQ4fq1KlTYevHxcWpZcuW51z/4MGDGjhwoNq2bas33nhD0dHRodsOHz6sEydOhC37rppzMmrmIhAIhN0eFRWl1q1bn/NY2rRpEzpvR5JuvPFGtWnTRnfddZcWLFigMWPGKBgMKicnRwcPHtRjjz2mzMxMxcfHKxgMqlevXqfNhySlpaWdtiw1NVUjRozQvHnz9PDDD2v79u16//33NW/evNA6NefBtG3b9oxjrtlepk6dqqlTp9a6zvfPXfn+nNQcdqlt7LU93q5du047/HSmx6rtuQMNAUEEzVK3bt1C/+j69eun6upq/ed//qdee+01DRs2TJKUnJwsSXrttddC77JrU/PPpLYTXQsKCsKu13za4/vnjnz/n0ZtVq9erYMHD2rt2rWhvSCSdOLEiVrXr8t3kBQVFWnQoEEKBoN6++235ff7w26vOZly5cqVtd4/ISFB0t/noqCgQJdcckno9qqqqlBIOV81ewu2bNkiSdq2bZu2bNmiBQsWKC8vL7Terl27zljjTHMyceJEvfLKK1q+fLlWrlypVq1ahe3VadOmjSTpwIEDZ6xds71MmzZNQ4cOrXWdyy677Iz3r6vk5GTFxsbqpZde+sHx1KiP76QBbCCIAJKeeeYZvf7663r88cc1dOhQRUREaMCAAYqKitIXX3zxg4c4LrvsMgUCAb366quaPHlyaPm+ffu0fv360MmOkkKfePjkk080YMCA0PK33nrrrGOs+UdS8665xnffuZ+PiooK3XHHHdqzZ4/WrVtX67v+wYMHa8mSJaqurlbPnj3PWCs7O1uStGjRIl1zzTWh5a+++uoZT5Y9Vx9//LEkhU6ytDkf11xzjXr37q2nn35a27Zt07/8y78oPj4+dHvv3r3l9/s1d+5cjRw5stZ/6pdddpk6d+6sLVu2aNasWXUew5n4fL5a95AMHjxYs2bNUuvWrdWxY0drjwdcbAQRQFJiYqKmTZumhx56SIsXL9a9996rDh066IknntAvfvELffnll6FzSQ4fPqwPP/xQ8fHxmjlzpiIiIjRz5kyNGTNGw4YN0z//8z/rxIkTmjlzptLS0kKHcKRvD1n0799fs2fPVmJiojIyMvSXv/xFb7zxxlnH2Lt3byUmJuqBBx7Q9OnT5fV6tWjRotAegvP14IMPavXq1Zo1a5ZKSkq0YcOG0G1t2rRRp06dNHLkSC1atEiDBg3SxIkT9Q//8A/yer06cOCA1qxZo9tuu0133HGHunXrpnvvvVe/+c1v5PV61b9/f23btk2/+tWv6nSY6PDhw6FxlJWV6eOPP9aTTz6pVq1ahb7npWvXrurUqZMefvhhGWOUlJSkP/3pT1q1atV5zcPEiRM1YsQIeTwejR07Nuy2Fi1a6Ne//rXuv/9+9e/fX6NHj1Zqaqp27dqlLVu2aM6cOZK+DUG5ubkaMGCARo0apUsuuUTHjh3Tjh079NFHH+m//uu/6jyuzMxMLVmyREuXLtWll16qmJgYZWZmatKkSXr99dfVt29fPfjgg7riiisUDAa1b98+/fnPf9aUKVN+MDQCDYbrs2WBi6nm0wrf/6SBMcacOnXKtG/f3nTu3NlUVVWFli9btsz069fPtGzZ0vh8PpORkWGGDRtm3n333bD7v/jii+ZHP/qRiY6ONl26dDEvvfSSue2228xVV10Vtt6hQ4fMsGHDTFJSkvH7/ebee+81mzZtOqdPzaxfv95cf/31Ji4uzrRp08bcf//95qOPPjrtvnl5eSY+Pr7WOfj+p2aysrLO+EmV737KpbKy0vzqV78yV155pYmJiTEtWrQwXbt2NWPGjDE7d+4MrVdeXm6mTJliUlJSTExMjOnVq5f561//ajIyMs7rUzNer9dceuml5r777jO7du0KW/fTTz81t9xyi0lISDCJiYnmJz/5idm3b5+RZKZPn37aXH799ddnfNzy8nLj8/nMwIEDz7jO22+/bbKyskx8fLyJi4szl19+uXn66afD1tmyZYsZPny4SUlJMV6v1wQCAXPTTTeZuXPnhtY503ZY26eq9uzZY3JyckxCQoKRFPa7KykpMY8++qi57LLLTHR0tPH7/SYzM9M8+OCDpqCgIGxOx40bd8bnBbjkMcaYixt9gObhxIkT6tKli26//Xa9+OKLroeDs/jTn/6kIUOGaMWKFRo0aJDr4QDNBkEEsKCgoEBPPfWU+vXrp9atW2vv3r167rnn9Nlnn2nTpk2hLwZDw/Ppp59q7969mjhxouLj4/XRRx9xYidwEXGOCGCBz+fTnj17NHbsWB07dkxxcXHq1auX5s6dSwhp4MaOHasPPvhAV199tRYuXEgIAS4y9ogAAABn+GZVAADgDEEEAAA4QxABAADONLiTVYPBoA4ePKiEhAROGgMAoJEwxqi4uFjp6elhX+R4Ng0uiBw8eFDt2rVzPQwAAHAe9u/f/4MNIr+vwQWRmuZZfTRIUaq9qyQAAGhYqlSpdXo79H/8XDW4IFJzOCZKXkV5CCIAADQK///LQOp6WgUnqwIAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMCZBvcV7wAuMk/Tfj/iibDTxdsEjZU6DZIJuh4BmrGm/QoEAAAaNIIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZwgiAADAGbrvAufCYodaW91gPZGRVuooyt7LgMfrtVMo0uJ8W/rdGVsdaqvtdbo1lZV2ClVVWSkTrLRTRxIdgZsR9ogAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZ+oURGbPnq3rrrtOCQkJSklJ0e23367PP/88bJ3s7Gx5PJ6wy8iRI60OGgAANA11CiL5+fkaN26cNmzYoFWrVqmqqko5OTkqLS0NW2/06NE6dOhQ6DJv3jyrgwYAAE1DnZpMrFy5Muz6/PnzlZKSos2bN6tv376h5XFxcQoEAnZGCAAAmqwL6nZVWFgoSUpKSgpbvmjRIv3hD39QamqqcnNzNX36dCUkJNRao7y8XOXl5aHrRUVFFzIkIJylhmcR0dFW6kiSJz7WTp0W8VbqmIQ4K3UkKdgixkqd6hiL/Tjt9BiUgsZKmcjyait1JCmiqMxKHc8JS6+7hfZev813/i9cUJ1qe/ON+nHef+3GGE2ePFl9+vRRjx49QsvvuecedezYUYFAQNu2bdO0adO0ZcsWrVq1qtY6s2fP1syZM893GAAAoBHzGGPOK+aPGzdOK1as0Lp169S2bdszrrd582Zde+212rx5s66++urTbq9tj0i7du2UrdsU5bHUUhzNF3tEzoo9IueIPSJnFWSPSLNWZSq1VstVWFioli1bnvP9zuuvfcKECXrrrbf03nvv/WAIkaSrr75aXq9XO3furDWI+Hw++Xy+8xkGAABo5OoURIwxmjBhgt58802tXbtWHTt2POt9tm/frsrKSqWlpZ33IAEAQNNUpyAybtw4LV68WMuXL1dCQoIKCgokSX6/X7Gxsfriiy+0aNEiDRo0SMnJyfr00081ZcoUXXXVVbrhhhvq5QkAAIDGq04H0F944QUVFhYqOztbaWlpocvSpUslSdHR0frLX/6iAQMG6LLLLtPPf/5z5eTk6N1331VkZGS9PAEAANB41fnQzA9p166d8vPzL2hAAACg+aDXDAAAcIYgAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcsdjQAbDIVo+YGDvtAyISW1mpI0nBtNZW6pReYqdHTGmqvZeBsmQ7dSpb2KkjScbSVxhFVNmpE1Vqp44kxR6xM1H+L2vvjl5Xvi/t9QcLfn3USh1zyk4/HpmgnTo4DXtEAACAMwQRAADgDEEEAAA4QxABAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzhBEAACAM3TfRYPkibTTMjWiRbyVOtVtLbWVlXSsu52OqScu81ipY9qfslJHktq1OWalTkbCcSt1JCk+qtxKnYqgnZfLr076rdSRpM/2B6zUKU+MsVInRSlW6kiSr9JSu+MqO3WCFRVW6uB07BEBAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzhBEAACAMwQRAADgDEEEAAA4Q9M7NEieCDsN3eTzWSlT3tpOUzBJKu5g57nFX26nwdxNbXdaqSNJ17XYbaVOB+/XVupIUuuIMit1fB5jpc7XQTvbpCStSLrSSp2F6mWlTuEJe38nbb5uaaWOp6jYSp3I6GhVl560UksmaKdOE8EeEQAAzsJaCMFpCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABn6L6Lps3Y6ZjqaYDNMqMi7QzK66m2UkeSyoJeK3WKg/a6uCZEVFip0zYy0kqd9lHRVupIklpusVLmg0sutVLnq+T2VupIUnULO9tAZBT/5ho69ogAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZ+oURGbPnq3rrrtOCQkJSklJ0e23367PP/88bJ3y8nJNmDBBycnJio+P15AhQ3TgwAGrgwYAAE1DnYJIfn6+xo0bpw0bNmjVqlWqqqpSTk6OSktLQ+tMmjRJb775ppYsWaJ169appKREgwcPVnW1va+RBgAATUOdvoR/5cqVYdfnz5+vlJQUbd68WX379lVhYaF+//vf65VXXlH//v0lSX/4wx/Url07vfvuuxowYIC9kQMAgEbvgs4RKSwslCQlJSVJkjZv3qzKykrl5OSE1klPT1ePHj20fv36WmuUl5erqKgo7AIAAJqH825LaIzR5MmT1adPH/Xo0UOSVFBQoOjoaCUmJoatm5qaqoKCglrrzJ49WzNnzjzfYaCJMpYO5ZlTp6zU8X1jp44kxR+w06H2WGKSlTorKrpbqSNJH7bIsFKndexJK3UkqWuL2l976mpgy0+s1LnWV2mljiS1stRZOMlnZ/ve77NSRpJkojx2CkXYqeOxVEeSDGcqhDnvPSLjx4/XJ598oj/+8Y9nXdcYI4+n9l/itGnTVFhYGLrs37//fIcEAAAamfMKIhMmTNBbb72lNWvWqG3btqHlgUBAFRUVOn78eNj6R44cUWpqaq21fD6fWrZsGXYBAADNQ52CiDFG48eP1xtvvKHVq1erY8eOYbdfc8018nq9WrVqVWjZoUOHtG3bNvXu3dvOiAEAQJNRp3NExo0bp8WLF2v58uVKSEgInffh9/sVGxsrv9+vn/70p5oyZYpat26tpKQkTZ06VZmZmaFP0QAAANSoUxB54YUXJEnZ2dlhy+fPn69Ro0ZJkp577jlFRUVp+PDhOnXqlG6++WYtWLBAkZGRVgYMAACajjoFEWPMWdeJiYnRb3/7W/32t78970EBAIDmgV4zAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHDmvJveAfXJBM/+UfFzqnOqzEqdiMPHz77SOUr8W7SVOhHVdjqMnTzkt1JHkr5qYafWnpZBK3Uk6bO2KVbq+DpVWakTiNpkpY4knQza+X6mCkt1PBabuXmq7bwGyNJrSURCgoLFxVZqIRx7RAAAOAtCSP0hiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJyh+y4aJmOn+6qptNMx1Zw8aaWOJEUXFFmp41dLK3V8hXa6AUtSWaKd9zalafbeI51KtPP8iqpjrdQpDtp72T0RtDOm42V26kTaaXYtSfJUWurAbCx1323RQtVFNL6rD+wRAQDgLAgh9YcgAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIbuu2iYPHYyssdrZxP3xMRYqSNJ1Ql2Op1WtPJaqXOqtb33IycDduqUZVTaKSSpe9sCK3Uy4/ZbqeP1WOoqK+mrykQrdY4Ut7BSx1dqpYwkKaKswkodU2lvW0L9YI8IAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZwgiAADAGZreoUGKsNSsLiLBTjOvqvYpVupI0vHL46zUKezksVKn8hI7zcUkqU2bQit1bkrdZ6WOJPX3b7dSJ9N32EqdSmPn9yZJX5SnWqlz8oSdRowJRcZKHUnylNlpVmeMnTFFxPgUPHXKSi2EY48IAABnQQipPwQRAADgDEEEAAA4QxABAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4Eydg8h7772nW2+9Venp6fJ4PFq2bFnY7aNGjZLH4wm79OrVy9qAAQBA01HnIFJaWqorr7xSc+bMOeM6AwcO1KFDh0KXt99++4IGCQAAmqY6dxbLzc1Vbm7uD67j8/kUCATOe1AAAKB5qJfuu2vXrlVKSopatWqlrKwsPfXUU0pJqb17aXl5ucrLy0PXi4qK6mNIuBg8Fk85ioy0U8dS992Tl9jpTipJx7vaqdO6x9dW6twQ+NJKHUnKjDtgpU5331dW6kjSj7x2urjGeXxW6myrtNehdmepna7QkUe9Vur4iqqt1JEklZWffZ1zUW1nTBHR0QpW2OtUjb+zfrJqbm6uFi1apNWrV+vXv/61Nm7cqJtuuiksbHzX7Nmz5ff7Q5d27drZHhIAABeEEFJ/rO8RGTFiROjnHj166Nprr1VGRoZWrFihoUOHnrb+tGnTNHny5ND1oqIiwggAAM1EvRya+a60tDRlZGRo586dtd7u8/nk89nZ5QkAABqXev8ekaNHj2r//v1KS0ur74cCAACNTJ33iJSUlGjXrl2h67t379bHH3+spKQkJSUlacaMGbrzzjuVlpamPXv26JFHHlFycrLuuOMOqwMHAACNX52DyKZNm9SvX7/Q9ZrzO/Ly8vTCCy9o69atevnll3XixAmlpaWpX79+Wrp0qRISEuyNGgAANAl1DiLZ2dky5swfP3vnnXcuaEAAAKD5oNcMAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJyp9694B86Lx2OnToSdrB2005xUkmSi7XRfjfPaacLli6iyUkeSyoydiSqo9lupI0lezzErdZIiyqzU+bq6tZU6kvRViZ15ij5h5+/NW2in07Ek6QyNUuvKBO11O0b9YI8IAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZwgiAADAGYIIAABwhu67sMcE7dWqstQRtvSklTJxBXY63UpSi90xVurs9aTZqZNkrxusL9ZO99WEWDudbiWpS9LXVurc0GqXlTqRsvd3Ul5t5yU8wtLmHVlmr5OzqbRUy+brEuoFe0QAAIAzBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzhBEAACAMwQRAADgDEEEAAA4QxABAADO0PQO9ngaXq41xSVW6kTvO2aljiSlVCdaqZNwwGelTkWCnTqSVO2z09DvpD/BSh1J+qBjkpU6J7rFWqlzdeJ+K3Ws8rgeQC2MsVMmaKeOPBE00KsnDe8/BwAADQ0hpN4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzhBEAACAMwQRAADgDEEEAAA4QxABAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM7QfbexstjpNsJraTMx18stAAAQKElEQVSIsrc5ebxea7VsMMdOWKvlPVVmp85XdjrdGl+0lTqSFLTUyfdkepyVOpJUHRNppc6R9i2s1Clpaa/bsTENrG0ufeFwHtgjAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJypcxB57733dOuttyo9PV0ej0fLli0Lu90YoxkzZig9PV2xsbHKzs7W9u3brQ0YAAA0HXUOIqWlpbryyis1Z86cWm9/5pln9Oyzz2rOnDnauHGjAoGAbrnlFhUXF1/wYAEAQNNS5+Ygubm5ys3NrfU2Y4x+85vf6Be/+IWGDh0qSVq4cKFSU1O1ePFijRkz5sJGCwAAmhSr54js3r1bBQUFysnJCS3z+XzKysrS+vXra71PeXm5ioqKwi4AAKB5sNp9t6CgQJKUmpoatjw1NVV79+6t9T6zZ8/WzJkzbQ6jYbPUNTci2l7H1Ah/gp1CiX47dSQFY+w8v4jySit1VFxqp44kVVTYqVNcYqWMp8JeN9gIr51OtybSXlfZyhbGSp3keDvbgNdTbaWOJJ2ssNOlOspOQ2h5Ku09N1VbrIUGrV4+NePxhL+IGGNOW1Zj2rRpKiwsDF32799fH0MCAAANkNU9IoFAQNK3e0bS0tJCy48cOXLaXpIaPp9PPp+9d2QAAKDxsLpHpGPHjgoEAlq1alVoWUVFhfLz89W7d2+bDwUAAJqAOu8RKSkp0a5du0LXd+/erY8//lhJSUlq3769Jk2apFmzZqlz587q3LmzZs2apbi4ON19991WBw4AABq/OgeRTZs2qV+/fqHrkydPliTl5eVpwYIFeuihh3Tq1CmNHTtWx48fV8+ePfXnP/9ZCQmWTogEAABNRp2DSHZ2tow581noHo9HM2bM0IwZMy5kXAAAoBmg1wwAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnLH6Fe84O0+knaZg1hrVSaruELBSp/jSeCt1JKm8pZ2M7C210/AsvqCllTqS5D1uqcNYtZ3nFoyz0zhNkk6mxVipc6yrvfdI0Z3sdPTu7i+wUudUtb2GlUXH46zUSTlhZ1uKOFlupY4kBW01vTNBO3VQb9gjAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnKH77kXm8Vqa8pb2uu+WtLfTNfdoD3u5tiK9wk6hMjvdjmMO+6zUkaToE3ZqeSw1Fa2y08BVknQq1U4X15gOhVbqSFKftrut1Iny2OkG+9ejHazUkSTfXjvbUlyBpb+3klI7dSSZyiprtdCwsUcEAAA4QxABAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzhBEAACAMwQRAADgDEEEAAA4Q/fdc+WxlNk8Hjt1oux0lZWkyng7Y6pIttcts1PGYSt1WvnKrNT5+pSdDsWSVHgqxkodY+z83lr6LHVeldS1RZGVOm1jj1upY9O6I5daqXNoe6qVOpKU8rmlbsf77fzejM3uu9V2uh2j4WOPCAAAcIYgAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBma3p0rE7RTp9pSnTJ7jcp8hXbG5D1qb3P6Jt1Ok7mOCces1OnaosBKHUnyRdhpDhjpsbQtWXSyOtpKnX2nkqzUkaSPD19ipc6pz1pZqZPyiZ1GdZLk315opY7nyFErdapP2WkyKcneay4aPPaIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJwhiAAAAGcIIgAAwBnrQWTGjBnyeDxhl0AgYPthAABAE1AvvWa6d++ud999N3Q9MjKyPh4GAAA0cvUSRKKiotgLAgAAzqpegsjOnTuVnp4un8+nnj17atasWbr00ktrXbe8vFzl5eWh60VFRfUxpAbDVFXaqWNxnuL32ul02zrOb6WOJB2rttN9dfWlsVbqpKecsFJHklLjSqzVsuFEeYy1WgVFLa3UKT1sZ5uUpPjddl7mAn+rtlKnxU5725IOHrFSJlhSaqWOrdc3NC/WzxHp2bOnXn75Zb3zzjv63e9+p4KCAvXu3VtHj9beZnr27Nny+/2hS7t27WwPCQAANFDWg0hubq7uvPNOZWZmqn///lqxYoUkaeHChbWuP23aNBUWFoYu+/fvtz0kAADQQNXLoZnvio+PV2Zmpnbu3Fnr7T6fTz6fr76HAQAAGqB6/x6R8vJy7dixQ2lpafX9UAAAoJGxHkSmTp2q/Px87d69W//7v/+rYcOGqaioSHl5ebYfCgAANHLWD80cOHBAd911l7755hu1adNGvXr10oYNG5SRkWH7oQAAQCNnPYgsWbLEdkkAANBE0WsGAAA4QxABAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM7U+1e8I5ypttPB01jqlilJEQfsdPBsVWnnuUmS77idLq7FB+x03z2WbKeOJB2Js1bKishT9mr5LDWWTfsmaKeQpPiDdp6g96CdJ2e+OWaljiQFS09aqUPXXLjEHhEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzhBEAACAMwQRAADgDEEEAAA4QxABAADOEEQAAIAzNL1rpIKVVdZqmeOFVupElJVbqSNJcceLrdSJ3dfCSp3qBHtN74K+SDuFjLFSJqLCXoO5yFI724Cn2GInvlI7DSKDlhpNmooKK3Uke000AZfYIwIAAJwhiAAAAGcIIgAAwBmCCAAAcIYgAgAAnCGIAAAAZwgiAADAGYIIAABwhiACAACcIYgAAABnCCIAAMAZgggAAHCGIAIAAJyh+25jZex1TDVVdmpVl9jrBBpRbqmT7/ETVspEer1W6khSZGQDy//VFrelyko7daosdpe21KHWBO10O7b5tws0BQ3sFREAADQnBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzhBEAACAMwQRAADgDEEEAAA4QxABAADONLimd8Z821iqSpWSpR5TuFjs5doI07AajHlsjUeSgg0s/9tsoGjsNL2TsddA0ViqZRrYNgk0NFX69u+/rn8rDS6IFBcXS5LW6W3HI0Gd2QyOFQ2sDgDgnBQXF8vv95/z+h5jLebbEQwGdfDgQSUkJMjj8bgeToNVVFSkdu3aaf/+/WrZsqXr4TR5zPfFxXxfXMz3xdOU59oYo+LiYqWnpysi4tz3/Da4PSIRERFq27at62E0Gi1btmxyG3NDxnxfXMz3xcV8XzxNda7rsiekRgM7WA0AAJoTgggAAHAmcsaMGTNcDwLnJzIyUtnZ2YqKanBH2Jok5vviYr4vLub74mGuwzW4k1UBAEDzwaEZAADgDEEEAAA4QxABAADOEEQAAIAzBBEAAOAMQaSRmTFjhjweT9glEAi4HlaT8d577+nWW29Venq6PB6Pli1bFna7MUYzZsxQenq6YmNjlZ2dre3btzsabeN2trkeNWrUadt6r169HI228Zs9e7auu+46JSQkKCUlRbfffrs+//zzsHXKy8s1YcIEJScnKz4+XkOGDNGBAwccjbhxO5f5zs7OPm0bHzlypKMRu0MQaYS6d++uQ4cOhS5bt251PaQmo7S0VFdeeaXmzJlT6+3PPPOMnn32Wc2ZM0cbN25UIBDQLbfcEmrWiHN3trmWpIEDB4Zt62+/TTPM85Wfn69x48Zpw4YNWrVqlaqqqpSTk6PS0tLQOpMmTdKbb76pJUuWaN26dSopKdHgwYNVXW2vG3JzcS7zLUmjR48O28bnzZvnaMQOGTQq06dPN1deeaXrYTQLksybb74Zuh4MBk0gEDD/9m//FlpWVlZm/H6/mTt3roshNhnfn2tjjMnLyzO33XaboxE1fUeOHDGSTH5+vjHGmBMnThiv12uWLFkSWuerr74yERERZuXKla6G2WR8f76NMSYrK8tMnDjR4agaBvaINEI7d+5Uenq6OnbsqJEjR+rLL790PaRmYffu3SooKFBOTk5omc/nU1ZWltavX+9wZE3X2rVrlZKSoi5dumj06NE6cuSI6yE1GYWFhZKkpKQkSdLmzZtVWVkZtn2np6erR48ebN8WfH++ayxatEjJycnq3r27pk6d2iz3rvL9so1Mz5499fLLL6tLly46fPiwnnzySfXu3Vvbt29X69atXQ+vSSsoKJAkpaamhi1PTU3V3r17XQypScvNzdVPfvITZWRkaPfu3Xrsscd00003afPmzfL5fK6H16gZYzR58mT16dNHPXr0kPTt9h0dHa3ExMSwdVNTU0PbPs5PbfMtSffcc486duyoQCCgbdu2adq0adqyZYtWrVrlcLQXH0GkkcnNzQ39nJmZqeuvv16dOnXSwoULNXnyZIcjaz48Hk/YdWPMactw4UaMGBH6uUePHrr22muVkZGhFStWaOjQoQ5H1viNHz9en3zyidatW3fWddm+L9yZ5nv06NGhn3v06KHOnTvr2muv1UcffaSrr776Yg/TGQ7NNHLx8fHKzMzUzp07XQ+lyav5dNL33x0eOXLktL0ksC8tLU0ZGRls6xdowoQJeuutt7RmzRq1bds2tDwQCKiiokLHjx8PW5/t+8Kcab5rc/XVV8vr9Ta7bZwg0siVl5drx44dSktLcz2UJq9mF+p3d5tWVFQoPz9fvXv3djiy5uHo0aPav38/2/p5MsZo/PjxeuONN7R69Wp17Ngx7PZrrrlGXq83bPs+dOiQtm3bxvZ9Hs4237XZvn27Kisrm902zqGZRmbq1Km69dZb1b59ex05ckRPPvmkioqKlJeX53poTUJJSYl27doVur579259/PHHSkpKUvv27TVp0iTNmjVLnTt3VufOnTVr1izFxcXp7rvvdjjqxumH5jopKUkzZszQnXfeqbS0NO3Zs0ePPPKIkpOTdccddzgcdeM1btw4LV68WMuXL1dCQkJoz57f71dsbKz8fr9++tOfasqUKWrdurWSkpI0depUZWZmqn///o5H3/icbb6/+OILLVq0SIMGDVJycrI+/fRTTZkyRVdddZVuuOEGx6O/yFx+ZAd1N2LECJOWlma8Xq9JT083Q4cONdu3b3c9rCZjzZo1RtJpl7y8PGPMtx/hnT59ugkEAsbn85m+ffuarVu3uh10I/VDc33y5EmTk5Nj2rRpY7xer2nfvr3Jy8sz+/btcz3sRqu2uZZk5s+fH1rn1KlTZvz48SYpKcnExsaawYMHM+fn6WzzvW/fPtO3b1+TlJRkoqOjTadOnczPf/5zc/ToUbcDd8BjjDEXM/gAAADU4BwRAADgDEEEAAA4QxABAADOEEQAAIAzBBEAAOAMQQQAADhDEAEAAM4QRAAAgDMEEQAA4AxBBAAAOEMQAQAAzvw/Ud5WCjL7HaQAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f7ace4bedd8>)"
      ]
     },
     "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** is 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
}