{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# OT distance on the Circle\n\nShows how to compute the Wasserstein distance on the circle\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Author: Cl\u00e9ment Bonet \n#\n# License: MIT License\n\n# sphinx_gallery_thumbnail_number = 2\n\nimport numpy as np\nimport matplotlib.pylab as pl\nimport ot\n\nfrom scipy.special import iv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot data\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def pdf_von_Mises(theta, mu, kappa):\n pdf = np.exp(kappa * np.cos(theta - mu)) / (2.0 * np.pi * iv(0, kappa))\n return pdf\n\n\nt = np.linspace(0, 2 * np.pi, 1000, endpoint=False)\n\nmu1 = 1\nkappa1 = 20\n\nmu_targets = np.linspace(mu1, mu1 + 2 * np.pi, 10)\n\n\npdf1 = pdf_von_Mises(t, mu1, kappa1)\n\n\npl.figure(1)\nfor k, mu in enumerate(mu_targets):\n pdf_t = pdf_von_Mises(t, mu, kappa1)\n if k == 0:\n label = \"Source distributions\"\n else:\n label = None\n pl.plot(t / (2 * np.pi), pdf_t, c=\"b\", label=label)\n\npl.plot(t / (2 * np.pi), pdf1, c=\"r\", label=\"Target distribution\")\npl.legend()\n\nmu2 = 0\nkappa2 = kappa1\n\nx1 = np.random.vonmises(mu1, kappa1, size=(10,)) + np.pi\nx2 = np.random.vonmises(mu2, kappa2, size=(10,)) + np.pi\n\nangles = np.linspace(0, 2 * np.pi, 150)\n\npl.figure(2)\npl.plot(np.cos(angles), np.sin(angles), c=\"k\")\npl.xlim(-1.25, 1.25)\npl.ylim(-1.25, 1.25)\npl.scatter(np.cos(x1), np.sin(x1), c=\"b\")\npl.scatter(np.cos(x2), np.sin(x2), c=\"r\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare the Euclidean Wasserstein distance with the Wasserstein distance on the circle\nThis examples illustrates the periodicity of the Wasserstein distance on the circle.\nWe choose as target distribution a von Mises distribution with mean $\\mu_{\\mathrm{target}}$\nand $\\kappa=20$. Then, we compare the distances with samples obtained from a von Mises distribution\nwith parameters $\\mu_{\\mathrm{source}}$ and $\\kappa=20$.\nThe Wasserstein distance on the circle takes into account the periodicity\nand attains its maximum in $\\mu_{\\mathrm{target}}+1$ (the antipodal point) contrary to the\nEuclidean version.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mu_targets = np.linspace(0, 2 * np.pi, 200)\nxs = np.random.vonmises(mu1 - np.pi, kappa1, size=(500,)) + np.pi\n\nn_try = 5\n\nxts = np.zeros((n_try, 200, 500))\nfor i in range(n_try):\n for k, mu in enumerate(mu_targets):\n # np.random.vonmises deals with data on [-pi, pi[\n xt = np.random.vonmises(mu - np.pi, kappa2, size=(500,)) + np.pi\n xts[i, k] = xt\n\n# Put data on S^1=[0,1[\nxts2 = xts / (2 * np.pi)\nxs2 = np.concatenate([xs[None] for k in range(200)], axis=0) / (2 * np.pi)\n\nL_w2_circle = np.zeros((n_try, 200))\nL_w2 = np.zeros((n_try, 200))\nL_lcot = np.zeros((n_try, 200))\n\nfor i in range(n_try):\n w2_circle = ot.wasserstein_circle(xs2.T, xts2[i].T, p=2)\n w2 = ot.wasserstein_1d(xs2.T, xts2[i].T, p=2)\n w_lcot = ot.linear_circular_ot(xs2.T, xts2[i].T)\n\n L_w2_circle[i] = w2_circle\n L_w2[i] = w2\n L_lcot[i] = w_lcot\n\nm_w2_circle = np.mean(L_w2_circle, axis=0)\nstd_w2_circle = np.std(L_w2_circle, axis=0)\n\nm_w2_lcot = np.mean(L_lcot, axis=0)\nstd_w2_lcot = np.std(L_lcot, axis=0)\n\nm_w2 = np.mean(L_w2, axis=0)\nstd_w2 = np.std(L_w2, axis=0)\n\npl.figure(1)\npl.plot(mu_targets / (2 * np.pi), m_w2_circle, label=\"Wasserstein circle\")\npl.fill_between(\n mu_targets / (2 * np.pi),\n m_w2_circle - 2 * std_w2_circle,\n m_w2_circle + 2 * std_w2_circle,\n alpha=0.5,\n)\npl.plot(mu_targets / (2 * np.pi), m_w2, label=\"Euclidean Wasserstein\")\npl.fill_between(\n mu_targets / (2 * np.pi), m_w2 - 2 * std_w2, m_w2 + 2 * std_w2, alpha=0.5\n)\npl.plot(mu_targets / (2 * np.pi), m_w2_lcot, label=\"Linear COT\")\npl.fill_between(\n mu_targets / (2 * np.pi),\n m_w2_lcot - 2 * std_w2_lcot,\n m_w2_lcot + 2 * std_w2_lcot,\n alpha=0.5,\n)\npl.vlines(\n x=[mu1 / (2 * np.pi)],\n ymin=0,\n ymax=np.max(w2),\n linestyle=\"--\",\n color=\"k\",\n label=r\"$\\mu_{\\mathrm{target}}$\",\n)\npl.legend()\npl.xlabel(r\"$\\mu_{\\mathrm{source}}$\")\npl.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wasserstein distance between von Mises and uniform for different kappa\nWhen $\\kappa=0$, the von Mises distribution is the uniform distribution on $S^1$.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kappas = np.logspace(-5, 2, 100)\nn_try = 20\n\nxts = np.zeros((n_try, 100, 500))\nfor i in range(n_try):\n for k, kappa in enumerate(kappas):\n # np.random.vonmises deals with data on [-pi, pi[\n xt = np.random.vonmises(0, kappa, size=(500,)) + np.pi\n xts[i, k] = xt / (2 * np.pi)\n\nL_w2 = np.zeros((n_try, 100))\nL_lcot = np.zeros((n_try, 100))\nfor i in range(n_try):\n L_w2[i] = ot.semidiscrete_wasserstein2_unif_circle(xts[i].T)\n L_lcot[i] = ot.linear_circular_ot(xts[i].T)\n\nm_w2 = np.mean(L_w2, axis=0)\nstd_w2 = np.std(L_w2, axis=0)\n\nm_lcot = np.mean(L_lcot, axis=0)\nstd_lcot = np.std(L_lcot, axis=0)\n\npl.figure(1)\npl.plot(kappas, m_w2, label=\"Wasserstein\")\npl.fill_between(kappas, m_w2 - std_w2, m_w2 + std_w2, alpha=0.5)\npl.plot(kappas, m_lcot, label=\"LCOT\")\npl.fill_between(kappas, m_lcot - std_lcot, m_lcot + std_lcot, alpha=0.5)\npl.legend()\npl.title(r\"Evolution of $W_2^2(vM(0,\\kappa), Unif(S^1))$\")\npl.xlabel(r\"$\\kappa$\")\npl.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.18" } }, "nbformat": 4, "nbformat_minor": 0 }