Tidy src/solution.jl (#582)

Author: odow

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

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])