summaryrefslogtreecommitdiff
path: root/notebooks/plot_stochastic.ipynb
diff options
context:
space:
mode:
authorRémi Flamary <remi.flamary@gmail.com>2018-08-29 14:10:40 +0200
committerRémi Flamary <remi.flamary@gmail.com>2018-08-29 14:10:40 +0200
commit1ec01f082a34c1086ef2df94b78a5d640a66db07 (patch)
treedc905fe366cb4b04187d3b4c760d8a03575ee97b /notebooks/plot_stochastic.ipynb
parentf12153c0c1be6f6377ace0050201409ec1b7e829 (diff)
add notebooks from new examples
Diffstat (limited to 'notebooks/plot_stochastic.ipynb')
-rw-r--r--notebooks/plot_stochastic.ipynb610
1 files changed, 610 insertions, 0 deletions
diff --git a/notebooks/plot_stochastic.ipynb b/notebooks/plot_stochastic.ipynb
new file mode 100644
index 0000000..e784e11
--- /dev/null
+++ b/notebooks/plot_stochastic.ipynb
@@ -0,0 +1,610 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "# Stochastic examples\n",
+ "\n",
+ "\n",
+ "This example is designed to show how to use the stochatic optimization\n",
+ "algorithms for descrete and semicontinous measures from the POT library.\n",
+ "\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Author: Kilian Fatras <kilian.fatras@gmail.com>\n",
+ "#\n",
+ "# License: MIT License\n",
+ "\n",
+ "import matplotlib.pylab as pl\n",
+ "import numpy as np\n",
+ "import ot\n",
+ "import ot.plot"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "COMPUTE TRANSPORTATION MATRIX FOR SEMI-DUAL PROBLEM\n",
+ "############################################################################\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "------------SEMI-DUAL PROBLEM------------\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"------------SEMI-DUAL PROBLEM------------\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "DISCRETE CASE\n",
+ "Sample two discrete measures for the discrete case\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Define 2 discrete measures a and b, the points where are defined the source\n",
+ "and the target measures and finally the cost matrix c.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "n_source = 7\n",
+ "n_target = 4\n",
+ "reg = 1\n",
+ "numItermax = 1000\n",
+ "\n",
+ "a = ot.utils.unif(n_source)\n",
+ "b = ot.utils.unif(n_target)\n",
+ "\n",
+ "rng = np.random.RandomState(0)\n",
+ "X_source = rng.randn(n_source, 2)\n",
+ "Y_target = rng.randn(n_target, 2)\n",
+ "M = ot.dist(X_source, Y_target)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Call the \"SAG\" method to find the transportation matrix in the discrete case\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Define the method \"SAG\", call ot.solve_semi_dual_entropic and plot the\n",
+ "results.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[2.55553509e-02 9.96395660e-02 1.76579142e-02 4.31178196e-06]\n",
+ " [1.21640234e-01 1.25357448e-02 1.30225078e-03 7.37891338e-03]\n",
+ " [3.56123975e-03 7.61451746e-02 6.31505947e-02 1.33831456e-07]\n",
+ " [2.61515202e-02 3.34246014e-02 8.28734709e-02 4.07550428e-04]\n",
+ " [9.85500870e-03 7.52288517e-04 1.08262628e-02 1.21423583e-01]\n",
+ " [2.16904253e-02 9.03825797e-04 1.87178503e-03 1.18391107e-01]\n",
+ " [4.15462212e-02 2.65987989e-02 7.23177216e-02 2.39440107e-03]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "method = \"SAG\"\n",
+ "sag_pi = ot.stochastic.solve_semi_dual_entropic(a, b, M, reg, method,\n",
+ " numItermax)\n",
+ "print(sag_pi)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "SEMICONTINOUS CASE\n",
+ "Sample one general measure a, one discrete measures b for the semicontinous\n",
+ "case\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Define one general measure a, one discrete measures b, the points where\n",
+ "are defined the source and the target measures and finally the cost matrix c.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "n_source = 7\n",
+ "n_target = 4\n",
+ "reg = 1\n",
+ "numItermax = 1000\n",
+ "log = True\n",
+ "\n",
+ "a = ot.utils.unif(n_source)\n",
+ "b = ot.utils.unif(n_target)\n",
+ "\n",
+ "rng = np.random.RandomState(0)\n",
+ "X_source = rng.randn(n_source, 2)\n",
+ "Y_target = rng.randn(n_target, 2)\n",
+ "M = ot.dist(X_source, Y_target)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Call the \"ASGD\" method to find the transportation matrix in the semicontinous\n",
+ "case\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Define the method \"ASGD\", call ot.solve_semi_dual_entropic and plot the\n",
+ "results.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[3.75309361 7.63288278 3.76418767 2.53747778 1.70389504 3.53981297\n",
+ " 2.67663944] [-2.49164966 -2.25281897 -0.77666675 5.52113539]\n",
+ "[[2.19699465e-02 1.03185982e-01 1.76983379e-02 2.87611188e-06]\n",
+ " [1.20688044e-01 1.49823131e-02 1.50635578e-03 5.68043045e-03]\n",
+ " [3.01194583e-03 7.75764779e-02 6.22686313e-02 8.78225379e-08]\n",
+ " [2.28707628e-02 3.52120795e-02 8.44977549e-02 2.76545693e-04]\n",
+ " [1.19721129e-02 1.10087991e-03 1.53333937e-02 1.14450756e-01]\n",
+ " [2.65247890e-02 1.33140544e-03 2.66861405e-03 1.12332334e-01]\n",
+ " [3.71512413e-02 2.86513804e-02 7.53932500e-02 1.66127118e-03]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "method = \"ASGD\"\n",
+ "asgd_pi, log_asgd = ot.stochastic.solve_semi_dual_entropic(a, b, M, reg, method,\n",
+ " numItermax, log=log)\n",
+ "print(log_asgd['alpha'], log_asgd['beta'])\n",
+ "print(asgd_pi)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Compare the results with the Sinkhorn algorithm\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Call the Sinkhorn algorithm from POT\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[2.55535622e-02 9.96413843e-02 1.76578860e-02 4.31043335e-06]\n",
+ " [1.21640742e-01 1.25369034e-02 1.30234529e-03 7.37715259e-03]\n",
+ " [3.56096458e-03 7.61460101e-02 6.31500344e-02 1.33788624e-07]\n",
+ " [2.61499607e-02 3.34255577e-02 8.28741973e-02 4.07427179e-04]\n",
+ " [9.85698720e-03 7.52505948e-04 1.08291770e-02 1.21418473e-01]\n",
+ " [2.16947591e-02 9.04086158e-04 1.87228707e-03 1.18386011e-01]\n",
+ " [4.15442692e-02 2.65998963e-02 7.23192701e-02 2.39370724e-03]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "sinkhorn_pi = ot.sinkhorn(a, b, M, reg)\n",
+ "print(sinkhorn_pi)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "PLOT TRANSPORTATION MATRIX\n",
+ "#############################################################################\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Plot SAG results\n",
+ "----------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAExZJREFUeJzt3X+wpQV93/H3hwUVAhHj3ihhxbVqd4pawdxgGqwa/IXE/JqYalJ/RdutrVhpTa0mmY7WaZImUyWd2qRbNcaIEo06k+aHhQYYytQfvasbhh8ygw66rCAXCAoUsCzf/vGcbe/c7u49u3vO+e6e837NnNl773nOeb7nwH3f5zznOeekqpAkzd5x3QNI0qIywJLUxABLUhMDLElNDLAkNTHAktTEAGvuJPn7SS47hOXfkOSaCa37liQvnsR1HYuSnJHkviSbumc5FhhgzZ2quqSqXto9x+FI8ookX0pyf5K7klySZMvovF8Zxe2+JA8m2bvm++tnMNuGf1yq6ptVdXJV7T2M639GksuS3J3kniQ7k1ywbpmnJHkkye/u5/JJcmGSa5P8ryS3J7kqyasPdZZZMcDSUSLJK4GPAxcDm4FnAA8B1yR5XFX9+ihuJwNvBj6/7/uqekbf5IMkxx/hVfwX4HLgicAPAv8U+O66ZV4H/DXwqiSPXnfevwcuAt4OPB44Hfg14PwjnGt6qsqTp5mdgH8J7AHuBW4CXjT6+XHAO4GvAXcBnwR+YHTeVqCAXwJ2M/wCvhn4EeBa4B7gP6xZxxuAaw4yw+OBP2H45f4S8N59y69Z1/Frlr8K+Aejr58KXDGa8U7gEuDUNcveArz4MO6XAN8A3rHu58cB1wH/et3PD3obD/N+O+BtA/4QeAR4ALgPeMea638T8E3g6rX3H/ADwK3AT46u42TgZuB1+5l18+hyp25wH30N+MfAt4FXrjnvbwJ7geXu/8cP5eQWsGYmyTbgQuBHquoU4GUMwQJ4K/AzwAuAH2KIxQfWXcVzgacDr2LYSvxV4MUMW4p/L8kLxhzlA8CDwGnAG0ensW8G8BujGf8W8CTg3WNdMPnFJNce4OxtwBnAp9b+sKoeAT4NvOQQZlxv3PvtgLetql7LENmfrGGL+7fWXP8LRsu/bN3sdzPct/85yQ8C7wd2VdVH9zPjXQxx/liSn0nyhP0s8zxgC3Apwx/o16857zxgd1WtbHhvHEUMsGZpL/Bo4MwkJ1TVLVX1tdF5bwZ+tapuraqHGH7xX7nuYe17q+rBqroMuB/4RFXdUVV7gP8OnL3RAKMnh34O+FdVdX9VXQf8wbg3oKpurqrLq+qhqloF3scQoHEu+/Gq+tsHOHvz6N/b9nPebWvOPxxj3W9HcNvePbovH1h/xmidnwL+ErgA+Ef7u4IaNmN/nOEP8r8DbktydZKnr1ns9cBfVNVfM+yqOX8Udhjun9vXXmeSW0f7kh9M8uQxbsfMGWDNTFXdzLCP7t3AHUkuTfJDo7OfDHx29AtzD3AjQ7DXbgl9e83XD+zn+5PXr3PdE1e/BywxPDzevWaxb4x7G5I8YTT3niTfBT7GkcVxnztH/562n/NOW3P+4RjrfjuC27Z7g/N3AM8EPlJVdx1oodEf3wur6qkM/z/cD3x0NNuJwM8z7Bahqj7PsEX+i6OL38W6+66qtozmfzTD1v1RxwBrpkZbgc9j+AUr4N+OztoNvLyqTl1zesxoK+1I1vd/n7iqqjcDq8DDDA+v9zljzdf3j/49ac3Pnrjm618fzf2sqvp+4DVM5pf7Job9pT+/9odJjmPYYv/LCaxjIxvdtgO9deIB31Jx9IhjB0NI/0mSp40zSFXtZthV9MzRj34W+H7gP46Obrid4Um2fbshrgC2JFke5/qPFgZYM5NkW5LzRs9eP8iw9fXI6OzfA/7NvoeKSZaS/PSkZ6jh8KjPAO9OclKSM1mzL3H00HsP8Jokm5K8keHJqX1OYXgS6jtJTgf+xYTmKuCXgV8b7St+TJInAh9kCM/7J7GeDWx0274N/I1DvM5fYQj0G4HfBj66v2OEkzwuyXuSPC3JcUk2jy7zhdEirwc+DDwLOGt0Ohd4dpJnVdVNwH8CLk3ykiQnjtbzY4c470wZYM3So4HfZHg4fTvDoUbvGp33OwxHJlyW5F6GX7znTmmOCxkedt8OfAT4/XXn/0OG+NzF8ETV/1hz3nuA5wDfAf6MIeZjGb1A5IDH61bVHwGvBf7ZaN03ACcC5x7sofsEbXTbfoPhD8Q9SX55oytL8sPAP2c46mEvw6OdYjjaZb3vMRxB8d8Yjk65juEQvDeM/hi8CLi4qm5fc9oJfI7/9wf0LQyHor0PuJvhEcV7GZ58/OZY98CMZXQIhyRpxtwClqQmBliSmhhgSWpigCWpyZG+eYaOcZs3b66tW7d2jyHNlZ07d95ZVUsbLWeAF9zWrVtZWTmmXj4vHfWSjPXqSndBSFITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNju8eQM1uugle+MLuKbTIzjoLLr64e4oWbgFLUhO3gBfdtm1w1VXdU0gLyS1gSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJqmq7hnUKMm9wE3dc0zYZuDO7iGmwNt17NhWVadstNDxs5hER7Wbqmq5e4hJSrIyb7cJvF3HkiQr4yznLghJamKAJamJAdaO7gGmYB5vE3i7jiVj3SafhJOkJm4BS1ITAyxJTQzwgkpyfpKbktyc5J3d80xCkg8nuSPJdd2zTEqSJyW5MskNSa5P8rbumSYhyWOSfCnJX41u13u6Z5qUJJuSfCXJn260rAFeQEk2AR8AXg6cCfxCkjN7p5qIjwDndw8xYQ8Db6+qM4EfBd4yJ/+tHgLOq6pnA2cB5yf50eaZJuVtwI3jLGiAF9M5wM1V9fWq+h5wKfDTzTMdsaq6Gri7e45JqqrbqurLo6/vZfjFPr13qiNXg/tG354wOh3zRwQk2QL8BPDBcZY3wIvpdGD3mu9vZQ5+qeddkq3A2cAXeyeZjNFD9V3AHcDlVTUPt+ti4B3AI+MsbIClY0CSk4FPAxdV1Xe755mEqtpbVWcBW4Bzkjyze6YjkeQVwB1VtXPcyxjgxbQHeNKa77eMfqajUJITGOJ7SVV9pnueSauqe4ArOfb3358L/FSSWxh2652X5GMHu4ABXkz/E3h6kqckeRTwauBPmmfSfiQJ8CHgxqp6X/c8k5JkKcmpo69PBF4CfLV3qiNTVe+qqi1VtZXhd+qKqnrNwS5jgBdQVT0MXAj8V4YndT5ZVdf3TnXkknwC+DywLcmtSd7UPdMEnAu8lmFratfodEH3UBNwGnBlkmsZNggur6oND9uaN74UWZKauAUsSU2m8obsmzdvrq1bt07jqjVhO3fuvLOqlrrnOFIvfOlvHtZDuZe9/+pJj3JQV77unJmuD6C+Mtu9S5c/8qnMdIXHsKkEeOvWraysjPWG8GqW5BvdM0iLyl0QktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUZKwAz+MHOEpStw0DPMcf4ChJrcbZAp7LD3A8FBddNJwkaZLGeTOe/X2A43PXL5RkO7Ad4IwzzpjIcEeLXbu6J5A0jyb2JFxV7aiq5apaXlo65t/dUJKmbpwA+wGOkjQF4wTYD3CUpCnYcB9wVT2cZN8HOG4CPjwPH+AoSd3G+kSMqvpz4M+nPIskLRRfCSdJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSk7FeiCEd7a746IcO63IXPP9nJzzJwdXXvzrT9QFs8s2xjlpuAUtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNNgxwkg8nuSPJdbMYSJIWxThbwB8Bzp/yHJK0cDYMcFVdDdw9g1kkaaG4D1iSmkwswEm2J1lJsrK6ujqpq5WkuTWxAFfVjqparqrlJd9/VJI25C4ISWoyzmFonwA+D2xLcmuSN01/LEmafxt+JFFV/cIsBpGkReMuCElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJarLhK+GkY8HLn/Zjh3W5b/7hSROe5OAe+NbyTNcH8PS3fnHm69R43AKWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmozzqchPSnJlkhuSXJ/kbbMYTJLm3TjvBfEw8Paq+nKSU4CdSS6vqhumPJskzbUNt4Cr6raq+vLo63uBG4HTpz2YJM27Q9oHnGQrcDbw/729UpLtSVaSrKyurk5mOkmaY2MHOMnJwKeBi6rqu+vPr6odVbVcVctLS0uTnFGS5tJYAU5yAkN8L6mqz0x3JElaDOMcBRHgQ8CNVfW+6Y8kSYthnC3gc4HXAucl2TU6XTDluSRp7m14GFpVXQNkBrNI0kLxlXCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNRnn/YClo96Df/fMw7rcYz8521+Bx7/x2zNdn45ubgFLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTcb5VOTHJPlSkr9Kcn2S98xiMEmad+O8EP4h4Lyqui/JCcA1Sf6iqr4w5dkkaa6N86nIBdw3+vaE0ammOZQkLYKx9gEn2ZRkF3AHcHlVfXE/y2xPspJkZXV1ddJzStLcGSvAVbW3qs4CtgDnJHnmfpbZUVXLVbW8tLQ06Tklae4c0lEQVXUPcCVw/nTGkaTFMc5REEtJTh19fSLwEuCr0x5MkubdOEdBnAb8QZJNDMH+ZFX96XTHkqT5N85RENcCZ89gFklaKL4STpKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQm47wSTjrqnXT9bYd1uUft+daEJzm447+wZabrA/izb+2a+To1HreAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCZjBzjJpiRfSeIHckrSBBzKFvDbgBunNYgkLZqxApxkC/ATwAenO44kLY5xt4AvBt4BPHKgBZJsT7KSZGV1dXUiw0nSPNswwEleAdxRVTsPtlxV7aiq5apaXlpamtiAkjSvxtkCPhf4qSS3AJcC5yX52FSnkqQFsGGAq+pdVbWlqrYCrwauqKrXTH0ySZpzHgcsSU0O6SOJquoq4KqpTCJJC8YtYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaHNILMaSj1f9+8uG9AVT2fGvCkxzc3j23zXR9AN955IGZru9xM13bsc0tYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJanJWC9FHn0k/b3AXuDhqlqe5lCStAgO5b0gfryq7pzaJJK0YNwFIUlNxg1wAZcl2Zlk+/4WSLI9yUqSldXV1clNKElzatwAP6+qngO8HHhLkuevX6CqdlTVclUtLy0d3lsDStIiGSvAVbVn9O8dwGeBc6Y5lCQtgg0DnOT7kpyy72vgpcB10x5MkubdOEdBPAH4bJJ9y3+8qj431akkaQFsGOCq+jrw7BnMIkkLxcPQJKmJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpyaG8Ibt01LrzWSce1uUee/IPT3iSg9v9+odnuj6AVz1100zXd9kDM13dMc0tYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJajJWgJOcmuSPk3w1yY1J/s60B5OkeTfuS5F/B/hcVb0yyaOAk6Y4kyQthA0DnOSxwPOBNwBU1feA7013LEmaf+PsgngKsAr8fpKvJPlgku+b8lySNPfGCfDxwHOA362qs4H7gXeuXyjJ9iQrSVZWV1cnPGavs84aTpI0SePsA74VuLWqvjj6/o/ZT4CragewA2B5ebkmNuFR4OKLuyeQNI823AKuqtuB3Um2jX70IuCGqU4lSQtg3KMg3gpcMjoC4uvAL01vJElaDGMFuKp2ActTnkWSFoqvhJOkJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCapmvz75iRZBb4x8SvWNDy5qpa6h5AW0VQCLEnamLsgJKmJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCb/B6HXs8MRx/3SAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 360x360 with 3 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "pl.figure(4, figsize=(5, 5))\n",
+ "ot.plot.plot1D_mat(a, b, sag_pi, 'semi-dual : OT matrix SAG')\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Plot ASGD results\n",
+ "-----------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE3lJREFUeJzt3X+wpQdd3/H3h81PSCTgXiFmExYKbom0JuQSsLGggWgSIjojVVCI/Gi3TAmTtFga1HagjtofMzR2ZNRtxEgJpInA6Fhqk5FkMFOE3oU1zQ/WiUzCbiTkBoz5BWGy+faP56xzve7uPbt7zv3unvN+zZzZe+/znPN8z9nsO899znPOSVUhSVp/T+seQJLmlQGWpCYGWJKaGGBJamKAJamJAZakJgZYR7UkP5PkxoNY/y1Jbp3Qtu9J8ppJ3NbRKMkZSR5NsqF7lqOVAdZRraquraof7p7jUCS5JMnnkzyW5OtJrk2yabTs50dxezTJt5LsWfH9Hesw25r/c6mqr1TVSVW15zC2c02SJ5OcuurnpyT5UJL7kzyS5M+TXLlieZJcluS2JI+P1rslyRtWrHPL6LF7JMnDSbYnuTLJ8Yc676QZYKlBktcDHwWuAjYC3ws8Adya5FlV9SujuJ0EvAP47N7vq+p7+yYfJDlmArfxDOAngL8G3rRq8X8BTgJeDDwTeB1w94rl/xW4Ang38J3AacAvAheuup3Lqupk4NTRum8APpUkhzv/RFSVFy9TuQD/BrgPeATYCbx69POnAVcCfwF8HbgeePZo2WaggLcCu4C/YgjQy4DbgIeAX1+xjbcAtx5ghu8E/gB4GPg88Et711+xrWNWrH8L8E9HX/894NOjGR8ErgVOWbHuPcBrDuFxCXAv8J5VP38acDvw71f9/ID38RAft/3eN+C/A08B3wQeBd6z4vbfDnwF+MzKxw94NrAb+NHRbZzEEMxLDzDzpaNZLwduX7XsduDH93O97wH2AItrPCZ/83e54mdnAI8Dl3T/+6gq94A1HUm2AJcBL6thD+RHGIIF8C7gx4FXAd/NEIsPrrqJlwMvAn6KYS/xF4DXMOwp/mSSV405ygeBbzHsAb1tdBn7bgC/OprxxcDpwPvGumLy00lu28/iLQwhuGHlD6vqKeDjwAUHMeNq4z5u+71vVfVmhsj+aA173P9pxe2/arT+j6ya/RsMj+1/S/JdDHuwO6rqwweY9WeBjwHXAX8/yTkrlv0p8MtJ3prkRauudz6wq6qW1ngs/o6q+gqwBPzjg73uNBhgTcse4HjgzCTHVtU9VfUXo2XvAH6hqnZX1RMM//Bfv+rX2l+qqm9V1Y3AY8DHquqBqroP+BPg7LUGGD059BPAv6uqx6rqduB3x70DVXV3Vd1UVU9U1TLwAYYAjXPdj1bVP9zP4o2jP7+6j2VfXbH8UIz1uB3GfXvf6LH85uoFo23eAPwxcDHwz/d3I0nOAH4I+GhVfW10nUtXrPIuhr3yy4A7k9yd5KLRso3A/atub3eSh0bHfJ+3xn34S4Y99nYGWFNRVXczHKN7H/BAkuuSfPdo8fOAT47+wTwE3MUQ7OesuImvrfj6m/v4/qTV21z1xNVvAgsMvx7vWrHavePehyTPGc19X5KHgY9weHHc68HRn6fuY9mpK5YfirEet8O4b7vWWL4NeAlwTVV9/QDrvRm4q6p2jL6/FvjpJMcCVNU3azgOfg7DYaTrgRuSPJvhsMnfeuyqatNo/uMZ9u4P5DTgG2ussy4MsKZmtBf4AwzBLeA/jhbtAi6qqlNWXE4Y7aUdzvb+5omrqnoHsAw8yfDr9V5nrPj6sdGfT1/xs+eu+PpXRnP/g6r6DoYniibx5M1OhuOl/2TlD5M8jWGP/Y8nsI21rHXf9vc2ift9+8TRbxzbgA8D/yLJCw+w/UuBF4zOXrifYQ98I8Oe89/eYNXDo3mfATyf4dj1piSLB7j9/c14OnAOw28D7QywpiLJliTnj075+RbD3tdTo8W/yXB873mjdReS/NikZ6jh9KhPAO9L8vQkZzIcd9y7fJnhScI3JdmQ5G0MT07tdTLDk1B/neQ04F9PaK4Cfg74xdGx4hOSPBe4GvgOhuOn07bWffsa8IKDvM2fZwj024D/DHx4X+cIJ/l+hsf5XOCs0eUlDGeFXDpa598meVmS45KcwPBE3UPAzqraCfwWcF2SC5KcONrOP9rfYKO//1cBv8/wZOynDvK+TYUB1rQcD/wHhl+n7we+C3jvaNmvMZyZcGOSRxiecHn5lOa4jOHX7vuBa4DfWbX8nzHE5+sMT1T9nxXL3g+8lOE0qf/JEPOxZHiByH7P162q/8Hwa/i/HG37TuBE4Lw1fnWflLXu268y/A/ioSQ/t9aNjZ5A+1cMZz3sYfhtpxjOdlntZ4Hfr6r/V1X3770w/HdxyegwQzH8XT3IcMz2AuC1VfXo6DbeyXAq2gcYDifsZjjD5acYnkDc69dH/419jeFJyY8DF46e8GyX0akZkqR15h6wJDUxwJLUxABLUhMDLElNDvsNNXR027hxY23evLl7DGmmbN++/cGqWlhrPQM85zZv3szS0kG/pF7SASQZ6xWXHoKQpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoc0z2Amu3cCT/4g91TaJ6ddRZcdVX3FC3cA5akJu4Bz7stW+CWW7qnkOaSe8CS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNUlXdM6hRkkeAnd1zTNhG4MHuIabA+3X02FJVJ6+10jHrMYmOaDurarF7iElKsjRr9wm8X0eTJEvjrOchCElqYoAlqYkB1rbuAaZgFu8TeL+OJmPdJ5+Ek6Qm7gFLUhMDLElNDPCcSnJhkp1J7k5yZfc8k5DkQ0keSHJ79yyTkuT0JDcnuTPJHUku755pEpKckOTzSf5sdL/e3z3TpCTZkOSLSf5wrXUN8BxKsgH4IHARcCbwxiRn9k41EdcAF3YPMWFPAu+uqjOBVwDvnJG/qyeA86vq+4CzgAuTvKJ5pkm5HLhrnBUN8Hw6F7i7qr5cVd8GrgN+rHmmw1ZVnwG+0T3HJFXVV6vqC6OvH2H4h31a71SHrwaPjr49dnQ56s8ISLIJeC1w9TjrG+D5dBqwa8X3u5mBf9SzLslm4Gzgc72TTMboV/UdwAPATVU1C/frKuA9wFPjrGyApaNAkpOAjwNXVNXD3fNMQlXtqaqzgE3AuUle0j3T4UhyCfBAVW0f9zoGeD7dB5y+4vtNo5/pCJTkWIb4XltVn+ieZ9Kq6iHgZo7+4/fnAa9Lcg/DYb3zk3zkQFcwwPPp/wIvSvL8JMcBbwD+oHkm7UOSAL8N3FVVH+ieZ1KSLCQ5ZfT1icAFwJd6pzo8VfXeqtpUVZsZ/k19uqredKDrGOA5VFVPApcB/5vhSZ3rq+qO3qkOX5KPAZ8FtiTZneTt3TNNwHnAmxn2pnaMLhd3DzUBpwI3J7mNYYfgpqpa87StWeNLkSWpiXvAktRkKm/IvnHjxtq8efM0bloTtn379geraqF7jsP16lf+8iH9KvczV39q0qMc0HVvvGBdtwdQX1zfo0s3PXVD1nWDR7GpBHjz5s0sLY31hvBqluTe7hmkeeUhCElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJajJWgGfxAxwlqduaAZ7hD3CUpFbj7AHP5Ac4HowrrhgukjRJ47wZz74+wPHlq1dKshXYCnDGGWdMZLgjxY4d3RNImkUTexKuqrZV1WJVLS4sHPXvbihJUzdOgP0AR0magnEC7Ac4StIUrHkMuKqeTLL3Axw3AB+ahQ9wlKRuY30iRlV9Cljfz26RpBnnK+EkqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJajLWCzGkI91N119zSNe7+DU/OdlB1vLnO9d3e8CGZz1r3bep8bgHLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDVZM8BJPpTkgSS3r8dAkjQvxtkDvga4cMpzSNLcWTPAVfUZ4BvrMIskzRWPAUtSk4kFOMnWJEtJlpaXlyd1s5I0syYW4KraVlWLVbW4sLAwqZuVpJnlIQhJajLOaWgfAz4LbEmyO8nbpz+WJM2+NT+SqKreuB6DSNK88RCEJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1WfOVcNLR4KIXvOKQrveX122Y8CQH9siuc9Z1ewAvetfn1n2bGo97wJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1KTcT4V+fQkNye5M8kdSS5fj8EkadaN814QTwLvrqovJDkZ2J7kpqq6c8qzSdJMW3MPuKq+WlVfGH39CHAXcNq0B5OkWXdQx4CTbAbOBv7O2ysl2ZpkKcnS8vLyZKaTpBk2doCTnAR8HLiiqh5evbyqtlXVYlUtLiwsTHJGSZpJYwU4ybEM8b22qj4x3ZEkaT6McxZEgN8G7qqqD0x/JEmaD+PsAZ8HvBk4P8mO0eXiKc8lSTNvzdPQqupWIOswiyTNFV8JJ0lNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1KTcd4PWDriPfnyFx/S9U684fgJT3Jgz33b7nXdno5s7gFLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTcb5VOQTknw+yZ8luSPJ+9djMEmadeO8F8QTwPlV9WiSY4Fbk/yvqvrTKc8mSTNtnE9FLuDR0bfHji41zaEkaR6MdQw4yYYkO4AHgJuq6nP7WGdrkqUkS8vLy5OeU5JmzlgBrqo9VXUWsAk4N8lL9rHOtqparKrFhYWFSc8pSTPnoM6CqKqHgJuBC6czjiTNj3HOglhIcsro6xOBC4AvTXswSZp145wFcSrwu0k2MAT7+qr6w+mOJUmzb5yzIG4Dzl6HWSRprvhKOElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKajPNKOOmId9yuvzqk6z3zT3ZNeJID27Djheu6PYAP3nvrum9T43EPWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmBliSmhhgSWoydoCTbEjyxSR+IKckTcDB7AFfDtw1rUEkad6MFeAkm4DXAldPdxxJmh/j7gFfBbwHeGp/KyTZmmQpydLy8vJEhpOkWbZmgJNcAjxQVdsPtF5VbauqxapaXFhYmNiAkjSrxtkDPg94XZJ7gOuA85N8ZKpTSdIcWDPAVfXeqtpUVZuBNwCfrqo3TX0ySZpxngcsSU0O6iOJquoW4JapTCJJc8Y9YElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaHNQLMaQj1ePfc2hvAHXcPbsmPMmBPXX3Peu6PYDHa8O6b1PjcQ9YkpoYYElqYoAlqYkBlqQmBliSmhhgSWpigCWpiQGWpCYGWJKaGGBJajLWS5FHH0n/CLAHeLKqFqc5lCTNg4N5L4gfqqoHpzaJJM0ZD0FIUpNxA1zAjUm2J9m6rxWSbE2ylGRpeXl5chNK0owaN8A/UFUvBS4C3pnklatXqKptVbVYVYsLC4f21oCSNE/GCnBV3Tf68wHgk8C50xxKkubBmgFO8owkJ+/9Gvhh4PZpDyZJs26csyCeA3wyyd71P1pVfzTVqSRpDqwZ4Kr6MvB96zCLJM0VT0OTpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqYkBlqQmB/OG7NIR67HnHNp/yk9dfM6EJzmw5bc8vq7bA3j3C/es6/Zu/Pa6bu6o5h6wJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1GSvASU5J8ntJvpTkriTfP+3BJGnWjfv6zV8D/qiqXp/kOODpU5xJkubCmgFO8kzglcBbAKrq24Cv9pakwzTOIYjnA8vA7yT5YpKrkzxjynNJ0swbJ8DHAC8FfqOqzgYeA65cvVKSrUmWkiwtLy9PeMxeZ501XCRpksY5Brwb2F1Vnxt9/3vsI8BVtQ3YBrC4uFgTm/AIcNVV3RNImkVr7gFX1f3AriRbRj96NXDnVKeSpDkw7lkQ7wKuHZ0B8WXgrdMbSZLmw1gBrqodwOKUZ5GkueIr4SSpiQGWpCYGWJKaGGBJamKAJamJAZakJgZYkpoYYElqYoAlqUmqJv++OUmWgXsnfsOahudV1UL3ENI8mkqAJUlr8xCEJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1MQAS1ITAyxJTQywJDUxwJLUxABLUhMDLElNDLAkNTHAktTEAEtSEwMsSU0MsCQ1McCS1OT/A4Bsx8/mq+t1AAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ "<Figure size 360x360 with 3 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "pl.figure(4, figsize=(5, 5))\n",
+ "ot.plot.plot1D_mat(a, b, asgd_pi, 'semi-dual : OT matrix ASGD')\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Plot Sinkhorn results\n",
+ "---------------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEc5JREFUeJzt3X2QXQV9xvHnMYRBDII0OxYIuBY1DmMl4IovKKUwYoIW246jUl+Ktc3YWgdaWt86baUztdY6No462AAqAkUpQsdBtGAJQ6kQu5FoCSGWUpHwYjalSFAEEp7+cU/sNi57Tzb33l/23u9nZofde84953fD7HfPnj33XicRAGDwnlI9AACMKgIMAEUIMAAUIcAAUIQAA0ARAgwARQgwsBew/Urbm/qw3TfbvqblumfYvnF3l2HuCDCGQhOIf7f9Y9v32z7X9kHNsk/bfrj5eMz249O+/uoAZovt58y2TpJ/SbJ0jtt/he1v2P6h7Qds/6vtFzfbvSTJKXPZLvqPAGPes322pL+W9MeSDpT0UknPknSt7X2TvDPJoiSLJH1I0hd3fp1kRd3kHbb32YP7Pl3SVZI+IelgSYdJOkfSo72Zrvf25PEOGwKMea0J0DmS3p3ka0keT/I9SW+QNC7pLXPY5om2N9t+j+0ttu+z/au2T7X93eYo8wPT1j/O9k22H2zW/aTtfZtlNzSrfbs54n7jtO2/1/b9kj6787bmPkc2+zi2+fpQ21O2T5xh3OdJUpJLk+xI8kiSa5J8p7nv/zt10ByNv9P2fzTzfsq2n+Tf4W9s32j7wGm3fdT2/9j+L9srpt1+qO0vN3PfYft3pi37oO3LbV9s+yFJZzS3XWb787a32d5ge2I3/1fNewQY893LJe0n6YrpNyZ5WNLVkl41x+3+fLPdwyT9maTz1In5iyS9UtKf2n52s+4OSX8gabGkl0k6WdLvNXOc0KxzdHPE/cVp2z9YnSP1lbvM/p+S3ivpYtv7S/qspAuTXD/DnN+VtMP2hbZX2H5Gi8f2WkkvlvRCdX5QvXr6QttPsX1es/yUJD9sFr1E0qbmcX5E0gXT4v0FSZslHSrp9ZI+ZPukaZt9naTLJR0k6ZLmttOa+x0k6cuSPtli9qFCgDHfLZa0Ncn2GZbd1yyfi8cl/WWSx9WJxGJJH0+yLckGSbdJOlqSkqxLcnOS7c3R999J+qUu239C0p8neTTJI7suTHKepDskrZV0iKQ/mWkjSR6S9ApJUeeHxFRzJPrMWfb94SQPJvm+pDWSlk1btlDSper8cPiVJD+etuyuJOcl2SHpwmauZ9o+XNLxkt6b5CdJ1ks6X9Lbpt33piT/mOSJaY/3xiRXN9u7SM2/5yghwJjvtkpa/CTnFQ9pls/FfzdhkKSdwfjBtOWPSFokSbafZ/uq5o9/D6lznrlb+KeS/KTLOudJeoGkTyR50nO6STYmOSPJkmb9QyWtmmW790/7/Mc7H0fjOeocrZ6T5LEnu9+0MC9q9vdAkm3T1r1Lnd8edrq7xRz7jdr5YQKM+e4mdf7g9OvTb7S9SNIKSf88gBnOlXS7pOcmebqkD0ia8bzqNLO+DGEz/ypJF0j6oO2D2wyS5HZJn1MnxHOxUdLbJX3VdturMu6VdLDtA6bddoSke6aPNsd5hhoBxrzWnJ88R9InbC+3vdD2uKTL1DknedEAxjhA0kOSHrb9fEm/u8vyH0j6hd3c5sclTSb5bUlfkfTpmVay/XzbZ9te0nx9uKTTJd28m/v7qSSXqvND5Ou2j2yx/t2SviHpr2zvZ/uFkt4h6eK5zjAqCDDmvSQfUScYH1UnhGvV+ZX35Nl+de+hP5L0G5K2qXPa4Iu7LP+gpAubqw7e0G1jtl8nabn+L+R/KOlY22+eYfVt6vxxbK3tH6kT3lslnT2Hx/FTSS6U9BeSrmt+oHVzujpXndwr6Up1zm9/fU9mGAXmBdkBoAZHwABQhAADQBECDABFCDAAFBmpi57xsxYvXpzx8fHqMYChsm7duq1JxrqtR4BH3Pj4uCYnJ6vHAIaK7bvarMcpCAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgyD7VA6DYpk3SiSdWT4FRtmyZtGpV9RQlOAIGgCIcAY+6pUul66+vngIYSRwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFHGS6hlQyPY2SZuq5+ixxZK2Vg/RBzyu+WNpkgO6rbTPICbBXm1TkonqIXrJ9uSwPSaJxzWf2J5ssx6nIACgCAEGgCIEGKurB+iDYXxMEo9rPmn1mPgjHAAU4QgYAIoQYAAoQoBHlO3ltjfZvsP2+6rn6QXbn7G9xfat1bP0iu3Dba+xfZvtDbbPrJ6pF2zvZ/ubtr/dPK5zqmfqFdsLbN9i+6pu6xLgEWR7gaRPSVoh6ShJp9s+qnaqnvicpOXVQ/TYdklnJzlK0kslvWtI/l89KumkJEdLWiZpue2XFs/UK2dK2thmRQI8mo6TdEeSO5M8JukLkl5XPNMeS3KDpAeq5+ilJPcl+Vbz+TZ1vrEPq51qz6Xj4ebLhc3HvL8iwPYSSa+RdH6b9QnwaDpM0t3Tvt6sIfimHna2xyUdI2lt7SS90fyqvl7SFknXJhmGx7VK0nskPdFmZQIMzAO2F0n6kqSzkjxUPU8vJNmRZJmkJZKOs/2C6pn2hO3XStqSZF3b+xDg0XSPpMOnfb2kuQ17IdsL1YnvJUmuqJ6n15I8KGmN5v/5++MlnWb7e+qc1jvJ9sWz3YEAj6Z/k/Rc28+2va+kN0n6cvFMmIFtS7pA0sYkH6uep1dsj9k+qPn8qZJeJen22qn2TJL3J1mSZFyd76nrkrxltvsQ4BGUZLuk35f0T+r8UeeyJBtqp9pzti+VdJOkpbY3235H9Uw9cLykt6pzNLW++Ti1eqgeOETSGtvfUeeA4NokXS/bGjY8FRkAinAEDABF+vKC7IsXL874+Hg/No0eW7du3dYkY9Vz7KkTT/nwnH6Ve/Xf3tDrUWa15m3HDXR/kpRbBnt26don/sED3eE81pcAj4+Pa3Ky1QvCo5jtu6pnAEYVpyAAoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaBIqwAP4xs4AkC1rgEe4jdwBIBSbY6Ah/INHHfHWWd1PgCgl9q8GM9Mb+D4kl1Xsr1S0kpJOuKII3oy3N5i/frqCQAMo579ES7J6iQTSSbGxub9qxsCQN+1CTBv4AgAfdAmwLyBIwD0QddzwEm22975Bo4LJH1mGN7AEQCqtXpHjCRXS7q6z7MAwEjhmXAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFGn1RAxgb3fd5y+Y0/1OPeHXejzJ7HLn7QPdnyQt4MWx9locAQNAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFugbY9mdsb7F96yAGAoBR0eYI+HOSlvd5DgAYOV0DnOQGSQ8MYBYAGCmcAwaAIj0LsO2VtidtT05NTfVqswAwtHoW4CSrk0wkmRjj9UcBoCtOQQBAkTaXoV0q6SZJS21vtv2O/o8FAMOv61sSJTl9EIMAwKjhFAQAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABTp+kw4YD5Y8ZyXz+l+379o/x5PMrtH7p0Y6P4k6bnvXjvwfaIdjoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIm3eFflw22ts32Z7g+0zBzEYAAy7Nq8FsV3S2Um+ZfsASetsX5vktj7PBgBDresRcJL7knyr+XybpI2SDuv3YAAw7HbrHLDtcUnHSPqZl1eyvdL2pO3Jqamp3kwHAEOsdYBtL5L0JUlnJXlo1+VJVieZSDIxNjbWyxkBYCi1CrDtherE95IkV/R3JAAYDW2ugrCkCyRtTPKx/o8EAKOhzRHw8ZLeKukk2+ubj1P7PBcADL2ul6EluVGSBzALAIwUngkHAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQJE2rwcM7PV+8sqj5nS/Ay8b7LfAz/3WDwa6P+zdOAIGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAirR5V+T9bH/T9rdtb7B9ziAGA4Bh1+aJ8I9KOinJw7YXSrrR9leT3Nzn2QBgqLV5V+RIerj5cmHzkX4OBQCjoNU5YNsLbK+XtEXStUnWzrDOStuTtienpqZ6PScADJ1WAU6yI8kySUskHWf7BTOsszrJRJKJsbGxXs8JAENnt66CSPKgpDWSlvdnHAAYHW2ughizfVDz+VMlvUrS7f0eDACGXZurIA6RdKHtBeoE+7IkV/V3LAAYfm2ugviOpGMGMAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAkTbPhAP2evtvuG9O99v3nnt7PMns9rl5yUD3J0lfuXf9wPeJdjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIq0DrDtBbZvsc0bcgJAD+zOEfCZkjb2axAAGDWtAmx7iaTXSDq/v+MAwOhoewS8StJ7JD3xZCvYXml70vbk1NRUT4YDgGHWNcC2XytpS5J1s62XZHWSiSQTY2NjPRsQAIZVmyPg4yWdZvt7kr4g6STbF/d1KgAYAV0DnOT9SZYkGZf0JknXJXlL3ycDgCHHdcAAUGS33pIoyfWSru/LJAAwYjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaDIbj0RA9hbPf6sub0AlO+5t8eTzG7HPfcNdH+S9MMnHhno/p4x0L3NbxwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEVaPRW5eUv6bZJ2SNqeZKKfQwHAKNid14L45SRb+zYJAIwYTkEAQJG2AY6ka2yvs71yphVsr7Q9aXtyamqqdxMCwJBqG+BXJDlW0gpJ77J9wq4rJFmdZCLJxNjY3F4aEABGSasAJ7mn+e8WSVdKOq6fQwHAKOgaYNtPs33Azs8lnSLp1n4PBgDDrs1VEM+UdKXtnev/fZKv9XUqABgBXQOc5E5JRw9gFgAYKVyGBgBFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARXbnBdmBvdbWX3zqnO534KIX9XiS2d39m9sHuj9JeuORCwa6v2seGeju5jWOgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoEirANs+yPbltm+3vdH2y/o9GAAMu7ZPRf64pK8leb3tfSXt38eZAGAkdA2w7QMlnSDpDElK8pikx/o7FgAMvzanIJ4taUrSZ23fYvt820/r81wAMPTaBHgfScdKOjfJMZJ+JOl9u65ke6XtSduTU1NTPR6z1rJlnQ8A6KU254A3S9qcZG3z9eWaIcBJVktaLUkTExPp2YR7gVWrqicAMIy6HgEnuV/S3baXNjedLOm2vk4FACOg7VUQ75Z0SXMFxJ2S3t6/kQBgNLQKcJL1kib6PAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEWc9P51c2xPSbqr5xtGPzwryVj1EMAo6kuAAQDdcQoCAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKPK/bk07WnJikdoAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "<Figure size 360x360 with 3 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "pl.figure(4, figsize=(5, 5))\n",
+ "ot.plot.plot1D_mat(a, b, sinkhorn_pi, 'OT matrix Sinkhorn')\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "COMPUTE TRANSPORTATION MATRIX FOR DUAL PROBLEM\n",
+ "############################################################################\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "------------DUAL PROBLEM------------\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(\"------------DUAL PROBLEM------------\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "SEMICONTINOUS CASE\n",
+ "Sample one general measure a, one discrete measures b for the semicontinous\n",
+ "case\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Define one general measure a, one discrete measures b, the points where\n",
+ "are defined the source and the target measures and finally the cost matrix c.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "n_source = 7\n",
+ "n_target = 4\n",
+ "reg = 1\n",
+ "numItermax = 100000\n",
+ "lr = 0.1\n",
+ "batch_size = 3\n",
+ "log = True\n",
+ "\n",
+ "a = ot.utils.unif(n_source)\n",
+ "b = ot.utils.unif(n_target)\n",
+ "\n",
+ "rng = np.random.RandomState(0)\n",
+ "X_source = rng.randn(n_source, 2)\n",
+ "Y_target = rng.randn(n_target, 2)\n",
+ "M = ot.dist(X_source, Y_target)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Call the \"SGD\" dual method to find the transportation matrix in the\n",
+ "semicontinous case\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Call ot.solve_dual_entropic and plot the results.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[ 1.67648902 5.3770004 1.70385554 0.4276547 -0.77206786 1.0474898\n",
+ " 0.54202203] [-0.23723788 -0.20259434 1.30855788 8.06179985]\n",
+ "[[2.62451875e-02 1.00499531e-01 1.78515577e-02 4.57450829e-06]\n",
+ " [1.20510690e-01 1.21972758e-02 1.27002374e-03 7.55197481e-03]\n",
+ " [3.65708350e-03 7.67963231e-02 6.38381061e-02 1.41974930e-07]\n",
+ " [2.64286344e-02 3.31748063e-02 8.24445965e-02 4.25479786e-04]\n",
+ " [9.59295422e-03 7.19190875e-04 1.03739180e-02 1.22100712e-01]\n",
+ " [2.09087627e-02 8.55676046e-04 1.77617241e-03 1.17896019e-01]\n",
+ " [4.18792948e-02 2.63326297e-02 7.17598381e-02 2.49335733e-03]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "sgd_dual_pi, log_sgd = ot.stochastic.solve_dual_entropic(a, b, M, reg,\n",
+ " batch_size, numItermax,\n",
+ " lr, log=log)\n",
+ "print(log_sgd['alpha'], log_sgd['beta'])\n",
+ "print(sgd_dual_pi)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Compare the results with the Sinkhorn algorithm\n",
+ "---------------------------------------------\n",
+ "\n",
+ "Call the Sinkhorn algorithm from POT\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[2.55535622e-02 9.96413843e-02 1.76578860e-02 4.31043335e-06]\n",
+ " [1.21640742e-01 1.25369034e-02 1.30234529e-03 7.37715259e-03]\n",
+ " [3.56096458e-03 7.61460101e-02 6.31500344e-02 1.33788624e-07]\n",
+ " [2.61499607e-02 3.34255577e-02 8.28741973e-02 4.07427179e-04]\n",
+ " [9.85698720e-03 7.52505948e-04 1.08291770e-02 1.21418473e-01]\n",
+ " [2.16947591e-02 9.04086158e-04 1.87228707e-03 1.18386011e-01]\n",
+ " [4.15442692e-02 2.65998963e-02 7.23192701e-02 2.39370724e-03]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "sinkhorn_pi = ot.sinkhorn(a, b, M, reg)\n",
+ "print(sinkhorn_pi)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Plot SGD results\n",
+ "-----------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEgBJREFUeJzt3X2QXQV9xvHnaQhFAWHarBaT0KUjg0U6BLtiLKII1QFx1No3rdBibTMWqTB16mg7Ktp2OjojpvW1KSK2oIiiHccqBSGMhSK4gYDyphZfEhKbjZSaIAoJT/+4J51tTHZPNufe3+6938/MTvbec+49v7OZfPfk3DcnEQBg8H6megAAGFUEGACKEGAAKEKAAaAIAQaAIgQYAIoQYCxoti+1/dcd3M+47dg+oIu5FiLbr7Z9TfUco4QAAx1wz5/b/qbtR2x/z/bf2v7ZZvkXbW9vvh6z/ei0yx/u82ytfrkkuTzJi+a4jZfZXm/7h7a32r7e9lHTlh9t+wrbU80637T9PtvLmuWn2H582s9ko+0rbT9rLvMsFAQY6MbfS1ol6fclHSrpDEmnSbpSkpKckeSQJIdIulzSu3ddTvK6qqF32Z8jf9tPk/RPkt4o6TBJR0n6gKSd05bfImmTpBOSPEnSSZL+U9Jzp93Vpubnc6iklZLulfTvtk+b62zzHQHGgmL7BNu32d5m+5OSDpq27BzbN+62fpoAyPaZtm9vjsA22L6wo5mOlnSupFcnuTnJjiR3SfpNSafbPnUO93mO7Ztsv9f2Q7bvt/1rzfUbbG+x/QfT1p9p377c/PlQc3T5nN3u/weSLpz+82u2tdX28uby8bb/2/bT9zDuCknfTnJderYluSrJ95rlF0q6KcmfJdkoSUm2JFmd5Ird76y5j41J3ibpYknv2tef30JBgLFg2D5Q0r9I+mdJPyfpU+pFrq2H1TtCPVzSmZL+xPbLW277g7Y/uJfFp0namOTW6Vcm2SDpK5JeuA8zTvdsSXdK+nlJH5d0haRnSXqapLMkvd/2Ic26M+3b85o/D2+OuG+edv/3S3qKpL/Zbfb/kPQPkj5m+wmSLpP01iT37mHO2yQ9vYn5C6bNtMuvS7pqn/e+5zOSnmn74Dnefl4jwFhIVkpaLGl1kseSfFrSV9veOMkNSb6W5PEkd0r6hKTnt7ztuUnO3cviJZI272XZ5mb5XHw7yUeT7JT0SUnLJb0zyU+SXCPpUfViPNd925Tkfc0R+yN7WH6heqcUbpX0gHqnFX5KkvslnSJpqXqnXLY2D47uCvESSd/ftb7t85qj+u22/3G2GSVZvV8sQ4cAYyF5qqQH8v/fQeq7bW9s+9m21zYPBP2PpNdp7nGcbqukI/ay7Ihm+Vz817TvH5GkJLtfd4g0533bMNPCJI9JulTScZLes9vPffd1v5Lkd5KMSTpZvaPuv2wW/0DTfj5J3p/kcEmr1fuFOpOlkiLpoVnWW5AIMBaSzZKW2va0646c9v3Dkp6464LtX9jt9h+X9DlJy5McJunD6h1d7a/rJS23feL0K5vzpyslXdfBNmYz077tLZwzvhWi7aWS3i7po5Les+sZHbNJ8lX1Th0c11x1naRXtLntHvyGpNuSPDzH289rBBgLyc2Sdkh6g+3Ftl8haXr07pD0DNsrbB+k3n+hpztU0oNJftzE8ve6GCrJN9QL3uW2V9peZPsZ6p33/FKSL3WxnVnMtG9Tkh6X9Ett76z5JXeppI9Ieq16v/z+ai/rPtf2H9t+cnP56ZJeqt75b6n393Cy7YuaqMv2Ekm/vLdt215q++2S/kjSX7Sde6EhwFgwkjyq3pHUOZIelPS76h1p7Vr+DUnvlPQlSd+UdONud3GupHfa3ibpbWqeItaG7Q975ufrnqfeI/aXSdou6WpJN2jfHiTcH3vdtyQ/Uu9Btpuac68rW9zfGyQ9Wb0H3iLpNZJeY/vkPaz7kHrB/ZrtXfv+WUnvbrb/DfUe8Fsm6Y5mxpvUO7/71mn389Tm9tvVO7f/K5JOac53DyXzhuwAUIMjYAAoQoABoAgBBoAiBBgAiozsW++hZ8mSJRkfH68eAxgq69at29q8KGVGBHjEjY+Pa3JysnoMYKjYbvUKTU5BAEARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFDqgeAMXuu0865ZTqKTDKVqyQVq+unqIER8AAUIQj4FF3zDHSDTdUTwGMJI6AAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAijhJ9QwoZHubpPuq5+jYEklbq4foA/Zr4TgmyaGzrXTAICbBvHZfkonqIbpke3LY9klivxYS25Nt1uMUBAAUIcAAUIQAY031AH0wjPsksV8LSat94kE4ACjCETAAFCHAAFCEAI8o26fbvs/2t2y/uXqeLti+xPYW21+vnqUrtpfbXmv7btt32T6/eqYu2D7I9q2272j26x3VM3XF9iLbt9v+/GzrEuARZHuRpA9IOkPSsZJeZfvY2qk6camk06uH6NgOSW9McqyklZJePyR/Vz+RdGqS4yWtkHS67ZXFM3XlfEn3tFmRAI+mEyV9K8n9SR6VdIWklxXPtN+SfFnSg9VzdCnJ5iS3Nd9vU+8f9tLaqfZferY3Fxc3Xwv+GQG2l0k6U9LFbdYnwKNpqaQN0y5v1BD8ox52tsclnSDpltpJutH8V329pC2Srk0yDPu1WtKbJD3eZmUCDCwAtg+RdJWkC5L8sHqeLiTZmWSFpGWSTrR9XPVM+8P2SyRtSbKu7W0I8Gh6QNLyaZeXNddhHrK9WL34Xp7kM9XzdC3JQ5LWauGfvz9J0kttf0e903qn2r5sphsQ4NH0VUlH2z7K9oGSXinpc8UzYQ9sW9JHJN2T5KLqebpie8z24c33T5D0Qkn31k61f5K8JcmyJOPq/Zu6PslZM92GAI+gJDsknSfp39R7UOfKJHfVTrX/bH9C0s2SjrG90fZrq2fqwEmSzlbvaGp98/Xi6qE6cISktbbvVO+A4Noksz5ta9jwUmQAKMIRMAAU6csbsi9ZsiTj4+P9uGt0bN26dVuTjFXPsb+ef8a75vRfuRe/Z23Xo8zourMH/1qD3D7Ys0vXPv4pD3SDC1hfAjw+Pq7JyVZvCI9itr9bPQMwqjgFAQBFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARVoFeBg/wBEAqs0a4CH+AEcAKNXmCHgoP8BxX1xwQe8LALrU5s149vQBjs/efSXbqyStkqQjjzyyk+Hmi/XrqycAMIw6exAuyZokE0kmxsYW/LsbAkDftQkwH+AIAH3QJsB8gCMA9MGs54CT7LC96wMcF0m6ZBg+wBEAqrX6RIwkX5D0hT7PAgAjhVfCAUARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCk1QsxgPnu6ks+NKfbveLk3+54kpnlO/cOdHuStIg3x5q3OAIGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAiswaYNuX2N5i++uDGAgARkWbI+BLJZ3e5zkAYOTMGuAkX5b04ABmAYCRwjlgACjSWYBtr7I9aXtyamqqq7sFgKHVWYCTrEkykWRijPcfBYBZcQoCAIq0eRraJyTdLOkY2xttv7b/YwHA8Jv1I4mSvGoQgwDAqOEUBAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFJn1lXDAQvDyo58/p9ttuOzgjieZ2Y82TQx0e5J09Hm3DHybaIcjYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaBIm09FXm57re27bd9l+/xBDAYAw67Ne0HskPTGJLfZPlTSOtvXJrm7z7MBwFCb9Qg4yeYktzXfb5N0j6Sl/R4MAIbdPp0Dtj0u6QRJP/X2SrZX2Z60PTk1NdXNdAAwxFoH2PYhkq6SdEGSH+6+PMmaJBNJJsbGxrqcEQCGUqsA216sXnwvT/KZ/o4EAKOhzbMgLOkjku5JclH/RwKA0dDmCPgkSWdLOtX2+ubrxX2eCwCG3qxPQ0tyoyQPYBYAGCm8Eg4AihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIm3eDxiY93588rFzut1hVwz2n8Bhf7hloNvD/MYRMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCkzaciH2T7Vtt32L7L9jsGMRgADLs2L4T/iaRTk2y3vVjSjba/mOQrfZ4NAIZam09FjqTtzcXFzVf6ORQAjIJW54BtL7K9XtIWSdcmuWUP66yyPWl7cmpqqus5AWDotApwkp1JVkhaJulE28ftYZ01SSaSTIyNjXU9JwAMnX16FkSShyStlXR6f8YBgNHR5lkQY7YPb75/gqQXSrq334MBwLBr8yyIIyR9zPYi9YJ9ZZLP93csABh+bZ4FcaekEwYwCwCMFF4JBwBFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0CRNq+EA+a9J961eU63O/CBTR1PMrMDbl460O1J0r9uWj/wbaIdjoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIq0DbHuR7dtt84GcANCBfTkCPl/SPf0aBABGTasA214m6UxJF/d3HAAYHW2PgFdLepOkx/e2gu1VtidtT05NTXUyHAAMs1kDbPslkrYkWTfTeknWJJlIMjE2NtbZgAAwrNocAZ8k6aW2vyPpCkmn2r6sr1MBwAiYNcBJ3pJkWZJxSa+UdH2Ss/o+GQAMOZ4HDABF9ukjiZLcIOmGvkwCACOGI2AAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAiuzTCzGA+eqx5UvmdDs/sKnjSWa2c/P3B7o9Sdq68+GBbu/JA93awsYRMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCk1UuRm4+k3yZpp6QdSSb6ORQAjIJ9eS+IFyTZ2rdJAGDEcAoCAIq0DXAkXWN7ne1Ve1rB9irbk7Ynp6amupsQAIZU2wA/N8kzJZ0h6fW2n7f7CknWJJlIMjE2NtbpkAAwjFoFOMkDzZ9bJH1W0on9HAoARsGsAbZ9sO1Dd30v6UWSvt7vwQBg2LV5FsRTJH3W9q71P57k6r5OBQAjYNYAJ7lf0vEDmAUARgpPQwOAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCL78obswLy19fgnzul2hz3pVzueZGbfO3vnQLcnSWc9bdFAt3fNIwPd3ILGETAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABRpFWDbh9v+tO17bd9j+zn9HgwAhl3blyL/naSrk/yW7QMlze11nwCA/zNrgG0fJul5ks6RpCSPSnq0v2MBwPBrcwriKElTkj5q+3bbF9s+uM9zAcDQaxPgAyQ9U9KHkpwg6WFJb959JdurbE/anpyamup4zForVvS+AKBLbc4Bb5S0McktzeVPaw8BTrJG0hpJmpiYSGcTzgOrV1dPAGAYzXoEnOT7kjbYPqa56jRJd/d1KgAYAW2fBfGnki5vngFxv6TX9G8kABgNrQKcZL2kiT7PAgAjhVfCAUARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAESfdv2+O7SlJ3+38jtEPv5hkrHoIYBT1JcAAgNlxCgIAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAo8r9wCGj9yW4UbQAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ "<Figure size 360x360 with 3 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "pl.figure(4, figsize=(5, 5))\n",
+ "ot.plot.plot1D_mat(a, b, sgd_dual_pi, 'dual : OT matrix SGD')\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Plot Sinkhorn results\n",
+ "---------------------\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEc5JREFUeJzt3X2QXQV9xvHnMYRBDII0OxYIuBY1DmMl4IovKKUwYoIW246jUl+Ktc3YWgdaWt86baUztdY6No462AAqAkUpQsdBtGAJQ6kQu5FoCSGWUpHwYjalSFAEEp7+cU/sNi57Tzb33l/23u9nZofde84953fD7HfPnj33XicRAGDwnlI9AACMKgIMAEUIMAAUIcAAUIQAA0ARAgwARQgwsBew/Urbm/qw3TfbvqblumfYvnF3l2HuCDCGQhOIf7f9Y9v32z7X9kHNsk/bfrj5eMz249O+/uoAZovt58y2TpJ/SbJ0jtt/he1v2P6h7Qds/6vtFzfbvSTJKXPZLvqPAGPes322pL+W9MeSDpT0UknPknSt7X2TvDPJoiSLJH1I0hd3fp1kRd3kHbb32YP7Pl3SVZI+IelgSYdJOkfSo72Zrvf25PEOGwKMea0J0DmS3p3ka0keT/I9SW+QNC7pLXPY5om2N9t+j+0ttu+z/au2T7X93eYo8wPT1j/O9k22H2zW/aTtfZtlNzSrfbs54n7jtO2/1/b9kj6787bmPkc2+zi2+fpQ21O2T5xh3OdJUpJLk+xI8kiSa5J8p7nv/zt10ByNv9P2fzTzfsq2n+Tf4W9s32j7wGm3fdT2/9j+L9srpt1+qO0vN3PfYft3pi37oO3LbV9s+yFJZzS3XWb787a32d5ge2I3/1fNewQY893LJe0n6YrpNyZ5WNLVkl41x+3+fLPdwyT9maTz1In5iyS9UtKf2n52s+4OSX8gabGkl0k6WdLvNXOc0KxzdHPE/cVp2z9YnSP1lbvM/p+S3ivpYtv7S/qspAuTXD/DnN+VtMP2hbZX2H5Gi8f2WkkvlvRCdX5QvXr6QttPsX1es/yUJD9sFr1E0qbmcX5E0gXT4v0FSZslHSrp9ZI+ZPukaZt9naTLJR0k6ZLmttOa+x0k6cuSPtli9qFCgDHfLZa0Ncn2GZbd1yyfi8cl/WWSx9WJxGJJH0+yLckGSbdJOlqSkqxLcnOS7c3R999J+qUu239C0p8neTTJI7suTHKepDskrZV0iKQ/mWkjSR6S9ApJUeeHxFRzJPrMWfb94SQPJvm+pDWSlk1btlDSper8cPiVJD+etuyuJOcl2SHpwmauZ9o+XNLxkt6b5CdJ1ks6X9Lbpt33piT/mOSJaY/3xiRXN9u7SM2/5yghwJjvtkpa/CTnFQ9pls/FfzdhkKSdwfjBtOWPSFokSbafZ/uq5o9/D6lznrlb+KeS/KTLOudJeoGkTyR50nO6STYmOSPJkmb9QyWtmmW790/7/Mc7H0fjOeocrZ6T5LEnu9+0MC9q9vdAkm3T1r1Lnd8edrq7xRz7jdr5YQKM+e4mdf7g9OvTb7S9SNIKSf88gBnOlXS7pOcmebqkD0ia8bzqNLO+DGEz/ypJF0j6oO2D2wyS5HZJn1MnxHOxUdLbJX3VdturMu6VdLDtA6bddoSke6aPNsd5hhoBxrzWnJ88R9InbC+3vdD2uKTL1DknedEAxjhA0kOSHrb9fEm/u8vyH0j6hd3c5sclTSb5bUlfkfTpmVay/XzbZ9te0nx9uKTTJd28m/v7qSSXqvND5Ou2j2yx/t2SviHpr2zvZ/uFkt4h6eK5zjAqCDDmvSQfUScYH1UnhGvV+ZX35Nl+de+hP5L0G5K2qXPa4Iu7LP+gpAubqw7e0G1jtl8nabn+L+R/KOlY22+eYfVt6vxxbK3tH6kT3lslnT2Hx/FTSS6U9BeSrmt+oHVzujpXndwr6Up1zm9/fU9mGAXmBdkBoAZHwABQhAADQBECDABFCDAAFBmpi57xsxYvXpzx8fHqMYChsm7duq1JxrqtR4BH3Pj4uCYnJ6vHAIaK7bvarMcpCAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgyD7VA6DYpk3SiSdWT4FRtmyZtGpV9RQlOAIGgCIcAY+6pUul66+vngIYSRwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFHGS6hlQyPY2SZuq5+ixxZK2Vg/RBzyu+WNpkgO6rbTPICbBXm1TkonqIXrJ9uSwPSaJxzWf2J5ssx6nIACgCAEGgCIEGKurB+iDYXxMEo9rPmn1mPgjHAAU4QgYAIoQYAAoQoBHlO3ltjfZvsP2+6rn6QXbn7G9xfat1bP0iu3Dba+xfZvtDbbPrJ6pF2zvZ/ubtr/dPK5zqmfqFdsLbN9i+6pu6xLgEWR7gaRPSVoh6ShJp9s+qnaqnvicpOXVQ/TYdklnJzlK0kslvWtI/l89KumkJEdLWiZpue2XFs/UK2dK2thmRQI8mo6TdEeSO5M8JukLkl5XPNMeS3KDpAeq5+ilJPcl+Vbz+TZ1vrEPq51qz6Xj4ebLhc3HvL8iwPYSSa+RdH6b9QnwaDpM0t3Tvt6sIfimHna2xyUdI2lt7SS90fyqvl7SFknXJhmGx7VK0nskPdFmZQIMzAO2F0n6kqSzkjxUPU8vJNmRZJmkJZKOs/2C6pn2hO3XStqSZF3b+xDg0XSPpMOnfb2kuQ17IdsL1YnvJUmuqJ6n15I8KGmN5v/5++MlnWb7e+qc1jvJ9sWz3YEAj6Z/k/Rc28+2va+kN0n6cvFMmIFtS7pA0sYkH6uep1dsj9k+qPn8qZJeJen22qn2TJL3J1mSZFyd76nrkrxltvsQ4BGUZLuk35f0T+r8UeeyJBtqp9pzti+VdJOkpbY3235H9Uw9cLykt6pzNLW++Ti1eqgeOETSGtvfUeeA4NokXS/bGjY8FRkAinAEDABF+vKC7IsXL874+Hg/No0eW7du3dYkY9Vz7KkTT/nwnH6Ve/Xf3tDrUWa15m3HDXR/kpRbBnt26don/sED3eE81pcAj4+Pa3Ky1QvCo5jtu6pnAEYVpyAAoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaBIqwAP4xs4AkC1rgEe4jdwBIBSbY6Ah/INHHfHWWd1PgCgl9q8GM9Mb+D4kl1Xsr1S0kpJOuKII3oy3N5i/frqCQAMo579ES7J6iQTSSbGxub9qxsCQN+1CTBv4AgAfdAmwLyBIwD0QddzwEm22975Bo4LJH1mGN7AEQCqtXpHjCRXS7q6z7MAwEjhmXAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFGn1RAxgb3fd5y+Y0/1OPeHXejzJ7HLn7QPdnyQt4MWx9locAQNAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQBECDABFugbY9mdsb7F96yAGAoBR0eYI+HOSlvd5DgAYOV0DnOQGSQ8MYBYAGCmcAwaAIj0LsO2VtidtT05NTfVqswAwtHoW4CSrk0wkmRjj9UcBoCtOQQBAkTaXoV0q6SZJS21vtv2O/o8FAMOv61sSJTl9EIMAwKjhFAQAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARQgwABTp+kw4YD5Y8ZyXz+l+379o/x5PMrtH7p0Y6P4k6bnvXjvwfaIdjoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIm3eFflw22ts32Z7g+0zBzEYAAy7Nq8FsV3S2Um+ZfsASetsX5vktj7PBgBDresRcJL7knyr+XybpI2SDuv3YAAw7HbrHLDtcUnHSPqZl1eyvdL2pO3Jqamp3kwHAEOsdYBtL5L0JUlnJXlo1+VJVieZSDIxNjbWyxkBYCi1CrDtherE95IkV/R3JAAYDW2ugrCkCyRtTPKx/o8EAKOhzRHw8ZLeKukk2+ubj1P7PBcADL2ul6EluVGSBzALAIwUngkHAEUIMAAUIcAAUIQAA0ARAgwARQgwABQhwABQhAADQJE2rwcM7PV+8sqj5nS/Ay8b7LfAz/3WDwa6P+zdOAIGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAirR5V+T9bH/T9rdtb7B9ziAGA4Bh1+aJ8I9KOinJw7YXSrrR9leT3Nzn2QBgqLV5V+RIerj5cmHzkX4OBQCjoNU5YNsLbK+XtEXStUnWzrDOStuTtienpqZ6PScADJ1WAU6yI8kySUskHWf7BTOsszrJRJKJsbGxXs8JAENnt66CSPKgpDWSlvdnHAAYHW2ughizfVDz+VMlvUrS7f0eDACGXZurIA6RdKHtBeoE+7IkV/V3LAAYfm2ugviOpGMGMAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAkTbPhAP2evtvuG9O99v3nnt7PMns9rl5yUD3J0lfuXf9wPeJdjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIq0DrDtBbZvsc0bcgJAD+zOEfCZkjb2axAAGDWtAmx7iaTXSDq/v+MAwOhoewS8StJ7JD3xZCvYXml70vbk1NRUT4YDgGHWNcC2XytpS5J1s62XZHWSiSQTY2NjPRsQAIZVmyPg4yWdZvt7kr4g6STbF/d1KgAYAV0DnOT9SZYkGZf0JknXJXlL3ycDgCHHdcAAUGS33pIoyfWSru/LJAAwYjgCBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaDIbj0RA9hbPf6sub0AlO+5t8eTzG7HPfcNdH+S9MMnHhno/p4x0L3NbxwBA0ARAgwARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEVaPRW5eUv6bZJ2SNqeZKKfQwHAKNid14L45SRb+zYJAIwYTkEAQJG2AY6ka2yvs71yphVsr7Q9aXtyamqqdxMCwJBqG+BXJDlW0gpJ77J9wq4rJFmdZCLJxNjY3F4aEABGSasAJ7mn+e8WSVdKOq6fQwHAKOgaYNtPs33Azs8lnSLp1n4PBgDDrs1VEM+UdKXtnev/fZKv9XUqABgBXQOc5E5JRw9gFgAYKVyGBgBFCDAAFCHAAFCEAANAEQIMAEUIMAAUIcAAUIQAA0ARAgwARXbnBdmBvdbWX3zqnO534KIX9XiS2d39m9sHuj9JeuORCwa6v2seGeju5jWOgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoEirANs+yPbltm+3vdH2y/o9GAAMu7ZPRf64pK8leb3tfSXt38eZAGAkdA2w7QMlnSDpDElK8pikx/o7FgAMvzanIJ4taUrSZ23fYvt820/r81wAMPTaBHgfScdKOjfJMZJ+JOl9u65ke6XtSduTU1NTPR6z1rJlnQ8A6KU254A3S9qcZG3z9eWaIcBJVktaLUkTExPp2YR7gVWrqicAMIy6HgEnuV/S3baXNjedLOm2vk4FACOg7VUQ75Z0SXMFxJ2S3t6/kQBgNLQKcJL1kib6PAsAjBSeCQcARQgwABQhwABQhAADQBECDABFCDAAFCHAAFCEAANAEQIMAEWc9P51c2xPSbqr5xtGPzwryVj1EMAo6kuAAQDdcQoCAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKEKAAaAIAQaAIgQYAIoQYAAoQoABoAgBBoAiBBgAihBgAChCgAGgCAEGgCIEGACKEGAAKPK/bk07WnJikdoAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ "<Figure size 360x360 with 3 Axes>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "pl.figure(4, figsize=(5, 5))\n",
+ "ot.plot.plot1D_mat(a, b, sinkhorn_pi, 'OT matrix Sinkhorn')\n",
+ "pl.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.6.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}