{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Regularized OT with generic solver\n",
    "\n",
    "\n",
    "Illustrates the use of the generic solver for regularized OT with\n",
    "user-designed regularization term. It uses Conditional gradient as in [6] and\n",
    "generalized Conditional Gradient as proposed in [5][7].\n",
    "\n",
    "\n",
    "[5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, Optimal Transport for\n",
    "Domain Adaptation, in IEEE Transactions on Pattern Analysis and Machine\n",
    "Intelligence , vol.PP, no.99, pp.1-1.\n",
    "\n",
    "[6] Ferradans, S., Papadakis, N., Peyré, G., & Aujol, J. F. (2014).\n",
    "Regularized discrete optimal transport. SIAM Journal on Imaging Sciences,\n",
    "7(3), 1853-1882.\n",
    "\n",
    "[7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). Generalized\n",
    "conditional gradient: analysis of convergence and applications.\n",
    "arXiv preprint arXiv:1510.06567.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "import ot.plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data\n",
    "-------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%% parameters\n",
    "\n",
    "n = 100  # nb bins\n",
    "\n",
    "# bin positions\n",
    "x = np.arange(n, dtype=np.float64)\n",
    "\n",
    "# Gaussian distributions\n",
    "a = ot.datasets.get_1D_gauss(n, m=20, s=5)  # m= mean, s= std\n",
    "b = ot.datasets.get_1D_gauss(n, m=60, s=10)\n",
    "\n",
    "# loss matrix\n",
    "M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))\n",
    "M /= M.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmcVNWZ//HPQzeL7CDNvrQsxhhBxY6CGsaA4hINGXGUDI6MG+OKg8aoITFxMnH5uYcEFZcEE1xxDW4/g7hFQGkRFZVFQARlFUQQge4+88dzO91CA71U31vL9/163Vd13bpd9XQ1/eXUueeeYyEERESk/jVIugARkVyhwBURiYkCV0QkJgpcEZGYKHBFRGKiwBURiYkCVyRLmNkPzGx+0nXIrilwRSJm9p9m9p6ZfW1mK83sDjNrHT12p5ltirZtZra90v3nYqgtmFnv3R0TQngthPCdOrzGCDObZWabzWx19PUFZmbR42ZmN5jZumi7ofwxqR4FrghgZpcBNwCXA62AAUAP4EUzaxRCOC+E0DyE0By4Fni4/H4I4fjkKndmll/H778MuB24EegIdADOA44AGkWHjQZ+AhwI9ANOAv6rLq+baxS4kvPMrCVwDXBxCOH5EML2EMJS4FSgEDi9Fs95lJktN7OfR63Fz83sJ2Z2gpktMLMvzOwXlY4/1MxmmNmG6Ng/mFmj6LFXo8PmRi3q0yo9/xVmthL4U/m+6Ht6Ra/RP7rf2czWmNlRVdTaCvgf4IIQwpQQwlfBzQkhjAwhbI0OHQXcHEJYHkJYAdwM/GdN35tcpsAVgcOBJsDjlXeGEDYBzwLH1PJ5O0bP2wW4GrgbD+9DgB8AvzKzfaJjS4GxQDtgIDAEuCCqY1B0zIFRi/rhSs/fFm+Jj96h9o+BK4C/mllT4E/ApBDCy1XUORBoDDy1h5/ne8DcSvfnRvukmhS4Ih5ya0MIJVU89nn0eG1sB34XQtgOPBQ9z+1RC3Ie8AH+8ZwQQnEIYWYIoSRqXd8F/Msenr8M+HUIYWsIYcuOD4YQ7gYWAbOATsC4XTzPTj+/mb0Rtba3mFl54DcHvqz0fV8CzdWPW30KXBFYC7TbRT9op+jx2lgXQiiNvi4PxFWVHt+Chxhmtq+ZTY1O1m3E+4n3FPRrQgjf7OGYu4EDgPGVugZ2qpMdfv4QwuEhhNbRY+U5sQloWen7WgKbgmbAqjYFrgjMALYCJ1feaWbNgeOBaTHUcAfwEdAnhNAS+AWwp5bjboMuqv824F7gN2bWdheHlv/8w/bwevOIWuSRA6N9Uk0KXMl5IYQv8ZNm483sODNraGaFwCPAcuAvMZTRAtgIbDKz/YDzd3h8FdCzhs95OzA7hHAO8AxwZ1UHhRA24D//BDM7xcxamFkDMzsIaFbp0PuBS82si5l1Bi4D/lzDmnJanYaSiGSLEML/M7N1wE1ALzz8ngRG7uajeCr9DJgI/ByYAzwMDK70+G+ASWa2F36CbPXunszMhgHHAX2jXZcC75jZyBDC5B2Pj37+FdHr3w9sBhbjJ97eiA67Cw/996L790T7pJpM3S8iIvFQl4KISEwUuCIiMVHgiojERIErIhITjVIQANq1axcKCwuTLkMkoxQXF68NIRRU93gFrgBQWFjI7Nmzky5DJKOY2Sc1OV5dCiIiMVHgiuSKLVvgrbdg0SLQ+PtEKHBFst2cOTB4MLRsCYceCn36QNu2cOml8NVXSVeXUxS4ItmqrAyuvBKKimDePLj8cnjsMbj7bjj+eLjtNthvP5g+PelKc4ZOmolko7IyOP98mDgRzj4bbrwR2rSpePycc+CSS+DMM+GEE2DqVBgyJLl6c4RauCLZ6OKLPWx/8Qtv0VYO23KHHQavvOJdDCeeCK++uvMxklIKXJFsc//9MGECXHYZ/O//wu4WZCgogGnToHt3OO00WL3bScikjhS4Itlk/ny44AIYNAiuv373YVuuoAAefRTWr4dRo7w7QuqFAlckW5SWwsiR0KQJTJ4M+TU4RdOvH9x6Kzz/PPzxj/VXY45T4Ipki4kTobjYA7Nr15p//3nnwdCh8MtfwqpVez5eakyBK5IN1q6FceN8vO2pp9buOcxg/Hi/QOLKK1NbnwAKXJHsMG6cX8Qwfnz1+m13Zd99/WTbn/8MM2akrDxxClyRTLdgAdxzj58s23//uj/fuHHQsSNcdZUuAU4xBa5IprvmGj9RNm5cap6veXMfv/vKK/DSS6l5TgEUuCKZbd48ePBBv9ChffvUPe+55/qJt1/+Uq3cFFLgimSya67xFunll6f2eZs08bCdOdOHiklKKHBFMtXChTBlClx0Eey9d+qf/8wzoVs3v4BCUkKBK5Kpbr4ZGjWCMWPq5/kbNYKxY32OhVmz6uc1cowCVyQTrVrlQ7dGjfIRBfXlnHOgdWufbUzqTIErkonGj4dt23zMbH1q0cKHmz3+uHdhSJ0ocEUyzZYtcMcdMGyYX6hQ3y6+GBo2hN//vv5fK8spcEUyzQMPwBdf+ATicejYEUaM8C6MjRvjec0spcAVySQheHdC377wL/8S3+tefDFs2uShK7WmwBXJJK+9BnPnegDWZc6EmioqgoEDPew1X26tKXBFMsn48b5czsiR8b/2mDG+xLouhKg1Ba5IpvjsM3jiCTjrLGjaNP7XHz7c+3PvuCP+184SClyRTHHvvb6qw3/9VzKv37ChrwD8zDPwySfJ1JDhFLgimaCkxFd0OOYYX2U3KaNHe9/x3XcnV0MGU+CKZIJnn4Xly+H885Oto3t3+NGPfP7dbduSrSUDKXBFMsGdd0LnznDiiUlX4mufrVoFTz6ZdCUZR4Erku6WLvWRAeec4/2oSTv2WOjRw7s4pEYUuCLp7p57vN/07LOTrsTl5Xn4T5vmw8Sk2hS4Iuls+3a47z44/njvP00XZ53lwauTZzWiwBVJZ888A59/ntxQsF3p3BlOOgn+9CedPKsBBa5IOrvrLujSxVu46Wb0aFizRifPakCBK5Kuli6FF17w/tL8/KSr2dnQoTp5VkMKXJF0VX6y7Kyzkq6kajp5VmMKXJF0lK4ny3akk2c1osAVSUdTp6bnybId6eRZjShwRdJROp8s25FOnlWbAlck3Sxe7CfLzj03PU+W7WjoUCgs1LSN1aDAFUk3d91VcUIqE+TledfHyy/DRx8lXU1aU+CKpJOtW/1k2Y9/7F0KmeKss3yehzvvTLqStKbAFUknjz8Oa9cmPw1jTbVv7ytCTJoEX3+ddDVpS4Erkk4mTIBevWDIkKQrqbnzzoMNG+DBB5OuJG0pcEXSxTvvwOuvw4UXQoMM/NMcNAgOOAD+8Adfzl12koG/VZEs9Yc/+OKQZ56ZdCW1YwYXXeT/cbzxRtLVpCUFrkg6WLcOJk+G//gPaN066Wpq7/TToVUr/89DdqLAFUkH990H33zj3QmZrFkzH7EwZYpfKSffosAVSVpJCYwfD0cdBX37Jl1N3V14oS/nPmFC0pWkHQWuSNKmTIFPP4VLL026ktTo1QuGDfMrzzRE7FsUuCJJCgFuvhn69PHlx7PFZZd5v/T99yddSVpR4Iok6R//gNmzYezYzBwKtitHHAHf/z7ceiuUlSVdTdrIot+wSAa68UZo2xbOOCPpSlLLzLtIFiyAp59Oupq0ocAVScp773kYjRnjZ/ezzSmnQM+ecN11uhAiosAVScp110Hz5nDxxUlXUj/y8+GKK+DNN30ZHlHgiiRi0SJ4+GGfpKZt26SrqT+jRvmqENdem3QlaUGBK5KE667z6QzHjk26kvrVuDH87GcwfbrPE5HjFLgicVuwwKcxPO886NQp6Wrq3+jR0KED/PKXOd+Xq8AVidtvfuMtv6uuSrqSeDRrBuPGwSuv5HxfrgJXJE7vvQcPPQSXXOKtvlwxerQv9z5uXE63chW4InG64gpo2dL7NXNJ48Zw9dU+YuGxx5KuJjEKXJG4PPecb7/6VXaPTNiVUaN8cp7LL/eZ0XKQAlckDtu3+5VXvXtn77jbPcnPh9tug6VL/ZLfHKTAFYnDH//oS4jffDM0apR0NckZPNhnEvvd72DFiqSriZ0CV6S+LVvmQ6KOOw5OOinpapJ3yy0+oc2FF+bcCTQFrkh9CqFiyfM77/RJXXJdz55wzTXw1FO+LHwOUeCK1KcHH4Rnn/WP0D16JF1N+hg7Fg4+2Bed/OKLpKuJjQJXpL4sWeKt24EDPVikQn6+r+O2bh2cc07OdC0ocEXqQ0kJjBzpXz/wAOTlJVtPOjroIJ9T4okn4K67kq4mFgpckfpw1VUwY4YHSWFh0tWkr7Fj4dhj/ba4OOlq6p0CVyTVJk2Cm27y7oQRI5KuJr01aODrnrVv78PFsnxpdQWuSCq99prPGzBkCNx+e9LVZIb27X3liw0bPHQ3bUq6onqjwBVJlTff9JV399kHHnnE57uV6jnwQJg82bsVfvxj2LIl6YrqhQJXJBVmz/a+yIICn4IwF+dKqKthw7w75uWX4eST4euvk64o5RS4InX13HNw1FHQqpWHbZcuSVeUuU4/He6+G154wbtl1q5NuqKUUuCK1FYIPhnLSSfBvvv6qASNSKi7s8+GKVPgnXdgwAC/zRIKXJHaWL3aPwKPHQsnnuirGeTCcjlxOflkeOkl78sdMADGj/f5FzKcAlekJkpL/SPvfvvB88/7SIQnnoAWLZKuLPsMHOit28GDYcwYOPJImDs36arqRIErUh2lpT7yoF8/H/bVr5//8Y8Zowlp6lNBATzzDPz5z7Bwoc+/8O//Dh9+mHRltaLAFdmdJUt84pmePeG003zfI4/4st/f/W6yteUKM18tYv58X6Lo6adh//3hmGP8d7F5c9IVVpuFHJk0QnavqKgozJ49O+kykrdxo4+nnT7duwzeftv3DxkCF1zg/baaFyFZa9bAxIm+LVsGe+0FQ4fC0UfDoEEexvn5sZRiZsUhhKJqH6/AFcihwN261acDXLUKVq70P9glS7z1NG8eLFjgx+Xnw6GHwr/+Kwwf7hczSHopLfUr+6ZM8Skwlyzx/Xvt5V0+3/0u9Onj02J26+arJBcU+PC9FP2nqcCVWvlW4D75JEydWv8vuqt/e+X7K9/uaist9bPXJSX+9fbtvm3d6tuWLb599ZVvVV3BlJ/va41997vQvz98//tw+OE6EZZpliyBN96At97y5eg//LDquRnM/HfbogU0awZNm0KTJr6ycKNGfoVgfr5veXk+38OgQf4JZ6enqlngxtPulsyyaJF/nI7Drk44le+vfFvVVv4HUf7H0bChb40bQ/Pm3trZa6+KP7A2bfwqsPbtvcXTvTt07qxugmywzz6+lU+LCd6/u2yZr5+2apVfSPHFF/Dll/4f8ObNfkXbN9/4f9BffeX/YZeU+FZW5luKJo9XC1eAHOpSEEmhmrZwNUpBRCQmClwRkZioS0EAMLM1wCeVdrUD0nnmkHSuL51rg/SuL51rg53r6xFCKKjuNytwpUpmNrsmfVNxS+f60rk2SO/60rk2qHt96lIQEYmJAldEJCYKXNmViUkXsAfpXF861wbpXV861wZ1rE99uCIiMVELV0QkJgpcEZGYKHBlJ2Z2nJnNN7NFZnZlwrV0M7PpZvaBmc0zs0ui/W3N7EUzWxjdtkmwxjwzm2NmU6P7+5jZrOj9e9jMGiVYW2szm2JmH5nZh2Y2MM3eu7HR7/V9M3vQzJok+f6Z2X1mttrM3q+0r8r3y9zvozrfNbP+e3p+Ba58i5nlAX8Ejgf2B35qZvsnWFIJcFkIYX9gAHBhVM+VwLQQQh9gWnQ/KZcAlZcguAG4NYTQG1gPnJ1IVe524PkQwn7AgXidafHemVkXYAxQFEI4AMgDRpDs+/dn4Lgd9u3q/Toe6BNto4E79vjsIQRt2v65AQOBFyrdvwq4Kum6KtXzFHAMMB/oFO3rBMxPqJ6u0R/hYGAqYPiVSPlVvZ8x19YKWEJ0crzS/nR577oAnwJt8ZkLpwLHJv3+AYXA+3t6v4C7gJ9WddyuNrVwZUflfwTllkf7EmdmhcDBwCygQwihfLLTlUCHhMq6Dfg5UL6k7N7AhhBCSXQ/yfdvH2AN8Keoy+MeM2tGmrx3IYQVwE3AMuBz4EugmPR5/8rt6v2q8d+KAlcygpk1Bx4D/juEsLHyY8GbF7GPbzSzE4HVIYTiuF+7mvKB/sAdIYSDgc3s0H2Q1HsHEPWFDsP/Y+gMNGPnj/Nppa7vlwJXdrQC6FbpftdoX2LMrCEetpNDCI9Hu1eZWafo8U7A6gRKOwL4sZktBR7CuxVuB1qbWfnk/km+f8uB5SGEWdH9KXgAp8N7B3A0sCSEsCaEsB14HH9P0+X9K7er96vGfysKXNnRW0Cf6ExxI/wkxtNJFWNmBtwLfBhCuKXSQ08Do6KvR+F9u7EKIVwVQugaQijE36eXQggjgenAKUnWFtW3EvjUzL4T7RoCfEAavHeRZcAAM2sa/Z7L60uL96+SXb1fTwNnRKMVBgBfVup6qFoSneXa0nsDTgAWAB8D4xKu5Uj8I9y7wDvRdgLeVzoNWAj8HWibcJ1HAVOjr3sCbwKLgEeBxgnWdRAwO3r/ngTapNN7B1wDfAS8D/wFaJzk+wc8iPcnb8c/IZy9q/cLP0H6x+jv5D18tMVun79Ol/aa2XH4R6g84J4QwvW1fjIRkSxX68CNxmsuwIfoLMc/iv40hPBB6soTEckedVm191BgUQhhMYCZPYSfcdxl4LZr1y4UFhbW4SWlrjZs8IVJu3X79v7i4uK1IZq5/pgG/6YZjUR28GLZo7tYYrr66hK4VY1BO2zHg8xsNH4VBt27d0crwybr5z+H8eNhx1+DmX1S9XeISKrU+yiFEMLEEEJRCKGooKDaS/9IPdm+HRo2TLoKkdxUl8BNu/GasmfbtkGjxKZSEcltdQnctBqvKdXz9dfQrFnSVYjkplr34YYQSszsIuAFfFjYfSGEeSmrTOrFxo3QvHnSVYjkprqcNCOE8CzwbIpqkRisWwd77510FSK5SZf25pgVK6Bjx6SrEMlNCtwcUloKS5dCr15JV5IgM99EEqDAzSHvvw8lJXDAAUlXIpKb6tSHK5nllVf8dtCgZOtIVPml7OWtXIvaHKHs24+L1AO1cHPIX/8K+++/82W9IhIPtXBzxMyZ8NZbflmvUNGSDaV+2yAPAGvkfxJhe7TCS1lp3JVJFlMLNwds2wbnnw8dOsAZZyRdjUjuUgs3y4UA48bBO+/Ak09Cy5ZJV5SmopZs2Oq3lu9/GraXX5YXtm332+3bEihOsoVauFksBJ8d7KabvIU7bFjSFYnkNrVws9SXX8KYMXD//XDhhfD73yddUWYJJSXfum3QpInftvLL9MLXWwAo+/rrBKqTTKUWbhZ68kkfjfDXv8LVV/uJsgb6TYskTi3cLDJjBlx3Hfztb9Cvnwfv97+fdFXZoeybb/yL6DYv6gzP7+Fj7MLGTQCUrl8ff3GSMdTuyXClpfDEE3DEEXD44fD663Dttb6ig8JWJL2ohZuhPvwQHn4YJk+GRYugsND7ac88U9MvxqF040b/IrrNi1YzadBvP79dswGAks9Xxl+cpC0FbgZZsMBD9pFHfF4EM79M93e/g5NPhnz9NkXSmv5E09imTfDaa/D3v8OLL8J773nIHnmknwgbPhw6dUq6SgEoXbPGv4huGxR29/0/7A9Ao+Xe4i1duDj+4iRtKHDTyPbtMGuWB+y0aX45bkmJr0F2xBFw221wyinQpUvSlYpIbShwE7R6tQds+TZjBmze7K3YQw6Byy6Do4/2sN1rr6SrlZooWboMgLzolu99B4CvRgwAoOXCrwAIxVqVKpcocGOydatfXjtzpofrzJmwZIk/lpcHffvCqFEwZAgcdRS0bZtouSJSDxS49WDTJnj3XQ/YOXP89t13fRIZ8C6BAQP8ctvDDvPWrFbSzW6l8+YD0CJq0JYdfiAAq8YeDkC7uVsByH+pOP7iJDYK3DpaterbwTpnDixcWDH7X9u2cPDBcMklHq6HHQZduyZbs4gkQ4FbTVu2wAcf+EiBytvKSsMsCws9XEeOhIMO8q+7dtUSWrIze2MuAB3f8PtbT/CrVD65cSAAnWb4ChRNH58Vf3FSbxS4OygthcWLdw7WRYugLFqFpUkTn6vg2GMrgvXAA6F162RrF5H0ltOBu3r1zsE6bx6UTwBl5ivc9u0LI0b4bd++0Lu3n+gSSZXGz74FQK9n/f7Gf/fRDEse6gdAq//vnfxt75sRf3GSMjkRuF9/7UG6Y7iuXl1xTEGBh+m551YE6/e+p5NZIpI6WRe4X3zhJ67efrvidsGCipNYTZr4MuE/+lFFsPbt68vPiKSLlg/MjG79/uoLfDRDs1d9zoaPn+gDQMdb34i/OKm1jA3cEOCzz3YO12XLKo7p1g369/92d0CvXuoOEJFkZEzglpV5t8Crr/r8Aq++Cp9/XvH4vvvCwIG+usHBB/vWrl1y9YqkUvsJ3pLdPMHvl17pLdwT5/n8u3c89CMAuv1WLd50lraBu307FBdXhOs//gHlczt36eJXYw0Y4C3YAw+EFi0SLVdEZI/SKnBD8Am077/fpyAsn3J03319+sEf/MCnIyws1NhWyW1drveW7NTr2wBQ+ls/SXHDEh+3e+rD/w3APldqVEM62WPgmlk34H6gAxCAiSGE282sLfAwUAgsBU4NIdRqfZHFi+Evf/GgXbzYRwYMHw4nneRTEXbsWJtnFRFJL9Vp4ZYAl4UQ3jazFkCxmb0I/CcwLYRwvZldCVwJXFHTAm6+GX72M2+xDhkCv/mNt2Y1HEuk+gp/5S3ZK351GABlN/r+x5b7aIeDohZvr8tmxl+c/NMeAzeE8DnwefT1V2b2IdAFGAYcFR02CXiZGgbuxIketsOHw623+qgCEZFsVaM+XDMrBA4GZgEdojAGWIl3OVTbkiVw3nlQVAQPPOCTbItIavS63Fu8wy/3K9bCrb7/hc/e8ccfOg+A3peqxRunaq/aa2bNgceA/w4hbKz8WAgh4P27VX3faDObbWaz15QvQ4JP6vLDH/r42WefrV3xIiKZpFqBa2YN8bCdHEJ4PNq9ysw6RY93AlZX9b0hhIkhhKIQQlFBtLIpQMOG8OST3sIdPhxOOMEXSPzmmzr9PCJShd5jZ9J77EyO7XwQx3Y+CAtgwft4H1s+k49vHMjH0UxlUn/2GLhmZsC9wIchhFsqPfQ0MCr6ehTwVE1fvEULeP55uOIKn9tgxAgfkTB6tI+7LZ+dS0QkG1gIVfYEVBxgdiTwGvAeUB6Bv8D7cR8BugOf4MPCvtjdcxUVFYXZs2dX+VhpKbz8MkyaBI895hPOtGnjw8IGDfIxuP37e8tYUs/MikMIRQDHNPi33f+jkKyz5Hpv3T5y2m0AnPZXH9VQPvpB4MWyR+s8+r86oxReB3b1QkPqWkC5vDwfFjZkCEyYAE89BdOn+5Vmf/ubH9O0qV++W34BxKGHaviYiGSOPbZwU2l3LdzdWbnSr0Arn0dh7ly/Ks0MvvMdnzehf/+KORS0AGPNqYUrlS272mcnu+C0ZwC46y8+V0P5FW65KJYWbjro2BFOOcU3gA0b4I034K23fJTD66/Dgw9WHN+jR0UA9+/vqzJ07qzLgUUkWRkRuDtq3dpHNZxwQsW+tWs9fMu3t9/2URDlDfg2bSqmaOzXz28POECT3ohUpfv/RHM1/I/P1ZA31veXz8e79MHeABTcoT7emsiILoXa+uor736YO9dHQbz7Lrz/vu8vV1j47YnI+/b1yXJy7eScuhSkJr44y0+yfTl0MwB7P90UqJg4PRvlTJdCbbVo4aMcjjyyYl8I8MknOy+389xzUFLixzRqBPvtt3MQawVeEamLrG7h1sTWrTB/vreCKwfx8uUVx7Ru7d0QO3ZNtGyZXN2pohau1MXm4T5pzsoBPrS/63RvvZQvjpkN1MJNocaNPUD79fv2/vXrvRuicghPnlwxVy/4Kr7ly6WX33bsqNawiHybWri1EAJ8+qm3ht95x7c5c3wu33Lt21cE8EEHwSGHeDCnawirhSupVDL4EADWHtiYDrO8n9femJtkSXWmFm5CzKB7d99OPLFi/5dfegjPmVMRwrfc4ssFgY8PPuww3wYM8As32rRJ5mcQkfgpcFOoVSu/Cu4HP6jYt20bfPABzJ4Ns2bBzJk+f0T5B4t99/XwPewwOPxw79JoUO053ETSU/5LxQB0fAnskO8BsHGETxXZ+j1fGKZ03vxkikuQAreeNWpU0a1wzjm+b+NGD+CZMz2En3/elxcCX2n4hz+Eo4/2y5x79kzfbggRqRkFbgJatoTBg32DiqFqr74K06b59uij/lhhYcUcE8cdpy4IyTyheB4ALYqjHX16AlD6w/4ANF7s82SXfPJp7LXFTYGbBsw8WAsL4YwzPIDnz68I38ceg3vv9Ysxhg6FU0+FYcO8C0NEMocCNw2Z+YUX++0HF17oU1e+9ZYH7yOPwDPPeFfFccd5+A4fDk2aJF21SPWULvThPHkLox1dOgPQoN9+ANjn6/y4SivEZAudnskAeXl+Yu3GG2HpUpgxw4O4uBhOP90n6/ntb2HduqQrFZHdUeBmGDMP31tugWXL4O9/9zG+V1/tqx5fdJH3B4tkipIVn1Gy4jPK3v2Isnc/8mvsS0rI79GN/B7daNCiBQ2yZJYpBW4Ga9DAT6Y9+2zFEkUTJ8L++8Ntt3lXhIikDwVuljjgALjvPli4EI46CsaO9XG98+YlXZlIzZSuX0/p+vWUfPKpj1woK4OyMvLa7U1eu71p0KQJDTL0pIUCN8v06AFTp8IDD/ilxocf7gtyikjyFLhZyAx++lOfhL1jRx9KNm1a0lWJ1E7Z5s2Ubd5M6dp1lK5dRwiBEAINmjWjQbNmWH4+lp8ZA64UuFmsWze/mKJnT1+eaMWKpCsSyW0K3CzXoQM88YTP93vuuRVzOIhkqrB1K2Hr1n+2fMtZ48ZY48b+ES9Nr4dX4OaA3r3h2mt9VYvp05OuRiR3KXBzxHnn+fSQEyYkXYlIaoWSEt+ilu8/NcjzLY0ocHNEkyYwahQ89RRU+hQFYQepAAAGdElEQVQmIjFS4OaQoUP9Ip6Z2buwqoifqAgBykp9K+/TTYO+XQVuDhnoK1vz5pvJ1iGSqzJj8JqkRKtWPmrh44+TrkQkRmk0NEct3BxTWOiT3ohI/BS4OaZdO03jKJIUBW6OadMG1q9PugqR3FTtwDWzPDObY2ZTo/v7mNksM1tkZg+bWaP6K1NSpWlT2LIl6SpEclNNWriXAB9Wun8DcGsIoTewHjg7lYVJ/WjSRIErkpRqBa6ZdQV+BNwT3TdgMDAlOmQS8JP6KFBSq2FD2L496SpEclN1W7i3AT8HyqL7ewMbQggl0f3lQJeqvtHMRpvZbDObvSYLF4XLNPn5fvGDiMRvj4FrZicCq0MIxXs6tiohhIkhhKIQQlFBQUFtnkJSqEEDn0BfROJXnQsfjgB+bGYnAE2AlsDtQGszy49auV0BzbaaARS4IsnZYws3hHBVCKFrCKEQGAG8FEIYCUwHTokOGwU8VW9VSsqYpdWFNyI5pS7jcK8ALjWzRXif7r2pKUnqkwJXJDk1mkshhPAy8HL09WLg0NSXJPUpTSfCF8kJutJMRCQmCtwcoxauSHIUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMVHgiojERIErIhITBa6ISEwUuCIiMalW4JpZazObYmYfmdmHZjbQzNqa2YtmtjC6bVPfxYqIZLLqtnBvB54PIewHHAh8CFwJTAsh9AGmRfdFRGQX9hi4ZtYKGATcCxBC2BZC2AAMAyZFh00CflJfRYqIZIPqtHD3AdYAfzKzOWZ2j5k1AzqEED6PjlkJdKjqm81stJnNNrPZa9asSU3VIiIZqDqBmw/0B+4IIRwMbGaH7oMQQgBCVd8cQpgYQigKIRQVFBTUtV4RkYxVncBdDiwPIcyK7k/BA3iVmXUCiG5X10+JIiLZYY+BG0JYCXxqZt+Jdg0BPgCeBkZF+0YBT9VLhSIiWSK/msddDEw2s0bAYuBMPKwfMbOzgU+AU+unRBGR7FCtwA0hvAMUVfHQkNSWIyKSvXSlmYhITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxqVbgmtlYM5tnZu+b2YNm1sTM9jGzWWa2yMweNrNG9V2siEgm22PgmlkXYAxQFEI4AMgDRgA3ALeGEHoD64Gz67NQEZFMV90uhXxgLzPLB5oCnwODgSnR45OAn6S+PBGR7LHHwA0hrABuApbhQfslUAxsCCGURIctB7pU9f1mNtrMZpvZ7DVr1qSmahGRDFSdLoU2wDBgH6Az0Aw4rrovEEKYGEIoCiEUFRQU1LpQEZFMV50uhaOBJSGENSGE7cDjwBFA66iLAaArsKKeahQRyQrVCdxlwAAza2pmBgwBPgCmA6dEx4wCnqqfEkVEskN1+nBn4SfH3gbei75nInAFcKmZLQL2Bu6txzpFRDJe/p4PgRDCr4Ff77B7MXBoyisSEclSutJMRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQJXRCQmClwRkZgocEVEYqLAFRGJiQI3xwwYAJdcknQVIrnJQgjxvZjZGuCT2F5QaqJHCKEg6SJEslmsgSsiksvUpSAiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITBS4IiIxUeCKiMREgSsiEhMFrohITP4PrQ161dhEqnEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2bf8869240>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% EMD\n",
    "\n",
    "G0 = ot.emd(a, b, M)\n",
    "\n",
    "pl.figure(3, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD with Frobenius norm regularization\n",
    "--------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "    0|1.760578e-01|0.000000e+00\n",
      "    1|1.669467e-01|-5.457501e-02\n",
      "    2|1.665639e-01|-2.298130e-03\n",
      "    3|1.664378e-01|-7.572776e-04\n",
      "    4|1.664077e-01|-1.811855e-04\n",
      "    5|1.663912e-01|-9.936787e-05\n",
      "    6|1.663852e-01|-3.555826e-05\n",
      "    7|1.663814e-01|-2.305693e-05\n",
      "    8|1.663785e-01|-1.760450e-05\n",
      "    9|1.663767e-01|-1.078011e-05\n",
      "   10|1.663751e-01|-9.525192e-06\n",
      "   11|1.663737e-01|-8.396466e-06\n",
      "   12|1.663727e-01|-6.086938e-06\n",
      "   13|1.663720e-01|-4.042609e-06\n",
      "   14|1.663713e-01|-4.160914e-06\n",
      "   15|1.663707e-01|-3.823502e-06\n",
      "   16|1.663702e-01|-3.022440e-06\n",
      "   17|1.663697e-01|-3.181249e-06\n",
      "   18|1.663692e-01|-2.698532e-06\n",
      "   19|1.663687e-01|-3.258253e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   20|1.663682e-01|-2.741118e-06\n",
      "   21|1.663678e-01|-2.624135e-06\n",
      "   22|1.663673e-01|-2.645179e-06\n",
      "   23|1.663670e-01|-1.957237e-06\n",
      "   24|1.663666e-01|-2.261541e-06\n",
      "   25|1.663663e-01|-1.851305e-06\n",
      "   26|1.663660e-01|-1.942296e-06\n",
      "   27|1.663657e-01|-2.092896e-06\n",
      "   28|1.663653e-01|-1.924361e-06\n",
      "   29|1.663651e-01|-1.625455e-06\n",
      "   30|1.663648e-01|-1.641123e-06\n",
      "   31|1.663645e-01|-1.566666e-06\n",
      "   32|1.663643e-01|-1.338514e-06\n",
      "   33|1.663641e-01|-1.222711e-06\n",
      "   34|1.663639e-01|-1.221805e-06\n",
      "   35|1.663637e-01|-1.440781e-06\n",
      "   36|1.663634e-01|-1.520091e-06\n",
      "   37|1.663632e-01|-1.288193e-06\n",
      "   38|1.663630e-01|-1.123055e-06\n",
      "   39|1.663628e-01|-1.024487e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   40|1.663627e-01|-1.079606e-06\n",
      "   41|1.663625e-01|-1.172093e-06\n",
      "   42|1.663623e-01|-1.047880e-06\n",
      "   43|1.663621e-01|-1.010577e-06\n",
      "   44|1.663619e-01|-1.064438e-06\n",
      "   45|1.663618e-01|-9.882375e-07\n",
      "   46|1.663616e-01|-8.532647e-07\n",
      "   47|1.663615e-01|-9.930189e-07\n",
      "   48|1.663613e-01|-8.728955e-07\n",
      "   49|1.663612e-01|-9.524214e-07\n",
      "   50|1.663610e-01|-9.088418e-07\n",
      "   51|1.663609e-01|-7.639430e-07\n",
      "   52|1.663608e-01|-6.662611e-07\n",
      "   53|1.663607e-01|-7.133700e-07\n",
      "   54|1.663605e-01|-7.648141e-07\n",
      "   55|1.663604e-01|-6.557516e-07\n",
      "   56|1.663603e-01|-7.304213e-07\n",
      "   57|1.663602e-01|-6.353809e-07\n",
      "   58|1.663601e-01|-7.968279e-07\n",
      "   59|1.663600e-01|-6.367159e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   60|1.663599e-01|-5.610790e-07\n",
      "   61|1.663598e-01|-5.787466e-07\n",
      "   62|1.663596e-01|-6.937777e-07\n",
      "   63|1.663596e-01|-5.599432e-07\n",
      "   64|1.663595e-01|-5.813048e-07\n",
      "   65|1.663594e-01|-5.724600e-07\n",
      "   66|1.663593e-01|-6.081892e-07\n",
      "   67|1.663592e-01|-5.948732e-07\n",
      "   68|1.663591e-01|-4.941833e-07\n",
      "   69|1.663590e-01|-5.213739e-07\n",
      "   70|1.663589e-01|-5.127355e-07\n",
      "   71|1.663588e-01|-4.349251e-07\n",
      "   72|1.663588e-01|-5.007084e-07\n",
      "   73|1.663587e-01|-4.880265e-07\n",
      "   74|1.663586e-01|-4.931950e-07\n",
      "   75|1.663585e-01|-4.981309e-07\n",
      "   76|1.663584e-01|-3.952959e-07\n",
      "   77|1.663584e-01|-4.544857e-07\n",
      "   78|1.663583e-01|-4.237579e-07\n",
      "   79|1.663582e-01|-4.382386e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   80|1.663582e-01|-3.646051e-07\n",
      "   81|1.663581e-01|-4.197994e-07\n",
      "   82|1.663580e-01|-4.072764e-07\n",
      "   83|1.663580e-01|-3.994645e-07\n",
      "   84|1.663579e-01|-4.842721e-07\n",
      "   85|1.663578e-01|-3.276486e-07\n",
      "   86|1.663578e-01|-3.737346e-07\n",
      "   87|1.663577e-01|-4.282043e-07\n",
      "   88|1.663576e-01|-4.020937e-07\n",
      "   89|1.663576e-01|-3.431951e-07\n",
      "   90|1.663575e-01|-3.052335e-07\n",
      "   91|1.663575e-01|-3.500538e-07\n",
      "   92|1.663574e-01|-3.063176e-07\n",
      "   93|1.663573e-01|-3.576367e-07\n",
      "   94|1.663573e-01|-3.224681e-07\n",
      "   95|1.663572e-01|-3.673221e-07\n",
      "   96|1.663572e-01|-3.635561e-07\n",
      "   97|1.663571e-01|-3.527236e-07\n",
      "   98|1.663571e-01|-2.788548e-07\n",
      "   99|1.663570e-01|-2.727141e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  100|1.663570e-01|-3.127278e-07\n",
      "  101|1.663569e-01|-2.637504e-07\n",
      "  102|1.663569e-01|-2.922750e-07\n",
      "  103|1.663568e-01|-3.076454e-07\n",
      "  104|1.663568e-01|-2.911509e-07\n",
      "  105|1.663567e-01|-2.403398e-07\n",
      "  106|1.663567e-01|-2.439790e-07\n",
      "  107|1.663567e-01|-2.634542e-07\n",
      "  108|1.663566e-01|-2.452203e-07\n",
      "  109|1.663566e-01|-2.852991e-07\n",
      "  110|1.663565e-01|-2.165490e-07\n",
      "  111|1.663565e-01|-2.450250e-07\n",
      "  112|1.663564e-01|-2.685294e-07\n",
      "  113|1.663564e-01|-2.821800e-07\n",
      "  114|1.663564e-01|-2.237390e-07\n",
      "  115|1.663563e-01|-1.992842e-07\n",
      "  116|1.663563e-01|-2.166739e-07\n",
      "  117|1.663563e-01|-2.086064e-07\n",
      "  118|1.663562e-01|-2.435945e-07\n",
      "  119|1.663562e-01|-2.292497e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  120|1.663561e-01|-2.366209e-07\n",
      "  121|1.663561e-01|-2.138746e-07\n",
      "  122|1.663561e-01|-2.009637e-07\n",
      "  123|1.663560e-01|-2.386258e-07\n",
      "  124|1.663560e-01|-1.927442e-07\n",
      "  125|1.663560e-01|-2.081681e-07\n",
      "  126|1.663559e-01|-1.759123e-07\n",
      "  127|1.663559e-01|-1.890771e-07\n",
      "  128|1.663559e-01|-1.971315e-07\n",
      "  129|1.663558e-01|-2.101983e-07\n",
      "  130|1.663558e-01|-2.035645e-07\n",
      "  131|1.663558e-01|-1.984492e-07\n",
      "  132|1.663557e-01|-1.849064e-07\n",
      "  133|1.663557e-01|-1.795703e-07\n",
      "  134|1.663557e-01|-1.624087e-07\n",
      "  135|1.663557e-01|-1.689557e-07\n",
      "  136|1.663556e-01|-1.644308e-07\n",
      "  137|1.663556e-01|-1.618007e-07\n",
      "  138|1.663556e-01|-1.483013e-07\n",
      "  139|1.663555e-01|-1.708771e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  140|1.663555e-01|-2.013847e-07\n",
      "  141|1.663555e-01|-1.721217e-07\n",
      "  142|1.663554e-01|-2.027911e-07\n",
      "  143|1.663554e-01|-1.764565e-07\n",
      "  144|1.663554e-01|-1.677151e-07\n",
      "  145|1.663554e-01|-1.351982e-07\n",
      "  146|1.663553e-01|-1.423360e-07\n",
      "  147|1.663553e-01|-1.541112e-07\n",
      "  148|1.663553e-01|-1.491601e-07\n",
      "  149|1.663553e-01|-1.466407e-07\n",
      "  150|1.663552e-01|-1.801524e-07\n",
      "  151|1.663552e-01|-1.714107e-07\n",
      "  152|1.663552e-01|-1.491257e-07\n",
      "  153|1.663552e-01|-1.513799e-07\n",
      "  154|1.663551e-01|-1.354539e-07\n",
      "  155|1.663551e-01|-1.233818e-07\n",
      "  156|1.663551e-01|-1.576219e-07\n",
      "  157|1.663551e-01|-1.452791e-07\n",
      "  158|1.663550e-01|-1.262867e-07\n",
      "  159|1.663550e-01|-1.316379e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  160|1.663550e-01|-1.295447e-07\n",
      "  161|1.663550e-01|-1.283286e-07\n",
      "  162|1.663550e-01|-1.569222e-07\n",
      "  163|1.663549e-01|-1.172942e-07\n",
      "  164|1.663549e-01|-1.399809e-07\n",
      "  165|1.663549e-01|-1.229432e-07\n",
      "  166|1.663549e-01|-1.326191e-07\n",
      "  167|1.663548e-01|-1.209694e-07\n",
      "  168|1.663548e-01|-1.372136e-07\n",
      "  169|1.663548e-01|-1.338395e-07\n",
      "  170|1.663548e-01|-1.416497e-07\n",
      "  171|1.663548e-01|-1.298576e-07\n",
      "  172|1.663547e-01|-1.190590e-07\n",
      "  173|1.663547e-01|-1.167083e-07\n",
      "  174|1.663547e-01|-1.069425e-07\n",
      "  175|1.663547e-01|-1.217780e-07\n",
      "  176|1.663547e-01|-1.140754e-07\n",
      "  177|1.663546e-01|-1.160707e-07\n",
      "  178|1.663546e-01|-1.101798e-07\n",
      "  179|1.663546e-01|-1.114904e-07\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  180|1.663546e-01|-1.064022e-07\n",
      "  181|1.663546e-01|-9.258231e-08\n",
      "  182|1.663546e-01|-1.213120e-07\n",
      "  183|1.663545e-01|-1.164296e-07\n",
      "  184|1.663545e-01|-1.188762e-07\n",
      "  185|1.663545e-01|-9.394153e-08\n",
      "  186|1.663545e-01|-1.028656e-07\n",
      "  187|1.663545e-01|-1.115348e-07\n",
      "  188|1.663544e-01|-9.768310e-08\n",
      "  189|1.663544e-01|-1.021806e-07\n",
      "  190|1.663544e-01|-1.086303e-07\n",
      "  191|1.663544e-01|-9.879008e-08\n",
      "  192|1.663544e-01|-1.050210e-07\n",
      "  193|1.663544e-01|-1.002463e-07\n",
      "  194|1.663543e-01|-1.062747e-07\n",
      "  195|1.663543e-01|-9.348538e-08\n",
      "  196|1.663543e-01|-7.992512e-08\n",
      "  197|1.663543e-01|-9.558020e-08\n",
      "  198|1.663543e-01|-9.993772e-08\n",
      "  199|1.663543e-01|-8.588499e-08\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  200|1.663543e-01|-8.737134e-08\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEYCAYAAAAeWvJ8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XecVdW5//HPMzP0QelIH7qCxsKI2L3BguUGNbFfb6zEgiYxlhtL4r35aWI0iUZjQaPivZrYY4mJLZYogoIaEJSqNKkC0oRpz++PtY8zjAxTz1lnZr7v12u/Tttn72cOzPnOXnvttczdERERiSUndgEiItK8KYhERCQqBZGIiESlIBIRkagURCIiEpWCSEREolIQiUidmdnBZjY7C+r4zMwOj12H1I2CSCQiMzvLzGaY2WYzW25md5lZh+S1u81sY7IUmVlxhcd/y0BtbmaDdrSOu//T3YfWcfufmdlXFX6mjWbWs27VSmOmIBKJxMx+AtwEXAHsDIwC+gEvm1lLd7/A3fPdPR+4EXg09djdj45XeWBmeQ2wmX+v8DPlu/vnadpPjZlZbib3JwoikSjMbCfgv4FL3P3v7l7s7p8BJwMFwH/UYZuHmdkSM7vSzFaa2TIzO97MjjGzOWa2xsyurrD+SDN7x8zWJeveYWYtk9feTFb7V3KkckqF7V9lZsuBB1LPJe8ZmOxjn+RxTzNbZWaH1fLnKEiOxs41s0XAP5Lnv2NmM5N6Xzez3Sq9dV8zm2Vma83sATNrXcP9PZgcib5gZpuAfzOzVmZ2i5ktMrMVydFpmwrvuTL5zD43s/NqcvQoVVMQicRxANAaeKrik+6+EXgBOKKO290l2W4v4GfAvYRQGwEcDFxnZv2TdUuBHwNdgP2B0cBFSR2HJOvsmRypPFph+50IR27jKtU+H7gK+D8zaws8AEx099fr+LMcCuwGHGVmQ4A/AT8CuhI+o+dSwZk4AzgKGAgMAa6txb5OB24A2gNvAb9KtrEXMIjyzxMzGwNcBhyevHZYnX46+ZqCSCSOLsBqdy/ZzmvLktfrohi4wd2LgT8n27nN3Te4+0xgFrAngLtPc/fJ7l6SHI3dQ/jy35Ey4OfuvtXdv6r8orvfC8wDpgA9gGuq2d5fkiOcdWb2l0qvXe/um5L9nAL81d1fTn62W4A2hEBPucPdF7v7GkKonFbNvit6xt3fdvcyYCshZH/s7mvcfQOhafTUZN2TgQfcfaa7bwaur8V+ZDsy2vYqIl9bDXQxs7zthFGP5PW6+MLdS5P7qaBYUeH1r4B8gOQo47dAIdCW8H0wrZrtr3L3LdWscy/wLDDO3bdWs+7x7v5KFa8trnC/J7Aw9cDdy8xsMeFIZXvrL0zeU1MV39uV8HlMM7PUcwakzh31BKZW8V6pAx0RicTxDuEv7xMrPmlm+cDRwKsZqOEu4BNgsLvvBFxN+MLdkR0O15/UfyvwR+B6M+tUj/oq7utzQnNgaj8G9AGWVlinT4X7fZP31GVfqwmBPdzdOyTLzkmnEQhHrL2r2K/UgYJIJAJ3/5LQWeF2MxtjZi3MrAB4DFgC/G8GymgPrAc2mtmuwIWVXl8BDKjlNm8Dprr7ecBfgbvrXWXwGHCsmY02sxbATwhBPqnCOhebWe8k/K4BHt3OdqqVNM/dC/zOzLoBmFkvMzuqQi1nm9luybmw6+r2I0mKgkgkEnf/NeEo5BZCIEwhNPOMrkGTVkO4nHCSfgPhi7fyF/f1wMTk/M3J1W3MzMYCYygPtMuAfczsjPoW6u6zCZ0ubiccsfw7oet3UYXVHgFeAhYA84H/l9TVN+n517cWu7yKcK5rspmtB14Bhia1/A34PfBaap3kPZn4N2uSTBPjiYjUXdKN/COgVRWdT6QaOiISEaklMzshudaoI+Gi5OcUQnWnIBIRqb0fACsJTYClfPP8mtSCmuZERCQqHRGJiEhUuqBVMqJLly5eUFAQuwwRyaBp06atdveu1a2nIJKMKCgoYOrUqdWvKCJNhpktrH4tNc2JSDZyh8WLYeZMKFFntKZOQSQi2cMdHnkECgqgb1/YfXfo0AGuvBI2bYpdnaSJgkhEskNJCZx+OpxxBnTvDnfcAQ89BMcfDzffDHvtBYsWxa5S0kBBJCLxlZaGAPrzn+GGG+Cdd+Dii+HMM+H//g9efx1WrYLDDoMlS2JXKw1MQSQi8f3yl/DYY/DrX8PVV0Nupdm6Dz0UXnoJVq+GU06B4uI4dUpaKIhEJK633oKf/zw0y11+edXrjRwJEybApElhfWkyFEQiEs/WrXDOOaFzwl13gVUzHdKpp4b1b7oJPvwwIyVK+imIRCSe3/0O5s6FP/wBdtqpZu+55Rbo3BnGjw+97KTRUxCJSByffw6/+EXoFTdmTM3f17Ej/OpX8Pbboau3NHoKIhGJ48YboagoHOHU1llnhe7cP/uZOi40AQoiEcm8RYvg3nvD+Z6BA2v//pyccDS1YAFMnNjw9UlGKYhEJPNuvDHcXntt3bdx7LGw334hkHRU1KgpiEQks5YvhwcfDM1rffrUfTtmcN114ejq0UcbqjqJQEEkIpl1++3h3NCOrhmqqaOPhuHDw4Ww6kHXaCmIRCRzNm6EO++EE0+EwYPrv72cHLjiCpgxA158sf7bkygURCKSORMnwrp1DXM0lHLaadCzJ9x6a8NtUzJKQSQimVFWFkbU3ndfGDWq4bbbsiVccEE4Ipo9u+G2KxmjIBKRzHjlFfjkE7j00obf9rhxIZDuuKPhty1ppyASkcy44w7o1g1OOqnht929exiH7sEHYcOGht++pJWCSETSb+FCeP55OP98aNUqPfu46KLQGeLhh9OzfUkbBZGIpN+994brfsaNS98+Ro6EvfcOo3irK3ejoiASkfQqKoL77oNjjoG+fdO3H7PQaWH69DDDqzQaCiIRSa9nn4UVK+DCC9O/r9NPh/bt4e67078vaTAKIhFJr3vuCUdCRx2V/n3l58N//Ac8/jisXZv+/UmDUBCJSPrMnx+6bZ93HuTmZmaf48bBli3wv/+bmf1JvSmIRCR97rsvBNA552Run3vtFTou3HOPOi00EgoiEUmPoiK4/3447jjo1Suz+x43DmbNgkmTMrtfqRMFkYikxzPPwMqV8IMfZH7fp54aOi3cc0/m9y21piASkfSYMCF0UjjyyMzvu1270GnhscdgzZrM719qRUEkIg1v3rzQSeH88zPXSaGyH/wAtm6Fhx6Ks3+pMQWRiDS8CRMy30mhsj33DFOJq9NC1lMQiUjD2rIldFIYOzbMExTTBReEEb/feCNuHbJDCiIRaVhPPglffJGZkRSqc8op0KGDRlrIcgoiEWlYd90FgwbBt78duxJo0wbOOgueegqWL49djVRBQSQiDeeDD+Dtt8OUDDlZ8vVy4YVQXBzOW0lWypL/KSLSJNxxB7RtC2efHbuSckOGhHHu7r47BJJkHQWRiDSML76ARx6BM88M52WyySWXwLJloYlOso6CSEQaxoQJocfc+PGxK/mmo4+GgQPhtttiVyLboSASkforKoLbbw+jKOy+e+xqviknB370ozBhnibNyzoKIhGpvz//OTR9XXZZ7EqqdtZZocnwt7+NXYlUoiASkfpxh1tugWHD4owrV1P5+WHYn6eeCkMQSdZQEIlI/Tz/PMyYAVddBWaxq9mxH/4QWrSAX/86diVSgYJIROrOHW64AQoK4LTTYldTvR49wvh3Dz4IS5fGrkYSCiIRqbtXX4UpU+DKK8ORRmNwxRVQVqajoiyiIBKRunGHa6+F3r2z6wLW6vTvHzou3H03LFoUuxpBQSQidfXcc+Fo6Gc/g9atY1dTO9ddF4L0F7+IXYmgIBKRuigpgWuuCReJnnVW7Gpqr1+/0IPugQfg449jV9PsKYhEpPbuuw8++ghuuqnxnBuq7LrrwpTil18eu5JmT0EkIrWzbl34Ej/0UDjxxNjV1F23bqFZ8YUX4G9/i11Ns6YgEpHaueoqWLMGbr01+68bqs4ll4TRucePh82bY1fTbCmIRKTm3nwzDG764x/DXnvFrqb+WraEe+6BBQvgv/87djXNloJIRGpmwwY499xw8WpT+tI+7DA47zz4zW80IGokCiIRqZnx48ORw8SJ4SR/U3LLLdCnTxgdYt262NU0OwoiEanefffBQw+FTgqHHBK7moa3887wpz/BkiWhO3pZWeyKmhUFkYjs2JtvwkUXhZG1r702djXpM2pUaJ575pkQuJIxebELEJEs9uGHMHYsDBgAjz4KeU38K+PSS8P1UTfeGAZIzcbZZpugJv6/SkTq7MMP4YgjoH17+Pvfw6RyTZ0Z3HknrFoVunbn5sKFF8auqslT05yIfNNLL4VzQa1bhxG2CwpiV5Q5LVqEo7/jjgtNktdeq3NGaaYgEpFyJSVw/fUwZkwYpXryZBg8OHZVmdeqFTz9dOjWfcMNcMwxsGJF7KqaLAWRiASTJkFhYbhG6D//E95+G3r1il1VPHl54eLdu++GN96A3XYLj0tKYlfW5CiIRJqzsjJ45ZXwF/+BB8Lq1fDEE2EG0/z82NXFZxZG6X7/ffjWt8L9YcNCd/ZNm2JX12QoiESamy1b4LXXwphxgweHDgnvvQe//CV88gl897uxK8w+u+0WPrOnnw4X855/PvTsGY4cH388dG6QOjN3j12DNAOFhYU+derU2GU0faWl4S/19evDwKQrV8KyZbBwIcydG7omz5gBxcWh6Wn0aDjjDDjppMY3uV0s7vDWW3D//eGao7Vrw/MDBsCee8Kuu4bOHb16Qffu0LlzuGC2ffvGO2VGHZnZNHcvrHY9BZFkQpMJohkz4PTT07f9ir+P7uWP3UMzWuq2rCyETnFxWIqKwpFOUVHV2+7dO/xlP2JEaIY75BDYaaf0/SzNQUkJTJ0aziG99x7MnAnz5lV9Hik3F9q0CZ0hWrQIfwzk5UFODuyxB/zlL5mtP81qGkS6jkikNlq3DtMGpFPFqRXMyh/n5ITFLNymvsTy8sIXW+vW4UsuPz8ETMeOYc6dXXYJ46i1aZPeupujvLwwIsOoUeXPlZbC55+HZcWKcGT65Zdh0NjNm8MfDFu3hrAqLg63ZWXNq4t8JQoikdoYPBiefDJ2FZLNcnND8PfpE7uSRkOdFUREJCqdI5KMMLNVwMI6vLULsLqBy2kI2VhXNtYEqqu2srGuutbUz927VreSgkiymplNrcnJzkzLxrqysSZQXbWVjXWluyY1zYmISFQKIhERiUpBJNluQuwCqpCNdWVjTaC6aisb60prTTpHJCIiUemISEREolIQiYhIVAoiyTpmdrOZfWJm083saTPrkDxfYGZfmdmHyXJ3hNrGmNlsM5tnZv+V6f1XqKOPmb1mZrPMbKaZ/TB5/nozW1rhMzomQm2fmdmMZP9Tk+c6mdnLZjY3ue2YwXqGVvg8PjSz9Wb2oxiflZndb2YrzeyjCs9t97Ox4PfJ/7XpZrZPhuvK2O+hzhFJ1jGzI4F/uHuJmd0E4O5XmVkB8Ly77x6prlxgDnAEsAR4DzjN3WdFqKUH0MPd3zez9sA04HjgZGCju9+S6Zoq1PYZUOjuqys892tgjbv/Kgnwju5+VYTacoGlwH7A2WT4szKzQ4CNwEOp/8dVfTZJMF4CHJPUe5u775fBujL2e6gjIsk67v6Su6eGL54M9I5ZTwUjgXnuvsDdi4A/A2NjFOLuy9z9/eT+BuBjIJunUx0LTEzuTySEZgyjgfnuXpdRPurN3d8E1lR6uqrPZiwhGNzdJwMdkj9AMlJXJn8P6xVE2dJMIU3aOcDfKjzub2YfmNkbZnZwhmvpBSyu8HgJWfDln/yFujcwJXlqfNKccn8mm8AqcOAlM5tmZuOS57q7+7Lk/nKge4S6AE4F/lThcezPCqr+bLLp/1tafw/rHETJIe4fgKOBYcBpZjasvgVJ82Bmr5jZR9tZxlZY5xqgBHg4eWoZ0Nfd9wYuAx4xs2Y9oY6Z5QNPAj9y9/XAXcBAYC/C5/WbCGUd5O77EL4bLk6afb7m4XxAxs8JmFlL4DvA48lT2fBZbSPWZ7Mjmfg9rM80EF83UyTFppopMt5eLo2Pux++o9fN7CzgOGB08suJu28Ftib3p5nZfGAIkKkZ95YCFcf27508F4WZtSCE0MPu/hSAu6+o8Pq9wPOZrsvdlya3K83sacJ3xQoz6+Huy5LmpZWZrosQjO+nPqNs+KwSVX020f+/Zer3sM6dFczse8AYdz8veXwmsJ+7j6/qPV26dPGCZjz5U1P32Wdh7q899vjma9OmTVtdk1F4ITT5Ar8FDnX3VRWe70o4qVtqZgOAfwJ7uHvlNve0MLM8YM7h9r3+mdifNE+v+BM3u/uVZnYsMJ7yzgq/d/eR6dpv5U4Imfw9TPvEeEkb8TiAvn370iSmi5btOuUUmD49zJxcmZnV5uTwHUAr4GULs5NOdvcLgEOA/zGzYqAMuCBTIQSQ9B4aD/w1U/uUZulXye0LhBCaB2wm9PJLCzP7E3AY0MXMlgA/B35Khn4P6xNENTpsdPcJJOMUFRYWZlXbpzSsjRuhbdv6b8fdB1Xx/JOEpqho3P2FI3JOilnCN6WmEtelGE1C6ks9aQq7OEP7PG07T/+xinUb/PewPr3m3gMGm1n/5CTgqcCzDVOWNEZr1kDHWP2ORKTRqnMQJf3LxwMvEq5heMzdZzZUYdL4LFkCvbPlip/mxD0sZmHJyS2/L9II1Osckbu/QGjHlGZuzZoQREOHxq6kGUs1zXlpCCPA8nLx0tLwfFlppMJEdkwjK0iDeO21cLv//nHrEJHGJ+295qR5uO8+6NkTDjoodiUCfH3042WlWIuWAFi7tviWreH54qJopYlUpiCSenvmGfj73+HGGyFP/6OyTip0vLiInKRbY06nDvjmrwAo27AhWm0ioKY5qafJk+Hss2HPPeEnP4ldjYg0RgoiqRN3eOQRGD0aOnWCp5+Gli1jVyXVKdu8mbLNmyldsRLLy8Py8sgdMpDcrl3J7VqjgS9EGpwaUqTWpk+HSy+FN96AkSND09wuu8SuSmqrdO3acGftWvL6JP3uR+5B7pqNAJQtXKpzSZIROiKSGikqgscfhyOOCM1wM2bA3XfDpEkKIRGpHx0RSZW2boV//hNeeAEefhhWroS+feF//gcuugg6d45doTSUksVLALCVq2BQAQDFh+xB3qZiAHLnL6N01aqq3i5SLwoi+Zo7zJsHL74YesG99hps3hzO/RxzDIwbB0ceCbm5sSuVdPGtWymdORuA1l90Z+uuYR62zQcPwMrCoOP5s9dS+vHcaDVK06MgasZWrID33gvLu++G2y++CK8NHgznngtjxsChh0K7dnFrFZGmS0HUDBQVhSOdWbPg44/hX/8KobNoUXg9JweGD4exY0Png8MPh4ED49Ys8ZUsX0Hu8jB3XPvhQ1m7ZxjRduWBXeCALgDs/GkRLSd/DIQeeSJ1oSBqQr76CubMKQ+cWbPCMnculJSUrzdgQBiK59JLQ/DsvTfk58erW7Jf6czZ7DQrDKJa8m/7sHqPVgAs368V7LcXAO0+d7q8EyYXLZ0zP06h0igpiBoR99CctmBBWObP3/b+smXl6+bmhqOaYcPghBNgt93C/aFD1cwmItlFQZRF3EPPtEWLtl0+/bQ8dCq2fpiFaRcGDICjjw63gweHwBk8GFq1ivezSBOUjO6d949p9J4TOjF8cWgf1u4WjpRWF5ax5lvhotjWq7rT/b0wrl3eq9MiFCuNiYIogzZvhsWLvxk0qWXx4tBluqK2bUPADBhQfu5mwIBw268ftG4d52eR5q1kSZiMeeeHl7LzqG8BsPTQfDb1D23AZf028tnuLcLKx40CoPOHRpfnQ4+80i8yNsO7NAIKogayfn2Yj6fysnRp+f01lX73zMKI1X37wogRoQmtb99tl44dNb+ZiDRtCqJquMPatdsPmYrL9gYw7tYtNJ0VFITpEXr12jZkevbU+GzSBEyeDkCvybD5xP0A+PyQfNoVfAnAkH6LAVg9PJ85R/UFoNWMXen7QhhiqOxfH2e6YskyzTqIyspg1aqqwyV1NPPVV9u+LycnDGvTp084H3PkkSFwevcOYdO7dwgZnaOR5qbtU1MAGPpmZ5adHKbr/ddB4RdhRL9FHLlLCJ0vd2vD24cOAGDZJ6Po90Jo0mvx0tRMlyxZoEkHUVFROO+ycGE4B7Nw4bb3t3dOJi+vPEz22Qe+853ykEktu+yieXdERBpKo/46TfUymzOnfEmFzcKFsHz51x19vtajR2gW22efcE6mT5+wpEKmW7dwxCMidVe6+gu63TkJgO7vDAfgwxOGsuqAcMHasbt8xE8HvQBA8cA8nhhZCMCk00fQ46/hayn/8SmZLlsiaRRBtGFDuChzzhyYPXvb4Fm/vny9Fi1CT7K+fcPQNKn7qds+fdRcJpJp/sFMAPp9AEVHhcC5fezhHLz3JwD8oPvr3Nj7eQCKe8GEPcN888+fOpx2z+0EQMcH38l02ZJBWRdEqYE33347LJMmhdEBUsxCqAwZAmeeGW5TS79+GpBTRKSxMa/cdlV5BbM+wENAd8CBCe5+m5l1Ah4FCoDPgJPdfe2OtlVYWOhTp37zZGRZGUycGCZYmzQpdCAA6NAhDEWz//5hLLQhQ8L1M23a1PbHlNjMbJq7F8auoyEckXPSjn9ppEY2nBquL1r5nS2cMfw9AC7p9C7tc0JX0mIvZcK6YQA8/GkhOU+HeUc63d+0j45eLnu82V2wUZMg6gH0cPf3zaw9MA04HjgLWOPuvzKz/wI6uvtVO9rW9oJo3jw477ww2+eAAXDwwXDggXDAAWFYGp2vaRoURLIja7+/PwBbT1jH9weFc0OXdVrw9etbvZjfr90VgIlzRtH22aTJbuLkb54IbuSaYxBV2zTn7suAZcn9DWb2MdALGAsclqw2EXgd2GEQVbZpExQWwpdfwu9/D+PH6+JNEZHmptojom1WNisA3gR2Bxa5e4fkeQPWph5Xes84YBxA3759RyxcuPDr19zhiivgt78NF33+9KfhiGjoUAVSU6MjIqmpL84PR0fFx67j9IGhBeXijjNoY6HJbm3ZV9y+ZiQAD8/cl+5/CeNc5T82OUK1Da85HhHVOIjMLB94A7jB3Z8ys3UVg8fM1rp7xx1to6pzRG+9FZrnZodhqOjUKTTNpZro9t1X54UaOwWR1MXas0Iobfr39YwdMAOACzpPojT5F1hcms99Kw4B4J8f7ErBM2UAtHyx8V4Y2xyDqEa95sysBfAk8LC7P5U8vcLMerj7suQ80sq6FnHQQWH+nNmzQ2eFVI+550OPTnJzw/mjij3kUkuvXjp6EhFpzGrSWcEI54DWuPuPKjx/M/BFhc4Kndz9yh1tq6ojoqqsXg3vvANTppRfQzR37rZD7rRtG6Y8SAVT//7bXjuk0amzg46IpL42nhx62S07toiDhswD4MQu79PaigGYtaUXz34eRgJf/k5P+j+9DoCyD2dtZ2vZqzkeEdUkiA4C/gnMAMqSp68GpgCPAX2BhYTu2zsc2722QbQ9ZWVhDLiKF7Wmlk8/hdLSbdfv1i2EUuWLW1PPaXTrzFAQSUMqGT0CgM+ObUGv4WE68/26fkbHvDBh1/QNvZi2MAywmv9WW3o8Gi6eLVu/ES8uilBxzTXHIKpJr7m3gKo+mNENW071cnLKh+UZXWnvxcUhpCqPKbdwIcyYEZr6tmzZ9j3t2n1zLLnKS+fOCisRkXTJupEV6qNFi9D7rqBg+6+7h4tlUwGVuk2Nsv3qq/D55+Goq6JWraoPK41RJ5I5qVlfB70KOd8K1xe9cNz+FO+xCYA+XdcyvNcyABYc2YnZg4cA0OuNMtr85d0IFcuONKkgqo5ZCIxu3cL1S9tTUgIrVnxzUrvUMmlSuC0u3vZ9qVG7UyN3b2/p0UOjdouIVKavxUoqBkpVyspCR4qq5jH64AN47rntz2OUGv27qkXnrERqp2x6OP/Tezrkdu4EwOrjhvLpiNC0YZ2KaNVnIwCLj2pL/pADAOj5ejKb5bszMlyxVFarC1rrqyE6KzQWqZldKx5VLV4clkWLypeiSudN27X7Zjj17x+6rw8cCF27Ns6gUmcFiaX48BGs2DcMu7+1k1PWJgRUi3WhLb3TLKfTm2EW2ZIlS+MUWYE6K0iDMQsX5nbqBHvssf11UjPEVgymissHH4T5lirKzw+hlFoGDiy/7ddPU4+LSOOjIIooJwe6dw/Lvvtuf52vvgodKubPhwULym/nzoUXX9y2+S8nJ4TRsGFh2W238tuddsrMzySSbVq8Mo3er4T7uUMGsnr/bgBs6hkOPNYXGF916QdA5492Ie8f06LU2ZwpiLJcmzaw665hqcw9zEKbCqj588P1VLNmwSuvbDsNeq9e5QE1fHgIvuHDQ09DkeaidM58Os6ZD0Dndu0AKBq1K18WhKaEL/u3pMXp4cLZjh+uoXTWnDiFNjMKokbMLHR+6NEjjMtXUWlpuMB31qxtl3vvhc3hmj9atw5Tpu+7L4wcGW4HDWqc56BEpPFSZ4VmpqwsHEG9915Y3n0X3n+/vImvY8dwofCYMWHZUe/B2lBnBWkMcocNYeOQMJaz5xit1oTrNFp+vITSFXUeTrNW1FlBmrycnHDUM2gQnHZaeK6kBGbODME0aRK89BI88UR4bY89QiB997vhqElHS9KUlc6aQ5tkaLrc7t0o7b8LAFt370OrLmFygbI5n2b9MEGNjcYCEPLyYM89w1Qc998fupjPmAE33xwu/r31Vhg1CvbaC+68M0xkKCLSUBRE8g1msPvucPnlodPD6tVwzz0hsC6+OJyTuvpq2LgxdqUi6VO6YiVMng6Tp9PizRlQVAxFxfiIXcnr1ZO8Xj1jl9hkKIikWjvtBOPGwbRpMHUqnHAC/PKXoSdfqglPpCnz4iJK5y6gdO4CmDwdLynBS0rIHTKQ3A47k9th59glNmoKIqmVESPg4YfDxIXdusFJJ8ENN4Su5CIidaHOClInBxwQJiw85xy49toQRNdeG7sqkcz4ugfdipXkJleL53bvhn+5HoCyyvPNyA6LxHASAAAGWklEQVQpiKTOWrSAiRPDNUvXXx9611U1qrlIU1W6PoQP69eTk0wJndO+PZ5cE+ElJbFKazTUNCf1kpMTetJ17QrXXBO7GhFpjBREUm8dOsAFF8DLL4fBWkWaq7ItW8KyYcPXz1mLlqErqi7Cq5KCSBrEd78bzhO9/nrsSkSyQ6pnnRcXgeWEJSc3dllZSUEkDWLYsDCX0jQNXCwitaTOCtIgcnLCBH4LF8auRCQLlZWW30810emah6/V+IjIzHLN7AMzez553N/MppjZPDN71Mw0JVsz161bmOhPRHbAXSFUSW2a5n4IfFzh8U3A79x9ELAWOLchC5PGZ+edIdWTVUSkpmoURGbWGzgWuC95bMC3gdQALxOB49NRoDQebdpsO2OsiEhN1PSI6FbgSqAsedwZWOfuqSu1lgDbnbnGzMaZ2VQzm7pK7TZNWsuWUFwcuwoRaWyqDSIzOw5Y6e516g/l7hPcvdDdC7t27VqXTUgjkZsbRlkQEamNmvSaOxD4jpkdA7QGdgJuAzqYWV5yVNQbWJq+MqUxyMlREIlI7VV7ROTuP3X33u5eAJwK/MPdzwBeA76XrPZ94Jm0VSmNgi4cF5G6qM8FrVcBl5nZPMI5oz82TEnSmKlXqojUVq0uaHX314HXk/sLgJENX5I0VmYKIhGpPQ3xIw1GTXMiUhcKIhERiUpBJCIiUSmIREQkKgWRiIhEpSASEZGoFEQiIhKVgkhERKJSEImISFQKIhERiUpBJCIiUSmIREQkKgWRiIhEpSASEZGoFEQiIhKVgkhERKJSEImISFQKIhERiUpBJCIiUSmIREQkKgWRiIhEVaMgMrMOZvaEmX1iZh+b2f5m1snMXjazucltx3QXKyIiTU9Nj4huA/7u7rsCewIfA/8FvOrug4FXk8ciIiK1Um0QmdnOwCHAHwHcvcjd1wFjgYnJahOB49NVpIiINF01OSLqD6wCHjCzD8zsPjNrB3R392XJOsuB7tt7s5mNM7OpZjZ11apVDVO1iIg0GTUJojxgH+Aud98b2ESlZjh3d8C392Z3n+Duhe5e2LVr1/rWKyIiTUxNgmgJsMTdpySPnyAE0woz6wGQ3K5MT4kiItKUVRtE7r4cWGxmQ5OnRgOzgGeB7yfPfR94Ji0ViohIk5ZXw/UuAR42s5bAAuBsQog9ZmbnAguBk9NTooiINGU1CiJ3/xAo3M5Loxu2HBERaW40soKIiESlIBIRkagURCIiEpWCSEREolIQiYhIVAoiERGJSkEkIiJRKYhERCQqBZGIiESlIBIRkagURCIiEpWCSEREolIQiYhIVAoiERGJSkEkIiJRKYhERCQqBZGIiESlIBIRkagURCIiEpWCSEREolIQiYhIVDUKIjP7sZnNNLOPzOxPZtbazPqb2RQzm2dmj5pZy3QXKyIiTU+1QWRmvYBLgUJ33x3IBU4FbgJ+5+6DgLXAueksVEREmqaaNs3lAW3MLA9oCywDvg08kbw+ETi+4csTEZGmrtogcvelwC3AIkIAfQlMA9a5e0my2hKg1/beb2bjzGyqmU1dtWpVw1QtIiJNRk2a5joCY4H+QE+gHTCmpjtw9wnuXujuhV27dq1zoSIi0jTVpGnucOBTd1/l7sXAU8CBQIekqQ6gN7A0TTWKiEgTVpMgWgSMMrO2ZmbAaGAW8BrwvWSd7wPPpKdEERFpympyjmgKoVPC+8CM5D0TgKuAy8xsHtAZ+GMa6xQRkSYqr/pVwN1/Dvy80tMLgJENXpGIiDQrGllBRESiUhCJiEhUCiIREYlKQSQiIlEpiEREJCoFkYiIRKUgEhGRqBREIiISlYJIRESiUhCJiEhUCiIREYlKQSQiIlEpiEREJCoFkYiIRKUgEhGRqBREIiISlYJIRESiUhCJiEhUCiIREYlKQSQiIlEpiEREJCoFkYiIRKUgEhGRqBRE0mD69YO99opdhYg0NubumduZ2SpgYcZ2KNmkn7t3jV2EiGSfjAaRiIhIZWqaExGRqBREIiISlYJIRESiUhCJiEhUCiIREYlKQSQiIlEpiEREJCoFkYiIRKUgEhGRqBREIiISlYJIRESiUhCJiEhUCiIREYlKQSQiIlEpiEREJCoFkYiIRKUgEhGRqBREIiISlYJIRESiUhCJiEhUCiIREYlKQSQiIlH9f33V1Cl94bk5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2c0386eb70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Example with Frobenius norm regularization\n",
    "\n",
    "\n",
    "def f(G):\n",
    "    return 0.5 * np.sum(G**2)\n",
    "\n",
    "\n",
    "def df(G):\n",
    "    return G\n",
    "\n",
    "\n",
    "reg = 1e-1\n",
    "\n",
    "Gl2 = ot.optim.cg(a, b, M, reg, f, df, verbose=True)\n",
    "\n",
    "pl.figure(3)\n",
    "ot.plot.plot1D_mat(a, b, Gl2, 'OT matrix Frob. reg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD with entropic regularization\n",
    "--------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "    0|1.692289e-01|0.000000e+00\n",
      "    1|1.617643e-01|-4.614437e-02\n",
      "    2|1.612639e-01|-3.102965e-03\n",
      "    3|1.611291e-01|-8.371098e-04\n",
      "    4|1.610468e-01|-5.110558e-04\n",
      "    5|1.610198e-01|-1.672927e-04\n",
      "    6|1.610130e-01|-4.232417e-05\n",
      "    7|1.610090e-01|-2.513455e-05\n",
      "    8|1.610002e-01|-5.443507e-05\n",
      "    9|1.609996e-01|-3.657071e-06\n",
      "   10|1.609948e-01|-2.998735e-05\n",
      "   11|1.609695e-01|-1.569217e-04\n",
      "   12|1.609533e-01|-1.010779e-04\n",
      "   13|1.609520e-01|-8.043897e-06\n",
      "   14|1.609465e-01|-3.415246e-05\n",
      "   15|1.609386e-01|-4.898605e-05\n",
      "   16|1.609324e-01|-3.837052e-05\n",
      "   17|1.609298e-01|-1.617826e-05\n",
      "   18|1.609184e-01|-7.080015e-05\n",
      "   19|1.609083e-01|-6.273206e-05\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   20|1.608988e-01|-5.940805e-05\n",
      "   21|1.608853e-01|-8.380030e-05\n",
      "   22|1.608844e-01|-5.185045e-06\n",
      "   23|1.608824e-01|-1.279113e-05\n",
      "   24|1.608819e-01|-3.156821e-06\n",
      "   25|1.608783e-01|-2.205746e-05\n",
      "   26|1.608764e-01|-1.189894e-05\n",
      "   27|1.608755e-01|-5.474607e-06\n",
      "   28|1.608737e-01|-1.144227e-05\n",
      "   29|1.608676e-01|-3.775335e-05\n",
      "   30|1.608638e-01|-2.348020e-05\n",
      "   31|1.608627e-01|-6.863136e-06\n",
      "   32|1.608529e-01|-6.110230e-05\n",
      "   33|1.608487e-01|-2.641106e-05\n",
      "   34|1.608409e-01|-4.823638e-05\n",
      "   35|1.608373e-01|-2.256641e-05\n",
      "   36|1.608338e-01|-2.132444e-05\n",
      "   37|1.608310e-01|-1.786649e-05\n",
      "   38|1.608260e-01|-3.103848e-05\n",
      "   39|1.608206e-01|-3.321265e-05\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   40|1.608201e-01|-3.054747e-06\n",
      "   41|1.608195e-01|-4.198335e-06\n",
      "   42|1.608193e-01|-8.458736e-07\n",
      "   43|1.608159e-01|-2.153759e-05\n",
      "   44|1.608115e-01|-2.738314e-05\n",
      "   45|1.608108e-01|-3.960032e-06\n",
      "   46|1.608081e-01|-1.675447e-05\n",
      "   47|1.608072e-01|-5.976340e-06\n",
      "   48|1.608046e-01|-1.604130e-05\n",
      "   49|1.608020e-01|-1.617036e-05\n",
      "   50|1.608014e-01|-3.957795e-06\n",
      "   51|1.608011e-01|-1.292411e-06\n",
      "   52|1.607998e-01|-8.431795e-06\n",
      "   53|1.607964e-01|-2.127054e-05\n",
      "   54|1.607947e-01|-1.021878e-05\n",
      "   55|1.607947e-01|-3.560621e-07\n",
      "   56|1.607900e-01|-2.929781e-05\n",
      "   57|1.607890e-01|-5.740229e-06\n",
      "   58|1.607858e-01|-2.039550e-05\n",
      "   59|1.607836e-01|-1.319545e-05\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   60|1.607826e-01|-6.378947e-06\n",
      "   61|1.607808e-01|-1.145102e-05\n",
      "   62|1.607776e-01|-1.941743e-05\n",
      "   63|1.607743e-01|-2.087422e-05\n",
      "   64|1.607741e-01|-1.310249e-06\n",
      "   65|1.607738e-01|-1.682752e-06\n",
      "   66|1.607691e-01|-2.913936e-05\n",
      "   67|1.607671e-01|-1.288855e-05\n",
      "   68|1.607654e-01|-1.002448e-05\n",
      "   69|1.607641e-01|-8.209492e-06\n",
      "   70|1.607632e-01|-5.588467e-06\n",
      "   71|1.607619e-01|-8.050388e-06\n",
      "   72|1.607618e-01|-9.417493e-07\n",
      "   73|1.607598e-01|-1.210509e-05\n",
      "   74|1.607591e-01|-4.392914e-06\n",
      "   75|1.607579e-01|-7.759587e-06\n",
      "   76|1.607574e-01|-2.760280e-06\n",
      "   77|1.607556e-01|-1.146469e-05\n",
      "   78|1.607550e-01|-3.689456e-06\n",
      "   79|1.607550e-01|-4.065631e-08\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "   80|1.607539e-01|-6.555681e-06\n",
      "   81|1.607528e-01|-7.177470e-06\n",
      "   82|1.607527e-01|-5.306068e-07\n",
      "   83|1.607514e-01|-7.816045e-06\n",
      "   84|1.607511e-01|-2.301970e-06\n",
      "   85|1.607504e-01|-4.281072e-06\n",
      "   86|1.607503e-01|-7.821886e-07\n",
      "   87|1.607480e-01|-1.403013e-05\n",
      "   88|1.607480e-01|-1.169298e-08\n",
      "   89|1.607473e-01|-4.235982e-06\n",
      "   90|1.607470e-01|-1.717105e-06\n",
      "   91|1.607470e-01|-6.148402e-09\n",
      "   92|1.607462e-01|-5.396481e-06\n",
      "   93|1.607461e-01|-5.194954e-07\n",
      "   94|1.607450e-01|-6.525707e-06\n",
      "   95|1.607442e-01|-5.332060e-06\n",
      "   96|1.607439e-01|-1.682093e-06\n",
      "   97|1.607437e-01|-1.594796e-06\n",
      "   98|1.607435e-01|-7.923812e-07\n",
      "   99|1.607420e-01|-9.738552e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  100|1.607419e-01|-1.022448e-07\n",
      "  101|1.607419e-01|-4.865999e-07\n",
      "  102|1.607418e-01|-7.092012e-07\n",
      "  103|1.607408e-01|-5.861815e-06\n",
      "  104|1.607402e-01|-3.953266e-06\n",
      "  105|1.607395e-01|-3.969572e-06\n",
      "  106|1.607390e-01|-3.612075e-06\n",
      "  107|1.607377e-01|-7.683735e-06\n",
      "  108|1.607365e-01|-7.777599e-06\n",
      "  109|1.607364e-01|-2.335096e-07\n",
      "  110|1.607364e-01|-4.562036e-07\n",
      "  111|1.607360e-01|-2.089538e-06\n",
      "  112|1.607356e-01|-2.755355e-06\n",
      "  113|1.607349e-01|-4.501960e-06\n",
      "  114|1.607347e-01|-1.160544e-06\n",
      "  115|1.607346e-01|-6.289450e-07\n",
      "  116|1.607345e-01|-2.092146e-07\n",
      "  117|1.607336e-01|-5.990866e-06\n",
      "  118|1.607330e-01|-3.348498e-06\n",
      "  119|1.607328e-01|-1.256222e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  120|1.607320e-01|-5.418353e-06\n",
      "  121|1.607318e-01|-8.296189e-07\n",
      "  122|1.607311e-01|-4.381608e-06\n",
      "  123|1.607310e-01|-8.913901e-07\n",
      "  124|1.607309e-01|-3.808821e-07\n",
      "  125|1.607302e-01|-4.608994e-06\n",
      "  126|1.607294e-01|-5.063777e-06\n",
      "  127|1.607290e-01|-2.532835e-06\n",
      "  128|1.607285e-01|-2.870049e-06\n",
      "  129|1.607284e-01|-4.892812e-07\n",
      "  130|1.607281e-01|-1.760452e-06\n",
      "  131|1.607279e-01|-1.727139e-06\n",
      "  132|1.607275e-01|-2.220706e-06\n",
      "  133|1.607271e-01|-2.516930e-06\n",
      "  134|1.607269e-01|-1.201434e-06\n",
      "  135|1.607269e-01|-2.183459e-09\n",
      "  136|1.607262e-01|-4.223011e-06\n",
      "  137|1.607258e-01|-2.530202e-06\n",
      "  138|1.607258e-01|-1.857260e-07\n",
      "  139|1.607256e-01|-1.401957e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  140|1.607250e-01|-3.242751e-06\n",
      "  141|1.607247e-01|-2.308071e-06\n",
      "  142|1.607247e-01|-4.730700e-08\n",
      "  143|1.607246e-01|-4.240229e-07\n",
      "  144|1.607242e-01|-2.484810e-06\n",
      "  145|1.607238e-01|-2.539206e-06\n",
      "  146|1.607234e-01|-2.535574e-06\n",
      "  147|1.607231e-01|-1.954802e-06\n",
      "  148|1.607228e-01|-1.765447e-06\n",
      "  149|1.607228e-01|-1.620007e-08\n",
      "  150|1.607222e-01|-3.615783e-06\n",
      "  151|1.607222e-01|-8.668516e-08\n",
      "  152|1.607215e-01|-4.000673e-06\n",
      "  153|1.607213e-01|-1.774103e-06\n",
      "  154|1.607213e-01|-6.328834e-09\n",
      "  155|1.607209e-01|-2.418783e-06\n",
      "  156|1.607208e-01|-2.848492e-07\n",
      "  157|1.607207e-01|-8.836043e-07\n",
      "  158|1.607205e-01|-1.192836e-06\n",
      "  159|1.607202e-01|-1.638022e-06\n",
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "  160|1.607202e-01|-3.670914e-08\n",
      "  161|1.607197e-01|-3.153709e-06\n",
      "  162|1.607197e-01|-2.419565e-09\n",
      "  163|1.607194e-01|-2.136882e-06\n",
      "  164|1.607194e-01|-1.173754e-09\n",
      "  165|1.607192e-01|-8.169238e-07\n",
      "  166|1.607191e-01|-9.218755e-07\n",
      "  167|1.607189e-01|-9.459255e-07\n",
      "  168|1.607187e-01|-1.294835e-06\n",
      "  169|1.607186e-01|-5.797668e-07\n",
      "  170|1.607186e-01|-4.706272e-08\n",
      "  171|1.607183e-01|-1.753383e-06\n",
      "  172|1.607183e-01|-1.681573e-07\n",
      "  173|1.607183e-01|-2.563971e-10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XecVdW5//HPM42hF2nSHMFKxEKIYr1G1KjRkBu9xkSvxMY1MeqPGKOGJFfvvcbkRqOG2LAkmhijwRpbriH2CAr2CiiIKL03hynr98ez98xwGJgZZmbvU77v1+u8DqfMOQ8H5jvPrL32WhZCQERE2l9R2gWIiBQKBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCLtzMwONbMP0q5D0qfAlaxnZt8xs7fMbIOZLTKzm8ysR/TYzWa2LrpsMrOqBrefSKC2YGa7bOs5IYTnQwi7b+frzzOzjQ3+TuvM7LfN/NpnzOzs7XlfaR8KXMlqZnYR8EvgYqA7MBrYCXjKzMpCCOeGELqEELoAPwfujW+HEI5Nr3JnZiVt8DInNPg7dQkhfL8NXrOtasva98tGClzJWmbWDbgCOD+E8GQIoSqEMA84GagATtuO1zzczBaY2Y/MbImZLTSzr5vZcWY2y8xWmNmPGzx/fzN7ycxWRc/9rZmVRY89Fz3tjajz/GaD17/EzBYBv4vvi75mWPQeI6PbA8xsqZkdvh1/l++Y2QtmdrWZrTSzuWZ2bPTYlcChwG8bdsVRR36emc0GZkf3HWRmr5jZ6uj6oAbv8YyZXWVmL5vZGjN72Mx6NbO+y81sipn90czWAN8xsyIzu9TMPjSz5WZ2X8PXM7PTzezj6LGfRh3+kS39bLKVAley2UFAOfBAwztDCOuAx4GjtvN1+0evOxD4GXArHt5fxEPqp2a2c/TcGmAC0Bs4EBgDfC+q47DoOftEnee9DV6/F96Jj8+o/UPgEuCPZtYJ+B1wZwjhme38uxwAfBDV97/A7WZmIYSJwPPA9xvpir8efd3wKOweA34D7AD8GnjMzHZo8PzTgTOBHYHq6LnNNRaYAvQA7gbOj97/X4ABwErgBgAzGw7cCJwavVd3/N8obyhwJZv1BpaFEKobeWxh9Pj2qAKuDCFUAX+OXuf6EMLaEMI7wLvAPgAhhJkhhGkhhOqou74FD4ttqQX+M4RQGULYmPlgCOFWYA4wHQ+WiU283kNRhx1fzmnw2MchhFtDCDXAndHr9Wvi9a4KIayIavsqMDuE8Ifo73gP8D5wQoPn/yGE8HYIYT3wU+BkMytu4j1iL4UQHgoh1Ebvdy4wMYSwIIRQCVwOnBQNN5wE/DWE8EIIYRP+wzCvFnsp+DEVyWrLgN5mVtJI6O4YPb49lkcBBRAH4uIGj28EugCY2W541zcK6IR/z8xs4vWXhhA+b+I5twKPAOOj4NmWr4cQ/r6VxxbFfwghbDAz4tq34ZMGfx4AfJzx+Mds3ll+kvFYKf5DquFn1pz3Au/6HzSz2gb31eA/JAY0fH7091nejPfIGepwJZu9BFQC32h4p5l1AY4FpiZQw014x7drCKEb8GPAmviabXZlUf3XAbcDlzd3THQ7bK2Ohvd/hodgQ0OATxvcHpzxWBXN/2GXWcMnwLEhhB4NLuUhhE/x31oGxU80s474MEfeUOBK1gohrMYPmk0ys2PMrNTMKoD7gAXAHxIooyuwBlhnZnsA3814fDEwtIWveT0wI4RwNj5+enOrq2xcc2p7HNjNzL5tZiVm9k1gOPBog+ecZmbDozHn/wKmNPgNoaVuBq40s50AzKyPmY2NHpsCnBAdxCvDhxua+uGWUxS4ktVCCP+Ld5VX48E3He+SxjTjV/G28EPg28BafBjg3ozHLwfujMZWT27qxaJwOYb64P4BMNLMTt3Gl/01Yx7ug82s/Xp8fHSlmTV6oCuEsBw4HrgIWA78CDg+hNCwg/0D8Ht8+KIcuKDB32edmR3azHrimh4B/s/M1gLT8AN4ROPn5+Pj6guBdcAS/LecvGBagFxEtsbMngH+GEK4LYX37gKswodz5ib9/u1BHa6IZA0zO8HMOplZZ/y3mreAeelW1XYUuCKSTcbiB/I+A3YFTgl59Gu4hhRERBKiDldEJCE68UEA6N27d6ioqEi7DJGcMnPmzGUhhD7Nfb4CVwCoqKhgxowZaZchklPMLPMsvW3SkIKISEIUuCKFYuNGeOUVmDMHdLA8FQpckXz32mtwxBHQrRvsvz/suiv06gU/+AGsXZt2dQVFgSuSr2pr4dJLYdQoeOcduPhiuP9+uPVWOPZYuO462GMPePrptCstGDpoJpKPamvhu9+FyZPhrLPgV7+Cnj3rHz/7bLjwQjjjDDjuOHj0URgzJr16C4Q6XJF8dP75HrY//rF3tA3DNnbAAfDssz7EcPzx8NxzWz5H2pQCVyTf3HUX3HgjXHQR/M//gG1jhcM+fWDqVBgyBL75TViyJLk6C5ACVySffPABfO97cNhh8ItfbDtsY336wF/+AitXwrhxPhwh7UKBK5Ivamrg1FOhvBzuvhtKWnCIZu+94dpr4ckn4YYb2q/GAqfAFckXkyfDzJkemIMGNf38TOeeC0cfDT/5CSxuznZl0lIKXJF8sGwZTJzo821PbnLjicaZwaRJfoLEpZe2bX0CKHBF8sPEiX4Sw6RJzRu33ZrddvODbb//Pbz0UpuVJ06BK5LrZs2C227zg2XDh7f+9SZOhP794bLLdApwG1PgiuS6K67wA2UTJ7bN63Xp4vN3n30W/vGPtnlNARS4IrntnXfgnnv8RIe+fdvudc85xw+8/eQn6nLbkAJXJJddcYV3pBdf3LavW17uYTttmk8VkzahwBXJVbNnw5Qp8P3vww47tP3rn3EGDB7sJ1BIm1DgiuSqa66BsjK44IL2ef2yMpgwwddYmD69fd6jwChwRXLR4sU+dWvcOJ9R0F7OPht69PDVxqTVFLgiuWjSJNi0yefMtqeuXX262QMP+BCGtIoCVyTXbNwIN90EY8f6iQrt7fzzobQUfvOb9n+vPKfAFck1f/oTrFjhC4gnoX9/OOUUH8JYsyaZ98xTClyRXBKCDyeMGAH/8i/Jve/558O6dR66st0UuCK55Pnn4Y03PABbs2ZCS40aBQce6GGv9XK3mwJXJJdMmuTb5Zx6avLvfcEFvsW6ToTYbgpckVzx2Wfw4INw5pnQqVPy73/iiT6ee9NNyb93nlDgiuSK22/3XR3+4z/Sef/SUt8B+LHH4OOP06khxylwRXJBdbXv6HDUUb7LblrGj/ex41tvTa+GHKbAFckFjz8OCxbAd7+bbh1DhsBXv+rr727alG4tOUiBK5ILbr4ZBgyA449PuxLf+2zxYnjoobQryTkKXJFsN2+ezww4+2wfR03bV74CO+3kQxzSIgpckWx3220+bnrWWWlX4oqLPfynTvVpYtJsClyRbFZVBXfcAcce6+On2eLMMz14dfCsRRS4Itnsscdg4cL0poJtzYABcMIJ8Lvf6eBZCyhwRbLZLbfAwIHe4Wab8eNh6VIdPGsBBa5Itpo3D/72Nx8vLSlJu5otHX20Dp61kAJXJFvFB8vOPDPtShqng2ctpsAVyUbZerAskw6etYgCVyQbPfpodh4sy6SDZy2iwBXJRtl8sCyTDp41mwJXJNt89JEfLDvnnOw8WJbp6KOhokLLNjaDAlck29xyS/0BqVxQXOxDH888A++/n3Y1WU2BK5JNKiv9YNnXvuZDCrnizDN9nYebb067kqymwBXJJg88AMuWpb8MY0v17es7Qtx5J2zYkHY1WUuBK5JNbrwRhg2DMWPSrqTlzj0XVq2Ce+5Ju5KspcAVyRavvw4vvADnnQdFOfitedhhsNde8Nvf+nbusoUc/FcVyVO//a1vDnnGGWlXsn3M4Pvf9x8c//xn2tVkJQWuSDZYvhzuvhv+/d+hR4+0q9l+p50G3bv7Dw/ZggJXJBvccQd8/rkPJ+Syzp19xsKUKX6mnGxGgSuStupqmDQJDj8cRoxIu5rWO+883879xhvTriTrKHBF0jZlCnzyCfzgB2lX0jaGDYOxY/3MM00R24wCVyRNIcA118Cuu/r24/nioot8XPquu9KuJKsocEXS9OKLMGMGTJiQm1PBtubgg+FLX4Jrr4Xa2rSryRp59C8skoN+9Svo1QtOPz3tStqWmQ+RzJoFjzySdjVZQ4Erkpa33vIwuuACP7qfb046CYYOhauu0okQEQWuSFquugq6dIHzz0+7kvZRUgKXXAIvv+zb8IgCVyQVc+bAvff6IjW9eqVdTfsZN853hfj5z9OuJCsocEXScNVVvpzhhAlpV9K+OnSAH/4Qnn7a14kocApckaTNmuXLGJ57Luy4Y9rVtL/x46FfP/jJTwp+LFeBK5K0yy/3zu+yy9KuJBmdO8PEifDsswU/lqvAFUnSW2/Bn/8MF17oXV+hGD/et3ufOLGgu1wFrkiSLrkEunXzcc1C0qED/OxnPmPh/vvTriY1ClyRpDzxhF9++tP8npmwNePG+eI8F1/sK6MVIAWuSBKqqvzMq112yd95t00pKYHrroN58/yU3wKkwBVJwg03+Bbi11wDZWVpV5OeI47wlcSuvBI+/TTtahKnwBVpb/Pn+5SoY46BE05Iu5r0/frXvqDNeecV3AE0Ba5Iewqhfsvzm2/2RV0K3dChcMUV8PDDvi18AVHgirSne+6Bxx/3X6F32intarLHhAmw336+6eSKFWlXkxgFrkh7mTvXu9sDD/RgkXolJb6P2/LlcPbZBTO0oMAVaQ/V1XDqqf7nP/0JiovTrScb7buvrynx4INwyy1pV5MIBa5Ie7jsMnjpJQ+Sioq0q8leEybAV77i1zNnpl1Nu1PgirS1O++Eq6/24YRTTkm7muxWVOT7nvXt69PF8nxrdQWuSFt6/nlfN2DMGLj++rSryQ19+/rOF6tWeeiuW5d2Re1GgSvSVl5+2Xfe3XlnuO8+X+9WmmeffeDuu31Y4Wtfg40b066oXShwRdrCjBk+Ftmnjy9BWIhrJbTW2LE+HPPMM/CNb8CGDWlX1OYUuCKt9cQTcPjh0L27h+3AgWlXlLtOOw1uvRX+9jcfllm2LO2K2pQCV2R7heCLsZxwAuy2m89K0IyE1jvrLJgyBV5/HUaP9us8ocAV2R5LlvivwBMmwPHH+24GhbBdTlK+8Q34xz98LHf0aJg0yddfyHEKXJGWqKnxX3n32AOefNJnIjz4IHTtmnZl+efAA727PeIIuOACOOQQeOONtKtqFQWuSHPU1PjMg7339mlfe+/t3/wXXKAFadpTnz7w2GPw+9/D7Nm+/sK3vw3vvZd2ZdtFgSuyLXPn+sIzQ4fCN7/p9913n2/7veee6dZWKMx8t4gPPvAtih55BIYPh6OO8n+L9evTrrDZLBTIohGybaNGjQozZsxIu4z0rVnj82mfftqHDF591e8fMwa+9z0ft9W6COlauhQmT/bL/PnQsSMcfTQceSQcdpiHcUlJIqWY2cwQwqhmP1+BK1BAgVtZ6csBLl4Mixb5N+zcud49vfMOzJrlzyspgf33h3/9VzjxRD+ZQbJLTY2f2Tdlii+BOXeu39+xow/57Lkn7LqrL4s5eLDvktynj0/fa6Mfmgpc2S6bBe5DD8Gjj7b/m27t/158f8PrrV1qavzodXW1/7mqyi+VlX7ZuNEva9f6pbEzmEpKfK+xPfeEkSPhS1+Cgw7SgbBcM3cu/POf8Morvh39e+81vjaDmf/bdu0KnTtDp05QXu47C5eV+RmCJSV+KS729R4OO8x/w9nipVoWuMn03ZJb5szxX6eTsLUDTvH9Da8bu8TfEPE3R2mpXzp0gC5dvNvp2LH+G6xnTz8LrG9f73iGDIEBAzRMkA923tkv8bKY4OO78+f7/mmLF/uJFCtWwOrV/gN4/Xo/o+3zz/0H9Nq1/gO7utovtbV+aaPF49XhClBAQwoibailHa5mKYiIJESBKyKSEA0pCABmthT4uMFdvYFsXjkkm+vL5togu+vL5tpgy/p2CiH0ae4XK3ClUWY2oyVjU0nL5vqyuTbI7vqyuTZofX0aUhARSYgCV0QkIQpc2ZrJaRfQhGyuL5trg+yuL5trg1bWpzFcEZGEqMMVEUmIAldEJCEKXNmCmR1jZh+Y2RwzuzTlWgab2dNm9q6ZvWNmF0b39zKzp8xsdnTdM8Uai83sNTN7NLq9s5lNjz6/e82sLMXaepjZFDN738zeM7MDs+yzmxD9u75tZveYWXman5+Z3WFmS8zs7Qb3Nfp5mftNVOebZjayqddX4MpmzKwYuAE4FhgOfMvMhqdYUjVwUQhhODAaOC+q51JgaghhV2BqdDstFwINtyD4JXBtCGEXYCVwVipVueuBJ0MIewD74HVmxWdnZgOBC4BRIYS9gGLgFNL9/H4PHJNx39Y+r2OBXaPLeOCmJl89hKCLLnUX4EDgbw1uXwZclnZdDep5GDgK+ADYMbpvR+CDlOoZFH0THgE8Chh+JlJJY59nwrV1B+YSHRxvcH+2fHYDgU+AXvjKhY8CX0n78wMqgLeb+ryAW4BvNfa8rV3U4Uqm+JsgtiC6L3VmVgHsB0wH+oUQ4sVOFwH9UirrOuBHQLyl7A7AqhBCdXQ7zc9vZ2Ap8LtoyOM2M+tMlnx2IYRPgauB+cBCYDUwk+z5/GJb+7xa/L2iwJWcYGZdgPuB/xdCWNPwseDtReLzG83seGBJCGFm0u/dTCXASOCmEMJ+wHoyhg/S+uwAorHQsfgPhgFAZ7b8dT6rtPbzUuBKpk+BwQ1uD4ruS42ZleJhe3cI4YHo7sVmtmP0+I7AkhRKOxj4mpnNA/6MDytcD/Qws3hx/zQ/vwXAghDC9Oj2FDyAs+GzAzgSmBtCWBpCqAIewD/TbPn8Ylv7vFr8vaLAlUyvALtGR4rL8IMYj6RVjJkZcDvwXgjh1w0eegQYF/15HD62m6gQwmUhhEEhhAr8c/pHCOFU4GngpDRri+pbBHxiZrtHd40B3iULPrvIfGC0mXWK/p3j+rLi82tga5/XI8Dp0WyF0cDqBkMPjUtjsFyX7L4AxwGzgA+BiSnXcgj+K9ybwOvR5Th8rHQqMBv4O9Ar5ToPBx6N/jwUeBmYA/wF6JBiXfsCM6LP7yGgZzZ9dsAVwPvA28AfgA5pfn7APfh4chX+G8JZW/u88AOkN0TfJ2/hsy22+fqtOrXXzI7Bf4UqBm4LIfxiu19MRCTPbXfgRvM1Z+FTdBbgv4p+K4TwbtuVJyKSP1qza+/+wJwQwkcAZvZn/IjjVgO3d+/eoaKiohVvKa21apVvTDp48Ob3z5w5c1mIVq4/qujftKKRSIanav+ylS2mm681gdvYHLQDMp9kZuPxszAYMmQI2hk2XT/6EUyaBJn/DGb2ceNfISJtpd1nKYQQJocQRoUQRvXp0+ytf6SdVFVBaWnaVYgUptYEbtbN15SmbdoEZaktpSJS2FoTuFk1X1OaZ8MG6Nw57SpECtN2j+GGEKrN7PvA3/BpYXeEEN5ps8qkXaxZA126pF2FSGFqzUEzQgiPA4+3US2SgOXLYYcd0q5CpDDp1N4C8+mn0L9/2lWIFCYFbgGpqYF582DYsLQrSYFZyy4i7UCBW0Defhuqq2GvvdKuRKQwtWoMV3LLs8/69WGHpVtHVrAmeo24yQ21m9/firVHRNThFpA//hGGD9/ytF4RSYY63AIxbRq88oqf1luQ4s40Hp+NO9eo07Ui2+w28e3ajI42+roQ368OWFpAHW4B2LQJvvtd6NcPTj897WpECpc63DwXAkycCK+/Dg89BN26pV1RyjI73UxRZ2vx46XFjT7fotcJNTV+R23GbXW+0gh1uHksBF8d7OqrvcMdOzbtikQKmzrcPLV6NVxwAdx1F5x3HvzmN2lXlGXqOs5oTDZqTOsmJxR7Zxt3uhbdpu46Gvut2byTDZs2+R+iTjfEj8djv3EHvEUdUgjU4eahhx7y2Qh//CP87Gd+oKxI/9IiqVOHm0deegmuugr++lfYe28P3i99Ke2qslxGhxl3oHUjtlGHG4o2n81gJdG3TnnJZs+zomhloErvdENVld+uqvbb1dF13e2qbdYj+UV9T46rqYEHH4SDD4aDDoIXXoCf/9x3dFDYimQXdbg56r334N574e67Yc4cqKjwcdozztDyi62SMfsgnm9rmZ1nUTTGWxp1wOW+qnvo4NfWpdPmz9/knWzR51HnG431hnXr/Tqj862b5aCON68ocHPIrFkesvfd5+simPlpuldeCd/4BpToX1Mkq+lbNIutWwfPPw9//zs89RS89ZaH7CGH+IGwE0+EHXdMu8o8FXeW0fSFEA+1xrMM6ubfZoz5lnjnW9O13B/u4N9ioSQa46321y1e5x1u0fqufl2Z0fFu/NxfLxobrt2wwW9GsyTijlhyiwI3i1RVwfTpHrBTp/rpuNXVvgfZwQfDddfBSSfBwIFpVyoi20OBm6IlSzxg48tLL8H69d7UfPGLcNFFcOSRHrYdO6ZdbYGLOtkQol42nm9b9/jmY71F0RlqcYdb3dFvV3WJxn5rojHfaFZDh5XesZatqvSvX+WdLmv9Oj66HXe2lnF7C3VrRmgMOJsocBNSWemn106b5uE6bRrMneuPFRfDiBEwbhyMGQOHHw69eqVaroi0AwVuO1i3Dt580wP2tdf8+s0365uigQNh9Gg/3faAA7yb1U66OSJzDYV4tkE0u6Aouj++Lq2MO1Afq62NZjVs6O2dbnUnv71uoN8uXV8WXXvn23GJDx53WOKdbvGKtf5+Gzf668VjvpmzG+JVz0LGmW2SKgVuKy1evHmwvvYazJ5d/5tcr16w335w4YUergccAIMGpVuziKRDgdtMGzfCu+/6TIGGl0WL6p9TUeHheuqpsO++/udBg7RFVl6KO914DDX6R65ZF81aqIzGYj/367JoFkJRZXf/uiKfp7sx6kTXD/Kv39jPX67WG11K1nUAoNMiv6PLQl/urcMyf73SRav9iSv9Oqz32QzxGg5Bq5ZlFQVuhpoa+OijLYN1zhyojf7vlpf7WgVf+Up9sO6zD/TokW7tIpLdCjpwlyzZMljfeQeiKY+Y+Q63I0bAKaf49YgRsMsu9YtGiQBbnbdbu9bHXC0a6y1e72OvXVd7p1s62K+Lqr2DXbOzd7pVvaPOub93yKuH+s1Vq/15Zcu9Q+78mU9f6flBT79/yTp/v+Wr/AuiTrs2et+6tRviMd5ajfEmqSACd8MGD9LMcF2ypP45ffp4mJ5zTn2wfuELOpglIm0n7wJ3xQo/cPXqq/XXs2bVNyDl5b5N+Fe/Wh+sI0b49jMibSZj3m6IzxSL59F+7meSdYzWWChb6bMS4tkJ61aVArB2mHeinQZ457rLkM8A2FTrv2LN+sz/467a08d6u8zrDUC3+d7xdp7rX1e80jvt2uUrvb54rYh41bKqjHnF0i5yNnBDgM8+2zJc58+vf87gwTBy5ObDAcOGaThARNKRM4FbW+vDAs895+sLPPccLFxY//huu8GBB/ruBvvt55fevdOrVwpc5tSUeFZDNKYa7yRRu3AxAEWr1gDQfY13pp2i2Qjly32Mdtm+Po/3vSrvFnbv7+Nh/77XdAB6lfh83EcW7Q3AnA/7A9DtPR8j7jYv6qDX9vHXnbfC61iy3OuIZzPEazVEdUrbytrAraqCmTPrw/XFF2Fl9NvQwIF+Ntbo0d7B7rMPdO2aarkiIk3KqsANwRfQvusuX4Jwjf/QZ7fdfPnBQw/15QgrKjS3VbJcE/Nda+MOMpotUDe2G12XrvGx1x1WesfbYY13qut29KO4b40YAsCqYd4Bn7nTiwDcsMufAfikwjvkP+19IAAvzt8ZgOLXvTPp2aMvAF0+9s63+DPvdMN675Rro7Uh4l2M685k06yGVmkycM1sMHAX0A8IwOQQwvVm1gu4F6gA5gEnhxBWbk8RH30Ef/iDB+1HH/nMgBNPhBNO8KUI+/ffnlcVEckuzelwq4GLQgivmllXYKaZPQV8B5gaQviFmV0KXApc0tICrrkGfvhD71jHjIHLL/duVtOxJK9lztuNN3iIzhSzaK0E2xDN213lsww69/TOtesCv17yRV8Q+b8WHw/Av37hdQAO7TYLgNuHvADAm/3/DsAdww4B4PFZXwBgxdve8XZa7J1u7zeiWQ1LfB5vWO2/ZgaL5vNWaieK1mgycEMIC4GF0Z/Xmtl7wEBgLHB49LQ7gWdoYeBOnuxhe+KJcO21PqtARCRfWWjBTyozqwCeA/YC5ocQekT3G7Ayvr01o0aNCjNmzAB8acJhw3ylrBdf9EW2JT1mNjOEMArgqKJ/U/uSLeLdgKPZA8X9fOw13ve+aohPxdnYz+fhrtjDn7dxD5/ne+a+/wTgy13eBWB4qd9fhf8TX75oDADTFu4EQOVLOwDQ633vvDsu8ueXLo7WaojWbIjn8cY7UdStnpbHne9TtX9p9ZGjZu/aa2ZdgPuB/xdCWNPwseCp3egnbWbjzWyGmc1YunRp3f2DBsGXv+zzZx9/fPuKFxHJJc3qcM2sFHgU+FsI4dfRfR8Ah4cQFprZjsAzIYTdt/U6DTtcgLVr4aij4JVXfCGYceNg7Fg/G0ySpQ43y0WdblEH72Tj+bJWFu8S7Ac9avr7rIZ1Q3ythbjjtZHemR4/9B0Aju32JgCHd/Qx2f/b4Ge2vbR+VwDunT3S32a6jxXv8K6fEddxQbRWw6JoVkM0xhyiVdHq1mrIw043kQ43Gi64HXgvDtvII8C46M/jgIdb+uZdu8KTT8Ill/jaBqec4jMSxo/3YYba2qZfQ0QkVzTZ4ZrZIcDzwFtAHIE/BqYD9wFDgI/xaWErtvVamR1uQzU18MwzcOedcP/9vuBMz54+Leyww3wO7siRUFrakr+eNJc63BwRT0C3zXulovLNO1+G+ir38Z5qa4d6B7x6aDTvd5R3vHv189M1j9vBO94vlfu58U+s28uvF/tshg/fHwBA93f99XvM8bUXOs7z2Qy2OlqrIV6PN5plEWpD/RSMWI52v23R4TZnlsILNNgFOsOY1hZysepGAAANoElEQVQQKy72aWFjxsCNN8LDD8PTT/uZZn/9qz+nUyc/fTc+AWL//TV9TERyR4tmKbTWtjrcbVm0yM9Ai9dReOMN/yFpBrvv7usmjBxZv4aCNmBsOXW4OWZru/LGY72dfAzXos7Xots1/Xwi0YZBfnvpvt5zfd7fzySrGOZrOxzV731/XrT1xNJNPk/3lUV+htv6N/ybrFu0EWqPWd7Rli6PdhtetMzL27ixbleMEJ+9lqNnqyXS4WaD/v3hpJP8ArBqFfzzn36w7bXXPIzvuaf++TvtVB/AI0f6rgwDBuh0YBFJV04EbqYePeC44/wSW7bMwze+vPoqPPRQfQPQs2f9Eo177+3Xe+2lRW8kR8X/sTM73ei6NloTId6+xKLdfYuinSe6LvCv6/Spz+v9vI9PDVoz2MdqfzfMr0t2jtZ06OpfP7Snz06YP9LHZZf29lkRlT28Y+42P9qDLbpdung1YcXKqKZoXDdknK1WlLFeao52wM2Rk4HbmN69fYrZUUfV37d2rQ8/vPGGz4J4801fryHa9QTwhXAaLkQ+YoQvlqODc5ITmhoSjA6uhShoa5Z6YFqpf+vbmuiU4W7RKb4f+nWP2b5YzsrdfShhZU+//7MdfXnH0r7Rqcedfbhg7Z5eR3Un/8ap7BoFcKcSOnTw+4pWRCdPrPUQr1v8PJ5KVlfzVoZL8kDeBG5junb1WQ6HHFJ/Xwjw8cdbbrfzxBMQb8BaVgZ77LFlEGsHXhFpjbwO3MaYeVdbUeGrkcUqK+GDD7wLjkP42Wfh7rvrn9Ojhw9DZA5NdOuW9N9CpJkyfz2PF8upjK4zO4hV3oWWR9u69//EI6JqgHe8lT28W13fz6cHbRgQbXrZ1YcJqrt4V7qhf7TAemkZXTr6sEPHLj5sUbzMl5QMK31KWe3GjCUg4yGHPOx0Cy5wt6ZDBw/Qvffe/P6VK+Httzfvhu++u36tXvBdfOPt0uPr/v3VDYvI5hS4TejZ0+f9Hnpo/X0hwCefeDf8+ut+efVVmDKl/jl9+9YH8L77+iI9u+yiEJYsEx9kiw+ulfo0sOr5n0a3PSLK1vnjZSV+u/MO/mtdZW/vVjf18Ps39vL/4NU+hEtNuVHZw8eRa0v9zvJO/tzScn+v4mjpydpoPLluGtmm/NvYUoG7HcxgyBC/HH98/f2rV3sIv/aah/Brr8Gvf+3bBYHPDz7gAL+MHu0nbvTsmc7fQUSSp8BtQ927b9kNb9oE774LM2bA9OkwbZqvHxEPS+22m4fvAQfAQQf5kEZRs9dwE2lbmdulh6hZqF7sK/0VlUXTd6Kx3o6LfSy3vLN3ul26exe7qaeP11Z3KaaqYzRTIpr9VdU5ip2+0caW0X/4+LXrFz2PxoE3xQvi5P7i5wrcdlZWVj+scPbZft+aNR7A06Z5CD/5pE9XA5/e9uUvw5FH+mnOQ4dqGEIkXyhwU9CtGxxxhF+gfqrac8/B1Kl++ctf/LGKivo1Jo45RkMQkrCMWQ61n0e3o64zXqTGVvt4bNFK72w7Lo5mInQqp6aXd8G1pdEW7CVRB1ETjR93ir42eg+LTwGOFuIpKo6Wfozfs26x89zreBW4WaDhVLXTT/f/Px98UB++998Pt9/uJ2McfTScfLKvG9y9e9qVi0hLKHCzkJmfeLHHHnDeeb505SuvePDedx889pgPVRxzjIfviSdq0XZJWG3G5pdx1xktRF4ULdNIcTEla3yslmiGQ+gYLSUZj9kWbz7fNkSzF+pG0uIxtfjgRjR7oW7eLrnT6erwTA4oLvYDa7/6FcybBy+95EE8cyacdpov1vPf/w3Ll6ddqYhsS04szyiNq631NYOvucZPTe7YEc48Ey6+2EO4JbQ8o7SJjCO8Vlxcv1h6UbR0ZMfo17F4m6B44ZKoA97iDLN4Xm40vzIey43nW9bN223njSwT3URSsk9RkR9Me/zx+i2KJk+G4cPhuuvqNlYVkSyhwM0Te+0Fd9wBs2fD4YfDhAk+r/edd9KuTApKCJtdQnU1obrKVwSrqYGaGmrWrKNmzTpq48uq1X5ZuYralasIa9f6Zf0Gv1RXexdbG6A2YGaYmXfIDS4WXTDL2rmUCtw8s9NO8Oij8Kc/wUcfeei++GLaVYkIKHDzkhl861u+vkP//j6VbOrUtKuSgrVZt1vtMxxqawhVmwhVm6jdVOWXjZ/7Zf1Gatdv9O15Nm4krFvvl8pKv9TU+Hht9Lp1Ha8VgRVt2elmUcerwM1jgwf7yRRDh/r2RJ9+mnZFIoVNgZvn+vWDBx/09X7POScnpipKoYk63rrONxrzra2s9EvUAYdNm6JLlV+qqv0SAo3Otoo63vrb6Xe6CtwCsMsu8POf+9Sxp59OuxqRwqXALRDnnuvLQ954Y9qViDQhY6ZDXecbjd3GHXDdJe5048drQ/2W7FlGgVsgysth3Dh4+GGIN3QVkWQpcAvI0Uf7STvTpqVdich2yOx8MzpgQm3jl0wpjuUqcAvIgQf69csvp1uHSKHSamEFpHt3n7Xw4YdpVyLSDnJgCo463AJTUQHz56ddhUhhUuAWmN69tYyjSFoUuAWmZ09YuTLtKkQKU7MD18yKzew1M3s0ur2zmU03szlmdq+ZlbVfmdJWOnWCaBsqEUlYSzrcC4H3Gtz+JXBtCGEXYCVwVlsWJu2jvFyBK5KWZgWumQ0CvgrcFt024AhgSvSUO4Gvt0eB0rZKS+sWyheRhDW3w70O+BF1u7WxA7AqhBDv4rYAGNjYF5rZeDObYWYzli5d2qpipfVKSup2LBGRhDUZuGZ2PLAkhDBze94ghDA5hDAqhDCqT58+2/MS0oaKinwvNBFJXnNOfDgY+JqZHQeUA92A64EeZlYSdbmDAK22mgMUuCLpabLDDSFcFkIYFEKoAE4B/hFCOBV4Gjgpeto44OF2q1LajFlOnJAjkpdaMw/3EuAHZjYHH9O9vW1KkvakwBVJT4vWUgghPAM8E/35I2D/ti9J2lOWbO0kUpB0ppmISEIUuAVGHa5IehS4IiIJUeCKiCREgSsikhAFrohIQhS4IiIJUeCKiCREgSsikhAFrohIQhS4IiIJUeCKiCREgSsikhAFrohIQhS4IiIJUeCKiCREgSsikhAFrohIQhS4IiIJUeCKiCREgSsikhAFrohIQhS4IiIJUeCKiCREgSsikhAFrohIQhS4IiIJUeCKiCREgSsikhAFrohIQhS4IiIJaVbgmlkPM5tiZu+b2XtmdqCZ9TKzp8xsdnTds72LFRHJZc3tcK8Hngwh7AHsA7wHXApMDSHsCkyNbouIyFY0Gbhm1h04DLgdIISwKYSwChgL3Bk97U7g6+1VpIhIPmhOh7szsBT4nZm9Zma3mVlnoF8IYWH0nEVAv8a+2MzGm9kMM5uxdOnStqlaRCQHNSdwS4CRwE0hhP2A9WQMH4QQAhAa++IQwuQQwqgQwqg+ffq0tl4RkZzVnMBdACwIIUyPbk/BA3ixme0IEF0vaZ8SRUTyQ5OBG0JYBHxiZrtHd40B3gUeAcZF940DHm6XCkVE8kRJM593PnC3mZUBHwFn4GF9n5mdBXwMnNw+JYqI5IdmBW4I4XVgVCMPjWnbckRE8pfONBMRSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIc0KXDObYGbvmNnbZnaPmZWb2c5mNt3M5pjZvWZW1t7FiojksiYD18wGAhcAo0IIewHFwCnAL4FrQwi7ACuBs9qzUBGRXNfcIYUSoKOZlQCdgIXAEcCU6PE7ga+3fXkiIvmjycANIXwKXA3Mx4N2NTATWBVCqI6etgAY2NjXm9l4M5thZjOWLl3aNlWLiOSg5gwp9ATGAjsDA4DOwDHNfYMQwuQQwqgQwqg+ffpsd6EiIrmuOUMKRwJzQwhLQwhVwAPAwUCPaIgBYBDwaTvVKCKSF5oTuPOB0WbWycwMGAO8CzwNnBQ9ZxzwcPuUKCKSH5ozhjsdPzj2KvBW9DWTgUuAH5jZHGAH4PZ2rFNEJOeVNP0UCCH8J/CfGXd/BOzf5hWJiOQpnWkmIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIFbYEaPhgsvTLsKkcJkIYTk3sxsKfBxYm8oLbFTCKFP2kWI5LNEA1dEpJBpSEFEJCEKXBGRhChwRUQSosAVEUmIAldEJCEKXBGRhChwRUQSosAVEUmIAldEJCEKXBGRhChwRUQSosAVEUmIAldEJCEKXBGRhChwRUQSosAVEUmIAldEJCEKXBGRhChwRUQSosAVEUmIAldEJCEKXBGRhPx/pesvP/7u8tEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2c0386ec50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Example with entropic regularization\n",
    "\n",
    "\n",
    "def f(G):\n",
    "    return np.sum(G * np.log(G))\n",
    "\n",
    "\n",
    "def df(G):\n",
    "    return np.log(G) + 1.\n",
    "\n",
    "\n",
    "reg = 1e-3\n",
    "\n",
    "Ge = ot.optim.cg(a, b, M, reg, f, df, verbose=True)\n",
    "\n",
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Ge, 'OT matrix Entrop. reg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD with Frobenius norm + entropic regularization\n",
    "-------------------------------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Loss        |Delta loss\n",
      "--------------------------------\n",
      "    0|1.693084e-01|0.000000e+00\n",
      "    1|1.610121e-01|-5.152589e-02\n",
      "    2|1.609378e-01|-4.622297e-04\n",
      "    3|1.609284e-01|-5.830043e-05\n",
      "    4|1.609284e-01|-1.111407e-12\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAAFgCAYAAAD3rsH6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XecXVW5//HPMy2TXsiQkDqEjnQjBEEuErpgvICKFyXSclEEbkQEBL1wVbCAlCglAgJSBEI1tB9GigoJJPSAJCE9QArpfcr6/fGsk5lMJpm+9zlzvu/Xa7/OnH3O7L3mzJzvPGfttde2EAIiItL2CtJugIhIvlDgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIEr0gAzu9XMfpp2O9qSma02syEpt+EuM/tFmm1oawpcaRYz+66ZvWtma83sUzO7xcx6xMdujW/g1Wa20cwqat1/ppXbcbiZzW/NbdYVQjg3hPDzttxHWzGzF83s7IaeF0LoEkKY2Yzt3xV/x6trLd9sXmvbPwWuNJmZXQT8GrgY6A4MAwYDz5tZSQyoLiGELsDVwIOZ+yGE41Job1HS+6ynDVea2ZVpt6OuVnptflPr99slhPBgG+2n0bLhd14fBa40iZl1A64Czg8hPBtCqAghzAa+AZQD327mdk8ws7fMbLmZvWJm+9R6bLaZ/cjM3jGzFWb2oJmVmlln4BmgX63qql8Mt3Fmdq+ZrQS+a2YdzOwGM/s4LjeYWYe4/cPNbL6Z/cTMlsT9nVZr/5t91DWzEbGtK83sIzM7tjk/8zZei2Bm3zez6Wa2ysx+bmY7xddlpZk9ZGYl8bk9zWy8mS02s2Xx6wHxsV8CXwJ+H1+b39fa/nlmNh2YXmvdzmZWEn+28+P6QjP7l5n9rBk/x2wzu8TM3gHWmFmRme0Rq+7lZjbVzL5a59t6m9nz8ed+ycwGN3Jf343tvN7MPgOujOvPNLMP4mvzXO3tmdnRZvZh/Ju6Oe6vwU8DLRJC0KKl0QtwLFAJFNXz2N3AA3XWXQnc28A29wcWAQcBhcBIYDbQIT4+G3gN6Af0Aj4Azo2PHQ7Mr2efFcDX8KKiI/B/wERge6AMeAX4ea1tVAK/AzoA/wGsAXaLj98F/CJ+fSCwAjgqbrs/sHsjXrcrgSsb+RoH4AmgG/A5YAMwARiCf6J4HxgZn7sdcDLQCegKPAw8XmtbLwJn17P95+Nr2bHWup3j13sBy4A9gMvj61a4lbZuem3qeWw28BYwMP4OioEZwE+AEuAIYFWd13kVcFj8PdwI/LORr9l34+/wfKAo7m9E3N8ecd0VwCvx+b2BlcBJ8bEL49/M2Y3ZX3MXVbjSVL2BJSGEynoe+yQ+3lSjgNtCCJNCCFUhhLvxkBlW6zk3hRA+DiEsBf4K7NfANl8NITweQqgOIawDTgP+L4SwKISwGK/Sv1Pne34aQtgQQngJeAqv2us6C7gzhPB83PaCEMK/m/4jN+g3IYSVIYSpwHvA/wshzAwhrMCr+v0BQgifhRAeCSGsDSGsAn6J/8NoyDUhhKXxtdlMCOE94BfA48CPgO+EEKq2sa0fxYp1uZktqfPYTSGEeXE/w4AuwK9CCBtDCH8HxgPfqvX8p0IIL4cQNuBhf7CZDWzEzwPwcQhhTAihMu7v3PhzfhD/Xq8G9otV7vHA1BDCo/Gxm4BPG7mfZlPgSlMtwT/21ddHtkN8vKkGAxfVetMux6uifrWeU/vNsBZ/427LvDr3+wFzat2fU2f7y0IIa7bxeMZA4KMG9g1A/Hif+XkuBS6t9TOOb+DbF9b6el0997vEfXQys9vMbE7sPnkZ6GFmhQ1sv+7rU9fd+O/l6RDC9Aaee20IoUdc6v7Drb2ffsC8EEJ1rXVz8E8JWzw/hLAaWEr9v4f61P2ZBgM31vodLAUs7q9fnX0FoE0PvoICV5ruVbz6PKn2SjPrAhyHf/RtqnnAL2u9aXuEEDqFEB5oxPdubbq7uus/xt+AGYPiuoyesU94a4/XbutOjWgXIYQTMj8P8Cu8ssv8fCc0ZhuNcBGwG3BQCKEb/nEcPFig8a9PXTfj1ecxZnZoC9pXez8fAwPNrHbuDAIW1Lq/qZqNf1O9qP/30NC+wH9X/13n76pjCOEV/NPYgFr7str324oCV5okfqS9ChhjZseaWbGZlQMP4RXCn5ux2T8C55rZQeY6m9lXzKxrI753IbCdmXVv4HkPAFeYWZmZ9QZ+Btxb5zlXxYNGXwJOwPtD67oDOMPMhptZgZn1N7PdG9HOttIVr3iXm1kv4H/rPL4Q7/ttNDP7DvB5vF/0AuDuGH4tNQn/dPLj+HdzOHAi8JdazznezA6NBwV/DkwMITRUjW/NrcBlZvY5ADPrbmZfj489BextZl+Ln9bOA/o2cz+NpsCVJgsh/AY/8HEtfuBhEl5NDI99b03d3mTgHOD3+MGaGfibvTHf+288TGfGj45b+/j5C2Ay8A7wLvBGXJfxadz3x8B9+EG5LfpmQwivAWcA1+MHz15i88o5aTfgB4iW4Ae3nq3z+I3AKfEo/U0NbczMBsVtnh5CWB1CuB9/3a5vaUNDCBvxgD0utvfmuJ/ar/P9+D+NpXjobxr1Ekc1nEYjhRAew4cv/iV2t7wX900IYQnwdeA3wGfAnvjP2eS/36aweMROJG/FSuveEEKbf6SU7BS7OeYDp4UQXmir/ajCFZG8ZGbHmFkP8/HYP8H7vSe25T4VuCKSrw7GR5wswbs6vlbfMLnWpC4FEZGEqMIVEUlIVk7wIMnr3bt3KC8vT7sZIjllypQpS0IIZY19vgJXACgvL2fy5MlpN0Mkp5jZnIafVUNdCiIiCVHgiuSLdevg9ddhxgzQwfJUKHBF2rs334QjjoBu3eDAA2GXXaBXL/jhD2HVqrRbl1cUuCLtVXU1XHopDB0KU6fCxRfDI4/AH/8Ixx0HN9wAu+8OL7TZiVVShw6aibRH1dXwve/B2LFw1lnw299Cz541j599Nlx4IZxxBhx/PIwfD8OHp9fePKEKV6Q9Ov98D9uf/MQr2tphm3HQQfDSS97FcMIJ8PLLybczzyhwRdqbe+6Bm2+Giy6CX/wCzLb+3LIymDABBg2Cb34TFi1Krp15SIEr0p58+CF8//tw2GHwq19tO2wzysrg4Ydh2TIYOdK7I6RNKHBF2ouqKjjtNCgthfvug6ImHKLZZx+4/np49ln4wx/aro15ToEr0l6MHQtTpnhgDmjG1L7nngtHHw1XXAELFzb8fGkyBa5Ie7BkCVx+uY+3/UZ9FxtuBDMYM8ZPkLj00tZtnwAKXJH24fLL/SSGMWMa12+7Nbvu6gfb7roLXn211ZonToErkuumTYPbb/eDZXvu2fLtXX459O0Ll12mU4BbmQJXJNdddZUfKLv88tbZXpcuPn73pZfg739vnW0KoMAVyW1Tp8IDD/iJDttv33rbPeccP/B2xRWqcluRAlckl111lVekF1/cutstLfWwnTjRh4pJq1DgiuSq6dNh3Dj4wQ9gu+1af/tnnAEDB/oJFNIqFLgiueq666CkBC64oG22X1ICo0f7HAuTJrXNPvKMAlckFy1c6EO3Ro70EQVt5eyzoUcPn21MWkyBK5KLxoyBjRt9zGxb6trVh5s9+qh3YUiLKHBFcs26dXDLLTBihJ+o0NbOPx+Ki+Gmm9p+X+2cAlck19x/Pyxd6hOIJ6FvXzj1VO/CWLkymX22UwpckVwSgncn7L03/Md/JLff88+H1as9dKXZFLgiueQf/4C33/YAbMmcCU01dCgcfLCHvebLbTYFrkguGTPGL5dz2mnJ7/uCC/wS6zoRotkUuCK54uOP4bHH4MwzoVOn5Pd/8snen3vLLcnvu51Q4Irkijvu8Ks6/Pd/p7P/4mK/AvBTT8GcOem0IccpcEVyQWWlX9HhqKP8KrtpGTXK+47/+Mf02pDDFLgiueDpp2H+fPje99Jtx6BB8JWv+Py7Gzem25YcpMAVyQW33gr9+sEJJ6TdEr/22cKF8Pjjabck5yhwRbLd7Nk+MuDss70fNW3HHAODB3sXhzSJAlck291+u/ebnnVW2i1xhYUe/hMm+DAxaTQFrkg2q6iAO++E447z/tNsceaZHrw6eNYkClyRbPbUU/DJJ+kNBduafv3gxBPhT3/SwbMmUOCKZLPbboP+/b3CzTajRsHixTp41gQKXJFsNXs2PPec95cWFaXdmi0dfbQOnjWRAlckW2UOlp15ZtotqZ8OnjWZAlckG2XrwbK6dPCsSRS4Itlo/PjsPFhWlw6eNYkCVyQbZfPBsrp08KzRFLgi2WbmTD9Yds452XmwrK6jj4byck3b2AgKXJFsc9ttNQekckFhoXd9vPgi/PvfabcmqylwRbLJhg1+sOyrX/UuhVxx5pk+z8Ott6bdkqymwBXJJo8+CkuWpD8NY1Ntv71fEeLuu2Ht2rRbk7UUuCLZ5OabYaedYPjwtFvSdOeeC8uXwwMPpN2SrKXAFckWb70F//wnnHceFOTgW/Oww2CvveD3v/fLucsWcvC3KtJO/f73fnHIM85IuyXNYwY/+IH/43jllbRbk5UUuCLZ4LPP4L774DvfgR490m5N833729C9u//zkC0ocEWywZ13wvr13p2Qyzp39hEL48b5mXKyGQWuSNoqK2HMGDj8cNh777Rb03LnneeXc7/55rRbknUUuCJpGzcO5s2DH/4w7Za0jp12ghEj/MwzDRHbjAJXJE0hwHXXwS67+OXH24uLLvJ+6XvuSbslWUWBK5Kmf/0LJk+G0aNzcyjY1hxyCHzhC3D99VBdnXZrskY7+g2L5KDf/hZ69YLTT0+7Ja3LzLtIpk2DJ59MuzVZQ4ErkpZ33/UwuuACP7rf3pxyCgwZAtdcoxMhIgWuSFquuQa6dIHzz0+7JW2jqAguuQRee80vwyMKXJFUzJgBDz7ok9T06pV2a9rOyJF+VYirr067JVlBgSuShmuu8ekMR49OuyVtq0MH+NGP4IUXfJ6IPKfAFUnatGk+jeG558IOO6TdmrY3ahT06QNXXJH3fbkKXJGkXXmlV36XXZZ2S5LRuTNcfjm89FLe9+UqcEWS9O678Je/wIUXetWXL0aN8su9X355Xle5ClyRJF1yCXTr5v2a+aRDB/jZz3zEwiOPpN2a1ChwRZLyzDO+/PSn7XtkwtaMHOmT81x8sc+MlocUuCJJqKjwM6923rn9jrttSFER3HADzJ7tp/zmIQWuSBL+8Ae/hPh110FJSdqtSc8RR/hMYr/8JSxYkHZrEqfAFWlrc+f6kKhjj4UTT0y7Nen73e98Qpvzzsu7A2gKXJG2FELNJc9vvdUndcl3Q4bAVVfBE0/4ZeHziAJXpC098AA8/bR/hB48OO3WZI/Ro2H//f2ik0uXpt2axChwRdrKrFle3R58sAeL1Cgq8uu4ffYZnH123nQtKHBF2kJlJZx2mn99//1QWJhue7LRfvv5nBKPPQa33ZZ2axKhwBVpC5ddBq++6kFSXp52a7LX6NFwzDF+O2VK2q1pcwpckdZ2991w7bXenXDqqWm3JrsVFPh1z7bf3oeLtfNLqytwRVrTP/7h8wYMHw433ph2a3LD9tv7lS+WL/fQXb067Ra1GQWuSGt57TW/8u6OO8JDD/l8t9I4++4L993n3Qpf/SqsW5d2i9qEAlekNUye7H2RZWU+BWE+zpXQUiNGeHfMiy/CSSfB2rVpt6jVKXBFWuqZZ+Dww6F7dw/b/v3TblHu+va34Y9/hOee826ZJUvSblGrUuCKNFcIPhnLiSfCrrv6qASNSGi5s86CcePgrbdg2DC/bScUuCLNsWiRfwQePRpOOMGvZpAPl8tJykknwd//7n25w4bBmDE+/0KOU+CKNEVVlX/k3X13ePZZH4nw2GPQtWvaLWt/Dj7Yq9sjjoALLoBDD4W33067VS2iwBVpjKoqH3mwzz4+7GufffzNf8EFmpCmLZWVwVNPwV13wfTpPv/Cf/0XfPBB2i1rFgWuyLbMmuUTzwwZAt/8pq976CG/7Pcee6Tbtnxh5leL+PBDv0TRk0/CnnvCUUf572LNmrRb2GgW8mTSCNm2oUOHhsmTJ6fdjPStXOnjaV94wbsM3njD1w8fDt//vvfbal6EdC1eDGPH+jJ3LnTsCEcfDUceCYcd5mFcVJRIU8xsSghhaKOfr8AVyKPA3bDBpwNcuBA+/dTfsLNmefU0dSpMm+bPKyqCAw+E//xPOPlkP5lBsktVlZ/ZN26cT4E5a5av79jRu3z22AN22cWnxRw40K+SXFbmw/da6Z+mAleaZbPAffxxGD++7Xe6tb+9zPrat1tbqqr86HVlpX9dUeHLhg2+rFvny6pVvtR3BlNRkV9rbI894IAD4AtfgC9+UQfCcs2sWfDKK/D66345+g8+qH9uBjP/3XbtCp07Q6dOUFrqVxYuKfEzBIuKfCks9PkeDjvMP+FssammBW4ydbfklhkz/ON0ErZ2wCmzvvZtfUvmDZF5cxQX+9KhA3Tp4tVOx441b7CePf0ssO2394pn0CDo10/dBO3Bjjv6kpkWE7x/d+5cv37awoV+IsXSpbBihf8DXrPGz2hbv97/Qa9a5f+wKyt9qa72pZUmj1eFK0AedSmItKKmVrgapSAikhAFrohIQtSlIACY2WJgTq1VvYFsnjkkm9uXzW2D7G5fNrcNtmzf4BBCWWO/WYEr9TKzyU3pm0paNrcvm9sG2d2+bG4btLx96lIQEUmIAldEJCEKXNmasWk3oAHZ3L5sbhtkd/uyuW3QwvapD1dEJCGqcEVEEqLAFRFJiAJXtmBmx5rZh2Y2w8wuTbktA83sBTN738ymmtmFcX0vM3vezKbH254ptrHQzN40s/Hx/o5mNim+fg+aWUmKbethZuPM7N9m9oGZHZxlr93o+Ht9z8weMLPSNF8/M7vTzBaZ2Xu11tX7epm7KbbzHTM7oKHtK3BlM2ZWCPwBOA7YE/iWme2ZYpMqgYtCCHsCw4DzYnsuBSaEEHYBJsT7abkQqH0Jgl8D14cQdgaWAWel0ip3I/BsCGF3YF+8nVnx2plZf+ACYGgIYS+gEDiVdF+/u4Bj66zb2ut1HLBLXEYBtzS49RCCFi2bFuBg4Lla9y8DLku7XbXa8wRwFPAhsENctwPwYUrtGRDfhEcA4wHDz0Qqqu/1TLht3YFZxIPjtdZny2vXH5gH9MJnLhwPHJP26weUA+819HoBtwHfqu95W1tU4UpdmTdBxvy4LnVmVg7sD0wC+oQQMpOdfgr0SalZNwA/BjKXlN0OWB5CqIz303z9dgQWA3+KXR63m1lnsuS1CyEsAK4F5gKfACuAKWTP65extderye8VBa7kBDPrAjwC/E8IYWXtx4KXF4mPbzSzE4BFIYQpSe+7kYqAA4BbQgj7A2uo032Q1msHEPtCR+D/GPoBndny43xWaenrpcCVuhYAA2vdHxDXpcbMivGwvS+E8GhcvdDMdoiP7wAsSqFphwBfNbPZwF/wboUbgR5mlpncP83Xbz4wP4QwKd4fhwdwNrx2AEcCs0IIi0MIFcCj+GuaLa9fxtZerya/VxS4UtfrwC7xSHEJfhDjybQaY2YG3AF8EEL4Xa2HngRGxq9H4n27iQohXBZCGBBCKMdfp7+HEE4DXgBOSbNtsX2fAvPMbLe4ajjwPlnw2kVzgWFm1in+njPty4rXr5atvV5PAqfH0QrDgBW1uh7ql0ZnuZbsXoDjgWnAR8DlKbflUPwj3DvAW3E5Hu8rnQBMB/4G9Eq5nYcD4+PXQ4DXgBnAw0CHFNu1HzA5vn6PAz2z6bUDrgL+DbwH/BnokObrBzyA9ydX4J8Qztra64UfIP1DfJ+8i4+22Ob2W3Rqr5kdi3+EKgRuDyH8qtkbExFp55oduHG85jR8iM58/KPot0II77de80RE2o+WXLX3QGBGCGEmgJn9BT/iuNXA7d27dygvL2/BLqWlli/3C5MOHLj5+ilTpiwJceb6owq+rhmNROp4vvrhrVxiuvFaErj1jUE7qO6TzGwUfhYGgwYNQleGTdePfwxjxkDdX4OZzan/O0SktbT5KIUQwtgQwtAQwtCyskZf+kfaSEUFFBen3QqR/NSSwM268ZrSsI0boSS1qVRE8ltLAjerxmtK46xdC507p90KkfzU7D7cEEKlmf0AeA4fFnZnCGFqq7VM2sTKldClS9qtEMlPLTloRgjhaeDpVmqLJOCzz2C77dJuhUh+0qm9eWbBAujbN+1WiOQnBW4eqaqC2bNhp53SbkkKzJq2iLQBBW4eee89qKyEvfZKuyUi+alFfbiSW156yW8POyzddmQFK6hzd/OqNlTXOdkuVNe5r5PxpOlU4eaRe++FPffc8rReEUmGKtw8MXEivP66n9ablzIVaZ3+2U2VbabijfetsM73ZyreWOmGOve32I9IPVTh5oGNG+F734M+feD009NujUj+UoXbzoUAl18Ob70Fjz8O3bql3aKUbapAM5Wq1xxWkKlUvbS1wljiZirigjp9vtXx+VVVvp2quL14f1Plq4pXalGF246F4LODXXutV7gjRqTdIpH8pgq3nVqxAi64AO65B847D266Ke0WZZmtVbrECjVT2cap1awovlUKY41Sd5RDdfy+Cr+6d6iMt5n7qnwFVbjt0uOP+2iEe++Fn/3MD5QV6DctkjpVuO3Iq6/CNdfAX/8K++zjwfuFL6TdqiyXqTRD7IsNdUYxbPoifhUrXSuJkwoXF23+eNyeVcbtVVT4+g0b/P7GirhelW8+Ut2T46qq4LHH4JBD4ItfhH/+E66+2q/ooLAVyS6qcHPUBx/Agw/CfffBjBlQXu79tGecoekXWyRWmJsqz02rfX2mQgmZ8bqxwg2lPqt7KI6jGzJ9OJmKd6NXtAXrN/r6dev94fWZ21gBV9atfFXxticK3BwybZqH7EMP+bwIZn6a7i9/CSedtOnTrohkKb1Fs9jq1fCPf8Df/gbPPw/vvushe+ihfiDs5JNhhx3SbmU7lal0Y8WZOdOsOt5ancrTYkUbSvwtVdXJ+3irO3jFu6kijt9fsMEr2MLVXtkWrFrnz1u91m/rVr6Zird688pbcosCN4tUVMCkSR6wEyb46biVlX4NskMOgRtugFNOgf79026piDSHAjdFixZ5wGaWV1+FNWu8iv385+Gii+DIIz1sO3ZMu7V5rjoziiGOJli3+ZlmmduCTGUcK95MhVvRxW8rS2MlnOnqrfRfbNFaPwWwZJWPYiha5hVv4Yo1/vyVq/12na+vjqMdNLohtyhwE7Jhg59eO3Gih+vEiTBrlj9WWAh77w0jR8Lw4XD44dCrV6rNFZE2oMBtA6tXwzvveMC++abfvvOOTyID3iUwbJifbnvQQV7N6kq6OaJO325m1jDLVLhxfWEcZ2tVXf155r/gqg5e4W7s7H26FZ38NsQz2Ao3+luyeHUpAKXLvfItXex/PEVLvNItXL7Sv291rIA31OnrVcWblRS4LbRw4ebB+uabMH16zd97r16w//5w4YUergcdBAMGpNtmEUmHAreR1q2D99/3kQK1l08/rXlOebmH62mnwX77+dcDBugSWe1apm93Y3W86/ct9rFaPNOsw4Y4DjdWvNABgKoO/sexvruvXdc50wfs9ws3eGdvyYpOAJQu8T7fzgt7+v1PfFRD4ZIV/n0rvPKtzozzrcz09arizQYK3DqqqmDmzC2DdcYMyMzIV1rqcxUcc0xNsO67L/TokW7bRSS75XXgLlq0ZbBOnQprvWjAzK9wu/fecOqpfrv33rDzzn6gS2STun27mUo39uUWxIqzJI6rLVzr/50LN8bhJ8Hfimv9hDU29o59sf28Ql1nvv3Va3x87/LP/PmdFnpp3Pljr5w7L/D9FH/qFS/L/LY609dbkRlXrPG8aciLwF271oO0brguWlTznLIyD9NzzqkJ1s99TgezRKT1tLvAXbrUD1y98UbN7bRpNV1YpaV+mfCvfKUmWPfe2y8/I9JqMhVvhY8uqMqMYohDVQrW+MeoLiu9L7Z4lVeqRbHEXV3hb831hb6dXn29Uu2zg49SYGe/+WSlj2L4ZKFXuKULvK+3yzyvnLvO7RXXe98uS5Z7u1atAmqN51XFm4icDdwQ4OOPtwzXuXNrnjNwIBxwwObdATvtpO4AEUlHzgRudbV3C7z8ss8v8PLL8MknNY/vuiscfLBf3WD//X3p3Tu99opsJlaQ1Rv8yKvFvl6LfbulcQ6F4uVx9MFyn/Jt5Urvs122ztdXD/bhC0P7zgPgy2UfArBhR3/eWyt8zOE7C/oBsHyO94l1mVMGQPfZ3nfcca5XuIWLlwIQVnnlXB37mFXxto2sDdyKCpgypSZc//UvWLbMH+vf38/GGjbMK9h994WuXbe5ORGR1GVV4IbgE2jfc49PQbgydjvtuqtPP/ilL/l0hOXlGtsqOWoroxkyV4IoiJVmt6Wx0l0cb+P425WfeZ/shJ28r3bZjr7+hLJ3ADh6wLsALN3BK+QJu+4JwN/m7QbAvFle4XaZ5dvpMdP7gDvFirdooVe81StjH68q3lbVYOCa2UDgHqAPEICxIYQbzawX8CBQDswGvhFCWNacRsycCX/+swftzJk+MuDkk+HEE30qwr59m7NVEZHsYqGBM1DMbAdghxDCG2bWFZgCfA34LrA0hPArM7sU6BlCuGRb2xo6dGiYPHnyZuuuuw5+9COvWIcPh9NP92pWw7GSZWZTQghDAY4q+LpOS0pL/OhmJT5aoaBLfCNs55Xu+kFeoa4Y4o+v3Mkf7rCLfxw8rvx9AE7tOQmA3Yq9z/jjWEk/tnI/v523LwCLZmwHQNeZfiS5x0deeXea49uzWPFuOoMtj0c1PF/9cIs/VzdY4YYQPgE+iV+vMrMPgP7ACODw+LS7gReBbQZuXWPHetiefDJcf72PKhARaa+a1IdrZuXA/sAkoE8MY4BP8S6HRps1C849F4YOhfvv90m2RfJepo83zv5VlTlTLY7bLV3mlWbpJ94H23WB98Eun+/jeB/d6SAAXtrFB+qePOgtAL7R7U0ALtlu+mb3Hxq4PwDjhvjt/PK43Y9iH+9HfjRmlkPGAAAPNElEQVQ6U/EWLsqMatA43uZo9FV7zawL8AjwPyGElbUfC94vUe/HUDMbZWaTzWzy4sWLN60fMAC+/GUfP/v0081rvIhILmmwDxfAzIqB8cBzIYTfxXUfAoeHED6J/bwvhhB229Z26vbhrloFRx0Fr7/uE8GMHAkjRvjZYJIs9eHmiALvay0o8XG31t0r3NDHK9K1g2PFu5N/eF01xCvP7Xf+DID/HPi233bzyrdfPAvowwqvvf6yzCvkZ2b76IYN03173T7y3Xef6WfKlc71M9b4zI+Tb5qrITPpczucnaw1+nAbrHDNzIA7gA8yYRs9CYyMX48Enmjqzrt2hWefhUsu8bkNTj3VRySMGuXjbjOzc4mItAeNGaVwKPAP4F0gE4E/wftxHwIGAXPwYWFLt7Wt+kYpZFRVwYsvwt13wyOP+IQzPXv6sLDDDvMxuAccAMXFTfnxpLFU4eaoTMVb6vPrFnTzPtfqTMU7yO8vH+IV7+od/S3cfYhXpsMHTPPbbj66oVehjwN+f4NfqXT84n0AeGPWIAA6fOQfP7vN9D+RbrP8GmslC3x7IVPx1p6Pt51Uu0mNUvgnsLUdDW9pAzIKC31Y2PDhcPPN8MQT8MILfqbZX//qz+nUyU/fzZwAceCBGj4mIrmjUX24rWVbFe62fPqpn4GWmUfh7bf9n6YZ7Labz5twwAE1cyjoAoxNpwq3nahT8VpXP+MslPmbYt1gr3hXDPaPiqsHx2uyDfY+2H36fwzAft3nA9ChwEchTFvjZx9N/tTHbq6Y46Mius7y/XWb7X3FnWfHa64tipXu8hXt5npriVS42aBvXzjlFF8Ali+HV17xg21vvulh/MADNc8fPLgmgA84wK/K0K+fTgcWkXTlRODW1aMHHH+8LxlLlnj4ZpY33oDHH6/5Z9qzZ80Ujfvs47d77aVJb6SdycxKlrlsSZwLoSCeKdYpzpvbaaafsba+v49CWDXI++amDtwVgCn9dwSgRx8fb7tDt5Wb3RbEvuClHb3S3djNK+b1PeKohjne19thQScs9utmrjCcz2N3czJw69O7tw8xO+qomnWrVnn3w9tv+yiId97x+RrimG3AJ8KpPRH53nv7ZDk6OCftQiaA18dwy1zcMs4MVbrIg7Z0tgdnj75+u6a/B+aaft4VMb2Pn1pcuV28DHxn305hJ7+/fgf/+Fhd7F0MlR29S6Nrt150WuAT7BQu8qFklrnsTx5e6LLdBG59unb1UQ6HHlqzLgSYM2fLy+088wzECZwoKYHdd98yiHUFXhFpiXYduPUx86q2vNxnI8vYsAE+/NCr4EwIv/QS3HdfzXN69PBuiLpdE926Jf1TiDRT5rLuG/y2Kp6oYCv9YFfxIj9Bouec2DWwvVe86/v4dJBr+sRL//T2j4Abu8eqtEPcfEyU9ZnJ/wuKqCrxKrpTJ3+wODOEbWm8tPuazU+ayPWDa9uSd4G7NR06eIDus8/m65ctg/fe27wavu++mrl6wa/im7lceua2b19VwyKyOQVuA3r29HG/X/pSzboQYN48r4bfesuXN96AceNqnrP99jUBvN9+8PnPezArhCWr1L3YZexPtXjQzZb6Aa/OC3x4WafuftCtoszvb+jts06t7+F9txVd/A+8Kk5GVV0IG7pn/ui9si0t9BNcSzp4/BQs9SeH1fEyP5m+3XZ4SXcFbjOYwaBBvpxwQs36FSs8hN9800P4zTfhd7/zywWBjw8+6CBfhg3zEzd69kznZxCR5ClwW1H37ltWwxs3wvvvw+TJMGkSTJzo80dkuqd23dXD96CD4Itf9C6NgkbP4SbSyrZyCaBM1Vmw3Ptdixf5yIOSbl7pdu7u/bSVPX39xq7ex1vZqWBTv26IhW5FV6+GLXi/cHH8gy8ojhVvkVe6Ya2fNtye+nYVuG2spKSmW+Hss33dypUewBMnegg/+6wPVwMf3vblL8ORR/ppzkOGqBtCpL1Q4KagWzc44ghfoGao2ssvw4QJvjz8sD9WXl4zx8Sxx6oLQhKWqSbDVsbzxpMZbKn3zxYvjJVvZ69eq7uUUtXFH6vu4JVtKNi8gqjqFAe9x4q3IFNhxKkjiZVu5iSOTZVuDvbtKnCzQO2haqef7n/jH35YE76PPAJ33OEnYxx9NHzjGz5vcPfuabdcRJpCgZuFzPzEi913h/PO86krX3/dg/ehh+Cpp7yr4thjPXxPPlmTtkvCMuN5M7eZ0Q3rvBq1eCqxlXagOPPH2dFvQ2kclVAS4ydT0Wb6jzvEydWrN/+jtvi86sxkOPFgNKF6s+/PZjo8kwMKC/3A2m9/C7Nnw6uvehBPmQLf/rZP1vPzn8Nnn6XdUhHZlpyYnlHqV13tcwZfd52fmtyxI5x5Jlx8sYdwU2h6RmlVmarVCrDYF2txFILFywPRIU4hmZm4pKhw8+/NZFMcjxsyk95kKtzM6IU6IyraqtJN5BI7kr0KCvxg2tNP11yiaOxY2HNPuOEG74oQkeyhwG0n9toL7rwTpk+Hww+H0aN9XO/UqWm3TPJSCL5UVxEqNhIqNlK9bh3V69ZRtXI1VStXU710eVyWUb10GWHpcl9WrvJlzTpfNlZ4dRuqfSks9KW4GIqLsaIiXwoLvZouiItZ1o2pVOC2M4MHw/jxcP/9MHOmh+6//pV2q0QENEqhXTKDb33Lp6U88kgfSvbkk979IJKaOmN6txjhkBl3awXxfqwHM+vrnoKZuax3webPMzLbzzw/e0YxqMJtxwYO9JMphgzxyxMtWJB2i0TymwK3nevTBx57zA/snnNOVvyTF9lc7O8NlZVxqSBUVlC9YYMv69ZTvW49Yd06X9Zv8CXz/KqqmhEK4BWyFWAFhhXYpvvZ0KerwM0DO+8MV1/tQ8deeCHt1ojkLwVunjj3XJ8e8uab026JSAMyIxxqjXSgumpTJZupgENF5WYLVVW+xNEMoToQquv5SJdipavAzROlpTByJDzxBMQrmohIwhS4eeToo/1CmRMnpt0SkWbYSuW7RQWc6dPNjNvNLFlAgZtHDj7Yb197Ld12iOQrjcPNI927+6iFjz5KuyUibSAHhuCows0z5eUwd27arRDJTwrcPNO7t6ZxFEmLAjfP9OwJy5al3QqR/NTowDWzQjN708zGx/s7mtkkM5thZg+aWUnbNVNaS6dOECflF5GENaXCvRD4oNb9XwPXhxB2BpYBZ7Vmw6RtlJYqcEXS0qjANbMBwFeA2+N9A44AxsWn3A18rS0aKK2ruBgqKhp+noi0vsZWuDcAP2bTPGdsBywPIVTG+/OB/vV9o5mNMrPJZjZ58eLFLWqstFxRkZ/8ICLJazBwzewEYFEIYUpzdhBCGBtCGBpCGFpWVtacTUgrKiiomUZURJLVmBMfDgG+ambHA6VAN+BGoIeZFcUqdwCg2VZzgAJXJD0NVrghhMtCCANCCOXAqcDfQwinAS8Ap8SnjQSeaLNWSqsxy4kTckTapZaMw70E+KGZzcD7dO9onSZJW1LgiqSnSXMphBBeBF6MX88EDmz9JklbyrKLmIrkFZ1pJiKSEAVunlGFK5IeBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIaFbhm1sPMxpnZv83sAzM72Mx6mdnzZjY93vZs68aKiOSyxla4NwLPhhB2B/YFPgAuBSaEEHYBJsT7IiKyFQ0Grpl1Bw4D7gAIIWwMISwHRgB3x6fdDXytrRopItIeNKbC3RFYDPzJzN40s9vNrDPQJ4TwSXzOp0Cf+r7ZzEaZ2WQzm7x48eLWabWISA5qTOAWAQcAt4QQ9gfWUKf7IIQQgFDfN4cQxoYQhoYQhpaVlbW0vSIiOasxgTsfmB9CmBTvj8MDeKGZ7QAQbxe1TRNFRNqHBgM3hPApMM/MdourhgPvA08CI+O6kcATbdJCEZF2oqiRzzsfuM/MSoCZwBl4WD9kZmcBc4BvtE0TRUTah0YFbgjhLWBoPQ8Nb93miIi0XzrTTEQkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGENCpwzWy0mU01s/fM7AEzKzWzHc1skpnNMLMHzaykrRsrIpLLGgxcM+sPXAAMDSHsBRQCpwK/Bq4PIewMLAPOasuGiojkusZ2KRQBHc2sCOgEfAIcAYyLj98NfK31myci0n40GLghhAXAtcBcPGhXAFOA5SGEyvi0+UD/+r7fzEaZ2WQzm7x48eLWabWISA5qTJdCT2AEsCPQD+gMHNvYHYQQxoYQhoYQhpaVlTW7oSIiua4xXQpHArNCCItDCBXAo8AhQI/YxQAwAFjQRm0UEWkXGhO4c4FhZtbJzAwYDrwPvACcEp8zEniibZooItI+NKYPdxJ+cOwN4N34PWOBS4AfmtkMYDvgjjZsp4hIzitq+CkQQvhf4H/rrJ4JHNjqLRIRaad0ppmISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBa6ISEIUuCIiCVHgiogkRIErIpIQBW6eGTYMLrww7VaI5CcLISS3M7PFwJzEdihNMTiEUJZ2I0Tas0QDV0Qkn6lLQUQkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGEKHBFRBKiwBURSYgCV0QkIQpcEZGE/H+ovX+5d2RUpwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2bf660a438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Example with Frobenius norm + entropic regularization with gcg\n",
    "\n",
    "\n",
    "def f(G):\n",
    "    return 0.5 * np.sum(G**2)\n",
    "\n",
    "\n",
    "def df(G):\n",
    "    return G\n",
    "\n",
    "\n",
    "reg1 = 1e-3\n",
    "reg2 = 1e-1\n",
    "\n",
    "Gel2 = ot.optim.gcg(a, b, M, reg1, reg2, f, df, verbose=True)\n",
    "\n",
    "pl.figure(5, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Gel2, 'OT entropic + matrix Frob. reg')\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}