From d82d637ce8916674974298e1eb19a6bd4f39647e Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sat, 11 Dec 2021 14:39:28 -0500 Subject: [PATCH 1/9] add modify functions --- src/Solvers/Solvers.jl | 229 +++++++++++++++++++--------- src/Solvers/systemsolvers/common.jl | 15 +- test/nativeinstances.jl | 78 +++++++++- test/nativesets.jl | 5 + test/runnativetests.jl | 185 ++++++++++++---------- 5 files changed, 347 insertions(+), 165 deletions(-) diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 7966512a2..ea5485c0c 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -34,6 +34,7 @@ include("point.jl") @enum Status begin NotLoaded Loaded + Modified SolveCalled PrimalInconsistent DualInconsistent @@ -254,13 +255,60 @@ mutable struct Solver{T <: Real} end end +function load(solver::Solver{T}, model::Models.Model{T}) where {T <: Real} + solver.orig_model = model + solver.status = Loaded + return solver +end + + + + +# TODO in the MOI wrapper the b, c, h changes can be cached and called at the +# start of optimize if there is a modify flag. this matters most for h because +# it will be updated constraint by constraint +function modify_c(solver::Solver{T}, c_new::Vector{T}) where {T <: Real} + @assert length(c_new) == solver.orig_model.n + solver.orig_model.c = c_new + solver.status = Modified + return solver +end + +function modify_b(solver::Solver{T}, b_new::Vector{T}) where {T <: Real} + @assert length(b_new) == solver.orig_model.p + solver.orig_model.b = b_new + solver.status = Modified + return solver +end + +function modify_h(solver::Solver{T}, h_new::Vector{T}) where {T <: Real} + @assert length(h_new) == solver.orig_model.q + solver.orig_model.h = h_new + solver.status = Modified + return solver +end + + + + function solve(solver::Solver{T}) where {T <: Real} - (solver.status == Loaded) || return solver - setup_solver(solver) - setup_point(solver) + @show solver.status + if solver.status == Loaded + setup_solver(solver) + setup_loaded(solver) + elseif solver.status == Modified + # TODO easier to update these system solvers (just need an update to setup_point_sub): + @assert solver.syssolver isa Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}} + setup_solver(solver) + setup_modified(solver) + else + # return solver # TODO old + error("solve was called but solver was not just loaded or modified") + end - if solver.status == SolveCalled + if !in(solver.status, (PrimalInconsistent, DualInconsistent)) setup_stepping(solver) + solver.verbose && print_header(solver.stepper, solver) while true step_and_check(solver) && break @@ -284,41 +332,10 @@ function solve(solver::Solver{T}) where {T <: Real} return solver end -function setup_solver(solver::Solver{T}) where {T <: Real} - solver.status = SolveCalled - solver.num_iters = 0 - solver.solve_time = time() - - solver.time_rescale = 0 - solver.time_initx = 0 - solver.time_inity = 0 - solver.time_loadsys = 0 - solver.time_upsys = 0 - solver.time_upfact = 0 - solver.time_uprhs = 0 - solver.time_getdir = 0 - solver.time_search = 0 - - solver.res_norm_cutoff = 0 - solver.max_ref_steps = 5 - solver.x_norm_res_t = NaN - solver.y_norm_res_t = NaN - solver.z_norm_res_t = NaN - solver.x_norm_res = NaN - solver.y_norm_res = NaN - solver.z_norm_res = NaN - solver.primal_obj_t = NaN - solver.dual_obj_t = NaN - solver.primal_obj = NaN - solver.dual_obj = NaN - solver.gap = NaN - solver.x_feas = NaN - solver.y_feas = NaN - solver.z_feas = NaN - solver.tau_feas = NaN +function setup_loaded(solver::Solver{T}) where {T <: Real} orig_model = solver.orig_model if solver.use_dense_model Models.densify!(orig_model) @@ -326,16 +343,53 @@ function setup_solver(solver::Solver{T}) where {T <: Real} solver.result = Point(orig_model) # copy original model to solver.model, which may be modified - solver.model = Models.Model{T}( + model = solver.model = Models.Model{T}( orig_model.c, orig_model.A, orig_model.b, orig_model.G, orig_model.h, orig_model.cones, obj_offset = orig_model.obj_offset) + + solver.time_rescale = @elapsed solver.used_rescaling = rescale_data(solver) + + setup_point(solver) + if solver.status != SolveCalled + return + end + + solver.x_residual = zero(model.c) + solver.y_residual = zero(model.b) + solver.z_residual = zero(model.h) + solver.tau_residual = 0 + + load(solver.stepper, solver) + solver.time_loadsys = @elapsed load(solver.syssolver, solver) + return +end + +function setup_modified(solver::Solver{T}) where {T <: Real} + orig_model = solver.orig_model + model = solver.model + + # TODO needs fixing for rescale/preproc/reduce + @. model.c = orig_model.c + @. model.b = orig_model.b + @. model.h = orig_model.h + + # TODO don't redo factorize stuff etc in here: + setup_point(solver) + if solver.status != SolveCalled + return + end + + # TODO cleaner interface: call a update_syssolver function that calls this: + # (and errors if syssolver does not support update) + set_point_sub_rhs(solver.syssolver, solver.model) return end function setup_point(solver::Solver{T}) where {T <: Real} # preprocess and find initial point (init_z, init_s) = initialize_cone_point(solver.orig_model) - solver.time_rescale = @elapsed solver.used_rescaling = rescale_data(solver) + + if solver.reduce # TODO don't find point / unnecessary stuff before reduce @@ -345,48 +399,82 @@ function setup_point(solver::Solver{T}) where {T <: Real} solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z, false) end + if solver.status != SolveCalled # TODO due to inconsistent statuses + return + end - if solver.status == SolveCalled - model = solver.model - point = solver.point = Point(model) - point.x .= init_x - point.y .= init_y - point.z .= init_z - point.s .= init_s - point.tau[] = one(T) - point.kap[] = one(T) - - calc_mu(solver) - if isnan(solver.mu) || abs(one(T) - solver.mu) > sqrt(eps(T)) - @warn("initial mu is $(solver.mu) but should be 1 (this could " * - "indicate a problem with cone barrier oracles)") - end - Cones.load_point.(model.cones, point.primal_views) - Cones.load_dual_point.(model.cones, point.dual_views) + + + + model = solver.model + point = solver.point = Point(model) + point.x .= init_x + point.y .= init_y + point.z .= init_z + point.s .= init_s + point.tau[] = one(T) + point.kap[] = one(T) + + calc_mu(solver) + if isnan(solver.mu) || abs(one(T) - solver.mu) > sqrt(eps(T)) + @warn("initial mu is $(solver.mu) but should be 1 (this could " * + "indicate a problem with cone barrier oracles)") end + Cones.load_point.(model.cones, point.primal_views) + Cones.load_dual_point.(model.cones, point.dual_views) return end + + + + + function setup_stepping(solver::Solver{T}) where {T <: Real} model = solver.model - # setup iteration helpers - solver.x_residual = zero(model.c) - solver.y_residual = zero(model.b) - solver.z_residual = zero(model.h) - solver.tau_residual = 0 - solver.x_conv_tol = inv(1 + norm(model.c, Inf)) solver.y_conv_tol = inv(1 + norm(model.b, Inf)) solver.z_conv_tol = inv(1 + norm(model.h, Inf)) solver.prev_is_slow = false solver.prev2_is_slow = false solver.worst_dir_res = 0 + return +end - load(solver.stepper, solver) - solver.time_loadsys = @elapsed load(solver.syssolver, solver) +function setup_solver(solver::Solver{T}) where {T <: Real} + solver.status = SolveCalled + solver.num_iters = 0 + solver.solve_time = time() - solver.verbose && print_header(solver.stepper, solver) - flush(stdout) + solver.time_rescale = 0 + solver.time_initx = 0 + solver.time_inity = 0 + solver.time_loadsys = 0 + solver.time_upsys = 0 + solver.time_upfact = 0 + solver.time_uprhs = 0 + solver.time_getdir = 0 + solver.time_search = 0 + + solver.res_norm_cutoff = 0 + solver.max_ref_steps = 5 + + solver.x_norm_res_t = NaN + solver.y_norm_res_t = NaN + solver.z_norm_res_t = NaN + solver.x_norm_res = NaN + solver.y_norm_res = NaN + solver.z_norm_res = NaN + + solver.primal_obj_t = NaN + solver.dual_obj_t = NaN + solver.primal_obj = NaN + solver.dual_obj = NaN + solver.gap = NaN + solver.x_feas = NaN + solver.y_feas = NaN + solver.z_feas = NaN + solver.tau_feas = NaN return end @@ -396,7 +484,6 @@ function step_and_check(solver::Solver{T}) where {T <: Real} if solver.verbose print_iteration(stepper, solver) - flush(stdout) end check_converged(solver) && return true @@ -594,6 +681,8 @@ function initialize_cone_point(model::Models.Model{T}) where {T <: Real} return (init_z, init_s) end +get_model(solver::Solver) = solver.orig_model + get_status(solver::Solver) = solver.status get_solve_time(solver::Solver) = solver.solve_time get_num_iters(solver::Solver) = solver.num_iters @@ -610,12 +699,6 @@ get_tau(solver::Solver) = solver.point.tau[] get_kappa(solver::Solver) = solver.point.kap[] get_mu(solver::Solver) = solver.mu -function load(solver::Solver{T}, model::Models.Model{T}) where {T <: Real} - solver.orig_model = model - solver.status = Loaded - return solver -end - include("process.jl") include("search.jl") @@ -642,6 +725,7 @@ function print_header(stepper::Stepper, solver::Solver) print_header_more(stepper, solver) println() + flush(stdout) return end print_header_more(stepper::Stepper, solver::Solver) = nothing @@ -662,6 +746,7 @@ function print_iteration(stepper::Stepper, solver::Solver) print_iteration_more(stepper, solver) end println() + flush(stdout) return end print_iteration_more(stepper::Stepper, solver::Solver) = nothing diff --git a/src/Solvers/systemsolvers/common.jl b/src/Solvers/systemsolvers/common.jl index 7aac2cd62..2f8c39e4e 100644 --- a/src/Solvers/systemsolvers/common.jl +++ b/src/Solvers/systemsolvers/common.jl @@ -186,24 +186,29 @@ function setup_point_sub( model::Models.Model{T}, ) where {T <: Real} (n, p, q) = (model.n, model.p, model.q) - dim_sub = n + p + q - z_start = model.n + model.p - sol_sub = syssolver.sol_sub = Point{T}() rhs_sub = syssolver.rhs_sub = Point{T}() rhs_const = syssolver.rhs_const = Point{T}() sol_const = syssolver.sol_const = Point{T}() for point_sub in (sol_sub, rhs_sub, rhs_const, sol_const) - point_sub.vec = zeros(T, dim_sub) + point_sub.vec = zeros(T, n + p + q) @views point_sub.x = point_sub.vec[1:n] @views point_sub.y = point_sub.vec[n .+ (1:p)] @views point_sub.z = point_sub.vec[n + p .+ (1:q)] point_sub.z_views = [view(point_sub.z, idxs) for idxs in model.cone_idxs] end + set_point_sub_rhs(syssolver, model) + return +end + +function set_point_sub_rhs( + syssolver::Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}}, + model::Models.Model{T}, + ) where {T <: Real} + rhs_const = syssolver.rhs_const @. rhs_const.x = -model.c @. rhs_const.y = model.b @. rhs_const.z = model.h - return end diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl index 2e9e01060..638451569 100644 --- a/test/nativeinstances.jl +++ b/test/nativeinstances.jl @@ -37,14 +37,17 @@ function build_solve_check( G, h::Vector{T}, cones::Vector{Cone{T}}, - tol::Real, - ; + tol::Real; obj_offset::T = zero(T), solver::Solvers.Solver{T} = Solvers.Solver{T}(), + already_loaded::Bool = false, ) where {T <: Real} - model = Hypatia.Models.Model{T}(c, A, b, G, h, cones, obj_offset = obj_offset) - - Solvers.load(solver, model) + if already_loaded + model = Solvers.get_model(solver) + else + model = Hypatia.Models.Model{T}(c, A, b, G, h, cones, obj_offset = obj_offset) + Solvers.load(solver, model) + end Solvers.solve(solver) status = Solvers.get_status(solver) @@ -2754,3 +2757,68 @@ function indirect5(T; options...) r = build_solve_check(c, A, b, G, h, cones, tol; options...) @test r.status == Solvers.PrimalInfeasible end + + + + +function modify1(T; options...) + tol = test_tol(T) + c = T[1, 0] + A = T[1 1] + b = [-T(2)] + G = SparseMatrixCSC(-one(T) * I, 2, 2) + h = zeros(T, 2) + cones = Cone{T}[Cones.Nonnegative{T}(2)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.PrimalInfeasible + solver = r.solver + +# TODO somehow test that it didn't redo factorization etc +# eg by checking that those subtimings are zero +# but first get it working with symindef without any factorizations + + b = [-T(1)] + Solvers.modify_b(solver, b) + @test Solvers.get_status(solver) == Solvers.Modified + r = build_solve_check(c, A, b, G, h, cones, tol; + solver = solver, already_loaded = true, options...) + @test r.status == Solvers.PrimalInfeasible + + c = T[2, 1] + b = [T(2)] + Solvers.modify_c(solver, c) + Solvers.modify_b(solver, b) + @test Solvers.get_status(solver) == Solvers.Modified + r = build_solve_check(c, A, b, G, h, cones, tol; + solver = solver, already_loaded = true, options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 2 atol=tol rtol=tol + @test r.x ≈ [0, 2] atol=tol rtol=tol + @test iszero(solver.time_rescale) + @test iszero(solver.time_loadsys) + + c = T[3, 1] + h = T[1, -2] + Solvers.modify_c(solver, c) + Solvers.modify_h(solver, h) + @test Solvers.get_status(solver) == Solvers.Modified + r = build_solve_check(c, A, b, G, h, cones, tol; + solver = solver, already_loaded = true, options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test r.x ≈ [-1, 3] atol=tol rtol=tol +end + +function modify2(T; options...) + tol = test_tol(T) + c = T[-1, 0] + A = zeros(T, 0, 2) + b = T[] + G = T[-1 0; 0 0; 0 -1] + h = T[0, 1, 0] + cones = Cone{T}[Cones.EpiPerSquare{T}(3)] +error("") + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.DualInfeasible +end diff --git a/test/nativesets.jl b/test/nativesets.jl index 3ecd15131..b9c4e46f8 100644 --- a/test/nativesets.jl +++ b/test/nativesets.jl @@ -177,3 +177,8 @@ inst_indirect = [ "indirect4", "indirect5", ] + +inst_modify = [ + "modify1", + # "modify2", + ] diff --git a/test/runnativetests.jl b/test/runnativetests.jl index 2eb2853ce..097375d3a 100644 --- a/test/runnativetests.jl +++ b/test/runnativetests.jl @@ -11,8 +11,8 @@ include(joinpath(@__DIR__, "nativesets.jl")) # default options to solvers default_options = ( - # verbose = true, - verbose = false, + verbose = true, + # verbose = false, default_tol_relax = 10, ) @@ -23,7 +23,7 @@ all_reals = [ ] diff_reals = [ Float64, - BigFloat, + # BigFloat, ] string_nameof(T) = string(nameof(T)) @@ -46,96 +46,115 @@ end @testset "native tests" begin -@testset "default options tests" begin - println("starting default options tests") - inst_defaults = vcat( - inst_preproc, - inst_infeas, - inst_cones_many, - ) - for inst_name in inst_defaults - test_instance_solver(inst_name, Float64, default_options) - end -end +# @testset "default options tests" begin +# println("starting default options tests") +# inst_defaults = vcat( +# inst_preproc, +# inst_infeas, +# inst_cones_many, +# ) +# for inst_name in inst_defaults +# test_instance_solver(inst_name, Float64, default_options) +# end +# end -@testset "no preprocess tests" begin - println("\nstarting no preprocess tests") - for inst_name in inst_cones_few, T in diff_reals - options = (; default_options..., preprocess = false, reduce = false, - syssolver = Solvers.SymIndefDenseSystemSolver{T}()) - test_instance_solver(inst_name, T, options) - end -end +# @testset "no preprocess tests" begin +# println("\nstarting no preprocess tests") +# for inst_name in inst_cones_few, T in diff_reals +# options = (; default_options..., preprocess = false, reduce = false, +# syssolver = Solvers.SymIndefDenseSystemSolver{T}()) +# test_instance_solver(inst_name, T, options) +# end +# end -@testset "indirect solvers tests" begin - println("\nstarting indirect solvers tests") - for inst_name in inst_indirect, T in diff_reals - options = (; default_options..., init_use_indirect = true, - preprocess = false, reduce = false, - syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), - tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, - tol_infeas = 1e-6) - test_instance_solver(inst_name, T, options) - end -end +# @testset "indirect solvers tests" begin +# println("\nstarting indirect solvers tests") +# for inst_name in inst_indirect, T in diff_reals +# options = (; default_options..., init_use_indirect = true, +# preprocess = false, reduce = false, +# syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), +# tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, +# tol_infeas = 1e-6) +# test_instance_solver(inst_name, T, options) +# end +# end + +# @testset "system solvers tests" begin +# println("\nstarting system solvers tests") +# syssolvers = [ +# (Solvers.NaiveDenseSystemSolver, diff_reals), +# (Solvers.NaiveSparseSystemSolver, [Float64,]), +# (Solvers.NaiveElimDenseSystemSolver, diff_reals), +# (Solvers.NaiveElimSparseSystemSolver, [Float64,]), +# (Solvers.SymIndefDenseSystemSolver, all_reals), +# (Solvers.SymIndefSparseSystemSolver, [Float64,]), +# (Solvers.QRCholDenseSystemSolver, all_reals), +# ] +# for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, +# T in real_types +# options = (; default_options..., syssolver = syssolver{T}(), +# reduce = false) +# test_instance_solver(inst_name, T, options, string_nameof(syssolver)) +# end +# end + +# @testset "PredOrCentStepper tests" begin +# verbose = true +# println("\nstarting PredOrCentStepper tests (with printing)") + +# # adjustment and curve search +# use_adj_curv = [(false, false), (true, false), (true, true)] +# for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals +# stepper = Solvers.PredOrCentStepper{T}(; +# use_adjustment = adj, use_curve_search = curv) +# options = (; default_options..., verbose = verbose, stepper = stepper) +# test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") +# end -@testset "system solvers tests" begin - println("\nstarting system solvers tests") +# # other options +# for inst_name in inst_minimal +# T = Float64 +# stepper = Solvers.PredOrCentStepper{T}(; +# # stepper options +# use_adjustment = false, use_curve_search = false, +# max_cent_steps = 8, pred_prox_bound = 0.0332, +# # searcher options +# min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, +# alpha_sched = [0.9999 * 0.7^i for i in 0:22]) +# options = (; default_options..., verbose = verbose, stepper = stepper) +# test_instance_solver(inst_name, T, options, "other") +# end +# end + +# @testset "CombinedStepper tests" begin +# verbose = true +# println("\nstarting CombinedStepper tests (with printing)") +# shifts = [0, 2] +# for inst_name in inst_minimal, shift in shifts, T in diff_reals +# options = (; default_options..., verbose = verbose, +# stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) +# test_instance_solver(inst_name, T, options, "shift=$shift") +# end +# end + +@testset "model modification tests" begin + println("\nstarting model modification tests") + # for inst_name in inst_modify, T in diff_reals + # test_instance_solver(inst_name, T, default_options) + # end syssolvers = [ - (Solvers.NaiveDenseSystemSolver, diff_reals), - (Solvers.NaiveSparseSystemSolver, [Float64,]), - (Solvers.NaiveElimDenseSystemSolver, diff_reals), - (Solvers.NaiveElimSparseSystemSolver, [Float64,]), - (Solvers.SymIndefDenseSystemSolver, all_reals), - (Solvers.SymIndefSparseSystemSolver, [Float64,]), - (Solvers.QRCholDenseSystemSolver, all_reals), + (Solvers.SymIndefDenseSystemSolver, diff_reals), + # (Solvers.SymIndefSparseSystemSolver, [Float64,]), + # (Solvers.QRCholDenseSystemSolver, diff_reals), ] - for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, + for inst_name in inst_modify, (syssolver, real_types) in syssolvers, T in real_types + # TODO later also rescale/preproc/reduce true options = (; default_options..., syssolver = syssolver{T}(), - reduce = false) + rescale = false, preprocess = false, reduce = false) test_instance_solver(inst_name, T, options, string_nameof(syssolver)) end end -@testset "PredOrCentStepper tests" begin - verbose = true - println("\nstarting PredOrCentStepper tests (with printing)") - - # adjustment and curve search - use_adj_curv = [(false, false), (true, false), (true, true)] - for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals - stepper = Solvers.PredOrCentStepper{T}(; - use_adjustment = adj, use_curve_search = curv) - options = (; default_options..., verbose = verbose, stepper = stepper) - test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") - end - - # other options - for inst_name in inst_minimal - T = Float64 - stepper = Solvers.PredOrCentStepper{T}(; - # stepper options - use_adjustment = false, use_curve_search = false, - max_cent_steps = 8, pred_prox_bound = 0.0332, - # searcher options - min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, - alpha_sched = [0.9999 * 0.7^i for i in 0:22]) - options = (; default_options..., verbose = verbose, stepper = stepper) - test_instance_solver(inst_name, T, options, "other") - end -end - -@testset "CombinedStepper tests" begin - verbose = true - println("\nstarting CombinedStepper tests (with printing)") - shifts = [0, 2] - for inst_name in inst_minimal, shift in shifts, T in diff_reals - options = (; default_options..., verbose = verbose, - stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) - test_instance_solver(inst_name, T, options, "shift=$shift") - end -end - end ; From 2ad121518f4e443520e27f7715e2cc34756cf88a Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sat, 11 Dec 2021 17:46:03 -0500 Subject: [PATCH 2/9] start cleanup of processing --- src/Solvers/Solvers.jl | 43 +++--- src/Solvers/process.jl | 333 ++++++++++++++++++++++------------------ test/nativeinstances.jl | 36 ++--- test/runnativetests.jl | 196 ++++++++++++----------- test/runtests.jl | 4 +- 5 files changed, 328 insertions(+), 284 deletions(-) diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index ea5485c0c..77e3776c3 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -198,7 +198,7 @@ mutable struct Solver{T <: Real} if reduce @assert preprocess # cannot use reduction without preprocessing # TODO only need primal eq preprocessing end - # @assert !(init_use_indirect && preprocess) # cannot use preprocessing and indirect methods for initial point + @assert !(init_use_indirect && preprocess) # cannot use preprocessing and indirect methods for initial point @assert near_factor >= 1 # factor to relax tolerances by if fail to converge should be at least 1 if isnothing(default_tol_power) @@ -292,17 +292,13 @@ end function solve(solver::Solver{T}) where {T <: Real} - @show solver.status - if solver.status == Loaded - setup_solver(solver) + init_status = solver.status + setup_solver(solver) + if init_status == Loaded setup_loaded(solver) - elseif solver.status == Modified - # TODO easier to update these system solvers (just need an update to setup_point_sub): - @assert solver.syssolver isa Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}} - setup_solver(solver) + elseif init_status == Modified setup_modified(solver) else - # return solver # TODO old error("solve was called but solver was not just loaded or modified") end @@ -354,13 +350,13 @@ function setup_loaded(solver::Solver{T}) where {T <: Real} return end + load(solver.stepper, solver) + solver.time_loadsys = @elapsed load(solver.syssolver, solver) + solver.x_residual = zero(model.c) solver.y_residual = zero(model.b) solver.z_residual = zero(model.h) solver.tau_residual = 0 - - load(solver.stepper, solver) - solver.time_loadsys = @elapsed load(solver.syssolver, solver) return end @@ -368,10 +364,19 @@ function setup_modified(solver::Solver{T}) where {T <: Real} orig_model = solver.orig_model model = solver.model - # TODO needs fixing for rescale/preproc/reduce - @. model.c = orig_model.c - @. model.b = orig_model.b - @. model.h = orig_model.h + copyto!(model.c, orig_model.c) + copyto!(model.b, orig_model.b) + copyto!(model.h, orig_model.h) + if solver.used_rescaling + # rescale using scalings computed during first load + model.c ./= solver.c_scale + model.b ./= solver.b_scale + model.h ./= solver.h_scale + end + # TODO needs fixing for preproc/reduce + + + # TODO don't redo factorize stuff etc in here: setup_point(solver) @@ -379,6 +384,8 @@ function setup_modified(solver::Solver{T}) where {T <: Real} return end + # TODO easier to update these system solvers (just need an update to setup_point_sub): + @assert solver.syssolver isa Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}} # TODO cleaner interface: call a update_syssolver function that calls this: # (and errors if syssolver does not support update) set_point_sub_rhs(solver.syssolver, solver.model) @@ -393,11 +400,11 @@ function setup_point(solver::Solver{T}) where {T <: Real} if solver.reduce # TODO don't find point / unnecessary stuff before reduce - solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z, true) + solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z) solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) else solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) - solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z, false) + solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z) end if solver.status != SolveCalled # TODO due to inconsistent statuses return diff --git a/src/Solvers/process.jl b/src/Solvers/process.jl index 8be11cd97..cc66c18f3 100644 --- a/src/Solvers/process.jl +++ b/src/Solvers/process.jl @@ -14,7 +14,9 @@ function rescale_data(solver::Solver{T}) where {T <: Real} solver.rescale || return false model = solver.model (c, A, b, G, h) = (model.c, model.A, model.b, model.G, model.h) - (isa(A, MatrixyAG) && isa(G, MatrixyAG)) || return false + if !isa(A, MatrixyAG) || !isa(G, MatrixyAG) + return false + end minval = sqrt(eps(T)) maxabsmin(v::AbstractVecOrMat) = mapreduce(abs, max, v; init = minval) @@ -182,38 +184,32 @@ end function find_initial_y( solver::Solver{T}, init_z::Vector{T}, - reduce::Bool, ) where {T <: Real} if solver.status != SolveCalled return zeros(T, 0) end model = solver.model p = model.p - if iszero(p) # y is empty (no primal variables) + + if iszero(p) + # y is empty (no primal variables) solver.y_keep_idxs = Int[] solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) solver.Ap_Q = I return zeros(T, 0) end - n = model.n - q = model.q - A = model.A - solver.y_keep_idxs = 1:p - if !reduce || !isa(A, MatrixyAG) - # rhs = -c - G' * point.z - rhs = copy(model.c) - mul!(rhs, model.G', init_z, -1, -1) + solver.y_keep_idxs = 1:p + rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1) - # indirect method - if solver.init_use_indirect - # TODO pick lsqr or lsmr - init_y = IterativeSolvers.lsqr(A', rhs) - return init_y - end + if solver.init_use_indirect + # use indirect method + # TODO pick lsqr or lsmr + return IterativeSolvers.lsqr(model.A', rhs_y) end # factorize A' + A = model.A if issparse(A) if !(T <: Float64) @warn("using dense factorization of A' in preprocessing and initial " * @@ -226,28 +222,39 @@ function find_initial_y( else Ap_fact = qr!(Matrix(A'), ColumnNorm()) end - Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) - if !reduce && !solver.preprocess - if Ap_rank < p - @warn("some primal equalities appear to be dependent " * - "(possibly inconsistent); try using preprocess = true") - end - init_y = Ap_fact \ rhs - return init_y + if solver.preprocess + return process_primal_eq(solver, Ap_fact, rhs_y) end - # preprocess dual equalities + Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) + if Ap_rank < p + @warn("some primal equalities appear to be dependent " * + "(possibly inconsistent); try using preprocess = true") + end + + return Ap_fact \ rhs_y +end + +# preprocess primal equalities +function process_primal_eq( + solver::Solver{T}, + Ap_fact::Factorization{T}, + rhs_y::Vector{T} + ) where {T <: Real} + Ap_Q = Ap_fact.Q + Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) Ap_R = UpperTriangular(Ap_fact.R[1:Ap_rank, 1:Ap_rank]) col_piv = (Ap_fact isa QRPivoted{T, Matrix{T}}) ? Ap_fact.p : Ap_fact.pcol - y_keep_idxs = col_piv[1:Ap_rank] - Ap_Q = Ap_fact.Q + model = solver.model + y_keep_idxs = col_piv[1:Ap_rank] b_sub = model.b[y_keep_idxs] - if Ap_rank < p + + if Ap_rank < model.p # some dependent primal equalities, so check if they are consistent # x_sub = Ap_Q * vcat((Ap_R' \ b_sub), zeros(n - Ap_rank)) - x_sub = zeros(T, n) + x_sub = zeros(T, model.n) x_sub1 = view(x_sub, 1:Ap_rank) copyto!(x_sub1, b_sub) ldiv!(Ap_R', x_sub1) @@ -256,7 +263,7 @@ function find_initial_y( if !(Ap_fact isa QRPivoted{T, Matrix{T}}) x_sub = x_sub[Ap_fact.rpivinv] end - residual = norm(A * x_sub - model.b, Inf) + residual = norm(model.A * x_sub - model.b, Inf) if residual > solver.init_tol_qr if solver.verbose println("some primal equality constraints are inconsistent " * @@ -266,94 +273,27 @@ function find_initial_y( return zeros(T, 0) end if solver.verbose + p = model.p println("$(p - Ap_rank) of $p primal equality constraints " * "are dependent") end end - if reduce && isa(model.G, MatrixyAG) - # remove all primal equalities by making A and b empty with - # n = n0 - p0 and p = 0 - # TODO improve efficiency - # TODO avoid calculating GQ1 explicitly if possible - # recover original-space solution using: - # x0 = Q * [(R' \ b0), x] - # y0 = R \ (-cQ1' - GQ1' * z0) - if !(Ap_fact isa QRPivoted{T, Matrix{T}}) - row_piv = Ap_fact.prow - model.c = model.c[row_piv] - model.G = model.G[:, row_piv] - solver.reduce_row_piv_inv = Ap_fact.rpivinv - else - solver.reduce_row_piv_inv = Int[] - end - - Q1_idxs = 1:Ap_rank - Q2_idxs = (Ap_rank + 1):n - - # [cQ1 cQ2] = c0' * Q - cQ = model.c' * Ap_Q - cQ1 = solver.reduce_cQ1 = cQ[Q1_idxs] - cQ2 = cQ[Q2_idxs] - # c = cQ2 - model.c = cQ2 - model.n = length(model.c) - # offset = offset0 + cQ1 * (R' \ b0) - Rpib0 = solver.reduce_Rpib0 = ldiv!(Ap_R', b_sub) - # solver.Rpib0 = Rpib0 # TODO - model.obj_offset += dot(cQ1, Rpib0) - - # [GQ1 GQ2] = G0 * Q - # very inefficient method used for sparse G * QRSparseQ - # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818 - if model.G isa UniformScaling - side = size(Ap_Q, 1) - G_mul = Matrix{T}(model.G, side, side) - elseif !isa(model.G, Matrix) - G_mul = Matrix{T}(model.G) - else - G_mul = model.G - end - GQ = G_mul * Ap_Q - GQ1 = solver.reduce_GQ1 = GQ[:, Q1_idxs] - GQ2 = GQ[:, Q2_idxs] - # h = h0 - GQ1 * (R' \ b0) - mul!(model.h, GQ1, Rpib0, -1, true) - - # G = GQ2 - model.G = GQ2 - - # A and b empty - model.p = 0 - model.A = zeros(T, 0, model.n) - model.b = zeros(T, 0) - solver.reduce_Ap_R = Ap_R - solver.reduce_Ap_Q = Ap_Q - - solver.reduce_y_keep_idxs = y_keep_idxs - solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) - solver.Ap_Q = I - return zeros(T, 0) +# TODO refac some of reduce_primal_eq and below + if solver.reduce && isa(model.G, MatrixyAG) + return reduce_primal_eq(solver, Ap_fact, Ap_rank, Ap_R, y_keep_idxs, b_sub) end - # init_y = Ap_R \ ((Ap_fact.Q' * (-c - G' * point.z))[1:Ap_rank]) - temp = copy(model.c) - mul!(temp, model.G', init_z, true, true) - lmul!(Ap_fact.Q', temp) - init_y = temp[1:Ap_rank] - init_y .*= -1 - ldiv!(Ap_R, init_y) - # modify solver.model to remove/reorder some dual variables y if !(Ap_fact isa QRPivoted{T, Matrix{T}}) row_piv = Ap_fact.prow - model.A = A[y_keep_idxs, row_piv] + model.A = model.A[y_keep_idxs, row_piv] model.c = model.c[row_piv] model.G = model.G[:, row_piv] solver.x_keep_idxs = solver.x_keep_idxs[row_piv] else - model.A = A[y_keep_idxs, :] + model.A = model.A[y_keep_idxs, :] end model.b = b_sub model.p = Ap_rank @@ -361,9 +301,90 @@ function find_initial_y( solver.Ap_R = Ap_R solver.Ap_Q = Ap_Q + + + # init_y = Ap_R \ ((Ap_Q' * (-c - G' * z))[1:Ap_rank]) + lmul!(Ap_Q', rhs_y) + init_y = rhs_y[1:Ap_rank] + ldiv!(Ap_R, init_y) return init_y end +# remove all primal equalities i.e. let n = n0 - p0 and p = 0 +# TODO improve efficiency +# TODO avoid calculating GQ1 explicitly if possible +# recover original-space solution using: +# x0 = Q * [(R' \ b0), x] +# y0 = R \ (-cQ1' - GQ1' * z0) +function reduce_primal_eq( + solver::Solver{T}, + Ap_fact::Factorization{T}, + Ap_rank::Int, + Ap_R, + y_keep_idxs, + b_sub, + ) where {T <: Real} + model = solver.model + + if !(Ap_fact isa QRPivoted{T, Matrix{T}}) + row_piv = Ap_fact.prow + model.c = model.c[row_piv] + model.G = model.G[:, row_piv] + solver.reduce_row_piv_inv = Ap_fact.rpivinv + else + solver.reduce_row_piv_inv = Int[] + end + + Q1_idxs = 1:Ap_rank + Q2_idxs = (Ap_rank + 1):model.n + + # [cQ1 cQ2] = c0' * Q + Ap_Q = Ap_fact.Q + cQ = model.c' * Ap_Q + cQ1 = solver.reduce_cQ1 = cQ[Q1_idxs] + cQ2 = cQ[Q2_idxs] + # c = cQ2 + model.c = cQ2 + model.n = length(cQ2) + # offset = offset0 + cQ1 * (R' \ b0) + Rpib0 = solver.reduce_Rpib0 = ldiv!(Ap_R', b_sub) + # solver.Rpib0 = Rpib0 # TODO + model.obj_offset += dot(cQ1, Rpib0) + + # [GQ1 GQ2] = G0 * Q + # very inefficient method used for sparse G * QRSparseQ + # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818 + if model.G isa UniformScaling + side = size(Ap_Q, 1) + G_mul = Matrix{T}(model.G, side, side) + elseif !isa(model.G, Matrix) + G_mul = Matrix{T}(model.G) + else + G_mul = model.G + end + GQ = G_mul * Ap_Q + GQ1 = solver.reduce_GQ1 = GQ[:, Q1_idxs] + GQ2 = GQ[:, Q2_idxs] + # h = h0 - GQ1 * (R' \ b0) + mul!(model.h, GQ1, Rpib0, -1, true) + + # G = GQ2 + model.G = GQ2 + + # A and b empty + model.p = 0 + model.A = zeros(T, 0, model.n) + model.b = zeros(T, 0) + solver.reduce_Ap_R = Ap_R + solver.reduce_Ap_Q = Ap_Q + + solver.reduce_y_keep_idxs = y_keep_idxs + solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) + solver.Ap_Q = I + + return zeros(T, 0) +end + # (pivoted) QR factorizations are usually rank-revealing but may be unreliable # see http://www.math.sjsu.edu/~foster/rankrevealingcode.html # TODO could replace this with rank(qr_fact) when available for both dense and sparse @@ -385,65 +406,46 @@ end function postprocess(solver::Solver{T}) where {T <: Real} point = solver.point result = solver.result - if solver.status in infeas_statuses - tau = one(T) - else - tau = point.tau[] - if tau <= 0 - result.vec .= NaN - return - end + tau = point.tau[] + if tau <= 0 || any(isnan, point.vec) + result.vec .= NaN + return end - # finalize s,z - @. result.s = point.s / tau - @. result.z = point.z / tau + copyto!(result.s, point.s) + copyto!(result.z, point.z) + x = copy(point.x) + y = copy(point.y) + if !in(solver.status, infeas_statuses) + # rescale non-infeasible certificates by 1/tau + result.s ./= tau + result.z ./= tau + x ./= tau + y ./= tau + end # finalize x - if solver.preprocess && !iszero(solver.orig_model.n) && !any(isnan, point.x) + if solver.preprocess && !iszero(solver.orig_model.n) # unpreprocess solver's solution if solver.reduce && !iszero(solver.orig_model.p) - # unreduce solver's solution - # x = Q * [(R' \ b0), x] - xa = zeros(T, solver.orig_model.n - length(solver.reduce_Rpib0)) - @. @views xa[solver.x_keep_idxs] = point.x / tau - if solver.status in infeas_statuses - Rpib0 = zeros(T, length(solver.reduce_Rpib0)) - else - Rpib0 = solver.reduce_Rpib0 - end - xb = vcat(Rpib0, xa) - lmul!(solver.reduce_Ap_Q, xb) - if isempty(solver.reduce_row_piv_inv) - result.x .= xb - else - @. @views result.x = xb[solver.reduce_row_piv_inv] - end + unreduce_x(solver, x) else - @. @views result.x[solver.x_keep_idxs] = point.x / tau + @views copyto!(result.x[solver.x_keep_idxs], x) end else - @. result.x = point.x / tau + copyto!(result.x, x) end # finalize y - if solver.preprocess && !iszero(solver.orig_model.p) && !any(isnan, point.y) + if solver.preprocess && !iszero(solver.orig_model.p) # unpreprocess solver's solution if solver.reduce - # unreduce solver's solution - # y = R \ (-cQ1' - GQ1' * z) - ya = solver.reduce_GQ1' * result.z - if !(solver.status in infeas_statuses) - ya .+= solver.reduce_cQ1 - end - @views ldiv!(solver.reduce_Ap_R, - ya[1:length(solver.reduce_y_keep_idxs)]) - @. @views result.y[solver.reduce_y_keep_idxs] = -ya + unreduce_y(solver, y) else - @. @views result.y[solver.y_keep_idxs] = point.y / tau + @views copyto!(result.y[solver.y_keep_idxs], y) end else - @. result.y = point.y / tau + copyto!(result.y, y) end if solver.used_rescaling @@ -456,3 +458,42 @@ function postprocess(solver::Solver{T}) where {T <: Real} return end + +function unreduce_x( + solver::Solver{T}, + x::Vector{T}, + ) where {T <: Real} + # x = Q * [(R' \ b0), x] + Rpib0 = if solver.status in infeas_statuses + zero(solver.reduce_Rpib0) + else + solver.reduce_Rpib0 + end + xa = zeros(T, solver.orig_model.n - length(Rpib0)) + @views copyto!(xa[solver.x_keep_idxs], x) + xb = vcat(Rpib0, xa) + lmul!(solver.reduce_Ap_Q, xb) + result_x = solver.result.x + if isempty(solver.reduce_row_piv_inv) + copyto!(result_x, xb) + else + @views copyto!(result_x, xb[solver.reduce_row_piv_inv]) + end + return +end + +function unreduce_y( + solver::Solver{T}, + y::Vector{T}, + ) where {T <: Real} + # y = R \ (-cQ1' - GQ1' * z) + ya = solver.reduce_GQ1' * solver.result.z + if !in(solver.status, infeas_statuses) + ya .+= solver.reduce_cQ1 + end + @views ya_sub = ya[1:length(solver.reduce_y_keep_idxs)] + ldiv!(solver.reduce_Ap_R, ya_sub) + @views y_sub = solver.result.y[solver.reduce_y_keep_idxs] + @. y_sub = -ya + return +end diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl index 638451569..f628816c1 100644 --- a/test/nativeinstances.jl +++ b/test/nativeinstances.jl @@ -2758,9 +2758,6 @@ function indirect5(T; options...) @test r.status == Solvers.PrimalInfeasible end - - - function modify1(T; options...) tol = test_tol(T) c = T[1, 0] @@ -2774,16 +2771,14 @@ function modify1(T; options...) @test r.status == Solvers.PrimalInfeasible solver = r.solver -# TODO somehow test that it didn't redo factorization etc -# eg by checking that those subtimings are zero -# but first get it working with symindef without any factorizations - b = [-T(1)] Solvers.modify_b(solver, b) @test Solvers.get_status(solver) == Solvers.Modified r = build_solve_check(c, A, b, G, h, cones, tol; solver = solver, already_loaded = true, options...) @test r.status == Solvers.PrimalInfeasible + @test iszero(solver.time_rescale) + @test iszero(solver.time_loadsys) c = T[2, 1] b = [T(2)] @@ -2808,17 +2803,20 @@ function modify1(T; options...) @test r.status == Solvers.Optimal @test r.primal_obj ≈ 0 atol=tol rtol=tol @test r.x ≈ [-1, 3] atol=tol rtol=tol + @test iszero(solver.time_rescale) + @test iszero(solver.time_loadsys) end -function modify2(T; options...) - tol = test_tol(T) - c = T[-1, 0] - A = zeros(T, 0, 2) - b = T[] - G = T[-1 0; 0 0; 0 -1] - h = T[0, 1, 0] - cones = Cone{T}[Cones.EpiPerSquare{T}(3)] -error("") - r = build_solve_check(c, A, b, G, h, cones, tol; options...) - @test r.status == Solvers.DualInfeasible -end +# function modify2(T; options...) +# tol = test_tol(T) +# c = T[-1, -1, 0] +# A = zeros(T, 0, 3) +# b = T[] +# G = repeat(SparseMatrixCSC(-one(T) * I, 3, 3), 2, 1) +# h = zeros(T, 6) +# cones = Cone{T}[Cones.EpiNormInf{T, T}(3), +# Cones.EpiNormInf{T, T}(3, use_dual = true)] +# error() +# r = build_solve_check(c, A, b, G, h, cones, tol; options...) +# @test r.status == Solvers.DualInfeasible +# end diff --git a/test/runnativetests.jl b/test/runnativetests.jl index 097375d3a..5bbb26c95 100644 --- a/test/runnativetests.jl +++ b/test/runnativetests.jl @@ -11,8 +11,8 @@ include(joinpath(@__DIR__, "nativesets.jl")) # default options to solvers default_options = ( - verbose = true, - # verbose = false, + # verbose = true, + verbose = false, default_tol_relax = 10, ) @@ -23,7 +23,7 @@ all_reals = [ ] diff_reals = [ Float64, - # BigFloat, + BigFloat, ] string_nameof(T) = string(nameof(T)) @@ -46,113 +46,111 @@ end @testset "native tests" begin -# @testset "default options tests" begin -# println("starting default options tests") -# inst_defaults = vcat( -# inst_preproc, -# inst_infeas, -# inst_cones_many, -# ) -# for inst_name in inst_defaults -# test_instance_solver(inst_name, Float64, default_options) -# end -# end - -# @testset "no preprocess tests" begin -# println("\nstarting no preprocess tests") -# for inst_name in inst_cones_few, T in diff_reals -# options = (; default_options..., preprocess = false, reduce = false, -# syssolver = Solvers.SymIndefDenseSystemSolver{T}()) -# test_instance_solver(inst_name, T, options) -# end -# end - -# @testset "indirect solvers tests" begin -# println("\nstarting indirect solvers tests") -# for inst_name in inst_indirect, T in diff_reals -# options = (; default_options..., init_use_indirect = true, -# preprocess = false, reduce = false, -# syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), -# tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, -# tol_infeas = 1e-6) -# test_instance_solver(inst_name, T, options) -# end -# end - -# @testset "system solvers tests" begin -# println("\nstarting system solvers tests") -# syssolvers = [ -# (Solvers.NaiveDenseSystemSolver, diff_reals), -# (Solvers.NaiveSparseSystemSolver, [Float64,]), -# (Solvers.NaiveElimDenseSystemSolver, diff_reals), -# (Solvers.NaiveElimSparseSystemSolver, [Float64,]), -# (Solvers.SymIndefDenseSystemSolver, all_reals), -# (Solvers.SymIndefSparseSystemSolver, [Float64,]), -# (Solvers.QRCholDenseSystemSolver, all_reals), -# ] -# for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, -# T in real_types -# options = (; default_options..., syssolver = syssolver{T}(), -# reduce = false) -# test_instance_solver(inst_name, T, options, string_nameof(syssolver)) -# end -# end - -# @testset "PredOrCentStepper tests" begin -# verbose = true -# println("\nstarting PredOrCentStepper tests (with printing)") - -# # adjustment and curve search -# use_adj_curv = [(false, false), (true, false), (true, true)] -# for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals -# stepper = Solvers.PredOrCentStepper{T}(; -# use_adjustment = adj, use_curve_search = curv) -# options = (; default_options..., verbose = verbose, stepper = stepper) -# test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") -# end - -# # other options -# for inst_name in inst_minimal -# T = Float64 -# stepper = Solvers.PredOrCentStepper{T}(; -# # stepper options -# use_adjustment = false, use_curve_search = false, -# max_cent_steps = 8, pred_prox_bound = 0.0332, -# # searcher options -# min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, -# alpha_sched = [0.9999 * 0.7^i for i in 0:22]) -# options = (; default_options..., verbose = verbose, stepper = stepper) -# test_instance_solver(inst_name, T, options, "other") -# end -# end - -# @testset "CombinedStepper tests" begin -# verbose = true -# println("\nstarting CombinedStepper tests (with printing)") -# shifts = [0, 2] -# for inst_name in inst_minimal, shift in shifts, T in diff_reals -# options = (; default_options..., verbose = verbose, -# stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) -# test_instance_solver(inst_name, T, options, "shift=$shift") -# end -# end +@testset "default options tests" begin + println("starting default options tests") + inst_defaults = vcat( + inst_preproc, + inst_infeas, + # inst_cones_few, + inst_cones_many, + ) + for inst_name in inst_defaults + test_instance_solver(inst_name, Float64, default_options) + end +end + +@testset "no preprocess tests" begin + println("\nstarting no preprocess tests") + for inst_name in inst_cones_few, T in diff_reals + options = (; default_options..., preprocess = false, reduce = false, + syssolver = Solvers.SymIndefDenseSystemSolver{T}()) + test_instance_solver(inst_name, T, options) + end +end + +@testset "indirect solvers tests" begin + println("\nstarting indirect solvers tests") + for inst_name in inst_indirect, T in diff_reals + options = (; default_options..., init_use_indirect = true, + preprocess = false, reduce = false, + syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), + tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, + tol_infeas = 1e-6) + test_instance_solver(inst_name, T, options) + end +end + +@testset "system solvers tests" begin + println("\nstarting system solvers tests") + syssolvers = [ + (Solvers.NaiveDenseSystemSolver, diff_reals), + (Solvers.NaiveSparseSystemSolver, [Float64,]), + (Solvers.NaiveElimDenseSystemSolver, diff_reals), + (Solvers.NaiveElimSparseSystemSolver, [Float64,]), + (Solvers.SymIndefDenseSystemSolver, all_reals), + (Solvers.SymIndefSparseSystemSolver, [Float64,]), + (Solvers.QRCholDenseSystemSolver, all_reals), + ] + for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, + T in real_types + options = (; default_options..., syssolver = syssolver{T}(), + reduce = false) + test_instance_solver(inst_name, T, options, string_nameof(syssolver)) + end +end + +@testset "PredOrCentStepper tests" begin + verbose = true + println("\nstarting PredOrCentStepper tests (with printing)") + + # adjustment and curve search + use_adj_curv = [(false, false), (true, false), (true, true)] + for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals + stepper = Solvers.PredOrCentStepper{T}(; + use_adjustment = adj, use_curve_search = curv) + options = (; default_options..., verbose = verbose, stepper = stepper) + test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") + end + + # other options + for inst_name in inst_minimal + T = Float64 + stepper = Solvers.PredOrCentStepper{T}(; + # stepper options + use_adjustment = false, use_curve_search = false, + max_cent_steps = 8, pred_prox_bound = 0.0332, + # searcher options + min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, + alpha_sched = [0.9999 * 0.7^i for i in 0:22]) + options = (; default_options..., verbose = verbose, stepper = stepper) + test_instance_solver(inst_name, T, options, "other") + end +end + +@testset "CombinedStepper tests" begin + verbose = true + println("\nstarting CombinedStepper tests (with printing)") + shifts = [0, 2] + for inst_name in inst_minimal, shift in shifts, T in diff_reals + options = (; default_options..., verbose = verbose, + stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) + test_instance_solver(inst_name, T, options, "shift=$shift") + end +end @testset "model modification tests" begin println("\nstarting model modification tests") - # for inst_name in inst_modify, T in diff_reals - # test_instance_solver(inst_name, T, default_options) - # end syssolvers = [ (Solvers.SymIndefDenseSystemSolver, diff_reals), # (Solvers.SymIndefSparseSystemSolver, [Float64,]), # (Solvers.QRCholDenseSystemSolver, diff_reals), ] for inst_name in inst_modify, (syssolver, real_types) in syssolvers, - T in real_types + T in real_types, rescale in (false, true) # TODO later also rescale/preproc/reduce true options = (; default_options..., syssolver = syssolver{T}(), - rescale = false, preprocess = false, reduce = false) - test_instance_solver(inst_name, T, options, string_nameof(syssolver)) + rescale = rescale, preprocess = false, reduce = false) + test_instance_solver(inst_name, T, options) end end diff --git a/test/runtests.jl b/test/runtests.jl index 85f995a71..3f830c86f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,8 +7,8 @@ using Printf using Hypatia test_files = [ - "polyutils", - "cone", + # "polyutils", + # "cone", "native", "moi", # "examples", # CI runs examples separately From 342f99837dad219e65230f9e6fbe6542b927a8a1 Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 12 Dec 2021 00:16:48 -0500 Subject: [PATCH 3/9] more preproc cleanup and split out updating part --- src/Solvers/Solvers.jl | 31 +- src/Solvers/process.jl | 449 ++++++++++++++++++++-------- src/Solvers/systemsolvers/qrchol.jl | 6 +- test/runnativetests.jl | 158 +++++----- 4 files changed, 431 insertions(+), 213 deletions(-) diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 77e3776c3..9e282f668 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -111,6 +111,8 @@ mutable struct Solver{T <: Real} model::Models.Model{T} x_keep_idxs::AbstractVector{Int} y_keep_idxs::AbstractVector{Int} + Ap_fact::Factorization{T} + Ap_rank::Int Ap_R::UpperTriangular{T, <:AbstractMatrix{T}} Ap_Q::Union{UniformScaling, AbstractMatrix{T}} reduce_cQ1 @@ -293,13 +295,16 @@ end function solve(solver::Solver{T}) where {T <: Real} init_status = solver.status + if !in(init_status, (Loaded, Modified)) + @warn("solve called when solver status is $init_status") + return + end + setup_solver(solver) - if init_status == Loaded - setup_loaded(solver) - elseif init_status == Modified + if init_status == Modified setup_modified(solver) else - error("solve was called but solver was not just loaded or modified") + setup_loaded(solver) end if !in(solver.status, (PrimalInconsistent, DualInconsistent)) @@ -325,7 +330,7 @@ function solve(solver::Solver{T}) where {T <: Real} free_memory(solver.syssolver) flush(stdout) - return solver + return end @@ -399,13 +404,21 @@ function setup_point(solver::Solver{T}) where {T <: Real} if solver.reduce - # TODO don't find point / unnecessary stuff before reduce - solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z) - solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) + solver.time_inity = @elapsed handle_primal_eq(solver) + solver.status == PrimalInconsistent && return + init_y = update_primal_eq(solver, init_z) + + # solver.time_initx = @elapsed + init_x = find_initial_x(solver, init_s) else solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) - solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z) + + solver.time_inity = @elapsed handle_primal_eq(solver) + solver.status == PrimalInconsistent && return + init_y = update_primal_eq(solver, init_z) end + + # TODO delete if solver.status != SolveCalled # TODO due to inconsistent statuses return end diff --git a/src/Solvers/process.jl b/src/Solvers/process.jl index cc66c18f3..c9adc7620 100644 --- a/src/Solvers/process.jl +++ b/src/Solvers/process.jl @@ -111,17 +111,17 @@ function find_initial_x( else AG = vcat(A, G) end - if issparse(AG) + AG_fact = if issparse(AG) if !(T <: Float64) @warn("using dense factorization of [A; G] in preprocessing and " * "initial point finding because sparse factorization for number " * - "type $T is not supported by SuiteSparse packages") - AG_fact = qr!(Matrix(AG), ColumnNorm()) + "type $T is not supported by SuiteSparse packages", maxlog = 1) + qr!(Matrix(AG), ColumnNorm()) else - AG_fact = qr(AG, tol = solver.init_tol_qr) + qr(AG, tol = solver.init_tol_qr) end else - AG_fact = qr!(AG, ColumnNorm()) + qr!(AG, ColumnNorm()) end AG_rank = get_rank_est(AG_fact, solver.init_tol_qr) @@ -149,12 +149,13 @@ function find_initial_x( if !(AG_fact isa QRPivoted{T, Matrix{T}}) yz_sub = yz_sub[AG_fact.rpivinv] end - @views residual = norm(A' * yz_sub[1:p] + G' * - yz_sub[(p + 1):end] - model.c, Inf) - if residual > solver.init_tol_qr + residual = A' * yz_sub[1:p] + G' * yz_sub[(p + 1):end] - model.c + + @views res_norm = norm(residual, Inf) + if res_norm > solver.init_tol_qr if solver.verbose println("some dual equality constraints are inconsistent " * - "(residual $residual, tolerance $(solver.init_tol_qr))") + "(residual norm $res_norm, tolerance $(solver.init_tol_qr))") end solver.status = DualInconsistent return zeros(T, 0) @@ -179,15 +180,8 @@ function find_initial_x( return init_x end -# optionally preprocess primal equalities and solve for y as least squares -# solution to A'y = -c - G'z -function find_initial_y( - solver::Solver{T}, - init_z::Vector{T}, - ) where {T <: Real} - if solver.status != SolveCalled - return zeros(T, 0) - end + +function handle_primal_eq(solver::Solver{T}) where {T <: Real} model = solver.model p = model.p @@ -198,72 +192,53 @@ function find_initial_y( solver.Ap_Q = I return zeros(T, 0) end - solver.y_keep_idxs = 1:p - rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1) - if solver.init_use_indirect - # use indirect method - # TODO pick lsqr or lsmr - return IterativeSolvers.lsqr(model.A', rhs_y) - end + solver.init_use_indirect && return # factorize A' A = model.A - if issparse(A) + Ap_fact = solver.Ap_fact = if issparse(A) if !(T <: Float64) @warn("using dense factorization of A' in preprocessing and initial " * "point finding because sparse factorization for number type $T " * - "is not supported by SuiteSparse packages") - Ap_fact = qr!(Matrix(A'), ColumnNorm()) + "is not supported by SuiteSparse packages", maxlog = 1) + qr!(Matrix(A'), ColumnNorm()) else - Ap_fact = qr(sparse(A'), tol = solver.init_tol_qr) + qr(sparse(A'), tol = solver.init_tol_qr) end else - Ap_fact = qr!(Matrix(A'), ColumnNorm()) + qr!(Matrix(A'), ColumnNorm()) end + Ap_rank = solver.Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) - if solver.preprocess - return process_primal_eq(solver, Ap_fact, rhs_y) - end - - Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) - if Ap_rank < p - @warn("some primal equalities appear to be dependent " * - "(possibly inconsistent); try using preprocess = true") + if !solver.preprocess + if Ap_rank < p + @warn("some primal equalities appear to be dependent " * + "(possibly inconsistent); try using preprocess = true") + end + return end - return Ap_fact \ rhs_y -end - -# preprocess primal equalities -function process_primal_eq( - solver::Solver{T}, - Ap_fact::Factorization{T}, - rhs_y::Vector{T} - ) where {T <: Real} + # preprocess Ap_Q = Ap_fact.Q - Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) Ap_R = UpperTriangular(Ap_fact.R[1:Ap_rank, 1:Ap_rank]) col_piv = (Ap_fact isa QRPivoted{T, Matrix{T}}) ? Ap_fact.p : Ap_fact.pcol - - model = solver.model y_keep_idxs = col_piv[1:Ap_rank] - b_sub = model.b[y_keep_idxs] - if Ap_rank < model.p + if Ap_rank < p # some dependent primal equalities, so check if they are consistent # x_sub = Ap_Q * vcat((Ap_R' \ b_sub), zeros(n - Ap_rank)) x_sub = zeros(T, model.n) - x_sub1 = view(x_sub, 1:Ap_rank) - copyto!(x_sub1, b_sub) + @views x_sub1 = x_sub[1:Ap_rank] + @views copyto!(x_sub1, model.b[y_keep_idxs]) ldiv!(Ap_R', x_sub1) lmul!(Ap_Q, x_sub) - if !(Ap_fact isa QRPivoted{T, Matrix{T}}) x_sub = x_sub[Ap_fact.rpivinv] end residual = norm(model.A * x_sub - model.b, Inf) + if residual > solver.init_tol_qr if solver.verbose println("some primal equality constraints are inconsistent " * @@ -279,94 +254,51 @@ function process_primal_eq( end end + if !(solver.reduce && isa(model.G, MatrixyAG)) + # not reducing + # modify solver.model to remove/reorder some dual variables y + if Ap_fact isa QRPivoted{T, Matrix{T}} + model.A = model.A[y_keep_idxs, :] + else + row_piv = Ap_fact.prow + model.A = model.A[y_keep_idxs, row_piv] + model.G = model.G[:, row_piv] + solver.x_keep_idxs = solver.x_keep_idxs[row_piv] + end -# TODO refac some of reduce_primal_eq and below - if solver.reduce && isa(model.G, MatrixyAG) - return reduce_primal_eq(solver, Ap_fact, Ap_rank, Ap_R, y_keep_idxs, b_sub) + model.p = Ap_rank + solver.y_keep_idxs = y_keep_idxs + solver.Ap_R = Ap_R + solver.Ap_Q = Ap_Q + return end - # modify solver.model to remove/reorder some dual variables y - if !(Ap_fact isa QRPivoted{T, Matrix{T}}) - row_piv = Ap_fact.prow - model.A = model.A[y_keep_idxs, row_piv] - model.c = model.c[row_piv] - model.G = model.G[:, row_piv] - solver.x_keep_idxs = solver.x_keep_idxs[row_piv] + # reduce + if Ap_fact isa QRPivoted{T, Matrix{T}} + solver.reduce_row_piv_inv = Int[] else - model.A = model.A[y_keep_idxs, :] - end - model.b = b_sub - model.p = Ap_rank - solver.y_keep_idxs = y_keep_idxs - solver.Ap_R = Ap_R - solver.Ap_Q = Ap_Q - - - - # init_y = Ap_R \ ((Ap_Q' * (-c - G' * z))[1:Ap_rank]) - lmul!(Ap_Q', rhs_y) - init_y = rhs_y[1:Ap_rank] - ldiv!(Ap_R, init_y) - return init_y -end - -# remove all primal equalities i.e. let n = n0 - p0 and p = 0 -# TODO improve efficiency -# TODO avoid calculating GQ1 explicitly if possible -# recover original-space solution using: -# x0 = Q * [(R' \ b0), x] -# y0 = R \ (-cQ1' - GQ1' * z0) -function reduce_primal_eq( - solver::Solver{T}, - Ap_fact::Factorization{T}, - Ap_rank::Int, - Ap_R, - y_keep_idxs, - b_sub, - ) where {T <: Real} - model = solver.model - - if !(Ap_fact isa QRPivoted{T, Matrix{T}}) row_piv = Ap_fact.prow - model.c = model.c[row_piv] model.G = model.G[:, row_piv] solver.reduce_row_piv_inv = Ap_fact.rpivinv - else - solver.reduce_row_piv_inv = Int[] end - Q1_idxs = 1:Ap_rank - Q2_idxs = (Ap_rank + 1):model.n - - # [cQ1 cQ2] = c0' * Q + model.n -= Ap_rank Ap_Q = Ap_fact.Q - cQ = model.c' * Ap_Q - cQ1 = solver.reduce_cQ1 = cQ[Q1_idxs] - cQ2 = cQ[Q2_idxs] - # c = cQ2 - model.c = cQ2 - model.n = length(cQ2) - # offset = offset0 + cQ1 * (R' \ b0) - Rpib0 = solver.reduce_Rpib0 = ldiv!(Ap_R', b_sub) - # solver.Rpib0 = Rpib0 # TODO - model.obj_offset += dot(cQ1, Rpib0) # [GQ1 GQ2] = G0 * Q # very inefficient method used for sparse G * QRSparseQ # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818 - if model.G isa UniformScaling + G_mul = if model.G isa UniformScaling side = size(Ap_Q, 1) - G_mul = Matrix{T}(model.G, side, side) + Matrix{T}(model.G, side, side) elseif !isa(model.G, Matrix) - G_mul = Matrix{T}(model.G) + Matrix{T}(model.G) else - G_mul = model.G + model.G end GQ = G_mul * Ap_Q - GQ1 = solver.reduce_GQ1 = GQ[:, Q1_idxs] - GQ2 = GQ[:, Q2_idxs] - # h = h0 - GQ1 * (R' \ b0) - mul!(model.h, GQ1, Rpib0, -1, true) + solver.reduce_GQ1 = GQ[:, 1:Ap_rank] + GQ2 = GQ[:, (Ap_rank + 1):end] # G = GQ2 model.G = GQ2 @@ -374,17 +306,284 @@ function reduce_primal_eq( # A and b empty model.p = 0 model.A = zeros(T, 0, model.n) - model.b = zeros(T, 0) solver.reduce_Ap_R = Ap_R solver.reduce_Ap_Q = Ap_Q solver.reduce_y_keep_idxs = y_keep_idxs solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) solver.Ap_Q = I + return +end + +# update data for b, c, h from primal equality preprocessing/reduction +# and return initial y as least squares solution to A'y = -c - G'z +function update_primal_eq( + solver::Solver{T}, + init_z::Vector{T}, + ) where {T <: Real} + orig_model = solver.orig_model + iszero(orig_model.p) && return zeros(T, 0) + model = solver.model + + if solver.init_use_indirect + # use indirect method TODO pick lsqr or lsmr + rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1) + return IterativeSolvers.lsqr(model.A', rhs_y) + end + + Ap_fact = solver.Ap_fact + if !solver.preprocess + rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1) + return Ap_fact \ rhs_y + end + + # preprocess + Ap_rank = solver.Ap_rank + + if !(solver.reduce && isa(model.G, MatrixyAG)) + # not reducing + model.b = model.b[solver.y_keep_idxs] + if !(Ap_fact isa QRPivoted{T, Matrix{T}}) + model.c = model.c[Ap_fact.prow] + end + + # init_y = Ap_R \ ((Ap_Q' * (-c - G' * z))[1:Ap_rank]) + rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1) + lmul!(solver.Ap_Q', rhs_y) + init_y = rhs_y[1:Ap_rank] + ldiv!(solver.Ap_R, init_y) + return init_y + end + + # reduce + # [cQ1 cQ2] = c0' * Q + # c = cQ2 + cQ = model.c' * solver.reduce_Ap_Q + cQ1 = solver.reduce_cQ1 = cQ[1:Ap_rank] + model.c = cQ[(Ap_rank + 1):end] + + Rpib0 = solver.reduce_Rpib0 = model.b[solver.reduce_y_keep_idxs] + ldiv!(solver.reduce_Ap_R', Rpib0) + # offset = offset0 + cQ1 * (R' \ b0) + model.obj_offset += dot(cQ1, Rpib0) + + model.b = zeros(T, 0) + + GQ1 = solver.reduce_GQ1 + # h = h0 - GQ1 * (R' \ b0) + mul!(model.h, GQ1, Rpib0, -1, true) return zeros(T, 0) end + + +# # optionally preprocess primal equalities and solve for y as least squares +# # solution to A'y = -c - G'z +# function find_initial_y( +# solver::Solver{T}, +# init_z::Vector{T}, +# ) where {T <: Real} +# if solver.status != SolveCalled +# return zeros(T, 0) +# end + +# model = solver.model +# p = model.p +# if iszero(p) +# # y is empty (no primal variables) +# solver.y_keep_idxs = Int[] +# solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) +# solver.Ap_Q = I +# return zeros(T, 0) +# end +# solver.y_keep_idxs = 1:p + +# rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1) + +# if solver.init_use_indirect +# # use indirect method +# # TODO pick lsqr or lsmr +# return IterativeSolvers.lsqr(model.A', rhs_y) +# end + +# # factorize A' +# A = model.A +# Ap_fact = if issparse(A) +# if !(T <: Float64) +# @warn("using dense factorization of A' in preprocessing and initial " * +# "point finding because sparse factorization for number type $T " * +# "is not supported by SuiteSparse packages", maxlog = 1) +# qr!(Matrix(A'), ColumnNorm()) +# else +# qr(sparse(A'), tol = solver.init_tol_qr) +# end +# else +# qr!(Matrix(A'), ColumnNorm()) +# end + +# if solver.preprocess +# return process_primal_eq(solver, Ap_fact, rhs_y) +# end + +# Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) +# if Ap_rank < p +# @warn("some primal equalities appear to be dependent " * +# "(possibly inconsistent); try using preprocess = true") +# end + +# return Ap_fact \ rhs_y +# end + +# preprocess primal equalities +# function process_primal_eq( +# solver::Solver{T}, +# Ap_fact::Factorization{T}, +# rhs_y::Vector{T} +# ) where {T <: Real} +# Ap_Q = Ap_fact.Q +# Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) +# Ap_R = UpperTriangular(Ap_fact.R[1:Ap_rank, 1:Ap_rank]) +# col_piv = (Ap_fact isa QRPivoted{T, Matrix{T}}) ? Ap_fact.p : Ap_fact.pcol + +# model = solver.model +# y_keep_idxs = col_piv[1:Ap_rank] +# b_sub = model.b[y_keep_idxs] + +# if Ap_rank < model.p +# # some dependent primal equalities, so check if they are consistent +# # x_sub = Ap_Q * vcat((Ap_R' \ b_sub), zeros(n - Ap_rank)) +# x_sub = zeros(T, model.n) +# x_sub1 = view(x_sub, 1:Ap_rank) +# copyto!(x_sub1, b_sub) +# ldiv!(Ap_R', x_sub1) +# lmul!(Ap_Q, x_sub) + +# if !(Ap_fact isa QRPivoted{T, Matrix{T}}) +# x_sub = x_sub[Ap_fact.rpivinv] +# end +# residual = norm(model.A * x_sub - model.b, Inf) +# if residual > solver.init_tol_qr +# if solver.verbose +# println("some primal equality constraints are inconsistent " * +# "(residual $residual, tolerance $(solver.init_tol_qr))") +# end +# solver.status = PrimalInconsistent +# return zeros(T, 0) +# end +# if solver.verbose +# p = model.p +# println("$(p - Ap_rank) of $p primal equality constraints " * +# "are dependent") +# end +# end + + +# # TODO refac some of reduce_primal_eq and below +# if solver.reduce && isa(model.G, MatrixyAG) +# return reduce_primal_eq(solver, Ap_fact, Ap_rank, Ap_R, y_keep_idxs, b_sub) +# end + +# # modify solver.model to remove/reorder some dual variables y +# if !(Ap_fact isa QRPivoted{T, Matrix{T}}) +# row_piv = Ap_fact.prow +# model.A = model.A[y_keep_idxs, row_piv] +# model.c = model.c[row_piv] +# model.G = model.G[:, row_piv] +# solver.x_keep_idxs = solver.x_keep_idxs[row_piv] +# else +# model.A = model.A[y_keep_idxs, :] +# end +# model.b = b_sub +# model.p = Ap_rank +# solver.y_keep_idxs = y_keep_idxs +# solver.Ap_R = Ap_R +# solver.Ap_Q = Ap_Q + + + +# # init_y = Ap_R \ ((Ap_Q' * (-c - G' * z))[1:Ap_rank]) +# lmul!(Ap_Q', rhs_y) +# init_y = rhs_y[1:Ap_rank] +# ldiv!(Ap_R, init_y) +# return init_y +# end + +# remove all primal equalities i.e. let n = n0 - p0 and p = 0 +# TODO improve efficiency +# TODO avoid calculating GQ1 explicitly if possible +# recover original-space solution using: +# x0 = Q * [(R' \ b0), x] +# y0 = R \ (-cQ1' - GQ1' * z0) +# function reduce_primal_eq( +# solver::Solver{T}, +# Ap_fact::Factorization{T}, +# Ap_rank::Int, +# Ap_R, +# y_keep_idxs, +# b_sub, +# ) where {T <: Real} +# model = solver.model + +# if !(Ap_fact isa QRPivoted{T, Matrix{T}}) +# row_piv = Ap_fact.prow +# model.c = model.c[row_piv] +# model.G = model.G[:, row_piv] +# solver.reduce_row_piv_inv = Ap_fact.rpivinv +# else +# solver.reduce_row_piv_inv = Int[] +# end + +# Q1_idxs = 1:Ap_rank +# Q2_idxs = (Ap_rank + 1):model.n + +# # [cQ1 cQ2] = c0' * Q +# Ap_Q = Ap_fact.Q +# cQ = model.c' * Ap_Q +# cQ1 = solver.reduce_cQ1 = cQ[Q1_idxs] +# cQ2 = cQ[Q2_idxs] +# # c = cQ2 +# model.c = cQ2 +# model.n = length(cQ2) +# # offset = offset0 + cQ1 * (R' \ b0) +# Rpib0 = solver.reduce_Rpib0 = ldiv!(Ap_R', b_sub) +# # solver.Rpib0 = Rpib0 # TODO +# model.obj_offset += dot(cQ1, Rpib0) + +# # [GQ1 GQ2] = G0 * Q +# # very inefficient method used for sparse G * QRSparseQ +# # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818 +# if model.G isa UniformScaling +# side = size(Ap_Q, 1) +# G_mul = Matrix{T}(model.G, side, side) +# elseif !isa(model.G, Matrix) +# G_mul = Matrix{T}(model.G) +# else +# G_mul = model.G +# end +# GQ = G_mul * Ap_Q +# GQ1 = solver.reduce_GQ1 = GQ[:, Q1_idxs] +# GQ2 = GQ[:, Q2_idxs] +# # h = h0 - GQ1 * (R' \ b0) +# mul!(model.h, GQ1, Rpib0, -1, true) + +# # G = GQ2 +# model.G = GQ2 + +# # A and b empty +# model.p = 0 +# model.A = zeros(T, 0, model.n) +# model.b = zeros(T, 0) +# solver.reduce_Ap_R = Ap_R +# solver.reduce_Ap_Q = Ap_Q + +# solver.reduce_y_keep_idxs = y_keep_idxs +# solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) +# solver.Ap_Q = I + +# return zeros(T, 0) +# end + # (pivoted) QR factorizations are usually rank-revealing but may be unreliable # see http://www.math.sjsu.edu/~foster/rankrevealingcode.html # TODO could replace this with rank(qr_fact) when available for both dense and sparse diff --git a/src/Solvers/systemsolvers/qrchol.jl b/src/Solvers/systemsolvers/qrchol.jl index c13661d1e..c88c99148 100644 --- a/src/Solvers/systemsolvers/qrchol.jl +++ b/src/Solvers/systemsolvers/qrchol.jl @@ -151,7 +151,11 @@ function load( # very inefficient method used for sparse G * QRSparseQ # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818 - GQ = model.G * solver.Ap_Q + GQ = if solver.Ap_Q == I + model.G + else + model.G * solver.Ap_Q + end syssolver.GQ2 = GQ[:, (p + 1):end] syssolver.HGQ2 = zeros(T, q, nmp) diff --git a/test/runnativetests.jl b/test/runnativetests.jl index 5bbb26c95..88b4f6c1d 100644 --- a/test/runnativetests.jl +++ b/test/runnativetests.jl @@ -23,7 +23,7 @@ all_reals = [ ] diff_reals = [ Float64, - BigFloat, + # BigFloat, ] string_nameof(T) = string(nameof(T)) @@ -51,8 +51,8 @@ end inst_defaults = vcat( inst_preproc, inst_infeas, - # inst_cones_few, - inst_cones_many, + inst_cones_few, + # # inst_cones_many, ) for inst_name in inst_defaults test_instance_solver(inst_name, Float64, default_options) @@ -62,8 +62,10 @@ end @testset "no preprocess tests" begin println("\nstarting no preprocess tests") for inst_name in inst_cones_few, T in diff_reals - options = (; default_options..., preprocess = false, reduce = false, - syssolver = Solvers.SymIndefDenseSystemSolver{T}()) + # options = (; default_options..., preprocess = false, reduce = false, + # syssolver = Solvers.SymIndefDenseSystemSolver{T}()) + options = (; default_options..., preprocess = true, reduce = false, + syssolver = Solvers.SymIndefDenseSystemSolver{T}()) test_instance_solver(inst_name, T, options) end end @@ -80,79 +82,79 @@ end end end -@testset "system solvers tests" begin - println("\nstarting system solvers tests") - syssolvers = [ - (Solvers.NaiveDenseSystemSolver, diff_reals), - (Solvers.NaiveSparseSystemSolver, [Float64,]), - (Solvers.NaiveElimDenseSystemSolver, diff_reals), - (Solvers.NaiveElimSparseSystemSolver, [Float64,]), - (Solvers.SymIndefDenseSystemSolver, all_reals), - (Solvers.SymIndefSparseSystemSolver, [Float64,]), - (Solvers.QRCholDenseSystemSolver, all_reals), - ] - for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, - T in real_types - options = (; default_options..., syssolver = syssolver{T}(), - reduce = false) - test_instance_solver(inst_name, T, options, string_nameof(syssolver)) - end -end - -@testset "PredOrCentStepper tests" begin - verbose = true - println("\nstarting PredOrCentStepper tests (with printing)") - - # adjustment and curve search - use_adj_curv = [(false, false), (true, false), (true, true)] - for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals - stepper = Solvers.PredOrCentStepper{T}(; - use_adjustment = adj, use_curve_search = curv) - options = (; default_options..., verbose = verbose, stepper = stepper) - test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") - end - - # other options - for inst_name in inst_minimal - T = Float64 - stepper = Solvers.PredOrCentStepper{T}(; - # stepper options - use_adjustment = false, use_curve_search = false, - max_cent_steps = 8, pred_prox_bound = 0.0332, - # searcher options - min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, - alpha_sched = [0.9999 * 0.7^i for i in 0:22]) - options = (; default_options..., verbose = verbose, stepper = stepper) - test_instance_solver(inst_name, T, options, "other") - end -end - -@testset "CombinedStepper tests" begin - verbose = true - println("\nstarting CombinedStepper tests (with printing)") - shifts = [0, 2] - for inst_name in inst_minimal, shift in shifts, T in diff_reals - options = (; default_options..., verbose = verbose, - stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) - test_instance_solver(inst_name, T, options, "shift=$shift") - end -end - -@testset "model modification tests" begin - println("\nstarting model modification tests") - syssolvers = [ - (Solvers.SymIndefDenseSystemSolver, diff_reals), - # (Solvers.SymIndefSparseSystemSolver, [Float64,]), - # (Solvers.QRCholDenseSystemSolver, diff_reals), - ] - for inst_name in inst_modify, (syssolver, real_types) in syssolvers, - T in real_types, rescale in (false, true) - # TODO later also rescale/preproc/reduce true - options = (; default_options..., syssolver = syssolver{T}(), - rescale = rescale, preprocess = false, reduce = false) - test_instance_solver(inst_name, T, options) - end -end +# @testset "system solvers tests" begin +# println("\nstarting system solvers tests") +# syssolvers = [ +# (Solvers.NaiveDenseSystemSolver, diff_reals), +# (Solvers.NaiveSparseSystemSolver, [Float64,]), +# (Solvers.NaiveElimDenseSystemSolver, diff_reals), +# (Solvers.NaiveElimSparseSystemSolver, [Float64,]), +# (Solvers.SymIndefDenseSystemSolver, all_reals), +# (Solvers.SymIndefSparseSystemSolver, [Float64,]), +# (Solvers.QRCholDenseSystemSolver, all_reals), +# ] +# for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, +# T in real_types +# options = (; default_options..., syssolver = syssolver{T}(), +# reduce = false) +# test_instance_solver(inst_name, T, options, string_nameof(syssolver)) +# end +# end + +# @testset "PredOrCentStepper tests" begin +# verbose = true +# println("\nstarting PredOrCentStepper tests (with printing)") + +# # adjustment and curve search +# use_adj_curv = [(false, false), (true, false), (true, true)] +# for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals +# stepper = Solvers.PredOrCentStepper{T}(; +# use_adjustment = adj, use_curve_search = curv) +# options = (; default_options..., verbose = verbose, stepper = stepper) +# test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") +# end + +# # other options +# for inst_name in inst_minimal +# T = Float64 +# stepper = Solvers.PredOrCentStepper{T}(; +# # stepper options +# use_adjustment = false, use_curve_search = false, +# max_cent_steps = 8, pred_prox_bound = 0.0332, +# # searcher options +# min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, +# alpha_sched = [0.9999 * 0.7^i for i in 0:22]) +# options = (; default_options..., verbose = verbose, stepper = stepper) +# test_instance_solver(inst_name, T, options, "other") +# end +# end + +# @testset "CombinedStepper tests" begin +# verbose = true +# println("\nstarting CombinedStepper tests (with printing)") +# shifts = [0, 2] +# for inst_name in inst_minimal, shift in shifts, T in diff_reals +# options = (; default_options..., verbose = verbose, +# stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) +# test_instance_solver(inst_name, T, options, "shift=$shift") +# end +# end + +# @testset "model modification tests" begin +# println("\nstarting model modification tests") +# syssolvers = [ +# (Solvers.SymIndefDenseSystemSolver, diff_reals), +# # (Solvers.SymIndefSparseSystemSolver, [Float64,]), +# # (Solvers.QRCholDenseSystemSolver, diff_reals), +# ] +# for inst_name in inst_modify, (syssolver, real_types) in syssolvers, +# T in real_types, rescale in (false, true) +# # TODO later also rescale/preproc/reduce true +# options = (; default_options..., syssolver = syssolver{T}(), +# rescale = rescale, preprocess = false, reduce = false) +# test_instance_solver(inst_name, T, options) +# end +# end end ; From 6f900f33e34acf0c5d10f9d327e20c354b17d57c Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 12 Dec 2021 11:04:44 -0500 Subject: [PATCH 4/9] start cleanup of dual eq preproc --- src/Solvers/Solvers.jl | 13 +++- src/Solvers/process.jl | 140 ++++++++++++++++++++++++++++++++++++++--- test/runnativetests.jl | 46 +++++++------- 3 files changed, 166 insertions(+), 33 deletions(-) diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 9e282f668..6db010fe8 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -115,6 +115,8 @@ mutable struct Solver{T <: Real} Ap_rank::Int Ap_R::UpperTriangular{T, <:AbstractMatrix{T}} Ap_Q::Union{UniformScaling, AbstractMatrix{T}} + AG_fact::Factorization{T} + AG_R::UpperTriangular{T, <:AbstractMatrix{T}} reduce_cQ1 reduce_Rpib0 reduce_GQ1 @@ -408,10 +410,15 @@ function setup_point(solver::Solver{T}) where {T <: Real} solver.status == PrimalInconsistent && return init_y = update_primal_eq(solver, init_z) - # solver.time_initx = @elapsed - init_x = find_initial_x(solver, init_s) + solver.time_initx = @elapsed handle_dual_eq(solver) + solver.status == DualInconsistent && return + init_x = update_dual_eq(solver, init_s) + # init_x = find_initial_x(solver, init_s) else - solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) + # solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) + solver.time_initx = @elapsed handle_dual_eq(solver) + solver.status == DualInconsistent && return + init_x = update_dual_eq(solver, init_s) solver.time_inity = @elapsed handle_primal_eq(solver) solver.status == PrimalInconsistent && return diff --git a/src/Solvers/process.jl b/src/Solvers/process.jl index c9adc7620..46d8f5b6c 100644 --- a/src/Solvers/process.jl +++ b/src/Solvers/process.jl @@ -135,9 +135,9 @@ function find_initial_x( end # preprocess dual equalities + AG_R = UpperTriangular(AG_fact.R[1:AG_rank, 1:AG_rank]) col_piv = (AG_fact isa QRPivoted{T, Matrix{T}}) ? AG_fact.p : AG_fact.pcol x_keep_idxs = col_piv[1:AG_rank] - AG_R = UpperTriangular(AG_fact.R[1:AG_rank, 1:AG_rank]) c_sub = model.c[x_keep_idxs] # yz_sub = AG_fact.Q * vcat((AG_R' \ c_sub), zeros(p + q - AG_rank)) @@ -164,7 +164,7 @@ function find_initial_x( println("$(n - AG_rank) of $n dual equality constraints are dependent") end - # modify solver.model to remove/reorder some primal variables x + # modify model to remove/reorder some primal variables x model.c = c_sub model.A = A[:, x_keep_idxs] model.G = G[:, x_keep_idxs] @@ -181,6 +181,133 @@ function find_initial_x( end +function handle_dual_eq(solver::Solver{T}) where {T <: Real} + model = solver.model + n = model.n + + if iszero(n) # x is empty (no primal variables) + solver.x_keep_idxs = Int[] + return zeros(T, 0) + end + solver.x_keep_idxs = 1:n + + A = model.A + G = model.G + if solver.init_use_indirect || !isa(A, MatrixyAG) || !isa(G, MatrixyAG) + return + end + + # direct method + AG = if iszero(model.p) + # A is empty + if issparse(G) + G + elseif G isa Matrix{T} + copy(G) + else + Matrix(G) + end + else + vcat(A, G) + end + AG_fact = solver.AG_fact = if issparse(AG) + if !(T <: Float64) + @warn("using dense factorization of [A; G] in preprocessing and " * + "initial point finding because sparse factorization for number " * + "type $T is not supported by SuiteSparse packages", maxlog = 1) + qr!(Matrix(AG), ColumnNorm()) + else + qr(AG, tol = solver.init_tol_qr) + end + else + qr!(AG, ColumnNorm()) + end + AG_rank = get_rank_est(AG_fact, solver.init_tol_qr) + + if !solver.preprocess || (AG_rank == n) + if AG_rank < n + @warn("some dual equalities appear to be dependent " * + "(possibly inconsistent); try using preprocess = true") + end + return + end + + # preprocess + AG_R = solver.AG_R = UpperTriangular(AG_fact.R[1:AG_rank, 1:AG_rank]) + col_piv = (AG_fact isa QRPivoted{T, Matrix{T}}) ? AG_fact.p : AG_fact.pcol + x_keep_idxs = col_piv[1:AG_rank] + + # yz_sub = AG_fact.Q * vcat((AG_R' \ c_sub), zeros(p + q - AG_rank)) + yz_sub = zeros(T, p + q) + yz_sub1 = view(yz_sub, 1:AG_rank) + @views copyto!(yz_sub1, model.c[x_keep_idxs]) + ldiv!(AG_R', yz_sub1) + lmul!(AG_fact.Q, yz_sub) + if !(AG_fact isa QRPivoted{T, Matrix{T}}) + yz_sub = yz_sub[AG_fact.rpivinv] + end + residual = A' * yz_sub[1:p] + G' * yz_sub[(p + 1):end] - model.c + + @views res_norm = norm(residual, Inf) + if res_norm > solver.init_tol_qr + if solver.verbose + println("some dual equality constraints are inconsistent " * + "(residual norm $res_norm, tolerance $(solver.init_tol_qr))") + end + solver.status = DualInconsistent + return zeros(T, 0) + end + if solver.verbose + println("$(n - AG_rank) of $n dual equality constraints are dependent") + end + + # modify model to remove/reorder some primal variables x + model.A = A[:, x_keep_idxs] + model.G = G[:, x_keep_idxs] + model.n = AG_rank + solver.x_keep_idxs = x_keep_idxs + return +end + +# update data for b, c, h from dual equality preprocessing/reduction +# and return initial x as least squares solution to Ax = b, Gx = h - s +function update_dual_eq( + solver::Solver{T}, + init_s::Vector{T}, + ) where {T <: Real} + model = solver.model + + if solver.init_use_indirect + # use indirect method TODO pick lsqr or lsmr + rhs_x = vcat(model.b, model.h - init_s) + if iszero(model.p) + AG = model.G + else + linmap(M::UniformScaling) = LinearMaps.LinearMap(M, model.n) + linmap(M) = LinearMaps.LinearMap(M) + AG = vcat(linmap(model.A), linmap(model.G)) + end + return IterativeSolvers.lsqr(AG, rhs_x) + end + + AG_fact = solver.AG_fact + if !solver.preprocess || (AG_rank == model.n) + rhs_x = vcat(model.b, model.h - init_s) + return AG_fact \ rhs_x + end + + model.c = model.c[solver.x_keep_idxs] + + # init_x = AG_R \ ((AG_fact.Q' * vcat(b, h - point.s))[1:AG_rank]) + temp = vcat(model.b, model.h - init_s) + lmul!(AG_fact.Q', temp) + init_x = temp[1:model.n] + ldiv!(AG_R, init_x) + + return init_x +end + + function handle_primal_eq(solver::Solver{T}) where {T <: Real} model = solver.model p = model.p @@ -256,7 +383,7 @@ function handle_primal_eq(solver::Solver{T}) where {T <: Real} if !(solver.reduce && isa(model.G, MatrixyAG)) # not reducing - # modify solver.model to remove/reorder some dual variables y + # modify model to remove/reorder some dual variables y if Ap_fact isa QRPivoted{T, Matrix{T}} model.A = model.A[y_keep_idxs, :] else @@ -321,8 +448,7 @@ function update_primal_eq( solver::Solver{T}, init_z::Vector{T}, ) where {T <: Real} - orig_model = solver.orig_model - iszero(orig_model.p) && return zeros(T, 0) + iszero(solver.orig_model.p) && return zeros(T, 0) model = solver.model if solver.init_use_indirect @@ -365,7 +491,7 @@ function update_primal_eq( Rpib0 = solver.reduce_Rpib0 = model.b[solver.reduce_y_keep_idxs] ldiv!(solver.reduce_Ap_R', Rpib0) # offset = offset0 + cQ1 * (R' \ b0) - model.obj_offset += dot(cQ1, Rpib0) + model.obj_offset = solver.orig_model.obj_offset + dot(cQ1, Rpib0) model.b = zeros(T, 0) @@ -484,7 +610,7 @@ end # return reduce_primal_eq(solver, Ap_fact, Ap_rank, Ap_R, y_keep_idxs, b_sub) # end -# # modify solver.model to remove/reorder some dual variables y +# # modify model to remove/reorder some dual variables y # if !(Ap_fact isa QRPivoted{T, Matrix{T}}) # row_piv = Ap_fact.prow # model.A = model.A[y_keep_idxs, row_piv] diff --git a/test/runnativetests.jl b/test/runnativetests.jl index 88b4f6c1d..9879fbe8b 100644 --- a/test/runnativetests.jl +++ b/test/runnativetests.jl @@ -46,18 +46,18 @@ end @testset "native tests" begin -@testset "default options tests" begin - println("starting default options tests") - inst_defaults = vcat( - inst_preproc, - inst_infeas, - inst_cones_few, - # # inst_cones_many, - ) - for inst_name in inst_defaults - test_instance_solver(inst_name, Float64, default_options) - end -end +# @testset "default options tests" begin +# println("starting default options tests") +# inst_defaults = vcat( +# inst_preproc, +# # inst_infeas, +# # inst_cones_few, +# # # inst_cones_many, +# ) +# for inst_name in inst_defaults +# test_instance_solver(inst_name, Float64, default_options) +# end +# end @testset "no preprocess tests" begin println("\nstarting no preprocess tests") @@ -70,17 +70,17 @@ end end end -@testset "indirect solvers tests" begin - println("\nstarting indirect solvers tests") - for inst_name in inst_indirect, T in diff_reals - options = (; default_options..., init_use_indirect = true, - preprocess = false, reduce = false, - syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), - tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, - tol_infeas = 1e-6) - test_instance_solver(inst_name, T, options) - end -end +# @testset "indirect solvers tests" begin +# println("\nstarting indirect solvers tests") +# for inst_name in inst_indirect, T in diff_reals +# options = (; default_options..., init_use_indirect = true, +# preprocess = false, reduce = false, +# syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), +# tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, +# tol_infeas = 1e-6) +# test_instance_solver(inst_name, T, options) +# end +# end # @testset "system solvers tests" begin # println("\nstarting system solvers tests") From 577323c027dec0bf2a011a3aab4d4b5def0cbe67 Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 12 Dec 2021 11:38:48 -0500 Subject: [PATCH 5/9] fix preproc bugs --- src/Solvers/Solvers.jl | 10 -- src/Solvers/process.jl | 366 +++-------------------------------------- test/runnativetests.jl | 90 +++++----- 3 files changed, 64 insertions(+), 402 deletions(-) diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 6db010fe8..1b403f0c5 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -413,9 +413,7 @@ function setup_point(solver::Solver{T}) where {T <: Real} solver.time_initx = @elapsed handle_dual_eq(solver) solver.status == DualInconsistent && return init_x = update_dual_eq(solver, init_s) - # init_x = find_initial_x(solver, init_s) else - # solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s) solver.time_initx = @elapsed handle_dual_eq(solver) solver.status == DualInconsistent && return init_x = update_dual_eq(solver, init_s) @@ -425,14 +423,6 @@ function setup_point(solver::Solver{T}) where {T <: Real} init_y = update_primal_eq(solver, init_z) end - # TODO delete - if solver.status != SolveCalled # TODO due to inconsistent statuses - return - end - - - - model = solver.model point = solver.point = Point(model) point.x .= init_x diff --git a/src/Solvers/process.jl b/src/Solvers/process.jl index 46d8f5b6c..f585b04cd 100644 --- a/src/Solvers/process.jl +++ b/src/Solvers/process.jl @@ -61,135 +61,11 @@ function rescale_data(solver::Solver{T}) where {T <: Real} return true end -# optionally preprocess dual equalities and solve for x as least squares -# solution to Ax = b, Gx = h - s -function find_initial_x( - solver::Solver{T}, - init_s::Vector{T}, - ) where {T <: Real} - if solver.status != SolveCalled - return zeros(T, 0) - end - model = solver.model - n = model.n - if iszero(n) # x is empty (no primal variables) - solver.x_keep_idxs = Int[] - return zeros(T, 0) - end - p = model.p - q = model.q - A = model.A - G = model.G - solver.x_keep_idxs = 1:n - - rhs = vcat(model.b, model.h - init_s) - - # indirect method - if solver.init_use_indirect || !isa(A, MatrixyAG) || !isa(G, MatrixyAG) - # TODO pick lsqr or lsmr - if iszero(p) - AG = G - else - linmap(M::UniformScaling) = LinearMaps.LinearMap(M, n) - linmap(M) = LinearMaps.LinearMap(M) - AG = vcat(linmap(A), linmap(G)) - end - init_x = IterativeSolvers.lsqr(AG, rhs) - return init_x - end - - # direct method - if iszero(p) - # A is empty - if issparse(G) - AG = G - elseif G isa Matrix{T} - AG = copy(G) - else - AG = Matrix(G) - end - else - AG = vcat(A, G) - end - AG_fact = if issparse(AG) - if !(T <: Float64) - @warn("using dense factorization of [A; G] in preprocessing and " * - "initial point finding because sparse factorization for number " * - "type $T is not supported by SuiteSparse packages", maxlog = 1) - qr!(Matrix(AG), ColumnNorm()) - else - qr(AG, tol = solver.init_tol_qr) - end - else - qr!(AG, ColumnNorm()) - end - AG_rank = get_rank_est(AG_fact, solver.init_tol_qr) - - if !solver.preprocess || (AG_rank == n) - if AG_rank < n - @warn("some dual equalities appear to be dependent " * - "(possibly inconsistent); try using preprocess = true") - end - init_x = AG_fact \ rhs - return init_x - end - - # preprocess dual equalities - AG_R = UpperTriangular(AG_fact.R[1:AG_rank, 1:AG_rank]) - col_piv = (AG_fact isa QRPivoted{T, Matrix{T}}) ? AG_fact.p : AG_fact.pcol - x_keep_idxs = col_piv[1:AG_rank] - - c_sub = model.c[x_keep_idxs] - # yz_sub = AG_fact.Q * vcat((AG_R' \ c_sub), zeros(p + q - AG_rank)) - yz_sub = zeros(T, p + q) - yz_sub1 = view(yz_sub, 1:AG_rank) - copyto!(yz_sub1, c_sub) - ldiv!(AG_R', yz_sub1) - lmul!(AG_fact.Q, yz_sub) - if !(AG_fact isa QRPivoted{T, Matrix{T}}) - yz_sub = yz_sub[AG_fact.rpivinv] - end - residual = A' * yz_sub[1:p] + G' * yz_sub[(p + 1):end] - model.c - - @views res_norm = norm(residual, Inf) - if res_norm > solver.init_tol_qr - if solver.verbose - println("some dual equality constraints are inconsistent " * - "(residual norm $res_norm, tolerance $(solver.init_tol_qr))") - end - solver.status = DualInconsistent - return zeros(T, 0) - end - if solver.verbose - println("$(n - AG_rank) of $n dual equality constraints are dependent") - end - - # modify model to remove/reorder some primal variables x - model.c = c_sub - model.A = A[:, x_keep_idxs] - model.G = G[:, x_keep_idxs] - model.n = AG_rank - solver.x_keep_idxs = x_keep_idxs - - # init_x = AG_R \ ((AG_fact.Q' * vcat(b, h - point.s))[1:AG_rank]) - temp = vcat(model.b, model.h - init_s) - lmul!(AG_fact.Q', temp) - init_x = temp[1:model.n] - ldiv!(AG_R, init_x) - - return init_x -end - - function handle_dual_eq(solver::Solver{T}) where {T <: Real} model = solver.model n = model.n - - if iszero(n) # x is empty (no primal variables) - solver.x_keep_idxs = Int[] - return zeros(T, 0) - end solver.x_keep_idxs = 1:n + iszero(n) && return A = model.A G = model.G @@ -238,7 +114,8 @@ function handle_dual_eq(solver::Solver{T}) where {T <: Real} x_keep_idxs = col_piv[1:AG_rank] # yz_sub = AG_fact.Q * vcat((AG_R' \ c_sub), zeros(p + q - AG_rank)) - yz_sub = zeros(T, p + q) + p = model.p + yz_sub = zeros(T, p + model.q) yz_sub1 = view(yz_sub, 1:AG_rank) @views copyto!(yz_sub1, model.c[x_keep_idxs]) ldiv!(AG_R', yz_sub1) @@ -246,9 +123,13 @@ function handle_dual_eq(solver::Solver{T}) where {T <: Real} if !(AG_fact isa QRPivoted{T, Matrix{T}}) yz_sub = yz_sub[AG_fact.rpivinv] end - residual = A' * yz_sub[1:p] + G' * yz_sub[(p + 1):end] - model.c - @views res_norm = norm(residual, Inf) + # residual = A' * yz_sub[1:p] + G' * yz_sub[(p + 1):end] - model.c + residual = copy(model.c) + @views mul!(residual, G', yz_sub[(p + 1):end], true, -1) + @views mul!(residual, A', yz_sub[1:p], true, true) + res_norm = norm(residual, Inf) + if res_norm > solver.init_tol_qr if solver.verbose println("some dual equality constraints are inconsistent " * @@ -277,7 +158,11 @@ function update_dual_eq( ) where {T <: Real} model = solver.model - if solver.init_use_indirect + model.c = model.c[solver.x_keep_idxs] + iszero(model.n) && return zeros(T, 0) + + if solver.init_use_indirect || + !isa(model.A, MatrixyAG) || !isa(model.G, MatrixyAG) # use indirect method TODO pick lsqr or lsmr rhs_x = vcat(model.b, model.h - init_s) if iszero(model.p) @@ -290,24 +175,21 @@ function update_dual_eq( return IterativeSolvers.lsqr(AG, rhs_x) end + rhs_x = vcat(model.b, model.h - init_s) AG_fact = solver.AG_fact - if !solver.preprocess || (AG_rank == model.n) - rhs_x = vcat(model.b, model.h - init_s) + + if !solver.preprocess || (model.n == size(AG_fact, 2)) return AG_fact \ rhs_x end - model.c = model.c[solver.x_keep_idxs] - # init_x = AG_R \ ((AG_fact.Q' * vcat(b, h - point.s))[1:AG_rank]) - temp = vcat(model.b, model.h - init_s) - lmul!(AG_fact.Q', temp) - init_x = temp[1:model.n] - ldiv!(AG_R, init_x) + lmul!(AG_fact.Q', rhs_x) + init_x = rhs_x[1:model.n] + ldiv!(solver.AG_R, init_x) return init_x end - function handle_primal_eq(solver::Solver{T}) where {T <: Real} model = solver.model p = model.p @@ -502,214 +384,6 @@ function update_primal_eq( return zeros(T, 0) end - - -# # optionally preprocess primal equalities and solve for y as least squares -# # solution to A'y = -c - G'z -# function find_initial_y( -# solver::Solver{T}, -# init_z::Vector{T}, -# ) where {T <: Real} -# if solver.status != SolveCalled -# return zeros(T, 0) -# end - -# model = solver.model -# p = model.p -# if iszero(p) -# # y is empty (no primal variables) -# solver.y_keep_idxs = Int[] -# solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) -# solver.Ap_Q = I -# return zeros(T, 0) -# end -# solver.y_keep_idxs = 1:p - -# rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1) - -# if solver.init_use_indirect -# # use indirect method -# # TODO pick lsqr or lsmr -# return IterativeSolvers.lsqr(model.A', rhs_y) -# end - -# # factorize A' -# A = model.A -# Ap_fact = if issparse(A) -# if !(T <: Float64) -# @warn("using dense factorization of A' in preprocessing and initial " * -# "point finding because sparse factorization for number type $T " * -# "is not supported by SuiteSparse packages", maxlog = 1) -# qr!(Matrix(A'), ColumnNorm()) -# else -# qr(sparse(A'), tol = solver.init_tol_qr) -# end -# else -# qr!(Matrix(A'), ColumnNorm()) -# end - -# if solver.preprocess -# return process_primal_eq(solver, Ap_fact, rhs_y) -# end - -# Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) -# if Ap_rank < p -# @warn("some primal equalities appear to be dependent " * -# "(possibly inconsistent); try using preprocess = true") -# end - -# return Ap_fact \ rhs_y -# end - -# preprocess primal equalities -# function process_primal_eq( -# solver::Solver{T}, -# Ap_fact::Factorization{T}, -# rhs_y::Vector{T} -# ) where {T <: Real} -# Ap_Q = Ap_fact.Q -# Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr) -# Ap_R = UpperTriangular(Ap_fact.R[1:Ap_rank, 1:Ap_rank]) -# col_piv = (Ap_fact isa QRPivoted{T, Matrix{T}}) ? Ap_fact.p : Ap_fact.pcol - -# model = solver.model -# y_keep_idxs = col_piv[1:Ap_rank] -# b_sub = model.b[y_keep_idxs] - -# if Ap_rank < model.p -# # some dependent primal equalities, so check if they are consistent -# # x_sub = Ap_Q * vcat((Ap_R' \ b_sub), zeros(n - Ap_rank)) -# x_sub = zeros(T, model.n) -# x_sub1 = view(x_sub, 1:Ap_rank) -# copyto!(x_sub1, b_sub) -# ldiv!(Ap_R', x_sub1) -# lmul!(Ap_Q, x_sub) - -# if !(Ap_fact isa QRPivoted{T, Matrix{T}}) -# x_sub = x_sub[Ap_fact.rpivinv] -# end -# residual = norm(model.A * x_sub - model.b, Inf) -# if residual > solver.init_tol_qr -# if solver.verbose -# println("some primal equality constraints are inconsistent " * -# "(residual $residual, tolerance $(solver.init_tol_qr))") -# end -# solver.status = PrimalInconsistent -# return zeros(T, 0) -# end -# if solver.verbose -# p = model.p -# println("$(p - Ap_rank) of $p primal equality constraints " * -# "are dependent") -# end -# end - - -# # TODO refac some of reduce_primal_eq and below -# if solver.reduce && isa(model.G, MatrixyAG) -# return reduce_primal_eq(solver, Ap_fact, Ap_rank, Ap_R, y_keep_idxs, b_sub) -# end - -# # modify model to remove/reorder some dual variables y -# if !(Ap_fact isa QRPivoted{T, Matrix{T}}) -# row_piv = Ap_fact.prow -# model.A = model.A[y_keep_idxs, row_piv] -# model.c = model.c[row_piv] -# model.G = model.G[:, row_piv] -# solver.x_keep_idxs = solver.x_keep_idxs[row_piv] -# else -# model.A = model.A[y_keep_idxs, :] -# end -# model.b = b_sub -# model.p = Ap_rank -# solver.y_keep_idxs = y_keep_idxs -# solver.Ap_R = Ap_R -# solver.Ap_Q = Ap_Q - - - -# # init_y = Ap_R \ ((Ap_Q' * (-c - G' * z))[1:Ap_rank]) -# lmul!(Ap_Q', rhs_y) -# init_y = rhs_y[1:Ap_rank] -# ldiv!(Ap_R, init_y) -# return init_y -# end - -# remove all primal equalities i.e. let n = n0 - p0 and p = 0 -# TODO improve efficiency -# TODO avoid calculating GQ1 explicitly if possible -# recover original-space solution using: -# x0 = Q * [(R' \ b0), x] -# y0 = R \ (-cQ1' - GQ1' * z0) -# function reduce_primal_eq( -# solver::Solver{T}, -# Ap_fact::Factorization{T}, -# Ap_rank::Int, -# Ap_R, -# y_keep_idxs, -# b_sub, -# ) where {T <: Real} -# model = solver.model - -# if !(Ap_fact isa QRPivoted{T, Matrix{T}}) -# row_piv = Ap_fact.prow -# model.c = model.c[row_piv] -# model.G = model.G[:, row_piv] -# solver.reduce_row_piv_inv = Ap_fact.rpivinv -# else -# solver.reduce_row_piv_inv = Int[] -# end - -# Q1_idxs = 1:Ap_rank -# Q2_idxs = (Ap_rank + 1):model.n - -# # [cQ1 cQ2] = c0' * Q -# Ap_Q = Ap_fact.Q -# cQ = model.c' * Ap_Q -# cQ1 = solver.reduce_cQ1 = cQ[Q1_idxs] -# cQ2 = cQ[Q2_idxs] -# # c = cQ2 -# model.c = cQ2 -# model.n = length(cQ2) -# # offset = offset0 + cQ1 * (R' \ b0) -# Rpib0 = solver.reduce_Rpib0 = ldiv!(Ap_R', b_sub) -# # solver.Rpib0 = Rpib0 # TODO -# model.obj_offset += dot(cQ1, Rpib0) - -# # [GQ1 GQ2] = G0 * Q -# # very inefficient method used for sparse G * QRSparseQ -# # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818 -# if model.G isa UniformScaling -# side = size(Ap_Q, 1) -# G_mul = Matrix{T}(model.G, side, side) -# elseif !isa(model.G, Matrix) -# G_mul = Matrix{T}(model.G) -# else -# G_mul = model.G -# end -# GQ = G_mul * Ap_Q -# GQ1 = solver.reduce_GQ1 = GQ[:, Q1_idxs] -# GQ2 = GQ[:, Q2_idxs] -# # h = h0 - GQ1 * (R' \ b0) -# mul!(model.h, GQ1, Rpib0, -1, true) - -# # G = GQ2 -# model.G = GQ2 - -# # A and b empty -# model.p = 0 -# model.A = zeros(T, 0, model.n) -# model.b = zeros(T, 0) -# solver.reduce_Ap_R = Ap_R -# solver.reduce_Ap_Q = Ap_Q - -# solver.reduce_y_keep_idxs = y_keep_idxs -# solver.Ap_R = UpperTriangular(zeros(T, 0, 0)) -# solver.Ap_Q = I - -# return zeros(T, 0) -# end - # (pivoted) QR factorizations are usually rank-revealing but may be unreliable # see http://www.math.sjsu.edu/~foster/rankrevealingcode.html # TODO could replace this with rank(qr_fact) when available for both dense and sparse diff --git a/test/runnativetests.jl b/test/runnativetests.jl index 9879fbe8b..301af8367 100644 --- a/test/runnativetests.jl +++ b/test/runnativetests.jl @@ -23,7 +23,7 @@ all_reals = [ ] diff_reals = [ Float64, - # BigFloat, + BigFloat, ] string_nameof(T) = string(nameof(T)) @@ -46,60 +46,58 @@ end @testset "native tests" begin -# @testset "default options tests" begin -# println("starting default options tests") -# inst_defaults = vcat( -# inst_preproc, -# # inst_infeas, -# # inst_cones_few, -# # # inst_cones_many, -# ) -# for inst_name in inst_defaults -# test_instance_solver(inst_name, Float64, default_options) -# end -# end +@testset "default options tests" begin + println("starting default options tests") + inst_defaults = vcat( + inst_preproc, + inst_infeas, + # inst_cones_few, + inst_cones_many, + ) + for inst_name in inst_defaults + test_instance_solver(inst_name, Float64, default_options) + end +end @testset "no preprocess tests" begin println("\nstarting no preprocess tests") for inst_name in inst_cones_few, T in diff_reals - # options = (; default_options..., preprocess = false, reduce = false, - # syssolver = Solvers.SymIndefDenseSystemSolver{T}()) - options = (; default_options..., preprocess = true, reduce = false, - syssolver = Solvers.SymIndefDenseSystemSolver{T}()) + options = (; default_options..., preprocess = false, reduce = false, + syssolver = Solvers.SymIndefDenseSystemSolver{T}()) test_instance_solver(inst_name, T, options) end end -# @testset "indirect solvers tests" begin -# println("\nstarting indirect solvers tests") -# for inst_name in inst_indirect, T in diff_reals -# options = (; default_options..., init_use_indirect = true, -# preprocess = false, reduce = false, -# syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), -# tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, -# tol_infeas = 1e-6) -# test_instance_solver(inst_name, T, options) -# end -# end +@testset "indirect solvers tests" begin + println("\nstarting indirect solvers tests") + for inst_name in inst_indirect, T in diff_reals + options = (; default_options..., init_use_indirect = true, + preprocess = false, reduce = false, + syssolver = Solvers.SymIndefIndirectSystemSolver{T}(), + tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, + tol_infeas = 1e-6) + test_instance_solver(inst_name, T, options) + end +end -# @testset "system solvers tests" begin -# println("\nstarting system solvers tests") -# syssolvers = [ -# (Solvers.NaiveDenseSystemSolver, diff_reals), -# (Solvers.NaiveSparseSystemSolver, [Float64,]), -# (Solvers.NaiveElimDenseSystemSolver, diff_reals), -# (Solvers.NaiveElimSparseSystemSolver, [Float64,]), -# (Solvers.SymIndefDenseSystemSolver, all_reals), -# (Solvers.SymIndefSparseSystemSolver, [Float64,]), -# (Solvers.QRCholDenseSystemSolver, all_reals), -# ] -# for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, -# T in real_types -# options = (; default_options..., syssolver = syssolver{T}(), -# reduce = false) -# test_instance_solver(inst_name, T, options, string_nameof(syssolver)) -# end -# end +@testset "system solvers tests" begin + println("\nstarting system solvers tests") + syssolvers = [ + (Solvers.NaiveDenseSystemSolver, diff_reals), + (Solvers.NaiveSparseSystemSolver, [Float64,]), + (Solvers.NaiveElimDenseSystemSolver, diff_reals), + (Solvers.NaiveElimSparseSystemSolver, [Float64,]), + (Solvers.SymIndefDenseSystemSolver, all_reals), + (Solvers.SymIndefSparseSystemSolver, [Float64,]), + (Solvers.QRCholDenseSystemSolver, all_reals), + ] + for inst_name in inst_minimal, (syssolver, real_types) in syssolvers, + T in real_types + options = (; default_options..., syssolver = syssolver{T}(), + reduce = false) + test_instance_solver(inst_name, T, options, string_nameof(syssolver)) + end +end # @testset "PredOrCentStepper tests" begin # verbose = true From a4d28b5510959132e5972827e2ae9734e946c9b5 Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 12 Dec 2021 13:11:28 -0500 Subject: [PATCH 6/9] clean up after preproc refac --- src/Solvers/Solvers.jl | 60 ++++++++++++++++++++++-------------------- src/Solvers/process.jl | 19 +++++++------ test/runnativetests.jl | 38 ++++++++++++++------------ 3 files changed, 62 insertions(+), 55 deletions(-) diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 1b403f0c5..756ccc6a0 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -116,6 +116,7 @@ mutable struct Solver{T <: Real} Ap_R::UpperTriangular{T, <:AbstractMatrix{T}} Ap_Q::Union{UniformScaling, AbstractMatrix{T}} AG_fact::Factorization{T} + AG_rank::Int AG_R::UpperTriangular{T, <:AbstractMatrix{T}} reduce_cQ1 reduce_Rpib0 @@ -292,9 +293,6 @@ function modify_h(solver::Solver{T}, h_new::Vector{T}) where {T <: Real} return solver end - - - function solve(solver::Solver{T}) where {T <: Real} init_status = solver.status if !in(init_status, (Loaded, Modified)) @@ -310,6 +308,7 @@ function solve(solver::Solver{T}) where {T <: Real} end if !in(solver.status, (PrimalInconsistent, DualInconsistent)) + initialize_point(solver) setup_stepping(solver) solver.verbose && print_header(solver.stepper, solver) @@ -335,9 +334,6 @@ function solve(solver::Solver{T}) where {T <: Real} return end - - - function setup_loaded(solver::Solver{T}) where {T <: Real} orig_model = solver.orig_model if solver.use_dense_model @@ -348,6 +344,8 @@ function setup_loaded(solver::Solver{T}) where {T <: Real} # copy original model to solver.model, which may be modified model = solver.model = Models.Model{T}( orig_model.c, orig_model.A, orig_model.b, orig_model.G, orig_model.h, + # orig_model.c, orig_model.A, orig_model.b, orig_model.G, copy(orig_model.h), + # orig_model.c, copy(orig_model.A), orig_model.b, copy(orig_model.G), orig_model.h, orig_model.cones, obj_offset = orig_model.obj_offset) solver.time_rescale = @elapsed solver.used_rescaling = rescale_data(solver) @@ -371,25 +369,17 @@ function setup_modified(solver::Solver{T}) where {T <: Real} orig_model = solver.orig_model model = solver.model - copyto!(model.c, orig_model.c) - copyto!(model.b, orig_model.b) - copyto!(model.h, orig_model.h) + model.c = copy(orig_model.c) + model.b = copy(orig_model.b) + model.h = copy(orig_model.h) if solver.used_rescaling # rescale using scalings computed during first load model.c ./= solver.c_scale model.b ./= solver.b_scale model.h ./= solver.h_scale end - # TODO needs fixing for preproc/reduce - - - - # TODO don't redo factorize stuff etc in here: - setup_point(solver) - if solver.status != SolveCalled - return - end + update_initial_point(solver) # TODO easier to update these system solvers (just need an update to setup_point_sub): @assert solver.syssolver isa Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}} @@ -399,12 +389,10 @@ function setup_modified(solver::Solver{T}) where {T <: Real} return end +# preprocess and find initial point function setup_point(solver::Solver{T}) where {T <: Real} - # preprocess and find initial point (init_z, init_s) = initialize_cone_point(solver.orig_model) - - if solver.reduce solver.time_inity = @elapsed handle_primal_eq(solver) solver.status == PrimalInconsistent && return @@ -423,12 +411,31 @@ function setup_point(solver::Solver{T}) where {T <: Real} init_y = update_primal_eq(solver, init_z) end - model = solver.model - point = solver.point = Point(model) + point = solver.point = Point(solver.model) point.x .= init_x point.y .= init_y point.z .= init_z point.s .= init_s + return +end + +function update_initial_point(solver::Solver{T}) where {T <: Real} + (init_z, init_s) = initialize_cone_point(solver.orig_model) + point = solver.point + if solver.reduce + point.y .= update_primal_eq(solver, init_z) + point.x .= update_dual_eq(solver, init_s) + else + point.x .= update_dual_eq(solver, init_s) + point.y .= update_primal_eq(solver, init_z) + end + point.z .= init_z + point.s .= init_s + return +end + +function initialize_point(solver::Solver{T}) where {T <: Real} + point = solver.point point.tau[] = one(T) point.kap[] = one(T) @@ -437,16 +444,13 @@ function setup_point(solver::Solver{T}) where {T <: Real} @warn("initial mu is $(solver.mu) but should be 1 (this could " * "indicate a problem with cone barrier oracles)") end + + model = solver.model Cones.load_point.(model.cones, point.primal_views) Cones.load_dual_point.(model.cones, point.dual_views) return end - - - - - function setup_stepping(solver::Solver{T}) where {T <: Real} model = solver.model solver.x_conv_tol = inv(1 + norm(model.c, Inf)) diff --git a/src/Solvers/process.jl b/src/Solvers/process.jl index f585b04cd..506f03373 100644 --- a/src/Solvers/process.jl +++ b/src/Solvers/process.jl @@ -76,7 +76,9 @@ function handle_dual_eq(solver::Solver{T}) where {T <: Real} # direct method AG = if iszero(model.p) # A is empty - if issparse(G) + if G isa UniformScaling + sparse(G, n, n) + elseif issparse(G) G elseif G isa Matrix{T} copy(G) @@ -98,7 +100,7 @@ function handle_dual_eq(solver::Solver{T}) where {T <: Real} else qr!(AG, ColumnNorm()) end - AG_rank = get_rank_est(AG_fact, solver.init_tol_qr) + AG_rank = solver.AG_rank = get_rank_est(AG_fact, solver.init_tol_qr) if !solver.preprocess || (AG_rank == n) if AG_rank < n @@ -161,10 +163,11 @@ function update_dual_eq( model.c = model.c[solver.x_keep_idxs] iszero(model.n) && return zeros(T, 0) + rhs_x = vcat(model.b, model.h - init_s) + if solver.init_use_indirect || !isa(model.A, MatrixyAG) || !isa(model.G, MatrixyAG) # use indirect method TODO pick lsqr or lsmr - rhs_x = vcat(model.b, model.h - init_s) if iszero(model.p) AG = model.G else @@ -175,10 +178,8 @@ function update_dual_eq( return IterativeSolvers.lsqr(AG, rhs_x) end - rhs_x = vcat(model.b, model.h - init_s) AG_fact = solver.AG_fact - - if !solver.preprocess || (model.n == size(AG_fact, 2)) + if !solver.preprocess || (solver.AG_rank == size(AG_fact, 2)) return AG_fact \ rhs_x end @@ -307,10 +308,8 @@ function handle_primal_eq(solver::Solver{T}) where {T <: Real} end GQ = G_mul * Ap_Q solver.reduce_GQ1 = GQ[:, 1:Ap_rank] - GQ2 = GQ[:, (Ap_rank + 1):end] - # G = GQ2 - model.G = GQ2 + model.G = GQ[:, (Ap_rank + 1):end] # A and b empty model.p = 0 @@ -379,7 +378,7 @@ function update_primal_eq( GQ1 = solver.reduce_GQ1 # h = h0 - GQ1 * (R' \ b0) - mul!(model.h, GQ1, Rpib0, -1, true) + model.h = mul!(copy(model.h), GQ1, Rpib0, -1, true) return zeros(T, 0) end diff --git a/test/runnativetests.jl b/test/runnativetests.jl index 301af8367..18ffeb9fa 100644 --- a/test/runnativetests.jl +++ b/test/runnativetests.jl @@ -62,8 +62,8 @@ end @testset "no preprocess tests" begin println("\nstarting no preprocess tests") for inst_name in inst_cones_few, T in diff_reals - options = (; default_options..., preprocess = false, reduce = false, - syssolver = Solvers.SymIndefDenseSystemSolver{T}()) + options = (; default_options..., rescale = false, preprocess = false, + reduce = false, syssolver = Solvers.SymIndefDenseSystemSolver{T}()) test_instance_solver(inst_name, T, options) end end @@ -138,21 +138,25 @@ end # end # end -# @testset "model modification tests" begin -# println("\nstarting model modification tests") -# syssolvers = [ -# (Solvers.SymIndefDenseSystemSolver, diff_reals), -# # (Solvers.SymIndefSparseSystemSolver, [Float64,]), -# # (Solvers.QRCholDenseSystemSolver, diff_reals), -# ] -# for inst_name in inst_modify, (syssolver, real_types) in syssolvers, -# T in real_types, rescale in (false, true) -# # TODO later also rescale/preproc/reduce true -# options = (; default_options..., syssolver = syssolver{T}(), -# rescale = rescale, preprocess = false, reduce = false) -# test_instance_solver(inst_name, T, options) -# end -# end +@testset "model modification tests" begin + println("\nstarting model modification tests") + syssolvers = [ + (Solvers.SymIndefDenseSystemSolver, diff_reals), + (Solvers.SymIndefSparseSystemSolver, [Float64,]), + (Solvers.QRCholDenseSystemSolver, diff_reals), + ] + for inst_name in inst_modify, (syssolver, real_types) in syssolvers, + T in real_types, rescale in (false, true), preprocess in (false, true), + reduce in (false, true) + if !preprocess + reduce && continue + syssolver == Solvers.QRCholDenseSystemSolver && continue + end + options = (; default_options..., syssolver = syssolver{T}(), + rescale = rescale, preprocess = preprocess, reduce = reduce) + test_instance_solver(inst_name, T, options) + end +end end ; From d94e25d84f426785c1c66c77ac3ab9450c558c45 Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 12 Dec 2021 14:54:28 -0500 Subject: [PATCH 7/9] add MOI modification functions --- src/MathOptInterface/wrapper.jl | 46 ++++++++++++++++++++ src/Solvers/Solvers.jl | 62 +++++++++++++++++++------- src/Solvers/process.jl | 2 +- test/nativeinstances.jl | 7 ++- test/runmoitests.jl | 2 +- test/runnativetests.jl | 77 ++++++++++++++++----------------- test/runtests.jl | 4 +- 7 files changed, 139 insertions(+), 61 deletions(-) diff --git a/src/MathOptInterface/wrapper.jl b/src/MathOptInterface/wrapper.jl index 948d536e3..2a5bc6f1b 100644 --- a/src/MathOptInterface/wrapper.jl +++ b/src/MathOptInterface/wrapper.jl @@ -212,6 +212,52 @@ end MOI.optimize!(opt::Optimizer) = Solvers.solve(opt.solver) +function MOI.modify( + opt::Optimizer{T}, + ::MOI.ObjectiveFunction{SAF{T}}, + chg::MOI.ScalarConstantChange{T}, + ) where {T} + obj_offset = chg.new_constant + if opt.obj_sense == MOI.MAX_SENSE + obj_offset = -obj_offset + end + Solvers.modify_obj_offset(opt.solver, obj_offset) + return +end + +function MOI.modify( + opt::Optimizer{T}, + ::MOI.ObjectiveFunction{SAF{T}}, + chg::MOI.ScalarCoefficientChange{T}, + ) where {T} + new_c = chg.new_coefficient + if opt.obj_sense == MOI.MAX_SENSE + new_c = -new_c + end + Solvers.modify_c(opt.solver, [chg.variable.value], [new_c]) + return +end + +function MOI.modify( + opt::Optimizer{T}, + ci::MOI.ConstraintIndex{VAF{T}, MOI.Zeros}, + chg::MOI.VectorConstantChange{T}, + ) where {T} + idxs = opt.zeros_idxs[ci.value] + Solvers.modify_b(opt.solver, idxs, chg.new_constant) + return +end + +function MOI.modify( + opt::Optimizer{T}, + ci::MOI.ConstraintIndex{VAF{T}, <:SupportedCone{T}}, + chg::MOI.VectorConstantChange{T}, + ) where {T} + idxs = opt.moi_cone_idxs[ci.value] + Solvers.modify_h(opt.solver, idxs, chg.new_constant) + return +end + function MOI.get(opt::Optimizer, ::MOI.TerminationStatus) status = opt.solver.status if status in (Solvers.NotLoaded, Solvers.Loaded) diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 756ccc6a0..63d775e38 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -263,34 +263,67 @@ end function load(solver::Solver{T}, model::Models.Model{T}) where {T <: Real} solver.orig_model = model solver.status = Loaded - return solver + return end +function modify_obj_offset(solver::Solver{T}, offset_new::T) where {T <: Real} + solver.orig_model.obj_offset = offset_new + solver.status = Modified + return +end - - -# TODO in the MOI wrapper the b, c, h changes can be cached and called at the -# start of optimize if there is a modify flag. this matters most for h because -# it will be updated constraint by constraint -function modify_c(solver::Solver{T}, c_new::Vector{T}) where {T <: Real} +function modify_c(solver::Solver{T}, c_new::AbstractVector{T}) where {T <: Real} @assert length(c_new) == solver.orig_model.n solver.orig_model.c = c_new solver.status = Modified - return solver + return +end + +function modify_c( + solver::Solver{T}, + idxs::AbstractVector{Int}, + c_new::Vector{T}, + ) where {T <: Real} + @assert length(c_new) == length(idxs) <= solver.orig_model.n + solver.orig_model.c[idxs] .= c_new + solver.status = Modified + return end -function modify_b(solver::Solver{T}, b_new::Vector{T}) where {T <: Real} +function modify_b(solver::Solver{T}, b_new::AbstractVector{T}) where {T <: Real} @assert length(b_new) == solver.orig_model.p solver.orig_model.b = b_new solver.status = Modified - return solver + return end -function modify_h(solver::Solver{T}, h_new::Vector{T}) where {T <: Real} +function modify_b( + solver::Solver{T}, + idxs::AbstractVector{Int}, + b_new::Vector{T}, + ) where {T <: Real} + @assert length(b_new) == length(idxs) <= solver.orig_model.p + solver.orig_model.b[idxs] .= b_new + solver.status = Modified + return +end + +function modify_h(solver::Solver{T}, h_new::AbstractVector{T}) where {T <: Real} @assert length(h_new) == solver.orig_model.q solver.orig_model.h = h_new solver.status = Modified - return solver + return +end + +function modify_h( + solver::Solver{T}, + idxs::AbstractVector{Int}, + h_new::Vector{T}, + ) where {T <: Real} + @assert length(h_new) == length(idxs) <= solver.orig_model.q + solver.orig_model.h[idxs] .= h_new + solver.status = Modified + return end function solve(solver::Solver{T}) where {T <: Real} @@ -344,8 +377,6 @@ function setup_loaded(solver::Solver{T}) where {T <: Real} # copy original model to solver.model, which may be modified model = solver.model = Models.Model{T}( orig_model.c, orig_model.A, orig_model.b, orig_model.G, orig_model.h, - # orig_model.c, orig_model.A, orig_model.b, orig_model.G, copy(orig_model.h), - # orig_model.c, copy(orig_model.A), orig_model.b, copy(orig_model.G), orig_model.h, orig_model.cones, obj_offset = orig_model.obj_offset) solver.time_rescale = @elapsed solver.used_rescaling = rescale_data(solver) @@ -369,6 +400,7 @@ function setup_modified(solver::Solver{T}) where {T <: Real} orig_model = solver.orig_model model = solver.model + model.obj_offset = orig_model.obj_offset model.c = copy(orig_model.c) model.b = copy(orig_model.b) model.h = copy(orig_model.h) @@ -383,8 +415,6 @@ function setup_modified(solver::Solver{T}) where {T <: Real} # TODO easier to update these system solvers (just need an update to setup_point_sub): @assert solver.syssolver isa Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}} - # TODO cleaner interface: call a update_syssolver function that calls this: - # (and errors if syssolver does not support update) set_point_sub_rhs(solver.syssolver, solver.model) return end diff --git a/src/Solvers/process.jl b/src/Solvers/process.jl index 506f03373..f28baf8e2 100644 --- a/src/Solvers/process.jl +++ b/src/Solvers/process.jl @@ -372,7 +372,7 @@ function update_primal_eq( Rpib0 = solver.reduce_Rpib0 = model.b[solver.reduce_y_keep_idxs] ldiv!(solver.reduce_Ap_R', Rpib0) # offset = offset0 + cQ1 * (R' \ b0) - model.obj_offset = solver.orig_model.obj_offset + dot(cQ1, Rpib0) + model.obj_offset += dot(cQ1, Rpib0) model.b = zeros(T, 0) diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl index f628816c1..b35055d6e 100644 --- a/test/nativeinstances.jl +++ b/test/nativeinstances.jl @@ -2782,13 +2782,15 @@ function modify1(T; options...) c = T[2, 1] b = [T(2)] + obj_offset = T(1) Solvers.modify_c(solver, c) Solvers.modify_b(solver, b) + Solvers.modify_obj_offset(solver, obj_offset) @test Solvers.get_status(solver) == Solvers.Modified - r = build_solve_check(c, A, b, G, h, cones, tol; + r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = obj_offset, solver = solver, already_loaded = true, options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 2 atol=tol rtol=tol + @test r.primal_obj ≈ 3 atol=tol rtol=tol @test r.x ≈ [0, 2] atol=tol rtol=tol @test iszero(solver.time_rescale) @test iszero(solver.time_loadsys) @@ -2797,6 +2799,7 @@ function modify1(T; options...) h = T[1, -2] Solvers.modify_c(solver, c) Solvers.modify_h(solver, h) + Solvers.modify_obj_offset(solver, T(0)) @test Solvers.get_status(solver) == Solvers.Modified r = build_solve_check(c, A, b, G, h, cones, tol; solver = solver, already_loaded = true, options...) diff --git a/test/runmoitests.jl b/test/runmoitests.jl index 8c1f5c58d..77b26e21c 100644 --- a/test/runmoitests.jl +++ b/test/runmoitests.jl @@ -29,7 +29,7 @@ end test_T = [ (Float64, 2 * sqrt(sqrt(eps())), 4, String[]), # TODO add test_linear after MOI 0.10.7 is tagged: - (BigFloat, 2 * eps(BigFloat)^0.2, 1, String["test_conic"]), + (BigFloat, 2 * eps(BigFloat)^0.2, 1, String["test_conic", "test_modification"]), ] @testset "MOI.Test tests: $T" for (T, tol_test, tol_relax, includes) in test_T diff --git a/test/runnativetests.jl b/test/runnativetests.jl index 18ffeb9fa..36b18b71f 100644 --- a/test/runnativetests.jl +++ b/test/runnativetests.jl @@ -51,7 +51,6 @@ end inst_defaults = vcat( inst_preproc, inst_infeas, - # inst_cones_few, inst_cones_many, ) for inst_name in inst_defaults @@ -99,44 +98,44 @@ end end end -# @testset "PredOrCentStepper tests" begin -# verbose = true -# println("\nstarting PredOrCentStepper tests (with printing)") - -# # adjustment and curve search -# use_adj_curv = [(false, false), (true, false), (true, true)] -# for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals -# stepper = Solvers.PredOrCentStepper{T}(; -# use_adjustment = adj, use_curve_search = curv) -# options = (; default_options..., verbose = verbose, stepper = stepper) -# test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") -# end - -# # other options -# for inst_name in inst_minimal -# T = Float64 -# stepper = Solvers.PredOrCentStepper{T}(; -# # stepper options -# use_adjustment = false, use_curve_search = false, -# max_cent_steps = 8, pred_prox_bound = 0.0332, -# # searcher options -# min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, -# alpha_sched = [0.9999 * 0.7^i for i in 0:22]) -# options = (; default_options..., verbose = verbose, stepper = stepper) -# test_instance_solver(inst_name, T, options, "other") -# end -# end - -# @testset "CombinedStepper tests" begin -# verbose = true -# println("\nstarting CombinedStepper tests (with printing)") -# shifts = [0, 2] -# for inst_name in inst_minimal, shift in shifts, T in diff_reals -# options = (; default_options..., verbose = verbose, -# stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) -# test_instance_solver(inst_name, T, options, "shift=$shift") -# end -# end +@testset "PredOrCentStepper tests" begin + verbose = true + println("\nstarting PredOrCentStepper tests (with printing)") + + # adjustment and curve search + use_adj_curv = [(false, false), (true, false), (true, true)] + for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals + stepper = Solvers.PredOrCentStepper{T}(; + use_adjustment = adj, use_curve_search = curv) + options = (; default_options..., verbose = verbose, stepper = stepper) + test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv") + end + + # other options + for inst_name in inst_minimal + T = Float64 + stepper = Solvers.PredOrCentStepper{T}(; + # stepper options + use_adjustment = false, use_curve_search = false, + max_cent_steps = 8, pred_prox_bound = 0.0332, + # searcher options + min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false, + alpha_sched = [0.9999 * 0.7^i for i in 0:22]) + options = (; default_options..., verbose = verbose, stepper = stepper) + test_instance_solver(inst_name, T, options, "other") + end +end + +@testset "CombinedStepper tests" begin + verbose = true + println("\nstarting CombinedStepper tests (with printing)") + shifts = [0, 2] + for inst_name in inst_minimal, shift in shifts, T in diff_reals + options = (; default_options..., verbose = verbose, + stepper = Solvers.CombinedStepper{T}(shift_sched = shift)) + test_instance_solver(inst_name, T, options, "shift=$shift") + end +end @testset "model modification tests" begin println("\nstarting model modification tests") diff --git a/test/runtests.jl b/test/runtests.jl index 3f830c86f..85f995a71 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,8 +7,8 @@ using Printf using Hypatia test_files = [ - # "polyutils", - # "cone", + "polyutils", + "cone", "native", "moi", # "examples", # CI runs examples separately From 96f57a416d67ae8da5ab3f197ded605103c95206 Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 12 Dec 2021 15:24:01 -0500 Subject: [PATCH 8/9] fix bug in permute_affine --- examples/JuMP_utils.jl | 4 +-- src/MathOptInterface/cones.jl | 4 +-- test/nativeinstances.jl | 46 +++++++++++++++++++++++------------ test/nativesets.jl | 2 +- 4 files changed, 35 insertions(+), 21 deletions(-) diff --git a/examples/JuMP_utils.jl b/examples/JuMP_utils.jl index ab8cd21f0..f4cd7be70 100644 --- a/examples/JuMP_utils.jl +++ b/examples/JuMP_utils.jl @@ -184,8 +184,8 @@ function solve_check( Hypatia.rescale_affine(cone, z_k) end if Hypatia.needs_permute(cone) - Hypatia.permute_affine(cone, s_k) - Hypatia.permute_affine(cone, z_k) + s_k = Hypatia.permute_affine(cone, s_k) + z_k = Hypatia.permute_affine(cone, z_k) end append!(s, s_k) append!(z, z_k) diff --git a/src/MathOptInterface/cones.jl b/src/MathOptInterface/cones.jl index 6ee7293df..902202f99 100644 --- a/src/MathOptInterface/cones.jl +++ b/src/MathOptInterface/cones.jl @@ -83,8 +83,8 @@ end needs_permute(cone::SpecNucCone) = needs_untransform(cone) function permute_affine(cone::SpecNucCone, vals::AbstractVector{T}) where {T} - @views vals[2:end] = reshape(vals[2:end], cone.row_dim, cone.column_dim)' - return vals + w_vals = reshape(vals[2:end], cone.row_dim, cone.column_dim)' + return vcat(vals[1], vec(w_vals)) end function permute_affine(cone::SpecNucCone, func::VAF{T}) where {T} diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl index b35055d6e..7a3326d76 100644 --- a/test/nativeinstances.jl +++ b/test/nativeinstances.jl @@ -2797,29 +2797,43 @@ function modify1(T; options...) c = T[3, 1] h = T[1, -2] + obj_offset = T(-1) Solvers.modify_c(solver, c) Solvers.modify_h(solver, h) - Solvers.modify_obj_offset(solver, T(0)) + Solvers.modify_obj_offset(solver, obj_offset) @test Solvers.get_status(solver) == Solvers.Modified - r = build_solve_check(c, A, b, G, h, cones, tol; + r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = obj_offset, solver = solver, already_loaded = true, options...) @test r.status == Solvers.Optimal - @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test r.primal_obj ≈ -1 atol=tol rtol=tol @test r.x ≈ [-1, 3] atol=tol rtol=tol @test iszero(solver.time_rescale) @test iszero(solver.time_loadsys) end -# function modify2(T; options...) -# tol = test_tol(T) -# c = T[-1, -1, 0] -# A = zeros(T, 0, 3) -# b = T[] -# G = repeat(SparseMatrixCSC(-one(T) * I, 3, 3), 2, 1) -# h = zeros(T, 6) -# cones = Cone{T}[Cones.EpiNormInf{T, T}(3), -# Cones.EpiNormInf{T, T}(3, use_dual = true)] -# error() -# r = build_solve_check(c, A, b, G, h, cones, tol; options...) -# @test r.status == Solvers.DualInfeasible -# end +function modify2(T; options...) + tol = test_tol(T) + c = T[-1, -1, 0] + A = zeros(T, 0, 3) + b = T[] + G = repeat(SparseMatrixCSC(-one(T) * I, 3, 3), 2, 1) + h = zeros(T, 6) + cones = Cone{T}[Cones.EpiNormInf{T, T}(3), + Cones.EpiNormInf{T, T}(3, use_dual = true)] + + r = build_solve_check(c, A, b, G, h, cones, tol; options...) + @test r.status == Solvers.DualInfeasible + solver = r.solver + + c = zeros(T, 3) + Solvers.modify_c(solver, c) + @test Solvers.get_status(solver) == Solvers.Modified + r = build_solve_check(c, A, b, G, h, cones, tol; + solver = solver, already_loaded = true, options...) + @test r.status == Solvers.Optimal + @test r.primal_obj ≈ 0 atol=tol rtol=tol + @test r.dual_obj ≈ 0 atol=tol rtol=tol + @test norm(r.z) ≈ 0 atol=tol rtol=tol + @test iszero(solver.time_rescale) + @test iszero(solver.time_loadsys) +end \ No newline at end of file diff --git a/test/nativesets.jl b/test/nativesets.jl index b9c4e46f8..9de6c5a6e 100644 --- a/test/nativesets.jl +++ b/test/nativesets.jl @@ -180,5 +180,5 @@ inst_indirect = [ inst_modify = [ "modify1", - # "modify2", + "modify2", ] From 9ea5ad68185c9d497249c00d4559e14a03b30a8e Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 12 Dec 2021 22:14:40 -0500 Subject: [PATCH 9/9] fix bug in modify h thru MOI; adjust illposed check --- src/MathOptInterface/wrapper.jl | 14 ++++++++++++-- src/Solvers/Solvers.jl | 5 ++--- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/MathOptInterface/wrapper.jl b/src/MathOptInterface/wrapper.jl index 2a5bc6f1b..0e8cfd8b8 100644 --- a/src/MathOptInterface/wrapper.jl +++ b/src/MathOptInterface/wrapper.jl @@ -253,8 +253,18 @@ function MOI.modify( ci::MOI.ConstraintIndex{VAF{T}, <:SupportedCone{T}}, chg::MOI.VectorConstantChange{T}, ) where {T} - idxs = opt.moi_cone_idxs[ci.value] - Solvers.modify_h(opt.solver, idxs, chg.new_constant) + i = ci.value + idxs = opt.moi_cone_idxs[i] + set = opt.moi_cones[i] + new_h = chg.new_constant + if needs_permute(set) + @assert !needs_rescale(set) + new_h = permute_affine(set, new_h) + end + if needs_rescale(set) + rescale_affine(set, new_h) + end + Solvers.modify_h(opt.solver, idxs, new_h) return end diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl index 63d775e38..7a81c93ff 100644 --- a/src/Solvers/Solvers.jl +++ b/src/Solvers/Solvers.jl @@ -229,7 +229,7 @@ mutable struct Solver{T <: Real} tol_infeas = default_tol_tight end if isnothing(tol_illposed) - tol_illposed = default_tol_tight / 100 + tol_illposed = default_tol_tight end @assert min(tol_rel_opt, tol_abs_opt, tol_feas, tol_infeas, tol_illposed, tol_slow) >= 0 @@ -702,8 +702,7 @@ function check_converged( end # TODO experiment with ill-posedness check - max_illp = max(solver.mu, tau / min(one(T), solver.point.kap[])) - if max_illp <= near_factor * solver.tol_illposed + if max(tau, solver.point.kap[]) <= near_factor * solver.tol_illposed solver.verbose && println("ill-posedness detected; terminating") solver.status = (check_near ? NearIllPosed : IllPosed) return true