add test for quadratic OT and POT.smooth_ot_dual

Author: zsteve

Project.toml

@@ -8,7 +8,9 @@ Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]

src/pot.jl

@@ -154,4 +154,10 @@ function sinkhorn_unbalanced2(
     )[1]
 end
 
+function smooth_ot_dual(
+    mu, nu, C, eps, reg_type = "l2", method = "L-BFGS-B", tol = 1e-9, max_iter = 500, verbose = false
+)
+    return pot.smooth.smooth_ot_dual(nu, mu, PyReverseDims(C), eps, reg_type = reg_type, method = method, stopThr = tol, numItermax = max_iter)'
+end
+
 end

test/runtests.jl

@@ -5,6 +5,7 @@ using Distances
 using PyCall
 using Tulip
 using MathOptInterface
+using SparseArrays
 
 using LinearAlgebra
 using Random
@@ -171,3 +172,23 @@ end
         @test norm(γ - γ_pot, Inf) < 1e-9
     end
 end
+
+@testset "quadratic optimal transport" begin
+    M = 250
+    N = 200
+    @testset "example" begin
+        # create two uniform histograms
+        μ = fill(1 / M, M)
+        ν = fill(1 / N, N)
+
+        # create random cost matrix
+        C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)
+
+        # compute optimal transport map (Julia implementation + POT)
+        eps = 0.5
+        γ = quadreg(μ, ν, C, eps)
+        γ_pot = sparse(POT.smooth_ot_dual(μ, ν, C, eps))
+        # need to use a larger tolerance here because of a quirk with the POT solver 
+        @test norm(γ - γ_pot, Inf) < 0.5e-4 
+    end
+end