Loss matrix of EMD in Python optimal transports

cost_matrix.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "347a3f1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "30265d7f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4, 5, 6, 7]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.array(range(8)).reshape(1,-1)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f75b5486",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 2, 3, 4, 5, 6, 7],\n",
       "       [0, 1, 2, 3, 4, 5, 6, 7],\n",
       "       [0, 1, 2, 3, 4, 5, 6, 7],\n",
       "       [0, 1, 2, 3, 4, 5, 6, 7],\n",
       "       [0, 1, 2, 3, 4, 5, 6, 7],\n",
       "       [0, 1, 2, 3, 4, 5, 6, 7],\n",
       "       [0, 1, 2, 3, 4, 5, 6, 7],\n",
       "       [0, 1, 2, 3, 4, 5, 6, 7]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.repeat(x,8, axis=0)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "13eea9d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0],\n",
       "       [1],\n",
       "       [2],\n",
       "       [3],\n",
       "       [4],\n",
       "       [5],\n",
       "       [6],\n",
       "       [7]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = np.array(range(8)).reshape(-1,1)\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5271daf4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0, 0, 0, 0, 0, 0],\n",
       "       [1, 1, 1, 1, 1, 1, 1, 1],\n",
       "       [2, 2, 2, 2, 2, 2, 2, 2],\n",
       "       [3, 3, 3, 3, 3, 3, 3, 3],\n",
       "       [4, 4, 4, 4, 4, 4, 4, 4],\n",
       "       [5, 5, 5, 5, 5, 5, 5, 5],\n",
       "       [6, 6, 6, 6, 6, 6, 6, 6],\n",
       "       [7, 7, 7, 7, 7, 7, 7, 7]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = np.repeat(y,8, axis=1)\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4d12443f",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.empty((8,8), dtype=np.complex128) #c stands for coordinate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "654fa579",
   "metadata": {},
   "outputs": [],
   "source": [
    "c.real = x\n",
    "c.imag = y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d570c6ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j, 4.+0.j, 5.+0.j, 6.+0.j, 7.+0.j],\n",
       "       [0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j, 4.+1.j, 5.+1.j, 6.+1.j, 7.+1.j],\n",
       "       [0.+2.j, 1.+2.j, 2.+2.j, 3.+2.j, 4.+2.j, 5.+2.j, 6.+2.j, 7.+2.j],\n",
       "       [0.+3.j, 1.+3.j, 2.+3.j, 3.+3.j, 4.+3.j, 5.+3.j, 6.+3.j, 7.+3.j],\n",
       "       [0.+4.j, 1.+4.j, 2.+4.j, 3.+4.j, 4.+4.j, 5.+4.j, 6.+4.j, 7.+4.j],\n",
       "       [0.+5.j, 1.+5.j, 2.+5.j, 3.+5.j, 4.+5.j, 5.+5.j, 6.+5.j, 7.+5.j],\n",
       "       [0.+6.j, 1.+6.j, 2.+6.j, 3.+6.j, 4.+6.j, 5.+6.j, 6.+6.j, 7.+6.j],\n",
       "       [0.+7.j, 1.+7.j, 2.+7.j, 3.+7.j, 4.+7.j, 5.+7.j, 6.+7.j, 7.+7.j]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "012d465c",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist = np.abs(np.subtract.outer(c,c)) #歐氏距離 (Euclidean distance)，可採用其他距離，例如 Chebyshev distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7ee85dc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 8, 8, 8)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dist.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f835baf6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dist[3,2,7,5]  # 可以解讀成是 (3,2) 和 (7,5)的距離"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "addc7478",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4142135623730951"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dist[0,0,1,1] #根號 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3b3f1745",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGdCAYAAAAv9mXmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAZLElEQVR4nO3dfUxVh/3H8Q9w5agVKCgoTERtrVYtzIoaR7s+aDXEmnZ/OGNsRnVbUoNTS5o0/DObLPO6P7bYLgYfZrFJ57Rbhu2aKVNXMUtlIoZE24Rq6yrTAnUiT65Xx72/v3o3flbkHPlyPPh+JSfZvTnX8wlhvHvv5SEhFovFBADAAEv0ewAAYGgiMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwERosC8YjUZ16dIlpaSkKCEhYbAvDwC4A7FYTJ2dncrJyVFiYt/PUQY9MJcuXVJubu5gXxYAMICampo0fvz4Ps8Z9MCkpKRIkjIyMm5bv7vNxIkT/Z7gSZCDPnbsWL8neJKenu73BE8cx/F7gieRSMTvCZ61tbX5PcGV69ev680334x/Le/LoAfm65fFEhMTAxeYUGjQP1wDYtiwYX5P8CyoX/CGDx/u9wRPgvrxDvLL7UP5Yx6sr/AAgMAgMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMCEp8Bs3bpVEydO1PDhwzVv3jydOHFioHcBAALOdWD27dunsrIybdy4UadOnVJBQYEWL16s1tZWi30AgIByHZhf/epX+vGPf6xVq1Zp+vTp2rZtm0aOHKk333zTYh8AIKBcBeb69euqr6/XwoUL//sPJCZq4cKFOn78+Dc+JhKJqKOjo9cBABj6XAXm8uXL6unp0dixY3vdP3bsWDU3N3/jY8LhsNLS0uJHbm6u97UAgMAw/y6y8vJytbe3x4+mpibrSwIA7gIhNyePGTNGSUlJamlp6XV/S0uLxo0b942PcRxHjuN4XwgACCRXz2CSk5M1e/ZsHTlyJH5fNBrVkSNHNH/+/AEfBwAILlfPYCSprKxMJSUlKiws1Ny5c7VlyxZ1d3dr1apVFvsAAAHlOjDLly/Xl19+qZ/+9Kdqbm7Wt7/9bR08ePCmN/4BAPc214GRpLVr12rt2rUDvQUAMITwu8gAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACU9/D2YgpKenKykpya/Le5KZmen3BE+ys7P9nuBZXl6e3xM8GTdunN8TPBk5cqTfEzy5du2a3xM8a25u9nuCK1999VW/z+UZDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATrgNz7NgxLV26VDk5OUpISND+/fsNZgEAgs51YLq7u1VQUKCtW7da7AEADBEhtw8oLi5WcXGxxRYAwBDiOjBuRSIRRSKR+O2Ojg7rSwIA7gLmb/KHw2GlpaXFj9zcXOtLAgDuAuaBKS8vV3t7e/xoamqyviQA4C5g/hKZ4zhyHMf6MgCAuww/BwMAMOH6GUxXV5fOnTsXv33+/Hk1NDQoIyNDEyZMGNBxAIDgch2YkydP6qmnnorfLisrkySVlJRo9+7dAzYMABBsrgPz5JNPKhaLWWwBAAwhvAcDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATLj+ezADZdSoUUpKSvLr8p6kpqb6PcGTjIwMvyd4Nm7cOL8neDJp0iS/J3iSnp7u9wRP2tra/J5wz7h27Vq/z+UZDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATrgITDoc1Z84cpaSkKCsrS88//7waGxuttgEAAsxVYGpqalRaWqra2lodOnRIN27c0KJFi9Td3W21DwAQUCE3Jx88eLDX7d27dysrK0v19fX67ne/O6DDAADB5iow/197e7skKSMj45bnRCIRRSKR+O2Ojo47uSQAICA8v8kfjUa1YcMGFRUVaebMmbc8LxwOKy0tLX7k5uZ6vSQAIEA8B6a0tFRnzpzR3r17+zyvvLxc7e3t8aOpqcnrJQEAAeLpJbK1a9fq/fff17FjxzR+/Pg+z3UcR47jeBoHAAguV4GJxWL6yU9+oqqqKh09elSTJk2y2gUACDhXgSktLdWePXv07rvvKiUlRc3NzZKktLQ0jRgxwmQgACCYXL0HU1FRofb2dj355JPKzs6OH/v27bPaBwAIKNcvkQEA0B/8LjIAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEy4+oNjAykpKUmhkG+X92TYsGF+T/DEcRy/J3g2cuRIvyd4kp6e7vcET6ZOner3BE8aGxv9nuBZS0uL3xPM8AwGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMuApMRUWF8vPzlZqaqtTUVM2fP18HDhyw2gYACDBXgRk/frw2b96s+vp6nTx5Uk8//bSee+45ffTRR1b7AAABFXJz8tKlS3vd/vnPf66KigrV1tZqxowZAzoMABBsrgLzv3p6evT73/9e3d3dmj9//i3Pi0QiikQi8dsdHR1eLwkACBDXb/KfPn1ao0aNkuM4eumll1RVVaXp06ff8vxwOKy0tLT4kZube0eDAQDB4DowU6dOVUNDg/7+979rzZo1Kikp0ccff3zL88vLy9Xe3h4/mpqa7mgwACAYXL9ElpycrAcffFCSNHv2bNXV1en111/X9u3bv/F8x3HkOM6drQQABM4d/xxMNBrt9R4LAACSy2cw5eXlKi4u1oQJE9TZ2ak9e/bo6NGjqq6uttoHAAgoV4FpbW3VD37wA33xxRdKS0tTfn6+qqur9cwzz1jtAwAElKvA7Nq1y2oHAGCI4XeRAQBMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBgwtUfHBtIPT09fl3asxs3bvg9wZNIJOL3BM+uXbvm9wRP2tra/J7gSWNjo98TPAnqx1sK3ue4m708gwEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABN3FJjNmzcrISFBGzZsGKA5AIChwnNg6urqtH37duXn5w/kHgDAEOEpMF1dXVq5cqV27typ9PT0gd4EABgCPAWmtLRUS5Ys0cKFCwd6DwBgiAi5fcDevXt16tQp1dXV9ev8SCSiSCQSv93R0eH2kgCAAHL1DKapqUnr16/Xb3/7Ww0fPrxfjwmHw0pLS4sfubm5noYCAILFVWDq6+vV2tqqRx99VKFQSKFQSDU1NXrjjTcUCoXU09Nz02PKy8vV3t4eP5qamgZsPADg7uXqJbIFCxbo9OnTve5btWqVpk2bpldffVVJSUk3PcZxHDmOc2crAQCB4yowKSkpmjlzZq/77rvvPo0ePfqm+wEA9zZ+kh8AYML1d5H9f0ePHh2AGQCAoYZnMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmLjjPzjmVVdXl5KSkvy6vCcdHR1+T/DkypUrfk/wrLm52e8J95SWlha/J3hy7do1vyd4FrTP8a+++qrf5/IMBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJV4F57bXXlJCQ0OuYNm2a1TYAQICF3D5gxowZOnz48H//gZDrfwIAcA9wXYdQKKRx48ZZbAEADCGu34M5e/ascnJyNHnyZK1cuVIXLlzo8/xIJKKOjo5eBwBg6HMVmHnz5mn37t06ePCgKioqdP78eT3++OPq7Oy85WPC4bDS0tLiR25u7h2PBgDc/VwFpri4WMuWLVN+fr4WL16sP//5z7p69areeeedWz6mvLxc7e3t8aOpqemORwMA7n539A79/fffr4ceekjnzp275TmO48hxnDu5DAAggO7o52C6urr06aefKjs7e6D2AACGCFeBeeWVV1RTU6N//OMf+vDDD/W9731PSUlJWrFihdU+AEBAuXqJ7J///KdWrFihf/3rX8rMzNRjjz2m2tpaZWZmWu0DAASUq8Ds3bvXagcAYIjhd5EBAEwQGACACQIDADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAE67+HsxAamtrU2JisPr25Zdf+j3BkxEjRvg94Z7T2dnp9wRPHMfxe4InkUjE7wmeXblyxe8Jrrj5WAfrKzwAIDAIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABMEBgBggsAAAEwQGACACQIDADBBYAAAJggMAMAEgQEAmHAdmIsXL+qFF17Q6NGjNWLECD3yyCM6efKkxTYAQICF3Jzc1tamoqIiPfXUUzpw4IAyMzN19uxZpaenW+0DAASUq8D84he/UG5uriorK+P3TZo0acBHAQCCz9VLZO+9954KCwu1bNkyZWVladasWdq5c2efj4lEIuro6Oh1AACGPleB+eyzz1RRUaEpU6aourpaa9as0bp16/TWW2/d8jHhcFhpaWnxIzc3945HAwDufgmxWCzW35OTk5NVWFioDz/8MH7funXrVFdXp+PHj3/jYyKRiCKRSPx2R0eHcnNzNWbMGCUmBuub2CZPnuz3BE8mTJjg9wTPsrOz/Z7gSUZGht8TPHEcx+8Jnvzv15iguXLlit8TXIlEItq2bZva29uVmpra57muvsJnZ2dr+vTpve57+OGHdeHChVs+xnEcpaam9joAAEOfq8AUFRWpsbGx132ffPKJ8vLyBnQUACD4XAXm5ZdfVm1trTZt2qRz585pz5492rFjh0pLS632AQACylVg5syZo6qqKv3ud7/TzJkz9bOf/UxbtmzRypUrrfYBAALK1c/BSNKzzz6rZ5991mILAGAICda3cQEAAoPAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABOu/+DYQLl8+bJfl/bMcRy/J9xz/v3vf/s9wZMgfn5L0rBhw/ye4MmNGzf8nuBZR0eH3xNccfOx5hkMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYcBWYiRMnKiEh4aajtLTUah8AIKBCbk6uq6tTT09P/PaZM2f0zDPPaNmyZQM+DAAQbK4Ck5mZ2ev25s2b9cADD+iJJ54Y0FEAgOBzFZj/df36db399tsqKytTQkLCLc+LRCKKRCLx2x0dHV4vCQAIEM9v8u/fv19Xr17Viy++2Od54XBYaWlp8SM3N9frJQEAAeI5MLt27VJxcbFycnL6PK+8vFzt7e3xo6mpyeslAQAB4uklss8//1yHDx/WH//4x9ue6ziOHMfxchkAQIB5egZTWVmprKwsLVmyZKD3AACGCNeBiUajqqysVElJiUIhz98jAAAY4lwH5vDhw7pw4YJWr15tsQcAMES4fgqyaNEixWIxiy0AgCGE30UGADBBYAAAJggMAMAEgQEAmCAwAAATBAYAYILAAABMEBgAgAkCAwAwQWAAACYIDADABIEBAJggMAAAEwQGAGCCwAAATAz6n6QM8t+SiUajfk/w5D//+Y/fEzy7ceOG3xM8CeruoP7/k8/xwfP1x7o/nyuDHpjOzs7BvuSAuXjxot8TPAnqbgB3r87OTqWlpfV5TkJskP+TJRqN6tKlS0pJSVFCQsKA/tsdHR3Kzc1VU1OTUlNTB/TftsTuwcXuwRfU7ey+WSwWU2dnp3JycpSY2Pe7LIP+DCYxMVHjx483vUZqamqgPhm+xu7Bxe7BF9Tt7O7tds9cvsab/AAAEwQGAGBiSAXGcRxt3LhRjuP4PcUVdg8udg++oG5n950Z9Df5AQD3hiH1DAYAcPcgMAAAEwQGAGCCwAAATAyZwGzdulUTJ07U8OHDNW/ePJ04ccLvSbd17NgxLV26VDk5OUpISND+/fv9ntQv4XBYc+bMUUpKirKysvT888+rsbHR71m3VVFRofz8/PgPn82fP18HDhzwe5ZrmzdvVkJCgjZs2OD3lD699tprSkhI6HVMmzbN71n9cvHiRb3wwgsaPXq0RowYoUceeUQnT570e9ZtTZw48aaPeUJCgkpLS33ZMyQCs2/fPpWVlWnjxo06deqUCgoKtHjxYrW2tvo9rU/d3d0qKCjQ1q1b/Z7iSk1NjUpLS1VbW6tDhw7pxo0bWrRokbq7u/2e1qfx48dr8+bNqq+v18mTJ/X000/rueee00cffeT3tH6rq6vT9u3blZ+f7/eUfpkxY4a++OKL+PG3v/3N70m31dbWpqKiIg0bNkwHDhzQxx9/rF/+8pdKT0/3e9pt1dXV9fp4Hzp0SJK0bNkyfwbFhoC5c+fGSktL47d7enpiOTk5sXA47OMqdyTFqqqq/J7hSWtra0xSrKamxu8prqWnp8d+85vf+D2jXzo7O2NTpkyJHTp0KPbEE0/E1q9f7/ekPm3cuDFWUFDg9wzXXn311dhjjz3m94wBsX79+tgDDzwQi0ajvlw/8M9grl+/rvr6ei1cuDB+X2JiohYuXKjjx4/7uOze0d7eLknKyMjweUn/9fT0aO/everu7tb8+fP9ntMvpaWlWrJkSa/P9bvd2bNnlZOTo8mTJ2vlypW6cOGC35Nu67333lNhYaGWLVumrKwszZo1Szt37vR7lmvXr1/X22+/rdWrVw/4Lxbur8AH5vLly+rp6dHYsWN73T927Fg1Nzf7tOreEY1GtWHDBhUVFWnmzJl+z7mt06dPa9SoUXIcRy+99JKqqqo0ffp0v2fd1t69e3Xq1CmFw2G/p/TbvHnztHv3bh08eFAVFRU6f/68Hn/88bv+T3Z89tlnqqio0JQpU1RdXa01a9Zo3bp1euutt/ye5sr+/ft19epVvfjii75tGPTfpoyhpbS0VGfOnAnEa+uSNHXqVDU0NKi9vV1/+MMfVFJSopqamrs6Mk1NTVq/fr0OHTqk4cOH+z2n34qLi+P/Oz8/X/PmzVNeXp7eeecd/fCHP/RxWd+i0agKCwu1adMmSdKsWbN05swZbdu2TSUlJT6v679du3apuLhYOTk5vm0I/DOYMWPGKCkpSS0tLb3ub2lp0bhx43xadW9Yu3at3n//fX3wwQfmf4JhoCQnJ+vBBx/U7NmzFQ6HVVBQoNdff93vWX2qr69Xa2urHn30UYVCIYVCIdXU1OiNN95QKBRST0+P3xP75f7779dDDz2kc+fO+T2lT9nZ2Tf9B8fDDz8ciJf3vvb555/r8OHD+tGPfuTrjsAHJjk5WbNnz9aRI0fi90WjUR05ciQwr60HTSwW09q1a1VVVaW//vWvmjRpkt+TPItGo4pEIn7P6NOCBQt0+vRpNTQ0xI/CwkKtXLlSDQ0NSkpK8ntiv3R1denTTz9Vdna231P6VFRUdNO33X/yySfKy8vzaZF7lZWVysrK0pIlS3zdMSReIisrK1NJSYkKCws1d+5cbdmyRd3d3Vq1apXf0/rU1dXV67/mzp8/r4aGBmVkZGjChAk+LutbaWmp9uzZo3fffVcpKSnx97rS0tI0YsQIn9fdWnl5uYqLizVhwgR1dnZqz549Onr0qKqrq/2e1qeUlJSb3t+67777NHr06Lv6fa9XXnlFS5cuVV5eni5duqSNGzcqKSlJK1as8Htan15++WV95zvf0aZNm/T9739fJ06c0I4dO7Rjxw6/p/VLNBpVZWWlSkpKFAr5/CXel+9dM/DrX/86NmHChFhycnJs7ty5sdraWr8n3dYHH3wQk3TTUVJS4ve0Pn3TZkmxyspKv6f1afXq1bG8vLxYcnJyLDMzM7ZgwYLYX/7yF79neRKEb1Nevnx5LDs7O5acnBz71re+FVu+fHns3Llzfs/qlz/96U+xmTNnxhzHiU2bNi22Y8cOvyf1W3V1dUxSrLGx0e8pMX5dPwDARODfgwEA3J0IDADABIEBAJggMAAAEwQGAGCCwAAATBAYAIAJAgMAMEFgAAAmCAwAwASBAQCYIDAAABP/B9tD3+X3oYNZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(dist[3,5,:,:], \"binary\") #即 (3,5) 與各座標的距離\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0c6cd7e5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.+0.j, 1.+0.j, 2.+0.j, 3.+0.j, 4.+0.j, 5.+0.j, 6.+0.j, 7.+0.j,\n",
       "       0.+1.j, 1.+1.j, 2.+1.j, 3.+1.j, 4.+1.j, 5.+1.j, 6.+1.j, 7.+1.j,\n",
       "       0.+2.j, 1.+2.j, 2.+2.j, 3.+2.j, 4.+2.j, 5.+2.j, 6.+2.j, 7.+2.j,\n",
       "       0.+3.j, 1.+3.j, 2.+3.j, 3.+3.j, 4.+3.j, 5.+3.j, 6.+3.j, 7.+3.j,\n",
       "       0.+4.j, 1.+4.j, 2.+4.j, 3.+4.j, 4.+4.j, 5.+4.j, 6.+4.j, 7.+4.j,\n",
       "       0.+5.j, 1.+5.j, 2.+5.j, 3.+5.j, 4.+5.j, 5.+5.j, 6.+5.j, 7.+5.j,\n",
       "       0.+6.j, 1.+6.j, 2.+6.j, 3.+6.j, 4.+6.j, 5.+6.j, 6.+6.j, 7.+6.j,\n",
       "       0.+7.j, 1.+7.j, 2.+7.j, 3.+7.j, 4.+7.j, 5.+7.j, 6.+7.j, 7.+7.j])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.ravel() #此為下圖的x軸和y軸"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "824ccf50",
   "metadata": {},
   "source": [
    "Python optimal transports 的 ot.emd() 及 ot.emd2() 的 M 可以用這種方式生出來。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "894ba1df",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = dist.reshape(64,64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0a43500e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGfCAYAAAAZGgYhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1KklEQVR4nO3df3BV9Z3/8VeQJKCSi/gjgTWwdErFHwUVEbPYXYtpGcY6uDJd7VjLdt06WlABdyrs+KOV1rA6W8U2YnVZtEPdbNkZbO0irsWK0/JLogz+2OX3LinklyG5l4T8zvn+4Xi/hnw+6f0k59zPSXw+Zu6MfO65537O5/54e+555fPJCYIgEAAAWTbCdwcAAJ9NFCAAgBcUIACAFxQgAIAXFCAAgBcUIACAFxQgAIAXFCAAgBcUIACAFxQgAIAXI6PacXl5uZ544gnV1NRo+vTp+slPfqKrr776Tz6up6dHx48f15gxY5STkxNV9wAAEQmCQCdPntSECRM0YkQ/5zlBBCoqKoK8vLzgX//1X4MPPvgg+M53vhOMHTs2qK2t/ZOPraqqCiRx48aNG7chfquqqur3+z4nCMKfjHTWrFmaOXOmfvrTn0r6+KymuLhY99xzj5YvX97vY5PJpMaOHauqqioVFBT0uq+mpsb4mOrq6j5tLtu6bn/8+HHjtuXl5cb2lpYWY3tDQ4Oxvba2NuP+mbaV7H00bf/ggw8at+3s7DS2J5NJY7vteOrr6zPqR3/ttuP/7ne/a2zv6urq02Z7HZqamoztJ06cMLZ/9NFHfdps/TYduyTdeuutxvbu7m5je1tbW5+2kydPGre1HY/t9TG12/p9ww03GNt7enqM7aZ+S9KpU6f6tNneV7bjaWxsNLab+n7dddcZt7WNt+29b+q3ZH5vpVIp47a2ftten5kzZxrbTTo6Oozt7e3txnbbZ8LUd9vrY2uXPn7tEomE9f7Qf4Lr6OhQZWWlVqxYkW4bMWKESktLtX379j7bt7e39xqcTz5UBQUFfQqQbbBMH8QzzzzTuO3o0aON7fn5+cb2vLy8Pm25ubnGbceMGWNst52C2j6cpr6H0W/J3Pezzz7buK3tzez64Rw1alSfNtd+jxxpfqvaXmdTAbJ9Sdo+nGG8J2z9tr2eti9EE9vr4/LaS+Y+nnHGGcZtbWNiG1vb/9+aXh9b/2zttj6aPm+2MbGNt+3nf9t732UMbd8Htue0Hb+J7XWwHadLHwdySeRPPSb0EMJHH32k7u5uFRYW9movLCw0/l9sWVmZEolE+lZcXBx2lwAAMeQ9BbdixQolk8n0raqqyneXAABZEPpPcOedd57OOOOMPr+J19bWqqioqM/2+fn5xtP6mpqaPj+5jR8/PtzOhsz2E+FZZ52V5Z64cf0pZ+zYsRH2xo3ppxzJ/JOI7adGH1x/ErH9ZJdttp94bD8rmX5+9cF1vG3i/lm2/USabUEQ9Htt6BOhnwHl5eVpxowZ2rJlS7qtp6dHW7ZsUUlJSdhPBwAYoiL5O6Bly5Zp4cKFuuqqq3T11VfrqaeeUktLi7797W9H8XQAgCEokgJ0yy23qL6+Xg8//LBqamp0+eWXa/PmzX2CCQCAz67IZkJYvHixFi9eHNXuAQBDnPcUHADgsymyM6DBqq6utv6l9+niko6z/fW4TVwSNba/NLel3eKUjmtubja2mxJvtj8K9ZGOa21tNbbb0m6mtJaPZJztj6dtabe4pONsf0BqE6d0nOsfSg922zD09PT4ScEBAJAJChAAwAsKEADACwoQAMCL2IYQampqrDMdZ8JHMMG2ZICrbIcTTMsL9CdO4QRbgMLEFjbwEU6wTdVvYwoc+Ji2xzbdlI1LOCHKYEJY02S5hBPC+hzbZpk3cQ0bRBVOyHRWd86AAABeUIAAAF5QgAAAXlCAAABeUIAAAF7ENgVXXV0dSZonynRclKu5RpmMO3z4cCj7cUnHhZWMO3jw4KD34SMdd+DAgUHvw2Xanv62dxFGvyVz4i3KaXtc3ydxSsfV1dU5bW+S7XScbaHI03EGBADwggIEAPCCAgQA8IICBADwggIEAPAitim4mpqarC6iFEY6bteuXSH0xE0Y6bht27aF0BM7U+ItrHnjtm/fPoAeZcYlHeeajHOdf89FlOk4l7n3XEW5qF0mi6NlwuXzFtaidmH13SSqRe0yPXbOgAAAXlCAAABeUIAAAF5QgAAAXlCAAABexDYFV11dbU1KZZNLOu7YsWMR9sSNS1onrJVcXYS1qqqPvpsSb67zxkWZgrNxScfZtm1oaAi1T5kIIx134sSJUPt0uijTcVEmD20Gm46zvTZ9tsv4WQAACBEFCADgBQUIAOAFBQgA4EVsQwjHjx9Xbm6u724Y2YIJPi6IuzJdLD1+/LiHnpi5hhPCWKwrDK6L2sWl35I5cGC7UB6nfruEE3yEPsJa1M5HCMEkiqnROAMCAHhBAQIAeEEBAgB4QQECAHhBAQIAeJETBEHguxOflkqllEgktG/fPo0ZM6bXfS0tLcbH1NfX92mzJdKqq6ud2k3T69j2vWnTJqd929JnpnaXbftrNx3P2rVrjdu6jLfkNuYu493fvleuXGls7+jo6NNmSxPZElK2xJepL7b+2V6HJUuWGNu7urqM7c3NzX3aXI/H5XWzHfvtt99ubO/u7ja2t7a2GttTqVSftsbGRuO2tuNxed3mz59v3Lanp8fY3tbWZmy3fSZMr4Xt9bFNZ2Qb89LSUmO7acw7OzuN29r6bXpfSeZF8GzTGZlehyAIlEwmlUwmVVBQYHycxBkQAMATChAAwAsKEADACwoQAMALChAAwIvYzgXX0tLSZ04n17mVss2W7HJZ1M4HW0Im7uMtmdNuknnuONs8cz7Y0m6uC9tlmy3tZpvHzLawXbbZ0m4ui9r54jrmJtn+LPf09BiTdKfjDAgA4AUFCADgBQUIAOAFBQgA4AUFCADghXMK7q233tITTzyhyspKVVdXa+PGjbrpppvS9wdBoEceeUTPP/+8mpqaNHv2bK1Zs0ZTpkxxep6GhgbrfEyni0tay3Vl0bik42xzhNnEZbwl+3xbpsSbbVVVH+k42xxcLiur+kjG2eZ2s6Xd4pKOs32XuKyq2t/2UbLN72bikoyTovss25J7p3M+A2ppadH06dNVXl5uvP/xxx/X008/rWeffVY7d+7UWWedpblz52ZcTAAAnw3OZ0Dz5s3TvHnzjPcFQaCnnnpKDz74YHr22Z///OcqLCzUyy+/rFtvvbXPY9rb29Xe3p7+t2mWXADA8BPqNaAjR46opqam1/ThiURCs2bN0vbt242PKSsrUyKRSN+Ki4vD7BIAIKZCLUCfrCtSWFjYq72wsNC6XsqKFSvS60Ykk0lVVVWF2SUAQEx5n4onPz9f+fn5vrsBAMiyUAtQUVGRJKm2trZXwqu2tlaXX365075qa2t15plnDrgvPpJarik4m2yn42xnp658jLltVUwTW9rNRzrOlt6zMSXefMwb53qN1iUdF2UyzjbfoU2c0nFhzNWY7XScba7D04X6E9zkyZNVVFSkLVu2pNtSqZR27typkpKSMJ8KADDEOZ8BNTc36+DBg+l/HzlyRHv27NG4ceM0ceJELVmyRD/84Q81ZcoUTZ48WQ899JAmTJjQ62+FAABwLkC7d+/Wl7/85fS/ly1bJklauHChXnjhBX3ve99TS0uL7rzzTjU1Nenaa6/V5s2bYzW9OQDAP+cCdN111ykIAuv9OTk5evTRR/Xoo48OqmMAgOHNewrOpqamJpKLklFeKA8rhGASZTAh6uh7lGN++PDhQe/DRzjh0z9jD5TLtD39be/iwIEDg96HZA4cRDltT1j9dgknhPWrj8t7xfWzFtWidplOH8RkpAAALyhAAAAvKEAAAC8oQAAALyhAAAAvYpuCq62tzeoccWEktfbs2TP4jjgKIx23a9euEHriJqxk3LZt20LZj4lLOs41GWebHT4MUabjXKY+chXlonauUx+5MiXewpq2J5lMDqhPn5btdFxHR0dGj+UMCADgBQUIAOAFBQgA4AUFCADgBQUIAOBFbFNwx48ft87FlU0u6ZEo54Jz5ZKOO3bsWIQ9ceOa1glrMT0XpsSb67xxPvrtko6zbRtlCs4mjHRcQ0NDqH3KRFiL2p04cSK0Pp0uqkXt+puw+tM4AwIAeEEBAgB4QQECAHhBAQIAeEEBAgB4EdsUXG1trXJzc313w8iWHIlTmszGlI7zkchyZRvzuCQPXVdVrauri7A3bkyJN9u8cXHqt0s6Lk79dk3H+UgeDjYdl2lijjMgAIAXFCAAgBcUIACAFxQgAIAXOUGmcyZkSSqVUiKR0Jtvvtnn4qhtkSPTYlO2C3e2i5G2C/GmdtuF7/Xr1xvbW1pajO319fUZP2d1dbVxW1u7LRBh2vemTZuc9m07fpd2133Yjmft2rXGdtOYu4y35DbmLuMtSStXrjS2x/09vmTJEmN7V1eXsb25udnY7nI8rq+b6fhvv/1247bd3d3G9tbWVmN7KpUytjc2NvZpsx2P6+s2f/58Y3tPT0+ftra2NuO2tu8g20J9pnbbdEamfvf09Ojo0aNKJpMqKCgwPk7iDAgA4AkFCADgBQUIAOAFBQgA4AUFCADgRWyn4uns7OyTCHJd9CvbbEkT10XWss2W9nJZ1M6XoTrmtrRb3N/jtrSbbeoe28J22WZLu7ksaueLKe0mmafusU3zk23d3d06evTon9yOMyAAgBcUIACAFxQgAIAXFCAAgBcUIACAF7FNwSWTSXV2dvZqc1n0y0dqyDZnlU1cklqui7rFKR3nMuZxGW/JPgdX3N/jtrndbGm3uKTjbHO7uSxq19/2UbLN72ZKvNkWtct2Ou70724bzoAAAF5QgAAAXlCAAABeUIAAAF5QgAAAXsQ2BdfQ0KBTp05ltK0pDeRjTi3bCo2usp3Wck3B2fhIx4Ux5j7ScbZVMW3i8h63pfdsXNJxUSbjbCuZ2sQpHWeb79DElnbLdjqOFBwAINYoQAAALyhAAAAvKEAAAC+cClBZWZlmzpypMWPG6IILLtBNN92kffv29dqmra1NixYt0rnnnquzzz5bCxYsUG1tbaidBgAMfU4puK1bt2rRokWaOXOmurq69I//+I/66le/qg8//DCdJFq6dKn+8z//Uxs2bFAikdDixYt188036w9/+INTx+rr6weV0HCZU6u/7V3YVhYNQ5RJrbBScDZRpuOqqqoi23eUY3748OFB78PHe/zgwYOD3odkTrxFOW/cgQMHBr0PyS0dF1YyLoy+ZzsdZ1vx93ROBWjz5s29/v3CCy/oggsuUGVlpf7yL/9SyWRSa9eu1UsvvaQ5c+ZIktatW6eLL75YO3bs0DXXXOPydACAYWxQ14CSyaQkady4cZKkyspKdXZ2qrS0NL3N1KlTNXHiRG3fvt24j/b2dqVSqV43AMDwN+AC1NPToyVLlmj27Nm67LLLJH38R4F5eXl9TvULCwutfzBYVlamRCKRvhUXFw+0SwCAIWTABWjRokV6//33VVFRMagOrFixQslkMn2L8jd9AEB8DGgqnsWLF+s3v/mN3nrrLV144YXp9qKiInV0dKipqanXWVBtba2KioqM+8rPz1d+fn6f9traWmP7YEV54Xb//v0D6NHghHGhfM+ePYPviKOwggm7du0KZT8uwhjzbdu2hdATsyjf47af0sMQ5aJ2rlMfuTIFDsKatsd1+iMXLuEEl2CC7TXr8zwZ71FSEARavHixNm7cqDfeeEOTJ0/udf+MGTOUm5urLVu2pNv27duno0ePqqSkxOWpAADDnNMZ0KJFi/TSSy/pV7/6lcaMGZO+rpNIJDR69GglEgndcccdWrZsmcaNG6eCggLdc889KikpIQEHAOjFqQCtWbNGknTdddf1al+3bp3+9m//VpL05JNPasSIEVqwYIHa29s1d+5cPfPMM6F0FgAwfDgVoCAI/uQ2o0aNUnl5ucrLywfcKQDA8MdccAAAL2K7IF1tba01tRMFl+SQbdtjx46F2KPBcUlqRT0VjwvXdFxcxtw1GRfW4oUuwniP++h3GOm4qFNwJmEtatfQ0BBanzJlSry5TNuTk5OT0fNwBgQA8IICBADwggIEAPCCAgQA8IICBADwIrYpuJqamoznE4qSKQ1kS+f5SAi5MqW14pIk648tHRf3Mbel4+KUPHR5j9fV1UXcm8y5pOPi1G/XdFxc+u4yb5wtMddnu0H1CACAAaIAAQC8oAABALygAAEAvKAAAQC8yAkymeI6i1KplBKJhCoqKnTmmWf2uq+rq8v4mObm5j5ttlUEbXNC1dfXG9tNKStbKmX16tXG9o6ODmO7Sx9tz2lLgdnaTemr9evXG7dtaWkxtruMlSRVV1dn1CbZE3m2fW/atCnj57Qlz8Jot21rO561a9ca213G3GW8+2s39dG275UrVxrb4/4eX7JkiXFbl+8Uye14XD8ntuO//fbbje3d3d192lpbW43bplIpY3tjY6Ox3XQ8tu9OU7+7u7v1/vvvK5lMqqCgwPg4iTMgAIAnFCAAgBcUIACAFxQgAIAXFCAAgBf+J1uz6Orq6pNQcVkB0QdbEsg2r5Zt1clssyWvXFf59MGW7HJdWTXbhuqYD9X3uC3tFvfvFMmcdpPMc8fZ5pnLNtt4n44zIACAFxQgAIAXFCAAgBcUIACAF7ENIbS0tKinp6dXm8sCVD4uItqm6bBdiI3LhVvblCE2cbpQ7rKwW5yCCUN1zIfqe9w2tY7Ld0p/20fJNr2OKXBgW9Qu2+GEzs7OjLbjDAgA4AUFCADgBQUIAOAFBQgA4AUFCADgRWxTcE1NTWpvb89oW1MyxUeKxbZgk41LcijK1JBtgSxXPpJaLik4Gx/puDDG3Md4D9X3uC29ZxOndJxtMTkTW9ot2+k425RNp+MMCADgBQUIAOAFBQgA4AUFCADgBQUIAOBFbFNwJ06cUH5+/oAf7yPFUldXN+h9SOY0UJRzatkWdQtLlGmtMFJwNlGm46qqqiLbd5Tjffjw4VD2k+33+MGDBwe9D8nPfJQHDhwY9D6ynY7LNMHMGRAAwAsKEADACwoQAMALChAAwAsKEADAi9im4D766CNrKmYwokzH7d27d0B9ykSUK07u379/AD0anLCSWnv27AllPy7CSMft2rUrhJ64CWPMt23bFkJPzKJ8j2/fvn0APcpclPNRus6/58IlHeeSjGtra8toO86AAABeUIAAAF5QgAAAXlCAAABeOIUQ1qxZozVr1uh///d/JUmXXnqpHn74Yc2bN0/Sxxee7r//flVUVKi9vV1z587VM888o8LCQueO1dbWKjc31/lxAxXGFBthLezmIowLt8eOHQuxR4PjeqE8yql4XLgGE4bqmA/V97iPfocVeIoyhGBjChxEMW2P0xnQhRdeqFWrVqmyslK7d+/WnDlzNH/+fH3wwQeSpKVLl+qVV17Rhg0btHXrVh0/flw333zzgDsHABi+nM6Abrzxxl7//tGPfqQ1a9Zox44duvDCC7V27Vq99NJLmjNnjiRp3bp1uvjii7Vjxw5dc8014fUaADDkDfgaUHd3tyoqKtTS0qKSkhJVVlaqs7NTpaWl6W2mTp2qiRMn9pvBb29vVyqV6nUDAAx/zgXovffe09lnn638/Hzddddd2rhxoy655BLV1NQoLy+vz++vhYWF/f7+WlZWpkQikb4VFxc7HwQAYOhxLkAXXXSR9uzZo507d+ruu+/WwoUL9eGHHw64AytWrFAymUzfolwnBQAQH85T8eTl5enzn/+8JGnGjBl6++23tXr1at1yyy3q6OhQU1NTr7Og2tpaFRUVWfeXn59vXHiuvr7emhbJJpcpNuKSyJLckkM+EkKubEmtOKXJTGzpuLiPuW28h+p7PKzFIsPgmo6LS99dpu2xJeZON+i/A+rp6VF7e7tmzJih3NxcbdmyJX3fvn37dPToUZWUlAz2aQAAw4zTKcaKFSs0b948TZw4USdPntRLL72kN998U6+99poSiYTuuOMOLVu2TOPGjVNBQYHuuecelZSUkIADAPThVIDq6ur0rW99S9XV1UokEpo2bZpee+01feUrX5EkPfnkkxoxYoQWLFjQ6w9RAQA4nVMBWrt2bb/3jxo1SuXl5SovLx9UpwAAwx9zwQEAvMgJgiDw3YlPS6VSSiQSevLJJ/ukLrq7u42PaW1tNe7HpLGx0dhum2/J1G5LpaxcudLY3tXVZWxvbm42tjc1NWXcv/r6emO7LWVl6vvq1auN23Z0dGTcP8neR9Nz2vpna7elr9avX29sb2lp6dPmOlbV1dUZt9vSeLZ9b9q0yek5TcdvG5Mw2m3HY/sVxDTektuYu4y35Dbmts/mUHiPL1myxNhu+l5x+U6R3L5XXL5Turq6tHPnTiWTSRUUFBgfJ3EGBADwhAIEAPCCAgQA8IICBADwggIEAPDC/2RrFt3d3X1Sb1GsyBcmW9rNdQXEbLMlgVxWnPTFlr5yXVk122zJLteVVbNtqI73UH6Pu3yvxOU7pbOzM6PtOAMCAHhBAQIAeEEBAgB4QQECAHhBAQIAeBHbFFxbW1ufNpcV+Xwk42zzMLmugJjtJIttniiXFSf72z5KtrnGTOKU1HJdWTQu6TiX8ZbiM+ZD+T3u8r0Sl+8UW+rwdJwBAQC8oAABALygAAEAvKAAAQC8iG0I4eTJkxlfyDIFDnxM22O70GnjEk6I8iKibVEqmzhduLUtkuXCx4Vy1xCCiY9gQhjjLWV/zIfye9zleyUugaf29vaMtuMMCADgBQUIAOAFBQgA4AUFCADgBQUIAOBFbFNwTU1N1sRJJlym7elvexeuSRubbE+xUVdXN+h9SG7JobBSQ7aF3cIQZVIrjBScTZTpuKqqqsj2HeV4Hz58OJT9+HiPHzx4cND7yHY6rrW1NaPtOAMCAHhBAQIAeEEBAgB4QQECAHhBAQIAeBHbFFxDQ4Nyc3ND32+U6bhDhw4NqE+ZiDLFsnfv3gH1KVOmNFBYc2rt379/AD0anDDSWnv27Bl8RxyFkY7btWtXCD1xE8Z4b9u2LYSe2EX5Ht++ffsAepSZqOajzM/Pz2g7zoAAAF5QgAAAXlCAAABeUIAAAF5QgAAAXsQ6BWdLeEXBJR1n2zas1SJdhJGO89HvsFacPHbsWEg9GhzXpFaUc8G5cknHxWW8JbcxH8rv8bh8r7h8p4wYkdm5DWdAAAAvKEAAAC8oQAAALyhAAAAvYhtCqK+vt06Pk02mwIGtX2Et7BYGl3BCnC6Iu1649XGB1oXtQnmcLuab2IIJcR9vyTzmQ/k9HpfvFZfvlEwDZJwBAQC8oAABALygAAEAvKAAAQC8oAABALzICYIgGOiDV61apRUrVui+++7TU089JUlqa2vT/fffr4qKCrW3t2vu3Ll65plnVFhYmNE+U6mUEomEHnjggT6LGvX09Bgf09bW1qetpaXFuG1TU5NTe0NDQ582Wypl+fLlxvbu7m5je2trq7E9lUr1aWtsbDRu+9FHHzm1m/q+cuVK47ZdXV3G9ubmZmO7bQxNfamvrzdua0tZ2cZ89erVxvaOjo5B9a+/5zT10dZvW/pq/fr1xnbb+9Y0XrbnrK6udmo3JfJs+960aZPTvm3Hb2p32ba/dtPxrF271rity3hLbmPuMt797dv2+Yzze7yrq0tvvfWWksmkCgoKjI+TBnEG9Pbbb+tnP/uZpk2b1qt96dKleuWVV7RhwwZt3bpVx48f18033zzQpwEADFMDKkDNzc267bbb9Pzzz+ucc85JtyeTSa1du1Y//vGPNWfOHM2YMUPr1q3Ttm3btGPHjtA6DQAY+gZUgBYtWqQbbrhBpaWlvdorKyvV2dnZq33q1KmaOHGidV3z9vZ2pVKpXjcAwPDnPBNCRUWF3nnnHb399tt97qupqVFeXl6fv/QtLCy0/n5YVlamH/zgB67dAAAMcU5nQFVVVbrvvvv0i1/8QqNGjQqlAytWrFAymUzfqqqqQtkvACDenM6AKisrVVdXpyuvvDLd1t3drbfeeks//elP9dprr6mjo0NNTU29zoJqa2tVVFRk3Gd+fn6ftJv0ceLt9NSbbZGjsIrhYNnSbra542wL22WbLe3msgCVL6YkkGSeV8s2B5cPtvSV68J22WZLdrksaufDUB1vaWi+x219Pp1TAbr++uv13nvv9Wr79re/ralTp+qBBx5QcXGxcnNztWXLFi1YsECStG/fPh09elQlJSUuTwUAGOacCtCYMWN02WWX9Wo766yzdO6556bb77jjDi1btkzjxo1TQUGB7rnnHpWUlOiaa64Jr9cAgCEv9OUYnnzySY0YMUILFizo9YeoAAB82qAL0Jtvvtnr36NGjVJ5ebnKy8sHu2sAwDDGXHAAAC9iuyJqW1ubTp+mzpZ2M6XjfCTjbHO72dJucUnH2eZ2c1kBsb/to2Sb+8qUBrKtOOkjOWSba8wmLmkt15VF45KOG6rjLQ3N93h7e3tG23EGBADwggIEAPCCAgQA8IICBADwggIEAPAitim4U6dOWecoO50p8eZj3jjXpSRc0nFRJuNsKRubOKXjbCs9mtiSQD6SQ7bZ4V34SGq5puBssp2OC2O8JT9jPhTf46ZVqk04AwIAeEEBAgB4QQECAHhBAQIAeBHbEEIymVRubu6AH+8ybU9/27tobGwc9D4kc+Agyml7XC5y9sclnBBWMKGurm7Q+/Bx4da2sFsYorxQHlYIwSTKYELUKy1HOeaHDx8e9D6y/R63TUt2Os6AAABeUIAAAF5QgAAAXlCAAABeUIAAAF7ENgXX1NQ0qBScTZTpuGPHjg2oT5mIclG7Q4cODahPmTIl3sKatmfv3r0D6lMmXJJDrqmh/fv3D6BHgxNGUmvPnj2D74ijMNJxu3btCqEnbsJKxm3bti2U/ZhE9R4/depURttxBgQA8IICBADwggIEAPCCAgQA8IICBADwIrYpuMbGRmvCKwou6TjbtmHNqeYijHRcWIt1uQhrUTsffTelgVzn1IoyMenKJa0V5VxwrlzScUN1vKWh+R7PNMHMGRAAwAsKEADACwoQAMALChAAwAsKEADAi9im4Orr663zs2WTKfFm61cYq3OGxSUdF6d+u6bj4pLKcl1x0keyyYUtqRWnNJmNKR0X9/GW7GM+FN/jHR0dGe3T/zc8AOAziQIEAPCCAgQA8IICBADwIicIgsB3Jz4tlUopkUjoW9/6Vp+LW93d3cbHdHZ29mlraWkxbtvc3GxsTyaTxvYTJ070abNNuXPvvfca23t6eoztbW1txnZT35uamozb2tobGhqM7abAwfLly43b2sa7tbXV2J5KpYztjY2NfdpsY2hrtwUlVq5caWzv6urq02Z77W1jaOtLfX19nzbbRW5bv1evXm1st128NfXRdaxsfTS12y58r1+/3thu+7yZxsr2nNXV1cZtbe22QIRp35s2bXLat+34Xdpd92E7nrVr1xrbTWPuMt6S25i7jHdnZ6d++9vfKplMqqCgwPg4iTMgAIAnFCAAgBcUIACAFxQgAIAXFCAAgBexnYqnu7u7TwrLZYE610WfwmBLu9mm7rEtbJdttrSby6J2vpjSbpJ56h7bND8+2NJurgvbZZst7ebj8+bClvZyWdTOl6E45kzFAwCINQoQAMALChAAwAsKEADACwoQAMALpxTc97//ff3gBz/o1XbRRRfpf/7nfyR9PLfZ/fffr4qKCrW3t2vu3Ll65plnVFhY6Nyxzs5O5eTkZLRtXNJxtrndbGm3uKTjbHO7uSxq19/2UbLN72ZKvNkWtfORjrPNP+ey6JePZJxtrjGbuCS1XBd1i1M6zmXM4zLe7e3tGW3nfAZ06aWXqrq6On37/e9/n75v6dKleuWVV7RhwwZt3bpVx48f18033+z6FACAzwDnvwMaOXKkioqK+rQnk0mtXbtWL730kubMmSNJWrdunS6++GLt2LFD11xzjXF/7e3tvaqlbUZlAMDw4nwGdODAAU2YMEGf+9zndNttt+no0aOSpMrKSnV2dqq0tDS97dSpUzVx4kRt377dur+ysjIlEon0rbi4eACHAQAYapwK0KxZs/TCCy9o8+bNWrNmjY4cOaIvfelLOnnypGpqapSXl9fnd+nCwkLrWhSStGLFCiWTyfStqqpqQAcCABhanH6CmzdvXvq/p02bplmzZmnSpEn65S9/OeAL0Pn5+crPzx/QYwEAQ9eg5oIbO3asvvCFL+jgwYP6yle+oo6ODjU1NfU6C6qtrTVeM/pTTp061WelU5eEh0syznXfNrY5m2xc0nFRJuNcr7vFKR1nS5OZ2NJuPtJxttVMbUyJNx/zxvX3a4aLbKe1XFNwNj7ScWGMebbH25asPd2g/g6oublZhw4d0vjx4zVjxgzl5uZqy5Yt6fv37duno0ePqqSkZDBPAwAYhpzOgP7hH/5BN954oyZNmqTjx4/rkUce0RlnnKFvfOMbSiQSuuOOO7Rs2TKNGzdOBQUFuueee1RSUmJNwAEAPrucCtAf//hHfeMb31BDQ4POP/98XXvttdqxY4fOP/98SdKTTz6pESNGaMGCBb3+EBUAgNM5FaCKiop+7x81apTKy8tVXl4+qE4BAIY/5oIDAHgR2xVRW1parOmk08UlHeeSyOqPKfEW5bxxjY2Ng96H5JaOCysZ55omM/GRjqurqxv0Plzmjetvexe2lUXDEGVSK6wUnE2U6bgo/zYyqjE/depURttxBgQA8IICBADwggIEAPCCAgQA8CK2IYRUKuUcGPg014trYSxq19DQ4PScLqJc1O7YsWMD6lOmTIGDsKbtOXTo0ID6lAmXcIJrMGHv3r0D6lMmogwn7N+/fwA9GpwwLpTv2bNn8B1xFFYwYdeuXaHsx8VgxzzTack4AwIAeEEBAgB4QQECAHhBAQIAeEEBAgB4EdsUXGNjozXhNRhRpuOiTMHZhJGOC2M6G1dhLWoX1gJpLkyJN9dpe3z02yUdZ9s26sSkC5fPctRT8bhwTcfFZcxdxjvTFC5nQAAALyhAAAAvKEAAAC8oQAAALyhAAAAvYpuCa2hoUE5OTtaeL4xF7cJYZCwsLum4OPXbNR0Xl767LmoXp1SWKfFmmzfOR3rPlemzHJckWX9s6bi4j7lpvHt6ejJ6LGdAAAAvKEAAAC8oQAAALyhAAAAvKEAAAC9ygiAIfHfi01KplBKJhEpLS5Wbm5vRY9rb2zNqk6Tm5mZjezKZNLY3NTVl1CZJf//3f29s7+7uNrZ3dnYa202rCbr2+8SJE8Z207xv9957r3FbW5Klra3N2G5bBdFlDG3z6dnSbsuXLze2m8a8tbXVuG0qlTK2NzY2GttNY2ibT8/W75UrVxrbu7q6jO2m1982hra+1NfXG9tNKStbv1evXm1s7+joMLa79NH2nLYUmK3dlDBcv369cVvbe9ZlrCSpuro6ozbJnsiz7XvTpk0ZP6ctXRlGu21b0/F0dnbqt7/9rZLJpAoKCoyPkzgDAgB4QgECAHhBAQIAeEEBAgB4QQECAHgR27ngXOTn5/vugiR72s1lVVXJfdXWwbKl3VxWVfXFZcxt88z5YEu7ua6smm22tJtt7jjbyqrZZku7ZfuzNhC2NJ3ryqrZZHufnI4zIACAFxQgAIAXFCAAgBcUIACAF7ENIXR0dPS5OO4SNvARTLBNrWPjEk6I8mKpbWodl0Xt+ts+SrbpdUyBA9t4+wgn2KZWclnYzkcwwTa1ji1sEJdwgm1qHZs4hRNcFi+MSzDBNhXa6TgDAgB4QQECAHhBAQIAeEEBAgB4QQECAHgR2xRce3u7dZqV08UlHRfWdB/ZTsfZ+m0Tp3ScbTE5E1vazUc6zpYmszEl3nxM22Nb7M7GJR0XZTLOttibKx/pOJcUnE2203G2dOrpOAMCAHhBAQIAeEEBAgB4QQECAHjhXICOHTumb37zmzr33HM1evRoffGLX9Tu3bvT9wdBoIcffljjx4/X6NGjVVpaqgMHDoTaaQDA0OeUgmtsbNTs2bP15S9/Wa+++qrOP/98HThwQOecc056m8cff1xPP/20XnzxRU2ePFkPPfSQ5s6dqw8//NApDdXS0uK8kNunuabdwkjH2eb3snFJ1ES5qJ1rIsvGJR0XVjKusbFx0PvwkY5zTZOZuMwb19/2Lurq6ga9D8mceIty3jjbom5hiTIdF0YKziaqdNypU6cy2s6pAP3TP/2TiouLtW7dunTb5MmT0/8dBIGeeuopPfjgg5o/f74k6ec//7kKCwv18ssv69Zbb3V5OgDAMOb0E9yvf/1rXXXVVfr617+uCy64QFdccYWef/759P1HjhxRTU2NSktL022JREKzZs3S9u3bjftsb29XKpXqdQMADH9OBejw4cNas2aNpkyZotdee01333237r33Xr344ouS/v8fexUWFvZ6XGFhofUPwcrKypRIJNK34uLigRwHAGCIcSpAPT09uvLKK/XYY4/piiuu0J133qnvfOc7evbZZwfcgRUrViiZTKZvVVVVA94XAGDocCpA48eP1yWXXNKr7eKLL9bRo0clSUVFRZKk2traXtvU1tam7ztdfn6+CgoKet0AAMOfUwhh9uzZ2rdvX6+2/fv3a9KkSZI+DiQUFRVpy5YtuvzyyyV9PFfXzp07dffddzt1LJVKWecWG4wo03HJZNK1O0bZTsc1NDQ47cOVKfEW1rxxx44dG1CfMuGSjnNNxh06dGhAfcpElOm4vXv3DqhPmYhyVdX9+/cPoEeDE1Yybs+ePaHsx8Vg03GZJoKdCtDSpUv1F3/xF3rsscf0N3/zN9q1a5eee+45Pffcc5KknJwcLVmyRD/84Q81ZcqUdAx7woQJuummm5wPAgAwfDkVoJkzZ2rjxo1asWKFHn30UU2ePFlPPfWUbrvttvQ23/ve99TS0qI777xTTU1Nuvbaa7V58+ZIZ0QGAAw9zssxfO1rX9PXvvY16/05OTl69NFH9eijjw6qYwCA4Y254AAAXsR2QbpkMqmcnJysPV8Yi9p9kgaMQpSL2kUdQjAJa1G7MKa0cWUKHLhO2xPWAmkuXMIJtm199DuMcEKUYRVXrp/lKKficeESTDh58mRG23EGBADwggIEAPCCAgQA8IICBADwggIEAPAi1im4OHBJx/lIZIWRjgtrkbEwuKbj4tJ310Xt4tJvyZx4s03bE5dEluSWjvOR3nNl+yzHKcFnYkrHZfq9xBkQAMALChAAwAsKEADACwoQAMCL2IUQgiBwvq+np6dPW3d3t3Hbrq4uY7vtYrHLmkQu/ZPc+tjZ2WnctqOjw6kvpuN07Z+tLy59tF3ktk2/ZHsdohzD9vZ2Y3tbW5ux3cT2vrK9D136aOtfa2ursd0WqDGNre31sfU7jDG09fvUqVPG9tzc3Iz7EtZr79LHlpYW47a2oI3tcxjle8I2tqb1fGzT65gCB6lUSlL/3+eSlBP8qS2y7I9//KOKi4t9dwMAMEhVVVW68MILrffHrgD19PTo+PHjGjNmjE6ePKni4mJVVVUN66W6U6kUxzlMfBaOUeI4h5uwjzMIAp08eVITJkzo91ek2P0EN2LEiHTF/OTnmIKCgmH94n+C4xw+PgvHKHGcw02Yx5lIJP7kNoQQAABeUIAAAF7EugDl5+frkUcecZoOZyjiOIePz8IxShzncOPrOGMXQgAAfDbE+gwIADB8UYAAAF5QgAAAXlCAAABeUIAAAF7EugCVl5frz//8zzVq1CjNmjVLu3bt8t2lQXnrrbd04403asKECcrJydHLL7/c6/4gCPTwww9r/PjxGj16tEpLS3XgwAE/nR2gsrIyzZw5U2PGjNEFF1ygm266Sfv27eu1TVtbmxYtWqRzzz1XZ599thYsWKDa2lpPPR6YNWvWaNq0aem/HC8pKdGrr76avn84HOPpVq1apZycHC1ZsiTdNhyO8/vf/75ycnJ63aZOnZq+fzgc4yeOHTumb37zmzr33HM1evRoffGLX9Tu3bvT92f7Oyi2Bejf//3ftWzZMj3yyCN65513NH36dM2dOzdWSxm7amlp0fTp01VeXm68//HHH9fTTz+tZ599Vjt37tRZZ52luXPnOs2+7NvWrVu1aNEi7dixQ6+//ro6Ozv11a9+tdfMwEuXLtUrr7yiDRs2aOvWrTp+/Lhuvvlmj712d+GFF2rVqlWqrKzU7t27NWfOHM2fP18ffPCBpOFxjJ/29ttv62c/+5mmTZvWq324HOell16q6urq9O33v/99+r7hcoyNjY2aPXu2cnNz9eqrr+rDDz/UP//zP+ucc85Jb5P176Agpq6++upg0aJF6X93d3cHEyZMCMrKyjz2KjySgo0bN6b/3dPTExQVFQVPPPFEuq2pqSnIz88P/u3f/s1DD8NRV1cXSAq2bt0aBMHHx5Sbmxts2LAhvc1///d/B5KC7du3++pmKM4555zgX/7lX4bdMZ48eTKYMmVK8Prrrwd/9Vd/Fdx3331BEAyf1/KRRx4Jpk+fbrxvuBxjEATBAw88EFx77bXW+318B8XyDKijo0OVlZUqLS1Nt40YMUKlpaXavn27x55F58iRI6qpqel1zIlEQrNmzRrSx5xMJiVJ48aNkyRVVlaqs7Oz13FOnTpVEydOHLLH2d3drYqKCrW0tKikpGTYHeOiRYt0ww039DoeaXi9lgcOHNCECRP0uc99TrfddpuOHj0qaXgd469//WtdddVV+vrXv64LLrhAV1xxhZ5//vn0/T6+g2JZgD766CN1d3ersLCwV3thYaFqamo89SpanxzXcDrmnp4eLVmyRLNnz9Zll10m6ePjzMvL09ixY3ttOxSP87333tPZZ5+t/Px83XXXXdq4caMuueSSYXWMFRUVeuedd1RWVtbnvuFynLNmzdILL7ygzZs3a82aNTpy5Ii+9KUv6eTJk8PmGCXp8OHDWrNmjaZMmaLXXntNd999t+699169+OKLkvx8B8VuOQYMH4sWLdL777/f6/f04eSiiy7Snj17lEwm9R//8R9auHChtm7d6rtboamqqtJ9992n119/3bqK53Awb9689H9PmzZNs2bN0qRJk/TLX/5So0eP9tizcPX09Oiqq67SY489Jkm64oor9P777+vZZ5/VwoULvfQplmdA5513ns4444w+SZPa2loVFRV56lW0Pjmu4XLMixcv1m9+8xv97ne/67UiYlFRkTo6OtTU1NRr+6F4nHl5efr85z+vGTNmqKysTNOnT9fq1auHzTFWVlaqrq5OV155pUaOHKmRI0dq69atevrppzVy5EgVFhYOi+M83dixY/WFL3xBBw8eHDavpSSNHz9el1xySa+2iy++OP1zo4/voFgWoLy8PM2YMUNbtmxJt/X09GjLli0qKSnx2LPoTJ48WUVFRb2OOZVKaefOnUPqmIMg0OLFi7Vx40a98cYbmjx5cq/7Z8yYodzc3F7HuW/fPh09enRIHadJT0+P2tvbh80xXn/99Xrvvfe0Z8+e9O2qq67Sbbfdlv7v4XCcp2tubtahQ4c0fvz4YfNaStLs2bP7/EnE/v37NWnSJEmevoMiiTaEoKKiIsjPzw9eeOGF4MMPPwzuvPPOYOzYsUFNTY3vrg3YyZMng3fffTd49913A0nBj3/84+Ddd98N/u///i8IgiBYtWpVMHbs2OBXv/pVsHfv3mD+/PnB5MmTg9bWVs89z9zdd98dJBKJ4M033wyqq6vTt1OnTqW3ueuuu4KJEycGb7zxRrB79+6gpKQkKCkp8dhrd8uXLw+2bt0aHDlyJNi7d2+wfPnyICcnJ/iv//qvIAiGxzGafDoFFwTD4zjvv//+4M033wyOHDkS/OEPfwhKS0uD8847L6irqwuCYHgcYxAEwa5du4KRI0cGP/rRj4IDBw4Ev/jFL4IzzzwzWL9+fXqbbH8HxbYABUEQ/OQnPwkmTpwY5OXlBVdffXWwY8cO310alN/97neBpD63hQsXBkHwcQzyoYceCgoLC4P8/Pzg+uuvD/bt2+e3045MxycpWLduXXqb1tbW4Lvf/W5wzjnnBGeeeWbw13/910F1dbW/Tg/A3/3d3wWTJk0K8vLygvPPPz+4/vrr08UnCIbHMZqcXoCGw3Hecsstwfjx44O8vLzgz/7sz4JbbrklOHjwYPr+4XCMn3jllVeCyy67LMjPzw+mTp0aPPfcc73uz/Z3EOsBAQC8iOU1IADA8EcBAgB4QQECAHhBAQIAeEEBAgB4QQECAHhBAQIAeEEBAgB4QQECAHhBAQIAeEEBAgB48f8AiZAcaVnptiwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(dist.reshape(64,64),\"binary\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ddc2d90",
   "metadata": {},
   "source": [
    "# 實驗一下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "05ead127",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.zeros((8,8))\n",
    "B = np.zeros((8,8))\n",
    "\n",
    "A[2,3]=1\n",
    "B[5,7]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "554e7dc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "07ccf212",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ot.emd2(A.ravel(),B.ravel(),M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "2e730f0d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fc8d3e63b80>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGfCAYAAAAZGgYhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdBElEQVR4nO3df2yV9d3/8Vdr20MFekornLajZTWiBRGGBcoZuDno1nEbA6M6NJgxRySyggJb1C4KbnGWaRTE8UOdA81kTJYA4m5gpkqNriBUiSCzgjZrJ5yDLvac0tlDoZ/vH96er0dAPeXgmx6ej+RK6HVd5+r7E5LzzNVzTpvinHMCAOBrlmo9AADg/ESAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJtLO1oWXL1+uBx98UIFAQCNGjNCjjz6qMWPGfOnjurq6dOjQIfXt21cpKSlnazwAwFninFNbW5sKCgqUmvoF9znuLFi3bp3LyMhwf/zjH91bb73lbrnlFpedne2CweCXPralpcVJYmNjY2Pr4VtLS8sXPt+nOJf4X0ZaVlam0aNH6/e//72kT+5qCgsLNXfuXN11111f+NhQKKTs7GyN1/8oTemJHg0AcJYdV6de0f+qtbVVXq/3tOcl/Edwx44dU0NDg6qrq6P7UlNTVV5ervr6+pPOj0QiikQi0a/b2tr+b7B0paUQIADocf7vtubLXkZJ+JsQPvzwQ504cUI+ny9mv8/nUyAQOOn8mpoaeb3e6FZYWJjokQAA5yDzd8FVV1crFApFt5aWFuuRAABfg4T/CO6iiy7SBRdcoGAwGLM/GAwqLy/vpPM9Ho88Hk+ixwAAnOMSfgeUkZGh0tJS1dbWRvd1dXWptrZWfr8/0d8OANBDnZXPAS1YsEAzZszQqFGjNGbMGC1dulTt7e26+eabz8a3AwD0QGclQNOmTdMHH3yghQsXKhAI6Fvf+pa2bt160hsTAADnr7PyOaAzEQ6H5fV6dbUm8zZsAOiBjrtObdcmhUIhZWVlnfY883fBAQDOTwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADARNwBevnll3XttdeqoKBAKSkp2rhxY8xx55wWLlyo/Px8ZWZmqry8XAcOHEjUvACAJBF3gNrb2zVixAgtX778lMcfeOABLVu2TKtWrdLOnTvVu3dvVVRUqKOj44yHBQAkj7R4HzBp0iRNmjTplMecc1q6dKnuvvtuTZ48WZL09NNPy+fzaePGjbrhhhtOekwkElEkEol+HQ6H4x0JANADJfQ1oKamJgUCAZWXl0f3eb1elZWVqb6+/pSPqampkdfrjW6FhYWJHAkAcI5KaIACgYAkyefzxez3+XzRY59XXV2tUCgU3VpaWhI5EgDgHBX3j+ASzePxyOPxWI8BAPiaJfQOKC8vT5IUDAZj9geDwegxAACkBAeouLhYeXl5qq2tje4Lh8PauXOn/H5/Ir8VAKCHi/tHcEePHtXBgwejXzc1NWnPnj3KyclRUVGR5s2bp/vuu0+DBw9WcXGx7rnnHhUUFGjKlCmJnBsA0MPFHaDdu3fre9/7XvTrBQsWSJJmzJihNWvW6I477lB7e7tmzZql1tZWjR8/Xlu3blWvXr0SNzUAoMdLcc456yE+KxwOy+v16mpNVlpKuvU4AIA4HXed2q5NCoVCysrKOu15/C44AIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgIk06wEAIJG2Hdpz0r6Kgm997XPgy3EHBAAwQYAAACYIEADABAECAJiIK0A1NTUaPXq0+vbtqwEDBmjKlClqbGyMOaejo0NVVVXKzc1Vnz59VFlZqWAwmNChAQA9X4pzzn3Vk3/4wx/qhhtu0OjRo3X8+HH96le/0r59+7R//3717t1bkjR79mz97W9/05o1a+T1ejVnzhylpqbq1Vdf/UrfIxwOy+v16mpNVlpKevdWBQAwc9x1ars2KRQKKSsr67TnxRWgz/vggw80YMAA1dXV6Tvf+Y5CoZD69++vtWvX6rrrrpMkvf322xoyZIjq6+s1duzYL70mAQKAnu2rBuiMXgMKhUKSpJycHElSQ0ODOjs7VV5eHj2npKRERUVFqq+vP+U1IpGIwuFwzAYASH7dDlBXV5fmzZuncePGadiwYZKkQCCgjIwMZWdnx5zr8/kUCAROeZ2amhp5vd7oVlhY2N2RAAA9SLcDVFVVpX379mndunVnNEB1dbVCoVB0a2lpOaPrAQB6hm79Kp45c+bo+eef18svv6yBAwdG9+fl5enYsWNqbW2NuQsKBoPKy8s75bU8Ho88Hk93xgAA9GBx3QE55zRnzhxt2LBBL774ooqLi2OOl5aWKj09XbW1tdF9jY2Nam5ult/vT8zEAICkENcdUFVVldauXatNmzapb9++0dd1vF6vMjMz5fV6NXPmTC1YsEA5OTnKysrS3Llz5ff7v9I74AAA54+4ArRy5UpJ0tVXXx2zf/Xq1frpT38qSVqyZIlSU1NVWVmpSCSiiooKrVixIiHDAgCSxxl9Duhs4HNAANCzfS2fAwIAoLsIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgIm4ArRy5UoNHz5cWVlZysrKkt/v15YtW6LHOzo6VFVVpdzcXPXp00eVlZUKBoMJHxoA0PPFFaCBAwdq8eLFamho0O7duzVhwgRNnjxZb731liRp/vz52rx5s9avX6+6ujodOnRIU6dOPSuDAwB6thTnnDuTC+Tk5OjBBx/Uddddp/79+2vt2rW67rrrJElvv/22hgwZovr6eo0dO/YrXS8cDsvr9epqTVZaSvqZjAYAMHDcdWq7NikUCikrK+u053X7NaATJ05o3bp1am9vl9/vV0NDgzo7O1VeXh49p6SkREVFRaqvrz/tdSKRiMLhcMwGAEh+cQdo79696tOnjzwej2699VZt2LBBQ4cOVSAQUEZGhrKzs2PO9/l8CgQCp71eTU2NvF5vdCssLIx7EQCAnifuAF122WXas2ePdu7cqdmzZ2vGjBnav39/tweorq5WKBSKbi0tLd2+FgCg50iL9wEZGRm65JJLJEmlpaXatWuXHnnkEU2bNk3Hjh1Ta2trzF1QMBhUXl7eaa/n8Xjk8XjinxwA0KOd8eeAurq6FIlEVFpaqvT0dNXW1kaPNTY2qrm5WX6//0y/DQAgycR1B1RdXa1JkyapqKhIbW1tWrt2rbZv365t27bJ6/Vq5syZWrBggXJycpSVlaW5c+fK7/d/5XfAAQDOH3EF6MiRI/rJT36iw4cPy+v1avjw4dq2bZu+//3vS5KWLFmi1NRUVVZWKhKJqKKiQitWrDgrgwMAerYz/hxQovE5IADo2c7654AAADgTBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMDEGQVo8eLFSklJ0bx586L7Ojo6VFVVpdzcXPXp00eVlZUKBoNnOicAIMl0O0C7du3SY489puHDh8fsnz9/vjZv3qz169errq5Ohw4d0tSpU894UABAculWgI4eParp06friSeeUL9+/aL7Q6GQnnzyST388MOaMGGCSktLtXr1av3jH//Qjh07EjY0AKDn61aAqqqqdM0116i8vDxmf0NDgzo7O2P2l5SUqKioSPX19ae8ViQSUTgcjtkAAMkvLd4HrFu3Tq+//rp27dp10rFAIKCMjAxlZ2fH7Pf5fAoEAqe8Xk1NjX7961/HOwYAoIeL6w6opaVFt99+u5555hn16tUrIQNUV1crFApFt5aWloRcFwBwbosrQA0NDTpy5IiuvPJKpaWlKS0tTXV1dVq2bJnS0tLk8/l07Ngxtba2xjwuGAwqLy/vlNf0eDzKysqK2QAAyS+uH8FNnDhRe/fujdl38803q6SkRHfeeacKCwuVnp6u2tpaVVZWSpIaGxvV3Nwsv9+fuKkBAD1eXAHq27evhg0bFrOvd+/eys3Nje6fOXOmFixYoJycHGVlZWnu3Lny+/0aO3Zs4qYGAPR4cb8J4cssWbJEqampqqysVCQSUUVFhVasWJHobwMA6OFSnHPOeojPCofD8nq9ulqTlZaSbj0OACBOx12ntmuTQqHQF76uz++CAwCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAm4grQvffeq5SUlJitpKQkeryjo0NVVVXKzc1Vnz59VFlZqWAwmPChAQA9X9x3QJdffrkOHz4c3V555ZXosfnz52vz5s1av3696urqdOjQIU2dOjWhAwMAkkNa3A9IS1NeXt5J+0OhkJ588kmtXbtWEyZMkCStXr1aQ4YM0Y4dOzR27NhTXi8SiSgSiUS/DofD8Y4EAOiB4r4DOnDggAoKCnTxxRdr+vTpam5uliQ1NDSos7NT5eXl0XNLSkpUVFSk+vr6016vpqZGXq83uhUWFnZjGQCAniauAJWVlWnNmjXaunWrVq5cqaamJl111VVqa2tTIBBQRkaGsrOzYx7j8/kUCAROe83q6mqFQqHo1tLS0q2FAAB6lrh+BDdp0qTov4cPH66ysjINGjRIzz77rDIzM7s1gMfjkcfj6dZjAQA91xm9DTs7O1uXXnqpDh48qLy8PB07dkytra0x5wSDwVO+ZgQAOL+dUYCOHj2qd999V/n5+SotLVV6erpqa2ujxxsbG9Xc3Cy/33/GgwIAkktcP4L75S9/qWuvvVaDBg3SoUOHtGjRIl1wwQW68cYb5fV6NXPmTC1YsEA5OTnKysrS3Llz5ff7T/sOOADA+SuuAP373//WjTfeqP/85z/q37+/xo8frx07dqh///6SpCVLlig1NVWVlZWKRCKqqKjQihUrzsrgAICeLcU556yH+KxwOCyv16urNVlpKenW4wAA4nTcdWq7NikUCikrK+u05/G74AAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEzEHaD3339fN910k3Jzc5WZmakrrrhCu3fvjh53zmnhwoXKz89XZmamysvLdeDAgYQODQDo+eIK0EcffaRx48YpPT1dW7Zs0f79+/XQQw+pX79+0XMeeOABLVu2TKtWrdLOnTvVu3dvVVRUqKOjI+HDAwB6rrR4Tv7d736nwsJCrV69OrqvuLg4+m/nnJYuXaq7775bkydPliQ9/fTT8vl82rhxo2644YYEjQ0A6OniugN67rnnNGrUKF1//fUaMGCARo4cqSeeeCJ6vKmpSYFAQOXl5dF9Xq9XZWVlqq+vP+U1I5GIwuFwzAYASH5xBei9997TypUrNXjwYG3btk2zZ8/WbbfdpqeeekqSFAgEJEk+ny/mcT6fL3rs82pqauT1eqNbYWFhd9YBAOhh4gpQV1eXrrzySt1///0aOXKkZs2apVtuuUWrVq3q9gDV1dUKhULRraWlpdvXAgD0HHEFKD8/X0OHDo3ZN2TIEDU3N0uS8vLyJEnBYDDmnGAwGD32eR6PR1lZWTEbACD5xRWgcePGqbGxMWbfO++8o0GDBkn65A0JeXl5qq2tjR4Ph8PauXOn/H5/AsYFACSLuN4FN3/+fH3729/W/fffrx//+Md67bXX9Pjjj+vxxx+XJKWkpGjevHm67777NHjwYBUXF+uee+5RQUGBpkyZcjbmBwD0UHEFaPTo0dqwYYOqq6v1m9/8RsXFxVq6dKmmT58ePeeOO+5Qe3u7Zs2apdbWVo0fP15bt25Vr169Ej48AKDnSnHOOeshPiscDsvr9epqTVZaSrr1OACAOB13ndquTQqFQl/4uj6/Cw4AYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMBHXb8P+Onz6u1GPq1M6p35NKgDgqziuTkn///n8dM65ALW1tUmSXtH/Gk8CADgTbW1t8nq9pz1+zv05hq6uLh06dEh9+/ZVW1ubCgsL1dLSktR/qjscDrPOJHE+rFFinckm0et0zqmtrU0FBQVKTT39Kz3n3B1QamqqBg4cKOmTv7AqSVlZWUn9n/8p1pk8zoc1Sqwz2SRynV905/Mp3oQAADBBgAAAJs7pAHk8Hi1atEgej8d6lLOKdSaP82GNEutMNlbrPOfehAAAOD+c03dAAIDkRYAAACYIEADABAECAJggQAAAE+d0gJYvX65vfvOb6tWrl8rKyvTaa69Zj3RGXn75ZV177bUqKChQSkqKNm7cGHPcOaeFCxcqPz9fmZmZKi8v14EDB2yG7aaamhqNHj1affv21YABAzRlyhQ1NjbGnNPR0aGqqirl5uaqT58+qqysVDAYNJq4e1auXKnhw4dHPznu9/u1ZcuW6PFkWOPnLV68WCkpKZo3b150XzKs895771VKSkrMVlJSEj2eDGv81Pvvv6+bbrpJubm5yszM1BVXXKHdu3dHj3/dz0HnbID+8pe/aMGCBVq0aJFef/11jRgxQhUVFTpy5Ij1aN3W3t6uESNGaPny5ac8/sADD2jZsmVatWqVdu7cqd69e6uiokIdHR1f86TdV1dXp6qqKu3YsUMvvPCCOjs79YMf/EDt7e3Rc+bPn6/Nmzdr/fr1qqur06FDhzR16lTDqeM3cOBALV68WA0NDdq9e7cmTJigyZMn66233pKUHGv8rF27dumxxx7T8OHDY/Ynyzovv/xyHT58OLq98sor0WPJssaPPvpI48aNU3p6urZs2aL9+/froYceUr9+/aLnfO3PQe4cNWbMGFdVVRX9+sSJE66goMDV1NQYTpU4ktyGDRuiX3d1dbm8vDz34IMPRve1trY6j8fj/vznPxtMmBhHjhxxklxdXZ1z7pM1paenu/Xr10fP+ec//+kkufr6eqsxE6Jfv37uD3/4Q9Ktsa2tzQ0ePNi98MIL7rvf/a67/fbbnXPJ83+5aNEiN2LEiFMeS5Y1OufcnXfe6caPH3/a4xbPQefkHdCxY8fU0NCg8vLy6L7U1FSVl5ervr7ecLKzp6mpSYFAIGbNXq9XZWVlPXrNoVBIkpSTkyNJamhoUGdnZ8w6S0pKVFRU1GPXeeLECa1bt07t7e3y+/1Jt8aqqipdc801MeuRkuv/8sCBAyooKNDFF1+s6dOnq7m5WVJyrfG5557TqFGjdP3112vAgAEaOXKknnjiiehxi+egczJAH374oU6cOCGfzxez3+fzKRAIGE11dn26rmRac1dXl+bNm6dx48Zp2LBhkj5ZZ0ZGhrKzs2PO7Ynr3Lt3r/r06SOPx6Nbb71VGzZs0NChQ5NqjevWrdPrr7+umpqak44lyzrLysq0Zs0abd26VStXrlRTU5OuuuoqtbW1Jc0aJem9997TypUrNXjwYG3btk2zZ8/WbbfdpqeeekqSzXPQOffnGJA8qqqqtG/fvpifpyeTyy67THv27FEoFNJf//pXzZgxQ3V1ddZjJUxLS4tuv/12vfDCC+rVq5f1OGfNpEmTov8ePny4ysrKNGjQID377LPKzMw0nCyxurq6NGrUKN1///2SpJEjR2rfvn1atWqVZsyYYTLTOXkHdNFFF+mCCy446Z0mwWBQeXl5RlOdXZ+uK1nWPGfOHD3//PN66aWXon/fSfpknceOHVNra2vM+T1xnRkZGbrkkktUWlqqmpoajRgxQo888kjSrLGhoUFHjhzRlVdeqbS0NKWlpamurk7Lli1TWlqafD5fUqzz87Kzs3XppZfq4MGDSfN/KUn5+fkaOnRozL4hQ4ZEf9xo8Rx0TgYoIyNDpaWlqq2tje7r6upSbW2t/H6/4WRnT3FxsfLy8mLWHA6HtXPnzh61Zuec5syZow0bNujFF19UcXFxzPHS0lKlp6fHrLOxsVHNzc09ap2n0tXVpUgkkjRrnDhxovbu3as9e/ZEt1GjRmn69OnRfyfDOj/v6NGjevfdd5Wfn580/5eSNG7cuJM+EvHOO+9o0KBBkoyeg87KWxsSYN26dc7j8bg1a9a4/fv3u1mzZrns7GwXCASsR+u2trY298Ybb7g33njDSXIPP/ywe+ONN9y//vUv55xzixcvdtnZ2W7Tpk3uzTffdJMnT3bFxcXu448/Np78q5s9e7bzer1u+/bt7vDhw9Htv//9b/ScW2+91RUVFbkXX3zR7d692/n9fuf3+w2njt9dd93l6urqXFNTk3vzzTfdXXfd5VJSUtzf//5351xyrPFUPvsuOOeSY52/+MUv3Pbt211TU5N79dVXXXl5ubvooovckSNHnHPJsUbnnHvttddcWlqa++1vf+sOHDjgnnnmGXfhhRe6P/3pT9Fzvu7noHM2QM459+ijj7qioiKXkZHhxowZ43bs2GE90hl56aWXnKSTthkzZjjnPnkb5D333ON8Pp/zeDxu4sSJrrGx0XboOJ1qfZLc6tWro+d8/PHH7uc//7nr16+fu/DCC92PfvQjd/jwYbuhu+FnP/uZGzRokMvIyHD9+/d3EydOjMbHueRY46l8PkDJsM5p06a5/Px8l5GR4b7xjW+4adOmuYMHD0aPJ8MaP7V582Y3bNgw5/F4XElJiXv88cdjjn/dz0H8PSAAgIlz8jUgAEDyI0AAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYOL/AeSnBbD9CTXXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(ot.emd(A.ravel(),B.ravel(),M))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baba7763",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}