From be8b9fa23984a226c3dc91bfd199d22eba36f268 Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Mon, 23 Dec 2024 13:05:09 -0600
Subject: [PATCH 01/12] working on benchmark to PATH

---
 benchmark/Project.toml | 7 +++++++
 1 file changed, 7 insertions(+)
 create mode 100644 benchmark/Project.toml

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
new file mode 100644
index 0000000..4cc4dd2
--- /dev/null
+++ b/benchmark/Project.toml
@@ -0,0 +1,7 @@
+[deps]
+MixedComplementarityProblems = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"
+PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b"
+ParametricMCPs = "9b992ff8-05bb-4ea1-b9d2-5ef72d82f7ad"
+
+[sources]
+MixedComplementarityProblems = {path = ".."}

From 9b99741b1422cdcbbffa22c374c5e0549c0f0f05 Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 09:07:46 -0600
Subject: [PATCH 02/12] first cut of benchmark loop

---
 benchmark/Project.toml        |   3 +
 benchmark/SolverBenchmarks.jl |  13 +++++
 benchmark/path.jl             | 101 ++++++++++++++++++++++++++++++++++
 test/runtests.jl              |   2 +-
 4 files changed, 118 insertions(+), 1 deletion(-)
 create mode 100644 benchmark/SolverBenchmarks.jl
 create mode 100644 benchmark/path.jl

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 4cc4dd2..a62ff56 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -2,6 +2,9 @@
 MixedComplementarityProblems = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"
 PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b"
 ParametricMCPs = "9b992ff8-05bb-4ea1-b9d2-5ef72d82f7ad"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [sources]
 MixedComplementarityProblems = {path = ".."}
diff --git a/benchmark/SolverBenchmarks.jl b/benchmark/SolverBenchmarks.jl
new file mode 100644
index 0000000..36f8d94
--- /dev/null
+++ b/benchmark/SolverBenchmarks.jl
@@ -0,0 +1,13 @@
+"Module for benchmarking different solvers against one another."
+module SolverBenchmarks
+
+using MixedComplementarityProblems: MixedComplementarityProblems
+using ParametricMCPs: ParametricMCPs
+using Random: Random
+using Statistics: Statistics
+using PATHSolver: PATHSolver
+using ProgressMeter: @showprogress
+
+include("path.jl")
+
+end # module SolverBenchmarks
diff --git a/benchmark/path.jl b/benchmark/path.jl
new file mode 100644
index 0000000..218f933
--- /dev/null
+++ b/benchmark/path.jl
@@ -0,0 +1,101 @@
+""" Generate a random large (convex) quadratic problem of the form
+                               min_x 0.5 xᵀ M x - θᵀ x
+                               s.t.  Ax - b ≥ 0.
+using Base: parameter_upper_bound
+
+NOTE: the problem may not be feasible!
+"""
+function generate_test_problem(
+    rng = Random.MersenneTwister(1);
+    num_primals = 1000,
+    num_inequalities = 1000,
+)
+    M = let
+        P = randn(rng, num_primals, num_primals)
+        P' * P
+    end
+
+    A = randn(rng, num_inequalities, num_primals)
+    b = randn(rng, num_inequalities)
+
+    G(x, y; θ) = M * x - θ - A' * y
+    H(x, y; θ) = A * x - b
+    K(z, θ) = begin
+        x = z[1:size(M, 1)]
+        y = z[(size(M, 1) + 1):end]
+
+        [G(x, y; θ); H(x, y; θ)]
+    end
+
+    (; G, H, K)
+end
+
+"Benchmark interior point solver against PATH on a bunch of random QPs."
+function benchmark(;
+    num_problems = 100,
+    num_samples_per_problem = 10,
+    num_primals = 1000,
+    num_inequalities = 1000,
+)
+    rng = Random.MersenneTwister(1)
+
+    # Generate random problems and parameters.
+    @showprogress desc = "Generating test problems..." problems = map(1:num_problems) do _
+        problem = generate_test_problem(rng; num_primals, num_inequalities)
+    end
+
+    θs = map(1:num_samples_per_problem) do _
+        randn(rng, num_primals)
+    end
+
+    # Generate corresponding MCPs.
+    @showprogress desc = "Generating IP MCPs... " ip_mcps = map(problems) do p
+        MixedComplementarityProblems.PrimalDualMCP(
+            p.G,
+            p.H;
+            unconstrained_dimension = num_primals,
+            constrained_dimension = num_inequalities,
+            parameter_dimension = num_primals,
+        )
+    end
+
+    @showprogress desc = "Generating PATH MCPs..." path_mcps = map(problems) do p
+        lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
+        upper_bounds = fill(Inf, num_primals + num_inequalities)
+        ParametricMCPs.ParametricMCP(p.K, lower_bounds, upper_bounds, num_primals)
+    end
+
+    @showprogress desc = "Solving IP MCPs..." ip_times = map(ip_mcps) do mcp
+        # Make sure everything is compiled in a dry run.
+        MixedComplementarityProblems.solve(
+            MixedComplementarityProblems.InteriorPoint(),
+            mcp,
+            zeros(num_primals),
+        )
+
+        # Solve and time.
+        map(θs) do θ
+            elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
+                MixedComplementarityProblems.InteriorPoint(),
+                mcp,
+                θ,
+            )
+
+            (; elapsed_time, sol.status)
+        end
+    end
+
+    @showprogress desc = "Solving PATH MCPs..." path_times = map(path_mcps) do mcp
+        # Make sure everything is compiled in a dry run.
+        ParametricMCPs.solve(mcp, zeros(num_primals))
+
+        # Solve and time.
+        map(θs) do θ
+            elapsed_time = @elapsed sol = ParametricMCPs.solve(mcp, θ)
+
+            (; elapsed_time, sol.status)
+        end
+    end
+
+    (; ip_times, path_times)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 716815f..1c2fd84 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -18,7 +18,7 @@ using FiniteDiff: FiniteDiff
     b = [1; 1]
     θ = rand(2)
 
