@@ -2,6 +2,7 @@ using ExactOptimalTransport
22
33using Distances
44using FillArrays
5+ using GLPK
56using PythonOT: PythonOT
67using Tulip
78using MathOptInterface
@@ -32,30 +33,34 @@ Random.seed!(100)
3233 pot_P = POT. emd (μ, ν, C)
3334 pot_cost = POT. emd2 (μ, ν, C)
3435
35- # compute optimal transport map and cost with Tulip
36- lp = Tulip. Optimizer ()
37- P = emd (μ, ν, C, lp)
38- @test size (C) == size (P)
39- @test MOI. get (lp, MOI. TerminationStatus ()) == MOI. OPTIMAL
40- @test maximum (abs, P .- pot_P) < 1e-2
41-
42- lp = Tulip. Optimizer ()
43- cost = emd2 (μ, ν, C, lp)
44- @test dot (C, P) ≈ cost atol = 1e-5
45- @test MOI. get (lp, MOI. TerminationStatus ()) == MOI. OPTIMAL
46- @test cost ≈ pot_cost atol = 1e-5
36+ # compute optimal transport map and cost with Tulip and GLPK
37+ for T in (Tulip. Optimizer, GLPK. Optimizer)
38+ lp = T ()
39+ P = emd (μ, ν, C, lp)
40+ @test size (C) == size (P)
41+ @test MOI. get (lp, MOI. TerminationStatus ()) == MOI. OPTIMAL
42+ @test maximum (abs, P .- pot_P) < 1e-2
43+
44+ lp = T ()
45+ cost = emd2 (μ, ν, C, lp)
46+ @test dot (C, P) ≈ cost atol = 1e-5
47+ @test MOI. get (lp, MOI. TerminationStatus ()) == MOI. OPTIMAL
48+ @test cost ≈ pot_cost atol = 1e-5
49+ end
4750 end
4851
4952 @testset " pre-computed plan" begin
5053 # create random cost matrix
5154 C = pairwise (SqEuclidean (), rand (1 , M), rand (1 , N); dims= 2 )
5255
53- # compute optimal transport map
54- P = emd (μ, ν, C, Tulip. Optimizer ())
56+ # compute optimal transport map with Tulip and GLPK
57+ for T in (Tulip. Optimizer, GLPK. Optimizer)
58+ P = emd (μ, ν, C, T ())
5559
56- # do not use μ and ν to ensure that provided map is used
57- cost = emd2 (similar (μ), similar (ν), C, Tulip. Optimizer (); plan= P)
58- @test cost ≈ emd2 (μ, ν, C, Tulip. Optimizer ())
60+ # do not use μ and ν to ensure that provided map is used
61+ cost = emd2 (similar (μ), similar (ν), C, T (); plan= P)
62+ @test cost ≈ emd2 (μ, ν, C, T ())
63+ end
5964 end
6065
6166 # https://github.com/JuliaOptimalTransport/OptimalTransport.jl/issues/71
@@ -106,8 +111,10 @@ Random.seed!(100)
106111 xs = rand (μ, m)
107112 μdiscrete = fill (1 / m, m)
108113 C = pairwise (Euclidean (), xs' , (1 : length (νprobs))' ; dims= 2 )
109- c2 = emd2 (μdiscrete, νprobs, C, Tulip. Optimizer ())
110- @test c2 ≈ c rtol = 1e-1
114+ for optimizer in (Tulip. Optimizer (), GLPK. Optimizer ())
115+ c2 = emd2 (μdiscrete, νprobs, C, optimizer)
116+ @test c2 ≈ c rtol = 1e-1
117+ end
111118 end
112119
113120 @testset " discrete case" begin
@@ -164,8 +171,10 @@ Random.seed!(100)
164171 # DiscreteNonParametric sorts the support automatically, here we have to sort
165172 # manually
166173 C = pairwise (Euclidean (), μsupport' , νsupport' ; dims= 2 )
167- c2 = emd2 (μprobs, νprobs, C, Tulip. Optimizer ())
168- @test c2 ≈ c rtol = 1e-5
174+ for optimizer in (Tulip. Optimizer (), GLPK. Optimizer ())
175+ c2 = emd2 (μprobs, νprobs, C, optimizer)
176+ @test c2 ≈ c rtol = 1e-5
177+ end
169178
170179 # compare with POT
171180 # disabled currently since https://github.com/PythonOT/POT/issues/169 causes bounds
@@ -222,8 +231,10 @@ Random.seed!(100)
222231 μprobs = normalize! (pdf (μ, μsupp' ), 1 )
223232 νprobs = normalize! (pdf (ν, νsupp' ), 1 )
224233 C = pairwise (SqEuclidean (), μsupp' , νsupp' ; dims= 2 )
225- @test emd2 (μprobs, νprobs, C, Tulip. Optimizer ()) ≈ ot_cost (SqEuclidean (), μ, ν) rtol =
226- 1e-3
234+ for optimizer in (Tulip. Optimizer (), GLPK. Optimizer ())
235+ @test emd2 (μprobs, νprobs, C, optimizer) ≈ ot_cost (SqEuclidean (), μ, ν) rtol =
236+ 1e-3
237+ end
227238
228239 # Use hcubature integration to perform ``\\int c(x,T(x)) d\\mu``
229240 T = ot_plan (SqEuclidean (), μ, ν)
0 commit comments