ipm-prerequisite.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Empirical estimation of IPMs (integral probability metrics): prerequisite<span class=\"tocSkip\"></span></h1>\n",
    "\n",
    "Author: [Sylvain Combettes](https://github.com/sylvaincom).\n",
    "\n",
    "Last update: Feb 3, 2020.\n",
    "\n",
    "GitHub repository: [here](https://github.com/sylvaincom/comparison-distributions).\n",
    "\n",
    "---\n",
    "This notebook is the prerequisite of my notebook `ipm.ipynb` about the empirical estimation of IPMs (integral probability metrics) and completes my [report](https://sylvaincom.github.io/docs/comparison_distributions_report.pdf) on the _Comparison of Empirical Probability Distributions_. One can consider this notebook as an introduction to `ipm.ipynb`.\n",
    "\n",
    "As $f$-divergences input probability distributions, IPMs input samples drawn from the probability distributions. Actually, for IPMs, we only focus on the Kantorovich metric.\n",
    "\n",
    "<br/>\n",
    "\n",
    "<div class=\"alert alert-info\"><h4>README<span class=\"tocSkip\"></span></h4>\n",
    "\n",
    "<p> I recommend reading the report before reading this notebook, as most explanations are not duplicated here.</p>\n",
    "\n",
    "<p>\n",
    "The best way to open this Jupyter Notebook is to use the table of contents from the extensions called <code>nbextensions</code>. See <a href=\"https://towardsdatascience.com/4-awesome-tips-for-enhancing-jupyter-notebooks-4d8905f926c5\">4 Awesome Tips for Enhancing Jupyter Notebooks</a> by George Seif.\n",
    "    \n",
    "The Python version is 3.7.3.\n",
    "</p></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Preliminary:-Solving-a-linear-programming-problem-using-PuLP\" data-toc-modified-id=\"Preliminary:-Solving-a-linear-programming-problem-using-PuLP-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Preliminary: Solving a linear programming problem using <code>PuLP</code></a></span></li><li><span><a href=\"#Example-of-empirical-estimation-of-the-Kantorovich-metric-$W$\" data-toc-modified-id=\"Example-of-empirical-estimation-of-the-Kantorovich-metric-$W$-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Example of empirical estimation of the Kantorovich metric $W$</a></span><ul class=\"toc-item\"><li><span><a href=\"#Getting-the-canonical-form-of-the-linear-programming-problem\" data-toc-modified-id=\"Getting-the-canonical-form-of-the-linear-programming-problem-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Getting the canonical form of the linear programming problem</a></span></li><li><span><a href=\"#Using-PuLP-to-solve-the-obtained-linear-programming-problem\" data-toc-modified-id=\"Using-PuLP-to-solve-the-obtained-linear-programming-problem-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Using <code>PuLP</code> to solve the obtained linear programming problem</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Imports<span class=\"tocSkip\"></span></h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pip install pulp\n",
    "from pulp import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preliminary: Solving a linear programming problem using `PuLP`\n",
    "\n",
    "As I am new to `PuLP`, this section constitutes a short introduction of how `PuLP` works.\n",
    "This section is inspired from the official `PuLP` [documentation](https://pypi.org/project/PuLP/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `LpVariable()` to create new variables. To create a variable $0 \\leq x \\leq 3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = LpVariable(\"x\", 0, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a variable $0 \\leq y \\leq 1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = LpVariable(\"y\", 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `LpProblem()` to create new problems. Create `“myProblem”`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = LpProblem(\"myProblem\", LpMinimize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine variables to create expressions and constraints, then add them to the problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob += x + y <= 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you add an expression (not a constraint), it will become the objective:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob += -4*x + y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `prob` variable contains all the data of our linear programming problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "myProblem:\n",
      "MINIMIZE\n",
      "-4*x + 1*y + 0\n",
      "SUBJECT TO\n",
      "_C1: x + y <= 2\n",
      "\n",
      "VARIABLES\n",
      "x <= 3 Continuous\n",
      "y <= 1 Continuous\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(prob)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve with the default included solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "status = prob.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the status of the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Optimal'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LpStatus[status]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can get the value of the variables using `value()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The solution is (2.00, 0.00) and the objective function is -8.00\n"
     ]
    }
   ],
   "source": [
    "print('The solution is (%1.2f, %1.2f) and the objective function is %1.2f'\n",
    "      % (value(x), value(y), value(prob.objective)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check this result graphically using [GeoGebra](https://www.geogebra.org/) and this [tutorial](https://www.youtube.com/watch?v=4vvlmSJqJto). With GeoGebra, we obtained the same result as with `PuLP`:\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/sylvaincom/comparison-distributions/master/img/geogebra.png\" width=\"600\" height=\"539.6\">\n",
    "\n",
    "For more information on `PuLP`, see examples on [GitHub](https://github.com/coin-or/pulp/tree/master/examples)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example of empirical estimation of the Kantorovich metric $W$\n",
    "\n",
    "We take the example from the report at subsection II.3.2. Our linear programming problem in the canonical form is:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{cases}\n",
    "    \\max \\left(c^Ta\\right) \\\\\n",
    "    Ma \\leq b\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "\n",
    "With $N=4$, we have:\n",
    "\n",
    "\\begin{equation}\n",
    "c =\n",
    "\\begin{pmatrix}\n",
    "    \\widetilde{Y}_1 \\\\\n",
    "    \\widetilde{Y}_2 \\\\\n",
    "    \\widetilde{Y}_3 \\\\\n",
    "    \\widetilde{Y}_4\n",
    "\\end{pmatrix} \\quad \\quad\n",
    "M = \\begin{pmatrix}\n",
    "        1 & -1 & 0 & 0 \\\\\n",
    "        -1 & 1 & 0 & 0 \\\\\n",
    "        1 & 0 & -1 & 0 \\\\\n",
    "        -1 & 0 & 1 & 0 \\\\\n",
    "        1 & 0 & 0 & -1 \\\\\n",
    "        -1 & 0 & 0 & 1 \\\\\n",
    "        0 & 1 & -1 & 0 \\\\\n",
    "        0 & -1 & 1 & 0 \\\\\n",
    "        0 & 1 & 0 & -1 \\\\\n",
    "        0 & -1 & 0 & 1 \\\\\n",
    "        0 & 0 & 1 & -1 \\\\\n",
    "        0 & 0 & -1 & 1\n",
    "    \\end{pmatrix} \\quad \\quad\n",
    "b = \\begin{pmatrix}\n",
    "        \\rho\\left(X_1, X_2\\right) \\\\\n",
    "        \\rho\\left(X_1, X_2\\right) \\\\\n",
    "        \\rho\\left(X_1, X_3\\right) \\\\\n",
    "        \\rho\\left(X_1, X_3\\right) \\\\\n",
    "        \\rho\\left(X_1, X_4\\right) \\\\\n",
    "        \\rho\\left(X_1, X_4\\right) \\\\\n",
    "        \\rho\\left(X_2, X_3\\right) \\\\\n",
    "        \\rho\\left(X_2, X_3\\right) \\\\\n",
    "        \\rho\\left(X_2, X_4\\right) \\\\\n",
    "        \\rho\\left(X_2, X_4\\right) \\\\\n",
    "        \\rho\\left(X_3, X_4\\right) \\\\\n",
    "        \\rho\\left(X_3, X_4\\right)\n",
    "    \\end{pmatrix}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We take the same uniform distributions as in the report at subsection II.1.3.\n",
    "\n",
    "We choose the samples $\\left\\{X_1^{(1)}, X_2^{(1)}, \\ldots, X_m^{(1)}\\right\\}$ drawn randomly from $\\mathbb{P}$, noted `X_p`, to be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_p is: \n",
      " [0.417022   0.72032449]\n"
     ]
    }
   ],
   "source": [
    "m = 2\n",
    "a = 0\n",
    "b = 1\n",
    "np.random.seed(1) # random seed for reproducability\n",
    "X_p = np.random.uniform(a, b, m)\n",
    "print('X_p is: \\n', X_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose the samples $\\left\\{X_1^{(2)}, X_2^{(2)}, \\ldots, X_n^{(2)}\\right\\}$ drawn randomly from $\\mathbb{Q}$, noted `X_q`, to be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_q is: \n",
      " [0.917022   1.22032449]\n"
     ]
    }
   ],
   "source": [
    "n = 2\n",
    "r = 1/2\n",
    "s = 3/2\n",
    "np.random.seed(1) # random seed for reproducability\n",
    "X_q = np.random.uniform(r, s, n)\n",
    "print('X_q is: \\n', X_q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X is: \n",
      " [0.417022   0.72032449 0.917022   1.22032449] \n",
      "\n",
      "N is: \n",
      " 4\n"
     ]
    }
   ],
   "source": [
    "X = np.concatenate((X_p, X_q), axis=0)\n",
    "print('X is: \\n', X, '\\n')\n",
    "\n",
    "N = m+n\n",
    "print('N is: \\n', N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting the canonical form of the linear programming problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_c(X_p, X_q):\n",
    "    m = len(X_p)\n",
    "    n = len(X_q)\n",
    "    \n",
    "    Y = [1/m]*m + [-1/n]*n\n",
    "    \n",
    "    return np.asarray(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c is: \n",
      " [ 0.5  0.5 -0.5 -0.5]\n"
     ]
    }
   ],
   "source": [
    "c = construct_c(X_p, X_q)\n",
    "print('c is: \\n', c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rho(x,y):\n",
    "    return abs(x-y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_b(X_p, X_q):\n",
    "    \n",
    "    X = np.concatenate((X_p, X_q), axis=0)\n",
    "    N = len(X)\n",
    "    \n",
    "    b_part = []\n",
    "    for i in range(N):\n",
    "        for j in range(i+1, N):\n",
    "            b_part.append(rho(X[i], X[j]))\n",
    "    \n",
    "    # Now, we duplicate each row to obtain a list of size 2*N\n",
    "    b = []\n",
    "    for i in range(N):\n",
    "        b.append(b_part[i])\n",
    "        b.append(b_part[i])\n",
    "    \n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b is: \n",
      " [0.3033024887395841, 0.3033024887395841, 0.5, 0.5, 0.803302488739584, 0.803302488739584, 0.1966975112604159, 0.1966975112604159]\n"
     ]
    }
   ],
   "source": [
    "b = construct_b(X_p, X_q)\n",
    "print('b is: \\n', b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_M(X_p, X_q):\n",
    "    \n",
    "    X = np.concatenate((X_p, X_q), axis=0)\n",
    "    N = len(X)\n",
    "    \n",
    "    M = []\n",
    "    for i in range(N):\n",
    "        for j in range(i+1, N):\n",
    "            l_M_1 = [0]*N\n",
    "            l_M_1[i] = 1\n",
    "            l_M_1[j] = -1\n",
    "            M.append(l_M_1)\n",
    "            l_M_2 = [0]*N\n",
    "            l_M_2[i] = -1\n",
    "            l_M_2[j] = 1\n",
    "            M.append(l_M_2)\n",
    "    M = np.asarray(M)\n",
    "\n",
    "    return M.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M is: \n",
      " [[ 1 -1  0  0]\n",
      " [-1  1  0  0]\n",
      " [ 1  0 -1  0]\n",
      " [-1  0  1  0]\n",
      " [ 1  0  0 -1]\n",
      " [-1  0  0  1]\n",
      " [ 0  1 -1  0]\n",
      " [ 0 -1  1  0]\n",
      " [ 0  1  0 -1]\n",
      " [ 0 -1  0  1]\n",
      " [ 0  0  1 -1]\n",
      " [ 0  0 -1  1]]\n"
     ]
    }
   ],
   "source": [
    "M = construct_M(X_p, X_q)\n",
    "print('M is: \\n', M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using `PuLP` to solve the obtained linear programming problem\n",
    "\n",
    "We rely on [test5](https://github.com/coin-or/pulp/blob/master/examples/test5.py) of the official `PuLP` examples to define our LP (linear programming) problem with matrices.\n",
    "\n",
    "From the previous subsection, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N is: \n",
      " 4 \n",
      "\n",
      "c is: \n",
      " [ 0.5  0.5 -0.5 -0.5] \n",
      "\n",
      "M is: \n",
      " [[ 1 -1  0  0]\n",
      " [-1  1  0  0]\n",
      " [ 1  0 -1  0]\n",
      " [-1  0  1  0]\n",
      " [ 1  0  0 -1]\n",
      " [-1  0  0  1]\n",
      " [ 0  1 -1  0]\n",
      " [ 0 -1  1  0]\n",
      " [ 0  1  0 -1]\n",
      " [ 0 -1  0  1]\n",
      " [ 0  0  1 -1]\n",
      " [ 0  0 -1  1]] \n",
      "\n",
      "b is: \n",
      " [0.3033024887395841, 0.3033024887395841, 0.5, 0.5, 0.803302488739584, 0.803302488739584, 0.1966975112604159, 0.1966975112604159]\n"
     ]
    }
   ],
   "source": [
    "print('N is: \\n', N, '\\n')\n",
    "print('c is: \\n', c, '\\n')\n",
    "print('M is: \\n', M, '\\n')\n",
    "print('b is: \\n', b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define our LP problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = LpProblem(\"LP_problem_for_estimating_the_Kantorovich_metric\", LpMaximize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create our vector of $N$ continuous variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = LpVariable.matrix(\"a\", list(range(N)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define our objective function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob += lpDot(c, a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define our number of constraints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 2*N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define our constraints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(p):\n",
    "    prob += lpDot(M[i], a) <= b[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print our LP problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LP_problem_for_estimating_the_Kantorovich_metric:\n",
      "MAXIMIZE\n",
      "0.5*a_0 + 0.5*a_1 + -0.5*a_2 + -0.5*a_3 + 0.0\n",
      "SUBJECT TO\n",
      "_C1: a_0 - a_1 <= 0.30330248874\n",
      "\n",
      "_C2: - a_0 + a_1 <= 0.30330248874\n",
      "\n",
      "_C3: a_0 - a_2 <= 0.5\n",
      "\n",
      "_C4: - a_0 + a_2 <= 0.5\n",
      "\n",
      "_C5: a_0 - a_3 <= 0.80330248874\n",
      "\n",
      "_C6: - a_0 + a_3 <= 0.80330248874\n",
      "\n",
      "_C7: a_1 - a_2 <= 0.19669751126\n",
      "\n",
      "_C8: - a_1 + a_2 <= 0.19669751126\n",
      "\n",
      "VARIABLES\n",
      "a_0 free Continuous\n",
      "a_1 free Continuous\n",
      "a_2 free Continuous\n",
      "a_3 free Continuous\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(prob)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solve our problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print the status of our solved LP problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status: Optimal\n"
     ]
    }
   ],
   "source": [
    "print(\"Status:\", LpStatus[prob.status])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print the value of the variables at the optimum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a_0 = 0.0\n",
      "a_1 = 0.0\n",
      "a_2 = -0.19669751\n",
      "a_3 = -0.80330249\n"
     ]
    }
   ],
   "source": [
    "for v in prob.variables():\n",
    "    print(v.name, \"=\", v.varValue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print the value of the objective:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "objective= 0.5\n"
     ]
    }
   ],
   "source": [
    "print(\"objective=\", value(prob.objective))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After solving, the value of objective function is the estimation of our Kantorovich metric.\n",
    "\n",
    "As in formula (II.11) of the report, the value is close close to $1/2$ (because $d=1$)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "308.833px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}