-    G(x, y; θ) = M * x - A' * y - θ
+    G(x, y; θ) = M * x - θ - A' * y
     H(x, y; θ) = A * x - b
     K(z; θ) = begin
         x = z[1:size(M, 1)]

From b93078e6eb656b32fe54d105dceda21d30e62a11 Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 09:14:10 -0600
Subject: [PATCH 03/12] fixing syntax for progress meter calls

---
 benchmark/path.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/benchmark/path.jl b/benchmark/path.jl
index 218f933..443937d 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -40,7 +40,7 @@ function benchmark(;
     rng = Random.MersenneTwister(1)
 
     # Generate random problems and parameters.
-    @showprogress desc = "Generating test problems..." problems = map(1:num_problems) do _
+    problems = @showprogress desc = "Generating test problems..." map(1:num_problems) do _
         problem = generate_test_problem(rng; num_primals, num_inequalities)
     end
 
@@ -49,7 +49,7 @@ function benchmark(;
     end
 
     # Generate corresponding MCPs.
-    @showprogress desc = "Generating IP MCPs... " ip_mcps = map(problems) do p
+    ip_mcps = @showprogress desc = "Generating IP MCPs... " map(problems) do p
         MixedComplementarityProblems.PrimalDualMCP(
             p.G,
             p.H;
@@ -59,13 +59,13 @@ function benchmark(;
         )
     end
 
-    @showprogress desc = "Generating PATH MCPs..." path_mcps = map(problems) do p
+    path_mcps = @showprogress desc = "Generating PATH MCPs..." map(problems) do p
         lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
         upper_bounds = fill(Inf, num_primals + num_inequalities)
         ParametricMCPs.ParametricMCP(p.K, lower_bounds, upper_bounds, num_primals)
     end
 
-    @showprogress desc = "Solving IP MCPs..." ip_times = map(ip_mcps) do mcp
+    ip_times = @showprogress desc = "Solving IP MCPs..." map(ip_mcps) do mcp
         # Make sure everything is compiled in a dry run.
         MixedComplementarityProblems.solve(
             MixedComplementarityProblems.InteriorPoint(),
@@ -85,7 +85,7 @@ function benchmark(;
         end
     end
 
-    @showprogress desc = "Solving PATH MCPs..." path_times = map(path_mcps) do mcp
+    path_times = @showprogress desc = "Solving PATH MCPs..." map(path_mcps) do mcp
         # Make sure everything is compiled in a dry run.
         ParametricMCPs.solve(mcp, zeros(num_primals))
 

From 9e0a29a5dbb118bf6de562b42e5abd3c6d1969f8 Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 09:24:12 -0600
Subject: [PATCH 04/12] adjusting defaults

---
 benchmark/path.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/benchmark/path.jl b/benchmark/path.jl
index 443937d..2822bb0 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -32,16 +32,16 @@ end
 
 "Benchmark interior point solver against PATH on a bunch of random QPs."
 function benchmark(;
-    num_problems = 100,
-    num_samples_per_problem = 10,
-    num_primals = 1000,
-    num_inequalities = 1000,
+    num_problems = 10,
+    num_samples_per_problem = 100,
+    num_primals = 10,
+    num_inequalities = 10,
 )
     rng = Random.MersenneTwister(1)
 
     # Generate random problems and parameters.
     problems = @showprogress desc = "Generating test problems..." map(1:num_problems) do _
-        problem = generate_test_problem(rng; num_primals, num_inequalities)
+        generate_test_problem(rng; num_primals, num_inequalities)
     end
 
     θs = map(1:num_samples_per_problem) do _

From 8c7d91fac1183efeb385d4ab0ead64b0e8f85a4b Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 10:46:06 -0600
Subject: [PATCH 05/12] can compute summary statistics now. should be timing
 how long it takes to build mcps though. currently we are 10x slower than path

---
 benchmark/path.jl | 44 +++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 39 insertions(+), 5 deletions(-)

diff --git a/benchmark/path.jl b/benchmark/path.jl
index 2822bb0..445ca58 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -65,7 +65,7 @@ function benchmark(;
         ParametricMCPs.ParametricMCP(p.K, lower_bounds, upper_bounds, num_primals)
     end
 
-    ip_times = @showprogress desc = "Solving IP MCPs..." map(ip_mcps) do mcp
+    ip = @showprogress desc = "Solving IP MCPs..." map(ip_mcps) do mcp
         # Make sure everything is compiled in a dry run.
         MixedComplementarityProblems.solve(
             MixedComplementarityProblems.InteriorPoint(),
@@ -81,11 +81,11 @@ function benchmark(;
                 θ,
             )
 
-            (; elapsed_time, sol.status)
+            (; elapsed_time, success = sol.status == :solved)
         end
     end
 
-    path_times = @showprogress desc = "Solving PATH MCPs..." map(path_mcps) do mcp
+    path = @showprogress desc = "Solving PATH MCPs..." map(path_mcps) do mcp
         # Make sure everything is compiled in a dry run.
         ParametricMCPs.solve(mcp, zeros(num_primals))
 
@@ -93,9 +93,43 @@ function benchmark(;
         map(θs) do θ
             elapsed_time = @elapsed sol = ParametricMCPs.solve(mcp, θ)
 
-            (; elapsed_time, sol.status)
+            (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
         end
     end
 
-    (; ip_times, path_times)
+    (; ip, path)
+end
+
+"Compute summary statistics from solver benchmark data."
+function summary_statistics(data)
+    accumulate_stats(solver_data) = begin
+        (; success_rate = fraction_solved(solver_data), runtime_stats(solver_data)...)
+    end
+
+    (; ip = accumulate_stats(data.ip), path = accumulate_stats(data.path))
+end
+
+"Estimate mean and standard deviation of runtimes for all problems."
+function runtime_stats(solver_data)
+    stats = map(solver_data) do problem_data
+        filtered_times = map(
+            datum -> datum.elapsed_time,
+            filter(datum -> datum.success, problem_data),
+        )
+        μ = Statistics.mean(filtered_times)
+        σ = Statistics.stdm(filtered_times, μ)
+
+        (; μ, σ)
+    end
+
+    μ = map(datum -> datum.μ, stats)
+    σ = map(datum -> datum.σ, stats)
+    (; μ, σ, mean_μ = Statistics.mean(μ), mean_σ = Statistics.mean(σ))
+end
+
+"Compute fraction of problems solved."
+function fraction_solved(solver_data)
+    Statistics.mean(solver_data) do problem_data
+        Statistics.mean(datum -> datum.success, problem_data)
+    end
 end

From f22ef5c37d58777dd087cca2bb470a9c223d918b Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 16:37:37 -0600
Subject: [PATCH 06/12] now with 100 primal dimensions, 100 inequality
 constraints, and 90% random sparsity IP is only 5x slower than path

---
 benchmark/Project.toml        |  1 +
 benchmark/SolverBenchmarks.jl |  1 +
 benchmark/path.jl             | 18 +++++++++++++-----
 3 files changed, 15 insertions(+), 5 deletions(-)

diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index a62ff56..b781345 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 MixedComplementarityProblems = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"
 PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b"
 ParametricMCPs = "9b992ff8-05bb-4ea1-b9d2-5ef72d82f7ad"
diff --git a/benchmark/SolverBenchmarks.jl b/benchmark/SolverBenchmarks.jl
index 36f8d94..a93b40a 100644
--- a/benchmark/SolverBenchmarks.jl
+++ b/benchmark/SolverBenchmarks.jl
@@ -5,6 +5,7 @@ using MixedComplementarityProblems: MixedComplementarityProblems
 using ParametricMCPs: ParametricMCPs
 using Random: Random
 using Statistics: Statistics
+using Distributions: Distributions
 using PATHSolver: PATHSolver
 using ProgressMeter: @showprogress
 
diff --git a/benchmark/path.jl b/benchmark/path.jl
index 445ca58..33d06c5 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -9,13 +9,20 @@ function generate_test_problem(
     rng = Random.MersenneTwister(1);
     num_primals = 1000,
     num_inequalities = 1000,
+    sparsity_rate = 0.9,
 )
+    bernoulli = Distributions.Bernoulli(1 - sparsity_rate)
+
     M = let
-        P = randn(rng, num_primals, num_primals)
+        P =
+            randn(rng, num_primals, num_primals) .*
+            rand(rng, bernoulli, num_primals, num_primals)
         P' * P
     end
 
-    A = randn(rng, num_inequalities, num_primals)
+    A =
+        randn(rng, num_inequalities, num_primals) .*
+        rand(rng, bernoulli, num_inequalities, num_primals)
     b = randn(rng, num_inequalities)
 
     G(x, y; θ) = M * x - θ - A' * y
@@ -34,14 +41,15 @@ end
 function benchmark(;
     num_problems = 10,
     num_samples_per_problem = 100,
-    num_primals = 10,
-    num_inequalities = 10,
+    num_primals = 100,
+    num_inequalities = 100,
+    sparsity_rate = 0.9,
 )
     rng = Random.MersenneTwister(1)
 
     # Generate random problems and parameters.
     problems = @showprogress desc = "Generating test problems..." map(1:num_problems) do _
-        generate_test_problem(rng; num_primals, num_inequalities)
+        generate_test_problem(rng; num_primals, num_inequalities, sparsity_rate)
     end
 
     θs = map(1:num_samples_per_problem) do _

From 5a1b9ada793cd697b9117e7d8a0f3f2fe135897f Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 17:14:26 -0600
Subject: [PATCH 07/12] refactoring problem to avoid needing to construct
 separate mcps, by parameterizing M, A, b

---
 benchmark/path.jl | 174 +++++++++++++++++++++++++---------------------
 1 file changed, 95 insertions(+), 79 deletions(-)

diff --git a/benchmark/path.jl b/benchmark/path.jl
index 33d06c5..ba63434 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -1,16 +1,35 @@
 """ Generate a random large (convex) quadratic problem of the form
-                               min_x 0.5 xᵀ M x - θᵀ x
+                               min_x 0.5 xᵀ M x - ϕᵀ x
                                s.t.  Ax - b ≥ 0.
-using Base: parameter_upper_bound
 
 NOTE: the problem may not be feasible!
 """
-function generate_test_problem(
-    rng = Random.MersenneTwister(1);
-    num_primals = 1000,
-    num_inequalities = 1000,
-    sparsity_rate = 0.9,
-)
+function generate_test_problem(; num_primals, num_inequalities)
+    G(x, y; θ) =
+        let
+            (; M, A, ϕ) = unpack_parameters(θ; num_primals, num_inequalities)
+            M * x - ϕ - A' * y
+        end
+
+    H(x, y; θ) =
+        let
+            (; A, b) = unpack_parameters(θ; num_primals, num_inequalities)
+            A * x - b
+        end
+
+    K(z, θ) =
+        let
+            x = z[1:num_primals]
+            y = z[(num_primals + 1):end]
+
+            [G(x, y; θ); H(x, y; θ)]
+        end
+
+    (; G, H, K)
+end
+
+"Generate a random parameter vector Θ corresponding to a convex QP."
+function generate_random_parameter(rng; num_primals, num_inequalities, sparsity_rate)
     bernoulli = Distributions.Bernoulli(1 - sparsity_rate)
 
     M = let
@@ -24,85 +43,92 @@ function generate_test_problem(
         randn(rng, num_inequalities, num_primals) .*
         rand(rng, bernoulli, num_inequalities, num_primals)
     b = randn(rng, num_inequalities)
+    ϕ = randn(rng, num_primals)
 
-    G(x, y; θ) = M * x - θ - A' * y
-    H(x, y; θ) = A * x - b
-    K(z, θ) = begin
-        x = z[1:size(M, 1)]
-        y = z[(size(M, 1) + 1):end]
+    [reshape(M, length(M)); reshape(A, length(A)); b; ϕ]
+end
 
-        [G(x, y; θ); H(x, y; θ)]
-    end
+"Unpack a parameter vector θ into the components of a convex QP."
+function unpack_parameters(θ; num_primals, num_inequalities)
+    M = reshape(θ[1:(num_primals^2)], num_primals, num_primals)
+    A = reshape(
+        θ[(num_primals^2 + 1):(num_primals^2 + num_inequalities * num_primals)],
+        num_inequalities,
+        num_primals,
+    )
 
-    (; G, H, K)
+    b =
+        θ[(num_primals^2 + num_inequalities * num_primals + 1):(num_primals^2 + num_inequalities * (num_primals + 1))]
+    ϕ = θ[(num_primals^2 + num_inequalities * (num_primals + 1) + 1):end]
+
+    (; M, A, b, ϕ)
 end
 
 "Benchmark interior point solver against PATH on a bunch of random QPs."
 function benchmark(;
-    num_problems = 10,
-    num_samples_per_problem = 100,
+    num_samples = 1000,
     num_primals = 100,
     num_inequalities = 100,
     sparsity_rate = 0.9,
 )
     rng = Random.MersenneTwister(1)
 
-    # Generate random problems and parameters.
-    problems = @showprogress desc = "Generating test problems..." map(1:num_problems) do _
-        generate_test_problem(rng; num_primals, num_inequalities, sparsity_rate)
-    end
+    # Generate problem and random parameters.
+    @info "Generating random problems..."
+    problem = generate_test_problem(; num_primals, num_inequalities)
 
-    θs = map(1:num_samples_per_problem) do _
-        randn(rng, num_primals)
+    θs = map(1:num_samples) do _
+        generate_random_parameter(rng; num_primals, num_inequalities, sparsity_rate)
     end
 
     # Generate corresponding MCPs.
-    ip_mcps = @showprogress desc = "Generating IP MCPs... " map(problems) do p
-        MixedComplementarityProblems.PrimalDualMCP(
-            p.G,
-            p.H;
-            unconstrained_dimension = num_primals,
-            constrained_dimension = num_inequalities,
-            parameter_dimension = num_primals,
-        )
-    end
-
-    path_mcps = @showprogress desc = "Generating PATH MCPs..." map(problems) do p
-        lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
-        upper_bounds = fill(Inf, num_primals + num_inequalities)
-        ParametricMCPs.ParametricMCP(p.K, lower_bounds, upper_bounds, num_primals)
-    end
-
-    ip = @showprogress desc = "Solving IP MCPs..." map(ip_mcps) do mcp
-        # Make sure everything is compiled in a dry run.
-        MixedComplementarityProblems.solve(
+    @info "Generating IP MCP..."
+    parameter_dimension = length(first(θs))
+    ip_mcp = MixedComplementarityProblems.PrimalDualMCP(
+        problem.G,
+        problem.H;
+        unconstrained_dimension = num_primals,
+        constrained_dimension = num_inequalities,
+        parameter_dimension,
+    )
+
+    @info "Generating PATH MCP..."
+    lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
+    upper_bounds = fill(Inf, num_primals + num_inequalities)
+    path_mcp = ParametricMCPs.ParametricMCP(
+        problem.K,
+        lower_bounds,
+        upper_bounds,
+        parameter_dimension,
+    )
+
+    # Warm up the solvers.
+    @info "Warming up IP solver..."
+    MixedComplementarityProblems.solve(
+        MixedComplementarityProblems.InteriorPoint(),
+        ip_mcp,
+        first(θs),
+    )
+
+    @info "Warming up PATH solver..."
+    ParametricMCPs.solve(path_mcp, first(θs))
+
+    # Solve and time.
+    ip = @showprogress desc = "Solving IP MCPs..." map(θs) do θ
+        elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
             MixedComplementarityProblems.InteriorPoint(),
-            mcp,
-            zeros(num_primals),
+            ip_mcp,
+            θ,
         )
 
-        # Solve and time.
-        map(θs) do θ
-            elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
-                MixedComplementarityProblems.InteriorPoint(),
-                mcp,
-                θ,
-            )
-
-            (; elapsed_time, success = sol.status == :solved)
-        end
+        (; elapsed_time, success = sol.status == :solved)
     end
 
-    path = @showprogress desc = "Solving PATH MCPs..." map(path_mcps) do mcp
-        # Make sure everything is compiled in a dry run.
-        ParametricMCPs.solve(mcp, zeros(num_primals))
-
+    path = @showprogress desc = "Solving PATH MCPs..." map(θs) do θ
         # Solve and time.
-        map(θs) do θ
-            elapsed_time = @elapsed sol = ParametricMCPs.solve(mcp, θ)
+        elapsed_time = @elapsed sol = ParametricMCPs.solve(path_mcp, θ)
 
-            (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
-        end
+        (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
     end
 
     (; ip, path)
@@ -119,25 +145,15 @@ end
 
 "Estimate mean and standard deviation of runtimes for all problems."
 function runtime_stats(solver_data)
-    stats = map(solver_data) do problem_data
-        filtered_times = map(
-            datum -> datum.elapsed_time,
-            filter(datum -> datum.success, problem_data),
-        )
-        μ = Statistics.mean(filtered_times)
-        σ = Statistics.stdm(filtered_times, μ)
+    filtered_times =
+        map(datum -> datum.elapsed_time, filter(datum -> datum.success, solver_data))
+    μ = Statistics.mean(filtered_times)
+    σ = Statistics.stdm(filtered_times, μ)
 
-        (; μ, σ)
-    end
-
-    μ = map(datum -> datum.μ, stats)
-    σ = map(datum -> datum.σ, stats)
-    (; μ, σ, mean_μ = Statistics.mean(μ), mean_σ = Statistics.mean(σ))
+    (; μ, σ)
 end
 
 "Compute fraction of problems solved."
 function fraction_solved(solver_data)
-    Statistics.mean(solver_data) do problem_data
-        Statistics.mean(datum -> datum.success, problem_data)
-    end
+    Statistics.mean(datum -> datum.success, solver_data)
 end

From c5ae1fbc6a1f3d383c0bc5ee81b9c41c9c392250 Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 17:40:47 -0600
Subject: [PATCH 08/12] saving out mcps + adding support for solver kwargs

---
 benchmark/path.jl | 48 +++++++++++++++++++++++++++--------------------
 src/solver.jl     | 11 ++++++++---
 2 files changed, 36 insertions(+), 23 deletions(-)

diff --git a/benchmark/path.jl b/benchmark/path.jl
index ba63434..e12bb64 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -67,9 +67,12 @@ end
 "Benchmark interior point solver against PATH on a bunch of random QPs."
 function benchmark(;
     num_samples = 1000,
-    num_primals = 100,
-    num_inequalities = 100,
+    num_primals = 200,
+    num_inequalities = 200,
     sparsity_rate = 0.9,
+    ip_mcp = nothing,
+    path_mcp = nothing,
+    ip_kwargs = (;),
 )
     rng = Random.MersenneTwister(1)
 
@@ -84,37 +87,42 @@ function benchmark(;
     # Generate corresponding MCPs.
     @info "Generating IP MCP..."
     parameter_dimension = length(first(θs))
-    ip_mcp = MixedComplementarityProblems.PrimalDualMCP(
-        problem.G,
-        problem.H;
-        unconstrained_dimension = num_primals,
-        constrained_dimension = num_inequalities,
-        parameter_dimension,
-    )
+    ip_mcp =
+        !isnothing(ip_mcp) ? ip_mcp :
+        MixedComplementarityProblems.PrimalDualMCP(
+            problem.G,
+            problem.H;
+            unconstrained_dimension = num_primals,
+            constrained_dimension = num_inequalities,
+            parameter_dimension,
+        )
 
     @info "Generating PATH MCP..."
     lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
     upper_bounds = fill(Inf, num_primals + num_inequalities)
-    path_mcp = ParametricMCPs.ParametricMCP(
-        problem.K,
-        lower_bounds,
-        upper_bounds,
-        parameter_dimension,
-    )
+    path_mcp =
+        !isnothing(path_mcp) ? path_mcp :
+        ParametricMCPs.ParametricMCP(
+            problem.K,
+            lower_bounds,
+            upper_bounds,
+            parameter_dimension,
+        )
 
     # Warm up the solvers.
     @info "Warming up IP solver..."
     MixedComplementarityProblems.solve(
         MixedComplementarityProblems.InteriorPoint(),
         ip_mcp,
-        first(θs),
+        first(θs);
+        ip_kwargs...,
     )
 
     @info "Warming up PATH solver..."
     ParametricMCPs.solve(path_mcp, first(θs))
 
     # Solve and time.
-    ip = @showprogress desc = "Solving IP MCPs..." map(θs) do θ
+    ip_data = @showprogress desc = "Solving IP MCPs..." map(θs) do θ
         elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
             MixedComplementarityProblems.InteriorPoint(),
             ip_mcp,
@@ -124,14 +132,14 @@ function benchmark(;
         (; elapsed_time, success = sol.status == :solved)
     end
 
-    path = @showprogress desc = "Solving PATH MCPs..." map(θs) do θ
+    path_data = @showprogress desc = "Solving PATH MCPs..." map(θs) do θ
         # Solve and time.
         elapsed_time = @elapsed sol = ParametricMCPs.solve(path_mcp, θ)
 
         (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
     end
 
-    (; ip, path)
+    (; ip_mcp, path_mcp, ip_data, path_data)
 end
 
 "Compute summary statistics from solver benchmark data."
@@ -140,7 +148,7 @@ function summary_statistics(data)
         (; success_rate = fraction_solved(solver_data), runtime_stats(solver_data)...)
     end
 
-    (; ip = accumulate_stats(data.ip), path = accumulate_stats(data.path))
+    (; ip = accumulate_stats(data.ip_data), path = accumulate_stats(data.path_data))
 end
 
 "Estimate mean and standard deviation of runtimes for all problems."
diff --git a/src/solver.jl b/src/solver.jl
index 1ba0410..64f2d09 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -25,6 +25,9 @@ function solve(
     tol = 1e-4,
     max_inner_iters = 20,
     max_outer_iters = 50,
+    tightening_rate = 1,
+    loosening_rate = 1,
+    verbose = false
 )
     # Set up common memory.
     ∇F = mcp.∇F_z!.result_buffer
@@ -52,14 +55,14 @@ function solve(
             # TODO! Can add some adaptive regularization.
             mcp.F!(F, x, y, s; θ, ϵ)
             mcp.∇F_z!(∇F, x, y, s; θ, ϵ)
-            δz .= (∇F \ F) .* -1
+            δz .= ((∇F + tol * I) \ F) .* -1
 
             # Fraction to the boundary linesearch.
             α_s = fraction_to_the_boundary_linesearch(s, δs; tol)
             α_y = fraction_to_the_boundary_linesearch(y, δy; tol)
 
             if isnan(α_s) || isnan(α_y)
-                @warn "Linesearch failed. Exiting prematurely."
+                verbose && @warn "Linesearch failed. Exiting prematurely."
                 status = :failed
                 break
             end
@@ -73,7 +76,9 @@ function solve(
             inner_iters += 1
         end
 
-        ϵ *= (status == :solved) ? 1 - exp(-inner_iters) : 1 + exp(-inner_iters)
+        ϵ *=
+            (status == :solved) ? 1 - exp(-tightening_rate * inner_iters) :
+            1 + exp(-loosening_rate * inner_iters)
         outer_iters += 1
     end
 

From 00a4022e0ccb8a3df752a996486b37ff5e23d74d Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 17:46:27 -0600
Subject: [PATCH 09/12] [skip ci] formatting

---
 src/solver.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solver.jl b/src/solver.jl
index 64f2d09..1db999d 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -27,7 +27,7 @@ function solve(
     max_outer_iters = 50,
     tightening_rate = 1,
     loosening_rate = 1,
-    verbose = false
+    verbose = false,
 )
     # Set up common memory.
     ∇F = mcp.∇F_z!.result_buffer

From 733d4afd5eac2e84c9bc4b59605efb3534f7123f Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 18:15:07 -0600
Subject: [PATCH 10/12] [skip ci] smaller test problems

---
 benchmark/path.jl | 4 ++--
 src/solver.jl     | 5 +++--
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/benchmark/path.jl b/benchmark/path.jl
index e12bb64..1c7cc06 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -67,8 +67,8 @@ end
 "Benchmark interior point solver against PATH on a bunch of random QPs."
 function benchmark(;
     num_samples = 1000,
-    num_primals = 200,
-    num_inequalities = 200,
+    num_primals = 100,
+    num_inequalities = 100,
     sparsity_rate = 0.9,
     ip_mcp = nothing,
     path_mcp = nothing,
diff --git a/src/solver.jl b/src/solver.jl
index 1db999d..f02d3d6 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -27,6 +27,7 @@ function solve(
     max_outer_iters = 50,
     tightening_rate = 1,
     loosening_rate = 1,
+    min_stepsize = 1e-2,
     verbose = false,
 )
     # Set up common memory.
@@ -58,8 +59,8 @@ function solve(
             δz .= ((∇F + tol * I) \ F) .* -1
 
             # Fraction to the boundary linesearch.
-            α_s = fraction_to_the_boundary_linesearch(s, δs; tol)
-            α_y = fraction_to_the_boundary_linesearch(y, δy; tol)
+            α_s = fraction_to_the_boundary_linesearch(s, δs; tol = min_stepsize)
+            α_y = fraction_to_the_boundary_linesearch(y, δy; tol = min_stepsize)
 
             if isnan(α_s) || isnan(α_y)
                 verbose && @warn "Linesearch failed. Exiting prematurely."

From 8b5d018a6414fef4568942ab25ad217e76e052cd Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 18:19:30 -0600
Subject: [PATCH 11/12] [skip ci] hiding path warnings

---
 benchmark/path.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/benchmark/path.jl b/benchmark/path.jl
index 1c7cc06..9f5e075 100644
--- a/benchmark/path.jl
+++ b/benchmark/path.jl
@@ -119,7 +119,7 @@ function benchmark(;
     )
 
     @info "Warming up PATH solver..."
-    ParametricMCPs.solve(path_mcp, first(θs))
+    ParametricMCPs.solve(path_mcp, first(θs); warn_on_convergence_failure = false)
 
     # Solve and time.
     ip_data = @showprogress desc = "Solving IP MCPs..." map(θs) do θ
@@ -134,7 +134,8 @@ function benchmark(;
 
     path_data = @showprogress desc = "Solving PATH MCPs..." map(θs) do θ
         # Solve and time.
-        elapsed_time = @elapsed sol = ParametricMCPs.solve(path_mcp, θ)
+        elapsed_time = @elapsed sol =
+            ParametricMCPs.solve(path_mcp, θ; warn_on_convergence_failure = false)
 
         (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
     end

From d9842fe6a97d091ccbf77f1726ce6230d8d2d912 Mon Sep 17 00:00:00 2001
From: dfridovi <dfk@utexas.edu>
Date: Tue, 24 Dec 2024 18:23:07 -0600
Subject: [PATCH 12/12] tuning rates a bit

---
 src/solver.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/solver.jl b/src/solver.jl
index f02d3d6..e43db3f 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -25,8 +25,8 @@ function solve(
     tol = 1e-4,
     max_inner_iters = 20,
     max_outer_iters = 50,
-    tightening_rate = 1,
-    loosening_rate = 1,
+    tightening_rate = 0.1,
+    loosening_rate = 0.5,
     min_stepsize = 1e-2,
     verbose = false,
 )