diff --git a/src/solution.jl b/src/solution.jl index 5ec64ec87..a257d5864 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -1,123 +1,88 @@ function add_variables!(model, var::AbstractVariable) - var.id_hash == objectid(:constant) && - error("Internal error: constant used as variable") - return if sign(var) == ComplexSign() + if sign(var) == ComplexSign() L = MOI.add_variables(model, length(var)) R = MOI.add_variables(model, length(var)) return (L, R) - else - return MOI.add_variables(model, length(var)) end + return MOI.add_variables(model, length(var)) end """ - solve!(problem::Problem, optimizer_factory; + solve!( + problem::Problem, + optimizer_factory; silent_solver = false, ) -Solves the problem, populating `problem.optval` with the optimal value, -as well as the values of the variables (accessed by [`evaluate`](@ref)) -and constraint duals (accessed by `cons.dual`), where applicable. +Solves the problem, populating `problem.optval` with the optimal value, as well +as the values of the variables (accessed by [`evaluate`](@ref)) and constraint +duals (accessed by `cons.dual`), where applicable. + Optional keyword arguments: -* `silent_solver`: whether the solver should be silent (and not emit output or logs) during the solution process. + + * `silent_solver`: whether the solver should be silent (and not emit output or + logs) during the solution process. """ function solve!(p::Problem, optimizer_factory; silent_solver = false) - context = Context(p, optimizer_factory) - model = context.model - - if silent_solver - MOI.set(model, MOI.Silent(), true) - end - - # check DCPness if problem_vexity(p) in (ConcaveVexity(), NotDcp()) throw(DCPViolationError()) end - + context = Context(p, optimizer_factory) + if silent_solver + MOI.set(context.model, MOI.Silent(), true) + end if context.detected_infeasible_during_formulation[] p.status = MOI.INFEASIBLE else - MOI.optimize!(model) - p.status = MOI.get(model, MOI.TerminationStatus()) + MOI.optimize!(context.model) + p.status = MOI.get(context.model, MOI.TerminationStatus()) end - - p.model = model - + p.model = context.model if p.status != MOI.OPTIMAL @warn "Problem wasn't solved optimally" status = p.status end - - dual_status = MOI.get(model, MOI.DualStatus()) - primal_status = MOI.get(model, MOI.PrimalStatus()) - - var_to_indices = context.var_id_to_moi_indices - id_to_variables = context.id_to_variables - - if primal_status != MOI.NO_SOLUTION - for (id, var_indices) in var_to_indices - var = id_to_variables[id] - vexity(var) == ConstVexity() && continue - if var_indices isa Tuple - vectorized_value_re = - MOI.get(model, MOI.VariablePrimal(), var_indices[1]) - vectorized_value_im = - MOI.get(model, MOI.VariablePrimal(), var_indices[2]) - set_value!( - var, - unpackvec(vectorized_value_re, size(var), false) + - im * unpackvec(vectorized_value_im, size(var), false), - ) - else - vectorized_value = - MOI.get(model, MOI.VariablePrimal(), var_indices) - set_value!( - var, - unpackvec( - vectorized_value, - size(var), - iscomplex(sign(var)), - ), - ) - end - end - else - for (id, var_indices) in var_to_indices - var = id_to_variables[id] - vexity(var) == ConstVexity() && continue + primal_status = MOI.get(context.model, MOI.PrimalStatus()) + for (id, indices) in context.var_id_to_moi_indices + var = context.id_to_variables[id] + if vexity(var) == ConstVexity() + continue # Fixed variable + elseif primal_status == MOI.NO_SOLUTION set_value!(var, nothing) + elseif indices isa Tuple # Complex-valued variable + primal_re = MOI.get(p.model, MOI.VariablePrimal(), indices[1]) + primal_im = MOI.get(p.model, MOI.VariablePrimal(), indices[2]) + set_value!( + var, + _unpack_vector(primal_re, size(var)) + + im * _unpack_vector(primal_im, size(var)), + ) + else + @assert !iscomplex(sign(var)) + primal = MOI.get(p.model, MOI.VariablePrimal(), indices) + set_value!(var, _unpack_vector(primal, size(var))) end end - - if dual_status != MOI.NO_SOLUTION - for (constr, MOI_constr_indices) in pairs(context.constr_to_moi_inds) - populate_dual!(model, constr, MOI_constr_indices) - end - else - # Empty duals - for constr in keys(context.constr_to_moi_inds) - constr.dual = nothing + dual_status = MOI.get(context.model, MOI.DualStatus()) + for (c, indices) in context.constr_to_moi_inds + if dual_status == MOI.NO_SOLUTION + c.dual = nothing + else + populate_dual!(context.model, c, indices) end end - - return nothing + return end -function populate_dual!(model::MOI.ModelLike, constr::Constraint, ::Nothing) - constr.dual = nothing - return nothing +function populate_dual!(::MOI.ModelLike, c::Constraint, ::Nothing) + # If the Constraint does not support `populate_dual!` it must return + # `nothing` from `_add_constraint!`. + c.dual = nothing + return end -# This is type unstable! -function unpackvec(v::AbstractVector, size::Tuple{Int,Int}, iscomplex::Bool) - if iscomplex && length(v) == 2 - return v[1] + im * v[2] - elseif iscomplex - l = length(v) รท 2 - # use views? - return reshape(v[1:l], size) + im * reshape(v[l+1:end], size) - elseif length(v) == 1 +function _unpack_vector(v::AbstractVector, size::Tuple{Int,Int}) + if length(v) == 1 return v[] - else - return reshape(v, size) end + return reshape(v, size) end diff --git a/test/test_utilities.jl b/test/test_utilities.jl index fceaee202..2504568b9 100644 --- a/test/test_utilities.jl +++ b/test/test_utilities.jl @@ -1136,6 +1136,26 @@ function test_tree_interface() return end +function test_multiple_constraint_dual() + x = Variable() + c = x >= 1 + p = minimize(x, [c, c]) + solve!(p, SCS.Optimizer) + @test isapprox(c.dual, 1.0; atol = 1e-5) + return +end + +function test_fixed_variable_value() + x = Variable() + y = Variable() + fix!(x, 2.0) + p = minimize(y, x + y >= 1) + solve!(p, SCS.Optimizer) + @test isapprox(x.value, 2.0; atol = 1e-5) + @test isapprox(y.value, -1.0; atol = 1e-5) + return +end + function test_scalar_fn_constant_objective() x = Variable() p = minimize(2.1, [x >= 1])