|
| 1 | +""" |
| 2 | +========================= |
| 3 | +Partial Wasserstein in 1D |
| 4 | +========================= |
| 5 | +
|
| 6 | +This script demonstrates how to compute and visualize the Partial Wasserstein distance between two 1D discrete distributions using `ot.partial.partial_wasserstein_1d`. |
| 7 | +
|
| 8 | +We illustrate the intermediate transport plans for all `k = 1...n`, where `n = min(len(x_a), len(x_b))`. |
| 9 | +""" |
| 10 | + |
| 11 | +# sphinx_gallery_thumbnail_number = 5 |
| 12 | + |
| 13 | +import numpy as np |
| 14 | +import matplotlib.pyplot as plt |
| 15 | +from ot.partial import partial_wasserstein_1d |
| 16 | + |
| 17 | + |
| 18 | +def plot_partial_transport( |
| 19 | + ax, x_a, x_b, indices_a=None, indices_b=None, marginal_costs=None |
| 20 | +): |
| 21 | + y_a = np.ones_like(x_a) |
| 22 | + y_b = -np.ones_like(x_b) |
| 23 | + min_min = min(x_a.min(), x_b.min()) |
| 24 | + max_max = max(x_a.max(), x_b.max()) |
| 25 | + |
| 26 | + ax.plot([min_min - 1, max_max + 1], [1, 1], "k-", lw=0.5, alpha=0.5) |
| 27 | + ax.plot([min_min - 1, max_max + 1], [-1, -1], "k-", lw=0.5, alpha=0.5) |
| 28 | + |
| 29 | + # Plot transport lines |
| 30 | + if indices_a is not None and indices_b is not None: |
| 31 | + subset_a = np.sort(x_a[indices_a]) |
| 32 | + subset_b = np.sort(x_b[indices_b]) |
| 33 | + |
| 34 | + for x_a_i, x_b_j in zip(subset_a, subset_b): |
| 35 | + ax.plot([x_a_i, x_b_j], [1, -1], "k--", alpha=0.7) |
| 36 | + |
| 37 | + # Plot all points |
| 38 | + ax.plot(x_a, y_a, "o", color="C0", label="x_a", markersize=8) |
| 39 | + ax.plot(x_b, y_b, "o", color="C1", label="x_b", markersize=8) |
| 40 | + |
| 41 | + if marginal_costs is not None: |
| 42 | + k = len(marginal_costs) |
| 43 | + ax.set_title( |
| 44 | + f"Partial Transport - k = {k}, Cumulative Cost = {sum(marginal_costs):.2f}", |
| 45 | + fontsize=16, |
| 46 | + ) |
| 47 | + else: |
| 48 | + ax.set_title("Original 1D Discrete Distributions", fontsize=16) |
| 49 | + ax.legend(loc="upper right", fontsize=14) |
| 50 | + ax.set_yticks([]) |
| 51 | + ax.set_xticks([]) |
| 52 | + ax.set_ylim(-2, 2) |
| 53 | + ax.set_xlim(min(x_a.min(), x_b.min()) - 1, max(x_a.max(), x_b.max()) + 1) |
| 54 | + ax.axis("off") |
| 55 | + |
| 56 | + |
| 57 | +# Simulate two 1D discrete distributions |
| 58 | +np.random.seed(0) |
| 59 | +n = 6 |
| 60 | +x_a = np.sort(np.random.uniform(0, 10, size=n)) |
| 61 | +x_b = np.sort(np.random.uniform(0, 10, size=n)) |
| 62 | + |
| 63 | +# Plot original distributions |
| 64 | +plt.figure(figsize=(6, 2)) |
| 65 | +plot_partial_transport(plt.gca(), x_a, x_b) |
| 66 | +plt.show() |
| 67 | + |
| 68 | +# %% |
| 69 | +indices_a, indices_b, marginal_costs = partial_wasserstein_1d(x_a, x_b) |
| 70 | + |
| 71 | +# Compute cumulative cost |
| 72 | +cumulative_costs = np.cumsum(marginal_costs) |
| 73 | + |
| 74 | +# Visualize all partial transport plans |
| 75 | +for k in range(n): |
| 76 | + plt.figure(figsize=(6, 2)) |
| 77 | + plot_partial_transport( |
| 78 | + plt.gca(), |
| 79 | + x_a, |
| 80 | + x_b, |
| 81 | + indices_a[: k + 1], |
| 82 | + indices_b[: k + 1], |
| 83 | + marginal_costs[: k + 1], |
| 84 | + ) |
| 85 | + plt.show() |
0 commit comments