Use an interface instead of accessing internal fields

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"

README.md

@@ -68,11 +68,11 @@ includes a few examples of how to use Loraine via JuMP; in particular,
 A few sample problems can be found in folder `examples/data`.
 
 ## Rank-one data
-If the solution does not have low rank, it is recommended to use a direct 
-solver `kit = 0`. However, if you know that your data matrices are all rank-one, 
-use the option `datarank = -1` to get a significant reduction in the complexity 
-(and CPU time). Examples of such problems are `maxG11` and `thetaG11` from the 
-`SDPLIB` collection.
+If the solution does not have low rank, it is recommended to use a direct
+solver `kit = 0`. If your data matrices are all rank-one, the factorization is
+detected by [LowRankOpt](https://github.com/blegat/LowRankOpt.jl/) and used to
+get a significant reduction in the complexity (and CPU time). Examples of such
+problems are `maxG11` and `thetaG11` from the `SDPLIB` collection.
 
 ## Documentation
 [Loraine documentation](https://kocvara.github.io/Loraine.jl/)
@@ -90,9 +90,6 @@ preconditioner  # 0...no; 1...H_alpha; 2...H_beta; 4...hybrid [1]
 erank           # estimated rank [1]
 aamat           # 0..A^TA; 1..diag(A^TA); 2..identity [2]
 verb            # 2..full output; 1..short output; 0..no output [1]
-datarank        # 0..full rank matrices expected [0]
-                # -1..rank-1 matrices expected, converted to vectors, if possible
-                # (TBD) 1..vectors expected for low-rank data matrices
 initpoint       # 0..Loraine heuristics, 1..SDPT3-like heuristics [0]
 timing          # 1..yes, 0..no
 maxit           # maximal number of global iterations [200]

docs/src/Loraine_options.md

@@ -11,9 +11,6 @@ preconditioner  # 0...no; 1...H_alpha; 2...H_beta; 4...hybrid [1]
 erank           # estimated rank [1]
 aamat           # 0..A^TA; 1..diag(A^TA); 2..identity [2]
 verb            # 2..full output; 1..short output; 0..no output [1]
-datarank        # 0..full rank matrices expected [0]
-                # -1..rank-1 matrices expected, converted to vectors, if possible
-                # (TBD) 1..vectors expected for low-rank data matrices
 initpoint       # 0..Loraine heuristics, 1..SDPT3-like heuristics [0]
 timing          # 1..yes, 0..no
 maxit           # maximal number of global iterations [200]
@@ -38,11 +35,6 @@ datasparsity    # data matrices treated as sparse when number of their
     - if you are not sure what the actual rank of the solution is, **always choose** `erank = 1`; with inreasing value of `erank`, the complexity of the preconditioner grows and the whole code could be slower, despite needing fewer CG iterations
     - only if you are sure about the rank of the solution, set `erank` to this value (but you should always compare it to `erank = 1`)
 
-- `datarank` (only used with the direct solver `kit = 0`)
-    - choose `datarank = -1` if you know (or suspect) that all the data matrices ``A_i`` have rank one; in this case, the matrices will be factorized as ``A_i = b_i b_i^T`` and vectors ``b_i`` will be used when constructing the Schur complement matrix
-    - if you are not sure about the rank of the data matrices, you can always try to set `datarank = -1`; if the factorization of any matrix fails, Lorain will switch to the default option `datarank = 0`
-    - for rank-one data matrices, option `datarank = -1` will result in a much faster code than the default `datarank = 0`
-
 - `timing` is not used when Loraine is called from JuMP
 
 - `tol_cg, tol_cg_up, aamat`: it is not recommended to change values of these options, unless you really want to

docs/src/index.md

@@ -48,6 +48,19 @@ set_optimizer(model, Loraine.Optimizer)
 optimize!(model)
 ```
 
+## Use without JuMP
+
+An SDP problem stored in SDPA format can also be solved without JuMP with
+[`loraine`](@ref):
+```julia
+using Loraine
+model = loraine("examples/data/theta1.dat-s", Dict("kit" => 1))
+```
+
+```@docs
+loraine
+```
+
 ## License and Original Contributors
 
 Loraine is licensed under the [MIT License](https://github.com/kocvara/Loraine.jl/blob/main/LICENSE.md).

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`).

docs/src/overview.md

@@ -69,4 +69,4 @@ Further assumptions
 - *At the moment, Loraine can only handle rank-one data.*
 - *This feature is only relevant for Loraine used with the **direct solver**.*
 
-If you know (or strongly suspect) that *all* data matrices ``A_i`` have rank one, select the option `datarank = -1`. Loraine will factorize the matrices as ``A_i = b_i b_i^\top`` and use only the vectors ``b_i`` in the interior point algorithm. This will gravely reduce the complexity (and the elapsed time) of Loraine.
+If *all* data matrices ``A_i`` have rank one, the factorization ``A_i = b_i b_i^\top`` is detected by [LowRankOpt](https://github.com/blegat/LowRankOpt.jl/) and only the vectors ``b_i`` are used in the interior point algorithm. This gravely reduces the complexity (and the elapsed time) of Loraine. See [Low-rank data, more details](@ref) for details.

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