From 39e1424a4370313500028a85d6b0a97285f46861 Mon Sep 17 00:00:00 2001 From: fpacaud Date: Tue, 23 Jul 2024 09:04:06 +0200 Subject: [PATCH 1/5] Revisit fixed variables treatment * Change the behavior of MakeParameter, depending on the callback: - `SparseCallback`: remove the fixed variables and reduce problem's dimension in MadNLP - `DenseCallback`: keep the fixed variables frozen (as done previously) * Add a new `fixed_handler` if the problem has no fixed variable: `NoFixedVariable`. * Ensure we always interact to the nonlinear program `nlp` through the callback wrapper `cb` to avoid side effect. * Remove the function `get_index_constraint` and stores indexes directly in the callback wrapper `cb`. * Fix the computation of the bound's multipliers for fixed variables. --- ext/MadNLPMOI/MadNLPMOI.jl | 2 +- lib/MadNLPTests/src/wrapper.jl | 52 +- src/IPM/IPM.jl | 42 +- src/IPM/kernels.jl | 13 - src/IPM/solver.jl | 8 +- src/IPM/utils.jl | 52 +- src/KKT/Dense/augmented.jl | 13 +- src/KKT/Dense/condensed.jl | 15 +- src/KKT/Sparse/augmented.jl | 12 +- src/KKT/Sparse/condensed.jl | 9 +- src/KKT/Sparse/unreduced.jl | 11 +- src/nlpmodels.jl | 1075 ++++++++++++++++++++------------ src/options.jl | 14 +- test/MOI_interface_test.jl | 5 - test/kkt_test.jl | 5 - test/madnlp_dense.jl | 21 +- test/madnlp_test.jl | 27 +- 17 files changed, 850 insertions(+), 526 deletions(-) diff --git a/ext/MadNLPMOI/MadNLPMOI.jl b/ext/MadNLPMOI/MadNLPMOI.jl index aada2d4e..f40a0c77 100644 --- a/ext/MadNLPMOI/MadNLPMOI.jl +++ b/ext/MadNLPMOI/MadNLPMOI.jl @@ -1142,7 +1142,7 @@ function MOI.get( MOI.check_result_index_bounds(model, attr) MOI.throw_if_not_valid(model, ci) rc = model.result.multipliers_L[ci.value] - model.result.multipliers_U[ci.value] - return rc + return -rc end function MOI.get( diff --git a/lib/MadNLPTests/src/wrapper.jl b/lib/MadNLPTests/src/wrapper.jl index 86098c05..8e478e33 100644 --- a/lib/MadNLPTests/src/wrapper.jl +++ b/lib/MadNLPTests/src/wrapper.jl @@ -2,38 +2,38 @@ abstract type AbstractWrapperModel{T,VT} <: NLPModels.AbstractNLPModel{T,VT} end struct DenseWrapperModel{T,VT,T2,VT2,MT2, I <: NLPModels.AbstractNLPModel{T2,VT2}} <: AbstractWrapperModel{T,VT} inner::I - + x::VT2 y::VT2 - + con::VT2 grad::VT2 jac::MT2 hess::MT2 - + meta::NLPModels.NLPModelMeta{T, VT} - counters::NLPModels.Counters + counters::NLPModels.Counters end struct SparseWrapperModel{T,VT,T2,VI2,VT2,I <: NLPModels.AbstractNLPModel{T2,VT2}} <: AbstractWrapperModel{T,VT} inner::I - + jrows::VI2 jcols::VI2 hrows::VI2 hcols::VI2 - + x::VT2 y::VT2 - + con::VT2 grad::VT2 jac::VT2 hess::VT2 - + meta::NLPModels.NLPModelMeta{T, VT} - counters::NLPModels.Counters + counters::NLPModels.Counters end @@ -121,10 +121,10 @@ function NLPModels.cons!( g::V ) where {M <: AbstractWrapperModel, V <: AbstractVector} - copyto!(m.x, x) + copyto!(m.x, x) NLPModels.cons!(m.inner, m.x, m.con) copyto!(g, m.con) - return + return end function NLPModels.grad!( m::M, @@ -143,7 +143,7 @@ function NLPModels.jac_structure!( rows::V, cols::V ) where {M <: SparseWrapperModel, V <: AbstractVector} - + NLPModels.jac_structure!(m.inner, m.jrows, m.jcols) copyto!(rows, m.jrows) copyto!(cols, m.jcols) @@ -159,13 +159,26 @@ function NLPModels.hess_structure!( copyto!(rows, m.hrows) copyto!(cols, m.hcols) end +function NLPModels.jtprod!( + m::SparseWrapperModel, + x::AbstractVector, + v::AbstractVector, + jtv::AbstractVector, +) + copyto!(m.x, x) + copyto!(m.grad, jtv) + copyto!(m.con, v) + NLPModels.jtprod!(m.inner, m.x, m.con, m.grad) + copyto!(jtv, m.grad) + return +end function NLPModels.jac_coord!( - m::M, - x::V, - jac::V - ) where {M <: SparseWrapperModel, V <: AbstractVector} + m::SparseWrapperModel, + x::AbstractVector, + jac::AbstractVector, +) - copyto!(m.x, x) + copyto!(m.x, x) NLPModels.jac_coord!(m.inner, m.x, m.jac) copyto!(jac, m.jac) return @@ -185,19 +198,18 @@ function NLPModels.hess_coord!( return end - - function MadNLP.jac_dense!( m::Model, x::V, jac::M ) where {Model <: DenseWrapperModel, V <: AbstractVector, M <: AbstractMatrix} - copyto!(m.x, x) + copyto!(m.x, x) MadNLP.jac_dense!(m.inner, m.x, m.jac) copyto!(jac, m.jac) return end + function MadNLP.hess_dense!( m::Model, x::V, diff --git a/src/IPM/IPM.jl b/src/IPM/IPM.jl index 65fa4e74..ac57d2be 100644 --- a/src/IPM/IPM.jl +++ b/src/IPM/IPM.jl @@ -134,20 +134,13 @@ function MadNLPSolver(nlp::AbstractNLPModel{T,VT}; kwargs...) where {T, VT} set_blas_num_threads(ipm_opt.blas_num_threads; permanent=true) @trace(logger,"Initializing variables.") - ind_cons = get_index_constraints( - get_lvar(nlp), get_uvar(nlp), - get_lcon(nlp), get_ucon(nlp); - fixed_variable_treatment=ipm_opt.fixed_variable_treatment, - equality_treatment=ipm_opt.equality_treatment - ) - - ind_lb = ind_cons.ind_lb - ind_ub = ind_cons.ind_ub + ind_lb = cb.ind_lb + ind_ub = cb.ind_ub - ns = length(ind_cons.ind_ineq) - nx = get_nvar(nlp) + ns = length(cb.ind_ineq) + nx = n_variables(cb) n = nx+ns - m = get_ncon(nlp) + m = n_constraints(cb) nlb = length(ind_lb) nub = length(ind_ub) @@ -155,7 +148,6 @@ function MadNLPSolver(nlp::AbstractNLPModel{T,VT}; kwargs...) where {T, VT} kkt = create_kkt_system( ipm_opt.kkt_system, cb, - ind_cons, ipm_opt.linear_solver; hessian_approximation=ipm_opt.hessian_approximation, opt_linear_solver=options.linear_solver, @@ -185,17 +177,17 @@ function MadNLPSolver(nlp::AbstractNLPModel{T,VT}; kwargs...) where {T, VT} c = VT(undef, m) rhs = VT(undef, m) - c_slk = view(c,ind_cons.ind_ineq) - x_lr = view(full(x), ind_cons.ind_lb) - x_ur = view(full(x), ind_cons.ind_ub) - xl_r = view(full(xl), ind_cons.ind_lb) - xu_r = view(full(xu), ind_cons.ind_ub) - zl_r = view(full(zl), ind_cons.ind_lb) - zu_r = view(full(zu), ind_cons.ind_ub) - x_trial_lr = view(full(x_trial), ind_cons.ind_lb) - x_trial_ur = view(full(x_trial), ind_cons.ind_ub) - dx_lr = view(d.xp, ind_cons.ind_lb) # TODO - dx_ur = view(d.xp, ind_cons.ind_ub) # TODO + c_slk = view(c,cb.ind_ineq) + x_lr = view(full(x), cb.ind_lb) + x_ur = view(full(x), cb.ind_ub) + xl_r = view(full(xl), cb.ind_lb) + xu_r = view(full(xu), cb.ind_ub) + zl_r = view(full(zl), cb.ind_lb) + zu_r = view(full(zu), cb.ind_ub) + x_trial_lr = view(full(x_trial), cb.ind_lb) + x_trial_ur = view(full(x_trial), cb.ind_ub) + dx_lr = view(d.xp, cb.ind_lb) # TODO + dx_ur = view(d.xp, cb.ind_ub) # TODO inertia_correction_method = if ipm_opt.inertia_correction_method == InertiaAuto is_inertia(kkt.linear_solver)::Bool ? InertiaBased : InertiaFree @@ -221,7 +213,7 @@ function MadNLPSolver(nlp::AbstractNLPModel{T,VT}; kwargs...) where {T, VT} d, p, _w1, _w2, _w3, _w4, x_trial, c_trial, zero(T), c_slk, rhs, - ind_cons.ind_ineq, ind_cons.ind_fixed, ind_cons.ind_llb, ind_cons.ind_uub, + cb.ind_ineq, cb.ind_fixed, cb.ind_llb, cb.ind_uub, x_lr, x_ur, xl_r, xu_r, zl_r, zu_r, dx_lr, dx_ur, x_trial_lr, x_trial_ur, iterator, zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), diff --git a/src/IPM/kernels.jl b/src/IPM/kernels.jl index a89c82b2..56ebe4ab 100644 --- a/src/IPM/kernels.jl +++ b/src/IPM/kernels.jl @@ -814,19 +814,6 @@ function get_ftype(filter,theta,theta_trial,varphi,varphi_trial,switching_condit return " " end -# fixed variable treatment ---------------------------------------------------- -function _get_fixed_variable_index( - mat::SparseMatrixCSC{Tv,Ti1}, ind_fixed::Vector{Ti2} -) where {Tv,Ti1,Ti2} - fixed_aug_index = Int[] - for i in ind_fixed - append!(fixed_aug_index,append!(collect(mat.colptr[i]+1:mat.colptr[i+1]-1))) - end - append!(fixed_aug_index,setdiff!(Base._findin(mat.rowval,ind_fixed),mat.colptr)) - - return fixed_aug_index -end - function dual_inf_perturbation!(px, ind_llb, ind_uub, mu, kappa_d) px[ind_llb] .-= mu*kappa_d px[ind_uub] .+= mu*kappa_d diff --git a/src/IPM/solver.jl b/src/IPM/solver.jl index 4cd65ad3..bf5e0010 100644 --- a/src/IPM/solver.jl +++ b/src/IPM/solver.jl @@ -152,6 +152,12 @@ function solve!( set_options!(solver.opt, kwargs) end + # If the problem has no free variable, do nothing + if solver.n == 0 + update!(stats, solver) + return stats + end + try if solver.status == INITIAL @notice(solver.logger,"This is $(introduce()), running with $(introduce(solver.kkt.linear_solver))\n") @@ -216,7 +222,7 @@ function regular!(solver::AbstractMadNLPSolver{T}) where T if (solver.cnt.k!=0 && !solver.opt.jacobian_constant) eval_jac_wrapper!(solver, solver.kkt, solver.x) end - + jtprod!(solver.jacl, solver.kkt, solver.y) sd = get_sd(solver.y,solver.zl_r,solver.zu_r,T(solver.opt.s_max)) sc = get_sc(solver.zl_r,solver.zu_r,T(solver.opt.s_max)) diff --git a/src/IPM/utils.jl b/src/IPM/utils.jl index 834a67af..27399fae 100644 --- a/src/IPM/utils.jl +++ b/src/IPM/utils.jl @@ -20,27 +20,38 @@ mutable struct MadNLPExecutionStats{T, VT} <: AbstractExecutionStats counters::MadNLPCounters end -MadNLPExecutionStats(solver::MadNLPSolver) =MadNLPExecutionStats( - solver.opt, - solver.status, - primal(solver.x)[1:get_nvar(solver.nlp)], - solver.obj_val / solver.cb.obj_scale[], - solver.c ./ solver.cb.con_scale, - solver.inf_du, - solver.inf_pr, - copy(solver.y), - primal(solver.zl)[1:get_nvar(solver.nlp)], - primal(solver.zu)[1:get_nvar(solver.nlp)], - 0, - solver.cnt, -) +function MadNLPExecutionStats(solver::MadNLPSolver{T, VT}) where {T, VT} + n, m = get_nvar(solver.nlp), get_ncon(solver.nlp) + x = similar(VT, n) + zl = similar(VT, n) + zu = similar(VT, n) + c = similar(VT, m) + unpack_cons!(c, solver.cb, solver.c) + unpack_x!(x, solver.cb, variable(solver.x)) + unpack_z!(zl, solver.cb, variable(solver.zl)) + unpack_z!(zu, solver.cb, variable(solver.zu)) + return MadNLPExecutionStats( + solver.opt, + solver.status, + x, + unpack_obj(solver.cb, solver.obj_val), + c, + solver.inf_du, + solver.inf_pr, + copy(solver.y), + zl, + zu, + 0, + solver.cnt, + ) +end function update!(stats::MadNLPExecutionStats, solver::MadNLPSolver) stats.status = solver.status - stats.solution .= @view(primal(solver.x)[1:get_nvar(solver.nlp)]) + unpack_x!(stats.solution, solver.cb, variable(solver.x)) + unpack_z!(stats.multipliers_L, solver.cb, variable(solver.zl)) + unpack_z!(stats.multipliers_U, solver.cb, variable(solver.zu)) stats.multipliers .= solver.y - stats.multipliers_L .= @view(primal(solver.zl)[1:get_nvar(solver.nlp)]) - stats.multipliers_U .= @view(primal(solver.zu)[1:get_nvar(solver.nlp)]) # stats.solution .= min.( # max.( # @view(primal(solver.x)[1:get_nvar(solver.nlp)]), @@ -48,12 +59,13 @@ function update!(stats::MadNLPExecutionStats, solver::MadNLPSolver) # ), # get_uvar(solver.nlp) # ) - stats.objective = solver.obj_val / solver.cb.obj_scale[] - stats.constraints .= solver.c ./ solver.cb.con_scale .+ solver.rhs + stats.objective = unpack_obj(solver.cb, solver.obj_val) + unpack_cons!(stats.constraints, solver.cb, solver.c) + stats.constraints .+= solver.rhs stats.constraints[solver.ind_ineq] .+= slack(solver.x) stats.dual_feas = solver.inf_du stats.primal_feas = solver.inf_pr - update_z!(solver.cb, stats.multipliers_L, stats.multipliers_U, solver.jacl) + update_z!(solver.cb, stats.solution, stats.multipliers, stats.multipliers_L, stats.multipliers_U, solver.jacl) stats.iter = solver.cnt.k return stats end diff --git a/src/KKT/Dense/augmented.jl b/src/KKT/Dense/augmented.jl index 5ecee833..cca5b7d1 100644 --- a/src/KKT/Dense/augmented.jl +++ b/src/KKT/Dense/augmented.jl @@ -42,21 +42,20 @@ end function create_kkt_system( ::Type{DenseKKTSystem}, cb::AbstractCallback{T,VT}, - ind_cons, linear_solver::Type; opt_linear_solver=default_options(linear_solver), hessian_approximation=ExactHessian, ) where {T, VT} - ind_ineq = ind_cons.ind_ineq - ind_lb = ind_cons.ind_lb - ind_ub = ind_cons.ind_ub + ind_ineq = cb.ind_ineq + ind_lb = cb.ind_lb + ind_ub = cb.ind_ub n = cb.nvar m = cb.ncon ns = length(ind_ineq) - nlb = length(ind_cons.ind_lb) - nub = length(ind_cons.ind_ub) + nlb = length(cb.ind_lb) + nub = length(cb.ind_ub) hess = create_array(cb, n, n) jac = create_array(cb, m, n) @@ -87,7 +86,7 @@ function create_kkt_system( hess, jac, quasi_newton, reg, pr_diag, du_diag, l_diag, u_diag, l_lower, u_lower, diag_hess, aug_com, - ind_ineq, ind_cons.ind_lb, ind_cons.ind_ub, + ind_ineq, cb.ind_lb, cb.ind_ub, _linear_solver, Dict{Symbol, Any}(), ) diff --git a/src/KKT/Dense/condensed.jl b/src/KKT/Dense/condensed.jl index f6ae1a17..6a3ab7a3 100644 --- a/src/KKT/Dense/condensed.jl +++ b/src/KKT/Dense/condensed.jl @@ -52,7 +52,6 @@ end function create_kkt_system( ::Type{DenseCondensedKKTSystem}, cb::AbstractCallback{T,VT}, - ind_cons, linear_solver::Type; opt_linear_solver=default_options(linear_solver), hessian_approximation=ExactHessian, @@ -60,10 +59,10 @@ function create_kkt_system( n = cb.nvar m = cb.ncon - ns = length(ind_cons.ind_ineq) + ns = length(cb.ind_ineq) n_eq = m - ns - nlb = length(ind_cons.ind_lb) - nub = length(ind_cons.ind_ub) + nlb = length(cb.ind_lb) + nub = length(cb.ind_ub) aug_com = create_array(cb, n+m-ns, n+m-ns) hess = create_array(cb, n, n) @@ -90,8 +89,8 @@ function create_kkt_system( fill!(du_diag, zero(T)) # Shift indexes to avoid additional allocation in views - ind_eq_shifted = ind_cons.ind_eq .+ n .+ ns - ind_ineq_shifted = ind_cons.ind_ineq .+ n .+ ns + ind_eq_shifted = cb.ind_eq .+ n .+ ns + ind_ineq_shifted = cb.ind_ineq .+ n .+ ns quasi_newton = create_quasi_newton(hessian_approximation, cb, n) _linear_solver = linear_solver(aug_com; opt = opt_linear_solver) @@ -101,9 +100,9 @@ function create_kkt_system( reg, pr_diag, du_diag, l_diag, u_diag, l_lower, u_lower, pd_buffer, diag_buffer, buffer, aug_com, - n_eq, ind_cons.ind_eq, ind_eq_shifted, + n_eq, cb.ind_eq, ind_eq_shifted, ns, - ind_cons.ind_ineq, ind_cons.ind_lb, ind_cons.ind_ub, + cb.ind_ineq, cb.ind_lb, cb.ind_ub, ind_ineq_shifted, _linear_solver, Dict{Symbol, Any}(), diff --git a/src/KKT/Sparse/augmented.jl b/src/KKT/Sparse/augmented.jl index a54d0464..ff67fa7e 100644 --- a/src/KKT/Sparse/augmented.jl +++ b/src/KKT/Sparse/augmented.jl @@ -40,13 +40,12 @@ end function create_kkt_system( ::Type{SparseKKTSystem}, cb::SparseCallback{T,VT}, - ind_cons, linear_solver::Type; opt_linear_solver=default_options(linear_solver), hessian_approximation=ExactHessian, ) where {T,VT} - n_slack = length(ind_cons.ind_ineq) + n_slack = length(cb.ind_ineq) # Deduce KKT size. n = cb.nvar @@ -59,19 +58,18 @@ function create_kkt_system( quasi_newton = create_quasi_newton(hessian_approximation, cb, n) hess_sparsity_I, hess_sparsity_J = build_hessian_structure(cb, hessian_approximation) - nlb = length(ind_cons.ind_lb) - nub = length(ind_cons.ind_ub) + nlb = length(cb.ind_lb) + nub = length(cb.ind_ub) force_lower_triangular!(hess_sparsity_I,hess_sparsity_J) - ind_ineq = ind_cons.ind_ineq + ind_ineq = cb.ind_ineq n_slack = length(ind_ineq) n_jac = length(jac_sparsity_I) n_hess = length(hess_sparsity_I) n_tot = n + n_slack - aug_vec_length = n_tot+m aug_mat_length = n_tot+m+n_hess+n_jac+n_slack @@ -136,7 +134,7 @@ function create_kkt_system( hess_raw, hess_com, hess_csc_map, jac_raw, jac_com, jac_csc_map, _linear_solver, - ind_ineq, ind_cons.ind_lb, ind_cons.ind_ub, + ind_ineq, cb.ind_lb, cb.ind_ub, ) end diff --git a/src/KKT/Sparse/condensed.jl b/src/KKT/Sparse/condensed.jl index b4697b5d..0ae11aef 100644 --- a/src/KKT/Sparse/condensed.jl +++ b/src/KKT/Sparse/condensed.jl @@ -55,12 +55,11 @@ end function create_kkt_system( ::Type{SparseCondensedKKTSystem}, cb::SparseCallback{T,VT}, - ind_cons, linear_solver::Type; opt_linear_solver=default_options(linear_solver), hessian_approximation=ExactHessian, ) where {T, VT} - ind_ineq = ind_cons.ind_ineq + ind_ineq = cb.ind_ineq n = cb.nvar m = cb.ncon n_slack = length(ind_ineq) @@ -82,8 +81,8 @@ function create_kkt_system( n_jac = length(jac_sparsity_I) n_hess = length(hess_sparsity_I) n_tot = n + n_slack - nlb = length(ind_cons.ind_lb) - nub = length(ind_cons.ind_ub) + nlb = length(cb.ind_lb) + nub = length(cb.ind_ub) reg = VT(undef, n_tot) @@ -126,7 +125,7 @@ function create_kkt_system( buffer, buffer2, aug_com, diag_buffer, dptr, hptr, jptr, _linear_solver, - ind_ineq, ind_cons.ind_lb, ind_cons.ind_ub, + ind_ineq, cb.ind_lb, cb.ind_ub, ext ) end diff --git a/src/KKT/Sparse/unreduced.jl b/src/KKT/Sparse/unreduced.jl index 8cd48ae1..97f49403 100644 --- a/src/KKT/Sparse/unreduced.jl +++ b/src/KKT/Sparse/unreduced.jl @@ -47,18 +47,17 @@ end function create_kkt_system( ::Type{SparseUnreducedKKTSystem}, cb::SparseCallback{T,VT}, - ind_cons, linear_solver::Type; opt_linear_solver=default_options(linear_solver), hessian_approximation=ExactHessian, ) where {T, VT} - ind_ineq = ind_cons.ind_ineq - ind_lb = ind_cons.ind_lb - ind_ub = ind_cons.ind_ub + ind_ineq = cb.ind_ineq + ind_lb = cb.ind_lb + ind_ub = cb.ind_ub n_slack = length(ind_ineq) - nlb = length(ind_cons.ind_lb) - nub = length(ind_cons.ind_ub) + nlb = length(cb.ind_lb) + nub = length(cb.ind_ub) # Deduce KKT size. n = cb.nvar m = cb.ncon diff --git a/src/nlpmodels.jl b/src/nlpmodels.jl index 8336995a..16c21dde 100644 --- a/src/nlpmodels.jl +++ b/src/nlpmodels.jl @@ -1,3 +1,37 @@ +#= + MadNLP wrappers + + MadNLP adapts any AbstractNLPModel to avoid numerical issues. + The `AbstractNLPModel` is wrapped as a `SparseCallback` or + as a `DenseCallback` (if the dense callbacks `jac_dense!` and `hess_dense!` are specified). + + The wrapper reformulates the model by: + 1. scaling the objective and the constraints. + 2. removing the fixed variables from the formulation. + + The scaling can be switched off by setting `nlp_scaling=false` in the options. + + Four cases can occur for the fixed variables: + + Case #1. The problem doesn't have any fixed variable + MadNLP sets `fixed_variable_treatment=NoFixedVariables()` + and fallbacks to the default callbacks (that just apply the scaling). + + Case #2. `fixed_variable_treatment=RelaxBound` + MadNLP relax slightly the bounds of the fixed variables during + the initialization. The wrapper fallbacks to the default callbacks, similarly as with Case #1. + + Case #3. `callback=DenseCallback` and `fixed_variable_treatment=MakeParameter` + Reformulate the fixed variables as dummy variables that are kept + at their bounds. The problem's dimension is not modified. The wrapper + modifies the Jacobian by filling the fixed columns with 0, + and the Hessian by filling the fixed columns and rows with 0 for + the non-diagonal elements, 1 for the diagonal ones. + + Case #4. `callback=SparseCallback` and `fixed_variable_treatment=MakeParameter` + Remove the fixed variables from the model. As a consequence, the problem's + dimension is reduced after the reformulation. +=# """ AbstractFixedVariableTreatment @@ -6,17 +40,28 @@ Abstract type to define the reformulation of the fixed variables inside MadNLP. """ abstract type AbstractFixedVariableTreatment end +""" + NoFixedVariables <: AbstractFixedVariableTreatment + +Do nothing if the problem has no fixed variables. +""" +struct NoFixedVariables <: AbstractFixedVariableTreatment end + """ MakeParameter{VT, VI} <: AbstractFixedVariableTreatment Remove the fixed variables from the optimization variables and define them as problem's parameters. + """ -struct MakeParameter{VT,VI} <: AbstractFixedVariableTreatment +struct MakeParameter{T, VT,VI} <: AbstractFixedVariableTreatment + free::VI fixed::VI - fixedj::VI - fixedh::VI - grad_storage::VT + ind_jac_free::VI + ind_hess_free::VI + hash_x::Ref{T} + x_full::VT + g_full::VT end """ @@ -25,10 +70,10 @@ end Relax the fixed variables ``x = x_{fixed}`` as bounded variables ``x_{fixed} - ϵ ≤ x ≤ x_{fixed} + ϵ``, with ``ϵ`` a small-enough parameter. + """ struct RelaxBound <: AbstractFixedVariableTreatment end - """ AbstractEqualityTreatment @@ -59,82 +104,6 @@ constraints only up to a tolerance ``ϵ``. """ struct RelaxEquality <: AbstractEqualityTreatment end - -""" - get_index_constraints(nlp::AbstractNLPModel) - -Analyze the bounds of the variables and the constraints in the `AbstractNLPModel` `nlp`. -Return a named-tuple witht the following keys:return ( - -* `ind_eq`: indices of equality constraints. -* `ind_ineq`: indices of inequality constraints. -* `ind_fixed`: indices of fixed variables. -* `ind_lb`: indices of variables with a lower-bound. -* `ind_ub`: indices of variables with an upper-bound. -* `ind_llb`: indices of variables with *only* a lower-bound. -* `ind_uub`: indices of variables with *only* an upper-bound. - -""" -function get_index_constraints( - nlp::AbstractNLPModel; options... -) - get_index_constraints( - get_lvar(nlp), get_uvar(nlp), - get_lcon(nlp), get_ucon(nlp); - options... - ) -end - -function get_index_constraints( - lvar, uvar, - lcon, ucon; - fixed_variable_treatment=MakeParameter, - equality_treatment=EnforceEquality, -) - ncon = length(lcon) - - if ncon > 0 - if equality_treatment == EnforceEquality - ind_eq = findall(lcon .== ucon) - ind_ineq = findall(lcon .!= ucon) - else - ind_eq = similar(lvar, Int, 0) - ind_ineq = similar(lvar, Int, ncon) .= 1:ncon - end - xl = [lvar;view(lcon,ind_ineq)] - xu = [uvar;view(ucon,ind_ineq)] - else - ind_eq = similar(lvar, Int, 0) - ind_ineq = similar(lvar, Int, 0) - xl = lvar - xu = uvar - end - - if fixed_variable_treatment == MakeParameter - ind_fixed = findall(xl .== xu) - ind_lb = findall((xl .!= -Inf) .* (xl .!= xu)) - ind_ub = findall((xu .!= Inf) .* (xl .!= xu)) - else - ind_fixed = similar(xl, Int, 0) - ind_lb = findall(xl .!=-Inf) - ind_ub = findall(xu .!= Inf) - end - - ind_llb = findall((lvar .!= -Inf).*(uvar .== Inf)) - ind_uub = findall((lvar .== -Inf).*(uvar .!= Inf)) - - # Return named tuple - return ( - ind_eq = ind_eq, - ind_ineq = ind_ineq, - ind_fixed = ind_fixed, - ind_lb = ind_lb, - ind_ub = ind_ub, - ind_llb = ind_llb, - ind_uub = ind_uub, - ) -end - """ AbstractCallback{T, VT} @@ -144,7 +113,7 @@ An `AbstractCallback` handles the scaling of the problem and the reformulations of the equality constraints and fixed variables. """ -abstract type AbstractCallback{T,VT} end +abstract type AbstractCallback{T, VT, FH} end """ create_callback( @@ -172,13 +141,12 @@ Wrap an `AbstractNLPModel` using sparse structures. """ struct SparseCallback{ T, - VT <: AbstractVector{T}, - VI <: AbstractVector{Int}, - I <: AbstractNLPModel{T, VT}, - FH <: AbstractFixedVariableTreatment, - EH <: AbstractEqualityTreatment, - } <: AbstractCallback{T, VT} - + VT<:AbstractVector{T}, + VI<:AbstractVector{Int}, + I<:AbstractNLPModel{T,VT}, + FH<:AbstractFixedVariableTreatment, + EH<:AbstractEqualityTreatment, +} <: AbstractCallback{T,VT, FH} nlp::I nvar::Int ncon::Int @@ -201,6 +169,14 @@ struct SparseCallback{ fixed_handler::FH equality_handler::EH + + ind_eq::VI + ind_ineq::VI + ind_fixed::VI + ind_lb::VI + ind_ub::VI + ind_llb::VI + ind_uub::VI end """ @@ -211,19 +187,17 @@ Wrap an `AbstractNLPModel` using dense structures. """ struct DenseCallback{ T, - VT <: AbstractVector{T}, - MT <: AbstractMatrix{T}, - I <: AbstractNLPModel{T, VT}, - FH <: AbstractFixedVariableTreatment, - EH <: AbstractEqualityTreatment, - } <: AbstractCallback{T, VT} - + VT<:AbstractVector{T}, + VI<:AbstractVector{Int}, + I<:AbstractNLPModel{T,VT}, + FH<:AbstractFixedVariableTreatment, + EH<:AbstractEqualityTreatment, +} <: AbstractCallback{T,VT, FH} nlp::I nvar::Int ncon::Int con_buffer::VT - jac_buffer::MT grad_buffer::VT obj_scale::Base.RefValue{T} @@ -231,19 +205,34 @@ struct DenseCallback{ fixed_handler::FH equality_handler::EH -end + ind_eq::VI + ind_ineq::VI + ind_fixed::VI + ind_lb::VI + ind_ub::VI + ind_llb::VI + ind_uub::VI +end create_array(cb::AbstractCallback, args...) = similar(get_x0(cb.nlp), args...) -function set_obj_scale!(obj_scale, f::VT, max_gradient) where {T, VT <: AbstractVector{T}} - obj_scale[] = min(one(T), max_gradient / norm(f, Inf)) +#= + Scaling +=# +function set_obj_scale!(obj_scale, f::VT, max_gradient) where {T,VT<:AbstractVector{T}} + return obj_scale[] = min(one(T), max_gradient / norm(f, Inf)) end -function set_con_scale_sparse!(con_scale::VT, jac_I,jac_buffer, max_gradient) where {T, VT <: AbstractVector{T}} +function set_con_scale_sparse!( + con_scale::VT, + jac_I, + jac_buffer, + max_gradient, +) where {T,VT<:AbstractVector{T}} fill!(con_scale, one(T)) _set_con_scale_sparse!(con_scale, jac_I, jac_buffer) - map!(x-> min(one(T), max_gradient / x), con_scale, con_scale) + return map!(x -> min(one(T), max_gradient / x), con_scale, con_scale) end function _set_con_scale_sparse!(con_scale, jac_I, jac_buffer) @inbounds @simd for i in 1:length(jac_I) @@ -252,29 +241,24 @@ function _set_con_scale_sparse!(con_scale, jac_I, jac_buffer) end end -function set_jac_scale_sparse!(jac_scale::VT, con_scale, jac_I) where {T, VT <: AbstractVector{T}} - copyto!(jac_scale, @view(con_scale[jac_I])) -end - -function set_con_scale_dense!(con_scale::VT, jac_buffer, max_gradient) where {T, VT <: AbstractVector{T}} - con_scale .= min.(one(T), max_gradient ./ mapreduce(abs, max, jac_buffer, dims=2, init=one(T))) +function set_jac_scale_sparse!( + jac_scale::VT, + con_scale, + jac_I, +) where {T,VT<:AbstractVector{T}} + return copyto!(jac_scale, @view(con_scale[jac_I])) end - -function create_dense_fixed_handler( - fixed_variable_treatment::Type{MakeParameter}, - nlp, -) - lvar = get_lvar(nlp) - uvar = get_uvar(nlp) - isfixed = (lvar .== uvar) - fixed = findall(isfixed) - return MakeParameter( - fixed, - similar(fixed,0), - similar(fixed,0), - similar(lvar, length(fixed)) - ) +function set_con_scale_dense!( + con_scale::VT, + jac_buffer, + max_gradient, +) where {T,VT<:AbstractVector{T}} + return con_scale .= + min.( + one(T), + max_gradient ./ mapreduce(abs, max, jac_buffer, dims = 2, init = one(T)), + ) end function create_sparse_fixed_handler( @@ -286,33 +270,93 @@ function create_sparse_fixed_handler( hess_J, hess_buffer, ) + n = get_nvar(nlp) lvar = get_lvar(nlp) uvar = get_uvar(nlp) nnzj = get_nnzj(nlp.meta) nnzh = get_nnzh(nlp.meta) - isfixed = (lvar .== uvar) + isfixed = (lvar .== uvar) + isfree = (lvar .< uvar) - fixed = findall(isfixed) - fixedj = findall(@view(isfixed[jac_J])) - fixedh = findall(@view(isfixed[hess_I]) .|| @view(isfixed[hess_J])) + fixed = findall(isfixed) nfixed = length(fixed) - nnzh = nnzh + nfixed + if nfixed == 0 + return NoFixedVariables(), n, nnzj, nnzh + end + + free = findall(isfree) + nx = length(free) + map_full_to_free = similar(jac_I, n) ; fill!(map_full_to_free, -1) + map_full_to_free[free] .= 1:nx + + ind_jac_free = findall(@view(isfree[jac_J])) + ind_hess_free = findall(@view(isfree[hess_I]) .&& @view(isfree[hess_J])) + + nnzh = length(ind_hess_free) + Hi, Hj = similar(hess_I, nnzh), similar(hess_J, nnzh) + # TODO: vectorize + for k in 1:nnzh + cnt = ind_hess_free[k] + i, j = hess_I[cnt], hess_J[cnt] + Hi[k] = map_full_to_free[i] + Hj[k] = map_full_to_free[j] + end resize!(hess_I, nnzh) resize!(hess_J, nnzh) - resize!(hess_buffer, nnzh) - copyto!(@view(hess_I[end-nfixed+1:end]), fixed) - copyto!(@view(hess_J[end-nfixed+1:end]), fixed) + hess_I .= Hi + hess_J .= Hj + + nnzj = length(ind_jac_free) + Ji, Jj = similar(jac_I, nnzj), similar(jac_J, nnzj) + for k in 1:nnzj + cnt = ind_jac_free[k] + i, j = jac_I[cnt], jac_J[cnt] + Ji[k] = i + Jj[k] = map_full_to_free[j] + end + resize!(jac_I, nnzj) + resize!(jac_J, nnzj) + jac_I .= Ji + jac_J .= Jj + + x_full = copy(lvar) fixed_handler = MakeParameter( + free, fixed, - fixedj, - fixedh, - similar(lvar, length(fixed)) + ind_jac_free, + ind_hess_free, + Ref(NaN), + x_full, + similar(lvar, n), ) - return fixed_handler, nnzj, nnzh + return fixed_handler, nx, nnzj, nnzh +end + + +function create_dense_fixed_handler(fixed_variable_treatment::Type{MakeParameter}, nlp) + n = get_nvar(nlp) + lvar = get_lvar(nlp) + uvar = get_uvar(nlp) + isfixed = (lvar .== uvar) + fixed = findall(lvar .== uvar) + if length(fixed) == 0 + return NoFixedVariables() + else + free = findall(lvar .< uvar) + return MakeParameter( + free, + fixed, + similar(fixed, 0), + similar(fixed, 0), + Ref(NaN), + similar(lvar, 0), + similar(lvar, n), + ) + end end function create_sparse_fixed_handler( @@ -324,54 +368,107 @@ function create_sparse_fixed_handler( hess_J, hess_buffer, ) + n = get_nvar(nlp) fixed_handler = RelaxBound() - return fixed_handler, get_nnzj(nlp.meta), get_nnzh(nlp.meta) + return fixed_handler, n, get_nnzj(nlp.meta), get_nnzh(nlp.meta) end +#= + Constructors +=# + function create_callback( ::Type{SparseCallback}, - nlp::AbstractNLPModel{T, VT}; - fixed_variable_treatment=MakeParameter, - equality_treatment=EnforceEquality, - ) where {T, VT} - + nlp::AbstractNLPModel{T,VT}; + fixed_variable_treatment = MakeParameter, + equality_treatment = EnforceEquality, +) where {T,VT} n = get_nvar(nlp) m = get_ncon(nlp) nnzj = get_nnzj(nlp.meta) nnzh = get_nnzh(nlp.meta) - x0 = get_x0(nlp) - - con_buffer = similar(x0, m) ; fill!(con_buffer, zero(T)) - grad_buffer = similar(x0, n) ; fill!(grad_buffer, zero(T)) - jac_buffer = similar(x0, nnzj) ; fill!(jac_buffer, zero(T)) - hess_buffer = similar(x0, nnzh) ; fill!(hess_buffer, zero(T)) + x0 = get_x0(nlp) jac_I = similar(x0, Int, nnzj) jac_J = similar(x0, Int, nnzj) hess_I = similar(x0, Int, nnzh) hess_J = similar(x0, Int, nnzh) + jac_buffer = similar(x0, nnzj) + fill!(jac_buffer, zero(T)) obj_scale = Ref(one(T)) - con_scale = similar(jac_buffer, m) ; fill!(con_scale, one(T)) - jac_scale = similar(jac_buffer, nnzj) ; fill!(jac_scale, one(T)) + con_scale = similar(jac_buffer, m) + fill!(con_scale, one(T)) NLPModels.jac_structure!(nlp, jac_I, jac_J) if nnzh > 0 NLPModels.hess_structure!(nlp, hess_I, hess_J) end - fixed_handler, nnzj, nnzh = create_sparse_fixed_handler( + hess_buffer = similar(x0, nnzh) + fill!(hess_buffer, zero(T)) + fixed_handler, nx, nnzj, nnzh = create_sparse_fixed_handler( fixed_variable_treatment, nlp, - jac_I, jac_J, hess_I, hess_J, + jac_I, + jac_J, + hess_I, + hess_J, hess_buffer, ) equality_handler = equality_treatment() + jac_scale = similar(jac_buffer, nnzj) + fill!(jac_scale, one(T)) + con_buffer = similar(x0, m) + fill!(con_buffer, zero(T)) + grad_buffer = similar(x0, nx) + fill!(grad_buffer, zero(T)) + + # Get indexing + lvar = get_lvar(nlp) + uvar = get_uvar(nlp) + lcon = get_lcon(nlp) + ucon = get_ucon(nlp) + + # Get fixed variables + ind_fixed = findall(lvar .== uvar) + if length(ind_fixed) > 0 && fixed_variable_treatment == MakeParameter + ind_free = findall(lvar .< uvar) + # Remove fixed variables from problem's formulation + lvar = lvar[ind_free] + uvar = uvar[ind_free] + end + + if m > 0 + if equality_treatment == EnforceEquality + ind_eq = findall(lcon .== ucon) + ind_ineq = findall(lcon .!= ucon) + else + ind_eq = similar(lvar, Int, 0) + ind_ineq = similar(lvar, Int, m) .= 1:m + end + xl = [lvar; view(lcon, ind_ineq)] + xu = [uvar; view(ucon, ind_ineq)] + else + ind_eq = similar(lvar, Int, 0) + ind_ineq = similar(lvar, Int, 0) + xl = lvar + xu = uvar + end + + ind_llb = findall((lvar .!= -Inf) .* (uvar .== Inf)) + ind_uub = findall((lvar .== -Inf) .* (uvar .!= Inf)) + ind_lb = findall(xl .!= -Inf) + ind_ub = findall(xu .!= Inf) + return SparseCallback( nlp, - n,m,nnzj,nnzh, + nx, + m, + nnzj, + nnzh, con_buffer, jac_buffer, grad_buffer, @@ -384,152 +481,246 @@ function create_callback( con_scale, jac_scale, fixed_handler, - equality_handler + equality_handler, + ind_eq, + ind_ineq, + ind_fixed, + ind_lb, + ind_ub, + ind_llb, + ind_uub, ) end function create_callback( ::Type{DenseCallback}, - nlp::AbstractNLPModel{T, VT}; - fixed_variable_treatment=MakeParameter, - equality_treatment=EnforceEquality, - ) where {T, VT} - + nlp::AbstractNLPModel{T,VT}; + fixed_variable_treatment = MakeParameter, + equality_treatment = EnforceEquality, +) where {T,VT} n = get_nvar(nlp) m = get_ncon(nlp) - x0 = similar(get_x0(nlp)) - con_buffer = similar(x0, m) ; fill!(con_buffer, zero(T)) - jac_buffer = similar(x0, m, n) ; fill!(jac_buffer, zero(T)) - grad_buffer = similar(x0, n) ; fill!(grad_buffer, zero(T)) + x0 = similar(get_x0(nlp)) + con_buffer = similar(x0, m) + fill!(con_buffer, zero(T)) + jac_buffer = similar(x0, m, n) + fill!(jac_buffer, zero(T)) + grad_buffer = similar(x0, n) + fill!(grad_buffer, zero(T)) obj_scale = Ref(one(T)) - con_scale = similar(x0, m) ; fill!(con_scale, one(T)) + con_scale = similar(x0, m) + fill!(con_scale, one(T)) - fixed_handler = create_dense_fixed_handler( - fixed_variable_treatment, - nlp, - ) + fixed_handler = create_dense_fixed_handler(fixed_variable_treatment, nlp) equality_handler = equality_treatment() + # Get indexing + lvar = get_lvar(nlp) + uvar = get_uvar(nlp) + lcon = get_lcon(nlp) + ucon = get_ucon(nlp) + + if m > 0 + if equality_treatment == EnforceEquality + ind_eq = findall(lcon .== ucon) + ind_ineq = findall(lcon .!= ucon) + else + ind_eq = similar(lvar, Int, 0) + ind_ineq = similar(lvar, Int, ncon) .= 1:ncon + end + xl = [lvar; view(lcon, ind_ineq)] + xu = [uvar; view(ucon, ind_ineq)] + else + ind_eq = similar(lvar, Int, 0) + ind_ineq = similar(lvar, Int, 0) + xl = lvar + xu = uvar + end + + ind_llb = findall((lvar .!= -Inf) .* (uvar .== Inf)) + ind_uub = findall((lvar .== -Inf) .* (uvar .!= Inf)) + + # Get fixed variables + ind_fixed = findall(lvar .== uvar) + if length(ind_fixed) > 0 && fixed_variable_treatment == MakeParameter + # Keep fixed variables but remove them from lb/ub + ind_lb = findall((xl .!= -Inf) .* (xl .!= xu)) + ind_ub = findall((xu .!= Inf) .* (xl .!= xu)) + else + ind_lb = findall(xl .!= -Inf) + ind_ub = findall(xu .!= Inf) + end + return DenseCallback( nlp, - n, m, + n, + m, con_buffer, - jac_buffer, grad_buffer, obj_scale, con_scale, fixed_handler, - equality_handler + equality_handler, + ind_eq, + ind_ineq, + ind_fixed, + ind_lb, + ind_ub, + ind_llb, + ind_uub, ) end -function _treat_fixed_variable_initialize!(fixed_handler::RelaxBound, x0, lvar, uvar) end -function _treat_fixed_variable_initialize!(fixed_handler::MakeParameter, x0, lvar, uvar) - fixed = fixed_handler.fixed - copyto!(@view(x0[fixed]), @view(lvar[fixed])) - fill!(@view(lvar[fixed]), -Inf) - fill!(@view(uvar[fixed]), Inf) -end +#= + Callback's initialization +=# function _treat_equality_initialize!(equality_handler::EnforceEquality, lcon, ucon, tol) end function _treat_equality_initialize!(equality_handler::RelaxEquality, lcon, ucon, tol) - set_initial_bounds!( - lcon, - ucon, - tol - ) + return set_initial_bounds!(lcon, ucon, tol) +end +# Initiate fixed variables. By default, do nothing. +function _treat_fixed_variable_initialize!(cb::AbstractCallback, x0, lvar, uvar) end +function _treat_fixed_variable_initialize!( + cb::DenseCallback{T, VT, VI, NLP, FH}, + x0, + lvar, + uvar, +) where {T, VT, VI, NLP, FH<:MakeParameter} + x0[cb.fixed_handler.fixed] .= lvar[cb.fixed_handler.fixed] + lvar[cb.fixed_handler.fixed] .= -Inf + uvar[cb.fixed_handler.fixed] .= Inf + return end function initialize!( cb::AbstractCallback, - x, xl, xu, y0, rhs, + x, + xl, + xu, + y0, + rhs, ind_ineq; - tol=1e-8, - bound_push=1e-2, - bound_fac=1e-2, + tol = 1e-8, + bound_push = 1e-2, + bound_fac = 1e-2, ) - - x0= variable(x) - lvar= variable(xl) - uvar= variable(xu) + x0 = variable(x) + lvar = variable(xl) + uvar = variable(xu) nlp = cb.nlp - con_buffer =cb.con_buffer - grad_buffer =cb.grad_buffer - + con_buffer = cb.con_buffer + grad_buffer = cb.grad_buffer - x0 .= get_x0(nlp) - y0 .= get_y0(nlp) - lvar .= get_lvar(nlp) - uvar .= get_uvar(nlp) + x0 .= get_x0(cb) + y0 .= get_y0(cb) + lvar .= get_lvar(cb) + uvar .= get_uvar(cb) lcon = copy(get_lcon(nlp)) ucon = copy(get_ucon(nlp)) - _treat_fixed_variable_initialize!(cb.fixed_handler, x0, lvar, uvar) _treat_equality_initialize!(cb.equality_handler, lcon, ucon, tol) + _treat_fixed_variable_initialize!(cb, x0, lvar, uvar) - set_initial_bounds!( - lvar, - uvar, - tol - ) - initialize_variables!( - x0, - lvar, - uvar, - bound_push, - bound_fac - ) + set_initial_bounds!(lvar, uvar, tol) + initialize_variables!(x0, lvar, uvar, bound_push, bound_fac) - NLPModels.cons!(nlp,x0,con_buffer) + _eval_cons_wrapper!(cb, x0, con_buffer) slack(xl) .= view(lcon, ind_ineq) slack(xu) .= view(ucon, ind_ineq) - rhs .= (lcon.==ucon).*lcon + rhs .= (lcon .== ucon) .* lcon copyto!(slack(x), @view(con_buffer[ind_ineq])) - set_initial_bounds!( - slack(xl), - slack(xu), - tol - ) - initialize_variables!( - slack(x), - slack(xl), - slack(xu), - bound_push, - bound_fac - ) + set_initial_bounds!(slack(xl), slack(xu), tol) + initialize_variables!(slack(x), slack(xl), slack(xu), bound_push, bound_fac) + return +end + +#= + Getters +=# + +n_variables(cb::AbstractCallback) = get_nvar(cb.nlp) +n_constraints(cb::AbstractCallback) = get_ncon(cb.nlp) +get_x0(cb::AbstractCallback) = get_x0(cb.nlp) +get_y0(cb::AbstractCallback) = get_y0(cb.nlp) +get_lvar(cb::AbstractCallback) = get_lvar(cb.nlp) +get_uvar(cb::AbstractCallback) = get_uvar(cb.nlp) +# Getters to unpack solution +unpack_obj(cb::AbstractCallback, obj_val) = obj_val ./ cb.obj_scale[] +function unpack_cons!(c_full::AbstractVector, cb::AbstractCallback, c::AbstractVector) + c_full .= c ./ cb.con_scale +end +function unpack_x!(x_full, cb::AbstractCallback, x) + x_full .= x +end +function unpack_z!(z_full, cb::AbstractCallback, z) + z_full .= z end +# N.B.: Special getters if we use SparseCallback with a MakeParameter fixed_handler, +# as the dimension of the problem is reduced. +function n_variables(cb::SparseCallback{T, VT, VI, NLP, FH}) where {T, VT, VI, NLP, FH<:MakeParameter} + return length(cb.fixed_handler.free) +end +function get_x0(cb::SparseCallback{T, VT, VI, NLP, FH}) where {T, VT, VI, NLP, FH<:MakeParameter} + free = cb.fixed_handler.free + x0 = get_x0(cb.nlp) + return x0[free] +end +function get_lvar(cb::SparseCallback{T, VT, VI, NLP, FH}) where {T, VT, VI, NLP, FH<:MakeParameter} + free = cb.fixed_handler.free + xl = get_lvar(cb.nlp) + return xl[free] +end +function get_uvar(cb::SparseCallback{T, VT, VI, NLP, FH}) where {T, VT, VI, NLP, FH<:MakeParameter} + free = cb.fixed_handler.free + xu = get_uvar(cb.nlp) + return xu[free] +end +function unpack_x!(x_full, cb::SparseCallback{T, VT, VI, NLP, FH}, x) where {T, VT, VI, NLP, FH<:MakeParameter} + _update_x!(cb.fixed_handler, x) + x_full .= cb.fixed_handler.x_full +end +function unpack_z!(z_full, cb::SparseCallback{T, VT, VI, NLP, FH}, z) where {T, VT, VI, NLP, FH<:MakeParameter} + free = cb.fixed_handler.free + z_full[free] .= z +end + + function set_scaling!( cb::SparseCallback, - x, xl, xu, y0, rhs, + x, + xl, + xu, + y0, + rhs, ind_ineq, - nlp_scaling_max_gradient - ) - - x0= variable(x) + nlp_scaling_max_gradient, +) + x0 = variable(x) nlp = cb.nlp obj_scale = cb.obj_scale con_scale = cb.con_scale jac_scale = cb.jac_scale - con_buffer =cb.con_buffer - jac_buffer =cb.jac_buffer - grad_buffer =cb.grad_buffer + con_buffer = cb.con_buffer + grad_buffer = cb.grad_buffer # Set scaling - NLPModels.jac_coord!(nlp,x0,jac_buffer) - set_con_scale_sparse!(con_scale, cb.jac_I, jac_buffer, nlp_scaling_max_gradient) + jac = similar(con_buffer, cb.nnzj) + _eval_jac_wrapper!(cb, x0, jac) + set_con_scale_sparse!(con_scale, cb.jac_I, jac, nlp_scaling_max_gradient) set_jac_scale_sparse!(jac_scale, con_scale, cb.jac_I) - NLPModels.grad!(nlp,x0,grad_buffer) + _eval_grad_f_wrapper!(cb, x0, grad_buffer) set_obj_scale!(obj_scale, grad_buffer, nlp_scaling_max_gradient) con_scale_slk = @view(con_scale[ind_ineq]) - y0 ./= con_scale + y0 ./= con_scale rhs .*= con_scale slack(x) .*= con_scale_slk slack(xl) .*= con_scale_slk @@ -539,29 +730,34 @@ end function set_scaling!( cb::DenseCallback, - x, xl, xu, y0, rhs, + x, + xl, + xu, + y0, + rhs, ind_ineq, - nlp_scaling_max_gradient - ) - + nlp_scaling_max_gradient, +) x0 = variable(x) nlp = cb.nlp obj_scale = cb.obj_scale con_scale = cb.con_scale - con_buffer =cb.con_buffer - jac_buffer =cb.jac_buffer - grad_buffer =cb.grad_buffer + con_buffer = cb.con_buffer + grad_buffer = cb.grad_buffer + # N.B.: Allocate jac_buffer here as it is use only once. + # GC should take care of it once the scaling has been initialized. + jac_buffer = similar(grad_buffer, cb.ncon, cb.nvar) # Set scaling - jac_dense!(nlp,x0,jac_buffer) + jac_dense!(nlp, x0, jac_buffer) set_con_scale_dense!(con_scale, jac_buffer, nlp_scaling_max_gradient) - NLPModels.grad!(nlp,x0,grad_buffer) + NLPModels.grad!(nlp, x0, grad_buffer) set_obj_scale!(obj_scale, grad_buffer, nlp_scaling_max_gradient) con_scale_slk = @view(con_scale[ind_ineq]) - y0 ./= con_scale + y0 ./= con_scale rhs .*= con_scale slack(x) .*= con_scale_slk slack(xl) .*= con_scale_slk @@ -569,217 +765,322 @@ function set_scaling!( return end -function _jac_sparsity_wrapper!( - cb::SparseCallback, - I::AbstractVector,J::AbstractVector - ) - - copyto!(I, cb.jac_I) - copyto!(J, cb.jac_J) - return -end +#= + Callbacks: default +=# -function _hess_sparsity_wrapper!( - cb::SparseCallback, - I::AbstractVector,J::AbstractVector - ) - copyto!(I, cb.hess_I) - copyto!(J, cb.hess_J) - return +function _eval_f_wrapper(cb::AbstractCallback, x::AbstractVector) + return NLPModels.obj(cb.nlp, x) * cb.obj_scale[] end - -function _eval_cons_wrapper!(cb::AbstractCallback,x::AbstractVector,c::AbstractVector) - NLPModels.cons!(cb.nlp, x,c) +function _eval_cons_wrapper!(cb::AbstractCallback, x::AbstractVector, c::AbstractVector) + NLPModels.cons!(cb.nlp, x, c) c .*= cb.con_scale return c end -function _eval_jac_wrapper!( - cb::SparseCallback, +function _eval_grad_f_wrapper!( + cb::AbstractCallback{T}, x::AbstractVector, - jac::AbstractVector - ) - - nnzj_orig = get_nnzj(cb.nlp.meta) - NLPModels.jac_coord!(cb.nlp, x, jac) - jac .*= cb.jac_scale - - _treat_fixed_variable_jac_coord!(cb.fixed_handler, cb, x, jac) + grad::AbstractVector, +) where {T} + NLPModels.grad!(cb.nlp, x, grad) + grad .*= cb.obj_scale[] + return grad end function _eval_jtprod_wrapper!( - cb::AbstractCallback{T}, + cb::AbstractCallback, x::AbstractVector, v::AbstractVector, jvt::AbstractVector, - ) where T - +) y = cb.con_buffer y .= v .* cb.con_scale NLPModels.jtprod!(cb.nlp, x, y, jvt) - _treat_fixed_variable_grad!(cb.fixed_handler, cb, x, jvt) return jvt end -function _treat_fixed_variable_jac_coord!(fixed_handler::RelaxBound, cb, x, jac) end -function _treat_fixed_variable_jac_coord!(fixed_handler::MakeParameter, cb::SparseCallback{T}, x, jac) where T - fill!(@view(jac[fixed_handler.fixedj]), zero(T)) -end +#= + Callbacks: SparseCallback +=# -function _eval_grad_f_wrapper!( - cb::AbstractCallback{T}, - x::AbstractVector, - grad::AbstractVector - ) where T +function _jac_sparsity_wrapper!(cb::SparseCallback, I::AbstractVector, J::AbstractVector) + copyto!(I, cb.jac_I) + copyto!(J, cb.jac_J) + return +end - NLPModels.grad!(cb.nlp, x, grad) - grad .*= cb.obj_scale[] - _treat_fixed_variable_grad!(cb.fixed_handler, cb, x, grad) +function _hess_sparsity_wrapper!(cb::SparseCallback, I::AbstractVector, J::AbstractVector) + copyto!(I, cb.hess_I) + copyto!(J, cb.hess_J) + return end -function _treat_fixed_variable_grad!(fixed_handler::RelaxBound, cb, x, grad) end -function _treat_fixed_variable_grad!(fixed_handler::MakeParameter, cb::AbstractCallback{T}, x, grad) where T - fixed_handler.grad_storage .= @view(grad[fixed_handler.fixed]) - map!( - (x,y)->x-y, - @view(grad[fixed_handler.fixed]), - @view(x[cb.fixed_handler.fixed]), - @view(get_lvar(cb.nlp)[cb.fixed_handler.fixed]) - ) + +function _eval_jac_wrapper!( + cb::SparseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, + jac::AbstractVector, +) where {T, VT, VI, NLP, FH} + NLPModels.jac_coord!(cb.nlp, x, jac) + jac .*= cb.jac_scale + return jac end -function _eval_f_wrapper(cb::AbstractCallback,x::AbstractVector) - return NLPModels.obj(cb.nlp,x)* cb.obj_scale[] +function _eval_jac_wrapper!( + cb::SparseCallback{T}, + x::AbstractVector, + jac::AbstractMatrix, +) where {T} + jac_buffer = view(cb.jac_buffer, 1:cb.nnzj) + _eval_jac_wrapper!(cb, x, jac_buffer) + fill!(jac, zero(T)) + @inbounds @simd for k in 1:length(cb.jac_I) + i, j = cb.jac_I[k], cb.jac_J[k] + jac[i, j] += jac_buffer[k] + end + return jac end function _eval_lag_hess_wrapper!( - cb::SparseCallback{T}, + cb::SparseCallback{T, VT, VI, NLP, FH}, x::AbstractVector, y::AbstractVector, hess::AbstractVector; - obj_weight = 1.0 - ) where T - - nnzh_orig = get_nnzh(cb.nlp.meta) - + obj_weight = 1.0, +) where {T, VT, VI, NLP, FH} cb.con_buffer .= y .* cb.con_scale + nnzh_orig = get_nnzh(cb.nlp.meta) NLPModels.hess_coord!( - cb.nlp, x, cb.con_buffer, view(hess, 1:nnzh_orig); - obj_weight=obj_weight * cb.obj_scale[] + cb.nlp, + x, + cb.con_buffer, + view(hess, 1:nnzh_orig); + obj_weight = obj_weight * cb.obj_scale[], ) - _treat_fixed_variable_hess_coord!(cb.fixed_handler, cb, hess) + return hess end -function _treat_fixed_variable_hess_coord!(fixed_handler::RelaxBound, cb, hess) end -function _treat_fixed_variable_hess_coord!(fixed_handler::MakeParameter, cb::SparseCallback{T}, hess) where T - nnzh_orig = get_nnzh(cb.nlp.meta) - fill!(@view(hess[fixed_handler.fixedh]), zero(T)) - fill!(@view(hess[nnzh_orig+1:end]), one(T)) +function _eval_lag_hess_wrapper!( + cb::SparseCallback{T}, + x::AbstractVector, + y::AbstractVector, + hess::AbstractMatrix; + obj_weight = one(T), +) where {T} + hess_buffer = view(cb.hess_buffer, 1:cb.nnzh) + _eval_lag_hess_wrapper!(cb, x, y, hess_buffer; obj_weight = obj_weight * cb.obj_scale[]) + fill!(hess, zero(T)) + @inbounds @simd for k in 1:length(cb.hess_I) + i, j = cb.hess_I[k], cb.hess_J[k] + hess[i, j] += hess_buffer[k] + end + return hess end +#= + Callback: DenseCallback +=# + function _eval_jac_wrapper!( - cb::SparseCallback{T}, + cb::DenseCallback{T, VT, VI, NLP, FH}, x::AbstractVector, - jac::AbstractMatrix - ) where T - - jac_buffer = cb.jac_buffer - _eval_jac_wrapper!(cb, x, jac_buffer) - fill!(jac, zero(T)) - @inbounds @simd for k=1:length(cb.jac_I) - i,j = cb.jac_I[k], cb.jac_J[k] - jac[i,j] += jac_buffer[k] - end + jac::AbstractMatrix, +) where {T, VT, VI, NLP, FH} + jac_dense!(cb.nlp, x, jac) + jac .*= cb.con_scale + return jac end function _eval_lag_hess_wrapper!( - cb::SparseCallback{T}, + cb::DenseCallback{T, VT, VI, NLP, FH}, x::AbstractVector, y::AbstractVector, hess::AbstractMatrix; - obj_weight = one(T) - ) where T + obj_weight = one(T), +) where {T, VT, VI, NLP, FH} + hess_dense!(cb.nlp, x, y, hess; obj_weight = obj_weight * cb.obj_scale[]) + return hess +end - hess_buffer = cb.hess_buffer - _eval_lag_hess_wrapper!(cb, x, y, hess_buffer; obj_weight=obj_weight * cb.obj_scale[]) - fill!(hess, zero(T)) - @inbounds @simd for k=1:length(cb.hess_I) - i,j = cb.hess_I[k], cb.hess_J[k] - hess[i,j] += hess_buffer[k] +#= + Callback for SparseCallback+MakeParameter +=# + +function _update_x!(fixed_handler::MakeParameter, x::AbstractVector) + idx = norm(x, 2) + if fixed_handler.hash_x[] == idx + return end - _treat_fixed_variable_hess_dense!(cb.fixed_handler, cb, hess) + fixed_handler.hash_x[] = idx + # Update x_full + free = fixed_handler.free + fixed_handler.x_full[free] .= x + return end -function _treat_fixed_variable_hess_dense!(fixed_handler::RelaxBound, cb, hess) end -function _treat_fixed_variable_hess_dense!(fixed_handler::MakeParameter, cb::SparseCallback{T}, hess) where T - nnzh_orig = get_nnzh(cb.nlp.meta) - fixed = fixed_handler.fixed - _set_diag!(hess, fixed, one(T)) +function _eval_f_wrapper( + cb::SparseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, +) where {T, VT, VI, NLP, FH<:MakeParameter} + _update_x!(cb.fixed_handler, x) + x_full = cb.fixed_handler.x_full + return NLPModels.obj(cb.nlp, x_full) * cb.obj_scale[] end -function _eval_jac_wrapper!( - cb::DenseCallback{T}, +function _eval_cons_wrapper!( + cb::SparseCallback{T, VT, VI, NLP, FH}, x::AbstractVector, - jac::AbstractMatrix - ) where T + c::AbstractVector, +) where {T, VT, VI, NLP, FH<:MakeParameter} + _update_x!(cb.fixed_handler, x) + x_full = cb.fixed_handler.x_full + NLPModels.cons!(cb.nlp, x_full, c) + c .*= cb.con_scale + return c +end - jac_dense!(cb.nlp, x, jac) - jac .*= cb.con_scale - _treat_fixed_variable_jac_dense!(cb.fixed_handler, cb, jac) +function _eval_grad_f_wrapper!( + cb::SparseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, + grad::AbstractVector, +) where {T, VT, VI, NLP, FH<:MakeParameter} + _update_x!(cb.fixed_handler, x) + x_full = cb.fixed_handler.x_full + g_full = cb.fixed_handler.g_full + NLPModels.grad!(cb.nlp, x_full, g_full) + grad .= g_full[cb.fixed_handler.free] + grad .*= cb.obj_scale[] + return grad +end + +function _eval_jac_wrapper!( + cb::SparseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, + jac::AbstractVector, +) where {T, VT, VI, NLP, FH<:MakeParameter} + _update_x!(cb.fixed_handler, x) + nnzj_orig = get_nnzj(cb.nlp.meta) + jac_full = cb.jac_buffer + x_full = cb.fixed_handler.x_full + @assert length(jac_full) == nnzj_orig + NLPModels.jac_coord!(cb.nlp, x_full, jac_full) + jac .= jac_full[cb.fixed_handler.ind_jac_free] + jac .*= cb.jac_scale + return jac end -function _treat_fixed_variable_jac_dense!(fixed_handler::RelaxBound, cb::DenseCallback, jac) end -function _treat_fixed_variable_jac_dense!(fixed_handler::MakeParameter, cb::DenseCallback{T}, jac) where T - jac[:,fixed_handler.fixed] .= zero(T) + +function _eval_jtprod_wrapper!( + cb::SparseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, + v::AbstractVector, + jvt::AbstractVector, +) where {T, VT, VI, NLP, FH<:MakeParameter} + _update_x!(cb.fixed_handler, x) + x_full = cb.fixed_handler.x_full + jvt_full = cb.fixed_handler.g_full + y = cb.con_buffer + y .= v .* cb.con_scale + NLPModels.jtprod!(cb.nlp, x_full, y, jvt_full) + jvt .= jvt_full[cb.fixed_handler.free] + return jvt end function _eval_lag_hess_wrapper!( - cb::DenseCallback{T}, + cb::SparseCallback{T, VT, VI, NLP, FH}, x::AbstractVector, y::AbstractVector, - hess::AbstractMatrix; - obj_weight = one(T) - ) where T - - hess_dense!( - cb.nlp, x, y, hess; - obj_weight=obj_weight * cb.obj_scale[] + hess::AbstractVector; + obj_weight = 1.0, +) where {T, VT, VI, NLP, FH<:MakeParameter} + _update_x!(cb.fixed_handler, x) + x_full = cb.fixed_handler.x_full + nnzh_orig = get_nnzh(cb.nlp.meta) + hess_full = cb.hess_buffer + @assert length(hess_full) == nnzh_orig + cb.con_buffer .= y .* cb.con_scale + NLPModels.hess_coord!( + cb.nlp, + x_full, + cb.con_buffer, + hess_full; + obj_weight = obj_weight * cb.obj_scale[], ) + hess .= hess_full[cb.fixed_handler.ind_hess_free] + return hess +end - _treat_fixed_variable_lag_hess_dense!(cb.fixed_handler, cb, hess) +#= + Callback for DenseCallback+MakeParameter +=# + +function _eval_grad_f_wrapper!( + cb::DenseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, + grad::AbstractVector, +) where {T, VT, VI, NLP, FH<:MakeParameter} + NLPModels.grad!(cb.nlp, x, grad) + grad .*= cb.obj_scale[] + map!( + (x, y) -> x - y, + @view(grad[cb.fixed_handler.fixed]), + @view(x[cb.fixed_handler.fixed]), + @view(get_lvar(cb.nlp)[cb.fixed_handler.fixed]) + ) + return grad end -function _treat_fixed_variable_lag_hess_dense!(fixed_handler::RelaxBound, cb::DenseCallback, hess) end -function _treat_fixed_variable_lag_hess_dense!(fixed_handler::MakeParameter, cb::DenseCallback{T}, hess) where T - fixed = fixed_handler.fixed - hess[:,fixed] .= zero(T) - hess[fixed,:] .= zero(T) - _set_diag!(hess, fixed, one(T)) + +function _eval_jac_wrapper!( + cb::DenseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, + jac::AbstractMatrix, +) where {T, VT, VI, NLP, FH<:MakeParameter} + jac_dense!(cb.nlp, x, jac) + jac .*= cb.con_scale + return jac[:, cb.fixed_handler.fixed] .= zero(T) end -function update_z!(cb, zl, zu, jacl) - _update_z!(cb.fixed_handler, zl, zu, jacl, get_minimize(cb.nlp) ? 1 : -1) +function _eval_lag_hess_wrapper!( + cb::DenseCallback{T, VT, VI, NLP, FH}, + x::AbstractVector, + y::AbstractVector, + hess::AbstractMatrix; + obj_weight = one(T), +) where {T, VT, VI, NLP, FH<:MakeParameter} + hess_dense!(cb.nlp, x, y, hess; obj_weight = obj_weight * cb.obj_scale[]) + fixed = cb.fixed_handler.fixed + hess[:, fixed] .= zero(T) + hess[fixed, :] .= zero(T) + _set_diag!(hess, fixed, one(T)) + return hess end -function _update_z!(fixed_handler::MakeParameter, zl, zu, jacl, sense) - zl_r = @view(zl[fixed_handler.fixed]) - zu_r = @view(zu[fixed_handler.fixed]) - jacl_r = @view(jacl[fixed_handler.fixed]) - map!( - (x,y)->sense * max(x+y,0), - zl_r, - fixed_handler.grad_storage, - jacl_r - ) - map!( - (x,y)->sense * max(-(x+y),0), - zu_r, - fixed_handler.grad_storage, - jacl_r, - ) +#= + Compute bounds' multipliers for fixed variables + + At a KKT solution, we have ∇f + ∇cᵀ y - zl + zu = 0 , (zl, zu) >= 0 +=# + +# N.B.: by default do nothing as the bounds' multipliers are computed by the algorithm +function update_z!(cb, x, y, zl, zu, jacl) end + +function update_z!(cb::AbstractCallback{T, VT, FH}, x, y, zl, zu, jacl) where {T, VT, FH<:MakeParameter} + fixed_handler = cb.fixed_handler::MakeParameter + sense = get_minimize(cb.nlp) ? 1 : -1 + ind_fixed = fixed_handler.fixed + g_full = fixed_handler.g_full + jtv = similar(g_full) ; fill!(jtv, zero(T)) + NLPModels.grad!(cb.nlp, x, g_full) # ∇f + NLPModels.jtprod!(cb.nlp, x, y, jtv) # ∇cᵀ y + g_full .+= jtv # ∇f + ∇cᵀ y + g_fixed = view(g_full, ind_fixed) + zl[ind_fixed] .= sense .* max.(zero(T), g_fixed) + zu[ind_fixed] .= sense .* max.(zero(T), .-g_fixed) + return end -function _update_z!(fixed_handler::RelaxBound, zl, zu, jacl, sense) end function _set_diag!(A, inds, a) @inbounds @simd for i in inds - A[i,i] = a + A[i, i] = a end end diff --git a/src/options.jl b/src/options.jl index 438b4207..78665e2c 100644 --- a/src/options.jl +++ b/src/options.jl @@ -45,8 +45,8 @@ end # NLP options kappa_d::Float64 = 1e-5 - fixed_variable_treatment::Type = kkt_system <: MadNLP.SparseCondensedKKTSystem ? MadNLP.RelaxBound : MadNLP.MakeParameter - equality_treatment::Type = kkt_system <: MadNLP.SparseCondensedKKTSystem ? MadNLP.RelaxEquality : MadNLP.EnforceEquality + fixed_variable_treatment::Type = kkt_system <: SparseCondensedKKTSystem ? RelaxBound : MakeParameter + equality_treatment::Type = kkt_system <: SparseCondensedKKTSystem ? RelaxEquality : EnforceEquality boudn_relax_factor::Float64 = 1e-8 jacobian_constant::Bool = false hessian_constant::Bool = false @@ -57,7 +57,7 @@ end # initialization options dual_initialized::Bool = false - dual_initialization_method::Type = kkt_system <: MadNLP.SparseCondensedKKTSystem ? DualInitializeSetZero : DualInitializeLeastSquares + dual_initialization_method::Type = kkt_system <: SparseCondensedKKTSystem ? DualInitializeSetZero : DualInitializeLeastSquares constr_mult_init_max::Float64 = 1e3 bound_push::Float64 = 1e-2 bound_fac::Float64 = 1e-2 @@ -102,8 +102,8 @@ end mu_linear_decrease_factor::Float64 = .2 end -is_dense_callback(nlp) = hasmethod(MadNLP.jac_dense!, Tuple{typeof(nlp), AbstractVector, AbstractMatrix}) && - hasmethod(MadNLP.hess_dense!, Tuple{typeof(nlp), AbstractVector, AbstractVector, AbstractMatrix}) +is_dense_callback(nlp) = hasmethod(jac_dense!, Tuple{typeof(nlp), AbstractVector, AbstractMatrix}) && + hasmethod(hess_dense!, Tuple{typeof(nlp), AbstractVector, AbstractVector, AbstractMatrix}) # smart option presets function MadNLPOptions( @@ -171,9 +171,9 @@ function _get_primary_options(options) end function load_options(nlp; options...) - + primary_opt, options = _get_primary_options(options) - + # Initiate interior-point options opt_ipm = MadNLPOptions(nlp; primary_opt...) linear_solver_options = set_options!(opt_ipm, options) diff --git a/test/MOI_interface_test.jl b/test/MOI_interface_test.jl index 750e4fc1..9ffa65be 100644 --- a/test/MOI_interface_test.jl +++ b/test/MOI_interface_test.jl @@ -37,11 +37,6 @@ function test_MOI_Test() ] ); exclude = [ - # TODO: MadNLP does not return the correct multiplier - # when a variable is fixed with MOI.EqualTo (Issue #229). - r"^test_linear_integration$", - "test_quadratic_constraint_GreaterThan", - "test_quadratic_constraint_LessThan", # MadNLP reaches maximum number of iterations instead # of returning infeasibility certificate. r"test_linear_DUAL_INFEASIBLE.*", diff --git a/test/kkt_test.jl b/test/kkt_test.jl index 8502b280..869cc6a3 100644 --- a/test/kkt_test.jl +++ b/test/kkt_test.jl @@ -10,7 +10,6 @@ using LinearAlgebra ind_lb = [2,3,4] ind_ub = [4,5,6] - rhs = KKTVector(VT, n, m, nlb, nub, ind_lb, ind_ub) @test length(rhs) == length(MadNLP.full(rhs)) @test MadNLP.number_primal(rhs) == length(MadNLP.primal(rhs)) == n @@ -33,11 +32,8 @@ end (MadNLP.DenseCondensedKKTSystem, MadNLP.DenseCallback), ] linear_solver = MadNLP.LapackCPUSolver - cnt = MadNLP.MadNLPCounters(; start_time=time()) nlp = MadNLPTests.HS15Model() - ind_cons = MadNLP.get_index_constraints(nlp) - cb = MadNLP.create_callback( Callback, nlp, ) @@ -45,7 +41,6 @@ end kkt = MadNLP.create_kkt_system( KKTSystem, cb, - ind_cons, linear_solver; ) MadNLPTests.test_kkt_system(kkt, cb) diff --git a/test/madnlp_dense.jl b/test/madnlp_dense.jl index 495bad6a..cb24da2a 100644 --- a/test/madnlp_dense.jl +++ b/test/madnlp_dense.jl @@ -18,7 +18,8 @@ function _compare_dense_with_sparse( :callback=>MadNLP.SparseCallback, :inertia_correction_method=>inertia, :linear_solver=>MadNLP.LapackCPUSolver, - :print_level=>MadNLP.ERROR, + :print_level=>MadNLP.INFO, + :dual_initialized=>true, :tol=>tol ) dense_options = Dict{Symbol, Any}( @@ -26,7 +27,8 @@ function _compare_dense_with_sparse( :callback=>MadNLP.DenseCallback, :inertia_correction_method=>inertia, :linear_solver=>MadNLP.LapackCPUSolver, - :print_level=>MadNLP.ERROR, + :print_level=>MadNLP.INFO, + :dual_initialized=>true, :tol=>tol ) @@ -44,10 +46,13 @@ function _compare_dense_with_sparse( @test result_sparse.objective ≈ result_dense.objective atol=atol @test result_sparse.solution ≈ result_dense.solution atol=atol @test result_sparse.multipliers ≈ result_dense.multipliers atol=atol - @test solver.kkt.jac_com[:, 1:n] == solverd.kkt.jac - if isa(solverd.kkt, MadNLP.AbstractReducedKKTSystem) - @test Symmetric(solver.kkt.aug_com, :L) ≈ solverd.kkt.aug_com atol=atol - end + ind_free = setdiff(1:n, ind_fixed) + n_free = length(ind_free) + @test solver.kkt.jac_com[:, 1:n_free] == solverd.kkt.jac[:, ind_free] + # @test solver.kkt.hess_com[1:n_free, 1:n_free] == solverd.kkt.hess[ind_free, ind_free] + # if isa(solverd.kkt, MadNLP.AbstractReducedKKTSystem) + # @test Symmetric(solver.kkt.aug_com[1:n_free, 1:n_free], :L) ≈ solverd.kkt.aug_com[ind_free, ind_free] atol=atol + # end end end @@ -68,7 +73,7 @@ end kkt = solverd.kkt @test isempty(kkt.jac) - @test solverd.kkt.linear_solver.A === kkt.aug_com + @test solverd.kkt.linear_solver.A === kkt.aug_com @test size(kkt.hess) == (n, n) @test length(kkt.pr_diag) == n @test length(kkt.du_diag) == m @@ -113,7 +118,7 @@ end _compare_dense_with_sparse(kkt, n, m, Int[], Int[1, 8]; inertia=MadNLP.InertiaFree) end @testset "Fixed variables" begin - n, m = 10, 5 + n, m = 10, 9 _compare_dense_with_sparse(kkt, n, m, Int[1, 2], Int[]) _compare_dense_with_sparse(kkt, n, m, Int[1, 2], Int[]; inertia=MadNLP.InertiaFree) end diff --git a/test/madnlp_test.jl b/test/madnlp_test.jl index a9863dbc..3852121e 100644 --- a/test/madnlp_test.jl +++ b/test/madnlp_test.jl @@ -60,6 +60,7 @@ testset = [ "SparseUnreducedKKTSystem", ()->MadNLP.Optimizer( kkt_system=MadNLP.SparseUnreducedKKTSystem, + linear_solver=UmfpackSolver, print_level=MadNLP.ERROR), [] ], @@ -111,7 +112,6 @@ if VERSION >= v"1.10" ) end - for (name,optimizer_constructor,exclude) in testset test_madnlp(name,optimizer_constructor,exclude) end @@ -137,6 +137,31 @@ end @test solver.status == MadNLP.SOLVE_SUCCEEDED end +@testset "Fixed variables" begin + nlp = MadNLPTests.HS15Model() + solver = MadNLPSolver(nlp; print_level=MadNLP.ERROR) + MadNLP.solve!(solver) + @test isa(solver.cb.fixed_handler, MadNLP.NoFixedVariables) + + # Fix first variable: + nlp.meta.lvar[1] = 0.5 + solver_sparse = MadNLP.MadNLPSolver(nlp; callback=MadNLP.SparseCallback, print_level=MadNLP.ERROR) + sol_sparse = MadNLP.solve!(solver_sparse) + @test length(solver_sparse.ind_fixed) == 1 + @test isa(solver_sparse.cb.fixed_handler, MadNLP.MakeParameter) + @test solver_sparse.n == 3 # fixed variables are removed + + solver_dense = MadNLP.MadNLPSolver(nlp; callback=MadNLP.DenseCallback, print_level=MadNLP.ERROR) + sol_dense = MadNLP.solve!(solver_dense) + @test length(solver_dense.ind_fixed) == 1 + @test isa(solver_dense.cb.fixed_handler, MadNLP.MakeParameter) + @test solver_dense.n == 4 # fixed variables are frozen + + @test sol_dense.iter == sol_sparse.iter + @test sol_dense.objective == sol_sparse.objective + @test sol_dense.solution == sol_sparse.solution +end + @testset "MadNLP warmstart" begin nlp = MadNLPTests.HS15Model() x0 = NLPModels.get_x0(nlp) From 78e8383961a3b44ac48b43570acaded2ce08f70b Mon Sep 17 00:00:00 2001 From: fpacaud Date: Thu, 3 Oct 2024 12:21:15 +0200 Subject: [PATCH 2/5] fix comparison of multipliers for problems with fixed variables --- test/madnlp_dense.jl | 23 ++++++++++------------- test/madnlp_test.jl | 3 +++ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/test/madnlp_dense.jl b/test/madnlp_dense.jl index cb24da2a..864eeef0 100644 --- a/test/madnlp_dense.jl +++ b/test/madnlp_dense.jl @@ -18,7 +18,7 @@ function _compare_dense_with_sparse( :callback=>MadNLP.SparseCallback, :inertia_correction_method=>inertia, :linear_solver=>MadNLP.LapackCPUSolver, - :print_level=>MadNLP.INFO, + :print_level=>MadNLP.ERROR, :dual_initialized=>true, :tol=>tol ) @@ -27,11 +27,13 @@ function _compare_dense_with_sparse( :callback=>MadNLP.DenseCallback, :inertia_correction_method=>inertia, :linear_solver=>MadNLP.LapackCPUSolver, - :print_level=>MadNLP.INFO, + :print_level=>MadNLP.ERROR, :dual_initialized=>true, :tol=>tol ) + # N.B.: depending on the variables we fix, the problem may not satisfy + # LICQ, implying the multipliers are non-unique at the solution. nlp = MadNLPTests.DenseDummyQP(zeros(T,n), m=m, fixed_variables=ind_fixed, equality_cons=ind_eq) solver = MadNLPSolver(nlp; sparse_options...) @@ -45,14 +47,9 @@ function _compare_dense_with_sparse( @test result_sparse.counters.k == result_dense.counters.k @test result_sparse.objective ≈ result_dense.objective atol=atol @test result_sparse.solution ≈ result_dense.solution atol=atol - @test result_sparse.multipliers ≈ result_dense.multipliers atol=atol ind_free = setdiff(1:n, ind_fixed) n_free = length(ind_free) @test solver.kkt.jac_com[:, 1:n_free] == solverd.kkt.jac[:, ind_free] - # @test solver.kkt.hess_com[1:n_free, 1:n_free] == solverd.kkt.hess[ind_free, ind_free] - # if isa(solverd.kkt, MadNLP.AbstractReducedKKTSystem) - # @test Symmetric(solver.kkt.aug_com[1:n_free, 1:n_free], :L) ≈ solverd.kkt.aug_com[ind_free, ind_free] atol=atol - # end end end @@ -117,11 +114,11 @@ end _compare_dense_with_sparse(kkt, n, m, Int[], Int[1, 8]) _compare_dense_with_sparse(kkt, n, m, Int[], Int[1, 8]; inertia=MadNLP.InertiaFree) end - @testset "Fixed variables" begin - n, m = 10, 9 - _compare_dense_with_sparse(kkt, n, m, Int[1, 2], Int[]) - _compare_dense_with_sparse(kkt, n, m, Int[1, 2], Int[]; inertia=MadNLP.InertiaFree) - end + # @testset "Fixed variables" begin + # n, m = 10, 5 + # _compare_dense_with_sparse(kkt, n, m, Int[1, 2], Int[]) + # _compare_dense_with_sparse(kkt, n, m, Int[1, 2], Int[]; inertia=MadNLP.InertiaFree) + # end end # Now we do not support custom KKT constructor @@ -131,7 +128,7 @@ end # end @testset "MadNLP: restart (PR #113)" begin - n, m = 10, 5 + n, m = 10, 6 nlp = MadNLPTests.DenseDummyQP(zeros(n); m=m) sparse_options = Dict{Symbol, Any}( :kkt_system=>MadNLP.SparseKKTSystem, diff --git a/test/madnlp_test.jl b/test/madnlp_test.jl index 3852121e..ff1d6f93 100644 --- a/test/madnlp_test.jl +++ b/test/madnlp_test.jl @@ -160,6 +160,9 @@ end @test sol_dense.iter == sol_sparse.iter @test sol_dense.objective == sol_sparse.objective @test sol_dense.solution == sol_sparse.solution + @test sol_dense.multipliers == sol_sparse.multipliers + @test sol_dense.multipliers_L == sol_sparse.multipliers_L + @test sol_dense.multipliers_U == sol_sparse.multipliers_U end @testset "MadNLP warmstart" begin From 22e19ff95e9d1ee2a390966a0f9542ecd38daf3a Mon Sep 17 00:00:00 2001 From: fpacaud Date: Thu, 3 Oct 2024 13:56:42 +0200 Subject: [PATCH 3/5] fix generation of documentation --- docs/src/lib/callbacks.md | 10 ++-------- docs/src/man/kkt.md | 2 -- docs/src/man/linear_solvers.md | 2 -- 3 files changed, 2 insertions(+), 12 deletions(-) diff --git a/docs/src/lib/callbacks.md b/docs/src/lib/callbacks.md index 65dbc11e..f3539f94 100644 --- a/docs/src/lib/callbacks.md +++ b/docs/src/lib/callbacks.md @@ -44,11 +44,5 @@ RelaxEquality ``` MadNLP has to keep in memory all the indexes associated to the equality -and inequality constraints. Similarly, MadNLP has to keep track -of the indexes of the bounded variables and the fixed variables. MadNLP -provides a utility [`get_index_constraints`](@ref) to import all the indexes required -by MadNLP. Each index vector is encoded as a `Vector{Int}`. -```@docs -get_index_constraints - -``` +and inequality constraints, as well as the indexes of the bounded variables and the fixed variables. +The indexes are stored explicitly as a `Vector{Int}` in the `AbstractCallback` structure used by MadNLP. diff --git a/docs/src/man/kkt.md b/docs/src/man/kkt.md index 5aac432e..6023846e 100644 --- a/docs/src/man/kkt.md +++ b/docs/src/man/kkt.md @@ -108,7 +108,6 @@ cb = MadNLP.create_callback( MadNLP.SparseCallback, nlp, ) -ind_cons = MadNLP.get_index_constraints(nlp) ``` @@ -130,7 +129,6 @@ the function [`create_kkt_system`](@ref): kkt = MadNLP.create_kkt_system( MadNLP.SparseKKTSystem, cb, - ind_cons, linear_solver, ) diff --git a/docs/src/man/linear_solvers.md b/docs/src/man/linear_solvers.md index eb59f8fb..8380e29c 100644 --- a/docs/src/man/linear_solvers.md +++ b/docs/src/man/linear_solvers.md @@ -13,12 +13,10 @@ cb = MadNLP.create_callback( MadNLP.SparseCallback, nlp, ) -ind_cons = MadNLP.get_index_constraints(nlp) linear_solver = LapackCPUSolver kkt = MadNLP.create_kkt_system( MadNLP.SparseKKTSystem, cb, - ind_cons, linear_solver, ) From 95d682dd0457c8d0a67be115bd9d8920eaba23bd Mon Sep 17 00:00:00 2001 From: fpacaud Date: Thu, 3 Oct 2024 14:08:05 +0200 Subject: [PATCH 4/5] fix wrong multipliers if nlp_scaling=true --- src/IPM/utils.jl | 9 +-------- src/nlpmodels.jl | 5 ++++- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/IPM/utils.jl b/src/IPM/utils.jl index 27399fae..4f749d97 100644 --- a/src/IPM/utils.jl +++ b/src/IPM/utils.jl @@ -49,16 +49,9 @@ end function update!(stats::MadNLPExecutionStats, solver::MadNLPSolver) stats.status = solver.status unpack_x!(stats.solution, solver.cb, variable(solver.x)) + unpack_y!(stats.multipliers, solver.cb, solver.y) unpack_z!(stats.multipliers_L, solver.cb, variable(solver.zl)) unpack_z!(stats.multipliers_U, solver.cb, variable(solver.zu)) - stats.multipliers .= solver.y - # stats.solution .= min.( - # max.( - # @view(primal(solver.x)[1:get_nvar(solver.nlp)]), - # get_lvar(solver.nlp) - # ), - # get_uvar(solver.nlp) - # ) stats.objective = unpack_obj(solver.cb, solver.obj_val) unpack_cons!(stats.constraints, solver.cb, solver.c) stats.constraints .+= solver.rhs diff --git a/src/nlpmodels.jl b/src/nlpmodels.jl index 16c21dde..98c92e79 100644 --- a/src/nlpmodels.jl +++ b/src/nlpmodels.jl @@ -658,7 +658,10 @@ function unpack_x!(x_full, cb::AbstractCallback, x) x_full .= x end function unpack_z!(z_full, cb::AbstractCallback, z) - z_full .= z + z_full .= z ./ cb.obj_scale[] +end +function unpack_y!(y_full, cb::AbstractCallback, y) + y_full .= (y .* cb.con_scale) ./ cb.obj_scale[] end # N.B.: Special getters if we use SparseCallback with a MakeParameter fixed_handler, From 5ac7e48836fade696d248444e92023f8a9171d01 Mon Sep 17 00:00:00 2001 From: fpacaud Date: Thu, 3 Oct 2024 18:07:19 +0200 Subject: [PATCH 5/5] add support for MakeParameter on GPU --- src/nlpmodels.jl | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/src/nlpmodels.jl b/src/nlpmodels.jl index 98c92e79..9fbc95ae 100644 --- a/src/nlpmodels.jl +++ b/src/nlpmodels.jl @@ -296,13 +296,8 @@ function create_sparse_fixed_handler( nnzh = length(ind_hess_free) Hi, Hj = similar(hess_I, nnzh), similar(hess_J, nnzh) - # TODO: vectorize - for k in 1:nnzh - cnt = ind_hess_free[k] - i, j = hess_I[cnt], hess_J[cnt] - Hi[k] = map_full_to_free[i] - Hj[k] = map_full_to_free[j] - end + Hi .= map_full_to_free[hess_I[ind_hess_free]] + Hj .= map_full_to_free[hess_J[ind_hess_free]] resize!(hess_I, nnzh) resize!(hess_J, nnzh) hess_I .= Hi @@ -310,12 +305,8 @@ function create_sparse_fixed_handler( nnzj = length(ind_jac_free) Ji, Jj = similar(jac_I, nnzj), similar(jac_J, nnzj) - for k in 1:nnzj - cnt = ind_jac_free[k] - i, j = jac_I[cnt], jac_J[cnt] - Ji[k] = i - Jj[k] = map_full_to_free[j] - end + Ji .= jac_I[ind_jac_free] + Jj .= map_full_to_free[jac_J[ind_jac_free]] resize!(jac_I, nnzj) resize!(jac_J, nnzj) jac_I .= Ji