notebooks/algo7.1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# アルゴリズム 7.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! pip install POT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import ot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f9e271a9b80>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAneklEQVR4nO2df5QdZZnnv093upNuYNNKx8l2d0I4LpsZDD1kCchMsutMAhMYTIiwtOCi7jgQOYyijichCBNDBklI9gyIDkezwHFmUaEdkpAxsJEfipMsOEkINirDrCOrpJtIgnZUumP/evaPe2/n3ttV99aP96236t7v55ycTtetrnqruutbz/u8zw9RVRBCCMkuDa4HQAghJB4UckIIyTgUckIIyTgUckIIyTgUckIIyTjTXJy0vb1d582b5+LUhBCSWQ4ePHhMVWeVb3ci5PPmzcOBAwdcnJoQQjKLiPzUaztdK4QQknEo5IQQknEo5IQQknEo5IQQknEo5IQQknEo5IT09QJ3LwA2tOW+9vW6HhEhoXASfkhIaujrBf7xJmB0OPf98ddy3wNAd4+7cTlm56F+bN3zCgYGh9HR1oI1y+dj1cJO18MiPtAiJ/XN0xtPiniB0eHc9jpl56F+3LL9JfQPDkMB9A8O45btL2HnoX7XQyM+UMhJfXP8cLjtdcDWPa9geHS8ZNvw6Di27nnF0YhINSjkpL6Z2RVuex0wMDgcajtxD4Wc1DfL1gNNLaXbmlpy2+uUjraWUNuJeyjkpL7p7gFW3AvMnANAcl9X3FvXC51rls9HS1NjybaWpkasWT7f0YhINRi1Qkh3T10LdzmF6BRGrWQHCnlc+npzEQ7HD+f8qsvWUxRI5lm1sJPCnSEo5HFgDDIhJAXE9pGLyBwR+baIvCwiPxSRT5gYWCZgDDIhJAWYsMjHAHxaVV8QkdMAHBSRJ1X1RwaOnW4Yg0wISQGxhVxVXwfwev7/vxaRlwF0Aqh9IZ/ZlXOneG2PCFOjCSFhMeojF5F5ABYC+J7HZ6sBrAaAuXPnmjytO5atL/WRA7FikAup0YWsukJqNIDExDwLL5IsjJGQJBFVNXMgkVMBPAvgc6q6vdK+ixYt0prp2WkwamXx5mfQ75E919nWgn3rlsYdaVXKXyRALn540xXnRBJKG4JreowkffBF7Y+IHFTVReXbjVjkItIE4FEAX60m4jWHwRhk16nRlWpshH2QbM0uTI6RpI80zEqziImoFQHwAICXVfVv4g+pfnGdGm3yRWKr8JLfWPoHh7F48zOs0JdxWLArGiZS9BcD+CCApSLyYv7fnxo4bt3hOjXa5IvE1uyi0lhYbjX7uJ6VZpXYQq6qe1VVVLVbVc/N/3vcxODqjVULO7HpinPQ2dYCQc43nqTv1+SLxNbswmuMxdB6yzauZ6VZhZmdKcNlarTJGhtrls/3XJSMO7soHqPXwjBgz3rjIpx9bP3d1DoUclKCqReJzcJLhTH6RfnYsN64CJcMLNgVDWPhh2GoqfBD4owkQxFdh4YSAlgOPyTEBbatt2JXip+5w0U4kgYo5CTT2FpT8LL2veAiHEkD7BBUb/T1AncvADa05b729boeUSrH5BXPXA4X4UhaoEVeT6Sxfnoax4TKLhMBuAhHUgWFvJ6oVD/dlWimcUzICTUXN0lWoGulnkhj/fQ0jgnus2wJCQOFvJ7wq5Meo356VHYe6sfizc/g8MTp3js4GFMxrrNsCQkDXSv1hOH66VEpjgjZ0tCD/9G0Dc0ydnKHxubEx+QFGxCTrECLvJ7o7gFW3AvMnANAcl9X3Ju4L7o8IkTLo7QdJKkRkmVokdcbBuunR6U4ImTttF5Ml7Iwv4lR54udxD22a9vUUu0cCjlJnOKIkA455r1TBhtY15IwuMZ2bZtaq51D1wpJnOKIkAFt997J8WJnWArC0J9P52dt9HjYbjBRaw0sKOQkcYojQraO9WAY00t3cLAAG5daEwbX2G4wUWsNLOhaIU44GRGyFOhbaKyBtStqTRhc45eQZaq2je3jJw2FnLgnBQuwcak1YXCN7QYTcY6fxrUQCjkhBggjDGkUgrRhu0Rx1OOndZGUjSUIMUQQgbbeDKOvN/NuqjTjusEIG0sQYpkgmaCVFkVjC3lKK0nWEmldC2HUCiEJYlUIKlWSJEbwW/NwvRZCISfhSGETCOBkEa4z1+3G4s3PpDZ+26oQpLSSZC2R1qqYFHISnMLU/fhrAPTk1N2xmGcpGceqEKSkumVWXqpRSGtVTC52kuDcvSAv4mXMnAN86gfJjyeP6wWosFiLWin3kQO55KoEC6NZX8ytc7jYSeKT0qm7Eb9zgtEe1srjFsYb5DosXa/VxVziC4WcBGdml49FfnLq7iJGOnYyTi1FewRJrrJ4vWmN6qh16CMnwVm2PjdVL6aoLoorX3Vsv3O9RXtYvN60RnXUOhRyEoidh/qx+PF2fOKtP8MRzIJ6NKZwVTgq9gKUr8vIY/ZRC5h0kZVFMd1z9v9NZVRHrUPXCqlK8QJWP5bgsRNLcgtYf3QOVnWfFEuX0+pYfmc/lxEkJ1RZc69UI4CLLBAeLprzX/os/v782/HJH53FEgQJYsQiF5EHReQNEXEXukCsEdTSzuy0etl6AOLxgQJP3Jz0aOxTxUUWGB8Xzfn/9gXsW7cUr26+DPvWLaWIJ4Ap18pXAFxi6FgkZQS1tNOaLFGV7h6gvG9ogeFfOI+TN46p3q0pjWKqR4y4VlT1uyIyz8SxSDKEiS4JGhViu2KdVWbO8feJR+wfmuoqhyZKB5ty0ZDY0Edeh4QtxRmmRGtoX3VaqvUtWw9sv977swgWZlrLnRpl2XrvBKSMdXdKCpsv9sSiVkRktYgcEJEDR48eTeq0xIOw0SXW0pLTlPLf3QO0vN37swgWZlIRPE7T4U25aBLC5b2yHZprLEU/71r5pqouqLYvU/Tdcua63Z4eYQFw9/vPTc4d4DLl32smABhLca90j1/dfFnkYRfDdPjguL5XpspI+KXoM468DvGLImlrbUo2ocfVYpnfTAAwZmEmEcFTcw2fLVbWdH2vbIfmmgo//DqA5wDMF5HDIvLnJo5L7OAVXSIAfjk0muwfu6tqfZUyG7t7crOBDYO5rxHdBElE8NRUOrxlN5vre2X7xW5EyFX1GlX996rapKpdqvqAieMSOxT7vIGciFdysFn7YzcVzxyQgo90YtAnOsXgTCCJcqdBxSETZWUtl0lwneNg+8XOqJU6pSAon+79PsbL1klWNuzF2mm96JBjGNB23N98LQAzft0SwlTri0mxj3SguR1dcmzqToZnAtaqHOYJEk2UmegZy262MJFXNrAdmkshr1MKD7iXiG9uuh+tMgIA6JJjuE2/BPS9y040gol45gAU+0i3jPWUXCOATIbNBRGHzJSVtRyTnoYcB5svdgp5neL1gAPA2mm9pQIHYNr4ichJMWmh2D20a2IJMIr8rONNNLRlt9t8NXFw7RsOTAIx6bZnSC6hkNcpfg9yh5fLAch82nV5duquiSXYNbIkF/71KXddhGxnf8au1Z4UCbrZahEKeZ3i94C/jnZ0wt9/bFt4bB3ftY/UCy//9d4d9+FPvvUoWoePRBOzsvj4e87+OD60/4xUXbcvCbnZapGaiCPPxKp8yvBbRR84b61vJInt7DSbx09j09xy99bKhr3YKNvQOvw6IoXgeYTw5crK/tTsdVuM9ybRyHzzZdcZW1nG1/r1qX9iu8lx1poo+xF0VlGe/bm3+SZ0NXjNhgJmuiaRKZuCBs/1TPabL/uIS2ZW5U1guMCU7+KPzxTX9sJZqOOnpdhWGVXD/YrG/dyMdtw5clVu8RUG1ieSyJStlkxFnJANIa/QLHZg8BTPH0ndqnxcUtAg2PbCWeDjp+BeTBlPXpwvRDsuHr8Ku7Bk8uNJw6JxX8m4Z+Mo7mq6HxjNLb4OaMz49iTKyrIGeSrJho+8ghXgOmMrMVLQINh2dlrg46fgXkxS5peejaPY3HQ/VjbsLdltYHDYc9wtMoLPNH8DAuD+5msx1jij9PhhQvCSyJR1VVaBVCQbQl7BCshsV5qwpMASMrFgWGlhOvDxU3AvJvEQ51YZwdpppQuAHW0tvuObjWN4dfNl2HDb7Zh2+ReiF+1KoqxswmUVSDCy4VqpMGVMQ8ZWIqSkG0ucpIog6eKBjp+SewHAV5w75M3J/08aFt8JMO64IXi2Q/gY751KsiHkVbK+ajlja5Ia6MZibGHa615AgLP+xMxAg1Dwi/uUG3tD2iFAqWHRmP3fIQDGe6eQbAh5rVsBFSIwToaynYIPn/pRrG15JHqyiGOMRb109wA/ex448CBOCqkC3/8aMPdC+/fEKwSvmKYWzF5xJ17tLis0Vut/x8QZmY8jzzwV4nJ3ji8OFCOf6ia/RRiNE0+wu1D5/X1Sbswn7Xgwcw7FmVgj+3HktUqFCIytv723qisiM2VKETNNvnzW4tfxPuaCZ7lo//HvzsKjB/tL7u+M6UdyRdzLmIDgnT+/Cx2Pt2DNeH/i9z8rL3Tb1ON9oEXumg1t8POzHp5ox5axnsmEkQLFfR/TmA1Z6UGK9JB5ujJ82mFII6ATkdwWXlnChbMU12ifQAOmycSUnz880Y4lI/dOft/W0oQNK9+ViIgwwzlHpPuQ0uQyL2iRp5UK1mVXwzFsLkoYKVAcI5+2MqXVZgiRFqa9Zi1QeIq55h/gCElCXouxBREvrl/egAmoAlJklQ9pM7aMlZ5ncHjU3uyoTHxefOtKDI9eULJLzWY4VyD0gnrakssiko048lrGKy63iPKY5HJXRNoSoqw0ufV1l+jJmGlpnPpxyCQhv5efV412EeTPKTg80Y51o9dNmTkB0a69ahE4j+JYa0fvm5KEVOmaytm/68sY3NAF/exM6IaZ+O3nzghVDCsthetCGzZpSi6LAYXcNSVJHN50NLzpmyCTtoQoKzME32zCOScbJetUVweAUD5z35eiTw0U1QksnrEdS0bu9RTxAmGuPVAFyIBJSECwF/r+XV/G7x+8BW34NURy85zpo4MY33FjIDG3XRUzDKENmzQll8WAQp4GCp3bfcS8YWYXXt18GfatWzplehgl29Km9WRlhhAkm9BA6rjfS/F48+947j+gp3uuT5QT5toDzWgCJCEBwV/oc17YimaZ2i2qUUcDWaZWZmERCW3Y1EjJAQp5moiY/rxqYSf2rVvqK/bF2LaerMwQgqSeG0gd93spvm3FHVOOPYzpuGu0ug817LUHmtH4iMyJ1tmRyie8Q4/6fxjAMk3TOk1ow6ZGSg5wsTNNJJAwYrvsr7WSCdWyCQ3dO+/F2KnHXnd0hac7ZWXDXnym+Rt4hx7FGzILr/2nNTh/4SWBzx+oAqRPlm/rpRuxrzt8pNIbMguz4SPmASzTtLWTC7WgXiNJWgw/rDPKmxkUKA5pJNXxCvtc2bAXdzU/gBb89uTGkE0XAofPGQyZK/jIy90r49KExvfdV/W4DH1MDoYfJk1KY1PTZj1lFa/kppubektFHAjddCHwjMZgvZPzV34U+wGc9cJfY6b+GhBgpKkN01dsDXQOJ4XrUvp8uYIWuQ1S3A6rLqynvl4MPbEeM4aPYGDidNzffC3OvWx1/D6VZcKxc3xxiXjtPXEFxG++s2Ew+rlJKSl+vmzjZ5FTyG2QYB2QKNR0CnNfL8Ye+zimjZ+Y3DSkzVivq7HkfTdGu86gwpGy33vh99w/OIxGEYyrorMWft8pu89JQtdKkqQ8NrVmy/729QI7bsA0LfX1tsoItuKL+PljjwCNd4a32oL2qbRZajikK6F85jWeN9jSXIsnMCl/vlzA8EMbJBGb2tebs0w2tOW+hsjCq0kKVrNOjYcGcpmYs3E0t0/YexVUOGx16PHI5Kx2HV7RSQUqxXinJUOzIjUS+20SWuQ2sN0EIsH6EGHdMM7cNp71WDyI0vE9TEciG00XInSurxbD7fV5Zipp1kCTFdPQIreB7d6JCdWHCJs8FCfZKIolWPwzEyGm1Rp2Cu46aSSCK6FaFJLX52nK0KxIEr1JM4YRi1xELgHweQCNAO5X1c0mjptpbLbDSshHGDZ5KGqyURRLsPxnBiZOR1eDd02Ucn6OdswOtGce10kjYXuU9vXiSVmPGdOPYEBPn1IK2S/bNE0ZmlVhu7kSYlvkItII4G8BXArgbADXiMjZcY9LKpCQjzDsgx1VCKJYguU/s2WsB0PaXLLPiE7Db7W0XMCQNmPTyFUVx+NJoR7OhsHc1yRFJMyMIO92ax1+HQ2ik6WQVzXuA1A5ZT1tlTRJcEy4Vi4A8GNV/YmqjgB4GMDlBo5L/Ehoqh/2wY4qBH5Cv+hXT/ou6Jb/zK6JJVg3eh0OT7SjMN2+s+ljWDP6URyeaMeEniw3e+DfXVxxPKkjjCvBpzLiPe278P+q1OJJWyXNzOEwAMGEa6UTQPG87zCAd5fvJCKrAawGgLlz5xo4bR2T0FQ/bGu2qK3cvLJNVzbsxebmB4Dj+UzJsgVdr5/ZNbEEB1svnuyMdG7e/bJrpNStsCmLwhTUlRDD7eYkQ7NWcNygwoSQe3QvnJrepqrbAGwDcglBBs5b3yTgIwz7YEcVgijp7kFeGrGFyUYauO3U8rD+9DJqNsfANhEii0xiQsgPAygupN0FYMDAcUkKCPtgRxECL8HtOPGm9855yzKoSEcWJhsWVpRjhhV+hua5wXGSkgkh3w/gLBE5E0A/gKsBfMDAcUmd4BV7Lt+pbllatR5tWFhhjxlF+Lt7gJ89Dxz8Si45ShqB3/9A4DHftvMlfP17r2FcFY0iuObdc3DHqnOCXV89E3MmFJfYi52qOgbgYwD2AHgZQK+q/jDucUl94Bd7vv+dH09t7Hbk7MewVpuP8B/Z/hn/c/b1At//2skMVx3PfR9g4e22nS/hoed/NpnOP66Kh57/GW7b+VLVn81ERqhNHOcaGEkIUtXHVfU/quo7VfVzJo5JTlLLD4lf6OEnf3SW26QPH0tqqGV29A5LYcNGfQT+HXoMn3rkRczz+nuIkSz29e95WJQVthcw3XUqk3/vjpOUmNmZctLU2NYGxpNQTIWA+VhYW0bfHz37MazV5iPwv9RT8E/NN+En0z+AR4aux94d9538e4jhqx33qYTqt72AyYzQTP+9O8w1oJCnnMykTUfEL8b8w6f+c+hCUVGKS/niY2H93W8u8Nw90IsnrNXmIfy/1UacJifQ1XAMDQJ0NRzDRtmGF3dvy+0QxOr3edk1ilcAmv/2AiZfxrX+924LCnnKyVTadAT8klDWNj0SzEVQLEo7bjBbg6a7JyemM7tyFu3TG3MvGA8CZz+Gsdrywn8EsyYTmt5CC5plrGS3VhnBdSMP5b6pZvVXeNld8+458MJvewGTGaG1/vduCwp5ynGSNp1ghppf1/PW4SPeP1DsIigXJZ8StoFDwMqv+5t/OUX0btMv4b82/5+SH5uMX7dx37p78Pzlz+JdEw9jyci9aMNvPHfraHhzcv+KVn8FH/odq87BtRfOnbTAG0Vw7YVzq0atmMwIZZmAaLCMbcqJmi0ZGQcZap5hhAHCDwOXrg0SAuZ13QceRHlu27TxE9h4yqN4rnVZafx64z5r960kZn6oHV0ytTjYiZbZaC18UylZrIoP/Y5V54QONzSZEZr433uNQCFPOYmnTUeJnzaQrVjeluwyWeHdkb54YTCApT2M6WgJEgLm+VLwXuRrHT6CfRuWlm68225m3+TLrm/TlFZ2Y40z0HppQPeRpXhnUzH9LBMQDQp5Bkg0bTps1IMBC74QqXDx+LN4pLkXHXIMA9qO3rH/jIsaX0SHvAnxekH4iNKYNqABigE9HVvHevD5IC8gH9H2xEv0ksrs6+7JPbRFL85pYV6cGcj8ZJmA8FDISSlhLTY/C/6JmwOLy9Y9r+Di8Wexuel+tMoIAKBLjuGD8hT+1/hF2HbqX2Dfp5ZO/UEPURrSZqwbvW6y/nann2/Vq6GyJ4ISkfcTPZ/7dgTteP5Qv1lhilNnx3VtdWIFLnaSUpatx1jjjJJNY40z/C02P4tz+BeBF/sGBoexdlrvpIgXaBDgg41P5crZ5ilJFnm8HfvPuR2YOQcKQb+2l4h4Rd9qEP96Uwuw6CPBwgU9okWGtBl3jlyVvjhol7XViRVokdcRQfpp7hxfjL2j1+GTeBgd8iYG9HTcM3E1lowvxiqvg/pZ8EBg/3BHWws6hr27+zQIcEvzNwBs8uwk9KH9Z2DTFXuwamEn9h/qx8E9r0CC+FYrujwkvKWa3+/I9s/gHXqstDPPRPUuSUaotFZhouqi7cqNJDIU8iRJ6kHwOM/O8cWB2qlt3fMK+kf+EP+APyw55HN+QrRsPbD9eu9xBPQPr1k+H6/vbEcnvMX8d3B0cmyVWsmF8q36upDm5KzUKHT34A++doqnt916HHSltQogfkSN43rbpDJ0rSSFyazDCOd5cfe26hlzfb14ZOh6/GT6B7C3+SasbNg7+ZGvEHX3AC1v9/4sRA3sgfPWYsLnc4EAfb1mk0UsFTlyFgddKdrIRLPuhBp+k2hQyJMiqQfB5zyTmX9lTIpg/gVQnPq9uen+STGvKESX3hVbFM+f9zY0NDT6fKrA0xvNiqSlIkfO2qVVipoxEVHjuN42qQxdK0mR1IPgc7zJzL/y7QUR9On1uHZaL57U91QWIhOREE9vBCZ8MjMB4PhhrLnccLKIhS5LzuKgq0UbxY0dd1xvm1SGQp4UST0IPuc50TIbLeON/iJY4QWw6XLvrusllInizkP92Lr5meBiVu2FNrMrM8kiVuKgq62vVIsPjxs7noH483qGQp4UST0IPudpvXQjNo2f4y+CPi+AhiIBDYpXdInXwmoJlaJfiu5TXSaLBFloDDIrijNjYvx5qhGtUmvYBosWLdIDBw4kfl7nOIxaqXoerwSZppZIfuPFm5+Z0uEeyCXnFDrcBzo/kFtIvfSu+haMuxeYj7AhmUREDqrqovLttMiTxIJP1th5DFpckaJLUmTxBYm3TxQL6yupu0YSCwo5OYmhF01HW4unRd4ggjPX7fYXjqRedBXYv+vLOP/gFvwTjmGguR1bftWDW7bnMk5tC52vuBpeX4nk+jIJE4uMw/BDYhyvEDwg1zIs1e27+nqx4IW/QqeUhmBePP6s9Q41FVucGY55d9qFJ6l8ijqDQk78idgoobxZhFerMNPCYaRh79MbS8vmIheCuaHp7/HI0PVWG21UFFfDMe9+Li6vWZRxmFhkBbpWiDcxU7KLo0vOXLfbcx9TaetVXQVBp/I+Pue34Td4e0O+M4+l1PSq6woG3U5+ri9B7l5ada8wscgKtMiJNwYtJ9tp6xWt2TBTeR+f85QJhQULMsnU/jXL58OrnbIC9t0rQZpDk9BQyIk3Bi0nL595U6Pgv5z4Ng6vfycmNrRh6K7fjeyyqGjNhnkhefiifYNzDVuQSab2r1rY6Xtd1ot7WapxU+9QyIk3Bi2ncp/521qb8F7Zi7/SL+Vqu0DROvw6xh77eCQxr2jNhnkhefiiJWZBsKD4NaG25ebwa7hhvbiXpRo39Q4Tgog3BhOEylm8+Rk8MnQ9uho8ytZGSHIp95EDOWt20xXnYNV3lsdLprF0H1zHcVe8Z4wnTy1MCCLhsJigMzA4jI7p3rXHo7gsKtZgaYxZGsHkfcgvuurxwzhfT8d5oz3ox5Lk47jBJse1Bi1ykjimLfKqpCEBxcOy9+ov6lvCgBD4W+T0kZPEWbN8Pu7B1RjS5pLtFXuDxiENPSorlAkuYH2hkdQsdK2QxMlN32/Elt3TcN3IQ+hoeBMnWmaj9dJgPT4ziV+ZYDlZJ976QmOd43pdwiaxhFxErgKwAcDvAbhAVekvIYHIJQzdDuB2AECr2+HYx6deyi/1FAAJdRGqY5zXl7FMXNfKDwBcAeC7Bsbinogp6YRUZdl6oKFpyubT5AT++6n/zGgRyzitL5MAsYRcVV9W1dq4EyzmQ2zS3QNMP23K5mYZw4ZTHqWIW8Zo4+4Ukthip4isFpEDInLg6NGjSZ02OCzmQ2wz/Evv7RmtM2KkUFlCJFkCAUDis/uqPnIReQrAbI+PblXVx4KeSFW3AdgG5MIPA48wKVjMh9jGx08+1DIbF+e7KjWKYFwVnSlfjMuaz3nNcsONuysRs+BcFKpa5Kp6kaou8PgXWMQzAYv51D6u10A86oyMNc7A+reunKxGOJ7P60htzfY8WfM5J1oCwcHsnuGHBdglvLZxYCVNwSNL9I63rsQ/jFzguXtBGNNo4WbR55xY424Hs/u44YfvA/AFALMA7BaRF1V1uZGRJU2KekYSfyLHAleykir8jo3HHpfVFf87n1rtBYwIo4XMVr+a5oyFh/HWfEGIJeSqugPADkNjcU8KekYSf2L5ZSNYSUn4gf0EsfjzWHzzL4EDD2KyIK+hmUiiPues4WB2zxR94p6AvutYftkIayBJ+IH9+psCBoSxr7dUxAsY8NeuWtiJK8/rnGzj1yiCK89LyHWRdhyU6qWPnLglhO86ll82gpWUhB+4uAqh8aiVpzfCtzVGTH/tzkP9ePRg/+Ti7LgqHj3Yj0VnvJ1iDiQ+u6eQE7eE8F3H8stGWANJyg9sbRGukljH9NdWmq1QyJOHrhXilhC+69jt0EJWQUyy/ZoVfMVaYvtrsxi1UstQyIlbQviuk26HlvT5jOPVHxMCLPqIkaiVMNuJXdhYol5IQ3MFv3FZailHYO33ntVWcVkvZctWb/VMGpJh/MZV8JFLI6DjuRX+mGKT9YfVKJYW3bLYKi5rZQXCQIu8Hrh7QbwGxDaw2NQ4i5Yisc/ifD2bcrLUYo+t3uoZQynDxqrd9fUCO26wUo8iazVASHLU8gItXSv1gIGUYWPT0oIlruPen8eMb67lh5XEI4lwUlduPVrk9YBX9ELIlGEvS/fi8Wdx4WPvCVZNsJC9uf36qZZ4MTHjmxlNQfywHU5aMHb6B4ehSLaCJYW8HjCQMlxu0a5s2IvNTfdjNo6iakelku5LFTBQjyLzsd/EGrbDSV269ehaqRdiRi+UT0vXTutFq4yU7uRXTdAre7McaTQScpjFaAo/GH1jHpulbF269SjkJBDl1e465Jj3jl4+7mp+b8Nx44nVnbZILYfK1SouS/vStUICUTwtBYABbffe0cvHXcnvnUBluCwSepruuvsRcerWo5CTwKxa2Il965ais60FW8Z6MKTNJZ8PY7q3j9tvsfWK/xmo5kloakDUgk7Tdx7qx6O3vx8Tj16fX4Oosl5BrOGypANdKyQ0A4PD6McSYDTnK++QNzGgp2PrWA8+7yXKSXZfSmsWa0iCTNN3HurH3h33YYv8bzRI2Y4Buh8R87hy61HISWgKIrNrYgl2jSyZ3N5ZyReYVH3miC3dQpFA3ZogHXi27nkFj+DhqSJewGKPSJIu6FohoUl1iJ9vFutrZtwsJaGU9twYQabpA4PD/ovOgNUekSRd0CInoUl1iJ9fFitgxs2ShMWfp9o0vaOtBQND7ejyFPP4NcdJdqCQk0iE9QUmFhPt1dKtmLiia6hujQnWLJ+Pe3ZcjY26rSSmXwGIgZrjJDtQyIl1Eo2JLllY9bPMY4iugbo1psjduxuxZfc0XDfyEDoa3sSJltlovZSLnPUGy9gS6zgrH2qjfC8bYRCHsIwtcYaz1GUDxcKmYKBuDSGmoWuFWMdZ6rKt+PWkQikJCQiFnFgnSEy0NSi6pA6gkBPrpDpckZAagEJOEqEWKhISkla42EkIIRmHQk4IIRknlmtFRLYCWAFgBMC/AfgzVR00MC6SYdjZhpBkiWuRPwlggap2A/hXALfEHxLJMi4b0BJSr8QSclX9lqqO5b99HgDLrdU5LhvQElKvmPSRfwTAEwaPRzKIywa0pEaogQ5PSVPVRy4iTwGY7fHRrar6WH6fWwGMAfhqheOsBrAaAObOnRtpsCT9uGxAS2qAGunwlDRVLXJVvUhVF3j8K4j4hwG8F8B/0woVuFR1m6ouUtVFs2bNMncFJFWkuukEST+V6r0TX+JGrVwC4GYA71HVITNDIlmGWZwkFimq954l4mZ2fhHAdABPiggAPK+qN8QeFck0zOIkkUlRvfcsEUvIVfU/mBoIIYR4dniKW3q4DmCtFZIqmExU59gqPVzjUMhJaki0JRxJLyw9HBrWWiGpgclEhESDQk5SA5OJCIkGhZykBr+kISYTEVIZCjlJDUwmIiQaXOwkqYHJRIREg0JOUgWTiQgJD10rhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhBCScSjkhKSFvl7g7gXAhrbc175e1yMiGYH1yAlJA329wD/eBIzm+5Mefy33PcCO8qQqtMgJSQNPbzwp4gVGh3PbCakChZyQNHD8cLjthBRBISckDczsCredkCIo5ISkgWXrgaaW0m1NLbnthFSBQk5IGujuAVbcC8ycA0ByX1fcy4VOEghGrRCSFrp7KNwkErTICSEk48QSchH5axHpE5EXReRbItJhamCEEEKCEdci36qq3ap6LoBvAuDKDCGEJEwsIVfVXxV9ewoAjTccQgghYYm92CkinwPwIQDHAfxxhf1WA1gNAHPnzo17WkIIIXlEtbIRLSJPAZjt8dGtqvpY0X63AJihqp+telKRowB+GnKsWaEdwDHXg0iQerteoP6umdebHs5Q1VnlG6sKeVBE5AwAu1V1gZEDZhQROaCqi1yPIynq7XqB+rtmXm/6iRu1clbRtysB/Eu84RBCCAlLXB/5ZhGZD2ACOVfJDfGHRAghJAyxhFxVrzQ1kBpim+sBJEy9XS9Qf9fM6005xnzkhBBC3MAUfUIIyTgUckIIyTgUcguIyFYR+Zd8HZodItLmekw2EZGrROSHIjIhIpkK2wqDiFwiIq+IyI9FZJ3r8dhGRB4UkTdE5Aeux5IEIjJHRL4tIi/n/54/4XpMQaGQ2+FJAAtUtRvAvwK4xfF4bPMDAFcA+K7rgdhCRBoB/C2ASwGcDeAaETnb7ais8xUAl7geRIKMAfi0qv4egAsB/EVWfscUcguo6rdUdSz/7fMAarpfl6q+rKqvuB6HZS4A8GNV/YmqjgB4GMDljsdkFVX9LoBfuB5HUqjq66r6Qv7/vwbwMoBOt6MKBoXcPh8B8ITrQZDYdAJ4rej7w8jIQ07CIyLzACwE8D3HQwkEOwRFJEgNGhG5Fbnp2leTHJsNgtbcqWHEYxtjd2sQETkVwKMAPllW4TW1UMgjoqoXVfpcRD4M4L0AlmkNBOtXu9464DCAOUXfdwEYcDQWYgkRaUJOxL+qqttdjycodK1YQEQuAXAzgJWqOuR6PMQI+wGcJSJnikgzgKsB7HI8JmIQEREADwB4WVX/xvV4wkAht8MXAZwG4Ml8G7wvuR6QTUTkfSJyGMAfANgtIntcj8k0+cXrjwHYg9wiWK+q/tDtqOwiIl8H8ByA+SJyWET+3PWYLLMYwAcBLM0/ty+KyJ+6HlQQmKJPCCEZhxY5IYRkHAo5IYRkHAo5IYRkHAo5IYRkHAo5IYRkHAo5IYRkHAo5IYRknP8PD5BL8G5gf2oAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# データ生成\n",
    "np.random.seed(0)\n",
    "n = 100\n",
    "x = np.random.randn(n, 2)\n",
    "y = np.random.randn(n, 2)\n",
    "a = np.ones(n) / n\n",
    "b = np.ones(n) / n\n",
    "plt.scatter(x[:, 0], x[:, 1])\n",
    "plt.scatter(y[:, 0], y[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.22545794341828113\n",
      "1 0.22701966725048606\n",
      "10 0.8318351456244244\n"
     ]
    }
   ],
   "source": [
    "# 通常の最適輸送（外れ値に鋭敏）\n",
    "for p in [0, 1, 10]:\n",
    "    x[0, 0] = p # 外れ値\n",
    "    C = ((x.reshape(-1, 1, 2) - y.reshape(1, -1, 2)) ** 2).sum(2)\n",
    "    P = ot.emd(a, b, C)\n",
    "    W = (C * P).sum()\n",
    "    print(p, W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.26815043097940594\n",
      "1 0.2693680277800869\n",
      "10 0.27817950004349623\n"
     ]
    }
   ],
   "source": [
    "# 不均衡最適輸送（外れ値に頑健）\n",
    "for p in [0, 1, 10]:\n",
    "    x[0, 0] = p # 外れ値\n",
    "    C = ((x.reshape(-1, 1, 2) - y.reshape(1, -1, 2)) ** 2).sum(2)\n",
    "\n",
    "    # アルゴリズム 7.1\n",
    "    eps = 0.1\n",
    "    lam = 1\n",
    "    K = np.exp(- C / eps) # ギブスカーネルの計算\n",
    "\n",
    "    v = np.ones(n) # すべて 1 で初期化\n",
    "    for i in range(100):\n",
    "        u = np.maximum(np.exp(-lam/eps) * (K @ v), np.minimum(a, np.exp(lam/eps) * (K @ v))) / (K @ v)\n",
    "        v = np.maximum(np.exp(-lam/eps) * (K.T @ u), np.minimum(a, np.exp(lam/eps) * (K.T @ u))) / (K.T @ u)\n",
    "\n",
    "    P = u.reshape(n, 1) * K * v.reshape(1, n) # 主解\n",
    "    W = (C * P).sum() + lam * np.abs(P.sum(1) - a).sum() + lam * np.abs(P.sum(0) - b).sum()\n",
    "    print(p, W)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.3 ('base')",
   "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.3"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "44a2b068c0f563b0ccee5ab3b9cf7ca27efea6d57d435efed0db203d351dddef"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}