Enh/options

README.md

@@ -54,7 +54,7 @@ solution = SCS_solve(m, n, A, ..., psize; max_iters=10, verbose=0);
 m = solve!(problem, SCSSolver(max_iters=10, verbose=0))
 ```
 
-Moreover, You may select one of the linear solvers to be used by `SCSSolver` via `linearsolver` keyword. The options available are `SCS.Indirect` (the default) and `SCS.Direct`.
+Moreover, You may select one of the linear solvers to be used by `SCSSolver` via `linear_solver` keyword. The options available are `SCS.Indirect` (the default) and `SCS.Direct`.
 
 ### High level wrapper
 

deps/.gitignore

@@ -1,2 +1,3 @@
 deps.jl
 usr
+build.log

src/MOIWrapper.jl

@@ -60,9 +60,11 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     maxsense::Bool
     data::Union{Nothing, ModelData} # only non-Void between MOI.copy_to and MOI.optimize!
     sol::MOISolution
-    function Optimizer()
-        new(ConeData(), false, nothing, MOISolution())
+    options
+    function Optimizer(; options...)
+        new(ConeData(), false, nothing, MOISolution(), options)
     end
+
 end
 
 function MOI.is_empty(optimizer::Optimizer)
@@ -319,10 +321,14 @@ function MOI.optimize!(optimizer::Optimizer)
     objconstant = optimizer.data.objconstant
     c = optimizer.data.c
     optimizer.data = nothing # Allows GC to free optimizer.data before A is loaded to SCS
-    sol = SCS_solve(SCS.Indirect, m, n, A, b, c,
+
+    linear_solver, options = sanatize_SCS_options(optimizer.options)
+
+    sol = SCS_solve(linear_solver, m, n, A, b, c, 
                     cone.f, cone.l, cone.qa, cone.sa, cone.ep, cone.ed, cone.p,
                     optimizer.sol.primal, optimizer.sol.dual,
-                    optimizer.sol.slack)
+                    optimizer.sol.slack; options...)
+
     ret_val = sol.ret_val
     primal = sol.x
     dual = sol.y

src/MPBWrapper.jl

@@ -93,17 +93,9 @@ end
 =#
 
 function optimize!(m::SCSMathProgModel)
+    linear_solver, options = sanatize_SCS_options(m.options)
     t = time()
-
-    T = SCS.Indirect # the default method
-    options = m.options
-    opts = Dict(m.options)
-    if :linearsolver in keys(opts)
-        T = opts[:linearsolver]
-        options = [(k,v) for (k,v) in options if k !=:linearsolver]
-    end
-
-    solution = SCS_solve(T, m.m, m.n, m.A, m.b, m.c, m.f, m.l, m.q,
+    solution = SCS_solve(linear_solver, m.m, m.n, m.A, m.b, m.c, m.f, m.l, m.q,
                          m.s, m.ep, m.ed, Float64[],
                          m.primal_sol, m.dual_sol, m.slack; options...)
     m.solve_time = time() - t

src/c_wrapper.jl

@@ -47,8 +47,9 @@ function SCS_solve(T::Union{Type{Direct}, Type{Indirect}},
 
     managed_matrix = ManagedSCSMatrix(m, n, A)
     matrix = Ref(SCSMatrix(managed_matrix))
-    settings = Ref(SCSSettings(;options...))
+    settings = Ref(SCSSettings(T; options...))
     data = Ref(SCSData(m, n, Base.unsafe_convert(Ptr{SCSMatrix}, matrix), pointer(b), pointer(c), Base.unsafe_convert(Ptr{SCSSettings},settings)))
+
     cone = Ref(SCSCone(f, l, q, s, ep, ed, p))
     info = Ref(SCSInfo())
 
@@ -86,6 +87,11 @@ end
 # Take Ref{}s because SCS might modify the structs
 for (T, lib) in zip([SCS.Direct, SCS.Indirect], [SCS.direct, SCS.indirect])
     @eval begin
+
+        function SCS_set_default_settings(::Type{$T}, data::Ref{SCSData})
+            ccall((:scs_set_default_settings, $lib), Nothing, (Ref{SCSData}, ), data)
+        end
+
         function SCS_init(::Type{$T}, data::Ref{SCSData}, cone::Ref{SCSCone}, info::Ref{SCSInfo})
 
             p_work = ccall((:scs_init, $lib), Ptr{Nothing},

src/types.jl

@@ -38,24 +38,49 @@ SCSMatrix(m::ManagedSCSMatrix) =
 
 
 struct SCSSettings
-    normalize::Int # boolean, heuristic data rescaling: 1
-    scale::Cdouble # if normalized, rescales by this factor: 1
-    rho_x::Cdouble # x equality constraint scaling: 1e-3
-    max_iters::Int # maximum iterations to take: 5000
-    eps::Cdouble # convergence tolerance: 1e-5
-    alpha::Cdouble # relaxation parameter: 1.5
-    cg_rate::Cdouble # for indirect, tolerance goes down like (1/iter)^cg_rate: 2
-    verbose::Int # boolean, write out progress: 1
-    warm_start::Int # boolean, warm start (put initial guess in Sol struct): 0
-    acceleration_lookback::Int # acceleration memory parameter: 20
+    normalize::Int # boolean, heuristic data rescaling
+    scale::Cdouble # if normalized, rescales by this factor
+    rho_x::Cdouble # x equality constraint scaling
+    max_iters::Int # maximum iterations to take
+    eps::Cdouble # convergence tolerance
+    alpha::Cdouble # relaxation parameter
+    cg_rate::Cdouble # for indirect, tolerance goes down like (1/iter)^cg_rate
+    verbose::Int # boolean, write out progress
+    warm_start::Int # boolean, warm start (put initial guess in Sol struct)
+    acceleration_lookback::Int # acceleration memory parameter
+
+    SCSSettings() = new()
+    SCSSettings(normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose, warm_start, acceleration_lookback) = new(normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose, warm_start, acceleration_lookback)
 end
 
-function SCSSettings(;normalize=1::Int, scale=convert(Cdouble, 1.0)::Cdouble, rho_x=convert(Cdouble,1e-3)::Cdouble,
-                        max_iters=5000::Int, eps=convert(Cdouble, 1e-5)::Cdouble, alpha=convert(Cdouble, 1.5)::Cdouble,
-                        cg_rate=convert(Cdouble,2)::Cdouble, verbose=1::Int, warm_start=0::Int, acceleration_lookback=20::Int)
-    return SCSSettings(normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose, warm_start, acceleration_lookback)
+struct Direct end
+struct Indirect end
+
+function _SCS_user_settings(default_settings::SCSSettings;
+        normalize=default_settings.normalize,
+        scale=default_settings.scale,
+        rho_x=default_settings.rho_x,
+        max_iters=default_settings.max_iters,
+        eps=default_settings.eps,
+        alpha=default_settings.alpha,
+        cg_rate=default_settings.cg_rate,
+        verbose=default_settings.verbose,
+        warm_start=default_settings.warm_start,
+        acceleration_lookback=default_settings.acceleration_lookback)
+    return SCSSettings(normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose,warm_start, acceleration_lookback)
 end
 
+function SCSSettings(linear_solver::Union{Type{Direct}, Type{Indirect}}; options...)
+
+    mmatrix = ManagedSCSMatrix(0,0,spzeros(1,1))
+    matrix = Ref(SCSMatrix(mmatrix))
+    default_settings = Ref(SCSSettings())
+    dummy_data = Ref(SCSData(0,0, Base.unsafe_convert(Ptr{SCSMatrix}, matrix),
+        pointer([0.0]), pointer([0.0]),
+        Base.unsafe_convert(Ptr{SCSSettings}, default_settings)))
+    SCS_set_default_settings(linear_solver, dummy_data)
+    return _SCS_user_settings(default_settings[]; options...)
+end
 
 struct SCSData
     # A has m rows, n cols
@@ -156,5 +181,25 @@ mutable struct Solution
     end
 end
 
-struct Direct end
-struct Indirect end
+function sanatize_SCS_options(options)
+    options = Dict(options)
+    if :linear_solver in keys(options)
+        linear_solver = options[:linear_solver]
+        if linear_solver == Direct || linear_solver == Indirect
+            nothing
+        else
+            throw(ArgumentError("Unrecognized linear_solver passed to SCS: $linear_solver;\nRecognized options are: $Direct, $Indirect."))
+        end
+        delete!(options, :linear_solver)
+    else
+        linear_solver = Indirect # the default linear_solver
+    end
+
+    SCS_options = fieldnames(SCSSettings)
+    unrecognized = setdiff(keys(options), SCS_options)
+    if length(unrecognized) > 0
+        plur = length(unrecognized) > 1 ? "s" : ""
+        throw(ArgumentError("Unrecognized option$plur passed to SCS: $(join(unrecognized, ", "));\nRecognized options are: $(join(SCS_options, ", ", " and "))."))
+    end
+    return linear_solver, options
+end

test/MOIWrapper.jl

@@ -14,17 +14,16 @@ MOIU.@model(SCSModelData,
             (MOI.ScalarAffineFunction,),
             (MOI.VectorOfVariables,),
             (MOI.VectorAffineFunction,))
-const optimizer = MOIU.CachingOptimizer(SCSModelData{Float64}(), SCS.Optimizer())
+for T in [SCS.Direct, SCS.Indirect]
+    optimizer = MOIU.CachingOptimizer(SCSModelData{Float64}(), SCS.Optimizer(linear_solver=T, eps=1e-8))
 
-# linear9test needs 1e-3 with SCS < 2.0 and 5e-1 with SCS 2.0
-# linear2test needs 1e-4
-const config = MOIT.TestConfig(atol=1e-4, rtol=1e-4)
+    config = MOIT.TestConfig(atol=1e-5)
 
-@testset "Continuous linear problems" begin
-    # AlmostSuccess for linear9 with SCS 2
-    MOIT.contlineartest(MOIB.SplitInterval{Float64}(optimizer), config, ["linear9"])
-end
+    @testset "Continuous linear problems with $T" begin
+        MOIT.contlineartest(MOIB.SplitInterval{Float64}(optimizer), config)
+    end
 
-@testset "Continuous conic problems" begin
-    MOIT.contconictest(MOIB.RootDet{Float64}(MOIB.LogDet{Float64}(optimizer)), config, ["rsoc", "geomean", "psds", "rootdet", "logdets"])
+    @testset "Continuous conic problems with $T" begin
+        MOIT.contconictest(MOIB.RootDet{Float64}(MOIB.LogDet{Float64}(optimizer)), config, ["rsoc", "geomean", "psds", "rootdet", "logdets"])
+    end
 end

test/MPBWrapper.jl

@@ -3,10 +3,12 @@
 end
 
 using Compat.Pkg: dir
-@testset "MathProgBase" begin
-    include(joinpath(dir("MathProgBase"),"test","conicinterface.jl"))
-    coniclineartest(SCS.SCSSolver(), duals=true, tol=1e-2)
-    conicSOCtest(SCS.SCSSolver(), duals=true, tol=1e-2)
-    conicEXPtest(SCS.SCSSolver(), duals=true, tol=1e-2)
-    conicSDPtest(SCS.SCSSolver(), duals=true, tol=1e-2)
+for T in [SCS.Direct, SCS.Indirect]
+    @testset "MathProgBase $T" begin
+        include(joinpath(dir("MathProgBase"),"test","conicinterface.jl"))
+        coniclineartest(SCS.SCSSolver(linear_solver=T, eps=1e-7), duals=true, tol=1e-5)
+        conicSOCtest(SCS.SCSSolver(linear_solver=T), duals=true, tol=1e-5)
+        conicEXPtest(SCS.SCSSolver(linear_solver=T), duals=true, tol=1e-5)
+        conicSDPtest(SCS.SCSSolver(linear_solver=T, eps=1e-7), duals=true, tol=1e-5)
+    end
 end

test/options.jl

@@ -20,20 +20,22 @@ s = SCSSolver()
 m = MathProgBase.ConicModel(s)
 MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
 MathProgBase.optimize!(m)
-@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-3)
 
-# With eps = 1e-8, solution should be far more accurate
-s = SCSSolver(eps=1e-8)
+@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-6, rtol=0.0)
+@test !isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-7, rtol=0.0)
+
+# With eps = 1e-12, solution should be far more accurate
+s = SCSSolver(eps=1e-12)
 m = MathProgBase.ConicModel(s)
 MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
 MathProgBase.optimize!(m)
-@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-5)
+@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-10, rtol=0.0)
 
 # With a warmstart from the eps = 1e-8 solution, solution should be extremely accurate even after 1 iteration
 SCS.addoption!(m, :warm_start, true)
 SCS.addoption!(m, :max_iters, 1)
 MathProgBase.optimize!(m)
-@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-5)
+@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-10, rtol=0.0)
 
 # Now let's do the same warmstart, but on a new instance of the same problem
 primal_sol = m.primal_sol
@@ -44,4 +46,34 @@ m = MathProgBase.ConicModel(s)
 MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
 MathProgBase.setwarmstart!(m, primal_sol; dual_sol = dual_sol, slack = slack)
 MathProgBase.optimize!(m)
-@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-5)
+@test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-10, rtol=0.0)
+
+# tests for incorrect options
+s = SCSSolver(eps=1e-12, epps=1.0)
+m = MathProgBase.ConicModel(s)
+MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
+
+@test_throws ArgumentError MathProgBase.optimize!(m)
+
+err = try
+    MathProgBase.optimize!(m)
+catch ex
+    ex
+end
+@test err.msg == "Unrecognized option passed to SCS: epps;
+Recognized options are: normalize, scale, rho_x, max_iters, eps, alpha, cg_rate, verbose, warm_start and acceleration_lookback."
+
+# tests for incorrect options
+s = SCSSolver(linear_solver="AAA", eps=1e-12, epps=1.0)
+m = MathProgBase.ConicModel(s)
+MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
+
+@test_throws ArgumentError MathProgBase.optimize!(m)
+
+err = try
+    MathProgBase.optimize!(m)
+catch ex
+    ex
+end
+
+@test err.msg == "Unrecognized linear_solver passed to SCS: AAA;\nRecognized options are: SCS.Direct, SCS.Indirect."