Support nonlinear constraints.

gen/generate_wrapper.jl

@@ -8,9 +8,10 @@ scip_headers = vcat(
     filter(h -> startswith(h, "scip_"), all_headers),
     "scipdefplugins.h",
     filter(h -> startswith(h, "cons_"), all_headers),
+    "intervalarith.h",                                 # for nlpi/pub_expr
 )
 lpi_headers = ["type_lpi.h"]
-nlpi_headers = ["type_expr.h", "type_nlpi.h"]
+nlpi_headers = ["type_expr.h", "type_nlpi.h", "type_exprinterpret.h", "pub_expr.h"]
 
 headers = vcat(
     [joinpath(HEADER_BASE, "scip", h) for h in scip_headers],

src/MOI_wrapper.jl

@@ -156,5 +156,6 @@ include(joinpath("MOI_wrapper", "quadratic_constraints.jl"))
 include(joinpath("MOI_wrapper", "soc_constraints.jl"))
 include(joinpath("MOI_wrapper", "sos_constraints.jl"))
 include(joinpath("MOI_wrapper", "abspower_constraints.jl"))
+include(joinpath("MOI_wrapper", "nonlinear_constraints.jl"))
 include(joinpath("MOI_wrapper", "objective.jl"))
 include(joinpath("MOI_wrapper", "results.jl"))

src/MOI_wrapper/nonlinear_constraints.jl

@@ -0,0 +1,22 @@
+# Nonlinear constraints (objective not supported)
+
+MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
+
+function MOI.set(o::Optimizer, ::MOI.NLPBlock, data::MOI.NLPBlockData)
+    # We don't actually store the data (to be queried later during the solve
+    # process). Instead, we extract the expression graphs and add the
+    # corresponding constraints to the model directly.
+    if data.has_objective
+        error("Nonlinear objective not supported by SCIP.jl!")
+    end
+
+    MOI.initialize(data.evaluator, [:ExprGraph])
+    for i in 1:length(data.constraint_bounds)
+        expr = MOI.constraint_expr(data.evaluator, i)
+        bounds = data.constraint_bounds[i]
+        cr = add_nonlinear_constraint(o.mscip, expr, bounds.lower, bounds.upper)
+        # Not registering or returning constraint reference!
+    end
+
+    return nothing
+end

src/MOI_wrapper/variable.jl

@@ -84,18 +84,21 @@ function MOI.delete(o::Optimizer, ci::CI{SVF,S}) where {S <: VAR_TYPES}
     vi = VI(ci.value)
     v = var(o, vi)
 
-    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
-        # Reset bounds from constraint.
-        bounds = o.binbounds[vi]
-        @SC SCIPchgVarLb(o, v, bounds.lower)
-        @SC SCIPchgVarUb(o, v, bounds.upper)
-    end
+    # Reset bounds from constraint if this was a binary, see below.
+    reset_bounds = SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
 
     # don't actually delete any SCIP constraint, just reset type
     infeasible = Ref{SCIP_Bool}()
     @SC SCIPchgVarType(o, var(o, vi), SCIP_VARTYPE_CONTINUOUS, infeasible)
     # TODO: warn if infeasible[] == TRUE?
 
+    # Can only change bounds after chaging the var type.
+    if reset_bounds
+        bounds = o.binbounds[vi]
+        @SC SCIPchgVarLb(o, v, bounds.lower)
+        @SC SCIPchgVarUb(o, v, bounds.upper)
+    end
+
     # but do delete the constraint reference
     delete!(o.constypes[SVF, S], ConsRef(ci.value))
 

src/SCIP.jl

@@ -9,6 +9,9 @@ include("wrapper.jl")
 # memory management
 include("managed_scip.jl")
 
+# constraints from nonlinear expressions
+include("nonlinear.jl")
+
 # implementation of MOI
 include("MOI_wrapper.jl")
 

src/managed_scip.jl

@@ -17,7 +17,7 @@ mutable struct ManagedSCIP
     cons_count::Int64
 
     function ManagedSCIP()
-        scip = Ref{Ptr{SCIP_}}()
+        scip = Ref{Ptr{SCIP_}}(C_NULL)
         @SC SCIPcreate(scip)
         @assert scip[] != C_NULL
         @SC SCIPincludeDefaultPlugins(scip[])
@@ -49,6 +49,12 @@ function free_scip(mscip::ManagedSCIP)
     @assert mscip.scip[] == C_NULL
 end
 
+"Set generic parameter."
+function set_parameter(mscip::ManagedSCIP, name::String, value)
+    @SC SCIPsetParam(mscip, name, Ptr{Cvoid}(value))
+    return nothing
+end
+
 "Return pointer to SCIP variable."
 var(mscip::ManagedSCIP, vr::VarRef) = mscip.vars[vr][]
 
@@ -73,7 +79,7 @@ end
 
 "Add variable to problem (continuous, no bounds), return var ref."
 function add_variable(mscip::ManagedSCIP)
-    var__ = Ref{Ptr{SCIP_VAR}}()
+    var__ = Ref{Ptr{SCIP_VAR}}(C_NULL)
     @SC SCIPcreateVarBasic(mscip, var__, "", -SCIPinfinity(mscip), SCIPinfinity(mscip),
                            0.0, SCIP_VARTYPE_CONTINUOUS)
     @SC SCIPaddVar(mscip, var__[])
@@ -117,7 +123,7 @@ Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
 function add_linear_constraint(mscip::ManagedSCIP, varrefs, coefs, lhs, rhs)
     @assert length(varrefs) == length(coefs)
     vars = [var(mscip, vr) for vr in varrefs]
-    cons__ = Ref{Ptr{SCIP_CONS}}()
+    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
     @SC SCIPcreateConsBasicLinear(
         mscip, cons__, "", length(vars), vars, coefs, lhs, rhs)
     @SC SCIPaddCons(mscip, cons__[])
@@ -148,7 +154,7 @@ function add_quadratic_constraint(mscip::ManagedSCIP, linrefs, lincoefs,
     quadvars1 = [var(mscip, vr) for vr in quadrefs1]
     quadvars2 = [var(mscip, vr) for vr in quadrefs2]
 
-    cons__ = Ref{Ptr{SCIP_CONS}}()
+    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
     @SC SCIPcreateConsBasicQuadratic(
         mscip, cons__, "", length(linvars), linvars, lincoefs,
         length(quadvars1), quadvars1, quadvars2, quadcoefs, lhs, rhs)
@@ -165,7 +171,7 @@ the right-hand side.
 """
 function add_second_order_cone_constraint(mscip::ManagedSCIP, varrefs)
     vars = [var(mscip, vr) for vr in varrefs]
-    cons__ = Ref{Ptr{SCIP_CONS}}()
+    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
     @SC SCIPcreateConsBasicSOC(mscip, cons__, "", length(vars) - 1,
                                vars[2:end], C_NULL, C_NULL, 0.0, vars[1], 1.0, 0.0)
     @SC SCIPaddCons(mscip, cons__[])
@@ -182,7 +188,7 @@ Add special-ordered-set of type 1 to problem, return cons ref.
 function add_special_ordered_set_type1(mscip::ManagedSCIP, varrefs, weights)
     @assert length(varrefs) == length(weights)
     vars = [var(mscip, vr) for vr in varrefs]
-    cons__ = Ref{Ptr{SCIP_CONS}}()
+    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
     @SC SCIPcreateConsBasicSOS1(mscip, cons__, "", length(vars), vars, weights)
     @SC SCIPaddCons(mscip, cons__[])
     return store_cons!(mscip, cons__)
@@ -198,14 +204,14 @@ Add special-ordered-set of type 2 to problem, return cons ref.
 function add_special_ordered_set_type2(mscip::ManagedSCIP, varrefs, weights)
     @assert length(varrefs) == length(weights)
     vars = [var(mscip, vr) for vr in varrefs]
-    cons__ = Ref{Ptr{SCIP_CONS}}()
+    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
     @SC SCIPcreateConsBasicSOS2(mscip, cons__, "", length(vars), vars, weights)
     @SC SCIPaddCons(mscip, cons__[])
     return store_cons!(mscip, cons__)
 end
 
 """
-Add special-ordered-set of type 2 to problem, return cons ref.
+Add abspower constraint to problem, return cons ref.
 
     lhs ≤ sign(x + a) * abs(x + a)^n + c*z ≤ rhs
 
@@ -221,15 +227,9 @@ Add special-ordered-set of type 2 to problem, return cons ref.
 Use `(-)SCIPinfinity(scip)` for one of the bounds if not applicable.
 """
 function add_abspower_constraint(mscip::ManagedSCIP, x, a, n, z, c, lhs, rhs)
-    cons__ = Ref{Ptr{SCIP_CONS}}()
+    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
     @SC SCIPcreateConsBasicAbspower(
         mscip, cons__, "", var(mscip, x), var(mscip, z), n, a, c, lhs, rhs)
     @SC SCIPaddCons(mscip, cons__[])
     return store_cons!(mscip, cons__)
 end
-
-"Set generic parameter."
-function set_parameter(mscip::ManagedSCIP, name::String, value)
-    @SC SCIPsetParam(mscip, name, Ptr{Cvoid}(value))
-    return nothing
-end

src/nonlinear.jl

@@ -0,0 +1,193 @@
+
+# Mapping from Julia (as given by MOI) to SCIP operators
+const OPMAP = Dict{Symbol, SCIP_ExprOp}(
+    :+ => SCIP_EXPR_SUM,
+    :* => SCIP_EXPR_PRODUCT,
+    :- => SCIP_EXPR_MINUS,
+    :/ => SCIP_EXPR_DIV,
+    :^ => SCIP_EXPR_REALPOWER,
+    :sqrt => SCIP_EXPR_SQRT,
+    :exp => SCIP_EXPR_EXP,
+    :log => SCIP_EXPR_LOG,
+)
+
+"""Subexpressions and variables referenced in an expression tree.
+
+Used to convert Julia expression to SCIP expression using recursive calls to the
+mutating push_expr!.
+"""
+mutable struct NonlinExpr
+    vars::Vector{Ptr{SCIP_VAR}}
+end
+
+"Extract operators from Julia expr recursively and convert to SCIP expressions."
+function push_expr!(nonlin::NonlinExpr, mscip::ManagedSCIP, expr::Expr)
+    # Storage for SCIP_EXPR*
+    expr__ = Ref{Ptr{SCIP_EXPR}}(C_NULL)
+    num_children = length(expr.args) - 1
+
+    if Meta.isexpr(expr, :comparison) # range constraint
+        # args: [lhs, <=, mid, <=, rhs], lhs and rhs constant
+        @assert length(expr.args) == 5
+        @assert expr.args[2][1] == expr.args[4][1] == :<=
+
+        # just call on middle expression, bounds are handled outside
+        return push_expr!(nonlin, mscip, expr.args[3])
+
+    elseif Meta.isexpr(expr, :call) # operator
+        op = expr.args[1]
+        if op in [:(==), :<=, :>=]
+            # args: [op, lhs, rhs], rhs constant
+            @assert length(expr.args) == 3
+
+            # just call on lhs expression, bounds are handled outside
+            return push_expr!(nonlin, mscip, expr.args[2])
+
+        elseif op == :^
+            # Special case: power with constant exponent.
+            # The Julia expression considers the base and exponent to be
+            # subexpressions. SCIP does in principle support constant
+            # expressions, but in the case of SCIP_EXPR_REALPOWER, the exponent
+            # value is stored directly, as the second child.
+            @assert num_children == 2
+
+            # Base (first child) is proper sub-expression.
+            base = push_expr!(nonlin, mscip, expr.args[2])
+
+            # Exponent (second child) is stored as value.
+            @assert isa(expr.args[3], Number)
+            exponent = Cdouble(expr.args[3])
+
+            # Create SCIP expression
+            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], base, exponent)
+
+        elseif op == :- && num_children == 1
+            # Special case: unary version of minus. SCIP only supports binary
+            # minus, so we will represent it as :(0 - child).
+
+            # First, insert constant 0 subexpression:
+            left = push_expr!(nonlin, mscip, 0.0)
+
+            # Then, insert the actual subexpression:
+            right = push_expr!(nonlin, mscip, expr.args[2])
+
+            # Finally, add the (binary) minus:
+            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], left, right)
+
+        elseif op in [:sqrt, :exp, :log]
+            # Unary operators
+            @assert num_children == 1
+
+            # Insert child expression:
+            child = push_expr!(nonlin, mscip, expr.args[2])
+
+            # Add this operator on top
+            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], child)
+
+        elseif op in [:-, :/]
+            # Binary operators
+            @assert num_children == 2
+
+            # Create left and right subexpression.
+            left = push_expr!(nonlin, mscip, expr.args[2])
+            right = push_expr!(nonlin, mscip, expr.args[3])
+
+            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], left, right)
+
+        elseif op in [:+, :*]
+            # N-ary operators
+            @assert num_children >= 2
+
+            # Create all children
+            children = Ptr{SCIP_EXPR}[]
+            for subexpr in expr.args[2:end]
+                child = push_expr!(nonlin, mscip, subexpr)
+                push!(children, child)
+            end
+            @assert length(children) == num_children
+
+            @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, OPMAP[op], Cint(num_children), children)
+
+        else
+            error("Operator $op (in $expr) not supported by SCIP.jl!")
+        end
+
+    elseif Meta.isexpr(expr, :ref) # variable
+        # It should look like this:
+        # :(x[MathOptInterface.VariableIndex(1)])
+        @assert expr.args[1] == :x
+        @assert num_children == 1
+        vi = expr.args[2] # MOI.VariableIndex
+        vr = VarRef(vi.value)
+        v = var(mscip, vr)
+
+        # 0-based indexing for SCIP
+        op = SCIP_EXPR_VARIDX
+        @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, op, Cint(length(nonlin.vars)))
+        push!(nonlin.vars, v)
+
+    else
+        error("Expression $expr not supported by SCIP.jl!")
+    end
+
+    # Return this expression to be referenced by parent expressions
+    return expr__[]
+end
+
+function push_expr!(nonlin::NonlinExpr, mscip::ManagedSCIP, expr::Number)
+    # Storage for SCIP_EXPR*
+    expr__ = Ref{Ptr{SCIP_EXPR}}(C_NULL)
+
+    op = SCIP_EXPR_CONST
+    value = Cdouble(expr)
+
+    @SC SCIPexprCreate(SCIPblkmem(mscip), expr__, op, value)
+
+    # double check whether value was saved correctly
+    value_stored = SCIPexprGetOpReal(expr__[])
+    if value != value_stored
+        error("Failed to create SCIP_EXPR / $(op):\n" *
+              "passed in constant value $(value),\n" *
+              "retrieved $(value_stored) instead,\n" *
+              "in expression stored at $(expr__[]).")
+    end
+
+    return expr__[]
+end
+
+
+"""
+Add nonlinear constraint to problem, return cons ref.
+
+    lhs ≤ expression ≤ rhs
+
+# Arguments
+- `expr::Expr`: Julia expression of nonlinear function, as given by MOI.
+- `lhs::Float64`: left-hand side for ranged constraint
+- `rhs::Float64`: right-hand side for ranged constraint
+
+"""
+function add_nonlinear_constraint(mscip::ManagedSCIP, expr::Expr, lhs::Float64, rhs::Float64)
+    nonlin = NonlinExpr([])
+
+    # convert expression recursively, extract root and variable pointers
+    root = push_expr!(nonlin, mscip, expr)
+    vars = nonlin.vars
+
+    # create expression graph object
+    tree__ = Ref{Ptr{SCIP_EXPRTREE}}(C_NULL)
+    @SC SCIPexprtreeCreate(SCIPblkmem(mscip), tree__, root, length(vars), 0, C_NULL)
+    @SC SCIPexprtreeSetVars(tree__[], length(vars), vars)
+
+    # create and add cons_nonlinear
+    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
+    @SC SCIPcreateConsBasicNonlinear(mscip, cons__, "", 0, C_NULL, C_NULL,
+                                     1, tree__, C_NULL, lhs, rhs)
+    @SC SCIPaddCons(mscip, cons__[])
+
+    # free memory
+    @SC SCIPexprtreeFree(tree__)
+
+    # register and return cons ref
+    return store_cons!(mscip, cons__)
+end

src/wrapper.jl

@@ -132,6 +132,10 @@ include(wrap("cons_symresack"))
 include(wrap("cons_varbound"))
 include(wrap("cons_xor"))
 
+# nonlinear expressions
+include(wrap("pub_expr"))
+include(wrap("expr_manual"))
+
 # SCIP_CALL: macro to check return codes, inspired by @assert
 macro SC(ex)
     return :(@assert $(esc(ex)) == SCIP_OKAY)