Decrease allocations

Project.toml

@@ -10,11 +10,15 @@ FameSVD = "9ba2d756-9ce3-11e9-1a71-0ffcb019784d"
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
@@ -25,8 +29,12 @@ DocumenterTools = "0.1.21"
 FameSVD = "0.1"
 GenericLinearAlgebra = "0.3"
 JuMP = "1"
+LowRankOpt = "0.2.1"
 MathOptInterface = "1"
 MultiFloats = "2"
+MutableArithmetics = "1.6.4"
+NLPModels = "0.21.5"
 Revise = "3"
+SolverCore = "0.3.8"
 TimerOutputs = "0.5"
 julia = "1"

docs/src/low-rank_data.md

@@ -9,9 +9,4 @@ with data matrices ``B_i\in{\mathbb R}^{m\times k}`` and with ``k \ll n``.
 The complexity of ``H`` assembly is then reduced from ``nm^3`` to ``knm^2``.
 
 *In particular:*
-If we know that ``A_i`` have rank one, the decomposition ``A_i = b_i b_i^\top`` is performed by Loraine automatically (`datarank = -1`).
-
-The model is represented internally via the following `struct`:
-```@docs
-Loraine.MyModel
-```
+If we know that ``A_i`` have rank one, the decomposition ``A_i = b_i b_i^\top`` is performed by [LowRankOpt](https://github.com/blegat/LowRankOpt.jl/) automatically (`detect_rank_1_tol = 1e-6`).

examples/benchmark_sdplib.jl

@@ -21,10 +21,7 @@
 #
 # Only options common to both `main` and `nlpmodel` are used, so the same
 # script produces comparable CSVs on either branch. The CSV is written with a
-# plain `print` (no CSV.jl dependency).
-#
-# The active project must have `Chairmarks` (used on the worker for a robust
-# minimum solve time) in addition to `JuMP` and `Loraine`.
+# plain `print` (no CSV.jl dependency), and it needs only `JuMP` and `Loraine`.
 
 using Distributed
 import Printf
@@ -63,15 +60,11 @@ end
 function start_worker()
     pid = addprocs(1; exeflags = `--project=$(Base.active_project())`)[1]
     ospid = remotecall_fetch(getpid, pid)
-    # Load in a separate eval from the function definition below: `@b` must be
-    # macro-expanded *after* `using Chairmarks` has taken effect.
     remotecall_wait(Core.eval, pid, Main, quote
         using JuMP
         import Loraine
-        using Chairmarks
-    end)
-    remotecall_wait(Core.eval, pid, Main, quote
-        function solve_problem(path, kit)
+        function solve_problem(path, kit, timeout)
+            t0 = time()
             model = read_from_file(path)
             set_optimizer(model, Loraine.Optimizer{Float64})
             set_attribute(model, "kit", kit)
@@ -80,28 +73,34 @@ function start_worker()
             set_attribute(model, "maxit", 100)
             set_attribute(model, "datasparsity", 8)
             optimize!(model) # warm up: absorb this problem's compilation
-            t = @elapsed optimize!(model)
-            # Millisecond solves are noise-dominated (GC/scheduler jitter), so
-            # take a robust minimum over repeated warm re-solves. Multi-second
-            # solves are measured once: they aren't noisy, and re-running them
-            # would risk blowing past the per-problem timeout.
-            if t < 0.5
-                t = (@b optimize!(model) seconds = 1).time
-            end
-            return (
-                t,
-                solve_time(model), # solver's own elapsed time (no harness)
+            sample() = (
+                solve_time(model), # solver's own `SolveTimeSec`
                 barrier_iterations(model),
                 objective_value(model),
                 string(termination_status(model)),
             )
+            # Re-solve for more samples (one row per solve, so `merge` can take
+            # the minimum), but only while another solve fits the budget (85% of
+            # the timeout). This keeps the whole call under the per-problem
+            # timeout while letting a larger `timeout` buy more samples; a solve
+            # slower than ~half the timeout just yields a single sample.
+            wall = @elapsed optimize!(model)
+            samples = [(wall, sample()...)]
+            # The 32-sample cap only binds for fast problems (slow ones are
+            # limited by the budget, so a larger `timeout` buys them more rows).
+            budget = 0.85 * timeout
+            while length(samples) < 32 && (time() - t0) + wall < budget
+                w = @elapsed optimize!(model)
+                push!(samples, (w, sample()...))
+            end
+            return samples
         end
     end)
     return pid, ospid
 end
 
-launch(pid, path, kit) =
-    remotecall((p, k) -> Main.solve_problem(p, k), pid, path, kit)
+launch(pid, path, kit, timeout) =
+    remotecall((p, k, t) -> Main.solve_problem(p, k, t), pid, path, kit, timeout)
 
 # Solve with a wall-clock cap. `isready` on a *remote* Future blocks while the
 # worker is busy (it has to query the worker), so we instead let an async task
@@ -110,7 +109,7 @@ launch(pid, path, kit) =
 function solve_capped(pid, path, kit, timeout)
     ch = Channel{Any}(1)
     @async put!(ch, try
-        fetch(launch(pid, path, kit))
+        fetch(launch(pid, path, kit, timeout))
     catch err
         err
     end)
@@ -121,16 +120,21 @@ function solve_capped(pid, path, kit, timeout)
     return (:timeout, nothing)
 end
 
-function bench(;
-    out = get(ARGS, 1, "sdplib_results.csv"),
-    label = get(ARGS, 2, "loraine"),
+function bench(
+    label,
+    only = get(ENV, "SDPLIB_ONLY", nothing);
+    out = label * ".csv",
     maxn = parse(Int, get(ENV, "SDPLIB_MAXN", "1000")),
     maxm = parse(Int, get(ENV, "SDPLIB_MAXM", "3000")),
     kit = parse(Int, get(ENV, "SDPLIB_KIT", "0")),
     timeout = parse(Float64, get(ENV, "SDPLIB_TIMEOUT", "60")),
 )
-    only = haskey(ENV, "SDPLIB_ONLY") ? Set(split(ENV["SDPLIB_ONLY"], ",")) :
-           nothing
+    # `only`: `nothing` (all problems within the size caps), a single name or a
+    # comma-separated string (e.g. "truss1" or "truss1,qap8"), or any iterable
+    # of names.
+    only =
+        isnothing(only) ? nothing :
+        only isa AbstractString ? Set(split(only, ",")) : Set(string.(only))
     ref = read_reference()
 
     problems = String[]
@@ -149,71 +153,80 @@ function bench(;
     pid, ospid = start_worker()
     # Warm up compilation on the smallest problem so timings exclude it.
     warmup = joinpath(DATA, problems[1] * ".dat-s")
-    fetch(launch(pid, warmup, kit))
+    fetch(launch(pid, warmup, kit, timeout))
 
-    open(out, "w") do io
-        println(
-            io,
-            "label,problem,m,n,status,time_s,solve_time_s,iterations,objective,optimal,rel_gap",
-        )
+    # Append so re-running accumulates more rows (one solve = one row); the
+    # minimum over rows is taken later, in `merge_results.jl`.
+    write_header = !isfile(out) || filesize(out) == 0
+    open(out, "a") do io
+        if write_header
+            println(
+                io,
+                "label,problem,m,n,status,time_s,solve_time_s,iterations,objective,optimal,rel_gap",
+            )
+        end
         ntot = length(problems)
         for (idx, name) in enumerate(problems)
             r = ref[name]
-            status, t, st, iters, obj = "ok", NaN, NaN, -1, NaN
             outcome, payload =
                 solve_capped(pid, joinpath(DATA, name * ".dat-s"), kit, timeout)
-            if outcome === :ok
-                t, st, iters, obj, status = payload
+            # `samples`: one `(wall, solve_time, iters, obj, status)` per solve.
+            # Error/timeout collapse to a single synthetic sample.
+            samples = if outcome === :ok
+                payload
             elseif outcome === :error
-                status = "ERROR: " * sprint(showerror, payload)[1:min(end, 60)]
+                msg = "ERROR: " * sprint(showerror, payload)[1:min(end, 60)]
+                [(NaN, NaN, -1, NaN, msg)]
             else
                 # Stuck: SIGKILL the worker's OS process, then respawn a fresh
                 # one (and re-warm it) for the remaining problems.
                 run(`kill -9 $ospid`)
                 rmprocs(pid; waitfor = 0)
-                status, t = "TIMEOUT", timeout
                 pid, ospid = start_worker()
-                fetch(launch(pid, warmup, kit))
+                fetch(launch(pid, warmup, kit, timeout))
+                [(timeout, NaN, -1, NaN, "TIMEOUT")]
             end
-            gap = isfinite(obj) && isfinite(r.opt) ?
-                  abs(obj - r.opt) / max(1, abs(r.opt)) : NaN
-            println(
-                io,
-                join(
-                    [
-                        label,
-                        name,
-                        r.m,
-                        r.n,
-                        "\"$status\"",
-                        round(t, sigdigits = 5),
-                        round(st, sigdigits = 5),
-                        iters,
-                        obj,
-                        r.opt,
-                        gap,
-                    ],
-                    ",",
-                ),
-            )
+            for (t, st, iters, obj, status) in samples
+                gap = isfinite(obj) && isfinite(r.opt) ?
+                      abs(obj - r.opt) / max(1, abs(r.opt)) : NaN
+                println(
+                    io,
+                    join(
+                        [
+                            label,
+                            name,
+                            r.m,
+                            r.n,
+                            "\"$status\"",
+                            round(t, sigdigits = 5),
+                            round(st, sigdigits = 5),
+                            iters,
+                            obj,
+                            r.opt,
+                            gap,
+                        ],
+                        ",",
+                    ),
+                )
+            end
+            flush(io)
+            # console: minimum solve time over the samples, and the count
+            fin = filter(isfinite, [s[2] for s in samples])
             Printf.@printf(
-                "[%2d/%2d] %-12s n=%-5d m=%-5d  %8.3fs (solve %8.3fs, %3d it)  obj=%-14.6g %s\n",
+                "[%2d/%2d] %-12s n=%-5d m=%-5d  min solve %8.3fs (%2d runs, %3d it)  obj=%-14.6g %s\n",
                 idx,
                 ntot,
                 name,
                 r.n,
                 r.m,
-                t,
-                st,
-                iters,
-                obj,
-                status,
+                isempty(fin) ? NaN : minimum(fin),
+                length(samples),
+                samples[end][3],
+                samples[end][4],
+                samples[end][5],
             )
-            flush(io)
         end
     end
     rmprocs(pid)
     println("\nWrote ", out)
 end
-
-bench()

examples/k.jl

@@ -5,7 +5,7 @@ import Loraine
 using MultiFloats
 
 # model = Model(Loraine.Optimizer)
-model = Model(Loraine.Optimizer{Float64x2})
+model = GenericModel{Float64x2}(Loraine.Optimizer{Float64x2})
 # model = JuMP.GenericModel{Float64x2}(Loraine.Optimizer{Float64x2})
 
 # @variable(model, x >= 0)
@@ -30,6 +30,7 @@ using Test
 @test primal_status(model) == MOI.FEASIBLE_POINT
 @test dual_status(model) == MOI.FEASIBLE_POINT
 @test objective_value(model) ≈ 4 rtol = 1e-6
+@test dual_objective_value(model) ≈ 4 rtol = 1e-6
 
 @test value(x) ≈ 2 rtol = 1e-6
 

examples/merge_results.jl

@@ -14,9 +14,18 @@ using CSV
 using DataFrames
 import Printf
 
+# Load a CSV that may hold several rows per problem (from repeated/appended
+# runs) and collapse each problem to its best run — the row with the smallest
+# finite `solve_time_s` (falling back to the first row if none finished).
 function load(path)
     df = CSV.read(path, DataFrame)
-    return df, isempty(df.label) ? basename(path) : first(df.label)
+    label = isempty(df.label) ? basename(path) : first(df.label)
+    best = combine(groupby(df, :problem)) do sub
+        fin = findall(x -> !ismissing(x) && isfinite(x), sub.solve_time_s)
+        i = isempty(fin) ? 1 : fin[argmin(sub.solve_time_s[fin])]
+        sub[i, :]
+    end
+    return best, label
 end
 
 function fmt_time(t)
@@ -43,22 +52,23 @@ function merge(baseline = ARGS[1], comparison = ARGS[2])
     base, base_label = load(baseline)
     comp, comp_label = load(comparison)
 
-    # Compare on `time_s`: the robust (Chairmarks minimum) warm-solve time,
-    # free of compilation and one-shot noise.
+    # Compare on `solve_time_s`: the solver's own `SolveTimeSec` (minimum over
+    # the repeated warm re-solves), free of the worker harness. `time_s` (wall)
+    # is kept in the CSV as a cross-check.
     a = select(
         base,
         :problem,
         :n,
         :m,
-        :time_s => :t_base,
+        :solve_time_s => :t_base,
         :iterations => :it_base,
         :objective => :obj_base,
         :status => :st_base,
     )
     b = select(
         comp,
         :problem,
-        :time_s => :t_comp,
+        :solve_time_s => :t_comp,
         :iterations => :it_comp,
         :objective => :obj_comp,
         :status => :st_comp,

examples/simple_linear.jl

@@ -0,0 +1,27 @@
+model = Model(Loraine.Optimizer)
+@variable(model, x)
+@constraint(model, con_ref, x <= 2)
+@objective(model, Max, 3x)
+optimize!(model)
+@test value(x) ≈ 2 rtol = 1e-6
+@test objective_value(model) ≈ 6 rtol = 1e-6
+@test termination_status(model) == MOI.OPTIMAL
+@test primal_status(model) == MOI.FEASIBLE_POINT
+@test dual_status(model) == MOI.FEASIBLE_POINT
+@test dual(con_ref) ≈ -3 rtol = 1e-6
+
+model = Model(Loraine.Optimizer)
+@variable(model, x)
+@variable(model, y)
+@constraint(model, cx, x >= 0)
+@constraint(model, cy, y <= 0)
+@objective(model, Min, x - y)
+optimize!(model)
+@test value(x) ≈ 0 atol = 1e-6
+@test value(y) ≈ 0 atol = 1e-6
+@test objective_value(model) ≈ 0 atol = 1e-6
+@test termination_status(model) == MOI.OPTIMAL
+@test primal_status(model) == MOI.FEASIBLE_POINT
+@test dual_status(model) == MOI.FEASIBLE_POINT
+@test dual(cx) ≈ 1 rtol = 1e-6
+@test dual(cy) ≈ -1 rtol = 1e-6

examples/solve_sdpa.jl

@@ -62,6 +62,19 @@ using Test
 @test objective_value(model) ≈ 23 rtol = 1e-6
 # # value.(X)
 
+# With CG now
+set_attribute(model, "kit", 1)
+optimize!(model)
+@test objective_value(model) ≈ 23 rtol = 1e-6
+
+set_attribute(model, "preconditioner", 0)
+optimize!(model)
+@test objective_value(model) ≈ 23 rtol = 1e-6
+
+set_attribute(model, "preconditioner", 2)
+optimize!(model)
+@test objective_value(model) ≈ 23 rtol = 1e-6
+
 # Mosek (CSDP, etc) for a comparison
 # Mosek must solve the dualized problem to be efficient