discrete/Sinkhorn algorithm and Birkhoff-von-Neumann decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "32d2f497",
   "metadata": {},
   "source": [
    "## Optimal transport formulations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be69ae82-22c5-4701-ae08-e9ef77902a02",
   "metadata": {},
   "source": [
    "Let $x = \\sum_{i=0}^{n} \\delta_{x_i}$ and $y = \\sum_{i=0}^{n} \\delta_{y_i}$ be two discrete measures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1834cb1-2ba4-4969-8520-7097fa9ef0e0",
   "metadata": {},
   "source": [
    "**Monge formulation**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46ccddda",
   "metadata": {},
   "source": [
    "the Monge formulation of the optimal transport problem is\n",
    "\n",
    "<font color=\"blue\">\n",
    "$$ \\begin{array}{c}\\text{min}\\\\ \\sigma \\in \\Pi_{N}\\end{array}\n",
    "\\left\\{ \\sum_{i=1}^{N} C_{i},\\sigma(i) \\right\\}$$\n",
    "</font>\n",
    "\n",
    "where\n",
    "- $\\Pi_{N}$: set of all the permutation matrices of $\\{1, \\ldots, N\\} \\text{x} \\{1, \\ldots, N\\}$\n",
    "- $\\sigma \\in \\Pi_{N}$: permutation matrix\n",
    "- $\\sigma(i)$: ith row of $\\sigma$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f96208c8",
   "metadata": {},
   "source": [
    "**Kantorovich formulation**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0bb2339-c30e-4f6d-a114-a98b9e80876d",
   "metadata": {},
   "source": [
    "<font color=\"blue\">\n",
    "$$M \\in \\textbf{B}_{N} = \\left\\{ P \\in \\mathbb{R}^{N\\text{x}N}_{+} : \\sum_{i}P_{i,j}=1 \\, \\sum_{j}P_{i,j}=1 \\right\\}$$\n",
    "</font>\n",
    "\n",
    "where\n",
    "- $\\textbf{B}_{N}$: set of all the bistochastic matrices in $R^{N\\text{x}N}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a430b95-c4aa-49e8-9ce5-9f781178b9b4",
   "metadata": {},
   "source": [
    "### Sinkhorn theorem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dbd4084-2cee-4a9a-bae6-6ad2633366e9",
   "metadata": {},
   "source": [
    "There exist diagonal matrices $U$ and $V$ such that $P = U M V$ is bistochastic matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c33a96e9-d93d-4360-8f4f-64924120d796",
   "metadata": {},
   "source": [
    "### Birkhoff-von-Neumann decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b28d23a3-05ba-4d27-b7c3-f58aa9503d9b",
   "metadata": {},
   "source": [
    "The Birkhoff-von-Neumann decomposition of a bistochastic matrix $P \\in \\mathbf{B}_N$ is \n",
    "\n",
    "$$\n",
    "P = \\sum_{i=0}^{m} \\lambda_i Pi\n",
    "$$\n",
    "\n",
    "where $\\lambda_i \\in [0, 1]$ and $P_i \\in \\Pi_n$ is a permutation matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7dea933-b177-436f-98f5-d2779942d3a1",
   "metadata": {},
   "source": [
    "## From Kantorovich to Monge"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db068c3e-302f-4341-ae05-c97fa5266b08",
   "metadata": {},
   "source": [
    "**Importing libraries**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d2edb000",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import networkx as nx\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from matplotlib.gridspec import GridSpec\n",
    "plt.style.use(\"science\")\n",
    "from scipy.optimize import linprog"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df6c94a5-4e7c-409e-8170-bc044fed6aed",
   "metadata": {},
   "source": [
    "**Discrete measures**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cfbb4af3",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 6\n",
    "x = np.array([24., 9., 3., 16., 21., 27.])\n",
    "y = np.array([3., 23., 18., 26., 16., 14.])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d465acf2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASsAAAD3CAYAAABFL3JUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAS1klEQVR4nO3dMW7jSL7H8V8/zAEI9WTeqDoZYKOh5LgD0+F0RG+jD2D5BtLMCTzSDdQ+gNE28YKOHiBO0LFmOA8P2FCM1tlA4gEW0AYG+SRblCyRlFj09xNZlERVt/74s1j1V9WbxWKxEADU3H8duwEA8BIkKwBWIFkBsALJCoAVvqv6A/7+97/r3bt3O73n4eFBJycnpXx+mecq+3x1Pde+55tOp/rnP/9ZWhvWIZ7sO9e+53sWT4uK/fTTTwd5zyHOVfb56nqufc9XdhvK+oym/T/bdq59z/f0PZXfBj48POjDhw+6vb2t+qPW+vTpU23PV9dz7er29lYfPnzQw8ND5Z9FPNl3rl3lxlNZmfOl2bGq96Bcdf3e6toubGZFz2ofx8zqeNSk76BJ/xZblfEdkKywVpO+gyb9W2zV2GQFAE+RrABYgWQFwAqVF4WmU82fPn1i7KBGfv3zr73f+/OP3688vr291e3t7UFLF15jPH379q3Q+9+/f19SS6qVF0+VJ6uTkxN9/fq16o/BEaWJ48OHD5V/FvHUfHnxxG0gACtU3rMCJCmKIoVhKEmaTCa6ubmR4zjZc5Lkuq7iOFaSJHJd91hNRU3Rs8JBhGGoXq+nXq+n09NTnZ2dZc+NRiO12229efNGV1dXMsYcsaWoK5IVKhdFka6vr7PHvu8riiLFcSxJarfbms/nms/nGo/HWY8LWMZtICrnuq5ubm6yx0mSSJJarVZ2jASFbUhWOAjf97O/v3z5Is/zsgSVJImCIJD0OJ7FrSDWoc4Khe1SZ5Umpj/++CM71u12s8RljNH5+bmm0+na96fxlCKumiONoxR1VijdLnVW/X7/2bhUHMfZ7J8xRnEcK47jtb0r4qm5nl54nsZTbrJiqhllGw6H6vf7MsZk41ZxHOvs7Ezz+XzltcvjWYC0YTaQqWaUKQgCua6bJaq7uzs5jiNjjAaDQfa6MAzl+z4D7nhmbc8qnWru9XqSHgdH+/1+1jVPp5olZnGwXRzHuri4WDnmOE42VtXpdDQcDuU4jqbTqe7v74/UUtTZ2mTFVDPKZIzRYrHIfd51XYYRsFXumBVTzQDqZOtsIFPNyLNtqrkKlMI0395LxDDVjDzbppqrQDw1X14pzMZkxVQzgLrILV1gqhlAnaztWTHVDKBu1iYrppoB1A3rWQGwwlFWXSiys4r0fHcVHNchd7fB68WqCyjskLvbHKLOqsiWV7Zsd1VnR9uKCygTF7/mYysuAFYjWQGwAskKgBVIVgCsUHmySmdvln+dj2a5vb3Vhw8fKF1ApSpPVunsDct5NNenT5/09etXnZycVP5ZXPyaL+/iR+kCrELpQvNRugDAaiQrAFYgWQGwAskKgBUoXUBhlC7gEFh1AYUdctUFvF7cBsIq9NSbjzorNAI99ebbaysuoCxRFCkMQ0mPu3jf3NxkOyLFcawgCLI9KJc30QVSJCscRBiG6vV6kh73ozw7O8t2+b64uMj+juNYl5eX7JiEZxizQuWiKNL19XX22Pd9RVGU7eS9zBiT9cCAZUfZMAL7KbLRRpWbbGzbMMJ1Xd3c3GSP0929W62W7u7unu3m3Wq1FEUR271hBaULKOwlpQu+72d/f/nyRZ7nyXGcLHE9NZvNym4mLJebrBgQRRWSJFEQBNkY1abXrZP21FP02Jsj7aGnXly6wIAoqtDv9zUej7OLm+M4z3pRs9ks9+JHT725nl54XrREDAOiqMJwOFS/35cxRkmSKEkSeZ639rWdTufArUPdrU1WmwZEwzDMHRAF8gRBINd1s0R1d3cnx3FkjFl5XRzH6nQ6DCvgmdzbQAZEUZY4jnVxcbFyzHEcdbtdSdL9/b36/b5OT081mUwYUsBaW2cDqxgQ1Q/nu7UStbRtQDRljNFiscg9jzFGg8FA0upFEli2NVlVMSBapF4I9bFtQBQo08YKdgZEUTesutB8easu5CYrBkRRR2zt1nx5W7utvQ1kQBRA3axNVgyIAqgbVl0AYAWSFQArsLsNCmN3GxwCS8SgMHa3wSGwrDGsYttijt++fdv7ve/fvy+xJfbIW8yRZAWr0FNvvryeOgPsAKxAsgJgBZIVACtQuoDCKF3AIVC6gMIoXcAhcBsIq9BTb768njqlC7AKPfXmo3QBgNVIVgCsQLICYAVKF1AYpQs4BEoXUBilCzgEbgMBWIHSBVjFtiVi6qyuy9ewRAwagWGF5qPOCkcVRZHa7fba41EUSXrcAi79G3iKZIXKBUEgSWsT0Wg0Urvd1ps3b3R1dfVsE10gVflt4GseY/j1z78Kvf/nH78vqSXVyhtjSG3aW7Ldbms+n0sSu3pjI0oXUFjR0gWSFF4i9zaQMQYcQpIkCoJAQRCo3+8rjuNjNwk1tbZnFQSBjDG5YwyfP3+WJHmep/v7+2pbiEbrdrtZz8oYo/Pzc02n09zXp8MKqdc4vNBU6XBC6kWlC4wx4FDiOJbrupIek1Ucx4rjOHegnWGF5np64SmldMFxHBIVCouiSGdnZ8+Ot1qtI7QGdbfzAHs6xiBJk8mE6WbsJEmSldu+wWCQPReGoXzf50KItXZOVmWMMeiH891bitrZNsaQCsNQ4/FYknR9fa3T09MsKXU6HQ2HQzmOo+l0yhgocu2crMoYYyhaf4R62DbGkPI8T57nrfSiUq7rZvEEbLLTmBVjDACOZWvPijEG1Mlr/kVEnZW5gsNOqy4wxoC6onSh+fJ+EbE2WTHGAKBuWHUBgBXYMAKFsWEEDoFVF1AYG0bgELgNBGAFkhUAK5CsYBXGQJsvbwyU3W1gFcZAm4/dbQBYjdIFFEbpAg6B0gUURukCDoHbQABWIFkBsALJClZhDLT5KF1AIzAG2nyULgCwGqULKIzSBRwCpQsojNIFHAK3gQCsQLICYAWSFQArkKxgFSZsmo86KzQCEzbNR50VjiqKIrXb7WfH4zjWcDhUEAQaDodKkuTwjYMVqLNCYdvqrIIgkPSYsJ66uLhQr9eT7/vyfV+Xl5eVthX2os4KhW2rs/J9f+3xOI5XHhtjFIZh6e1DM+T2rOi2o2phGKrVaq0ca7Vaa3tgwNqeVRAEMsbkdtv/+OMPSY+J6/LyUvf399W2Eo2Ud6GbzWaHbQissDZZ0W3HMW3qradjoKn0FhT2u729XRnbLlS6sKnb7rpugWbiNXIc51kvajabyXGc3PcwBtpcTy88hUoX6LajTJ7nrT3e6XQO3BLYoJTZwF277frhvIyPxZFt67avkyRJ1nMyxqw8F8exOp3Oxp4VXq+dklVZ3fZf//xrl49FTW3rtqfCMNR4PJYkXV9f6/T0NBsXvb+/V7/f1+npqSaTCZM1yLVTsvI8T6PR6Nlxuu3YxPM8eZ6nwWDw7DljTHY8b2IHkF4wZrV8i0e3HcCxrE1WYRiq3+9Leuy2pz+XkP6/2x4EgUajEd12HBQ/32q+nVZdoNuOuqJ0oflYdQGA1Sr/IXPaba+y0rjI7OLPP35fYktep7SEgd1tUCVWXUBh7G6DQ+A2EIAVSFYArECyAmAFkhWsQp1V87G7DRqBCZvmO1qdFVfC5tu2YQRQBkoXUBilCzgExqwAWIFkBcAKJCsAViBZwSpM2DQfpQtoBCZsmo/SBVSG0gUcAqULKIzSBRwCY1YArECyAmAFkhUAK5CsAFiBZAWrMLvcfHmzy5QuoLBDli6ks8tVbT6C4/v06ZO+fv2qk5OTleOULqCwoqULURRJklzXVRzHSpJEruuW2UQ0wF49qyiKsgCL4zj7G9jHaDRSu93WmzdvdHV1JWPMsZuEGtqrZzUajfT582dJj7s3s4U8imi325rP55Ikx3GO2xjU1l7JiuBC2YgjbLP3mBXBhbIkSaIgCCRJk8mEW0GstVeyIrhQpm63m138jDE6Pz/XdDpd+9p0djmVDu7Dfre3tytVA6UsEVM0uPTD+T4fi5rZFlwvFcdxNvtnjFEcx4rjeO0FkNnl5np64Xk6u7xXsioaXL/++dc+H4ua2RZcLxFFkc7OzrIx0FSr1SrcPjTLzqULaXA9RXBhH8YYDQaD7HEYhvJ9nzFRPLNzz4rgQpkcx1Gn09FwOJTjOJpOp5TCYK2dkxXBhbK5rkvFOrbaa8yK4AJwaKy6AKvww/jmY9UFVIZVF1AmVl1AZdgwAofAbSAAK5CsAFiBZAXACiQrAFYgWQGwQuWzgUCZ0lKY5R9Rf/v2rdA5379/X0bTUJJ0NY9SlojZxbrgqrMiK0L8/OP3JbbEHnnBVQVKYZovrxSGOisURp0VDoExKwBWIFkBsALJCoAVSFYArECyglVYxaP5WCIGlWGJGJSJJWJQGUoXcAjcBgKwAskKgBVIVgCsQLICYAWSFQArULqAwg5ZukA8NV9ePFG6gMIOWbpAPDVfqUvExHGsIAhkjFEcx+p2u3Icp4x24hUinvASe90GXlxcqNfryfd9+b6vy8vLUhv1v//z36WeD7s75HdQdTz99ttvpZ4PuyvjO9g5WcVxvPLYGKMwDAs3ZNn/kayO7lDfwSHiiWR1fEdJVmEYqtVqrRxrtVqKoqhwY6pQdg+hzPPV9VyHZFs8lZ34yjxfXc9Vlp2TVZIka4/PZrOibalE2T2EMs9X13Mdkm3xRLI6njeLxWKxyxuGw6HG47HG43F27N27dxoMBvJ9/9nr//a3v+nf//539vjk5OTZr6mfenh42PqalyrzXGWfr67neun5Hh4eVqaXv/vuO/3rX//a6XOIp3q2rY7xtPNsoOM4z656s9ksd/Zm1+DF60I84aV2vg30PG/t8U6nU7gxeH2IJ7zUzsnKGLPyOI5jdTod6mKwF+IJL7XzmJX0GFCj0Uinp6eaTCb65ZdfCC7sjXjCS+yVrKpCJfNxpeUCrusqjmMlSSLXdY/cqv0RT8dVejwtasR13ezv6XS68H3/iK15fbrd7kLSQtLC87zFfD4/dpMKIZ6Oq+x4qvyHzC91iEpmbNZutzWfzyXJ+h4I8XR8ZcdTbdazsq2Suakcx7E+UUnEU12UGU+16VnZVsncREmSKAgCSdJkMtHV1dWz2TpbEE/HV3Y81SZZ5ckLOpRveQDaGKPz83NNp9PjNqpkxNPhlB1PtbkN3LWSGeVbHudJZ9Cejv3Ygng6vrLjqTbJikrm44qiSGdnZ8+OPx33sQXxdFxVxFNtkhWVzMdljNFgMMgeh2Eo3/et/f8nno6riniqXVEolczHE0WRwjCU4ziaTqcrwWYj4um4yo6nWiUrAMhTm9tAANiEZAXACiQrAFYgWQGwAskKOKAkSRSGIZX0e3iVySoMQ7179y773dKh3rus3W4XPgfskiSJLi8vZYxZWzBZlqbGVu1/G1gFz/NyK5yrfO+ywWBANfUrE4ahTk9PZYzZutVVkiR714Q1NbZeZc9KKra+ThmFhZ7nUaD4ysRxnH3nm777OI51d3e39+c0NbZebbIC6sr2Xw5U5VXeBqbG47Ecx1GSJM9+jhEEgRzHURzHG38qcHV1pTAMNRgM5DiOoiiSMUaTyUSDwUBhGOrq6kr9fl+SNBqNdHNzo8vLS11dXanb7Up6vEV4+t4kSXR3dydjTNZGArkeln/KM5vN9I9//EOO4+TGTRRFKxu5drvdtd95GIb6/fffsxUjPM/LVjndN7by2mRdfBVfadlOvV5vMRgMssfj8XjheV72WNJiOp0uFovHtaTv7+9X3nt/f7+Yz+fZOabT6cqa36PRKHuu2+0uut3uYrFYZOcZDAaL0Wi08b2DwWAxHo9XjuP45vP5whiTrSm+HEub4ib9TheLzfHS6/VWvusisbWpTbbF16vuWS3f13uep/Pz82xgcz6fZ1ej2Wz2bB2eyWSi6+vrbKB0NBqp1WqtrPM9mUyyz3n79q0krd0SPe+9g8FA7XZbxhh9/Pgxu1LiuNLeSBo/v/zyS/bctrhJbYqXXV67LbY2tcn3favi61Unq02ur6/19u1b+b6/dinWt2/fajAY6OLiIuveu667MlO4/OVvW8513XuTJNF8PlcURfry5cvKZ+F4ns7ULf+9LW6WbYqX5c/a9tptn5PXplarZVV8veoB9uXCvCiKslmUdHyg1+tl9/OSVq5sxphsPGE4HOrjx4/Pdk9Zfrxu7e/0vHnvvb6+VhzHcl03GxPD8fm+/2zjiTAMXxQ3qW3xsnysSGxtapNt8fWqe1bpUqtJkuj333/X/f29JGWLtKUBcXFxodFolA10pgHguq7evXunfr+vXq+nwWCgfr+v09NTSY+3lsuvTxNceiVrtVryfT8LlqfvTdcDarVams1m+vjx43H+o7DCGKPRaPTs+5KUGzfL37kxRr7vr/3OpcdJm8FgoM+fP2cXxH1ja1Msv3371qr4Yj0rAFZ41beBAOxBsgJgBZIVACuQrABYgWQFwAokKwBWIFkBsALJCoAV/gNlTOJ34AFcfQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 350x262.5 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(1, 2);\n",
    "axes[0].bar(range(N), x, color = \"skyblue\")\n",
    "axes[0].set_xlabel(\"bakeries\")\n",
    "axes[1].bar(range(N), y, color = \"silver\")\n",
    "axes[1].set_xlabel(\"cafeterias\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a171e38c-7708-46c8-8072-81be29f13ee8",
   "metadata": {},
   "source": [
    "**Cost matrix**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "80bf65d1-f586-44ec-97b7-6ca1c9f47a24",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.zeros((N, N))\n",
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        c[i,j] = np.abs(x[i] - y[j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "75fa9506-0af4-4be9-a30c-430214d563d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[21.,  1.,  6.,  2.,  8., 10.],\n",
       "       [ 6., 14.,  9., 17.,  7.,  5.],\n",
       "       [ 0., 20., 15., 23., 13., 11.],\n",
       "       [13.,  7.,  2., 10.,  0.,  2.],\n",
       "       [18.,  2.,  3.,  5.,  5.,  7.],\n",
       "       [24.,  4.,  9.,  1., 11., 13.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c0d5456",
   "metadata": {},
   "source": [
    "**Solving the problem**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5641846d",
   "metadata": {},
   "source": [
    "The optimal transport problem can be solved by the an integer linear program \n",
    "\n",
    "<font color=\"blue\">\n",
    "$$\n",
    "\\begin{align}\n",
    "\\text{min} \\big\\{ <c, z> : p \\geq 0, Az = b\\big\\}\n",
    "\\end{align}\n",
    "$$\n",
    "</font>\n",
    "\n",
    "where c is the vectorized cost matrix C, z is the vectorized matrix Z,\n",
    " $A = \\left[\n",
    "\\begin{array}{c}\n",
    "I_{N} \\otimes 1_{1\\text{x}N}\\\\\n",
    "1_{1\\text{x}N} \\otimes I_{N}\n",
    "\\end{array}\n",
    "\\right]$\n",
    "and $b=\\left[\\begin{array}{c}x\\\\ y \\end{array}\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f075fa56",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solv_lp(x, y, N, c):\n",
    "    Ax = np.kron(np.identity(N), np.ones(N))\n",
    "    Ay = np.kron(np.ones(N), np.identity(N))\n",
    "    A = np.concatenate((Ax, Ay), 0)\n",
    "    b = np.concatenate((x, y))\n",
    "    res = linprog(c, A_eq = A, b_eq = b, method = 'highs')\n",
    "    z = res.x\n",
    "    Z = z.reshape((N, N))\n",
    "    return Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e009a278",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 6)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = solv_lp(x, y, N, c.reshape(N*N))\n",
    "z.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "87483699",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0., 23.,  0.,  0.,  0.,  1.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  9.],\n",
       "       [ 3.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0., 15.,  1.],\n",
       "       [ 0.,  0., 18.,  0.,  0.,  3.],\n",
       "       [ 0.,  0.,  0., 26.,  1.,  0.]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "006aa17d-c66d-450e-bb34-1c599c2a0716",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAFjCAYAAACwvfn4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAe7UlEQVR4nO3dP3Li2Pf38c88NQtQk3amDieT7RU0DnsiPC4voPEOoLwCF94B7gW4bBFNKn1XwA9lE0pZhwPawX0CjTBgdA1YCEn9flV1tZABnSOJoz9XV/rNGGMEANjq/506AACoM4okAFhQJAHAgiIJABYUSaBEaZqeOoTWq3oe17ZIpmmq4XCoyWSiyWSix8dHPT4+njqsg0VRpMvLS3358uXUobTaKefzcDiU4zhKkkQPDw/69OmTvnz5ooeHh7UfdhiG+u2333R2dqbJZFJ5nHWVL7uzszPr+x4fH5UkSUVRSTI1NJvNTLfbNYvFYm287/um2+0e/L3j8fjNuMFgYHq93sHfuY8gCIzruu++b1ucbVZ2vrvO5zKNx2MTx/HaOM/zTL/ff/PeOI7NYDA4aix1ss9vbNdlt22+Hkst9yS/fv263Cqv6vV6cl1Xt7e3B31vEARvxl1eXur6+vqg79tXp9PZ6X3b4myzsvPddT6XJUkSzWYzua77bhxpmioMQ41Go6PFU7f1Z5/f2K7L7urqSg8PDx8Ja2e1K5LD4VCu66rb7Rb+/ZDd7aLPdLtd9Xq9g2I9hsoPJU6sDfmOx2MNh8N335emqR4fH9Xv948WSx3n5zF+Y91uV8/Pz6V+Z5HfK5nKHiaTiTzPK/x7vrWeTCYaDAYKw1DD4VCdTkdXV1eSspXx33//XW6twzBUEATLc0WSNBgMFEWRhsOhkiRRHMfL9+Yr/I8fP5QkiebzuWazmcbjsR4fH9XpdPT8/Ky7u7u1WPMfgeu6CoJAt7e31lw2FcWZx5TvRed7CqPRyDrNbbkkSbI2byQtP5+mqZIkkeM4cl333fm66eHhYbl8kiTRYDBYi2Mz/svLy6352ubPvjHlyphPtrh2mf79/b31fR+NsWj9ye27fPLvLfrcZDLR/f29kiSR7/vqdrtKkkSXl5dyXVfj8XjZtrD6G8tttjGsbjyiKHp3ObiuqyiK9vqNHaSyA/sdSTKj0cj6Htd1185x+L5vJK2dExoMBmvnLYIgMJ7nvfmu2Wz25hxIfl4kCIK1aa6eR/J9/833DQaDtRhc1107r7ptWpuK4synFwSBmc1my1jem2ZRLrPZbPm9q+ew4jhevt5lvuZ6vd7aNOI4Xjt/XBR/Ub5Fdolp23z+6Hwqspnnqm63a/r9vlksFqbX6xnHcd6cZy87xqL5eejyee9z284hbv5+ty2P0Wj05vfk+76ZzWbGcZydlsN4PH63VpShdofbkvTvv//u9X7HceR53to5obu7u4MPPTqdjpIkWTvk3zzf5Hnem+9OkkRhGK59ZvX1RziOoyiK1O125Xnecsv63jSLclmN3ff9Zeur67o6Pz9fTnOX+RpFkcIwfDON+Xy+jKUo/kPmwyHLuoz5tE2apm/WjVXz+VyPj4/yfV+dTkffv3+vPMZDl88un+t2u5rP54qiaPmezbaETfne5d3d3XLc8/PzMo80TXfKsdPp7F0rDlG7IrnLQk+SRBcXF9b3OI6zXPCHxrH5fe9dVuL7vvr9vtI0VRRFms/nms/nB01/l5h2nea2XPL35OeKPn36pLOzMz08PFgPX7bN1//7v//bGlt+2GiLvwy7LOuPzqci+emJImmaLg9Pfd/XZDIp3HAeK8ZDl8+un+v3+xqPx5KyQ/e//vrr3XjyZZbzfX85n3bNcZdaUYbaFclut2td2fO/FTXs7OoYMzeKIl1dXenl5UWu65ZSFFbj3PZjLGOaQRBoNpvp+vpa4/F471bDXS/ufW8PQzrOcpGOs2x2sTodz/PU7/cLr8441vpz6PLZ9XO3t7d6eXlZTm+XPckyzOfzSq5kqF2RHI1Gms/nhRfZDodD9fv9d0/WpmmqNE0L33foHqZtel+/ftXd3Z36/b4cx1muDB/54dviLGOa+clzz/M0GAw0m82srYbb5mt+wn7TLnv8mw5ZLu8t62MtG2n/vZnxeKz5fP6mNfyY68+hy2fXz7muq06no8lkslPR8jxva6Hct3imaVpJp4HaFUnHceT7/rLVbFX+g952PiuKorWZfH9/r36/v9war67MSZLs3SL23gLMt9ir35sfIuzzw98nzkOnuZpL3qK6GUPuvfkqZSt9t9tdO4zMp//epR+HLJddYlpVxnzaJf5NRYfBP3780MPDw9q0y4px2/w8dPns87nb21t9//59pyM813XV6/XWjljSNF3ujb6XY26XvdZSHL1p6ECLxcIMBgMzHo/X/m2Tt+j5vm+CIDCj0Whrq9dgMFh+pzFZq1uv11trUd82bjQaGcdxTLfbXbb+5e8ZDAbLFsj8+4MgMEEQmDiOTa/XW7barX7GZjPOIAhMt9s1juOY0Wj0pmV3l2lu5pLPr3y+5q2Lo9Fomc+u83U1lvz7VnO0xb8tX5v3Yiqazx+dTzabrdtxHC8/77ru2jw1xiynlc+PsmMsmp8fWT7bPrdqsVhsvephW9yr3zsajZbr4b7LodfrWa8WKEtti+Q+9r2MBLup43ytY0yj0WjtkhVUo6ruxLU73AaaZjAYLFt3UY2Hh4eDuyfviyIJlOD6+po7+lQk72X10StcdtX4Ipl3CYuiqLIO77+COs7XOsaUyxsy6tZvuo0eHx+PeoOQTb8Zw9MSAaBI4/ckAeCYKJIAYFG7W6VtSpJEk8lkeZFs3huhiaIo0vfv3zWbzU4dSinyGyBI0nQ61Y8fPxq7bPI80jTVdDrV9fX18W/BVYH8RhJNXS7S6wXs+U1lbL2rjqKSC40+YPWauPzi2ibKLwhuwCzf2erFwaPRqHbXL+7DcZzl7bjG43Hlj384hnx9q+KC62Pq9/tGkpG09bEux1brw+3NlsIybz1WtV6v14o9k1wURbq/v1++7vV6yxulNpHv+2vLp8l7XrkkSSq7kccxnZ2dabFYaLFYKAiCypdNrYtkGIZvOsx3Op3Sb06B/Xmepx8/fixf5/1rq36+TFlWr7nzfb+yC5WPZTKZ1OqxJB+1eWu1KtX6nGRRx/Yy79GIw63+CJ+fn9Xtdhu9BxZFkZ6fn3V5eXnU59AcW5qmjV4Om9I0XV6oP51OdXt7W+kecq2LZJGy7keHcuQrcdMbpPI7nufPe2/qntjLy0uji/ym1cZa13V1eXn55nk5x1Trw+1tdySez+et2kq2wXA4PMm5omNwHEdXV1e6urpq5MZ4lzuDN83qee78Kpcqz33XukgW9c3Mn8GC03t4eFg+aS+/+W3ThGGoT58+LV+vPhmwiV5eXvT4+Lh87s/9/X1jz+NHUaSvX7++GV/lue9aH25vnndIkkTn5+eN32Npyzmj/PG/eYFs6mFep9NZ2yBHUbR84FjTbO5Y3N7eVn4Or0yu66710w7DUL1er9LfT+37bidJovF4rIuLC02n08ZeGJs/E/nh4UGDwUAXFxeNPeclZctl89b5juNosVicKKKPmUwmy1M7QRBoNBo1trBIr3eczx93su8z4Osk77TgOI7iOK705hZSA4okAJxSrc9JAsCpUSQBwIIiCQAWFEkAsKBIAoAFRRIALEq5mPyPP/54c83cz58/9fnz53c/W/b7ftVpE+OvE+Mpp/0rxBjHsf7555/XEWXclPLbt287jdv1sx953686bWKs9n2/6rR/hRg3xx3tcPvm5qbU951y2vvEWHY+p4zxVLns817Ws3KwnlnsXHL3rMZlq2IaVWpTPuRSX23Kp6pcKtuTLNsx9gROqU35kEt9tSmfU+VSSt/ts7Mzff78WTc3N61aKAB+HU9PT3p6etLPnz/XbiBdSpH8888/9ffff3/0awDg5DbrWWMOtwHgFCiSAGBBkQQAC4okAFhQJAHAgiIJABYUSQCwoEgCgEUpRfLnz5/6888/9fT0VMbXAUDlnp6e9Oeff+rnz59r4+lxAwAr6HEDAHugSAKABUUSACxKecbNMSVJoslkItd1lSSJ+v2+HMc5dVi7iyIpDLPh6VT68UPK48/Hp2n2t+tryfNOEeVBGr9sVpBLjUWR9P27tHL7skoV3Z03jmMzGo2M7/tmNBqZxWKx8518y+R53lpMvV7vaNM6itFofXglH+M4xsxm2fB4bIzrVhvbBzV+2awgl5ry/ew3Us5DFHayWc8Kp7zPjD5WkYzjeC0OY4xxHOco0zqK2SwrhLk4zhZ2HGevg+D1b+PxegGtucYvmxXk0gAnLJJbz0kmSbL22nVdhfmhYYXCMFSn01kb1+l0FEVR5bEcxPOyw+tcmmb/5zl1u69/833p9ray0D6q8ctmBbnAZmuRrMuMTvOismE+n1cax4f0eq/Dz89ZYVw9PxRF0nAoXV5K/X7l4R2qFcvmP+QCm61Fsu4zuii+WktTaTLJ9hhXeZ50dyfFcfb3hmvksilALpD2bN0umtF5t8RcWQ8EcxznTWGez+fNbKkbDqUgWN+LzDmOdHWV7U0uFtvfUzNtWjbk8mvLHwCW2+yWuPVs6Hg83nryN1htaLCc6CxL0UloW0t7LY1Gr401i0X2Lwi2N+rkrd0115plY8ilEerWcNNdbVBYcX5+Xl753oHrumuvkyTR+fl5s7aKk0l2SO262SH3y0u2p9jprDfcRFE2viHXSbZi2fyHXBriRKcMth5u12lG+76v4XCoi4sLTadT+Zvn9OosSbLD6FWOkzXQeF528fjjYzY+CE53seyBGr1sNpBLTYVh9tuQpPt76eJivTG0AoV3AUqSROPxeDmj7+7uCoskdwEC0Bab9ayw4cZ1XY1GI0lSr+LKDQB1wQ0uAMCCIgkAFhRJALCgSAKABQ8CAwDxIDAA2AkPAgOAPVAkAcCCIgkAFhRJALCgSAKABUUSACwokgBgQZEEAAt63ACA6HEDADuhxw0A7IEiCQAWFEkAsKBIAoAFRRIALCiSAGBBkQQAC4okAFjQ4wYARI8bANgJPW4AYA8USQCwoEgCgAVFEgAsKJIAYEGRBAALiiQAWFAkAcCCIgkAFnRLBADRLREAdkK3RADYA0USACx+P3UA70mSRJPJRK7rKkkS9ft9OY5z6rAO0qZcpHblQy41liTSeCx9+SLFsXR3J1WZjynBt2/fyviarTzPWw7HcWx6vd7RpnVsbcrFmHblQy415rrGLBbZ8GxmTL9/1Mlt1rPCw+0oinR2dlZdtd4iSZK1167rKgzDE0XzMW3KRWpXPuRSY3ns+Z6j50mPj5WGsLVITiYTSVmhPKUwDNXpdNbGdTqdk8d1iDblIrUrH3KpsTTdPr7CfLaek+z1epUFYJMWzKD5fF5tICVoUy5Su/IhlxrzvOycZC4vjhXm08jW7aIVoYnalIvUrnzIpQZcVxqNskPsNH0tmBt7y8dUSut23uMmd3Nzo5ubmw9/r+M4b7aA8/m8kS11bcpFalc+5FJzg0FWHJNE6nazca5b2tc/PT2t9Rbc7HFjbd1+58+FrUFlieN4raXOGGMcxzGLvKWrQdqUizHtyodcai6OX4dnM2M28ivbzq3bdeBubC2SJNH5+Xkjt4ptykVqVz7kUnNnZ68NOONxdvhdodpfTO77vobDoS4uLjSdTuX7/qlDOlibcpHalQ+51NholF0KNJ9LV1evh9wVsd7g4rffftNisXh3K8QNLgC0xU43uAjDUMPhUJJ0f3+/vG4SAH41Ww+3u92uut2uRhUf+wNA3dS64QYATo0iCQAWFEkAsOAZNwAgnnEDADvhGTcAsAeKJABYUCQBwIIiCQAWFEkAsKBIAoAFRRIALCiSAGBBkQQAC7olAoDolggAO6FbIgDsgSIJABYUSQCwoEgCgAVFEgAsKJIAYEGRBAALiiQAWNDjBgBEjxsA2Ak9bgBgDxRJALCgSAKABUUSACwokgBgQZEEAAuKJABYUCQBwIIeNwAgetwAwE7ocQMAe6BIAoAFRRIALH4/dQDvCsPs/zSVplPp+lryvJOGdKgkSTSZTOS6rpIkUb/fl+M4pw7rYG3Kh1zq6+T5mAKz2cyMRiMzGo1Mr9czi8Wi6K3m27dvhX/7MMcxZjbLhsdjY1z3eNM6Ms/zlsNxHJter3fCaD6uTfmQS31Vnc9mPSvckwzDUIPBQJL08PCgr1+/ajabVVa8l3x/fc+xoVvEJEnWXruuqzDfS26gNuVDLvVVh3y2npOMokj39/fL171eT1EUvQm4Et3u67DvS7e31cdQgjAM1el01sZ1Oh1FUXSiiD6mTfmQS33VIZ+te5Ke5+nHjx/L12maStKbYCsTRdLzs3R5KfX7p4nhg/J5uGk+n1cbSEnalA+51Fcd8ik83O71esvh5+dndbvdwpOleY+b3M3NjW5ubsqL0vMk15WGQ2kykVZia7qilaCp2pQPudRXmfk8PT2t9Rbc7HHzbut2mqaaTCbW85GfP38+fo8bx5GurrK9ycWicecmHcd5s/Wbz+eNbXVsUz7kUl9V5LO5U7e6wyftcJ3kcDhUEASnmclhKH369PradbP/T3Fu9IO6q+dWV5yfn1ccSTnalA+51Fcd8rEWyYeHBw2HQ7muqzRNq99l73TWG26iKNuDbOB1km5e4P+TJInOz88bu4VvUz7kUl91yKfwcHsymcjzvGWBfHl5Ub/qRhPPyy4ef3zMXgeBdIrLkEri+76Gw6EuLi40nU7l+/6pQ/qQNuVDLvV16ny23gUoSRJ9+fJlbZzjOFosFlu/hLsAAWiLzXq2dU/SdV1tqZ0A8MvhBhcAYEGRBAALiiQAWFAkAcCCB4EBgHgQGADshAeBAcAeKJIAYEGRBAALiiQAWFAkAcCCIgkAFhRJALCgSAKABT1uAED0uAGAndDjBgD2QJEEAAuKJABYUCQBwIIiCQAWFEkAsKBIAoAFRRIALOhxAwCixw0A7IQeNwCwB4okAFhQJAHAgiIJABYUSQCwoEgCgAVFEgAsKJIAYEGPGwAQPW4AYCf0uAGAPVAkAcCCIgkAFr+fOoD3JEmiyWQi13WVJIn6/b4cxzl1WAdpUy5Su/JpfC5RJH3/Ls1mb8dLkudJSSKlaTbcJEW5VcUUCILABEFgfN83g8HAzGazoreab9++Ff7tozzPWw7HcWx6vd7RpnVsbcrFmHbl0+hcfN+Y2cyYbT/nfj8bLxnT7RqzWFQe3ofYcjuSzXpWuCd5dXWl//3vf+p2u5rP57q6ulIcx1XWbyVJsvbadV2FYVhpDGVpUy5Su/JpfC69XvHfzs6kxSIbbtKecc6WW0UKz0n6vi9vZbf8FIceYRiq0+msjet0OoryQ4gGaVMuUrvyaVMuWzlOMwtkTRTuSXa73eWw7/u6vb2tJKBVaZpuHT+fz6sNpARtykVqVz5tyuWNNJUmk2x4OpVubyXXPWlITWNtuImiSM/Pz7q8vFS/368qpncVrdRN1KZcpHbl04pc+v3XvUjXlS4vpYpPmzWdtUh6nifXdTUcDjWZTNQrOD+Qd0vM3dzc6Obm5sPBOY7zZms+n8+b1er4nzblIrUrnzbl8kaSvLZmu272OknYm1zx9PS01qV6s1viTk1GQRAYSWZR0DJ2rNbtOI7XWh2NMcZxnMI46qxNuRjTrnxak8vmz3k2M8ZxXl8vFtl7mpaXMSdt3d7acBOGoT59+rR87f631dlsBTw2d2NrlySJzs/PG7mFb1MuUrvyaVMuWj1F4LrSaPT6Ogyz1uIm5iWt51ahrYfbnU5nreEmiiI5jrPW2l0V3/c1HA51cXGh6XQq3/crj6EsbcpFalc+jc4lDKUgyIbv76WLi9dieH4uPTxkw3EsNSkvqTi3ChXeBWgymSzP0wRBoNFo9GaLm+MuQADaYrOeFTbcrDbS1KllGwCqxA0uAMCCIgkAFhRJALCgSAKABQ8CAwDxIDAA2AkPAgOAPVAkAcCCIgkAFhRJALCgSAKABUUSACwokgBgQZEEAAt63ACA6HEDADuhxw0A7IEiCQAWFEkAsKBIAoAFRRIALCiSAGBBkQQAC4okAFjQ4wYARI8bANgJPW4AYA8USQCwoEgCgAVFEgAsKJIAYEGRBAALiiQAWFAkAcCCIgkAFnRLBADRLREAdkK3RADYA0USACx+P3UA70mSRJPJRK7rKkkS9ft9OY5z6rAO0vhcokj6/l2azdbHJ4kUhlKnkw33epLrnibGAzV+2axoUy6SsnVLktJUmk6l62vJ86qbvinBt2/fyviarTzPWw7HcWx6vd7RpnVsjc7F942ZzYzZtsqMRuuv+/1qYipRo5fNhjblYowxxnGydc8YY8ZjY1z3qJPbrGc7HW4Ph0OlaXrkcv1WkiRrr13XVZhvVRqm8bn0esVb7+fnamMpWeOXzYo25bLk++vrXsV7xe8WySiK9PDwUEUsb4RhqE6nszau0+koiqKTxPMRbcrljU5HOjt7Pey+vDx1RHtp07JpUy5L3e7rsO9Lt7eVTv7dIpkkidwTnV8q2nudz+fVBlKCNuXyhu9n/3/5kg33eqeNZ09tWjZtymVNFEnDYbYB7vcrnbS1SE4mE/VquMKf4tD/WFqRSxhKo5E0HkuPj5Vv6Y+lFcvmP43PxfOkuzspjqXJpNJJF7Zup2m6c4tY3uMmd3Nzo5ubmw8H5zjOmy3gfD5vZEtdm3JZkyRZi+NolL3udrND7+GwMS3cbVo2bcrlDceRrq6yvcnForRzk09PT2u9BTd73BS2bo/H4+Ww67pmsVjs3BpUljiO11rqjDHGcRxrLHXVmlw2Vxnfz/6tGo1eWyMboDXLxrQrF2OMMUGQtW7n4jhbB4+4fu3Uuh2Gof76669SqvRHbJ4LTZJE5+fnjdwqtikXrR66eV62J7nq33+rvY7tg9q0bNqUi6SsUXC14SaKsj3ICtevwsPtl5eX5XCSJLq/v9f19bW8ild+3/c1HA51cXGh6XQqP28kaKBG5xKGUhBkw/f30sXF60Xjl5fSw8Pr4U8Dz0k2etlsaFMu8rzs4vHHx+x1ELztzHBkO93g4rffflMcx4Wt3NzgAkBb7HWDizRNl9dIjkajZl9rBQAHsPbddhxHg8FAg8GgqngAoFa4CxAAWFAkAcCCIgkAFjzjBgDEM24AYCc84wYA9kCRBAALiiQAWFAkAcCCIgkAFhRJALCgSAKABUUSACwokgBgQbdEABDdEgFgJ3RLBIA9UCQBwIIiCQAWFEkAsKBIAoAFRRIALCiSAGBBkQQAC3rcAIDocQMAO6HHDQDsgSIJABYUSQCwoEgCgAVFEgAsKJIAYEGRBAALiiQAWNDjBgBEjxsA2Ak9bgBgDxRJALCgSAKAxe+nDuA9SZJoMpnIdV0lSaJ+vy/HcU4d1kHalIvUgnyiSArDbHg6lX78kFbjD0MpSSTXzV53u5WHeJAokr5/l2azU0dSipOvZ6bAbDYzs9nMGGNMHMfL4W2+fftW+LeP8jxvORzHsen1ekeb1rG1KRdjWpDPaLQ+vJKPCQJj+v1sOI6Ncd1qYzuU7xszmxlT/NNunKrXs816Vjgn+/2+kWQkmW63axaLxc5fWpY4jtdmkDHGOI5zlGkdW5tyMaYF+cxmxqzGG8dZYYnj7LXrGrO6zufjm6IlRfIU69lmPSs8J3l2dqbFYqHFYqEgCE5yGBWGoTqdztq4TqejKIoqj+Wj2pSL1IJ8PC87vM6lafZ/p5MdYs/n2aF3FGV/yw+5Uak6rGfWhhvHcU56jinNV9wN8/m82kBK0KZcpJbk0+u9Dj8/Z+cc88LY6UiTSVYcHx+zYVSuDutZYcNNmqaa/LdiTKdT3d7eyi3YmuY9bnI3Nze6ubkpOdT12NqiTblIDc0nTbMimDd0zOfZ3mReNPt96dMn6eP9LlCSMtezp6entd6Cmz1uCovkaguS67q6vLxUHMdb3/v58+ej9LhxHOfNFmM+nzerBfU/bcpFalk+w6EUBK8t266bDeev8/+jKDtMR2WqWM82d+pWd/gky+F2kiTL4bzpfXVcFboFl1ycn59XGkcZ2pSL1KJ8Hh6yIum62R4l5x9rpQ7r2dYiGUWRvn79+mb85gnUY9s8vE+SROfn543cW2lTLlJL8plMsj3DvEC+vGR7ja4rnZ+/Nubk10o2bS+yiac+NtRhPdt6uO26rkaj0fJ1GIbq9Xon+QH4vq/hcKiLiwtNp1P5vl95DGVpUy5Sw/NJEunqan1cfv5Rknw/28M8O8vOVQZB5SEeJAxfY72/ly4u1huoGujU61nhXYCiKFIYhnIcR3EcrxXNTdwFCEBbbNazwoYbz/PkNe3wAgBKxg0uAMCCIgkAFhRJALCgSAKABQ8CAwDxIDAA2AkPAgOAPVAkAcCCIgkAFhRJALCgSAKABUUSACwokgBgQZEEAAt63ACAinvclFIk8weBHfMJiW0rwG3Kh1zqq035HDuXm5sb/f333/r8+fPa+MYcbrdpYUvtyodc6qtN+Zwql6MVyV0TOkbiZU97nxjLzueUMZ4ql33ey3pWDtazYhTJEt+373vL/L5fceU9xvv2wXpW/XeW/X27vLeUuwD98ccf+vLly9q4nz9/vjm236bs9/2q0ybGXyfGU077V4gxjmP9888/y9elFEkAaKvGNNwAwClQJAHAgiIJABb/H2KSon5s2S/EAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (4, 4))\n",
    "ax.matshow(c, cmap = mpl.colors.ListedColormap([\"white\"]))\n",
    "ax.set_title(\"Optimal transport plan (Kantorovich)\")\n",
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        color = 'red' if z[i,j] > 0 else 'black'\n",
    "        ax.text(j, i, int(z[i, j]), horizontalalignment = 'center',\n",
    "                verticalalignment = 'center', weight = 'bold', color = color)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3139b825-f833-4101-8a68-dacc2ba641d0",
   "metadata": {},
   "source": [
    "**Sinkhorn algorithm**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "636512c7-7f49-413f-96b9-c983502dd3da",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sinkhorn_algorithm(z, N, itermax = 1000, eps = 1e-1):\n",
    "    u = np.ones(N)\n",
    "    u0 = u.copy()\n",
    "    v = np.ones(N)\n",
    "    v0 = v.copy()\n",
    "    for i in range(itermax):\n",
    "        v = np.ones(N) / (z.T @ u)\n",
    "        u = np.ones(N) / (z @ v)\n",
    "        if np.linalg.norm(u - u0) < eps and np.linalg.norm(v - v0) < eps:\n",
    "            print(\"Sinkhorn algorithm converged after {} iterations\".format(i))\n",
    "            break\n",
    "        else:\n",
    "            u0 = u.copy()\n",
    "            v0 = v.copy()\n",
    "    return u,v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "538169ac-6a36-4f5e-934b-89f1398cb4ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "u, v = sinkhorn_algorithm(z, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "940d332f-6d83-4880-8e04-c8958c2d5541",
   "metadata": {},
   "source": [
    "The matrix $P = Diag(u) M Diag(v)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f4827d42-0d50-45ad-98db-5b5505c6a93b",
   "metadata": {},
   "outputs": [],
   "source": [
    "P = np.diag(u) @ z @ np.diag(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08ba8f20-470b-4c75-bc29-bd62259805ab",
   "metadata": {},
   "source": [
    "We can check that P is indeed a bistochastic matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "220b5f12-1cf2-485a-bc25-2b0455cb547a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 0.99988094, 0.9998426 , 0.99915164, 1.00010294,\n",
       "       1.00102188])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.sum(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c2426291-e959-490b-a66a-867fbdec0e6c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1.])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.sum(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a754c46-1578-464a-8ee1-6981fab10e89",
   "metadata": {},
   "source": [
    "**Solution to the Monge formulation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "924dc1d5-c83c-4a9e-84e4-ebc3cddb39cd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUsAAAFRCAYAAAD9+o9jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAPuElEQVR4nO3dz2sc990H8I8eCsHQpOOF58E/iA/jU3syaxkKPVpqLoU6sIqbU09Z/Qde+hcI6T+wfcoptbSQprdUew0UhAZfetQcUuz6UKwtOZRAYZ9D2GU32rW/7s7+1Ot10kwmM5/RrN878535aDZ6vV4vAHij/1l0AQCrQFgCJBCWAAmEJUCCnyy6gLcpyzLa7XbkeR5lWUaz2YwsyxZd1n+tKIr47LPP4vT0dNGlTK0oiuh0OhERcXJyEk+fPl3pY9Pfl263GycnJ/Hw4cOo1+sLrqoarVYr/vCHP6zs8SmKIiIi6vV6lGUZ3W53/semt+Tq9frg57Ozs16j0VhgNdM5OjrqnZ6e9lbg155kf39/5OfhY7WKsizrnZ6e9nq9Xu/x48e9PM8XXFE1+p+58/PzRZfyX2s2m72I6EVEb2trayH7stSX4WVZjkzneT749l9FjUZjbc5UiqKIvb29wXSj0YiiKC4cs1VydHQ0cnxW9Szsx8qyjDzPF13GVO7evRvn5+dxfn4ex8fHCzk2Sx2WnU4narXayLxarTY4JWdx6vV6PH36dDDd7XYjIi4cr1WytbU1+Pno6Ch2d3cXWE012u12NBqNRZdRiSzLFvoFttRjlv1/gD/2+vXr+RbCWMP/CJ89exZbW1srfzZWFEU8e/Ystre3o9lsLrqcqXS73ZU/Hn3dbjfa7XZE/DA+vru7O/ez5aUOy0kmhSiL0f8gr8NNq3q9HnmeR6vVWvmzssPDw5UP/L7hG7t5nsf29nacnZ3NtYalvgzPsuzCWeTr16/X5ttyXbRarYWNI81ClmWxs7MTOzs7K/vF3Ol04pNPPll0GZUZHgvvPxkz7/HxpQ7L4TGkYZubm3OuhEkODg6i1WpFnufR7XZXOlyuXr06mO5f4q3yDavDw8N48uRJPHnyJMqyjL29vZUc7y+KIu7fv39h/rzHx5f6MvzHYxJlWcbm5uZanMGsw3hSu90eXLZ2u92Vvuyr1WojX85FUUSWZSv79MKPTzR2d3cXMs5XhTzPY39/fzDd6XSi0WjM/d/PRq+33H91qCzLePz4cdy7dy9OTk5W+sHaTqcTx8fHcXBwEI8ePYp79+6t7JhYWZZx+/btkXlZlsX5+fmCKppeu90eDPscHx/H/v7+SobLsG63G0+ePIlWqxXNZjN2d3dX8gug3wCRZVmcnZ2NhOe8LH1YAiyDpR6zBFgWwhIggbAESCAsARIIS4AEwhIgQSUPpf9044O4Ej8dmfd9/Dveiytv/X9Tl5vFOict98tf/+vCvBf/+E/cvP72X1fVy027zr/+5WcXllvG37lt2/aybfv/fnE1/va3vw2mKwnLK/HTuLPxq5F5z3vfXJg3Tupys1jnpOW++vz5hXm//f3L+OrzG29dZ9XLTbvOj27cubDcMv7Obdu2l23bV25/PzI9s8vwa/Fhpcstetu/e/B+pcu9i6q3/S77XfXxmcW21+2ztk77s06ftdmF5catSpdb9LY//TgtiFKXexdVb/td9rvq4zOLba/bZ22d9medPmsrdYNnFt+kizKLM9BFWqdjE7Fe+7NO+xKxuP2ppDf8g42r8V5ciWvx4Uy+7ebt65fPF11CZcaNWQKTvep9G6/i73Grfn3kD1pXcoPnvbiSPJAKsMyubdyKa3Errtyc0w0egHUiLAESCEuABMISIIGwBEggLAESCEuABMISIIGwBEhQSQfP9/HveN77Zm3aHYHLa9Du+OL6yHztjgBDtDsCTEFYAiQQlgAJJo5ZlmUZ7XY78jyPsiyj2WxGlmVzLA1geUwMy52dncEfvizLMj777LM4OjqaW2EAy2TsZXhZliPTeZ5Hp9OZS0EAy2hsWHY6najVaiPzarVaFEUxl6IAls3YsOx2u2MXfv369SxrAVha7/RQ+qQQ7Xfw9OnkAVZNv3OnL6mDJ8uyC2eRr1+/nng3XAcPsOr6nTt9SR08W1tbY1e2ublZYWkAq2NsWOZ5PjJdlmVsbm56zhK4tCaOWR4dHUWr1Yp79+7FycmJZyyBS21iWOZ5Hvv7+xER0Wg05lYQwDLSGw6QQFgCJBCWAAmEJUAC7+ABGOIdPAAJvIMHYArCEiCBsARIICwBEghLgATCEiCBsARIICwBEghLgATaHQGGaHcESKDdEWAKwhIggbAESFDJmOW6+ejGnUWXACvn65fPF11CpR40fz4y7cwSIIGwBEggLAESCEuABJWEZb+D51Xv2ypWB7AwX3z5Xfz29y/jxYsXI/MrCct+B49WR2DVffrx+/HV5zfi5s2bI/NdhgMkEJYACYQlQAJhCZBAWAIkEJYACYQlQAJhCZBABw/AEB08AAl08ABMQVgCJBCWAAkmhmVRFHH37t151gKwtMaGZbvdjogfAhOACW93bDQa864DYKkZswRIICwBEghLgARjxyzfVb/dse9afKibB1gpX3z5XfzxT98Npl/+84OR/15JWPbbHQFW1acfvx+ffvz+YPpB8x3bHbvdbuVFAayasWHZ6XSi1WpFRMTe3t7guUuAy2qj1+v1pl3J/27ccBkOl9zXL58vuoRKPWj+PP785z8Ppt0NB0ggLAESCEuABMISIIGwBEjghWUAQ7ywDCCBF5YBTEFYAiQQlgAJhCVAAmEJkEBYAiQQlgAJhCVAAh08AEN08AAk0MEDMAVhCZBAWAIkqOS94b/89b/iq8+fV7GqpfDRjTuLLgFWzrr9u7nym+9Hpp1ZAiQQlgAJhCVAAmEJkEBYAiSoJCxf/OM/8dvfv4wvvvyuitUBLMyr3rfxvPfNhXbHSh4dunn9J/HV5zeqWBXAQl3buBXX4lZcuenRIYB3JiwBEghLgATCEiCBsARIICwBEghLgATCEiCBDh6AITp4ABLo4AGYgrAESCAsARJMHLMsiiI6nU5ERJycnMTTp08jy7J51QWwVCaGZafTiUePHkVExMHBQdy/fz9OT0/nVhjAMhl7GV4URezt7Q2mG41GFEURZVnOrTCAZTI2LOv1ejx9+nQw3e12IyKiVqvNpSiAZTPxBk+j0Rj8/OzZs9ja2jJmCVxab30ovdvtRrvdfuN4Zb+Dp+93D96PTz9+v5oKAebgVe/beBV/H0zfenF95L+/NSxbrVYcHx+/8axSBw+w6vqdO33v1MFzcHAQrVYr8jyPbrc7GLsEuGwmhmW73Y56vT4IysPDQ2OWwKU19jK8LMvY2dkZmZdlWTSbzbkUBbBsxoZlnufR6/XmXQvA0tIbDpBAWAIkEJYACYQlQAJhCZDAC8sAhnhhGUACLywDmIKwBEggLAESCEuABMISIIGwBEggLAESCEuABDp4AIbo4AFIoIMHYArCEiCBsARIUMmY5V//8rP46MadKlYFl8rXL58vuoTKrHsGOLMESCAsARIIS4AEwhIgQSU3eL6Pf8fz3jdxLT6Maxu3qlglwEK86n0br+LvcevF9ZH5lYTle3El7mz8qopVASyUDh6AKQhLgATCEiCBsARIICwBEghLgATCEiCBsARIICwBEmh3BBii3REggXZHgCkIS4AEwhIgwcQxy06nExER3W43Tk5O4uHDh1Gv1+dWGMAymXhmubOzE7VaLRqNRty+fTt2dnbmWRfAUpkYlkdHRyNnklmWzaMegKU0MSy3trYGPx8dHcXu7u5cCgJYRm98zrIoinj27Flsb29Hs9mcV00AS+eNYVmv1yPP82i1WtFut6PRaIxdrt/B06eTB1g1/c6dvnfu4MmyLHZ2dmJ7ezvOz8/Hjl3q4AFWXb9zpy+pg6fT6cTVq1cH03meR0REWZazqBFg6Y0Ny1qtNnKDpyiKyLLMc5bApTX2Mrxer8fDhw/jyZMnERFxfHwcp6ency0MYJlMHLMcvpnjTjhw2ekNB0ggLAESCEuABMISIIF38AAM8Q4egATewQMwBWEJkEBYAiQQlgAJhCVAAmEJkEBYAiQQlgAJdPAADNHBA5BABw/AFIQlQAJhCZBAWAIkqOQGD8vr65fPF11CpT66cWfRJVRq3fZnnTmzBEggLAESCEuABMISIIF2R4Ah2h0BEmh3BJiCsARIICwBEghLgATCEiCBsARIICwBEghLgAQ6eACG6OABSKCDB2AKwhIggbAESJAUlq1WK7rd7oxLAVhebw3Loiji4OBgHrUALK23hmVZlpHn+TxqAVhabwzLdrsdjUZjXrUALK2JYdntdiPLsjmWArC8Jobl4eFhbG1tzbMWgKU1toOn0+nEJ598krySfrtjn7ZHYNX02xz7ktsdDw8PBz+XZRl7e3vx8OHDqNfrF5bV7gisun6bY9+P2x3HhuWPL793d3djd3fXXXHg0nrj3fButzt4xnJ/fz+KophLUQDL5o1/dSjLsnj06FE8evRoXvUALCW94QAJhCVAAmEJkEBYAiQQlgAJvLAMYIgXlgEk8MIygCkIS4AEwhIggbAESCAsARIIS4AEwhIggbAESKCDB2CIDh6ABDp4AKYgLAESCEuABJWMWbK8PrpxZ9ElVOrrl88XXUKl1u34rDNnlgAJhCVAAmEJkEBYAiQQlgAJtDsCDNHuCJBAuyPAFIQlQAJhCZBAWAIkEJYACYQlQAJhCZBAWAIk0MEDMEQHD0ACHTwAUxCWAAmEJUCCiWFZFEUURREREWVZDn4GuIwmhuXjx4/j7t27sbGxEbu7u5Hn+TzrAlgqE++G3717N87PzyMiIsuyedUDsJTe+OiQkAT4wcSw7Ha70W63IyLi5OTEpThwqU0My2azOTizzPM8tre34+zsbOyy/Q6ePp08wKrpd+70JXfwlGUZ9Xo9In4Iy7IsoyzLsWeXOniAVdfv3OlL6uApiiLu379/YX6tVqu4PIDVMDYs8zyP/f39wXSn04lGo+GGD3Bpjb0Mz7IsNjc34+DgILIsi7Ozszg6Opp3bQBLY+KYZb1eH4xZAlx2esMBEghLgATCEiCBsARIICwBEnhhGcAQLywDSOCFZQBTEJYACYQlQAJhCZBAWAIkEJYACYQlQAJhCZBABw/AkEkdPJWcWfY7eGYdlK963850/fO0TvsSsX7788WX3y26hMqs27GZ9f5c27gVdzZ+FTdv3hyZv1KX4cOvqVx167QvEeu3P3/80xqF5Zodm0Xtz8zCMjX9Z/EtMYttr9P+zGK/Z7G+1GVTzwJncbY4i237rFWj6m3PLiwT038W3xKz2PY67c8s9nsW60tdNvUscBZni7PYts9aNarediU3eP7vF1fjyu3Rv9Bx68X1C3+1Y5zU5WaxTttevW0/aP78wnIv//lBPGjevDD/v11uFuuctNyV31z8XSzb7/yybvvs7GxkeqPX6/WS1gZwia3UDR6ARRGWAAmEJUCC/weZrbwWe68c8wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.matshow(P)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a8541e8b-a4f1-4b70-b2f6-8439b244b9a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = np.zeros((N, N))\n",
    "Q[P > 0.1] = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "930aac0e-1cfb-4719-b643-b7b8f931352b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUsAAAFRCAYAAAD9+o9jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAPuElEQVR4nO3dz2sc990H8I8eCsHQpOOF58E/iA/jU3syaxkKPVpqLoU6sIqbU09Z/Qde+hcI6T+wfcoptbSQprdUew0UhAZfetQcUuz6UKwtOZRAYZ9D2GU32rW/7s7+1Ot10kwmM5/RrN878535aDZ6vV4vAHij/1l0AQCrQFgCJBCWAAmEJUCCnyy6gLcpyzLa7XbkeR5lWUaz2YwsyxZd1n+tKIr47LPP4vT0dNGlTK0oiuh0OhERcXJyEk+fPl3pY9Pfl263GycnJ/Hw4cOo1+sLrqoarVYr/vCHP6zs8SmKIiIi6vV6lGUZ3W53/semt+Tq9frg57Ozs16j0VhgNdM5OjrqnZ6e9lbg155kf39/5OfhY7WKsizrnZ6e9nq9Xu/x48e9PM8XXFE1+p+58/PzRZfyX2s2m72I6EVEb2trayH7stSX4WVZjkzneT749l9FjUZjbc5UiqKIvb29wXSj0YiiKC4cs1VydHQ0cnxW9Szsx8qyjDzPF13GVO7evRvn5+dxfn4ex8fHCzk2Sx2WnU4narXayLxarTY4JWdx6vV6PH36dDDd7XYjIi4cr1WytbU1+Pno6Ch2d3cXWE012u12NBqNRZdRiSzLFvoFttRjlv1/gD/2+vXr+RbCWMP/CJ89exZbW1srfzZWFEU8e/Ystre3o9lsLrqcqXS73ZU/Hn3dbjfa7XZE/DA+vru7O/ez5aUOy0kmhSiL0f8gr8NNq3q9HnmeR6vVWvmzssPDw5UP/L7hG7t5nsf29nacnZ3NtYalvgzPsuzCWeTr16/X5ttyXbRarYWNI81ClmWxs7MTOzs7K/vF3Ol04pNPPll0GZUZHgvvPxkz7/HxpQ7L4TGkYZubm3OuhEkODg6i1WpFnufR7XZXOlyuXr06mO5f4q3yDavDw8N48uRJPHnyJMqyjL29vZUc7y+KIu7fv39h/rzHx5f6MvzHYxJlWcbm5uZanMGsw3hSu90eXLZ2u92Vvuyr1WojX85FUUSWZSv79MKPTzR2d3cXMs5XhTzPY39/fzDd6XSi0WjM/d/PRq+33H91qCzLePz4cdy7dy9OTk5W+sHaTqcTx8fHcXBwEI8ePYp79+6t7JhYWZZx+/btkXlZlsX5+fmCKppeu90eDPscHx/H/v7+SobLsG63G0+ePIlWqxXNZjN2d3dX8gug3wCRZVmcnZ2NhOe8LH1YAiyDpR6zBFgWwhIggbAESCAsARIIS4AEwhIgQSUPpf9044O4Ej8dmfd9/Dveiytv/X9Tl5vFOict98tf/+vCvBf/+E/cvP72X1fVy027zr/+5WcXllvG37lt2/aybfv/fnE1/va3vw2mKwnLK/HTuLPxq5F5z3vfXJg3Tupys1jnpOW++vz5hXm//f3L+OrzG29dZ9XLTbvOj27cubDcMv7Obdu2l23bV25/PzI9s8vwa/Fhpcstetu/e/B+pcu9i6q3/S77XfXxmcW21+2ztk77s06ftdmF5catSpdb9LY//TgtiFKXexdVb/td9rvq4zOLba/bZ22d9medPmsrdYNnFt+kizKLM9BFWqdjE7Fe+7NO+xKxuP2ppDf8g42r8V5ciWvx4Uy+7ebt65fPF11CZcaNWQKTvep9G6/i73Grfn3kD1pXcoPnvbiSPJAKsMyubdyKa3Errtyc0w0egHUiLAESCEuABMISIIGwBEggLAESCEuABMISIIGwBEhQSQfP9/HveN77Zm3aHYHLa9Du+OL6yHztjgBDtDsCTEFYAiQQlgAJJo5ZlmUZ7XY78jyPsiyj2WxGlmVzLA1geUwMy52dncEfvizLMj777LM4OjqaW2EAy2TsZXhZliPTeZ5Hp9OZS0EAy2hsWHY6najVaiPzarVaFEUxl6IAls3YsOx2u2MXfv369SxrAVha7/RQ+qQQ7Xfw9OnkAVZNv3OnL6mDJ8uyC2eRr1+/nng3XAcPsOr6nTt9SR08W1tbY1e2ublZYWkAq2NsWOZ5PjJdlmVsbm56zhK4tCaOWR4dHUWr1Yp79+7FycmJZyyBS21iWOZ5Hvv7+xER0Wg05lYQwDLSGw6QQFgCJBCWAAmEJUAC7+ABGOIdPAAJvIMHYArCEiCBsARIICwBEghLgATCEiCBsARIICwBEghLgATaHQGGaHcESKDdEWAKwhIggbAESFDJmOW6+ejGnUWXACvn65fPF11CpR40fz4y7cwSIIGwBEggLAESCEuABJWEZb+D51Xv2ypWB7AwX3z5Xfz29y/jxYsXI/MrCct+B49WR2DVffrx+/HV5zfi5s2bI/NdhgMkEJYACYQlQAJhCZBAWAIkEJYACYQlQAJhCZBABw/AEB08AAl08ABMQVgCJBCWAAkmhmVRFHH37t151gKwtMaGZbvdjogfAhOACW93bDQa864DYKkZswRIICwBEghLgARjxyzfVb/dse9afKibB1gpX3z5XfzxT98Npl/+84OR/15JWPbbHQFW1acfvx+ffvz+YPpB8x3bHbvdbuVFAayasWHZ6XSi1WpFRMTe3t7guUuAy2qj1+v1pl3J/27ccBkOl9zXL58vuoRKPWj+PP785z8Ppt0NB0ggLAESCEuABMISIIGwBEjghWUAQ7ywDCCBF5YBTEFYAiQQlgAJhCVAAmEJkEBYAiQQlgAJhCVAAh08AEN08AAk0MEDMAVhCZBAWAIkqOS94b/89b/iq8+fV7GqpfDRjTuLLgFWzrr9u7nym+9Hpp1ZAiQQlgAJhCVAAmEJkEBYAiSoJCxf/OM/8dvfv4wvvvyuitUBLMyr3rfxvPfNhXbHSh4dunn9J/HV5zeqWBXAQl3buBXX4lZcuenRIYB3JiwBEghLgATCEiCBsARIICwBEghLgATCEiCBDh6AITp4ABLo4AGYgrAESCAsARJMHLMsiiI6nU5ERJycnMTTp08jy7J51QWwVCaGZafTiUePHkVExMHBQdy/fz9OT0/nVhjAMhl7GV4URezt7Q2mG41GFEURZVnOrTCAZTI2LOv1ejx9+nQw3e12IyKiVqvNpSiAZTPxBk+j0Rj8/OzZs9ja2jJmCVxab30ovdvtRrvdfuN4Zb+Dp+93D96PTz9+v5oKAebgVe/beBV/H0zfenF95L+/NSxbrVYcHx+/8axSBw+w6vqdO33v1MFzcHAQrVYr8jyPbrc7GLsEuGwmhmW73Y56vT4IysPDQ2OWwKU19jK8LMvY2dkZmZdlWTSbzbkUBbBsxoZlnufR6/XmXQvA0tIbDpBAWAIkEJYACYQlQAJhCZDAC8sAhnhhGUACLywDmIKwBEggLAESCEuABMISIIGwBEggLAESCEuABDp4AIbo4AFIoIMHYArCEiCBsARIUMmY5V//8rP46MadKlYFl8rXL58vuoTKrHsGOLMESCAsARIIS4AEwhIgQSU3eL6Pf8fz3jdxLT6Maxu3qlglwEK86n0br+LvcevF9ZH5lYTle3El7mz8qopVASyUDh6AKQhLgATCEiCBsARIICwBEghLgATCEiCBsARIICwBEmh3BBii3REggXZHgCkIS4AEwhIgwcQxy06nExER3W43Tk5O4uHDh1Gv1+dWGMAymXhmubOzE7VaLRqNRty+fTt2dnbmWRfAUpkYlkdHRyNnklmWzaMegKU0MSy3trYGPx8dHcXu7u5cCgJYRm98zrIoinj27Flsb29Hs9mcV00AS+eNYVmv1yPP82i1WtFut6PRaIxdrt/B06eTB1g1/c6dvnfu4MmyLHZ2dmJ7ezvOz8/Hjl3q4AFWXb9zpy+pg6fT6cTVq1cH03meR0REWZazqBFg6Y0Ny1qtNnKDpyiKyLLMc5bApTX2Mrxer8fDhw/jyZMnERFxfHwcp6ency0MYJlMHLMcvpnjTjhw2ekNB0ggLAESCEuABMISIIF38AAM8Q4egATewQMwBWEJkEBYAiQQlgAJhCVAAmEJkEBYAiQQlgAJdPAADNHBA5BABw/AFIQlQAJhCZBAWAIkqOQGD8vr65fPF11CpT66cWfRJVRq3fZnnTmzBEggLAESCEuABMISIIF2R4Ah2h0BEmh3BJiCsARIICwBEghLgATCEiCBsARIICwBEghLgAQ6eACG6OABSKCDB2AKwhIggbAESJAUlq1WK7rd7oxLAVhebw3Loiji4OBgHrUALK23hmVZlpHn+TxqAVhabwzLdrsdjUZjXrUALK2JYdntdiPLsjmWArC8Jobl4eFhbG1tzbMWgKU1toOn0+nEJ598krySfrtjn7ZHYNX02xz7ktsdDw8PBz+XZRl7e3vx8OHDqNfrF5bV7gisun6bY9+P2x3HhuWPL793d3djd3fXXXHg0nrj3fButzt4xnJ/fz+KophLUQDL5o1/dSjLsnj06FE8evRoXvUALCW94QAJhCVAAmEJkEBYAiQQlgAJvLAMYIgXlgEk8MIygCkIS4AEwhIggbAESCAsARIIS4AEwhIggbAESKCDB2CIDh6ABDp4AKYgLAESCEuABJWMWbK8PrpxZ9ElVOrrl88XXUKl1u34rDNnlgAJhCVAAmEJkEBYAiQQlgAJtDsCDNHuCJBAuyPAFIQlQAJhCZBAWAIkEJYACYQlQAJhCZBAWAIk0MEDMEQHD0ACHTwAUxCWAAmEJUCCiWFZFEUURREREWVZDn4GuIwmhuXjx4/j7t27sbGxEbu7u5Hn+TzrAlgqE++G3717N87PzyMiIsuyedUDsJTe+OiQkAT4wcSw7Ha70W63IyLi5OTEpThwqU0My2azOTizzPM8tre34+zsbOyy/Q6ePp08wKrpd+70JXfwlGUZ9Xo9In4Iy7IsoyzLsWeXOniAVdfv3OlL6uApiiLu379/YX6tVqu4PIDVMDYs8zyP/f39wXSn04lGo+GGD3Bpjb0Mz7IsNjc34+DgILIsi7Ozszg6Opp3bQBLY+KYZb1eH4xZAlx2esMBEghLgATCEiCBsARIICwBEnhhGcAQLywDSOCFZQBTEJYACYQlQAJhCZBAWAIkEJYACYQlQAJhCZBABw/AkEkdPJWcWfY7eGYdlK963850/fO0TvsSsX7788WX3y26hMqs27GZ9f5c27gVdzZ+FTdv3hyZv1KX4cOvqVx167QvEeu3P3/80xqF5Zodm0Xtz8zCMjX9Z/EtMYttr9P+zGK/Z7G+1GVTzwJncbY4i237rFWj6m3PLiwT038W3xKz2PY67c8s9nsW60tdNvUscBZni7PYts9aNarediU3eP7vF1fjyu3Rv9Bx68X1C3+1Y5zU5WaxTttevW0/aP78wnIv//lBPGjevDD/v11uFuuctNyV31z8XSzb7/yybvvs7GxkeqPX6/WS1gZwia3UDR6ARRGWAAmEJUCC/weZrbwWe68c8wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.matshow(Q)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bc805b06-05b6-443a-8ba8-afb49fbd18c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAFjCAYAAACwvfn4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkJklEQVR4nO3dP3KjyNsH8O++tQfASifrCTfD9gkGhbMRsssHGHQDUXMCFbqB5AO4bIg2hRtoIduQ/kUTjuAG/QZeWCGJtmQjodZ8P1Wu0b8RT0PzqOmm4TellAIREe30f30HQER0zpgkiYg0mCSJiDSYJImINJgkic5YWZZ9h/Bhppfhl0mSZVnC931EUYQoirBYLLBYLPoO692yLMNwOMTnz5/7DuWi9bmefd+HZVmQUsL3fVxdXb0Zx+fPn3F1dQXf988mOS0WC0gp+w7j/dQvIE1T5TiOKoqi8XoYhspxnHd/73w+33ptMpko13Xf/Z2HiONYCSHe/NyuOC9Z1+Xddz13aT6fqzzPG68FQaAsy1Jpmu78P3EcK9d1led5pwjxIOcY075+iZbkly9f6l/lda7rQgiB8Xj8ru+N43jrteFwiPv7+3d936EGg8Fen9sV5yXrurz7rueuSCmRpimEEI3XLcvC3d0d5vP5zv93Li3HXUajEWazWd9hvMvFJ0nf9yGEgOM4re+/53Cg7f84jgPXdd8V6zEYf6hzoEso73w+h+/7O98bj8d4eXnZer0sy62kek4cx8Hz83PfYbzL730HcGxRFMG27db3q4oVRREmkwmSJIHv+xgMBhiNRgBeK+DPnz8RBAEAIEkSxHEMKWX96ziZTJBlGXzfh5QSeZ7Xn60q/OPjI6SUWK1WSNMU8/kci8UCg8EAz8/P+P79eyPWsiyxWCwghEAcxxiPx9qybGqLs4qpakVXLa8gCLTL3FUWKWVj3QCo/39ZlpBSwrIsCCHeXK+bZrNZvX2klJhMJo04NuMfDoc7y6tbP4fGVOliPeniavucbdsYDAaIoqjxY/z333+3NgSAt9flW7FWffrX19dI0xSfP3+GbduYz+cIw7CxDCklhBBbjQUhBLIsO6gOn4W+j/ePDYAKgkD7GSFEox8xDEMFoNEnNJlMGv0qcRwr27a3vitN063+q6pPK47jxjInk0ljmZvfN5lMGjEIIRr9qruWtaktzmp5cRyrNE3rWN5aZltZqn6yMAwbfYJ5ntfP91mvFdd1G8vI87zRf9wWf1t52+wT0671/NH11GaznOuq9bir3zsMQ6WU2tkn+da63CdW13XrZRRFsbU+1t9XSinHcbbKOp/P39wXz9HFH24DwM+fPw/6vGVZsG27cfjy/fv3dx/KDQYDSCkbv/Sbh0a2bW99t5QSSZI0/s/684+wLAtZlsFxHNi2Xbca3lpmW1nWYw/DsO4fE0Lg5uamXuY+6zXLMiRJsrWM1WpVx9IW/3vWw3u2dRfraZd9DpvH4zGiKKrXsZSytXW2z7rcJ9b1I7JqxL16X0q51bIdjUZbfaeDweDgffEcXHyS3KdiSilxe3ur/YxlWfWO+d44Nr/vrdM5wjCE53koyxJZlmG1WmG1Wr1r+fvEtO8yd5Wl+ky1o1xdXeH6+hqz2Ux7eLVrvf799987Y6sObXXxd2Gfbf3R9dSm6p7QEULAtu26bzLLstZ18d51uRmrbdtbA0PV/0mSBJZlIUmS+i/P8639bp998RxdfJJ0HEdb2av3dP05+zjGxs+yDKPRCC8vLxBCdJIU1uPctTN2scw4jpGmKe7v7zGfzw8e1dx3lPatZAIcZ7sAx9k2hxiPx3VLTbceuhrxHo/HmE6nKMsSs9ms0UqsWr+O49R/QRBsnWWwWq1OfqZAFy4+SQZBgNVqhSiKdr7v+z48z3uzM7ksS5RlqT2s6VJZlvjy5Qu+f/8Oz/NgWVbj8Oq9dHF2sczqBH3btjGZTJCmqXZUc9d6dRxn5/L2afFves92eWtbH2vbAPu3tu7u7pBlGaIoqrszdulyXT4+PiJJEriuC8/z6td3dRUB2wm6LEsjJz9cfJK0LAthGGI6nW5tyGqH3tWflWVZYyNPp1N4nle3GNYrs65PqM1bv/BSyq0dtTr8OWTHPyTO9y5zvSzVqO9mDJW31ivwutM5jtPo46uW/9bpVe/ZLvvEtK6L9bRP/JuqMyaA13pdnVaja0m+d11uxpqmKSzLqs8tXuc4Dm5ubrYaIpunKu3TlXCW+h45OpWiKNRkMlHz+bzxt0s1QhqGoYrjWAVBsHNUbjKZ1N+p1OsoqOu6jRH1Xa9VMyccx6lHZ6vPTCaTepS0+v44jlUcxyrP83oUcfP/6GzGGcexchxHWZalgiDYGtndZ5mbZanWV7VewzBUYRiqIAjq8uy7Xtdjqb5vvYy6+HeVV+etmNrW80fXk87m6Hae58rzPAVAeZ5Xl7dax0q91u9qGUKIxnrXrct9Y63OAqi+37btrfVbLWM9rnWu627NejPBL5MkD3HoaSS0n3Ncr+cYUxAEjdNx+lYURSPBFUVRT/U9ZAroqabrdu3iD7eJTDOZTFqnHvZhsVhgPB7Xh8rVaVNBECBN072+YzabvXv6b9+YJInO0P39fetg46nZtl3PqlmXJAmGw+Gb/7+axfTRM0h603dT9txU/V3YY6YO7e8c1+s5xrQuDMOt/ta+VP21633Nb/WtVs5x3R7iN6V4t0QiojY83CYi0mCSJCLSOPtLpVWT56uTbKsZDibKsgzfvn3be0Tw3FUXTwCA5XKJx8dHY7dNVY6yLLFcLnF/f2/eJb128H0f379/N3a7AP+d/F7N7NHNhjqKvjtF37J+Dlt1wq6JqpOMDVjle1vvkA+C4OzONzzE+m0R5vP5yW/XcAxVfTPxBO511Yn0AHbehuXYzvpwe9dVRLq6VNipua57ES2TSpZlmE6n9XPXdZFlmZFXeQFer+qzvn1MbnlVqovfmu76+hpFUaAoCsRxfPJtc9ZJMkmSrauGDAaDzi8mQYezbRuPj4/182qur4lXeQGaV4EKw9DYE58rm9d3NF11+bo+nHWfZNsFAbq8piK93/pO+Pz8DMdxjG6BZVmG5+dnDIfDxlVuTFOWpdHbYVNZlvWJ9cvlEuPx+KQt5LNOkm3O+a5wv6KqEps+IFVdoby6P7upLbGXlxejk/ym9cFaIQSGw2HjikjHdtaH27uu5LxarS7qV/IS+L7fS1/RMViWhdFohNFoZOSPcZIkuLu76zuMTq33c1dnuZyy7/usk2TbXE/dRUbptGazWX3nwupitaZJkgRXV1f18/W7Cpro5eUFi8Wivk/PdDo1th8/yzJ8+fJl6/VT9n2f9eH2Zr+DlBI3NzfGt1gupc+oujlUlSBNPcwbDAaNH+Qsy+or3Zhms2ExHo9P3ofXJSFE46LY1ZXRT7n/nP3cbSkl5vM5bm9vsVwujT0xtroH9mw2w2Qywe3trbF9XsDrdtm8FL9lWSiKoqeIPiaKorprJ45jBEFgbGIB/rtCfHV7kkPv2X5OqkkLlmUhz/N33xnzvc4+SRIR9ems+ySJiPrGJElEpMEkSUSkwSRJRKTBJElEpMEkSUSk0cnJ5H/88cfWOXM/fvzAp0+f3vy/XX/uV102Y/x1Yuxz2b9CjHme459//vnvhS4uSvn169e9Xtv3/37kc7/qshnjaT/3qy77V4hx87WjHW4/PDx0+rk+l31IjF2Xp88Y+yrLIZ9lPesG65nG3in3wGzctVMs45QuqTwsy/m6pPKcqiwna0l27RgtgT5dUnlYlvN1SeXpqyydzN2+vr7Gp0+f8PDwcFEbhYh+HU9PT3h6esKPHz8aF5DuJEn++eef+Ouvvz76NUREvdvMZ8YcbhMR9YFJkohIg0mSiEiDSZKISINJkohIg0mSiEiDSZKISINJkohIo5Mk+ePHD/z55594enrq4uuIiE7u6ekJf/75J378+NF4nTNuiIjWcMYNEdEBmCSJiDSYJImINDq5x03XsixDkiQAgOVyicfHR1iW1Xj/27dvjcsZGSHLgG/fANPi1kiSBFJKCCEAAI7j9BzR+7xV50wipUSSJBgMBpBSwnXdevuYoG3/llIiiiIIISClhOd5p9lGbVfnzfNcBUGgwjBUQRCooij2vpLvRwVB0Hhs23b9PAxDlaap0oR+nsJQqTRVyrS4NeI4Vp7nKaVe64sQoueI3k9X50yzXhalVL2NTKDbv9e3SZ7nynXdo8Swmc9a99hDAuoySaZpqizLaiwbgMrzvPE545JkxdS4dxBCNH48N7eRKfatc6bYTPAmJcnK5v6d5/lWuda3WZf2un2DlLLxXAhRH4ocm23beHx8rJ+XZQkAGAwGJ1k+7UdKidVqBcuykGUZyrI06pBu3aXVucFggOvr6/qwezgc9h3Sh1XdB+sGgwGyLDv6sncmyT4DAgDXdevHz8/PcBzH2P6hS5VlGQaDQd1HtFgsEEVR32G92yXVuTAMAQCfP39GGIaNspmq+uHatFqtjr7snQM3fQa0GUcUReYN0PwCVqsVpJR1MvE8D1dXV1Afn5vQq0uoc0mSIAgCSCkxHo8BAPP5vOeojqMtV3XpoNHttoCqaYmVrm4I5vs+4jg29hf9kgkhYFlWvW2qf7Msg23b/QX2QabXOSkllsslgiAA8Hq2wfX1NXzfN7Y7BHitX5uNtKq756OqG4BVNqcl7kyShwb06dOnzqclzmazesNWydnUinuJTN7h2lxCncuyDLe3t/VzIQS+f/9+khbXMTmOs7M1fHNz8+Hv3mzUrTf4gJY+ybZz3boIaB9RFMG27bqyvry87Kysxm54U+NeI4TAzc1NvQ2qcyVNbUXuW+fOnW3bWC6Xjdd+/vxp5HZZ3783f5SllLi5uTnJNtrZkuwzICklRqNR47Wqzwt47W+J4xgAMJ1OcXt7a0bHdJIA/8aN6RS4vQVMiFsjDEP4vo/r62ukaVpvF9O8VedMIoTAcDjEbDar99eqX9IEuv27qm+3t7dYLpf1ANWxtV4FSEqJ+XxeB/T9+/fWJMmrABHRpdjMZ60DN0KIuvPXiJYaEdER8AIXREQaTJJERBpMkkREGkySREQavBEYERF4IzAior3wRmBERAdgkiQi0mCSJCLSYJIkItJgkiQi0mCSJCLSYJIkItJgkiQi0uCMGyIicMYNEdFeOOOGiOgATJJERBpMkkREGkySREQaTJJERBpMkkREGkySREQaTJJERBqccUNEBM64ISLaC2fcEBEdgEmSiEiDSZKISINJkohIg0mSiEiDSZKISINJkohIg0mSiEiDSZKISIPTEomIwGmJRER74bREIqIDMEkSEWkYkSSTJMFisUCSJEiSpO9wDpJlGa6vr7Wf8X0fZVmeJqAOSSnh+z4Wi4WRZWjbNlEUoSxL48pTMT3+LVn2+gcAUv73+FRUB75+/drF1+wUx7HyPE8ppVSe50oIcbRldS0MQ5WmqdKt5ur9oihOF1hHhBB13Gma1tvJBLptA2DrLwiCHqJ8H9Pj3+J5SgGvf46j1JH3lc181rr3pmmqbNt+15d2aX1HVOo1UZpGlyTDMNwqowniON76weroN/ekNmMuikKFYdh4zaQEY3r8O83nr4nxRPvIZj7bebgdRRGA18ORPkkpsVqtYFkWsixDWZYQQvQaU5eiKILrun2H8S5th3J915kurG8TE7eR6fHvZFmvfz34fdeL57JSsyzDYDBAFEVwHAeLxQJCiLOJ7yPKsoTV00bvgm3bkFLWz6vkuFqt+gqpE+vbpCxLrFYro36YTY9/p7IE/m24YbkExmPghGXamSTPxWq1gpQSjuPAsix4noerqyuoj5/a2buXlxd4ntd3GO8mhEAQBFgsFri7u6sT5mAw6Dmy7vi+jyAI+g7j3UyPv+Z5/7UihQCGQyDPT7b4TpJkNeOm8vDwgIeHhw9/rxAClmXVv47Vv1mWwbbtD39/X5Ikwd3dXd9hfNhkMoGUsv4hA2B+q+VfZVkiSRJjW/umx98gJVDt70K8Ppeys9bk09NTY7bg5oybTpLkp0+fjjLj5lJ2uF1eXl7qx1JKTKdT3N/fG5X8pZT1Nqp+uC5ipwTw999/G10W0+OvZRnw5QtQFM3XOzxi2WzUrTf4gDM/3BZC4Obmpu6/q3ZKkxJJZb0Psmp1VcbjMcbjsXE/CtfX1/jf//4Hy7Iwn8+NPbTb1T9c9YebyvT4a0IA6/UqSQDXPekgzlknSQAIwxC+7+P6+hppmiKO475D2luSJHW80+kUt7e3jUGnsiyxWCwAAEEQYDweG/UDEAQBkiTBarXCaDTaSv7n7K1tA5h/JGN6/ABek+HNDTCbvT7OcyAMTxqC9gIXv/32G4qieLPZzgtcENGl2OsCF0mSwPd9AK+/stV5k0REv5qdh9uO48BxHGP7mIiIumLEBS6IiPrCJElEpMEkSUSkwXvcEBGB97ghItoL73FDRHQAJkkiIg0mSSIiDSZJIiINJkkiIg0mSSIiDSZJIiINJkkiIg0mSSIiDU5LJCICpyUSEe2F0xKJiA7AJElEpMEkSUSkwSRJRKTBJElEpMEkSUSkwSRJRKTBJElEpMEZN0RE4IwbIqK9cMYNEdEBmCSJiDSYJImINJgkiYg0mCSJiDSYJImINJgkiYg0mCSJiDQ444aICJxxQ0S0F864ISI6AJMkEZEGkyQRkcbvfQfwJimBKAKEeH3seYBl9R3V3rIsQ5IkAIDlconHx0dY/8YvpUQURRBCQEoJz/Pq985RlmX49u0b0jTdeh0AbNuGlBJlWcK27T5CPIhu21Svl2WJ5XKJ+/v7sy5T27YxrRxAe1neeu9oVIs0TVUQBCoIAuW6riqKou2j6uvXr63vfZht//c4z5Vy3eMt6wiCIGg8ttfKs/44z3PlnnHZwjBUaZqqXVXG8zwFQAFQjuNo68o50W0by7JUmqZKKaXm87kSQpw8vn3pto1J5VBKXxbde13azGetS9NVoLe+tDN53kySSillWcdZ1hGkaaqstXjzPFcAVJ7nKs/zrXVqGVC2XRV0Pp+roiiMSY5K6beNUkrFcVy/N5/PtfX/XOzaNiaWQ6ndZdnnvS5s5rOdfZJZlmE6ndbPXddFlmWQUh61VbslSYDBoPnaYAD8e3h37mzbxuPjY/28LEsAwGAwQJIkGGyUbTAY1IeuprEs66y7Cjbptg0AOI5TvxeGIcbj8Unj68qllKNPO/sk36pAJ/PvcresVicN4yNc160fPz8/w3EcWJZVr9NNK4PKVinLElEUAXjt2xuPxxBC9BzV29q2TSXLMjw/P2M4HMLzvB4i7MallKMvrQM3b1WgddWMm8rDwwMeHh66i3JTW/I8Y1UieavDuS15nrP1ASchBIbDIfI87zeoA7RtG9u2IYSA7/uIoqixT5jkUspxLE9PT43Zgpszbt48uC+KQgkh+hm4mc9390mu9bOYwvO8ur9Lqd39Q5ZlNfqQztGuKlMNDCj1Wl+w1rdngs1tsymOYwXg7Ptc39qdTSmHUgb0Sa7zfR9xHPfT37TWn9Jwc3PaOD5oNpvB930IIVCWJcqybPQVrbsxrGxZluHLly9br5+8a+addm2bJElwdXVVf6bqOjh5n/wHXUo5+qZNkrsq0Elt9mtJ+ZogDRogiKKoPtwpyxIvLy+wLGurz05KiZubGyMGP9brgRACQRDUz5Mkgeu6RpSjbdsMBoPGj1iWZbAs6+zPLwSa28bkcgD6rqdT5qLWPsldFaiXTt8wBHwfuL0FlsvX54aQUmI0GjVesyyrXo9hGML3fdze3mK5XCI847IlSYI4jgEA0+kUt7e3dTK8ubnBbDaDZVnI8/ysy1HRbRvbtnF/f4/FYgEAiOP4tCcvH6ht25hWDqC9LG+9d0w7rwIkpcTnz58br1mWhaIodn4JrwJERJdiM5/tbEkKIbAjdxIR/XJ4gQsiIg0mSSIiDSZJIiINJkkiIg3eCIyICLwRGBHRXngjMCKiAzBJEhFpMEkSEWkwSRIRaTBJEhFpMEkSEWkwSRIRaTBJEhFpcMYNERE444aIaC+ccUNEdAAmSSIiDSZJIiINJkkiIg0mSSIiDSZJIiINJkkiIg0mSSIiDc64ISICZ9wQEe2FM26IiA7AJElEpMEkSUSkwSRJRKTBJElEpMEkSUSkwSRJRKTBJElEpMEZN0RE4IwbIqK9cMYNEdEBmCSJiDSYJImINH7vO4BdsizDt2/fkKZp4/UkSQAAZVliuVzi/v4etm33EeK7RFEEx3EAAJZl9RtMB7Isq7fJcrnE4+OjUeVqq2dSSkRRBCEEpJTwPM+ockFKIIoAIV4fex5gUvwbeq9nqkUcxyqOYxWGoZpMJipN07aPqq9fv7a+d6gwDFWapmpXaJZl1XHM53MlhOhsuacAYOsvCIK+w3q39diDIFC2bfcYzWF09Wy9HHmeK9d1Txnax61vhzxXyrT4N5y6nm3ms9YkeUhC6jJJ1oHtqLxxHNeP5/O5UTtlURQqDMPGayYnyDRNlWVZ9fM8zxUAled5j1EdbrOe5Xm+Va/Wy3n28ryZJJVSyqT4N/RRzzbzWWufZBiGjUPZczjcqA5Vgdf4xuNxj9EcznXd+nEURY3nprFtG4+Pj/XzsiwBAIPBoKeIupEkyVYZBoMBsizrKaIDJQmwuQ0GA8CU+DecQz1r7ZM814SUZRmen58xHA7heV7f4ext/UemLEusVisIIfoLqAPrSf75+RmO45zFj+lHVDvhptVqddpA3qslfpgS/w591zPt6HaWZfB9/6wSkm3b+P79O/I8RxRFfYfzLr7v4+7uru8wOlOWJaIoQhiGfYdyNG3J0ximx4/+6pl2dNu2bQgh4Pu+9vCwmpZYeXh4wMPDQ7eRrrEsC6PRCMPhEEVRGNV6KcsSSZIYFfNbfN9HHMcXUSbLsrZajavVypyyWdZ2q3G1Mnp0u3Ksevb09NSYUr05LbF14GZdHMcKgCqKYq+Ozi5shhbH8c4OXN2o+zmK49ioAae3BEFQd6IXRdFaR87VZj1rG7gxplxtAzemxN/ilPVsr4GbJElwdXVVP6/6zqSUnWbwt6wf4gwGg0Y/aZZlsCzLqPMkgde4TR/cqERRVB9tlGWJl5cXc1pca9br2WY/sZQSNzc35pRrs59bSuDmxuiWZN/1bOfhdp8JKUkSxHEMAJhOp7i9vYXrurBtG/f391gsFgCAOI63TgI2hekDNsBr8hiNRo3XLMs6m77rt7TVM+B1oNL3fdze3mK5XJrX1xqGgO8Dt7fAcvn63FDnUM9arwIURVHdNxPHMYIgaN25eRUgIroUm/msdeBmfZDGlNYBEVHXeIELIiINJkkiIg0mSSIiDSZJIiIN3giMiAi8ERgR0V54IzAiogMwSRIRaTBJEhFpMEkSEWkwSRIRaTBJEhFpMEkSEWkwSRIRaXDGDREROOOGiGgvnHFDRHQAJkkiIg0mSSIiDSZJIiINJkkiIg0mSSIiDSZJIiINJkkiIg3OuCEiAmfcEBHthTNuiIgOwCRJRKTBJElEpMEkSUSkwSRJRKTBJElEpMEkSUSkwSRJRKTBJElEpMFpiURE4LREIqK9cFoiEdEBmCSJiDR+7zuAXbIsw7dv35CmaeN1KSWSJMFgMICUEq7rQgjRU5SHy7IMSZIAAJbLJR4fH2FZVr9BfcS/ZUFZAsslcH8P2HavIb1XlmUAANu2IaVEWZawWZazEEURHMcBgH72F9WBr1+/dvE1SimlwjBUaZqqXaEFQdB47nleZ8s9hfX4gyBQtm33GE0HLEupNH19PJ8rJUS/8XyA53kKgAKgHMdRRVH0HdK7XVJZlFJ1Wdb/NnNBlzbz2V6H277voyzLoyTpTa7rtv7qPT8/nySGY8iyDNPptH7uui6yLIOUsseoPigMmy1Hg1vF19fXKIoCRVEgjmOjW/iXVJayLBGGIZRS9V8QBJhMJieL4c0kmWUZZrPZKWJ502AwwPX1dX3YPRwO+w5pb7Zt4/HxsX5e/egMBoOeIurAv4dAAF4T5njcXywdsCzL6ISy7pLK4rpu/TiKosbzU3izT1JKeTb9fmEY4suXL/j8+TM8z8N8Pu87pIOsb9zn52c4jmN+Rc4y4PkZGA4Bz+s7mncryxJRFAF47S8ej8dnU+8PdUllWd8/yrLEarU6fVl0x+ZhGCqllBJCaPs1uuyTrOwKLQxDFcexms/nCoBxfZKVoijeXKdGKQqlPE+pf+uLida3RZqmShjcv3pJZVnned5J9pm9+yTLsty7lVPNuKn+jjHzRkqJ5XIJx3HgeR7yPMfLy4uRfXq+7xvfV9RgWcBo9Pp3or7rrq3XIyEEpJRG1i3gsspSKcsSSZIcZZ+pZtpUf5szblpbkvP5vH58Di3JMAzrlm0lCAKVVqOrhgiCQOV5rpR6/cU3tjUZx6+j25U8Vwr4b7TbIGmaKmutLEVRKABGbptLKsu6OI5PdjbIXi3JJElwd3fXecY+1PqIum3bWC6Xjfd//vxp1PlfURTBtm0IIVCWJV5eXsxtTQ4GzYGbLHttURq0PSpCCARBUD9PkgSu6xq5bS6pLOuyLOttkLN14Obl5aV+LKXEdDrF/f390ZNSkiSI4xgAMJ1OcXt7W580PhwOMZvN6g0+Nmg0VUqJ0WjUeM2yLHimDnbY9uvJ44vF6/M4BjZO/jeFZVm4ubmp61ae5wjDsO+w3uWSyrKpr8GnvS5w8dtvvyHP89YgeYELIroUB13goizL+hzJIAjq6U5ERL8K7XmSlmVhMpmc9Ox2IqJzwqsAERFpMEkSEWkwSRIRafAeN0RE4D1uiIj2wnvcEBEdgEmSiEiDSZKISINJkohIg0mSiEiDSZKISINJkohIg0mSiEiDSZKISIPTEomIwGmJRER74bREIqIDMEkSEWkwSRIRaTBJEhFpMEkSEWkwSRIRaTBJEhFpMEkSEWlwxg0RETjjhohoL5xxQ0R0ACZJIiINJkkiIg0mSSIiDSZJIiINJkkiIg0mSSIiDSZJIiINzrghIgJn3BAR7YUzboiIDsAkSUSkwSRJRKRxlkkyyzLMZjPMZjOMRiOUZbnzc77vt7537kyOvSKlhO/7WCwW5pcny4Dr676jeLcsy3DdEr/uvXPUFm+SJEiSBFEUwfd9ZFl2moBUizRNVZqmSiml8jyvH+/y9evX1vfeIwiCxmPbtnfGB0AVRdHpsk/B5NjXCSHqMqRpqjzP6zeg9wpDpdJUqfbd4ayFYVjXqUPeO0e6eC3LqvPQfD5XQoijxLCZz1rXnOd5CoACoBzH0e7QXSbJNE2VZVn18zzPFQCV53njc2EYNnZSk5gceyWO461KasqO2Mrw+HXr37RtsyveOI7rx/P5fGfjqQub+az1cPv6+hpFUaAoCsRxDMuyjtymfWXbNh4fH+vn1SHcYDCoX4uiCK7rniSerpkc+7q2Q+uTHQLRL8dxnPpxGIYYj8cnWe7vujdPlRg3rSeR5+dnOI5Tx1KWZW9xfZTJsW+ybRtSyvp5lRxXq1VfIdEvIMsyPD8/YzgcwvO8kyyzNUmWZYkoigAAy+US4/EYQoidn61m3FQeHh7w8PDw4eCqGNI0rV97eXk52crpmsmxbxJCIAgCLBYL3N3d1QlzvcVP1DXbtiGEgO/7nR2VPT09NWYLbs64ae2oWO8vS9NU20na9cBNxfO8Rl9kHMeNuEzq1zM5dp1qUK8oCvMHowzrt9uk2Z0vok9yXRzHR6tvm/mstSUppYRt2wBeWw1SSkgpW1uTXZvNZvB9H0KIRv/Xy8tLI8bpdIr7+/s61nNmcuy7rNeHLMtg2/bFdCfQeUmSBKPRCEVRAEBd79bz1NHsyqSbI8xvtRK6bkmGYViPZBVFoebz+c7PYceotylMjr1iWVZdJzzPa4w+GglQyuCWsG4f1b13jjbjTdNUua5bPw/DsJGjurRXS7Lqb6okSQLXdU/SSpBSYjQaNV6zLKvRl1eWJRaLBQAgCAKMx2NjWmMmx74pCAIkSYLVaoXRaNQYfTRKkgBx/Pp4OgVubwGDzkBIkgTxv/FPp1Pc3t7WfXW6985RW7y2beP+/r7ed+I4boxVHFPrVYCyLEOSJLAsC3meN5LmJl4FiIguxWY+a+2TtG3b2BYOEVFXznLuNhHRuWCSJCLSYJIkItJgkiQi0uCNwIiIwBuBERHthTcCIyI6AJMkEZEGkyQRkQaTJBGRBpMkEZEGkyQRkQaTJBGRBpMkEZEGZ9wQEaF9xk0nSfLTp0/466+/OrlDYptLS8CXVB6W5XxdUnmOXZaHhwf89ddf+PTpU+N1Yw63L2ljA5dVHpblfF1Sefoqy9GS5L4FOkbBu172ITF2XZ4+Y+yrLId8lvWsG6xn7ZgkO/zcoZ/t8vt+xcp7jM8dgvXs9N/Z9fft89lOrgL0xx9/4PPnz43Xfvz4sXVsv0vXn/tVl80Yf50Y+1z2rxBjnuf4559/6uedJEkioktlzMANEVEfmCSJiDSYJImINP4fQKLX31ckDckAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (4, 4))\n",
    "ax.matshow(c, cmap = mpl.colors.ListedColormap([\"white\"]))\n",
    "ax.set_title(\"Optimal transport plan (Monge)\")\n",
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        if (Q[i,j] > 0):\n",
    "            color = 'red'\n",
    "        else:\n",
    "            color = 'black'\n",
    "        ax.text(j, i, int(c[i, j]), horizontalalignment = 'center',\n",
    "                verticalalignment = 'center', weight = 'bold', color = color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c05324e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}