Xpress NLP solver implementation (#205)

src/MOI/nlp.jl

@@ -0,0 +1,218 @@
+MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
+
+function walk_and_strip_variable_index!(expr::Expr)
+    for i in 1:length(expr.args)
+        if expr.args[i] isa MOI.VariableIndex
+            expr.args[i] = expr.args[i].value
+        end
+        walk_and_strip_variable_index!(expr.args[i])
+    end
+    return
+end
+
+walk_and_strip_variable_index!(not_expr) = nothing
+
+function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
+    if model.nlp_block_data !== nothing
+        model.nlp_block_data = nothing
+        # error("Nonlinear block already set; cannot overwrite. Create a new model instead.")
+    end
+    model.nlp_block_data = nlp_data
+
+    nlp_eval = nlp_data.evaluator
+
+    MOI.initialize(nlp_eval, [:ExprGraph])
+
+    if nlp_data.has_objective
+        # according to test: test_nonlinear_objective_and_moi_objective_test
+        # from MOI 0.10.9, linear objectives are just ignores if the noliena exists
+        # if model.inner.objective_expr !== nothing
+            # error("Two objectives set: One linear, one nonlinear.")
+        # end
+        obj = verify_support(MOI.objective_expr(nlp_eval))
+        walk_and_strip_variable_index!(obj)
+        if obj == :NaN
+            model.objective_expr = 0.0
+            model.termination_status = MOI.INVALID_MODEL
+        else
+            model.objective_expr = obj
+        end
+        model.objective_type = NLP_OBJECTIVE
+    else
+        model.objective_expr = 0.0
+    end
+
+    for i in 1:length(nlp_data.constraint_bounds)
+        expr = verify_support(MOI.constraint_expr(nlp_eval, i))
+        lb = nlp_data.constraint_bounds[i].lower
+        ub = nlp_data.constraint_bounds[i].upper
+        @assert expr.head == :call
+        if expr.args[1] == :(==)
+            @assert lb == ub == expr.args[3]
+        elseif expr.args[1] == :(<=)
+            @assert lb == -Inf
+            lb = nothing
+            @assert ub == expr.args[3]
+        elseif expr.args[1] == :(>=)
+            @assert lb == expr.args[3]
+            @assert ub == Inf
+            ub = nothing
+        else
+            error("Unexpected expression $expr.")
+        end
+        expr = expr.args[2]
+        walk_and_strip_variable_index!(expr)
+        push!(model.nlp_constraint_info, Xpress.NLPConstraintInfo(expr, lb, ub, nothing))
+    end
+    return
+end
+
+# Converting expressions in strings adapted to chgformulastring and chgobjformulastring
+wrap_with_parens(x::String) = string("( ", x, " )")
+
+to_str(x) = string(x)
+
+function to_str(c::Expr)
+    if c.head == :comparison
+        if length(c.args) == 3
+            return join([to_str(c.args[1]), " ", c.args[2], " ", c.args[3]])
+        elseif length(c.args) == 5
+            return join([c.args[1], " ", c.args[2], " ", to_str(c.args[3]), " ",
+                         c.args[4], " ", c.args[5]])
+        else
+            throw(UnrecognizedExpressionException("comparison", c))
+        end
+    elseif c.head == :call
+        if c.args[1] in (:<=,:>=,:(==))
+            if length(c.args) == 3
+                return join([to_str(c.args[2]), " ", to_str(c.args[1]), " ", to_str(c.args[3])])
+            elseif length(c.args) == 5
+                return join([to_str(c.args[1]), " ", to_str(c.args[2]), " ", to_str(c.args[3]), " ", to_str(c.args[4]), " ", to_str(c.args[5])])
+            else
+                throw(UnrecognizedExpressionException("comparison", c))
+            end
+        elseif c.args[1] in (:+,:-,:*,:/)
+            if all(d->isa(d, Real), c.args[2:end]) # handle unary case
+                return wrap_with_parens(string(eval(c)))
+            elseif c.args[1] == :- && length(c.args) == 2
+                return wrap_with_parens(string("( - $(to_str(c.args[2])) )"))
+            else
+                return wrap_with_parens(string(join([to_str(d) for d in c.args[2:end]], join([" ",string(c.args[1]), " "]))))
+            end
+        elseif c.args[1] == :^
+            if length(c.args) != 3
+                throw(UnrecognizedExpressionException("function call", c))
+            end 
+            return wrap_with_parens(join([to_str(c.args[2]), " ",to_str(c.args[1]), " ",to_str(c.args[3])]))
+        elseif c.args[1] in (:exp,:log,:sin,:cos,:abs,:tan,:sqrt)
+            if length(c.args) != 2
+                throw(UnrecognizedExpressionException("function call", c))
+            end
+            return wrap_with_parens(string(join([uppercase(string(c.args[1])), " "]), wrap_with_parens(to_str(c.args[2]))))
+        else
+            throw(UnrecognizedExpressionException("function call", c))
+        end
+    elseif c.head == :ref
+        if c.args[1] == :x
+            idx = c.args[2]
+            @assert isa(idx, Int)
+            # TODO decide is use use defined names
+            # might be messy becaus a use can call his variable "sin"
+            return "x$idx"
+        else
+            throw(UnrecognizedExpressionException("reference", c))
+        end
+    end
+end
+
+verify_support(c) = c
+
+function verify_support(c::Real)
+    if isfinite(c) # blocks NaN and +/-Inf
+        return c
+    end
+    error("Expected number but got $c")
+end
+
+function verify_support(c::Expr)
+    if c.head == :comparison
+        map(verify_support, c.args)
+        return c
+    end
+    if c.head == :call
+        if c.args[1] in (:+, :-, :*, :/, :exp, :log)
+            return c
+        elseif c.args[1] in (:<=, :>=, :(==))
+            map(verify_support, c.args[2:end])
+            return c
+        elseif c.args[1] == :^
+            @assert isa(c.args[2], Real) || isa(c.args[3], Real)
+            return c
+        else # TODO: do automatic transformation for x^y, |x|
+            error("Unsupported expression $c")
+        end
+    end
+    return c
+end
+
+function set_lower_bound(info::NLPConstraintInfo, value::Union{Number, Nothing})
+    if value !== nothing
+        info.lower_bound !== nothing && throw(ArgumentError("Lower bound has already been set"))
+        info.lower_bound = value
+    end
+    return
+end
+
+function set_upper_bound(info::NLPConstraintInfo, value::Union{Number, Nothing})
+    if value !== nothing
+        info.upper_bound !== nothing && throw(ArgumentError("Upper bound has already been set"))
+        info.upper_bound = value
+    end
+    return
+end
+
+function set_bounds(info::NLPConstraintInfo, set::MOI.EqualTo)
+    set_lower_bound(info, set.value)
+    set_upper_bound(info, set.value)
+end
+
+function set_bounds(info::NLPConstraintInfo, set::MOI.GreaterThan)
+    set_lower_bound(info, set.lower)
+end
+
+function set_bounds(info::NLPConstraintInfo, set::MOI.LessThan)
+    set_upper_bound(info, set.upper)
+end
+
+function set_bounds(info::NLPConstraintInfo, set::MOI.Interval)
+    set_lower_bound(info, set.lower)
+    set_upper_bound(info, set.upper)
+end
+
+# Transforming NLconstraints in constraint sets used for affine functions creation in optimize!
+function to_constraint_set(c::Xpress.NLPConstraintInfo)
+    if c.lower_bound !== nothing || c.upper_bound !== nothing
+        if c.upper_bound === nothing
+            return [MOI.GreaterThan(c.lower_bound)]
+        elseif c.lower_bound === nothing
+            return [MOI.LessThan(c.upper_bound)]
+        elseif c.lower_bound==c.upper_bound
+            return [MOI.EqualTo(c.lower_bound)]
+        else
+            return MOI.GreaterThan(c.lower_bound), MOI.LessThan(c.upper_bound)
+        end
+    end
+end
+
+# The problem is Nonlinear if a NLPBlockData has been defined
+function is_nlp(model)
+    return model.nlp_block_data !== nothing
+end
+
+MOI.supports(::Optimizer, ::MOI.NLPBlockDual) = true
+
+function MOI.get(model::Optimizer, attr::MOI.NLPBlockDual)
+    MOI.check_result_index_bounds(model, attr)
+    s = _dual_multiplier(model) 
+    return s .*model.cached_solution.linear_dual[(1:length(model.nlp_constraint_info))]
+end

src/api.jl

@@ -4964,3 +4964,4 @@ end
 function _bo_validate(obranch, p_status)
     @checked Lib.XPRS_bo_validate(obranch, p_status)
 end
+

src/common.jl

@@ -186,7 +186,6 @@ const XPRS_DUALIZE = 8144
 const XPRS_DUALGRADIENT = 8145
 const XPRS_SBITERLIMIT = 8146
 const XPRS_SBBEST = 8147
-const XPRS_MAXCUTTIME = 8149 # kept for compatibility (removed in v41)
 const XPRS_ACTIVESET = 8152 # kept for compatibility (removed in v41)
 const XPRS_BARINDEFLIMIT = 8153
 const XPRS_HEURSTRATEGY = 8154 # kept for compatibility (removed in v41)
@@ -297,7 +296,6 @@ const XPRS_TUNERMODE = 8359
 const XPRS_TUNERMETHOD = 8360
 const XPRS_TUNERTARGET = 8362
 const XPRS_TUNERTHREADS = 8363
-const XPRS_TUNERMAXTIME = 8364 # kept for compatibility (removed in v41)
 const XPRS_TUNERHISTORY = 8365
 const XPRS_TUNERPERMUTE = 8366
 const XPRS_TUNERROOTALG = 8367 # kept for compatibility (removed in v41)
@@ -418,8 +416,8 @@ const XPRS_BRANCHVAR = 1036
 const XPRS_MIPTHREADID = 1037
 const XPRS_ALGORITHM = 1049
 const XPRS_SOLSTATUS = 1053
-const XPRS_TIME = 1122 # kept for compatibility (removed in v41)
 const XPRS_ORIGINALROWS = 1124
+const NLPORIGINALROWS = 1125
 const XPRS_CALLBACKCOUNT_OPTNODE = 1136
 const XPRS_CALLBACKCOUNT_CUTMGR = 1137
 const XPRS_ORIGINALQELEMS = 1157

src/helper.jl

@@ -234,6 +234,7 @@ n_quadratic_elements(prob::XpressProblem) = @_invoke Lib.XPRSgetintattrib(prob,
 n_quadratic_row_coefficients(prob::XpressProblem) = @_invoke Lib.XPRSgetintattrib(prob, Lib.XPRS_ORIGINALQCELEMS, _)::Int
 n_entities(prob::XpressProblem) = @_invoke Lib.XPRSgetintattrib(prob, Lib.XPRS_ORIGINALMIPENTS, _)::Int
 n_setmembers(prob::XpressProblem) = @_invoke Lib.XPRSgetintattrib(prob, Lib.XPRS_ORIGINALSETMEMBERS, _)::Int
+n_nonlinear_coefs(prob::XpressProblem) = @_invoke Lib.XPRSgetintattrib(prob, Lib.XPRS_SLPCOEFFICIENTS, _)::Int
 
 n_original_variables(prob::XpressProblem) = @_invoke Lib.XPRSgetintattrib(prob, Lib.XPRS_ORIGINALCOLS, _)::Int
 n_original_constraints(prob::XpressProblem) = @_invoke Lib.XPRSgetintattrib(prob, Lib.XPRS_ORIGINALROWS, _)::Int
@@ -243,11 +244,17 @@ objective_sense(prob::XpressProblem) = obj_sense(prob)  == Lib.XPRS_OBJ_MINIMIZE
 
 # derived attribute functions
 
+"""
+    n_nonlinear_constraints(prob::XpressProblem)
+Return the number of nonlinear contraints in the XpressProblem
+"""
+n_nonlinear_constraints(prob::XpressProblem) = max(n_nonlinear_coefs(prob) - 1,0)
+
 """
     n_linear_constraints(prob::XpressProblem)
 Return the number of purely linear contraints in the XpressProblem
 """
-n_linear_constraints(prob::XpressProblem) = n_constraints(prob) - n_quadratic_constraints(prob)
+n_linear_constraints(prob::XpressProblem) =n_constraints(prob) - n_quadratic_constraints(prob)-n_nonlinear_constraints(prob)
 
 """
     is_qcp(prob::XpressProblem)
@@ -261,6 +268,12 @@ Return `true` if there are integer entities in the XpressProblem
 """
 is_mixedinteger(prob::XpressProblem) = (n_entities(prob) + n_special_ordered_sets(prob)) > 0
 
+"""
+    is_nonlinear(prob::XpressProblem)
+Return `true` if there are nonlinear strings in the XpressProblem
+"""
+is_nonlinear(prob::XpressProblem) = n_nonlinear_coefs(prob) > 0
+
 """
     is_quadratic_objective(prob::XpressProblem)
 Return `true` if there are quadratic terms in the objective in the XpressProblem
@@ -273,6 +286,7 @@ Return a symbol enconding the type of the problem.]
 Options are: `:LP`, `:QP` and `:QCP`
 """
 function problem_type(prob::XpressProblem)
+    is_nonlinear(prob) ? (:NLP) :
     is_quadratic_constraints(prob) ? (:QCP) :
     is_quadratic_objective(prob)  ? (:QP)  : (:LP)
 end
@@ -293,12 +307,31 @@ function Base.show(io::IO, prob::XpressProblem)
     println(io, "    number of linear constraints           = $(n_linear_constraints(prob))")
     println(io, "    number of quadratic constraints        = $(n_quadratic_constraints(prob))")
     println(io, "    number of sos constraints              = $(n_special_ordered_sets(prob))")
+    println(io, "    number of nonlinear constraints        = $(n_nonlinear_constraints(prob))")
     println(io, "    number of non-zero coeffs              = $(n_non_zero_elements(prob))")
     println(io, "    number of non-zero qp objective terms  = $(n_quadratic_elements(prob))")
     println(io, "    number of non-zero qp constraint terms = $(n_quadratic_row_coefficients(prob))")
     println(io, "    number of integer entities             = $(n_entities(prob))")
 end
 
+const NLPSTATUS_STRING = Dict{Int,String}(
+    Lib.XPRS_NLPSTATUS_UNSTARTED => "0 Unstarted ( XPRS_NLPSTATUS_UNSTARTED).",
+    Lib.XPRS_NLPSTATUS_SOLUTION => "1 Global search incomplete - an integer solution has been found ( XPRS_NLPSTATUS_SOLUTION).",
+    Lib.XPRS_NLPSTATUS_OPTIMAL => "2 Optimal ( XPRS_NLPSTATUS_OPTIMAL).",
+    Lib.XPRS_NLPSTATUS_NOSOLUTION => "3 Global search complete - No solution found ( XPRS_NLPSTATUS_NOSOLUTION).",
+    Lib.XPRS_NLPSTATUS_INFEASIBLE => "4 Infeasible ( XPRS_NLPSTATUS_INFEASIBLE).",
+    Lib.XPRS_NLPSTATUS_UNBOUNDED => "5 Unbounded ( XPRS_NLPSTATUS_UNBOUNDED).",
+    Lib.XPRS_NLPSTATUS_UNFINISHED => "6 Unfinished ( XPRS_NLPSTATUS_UNFINISHED).",
+    Lib.XPRS_NLPSTATUS_UNSOLVED => "7 Problem could not be solved due to numerical issues. ( XPRS_NLPSTATUS_UNSOLVED).",
+)
+
+function nlp_solve_complete(stat)
+    stat in [Lib.XPRS_NLPSTATUS_INFEASIBLE, Lib.XPRS_NLPSTATUS_OPTIMAL]
+end
+function nlp_solve_stopped(stat)
+    stat in [Lib.XPRS_NLPSTATUS_INFEASIBLE, Lib.XPRS_NLPSTATUS_OPTIMAL]
+end
+
 const MIPSTATUS_STRING = Dict{Int,String}(
     Lib.XPRS_MIP_NOT_LOADED => "0 Problem has not been loaded ( XPRS_MIP_NOT_LOADED).",
     Lib.XPRS_MIP_LP_NOT_OPTIMAL => "1 Global search incomplete - the initial continuous relaxation has not been solved and no integer solution has been found ( XPRS_MIP_LP_NOT_OPTIMAL).",

src/lib.jl

@@ -2419,3 +2419,12 @@ end
 function XPRSremovecbgloballog(prob, globallog, data)
     ccall((:XPRSremovecbgloballog, libxprs), Cint, (XPRSprob, Ptr{Cvoid}, Ptr{Cvoid}), prob, globallog, data)
 end
+
+function XSLPcreateprob(prob, _probholder)
+    ccall((:XSLPcreateprob, libxprs), Cint, (XPRSprob, Ptr{XPRSprob},), prob, _probholder)
+end
+
+function XSLPinit()
+    ccall((:XSLPinit, libxprs), Cint, ())
+end
+

src/utils.jl

@@ -144,4 +144,4 @@ function _check(prob, val::Cint)
         throw(XpressError(val, "Xpress internal error:\n\n$e.\n"))
     end
     return nothing
-end
+end

test/MathOptInterface/MOI_callbacks.jl

@@ -7,9 +7,9 @@ const MOI = MathOptInterface
 
 function callback_simple_model()
     model = Xpress.Optimizer(
-        HEURSTRATEGY = 0, # before v41
+        HEURSTRATEGY = 0,
         HEUREMPHASIS = 0,
-        OUTPUTLOG = 0,
+        OUTPUTLOG = 0
     )
 
     MOI.Utilities.loadfromstring!(model, """
@@ -28,7 +28,7 @@ end
 function callback_knapsack_model()
     model = Xpress.Optimizer(
         OUTPUTLOG    = 0,
-        HEURSTRATEGY = 0, # before v41
+        HEURSTRATEGY = 0,
         HEUREMPHASIS = 0,
         CUTSTRATEGY  = 0,
         PRESOLVE     = 0,