From 97d5e2247ab2d464d2a34045bffe8f8c9ef34637 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sun, 11 Apr 2021 15:27:49 +1200 Subject: [PATCH] Rewrite and overhaul of AmplNLWriter. (#125) * WIP: begin a rewrite and overhaul of AmplNLWriter. Ideally, some, or all of this would be migrated into MOI.FileFormats as an NL submodule. This needs many more tests before merging. This implementation is simplier than the current one, because it doesn't try to detect or simplify linear expressions. * Remove MathProgBase * More updates * More fixes * More fixes * More fixes * More fixes * Fix on Julia 1.0 and add comments * More fixes * More fixes * Delete rev_opcode.jl * More fixes * Typos * Add another test * Update MOI_wrapper.jl * Check ExprGraph available * features_available * Fix typo * Change Optimizer to the copy_to interface * Fix MINLPTests * Fix part II --- .solverdata/.gitignore | 2 - Project.toml | 5 +- README.md | 68 +- gen/gen.jl | 15 + src/AmplNLWriter.jl | 1721 ++++++++++++++++------------- src/MOI_wrapper.jl | 531 --------- src/NLExpr.jl | 477 ++++++++ src/nl_convert.jl | 111 -- src/nl_linearity.jl | 263 ----- src/nl_params.jl | 74 -- src/nl_write.jl | 266 ----- src/opcode.jl | 72 ++ test/MINLPTests/run_minlptests.jl | 13 +- test/MOI_wrapper.jl | 55 +- test/NLModel.jl | 692 ++++++++++++ test/nl_convert.jl | 87 -- test/nl_linearity.jl | 24 - test/nl_write.jl | 50 - test/runtests.jl | 15 +- test/sol_file_parser.jl | 32 - 20 files changed, 2281 insertions(+), 2292 deletions(-) delete mode 100644 .solverdata/.gitignore create mode 100644 gen/gen.jl delete mode 100644 src/MOI_wrapper.jl create mode 100644 src/NLExpr.jl delete mode 100644 src/nl_convert.jl delete mode 100644 src/nl_linearity.jl delete mode 100644 src/nl_params.jl delete mode 100644 src/nl_write.jl create mode 100644 src/opcode.jl create mode 100644 test/NLModel.jl delete mode 100644 test/nl_convert.jl delete mode 100644 test/nl_linearity.jl delete mode 100644 test/nl_write.jl delete mode 100644 test/sol_file_parser.jl diff --git a/.solverdata/.gitignore b/.solverdata/.gitignore deleted file mode 100644 index 2b1bf90..0000000 --- a/.solverdata/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -*.nl -*.sol diff --git a/Project.toml b/Project.toml index 6cadb1f..0a801ec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,17 +1,14 @@ name = "AmplNLWriter" uuid = "7c4d4715-977e-5154-bfe0-e096adeac482" -version = "0.6.1" +version = "0.7.0" [deps] MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" -MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] Ipopt = "0.6" Ipopt_jll = "3.13.2" MathOptInterface = "0.9.13" -MathProgBase = "~0.5.0, ~0.6, ~0.7" julia = "1" [extras] diff --git a/README.md b/README.md index df1eb57..3f0ee8d 100644 --- a/README.md +++ b/README.md @@ -1,20 +1,17 @@ # AmplNLWriter.jl -This [Julia](https://github.com/JuliaLang/julia) package is an interface between -[MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) and -[AMPL-enabled solvers](http://ampl.com/products/solvers/all-solvers-for-ampl/). +AmplNLWriter.jl is an interface between [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) +and [AMPL-enabled solvers](http://ampl.com/products/solvers/all-solvers-for-ampl/). -A list of AMPL-enabled solvers is available [here](http://ampl.com/products/solvers/all-solvers-for-ampl/). - -*Development of AmplNLWriter.jl is community driven and has no official connection with the AMPL modeling language or AMPL Optimization Inc.* +*Development of AmplNLWriter.jl is community driven and has no official +connection with the AMPL modeling language or AMPL Optimization Inc.* [![Build Status](https://github.com/jump-dev/AmplNLWriter.jl/workflows/CI/badge.svg?branch=master)](https://github.com/jump-dev/AmplNLWriter.jl/actions?query=workflow%3ACI) [![codecov](https://codecov.io/gh/jump-dev/AmplNLWriter.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/jump-dev/AmplNLWriter.jl) ## Installation -AmplNLWriter.jl can be installed using the Julia package manager with the -following command: +The package can be installed with `Pkg.add`. ```julia import Pkg @@ -23,25 +20,27 @@ Pkg.add("AmplNLWriter") ## Usage -AmplNLWriter.jl provides `AmplNLWriter.Optimizer` as a usable solver in JuMP. -The following Julia code uses the Bonmin solver in JuMP via AmplNLWriter.jl: +To call Ipopt via AmplNLWriter.jl, use: +```julia +using JuMP, AmplNLWriter, Ipopt_jll +model = Model(() -> AmplNLWriter.Optimizer(Ipopt_jll.amplexe)) +``` + +Replace `Ipopt_jll` with `Bonmin_jll` or `Couenne_jll` as appropriate. + +**Note: the `_jll` packages require Julia 1.3 or later.** +You can also pass a string pointing to an AMPL-compatible solver executable. For +example, if the `bonmin` executable is on the system path, it as follows: ```julia using JuMP, AmplNLWriter model = Model(() -> AmplNLWriter.Optimizer("bonmin")) ``` -You can then model and solve your optimization problem as usual. See -[JuMP.jl](https://github.com/jump-dev/JuMP.jl/blob/master/README.md) for more -details. - -The `AmplNLWriter.Optimizer()` constructor requires as the first argument the -name of the solver command needed to run the desired solver. +If the solver is not on the path, the full path to the solver will need to be +passed in. -For example, if the `bonmin` executable is on the system path, you can use -this solver using `AmplNLWriter.Optimizer("bonmin")`. If the solver is not on -the path, the full path to the solver will need to be passed in. This solver -executable must be an AMPL-compatible solver. +## Options The second (optional) argument to `AmplNLWriter.Optimizer()` is a `Vector{String}` of solver options. These options are appended to the solve @@ -50,7 +49,7 @@ you are using. Generally, they will be of the form `"key=value"`, where `key` is the name of the option to set and `value` is the desired value. For example, to set the NLP log level to 0 in Bonmin, you would run -`AmplNLWriter.Optimizer("bonmin", ["bonmin.nlp_log_level=0"])`. +`AmplNLWriter.Optimizer(Bonmin_jll.amplexe, ["bonmin.nlp_log_level=0"])`. For a list of options supported by your solver, check the solver's documentation, or run `/path/to/solver -=` at the command line, e.g., run @@ -78,28 +77,21 @@ A list of available options for the respective `.opt` files can be found here: - [Bonmin](https://github.com/coin-or/Bonmin/blob/master/Bonmin/test/bonmin.opt) (plus Ipopt options) - [Couenne](https://github.com/coin-or/Couenne/blob/master/Couenne/src/couenne.opt) (plus Ipopt and Bonmin options) -## Guides for specific solvers +## SCIP -### Bonmin/Couenne/Ipopt +To use SCIP, you must first compile the `scipampl` binary, which is a version of +SCIP with support for the AMPL .nl interface. -**NOTE: AmplNLWriter v0.4.0 introduced a breaking change by removing -`BonminNLSolver`, `CouenneNLSolver`, and `IpoptNLSolver`. Users are now expected -to pass the path of the solver executable to `AmplNLWriter.Optimizer`.** - -The easiest way to obtain a solver executable for Bonmin, Couenne, or Ipopt is -to download one from [AMPL](https://ampl.com/products/solvers/open-source/). - -### SCIP - -To use SCIP with AmplNLWriter.jl, you must first compile the `scipampl` binary -which is a version of SCIP with support for the AMPL .nl interface. To do this, -you can follow the instructions [here](http://zverovich.net/2012/08/07/using-scip-with-ampl.html), +To do this, you can follow the instructions [here](http://zverovich.net/2012/08/07/using-scip-with-ampl.html), which we have tested on OS X and Linux. After doing this, you can access SCIP through -`AmplNLWriter.Optimizer("/path/to/scipampl")`. Options can be specified for SCIP -using a `scip.set` file, where each line is of the form `key = value`. For -example, the following `scip.set` file will set the verbosity level to 0: +`AmplNLWriter.Optimizer("/path/to/scipampl")`. + +Options can be specified for SCIP using a `scip.set` file, where each line is of +the form `key = value`. + +For example, the following `scip.set` file will set the verbosity level to 0: ``` display/verblevel = 0 ``` diff --git a/gen/gen.jl b/gen/gen.jl new file mode 100644 index 0000000..5593767 --- /dev/null +++ b/gen/gen.jl @@ -0,0 +1,15 @@ +# This script builds the list of recognized ASL opcodes using the header files +# in ASL_jll. Only re-run it if AMPL adds new opcodes (which is unlikely). + +using ASL_jll + +open(joinpath(dirname(@__DIR__), "src", "opcode.jl"), "w") do io + println(""" + # Do not modify. This file is automatically created by the script in `gen.jl`. + """) + filename = joinpath(ASL_jll.artifact_dir, "include", "opcode.hd") + for line in readlines(filename) + items = split(line, c -> c == '\t' || c == ' '; keepempty = false) + println(io, "const ", items[2], " = ", items[3]) + end +end diff --git a/src/AmplNLWriter.jl b/src/AmplNLWriter.jl index 42f7894..452e853 100644 --- a/src/AmplNLWriter.jl +++ b/src/AmplNLWriter.jl @@ -1,596 +1,473 @@ module AmplNLWriter -using MathProgBase -using MathProgBase.SolverInterface -using SparseArrays - -debug = false -setdebug(b::Bool) = global debug = b - -const solverdata_dir = abspath(joinpath(@__DIR__, "..", ".solverdata")) - -include("nl_linearity.jl") -include("nl_params.jl") -include("nl_convert.jl") - -export AmplNLSolver, - BonminNLSolver, - CouenneNLSolver, - IpoptNLSolver, - getsolvername, - getsolveresult, - getsolveresultnum, - getsolvemessage, - getsolveexitcode - -struct AmplNLSolver <: AbstractMathProgSolver - solver_command::String - options::Vector{String} - filename::String -end - -function AmplNLSolver( - solver_command::String, - options::Vector{String} = String[]; - filename::String = "", -) - return AmplNLSolver(solver_command, options, filename) -end - -function BonminNLSolver(options = String[]; filename::String = "") - return error( - """BonminNLSolver is no longer available by default through AmplNLWriter. - - You should install CoinOptServices via - - Pkg.add("CoinOptServices") +import MathOptInterface +const MOI = MathOptInterface - and then replace BonminNLSolver(options=String[]; filename::String="") - with +include("NLExpr.jl") - AmplNLSolver(CoinOptServices.bonmin, options; filename=filename) +### ============================================================================ +### Nonlinear constraints +### ============================================================================ - alternatively, you can call - - AmplNLSolver("path/to/bonmin", options; filename=filename) - """, - ) +struct _NLConstraint + lower::Float64 + upper::Float64 + opcode::Int + expr::_NLExpr end -function CouenneNLSolver(options = String[]; filename::String = "") - return error( - """CouenneNLSolver is no longer available by default through AmplNLWriter. +""" + _NLConstraint(expr::Expr, bound::MOI.NLPBoundsPair) - You should install CoinOptServices via +Convert a constraint in the form of a `expr` into a `_NLConstraint` object. - Pkg.add("CoinOptServices") +See `MOI.constraint_expr` for details on the format. - and then replace CouenneNLSolver(options=String[]; filename::String="") - with +As a validation step, the right-hand side of each constraint must be a constant +term that is given by the `bound`. (If the constraint is an interval constraint, +both the left-hand and right-hand sides must be constants.) - AmplNLSolver(CoinOptServices.couenne, options; filename=filename) +The six NL constraint types are: - alternatively, you can call + l <= g(x) <= u : 0 + g(x) >= l : 1 + g(x) <= u : 2 + g(x) : 3 # We don't support this + g(x) == c : 4 + x ⟂ g(x) : 5 # TODO(odow): Complementarity constraints +""" +function _NLConstraint(expr::Expr, bound::MOI.NLPBoundsPair) + if expr.head == :comparison + @assert length(expr.args) == 5 + if !(expr.args[1] ≈ bound.lower && bound.upper ≈ expr.args[5]) + _warn_invalid_bound(expr, bound) + end + return _NLConstraint( + expr.args[1], + expr.args[5], + 0, + _NLExpr(expr.args[3]), + ) + else + @assert expr.head == :call + @assert length(expr.args) == 3 + if expr.args[1] == :(<=) + if !(-Inf ≈ bound.lower && bound.upper ≈ expr.args[3]) + _warn_invalid_bound(expr, bound) + end + return _NLConstraint(-Inf, expr.args[3], 1, _NLExpr(expr.args[2])) + elseif expr.args[1] == :(>=) + if !(expr.args[3] ≈ bound.lower && bound.upper ≈ Inf) + _warn_invalid_bound(expr, bound) + end + return _NLConstraint(expr.args[3], Inf, 2, _NLExpr(expr.args[2])) + else + @assert expr.args[1] == :(==) + if !(expr.args[3] ≈ bound.lower ≈ bound.upper) + _warn_invalid_bound(expr, bound) + end + return _NLConstraint( + expr.args[3], + expr.args[3], + 4, + _NLExpr(expr.args[2]), + ) + end + end +end - AmplNLSolver("path/to/couenne", options; filename=filename) - """, +function _warn_invalid_bound(expr::Expr, bound::MOI.NLPBoundsPair) + return @warn( + "Invalid bounds detected in nonlinear constraint. Expected " * + "`$(bound.lower) <= g(x) <= $(bound.upper)`, but got the constraint " * + "$(expr)", ) end -function IpoptNLSolver(options = String[]; filename::String = "") - return error( - """IpoptNLSolver is no longer available by default through AmplNLWriter. - - You should install Ipopt via - - Pkg.add("Ipopt") - - and then replace IpoptNLSolver(options=String[]; filename::String="") - with +### ============================================================================ +### Nonlinear models +### ============================================================================ + +@enum(_VariableType, _BINARY, _INTEGER, _CONTINUOUS) + +mutable struct _VariableInfo + # Variable lower bound. + lower::Float64 + # Variable upper bound. + upper::Float64 + # Whether variable is binary or integer. + type::_VariableType + # Primal start of the variable. + start::Union{Float64,Nothing} + # Number of constraints that the variable appears in. + jacobian_count::Int + # If the variable appears in the objective. + in_nonlinear_objective::Bool + # If the objetive appears in a nonlinear constraint. + in_nonlinear_constraint::Bool + # The 0-indexed column of the variable. Computed right at the end. + order::Int + function _VariableInfo(model::MOI.ModelLike, x::MOI.VariableIndex) + start = MOI.get(model, MOI.VariablePrimalStart(), x) + return new(-Inf, Inf, _CONTINUOUS, start, 0, false, false, 0) + end +end - AmplNLSolver(Ipopt.amplexe, options; filename=filename) +struct _NLResults + raw_status_string::String + termination_status::MOI.TerminationStatusCode + primal_status::MOI.ResultStatusCode + objective_value::Float64 + primal_solution::Dict{MOI.VariableIndex,Float64} +end - alternatively, you can call +""" + _solver_command(x::Union{Function,String}) - AmplNLSolver("path/to/ipopt", options; filename=filename) - """, - ) +Functionify the solver command so it can be called as follows: +```julia +foo = _solver_command(x) +foo() do path + run(`\$(path) args...`) end +``` +""" +_solver_command(x::String) = f -> f(x) +_solver_command(x::Function) = x -getsolvername(s::AmplNLSolver) = basename(s.solver_command) - -mutable struct AmplNLMathProgModel <: AbstractMathProgModel +mutable struct Optimizer <: MOI.AbstractOptimizer + optimizer::Function options::Vector{String} + results::_NLResults + # Store MOI.Name(). + name::String + # The objective expression. + f::_NLExpr + sense::MOI.OptimizationSense + # A vector of nonlinear constraints + g::Vector{_NLConstraint} + # A vector of linear constraints + h::Vector{_NLConstraint} + # A dictionary of info for the variables. + x::Dict{MOI.VariableIndex,_VariableInfo} + # A struct to help sort the mess that is variable ordering in NL files. + types::Vector{Vector{MOI.VariableIndex}} + # A vector of the final ordering of the variables. + order::Vector{MOI.VariableIndex} +end - solver_command::String - - x_l::Vector{Float64} - x_u::Vector{Float64} - g_l::Vector{Float64} - g_u::Vector{Float64} - - nvar::Int - ncon::Int - - obj::Any - constrs::Vector{Any} - - lin_constrs::Vector{Dict{Int,Float64}} - lin_obj::Dict{Int,Float64} - - r_codes::Vector{Int} - j_counts::Vector{Int} - - vartypes::Vector{Symbol} - varlinearities_con::Vector{Symbol} - varlinearities_obj::Vector{Symbol} - conlinearities::Vector{Symbol} - objlinearity::Symbol - - v_index_map::Dict{Int,Int} - v_index_map_rev::Dict{Int,Int} - c_index_map::Dict{Int,Int} - c_index_map_rev::Dict{Int,Int} - - sense::Symbol - - x_0::Vector{Float64} +""" + Optimizer( + solver_command::Union{String,Function}, + options::Vector{String} = String[], + ) - file_basename::String - probfile::String - solfile::String +Create a new Optimizer object. - objval::Float64 - solution::Vector{Float64} +`solver_command` should be one of two things: - status::Symbol - solve_exitcode::Int - solve_result_num::Int - solve_result::String - solve_message::String - solve_time::Float64 +* A `String` of the full path of an AMPL-compatible executable +* A function that takes takes a function as input, initializes any environment + as needed, calls the input function with a path to the initialized executable, + and then destructs the environment. - d::AbstractNLPEvaluator +# Examples - function AmplNLMathProgModel( - solver_command::String, - options::Vector{String}, - filename::String, - ) - return new( - options, - solver_command, - zeros(0), - zeros(0), - zeros(0), - zeros(0), - 0, - 0, - :(0), - [], - Dict{Int,Float64}[], - Dict{Int,Float64}(), - Int[], - Int[], - Symbol[], - Symbol[], - Symbol[], - Symbol[], - :Lin, - Dict{Int,Int}(), - Dict{Int,Int}(), - Dict{Int,Int}(), - Dict{Int,Int}(), - :Min, - zeros(0), - filename, - "", - "", - NaN, - zeros(0), - :NotSolved, - -1, - -1, - "?", - "", - NaN, - ) - end -end -mutable struct AmplNLLinearQuadraticModel <: AbstractLinearQuadraticModel - inner::AmplNLMathProgModel -end -mutable struct AmplNLNonlinearModel <: AbstractNonlinearModel - inner::AmplNLMathProgModel -end +A string to an executable: +```julia +Optimizer("/path/to/ipopt.exe", ["print_level=0"]) +``` -include("nl_write.jl") +A function or string provided by a package: +```julia +Optimizer(Ipopt.amplexe, ["print_level=0"]) +# or +Optimizer(Ipopt_jll.amplexe, ["print_level=0"]) +``` -function SolverInterface.NonlinearModel(s::AmplNLSolver) - return AmplNLNonlinearModel( - AmplNLMathProgModel(s.solver_command, s.options, s.filename), - ) +A custom function +```julia +function solver_command(f::Function) + # Create environment ... + ret = f("/path/to/ipopt") + # Destruct environment ... + return ret end -function SolverInterface.LinearQuadraticModel(s::AmplNLSolver) - return AmplNLLinearQuadraticModel( - AmplNLMathProgModel(s.solver_command, s.options, s.filename), +Optimizer(solver_command) +``` +""" +function Optimizer( + solver_command::Union{String,Function} = "", + options::Vector{String} = String[], +) + return Optimizer( + _solver_command(solver_command), + options, + _NLResults( + "Optimize not called.", + MOI.OPTIMIZE_NOT_CALLED, + MOI.NO_SOLUTION, + NaN, + Dict{MOI.VariableIndex,Float64}(), + ), + "", + _NLExpr(false, _NLTerm[], Dict{MOI.VariableIndex,Float64}(), 0.0), + MOI.FEASIBILITY_SENSE, + _NLConstraint[], + _NLConstraint[], + Dict{MOI.VariableIndex,_VariableInfo}(), + [MOI.VariableIndex[] for _ in 1:9], + MOI.VariableIndex[], ) end -function SolverInterface.loadproblem!( - outer::AmplNLNonlinearModel, - nvar::Integer, - ncon::Integer, - x_l, - x_u, - g_l, - g_u, - sense::Symbol, - d::AbstractNLPEvaluator, -) - m = outer.inner - - m.nvar, m.ncon = nvar, ncon - loadcommon!(m, x_l, x_u, g_l, g_u, sense) - - m.d = d - if !(:ExprGraph in features_available(m.d)) - error( - "Unable to solve problem using AmplNLWriter because the " * - "nonlinear evaluator does not supply expression graphs.", - ) - end - initialize(m.d, [:ExprGraph]) - - # Process constraints - m.constrs = map(1:m.ncon) do i - c = constr_expr(m.d, i) +Base.show(io::IO, ::Optimizer) = print(io, "An AMPL (.nl) model") - # Remove relations and bounds from constraint expressions - if length(c.args) == 3 - expected_head = :call - expr_index = 2 - rel_index = 1 +MOI.get(model::Optimizer, ::MOI.SolverName) = "AmplNLWriter" - @assert c.head == expected_head - # Single relation constraint: expr rel bound - rel = c.args[rel_index] - m.r_codes[i] = relation_to_nl[rel] - if rel == [:<=, :(==)] - m.g_u[i] = c.args[3] - end - if rel in [:>=, :(==)] - m.g_l[i] = c.args[3] - end - c = c.args[expr_index] - else - # Double relation constraint: bound <= expr <= bound - @assert c.head == :comparison - m.r_codes[i] = relation_to_nl[:multiple] - m.g_u[i] = c.args[5] - m.g_l[i] = c.args[1] - c = c.args[3] - end +MOI.supports(::Optimizer, ::MOI.Name) = true +MOI.get(model::Optimizer, ::MOI.Name) = model.name +MOI.set(model::Optimizer, ::MOI.Name, name::String) = (model.name = name) - # Convert non-linear expression to non-linear, linear and constant - c, constant, m.conlinearities[i] = - process_expression!(c, m.lin_constrs[i], m.varlinearities_con) - - # Update bounds on constraint - m.g_l[i] -= constant - m.g_u[i] -= constant - - # Update jacobian counts using the linear constraint variables - for j in keys(m.lin_constrs[i]) - m.j_counts[j] += 1 - end - return c +function MOI.empty!(model::Optimizer) + model.results = _NLResults( + "Optimize not called.", + MOI.OPTIMIZE_NOT_CALLED, + MOI.NO_SOLUTION, + NaN, + Dict{MOI.VariableIndex,Float64}(), + ) + model.f = _NLExpr(false, _NLTerm[], Dict{MOI.VariableIndex,Float64}(), 0.0) + empty!(model.g) + empty!(model.h) + empty!(model.x) + for i in 1:9 + empty!(model.types[i]) end + empty!(model.order) + return +end - # Process objective - m.obj = obj_expr(m.d) - if typeof(m.obj) <: Expr - if length(m.obj.args) < 2 - # TODO(odow): what does this case mean? - m.obj = 0 - else - # Convert non-linear expression to non-linear, linear and constant - m.obj, constant, m.objlinearity = - process_expression!(m.obj, m.lin_obj, m.varlinearities_obj) - - # Add constant back into non-linear expression - if constant != 0 - m.obj = add_constant(m.obj, constant) - end - end - end - return m +function MOI.is_empty(model::Optimizer) + return isempty(model.g) && isempty(model.h) && isempty(model.x) end -function SolverInterface.loadproblem!( - outer::AmplNLLinearQuadraticModel, - A::AbstractMatrix, - x_l, - x_u, - c, - g_l, - g_u, - sense, +const _SCALAR_FUNCTIONS = Union{ + MOI.SingleVariable, + MOI.ScalarAffineFunction{Float64}, + MOI.ScalarQuadraticFunction{Float64}, +} + +const _SCALAR_SETS = Union{ + MOI.LessThan{Float64}, + MOI.GreaterThan{Float64}, + MOI.EqualTo{Float64}, + MOI.Interval{Float64}, +} + +function MOI.supports_constraint( + ::Optimizer, + ::Type{<:_SCALAR_FUNCTIONS}, + ::Type{<:_SCALAR_SETS}, ) - m = outer.inner - m.ncon, m.nvar = size(A) - - loadcommon!(m, x_l, x_u, g_l, g_u, sense) - - # Load A into the linear constraints - @assert (m.ncon, m.nvar) == size(A) - load_A!(m, A) - m.constrs = zeros(m.ncon) # Dummy constraint expression trees - - # Load c - for (index, val) in enumerate(c) - m.lin_obj[index] = val - end - m.obj = 0 # Dummy objective expression tree - - # Process variables bounds - for j in 1:m.ncon - lower = m.g_l[j] - upper = m.g_u[j] - if lower == -Inf - if upper == Inf - error("Neither lower nor upper bound on constraint $j") - else - m.r_codes[j] = 1 - end - else - if lower == upper - m.r_codes[j] = 4 - elseif upper == Inf - m.r_codes[j] = 2 - else - m.r_codes[j] = 0 - end - end - end - return m + return true end -function load_A!(m::AmplNLMathProgModel, A::SparseMatrixCSC{Float64}) - for var in 1:A.n, k in A.colptr[var]:(A.colptr[var+1]-1) - m.lin_constrs[A.rowval[k]][var] = A.nzval[k] - m.j_counts[var] += 1 - end +function MOI.supports_constraint( + ::Optimizer, + ::Type{MOI.SingleVariable}, + ::Type{<:Union{MOI.ZeroOne,MOI.Integer}}, +) + return true end -function load_A!(m::AmplNLMathProgModel, A::Matrix{Float64}) - for con in 1:m.ncon, var in 1:m.nvar - val = A[con, var] - if val != 0 - m.lin_constrs[con][var] = val - m.j_counts[var] += 1 - end - end -end +MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true +MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{<:_SCALAR_FUNCTIONS}) = true -function loadcommon!(m::AmplNLMathProgModel, x_l, x_u, g_l, g_u, sense) - @assert m.nvar == length(x_l) == length(x_u) - @assert m.ncon == length(g_l) == length(g_u) +# ============================================================================= - m.x_l, m.x_u = x_l, x_u - m.g_l, m.g_u = g_l, g_u - setsense!(m, sense) +MOI.supports(::Optimizer, ::MOI.NLPBlock) = true - m.lin_constrs = [Dict{Int,Float64}() for _ in 1:m.ncon] - m.j_counts = zeros(Int, m.nvar) +function MOI.set( + model::Optimizer, + ::MOI.NLPBlock, + block::Union{Nothing,MOI.NLPBlockData}, +) + model.nlp_block = block + return +end - m.r_codes = Array{Int}(undef, m.ncon) +MOI.get(model::Optimizer, ::MOI.NLPBlock) = model.nlp_block - m.varlinearities_con = fill(:Lin, m.nvar) - m.varlinearities_obj = fill(:Lin, m.nvar) - m.conlinearities = fill(:Lin, m.ncon) - m.objlinearity = :Lin +# ============================================================================== - m.vartypes = fill(:Cont, m.nvar) - return m.x_0 = zeros(m.nvar) +function MOI.supports( + ::Optimizer, + ::MOI.VariablePrimalStart, + ::Type{MOI.VariableIndex}, +) + return true end -SolverInterface.getvartype(m::AmplNLMathProgModel) = copy(m.vartypes) -function SolverInterface.setvartype!( - m::AmplNLMathProgModel, - cat::Vector{Symbol}, -) - @assert all(x -> (x in [:Cont, :Bin, :Int]), cat) - return m.vartypes = copy(cat) +function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, x, v::Real) + model.x[x].start = Float64(v) + return end -SolverInterface.getsense(m::AmplNLMathProgModel) = m.sense -function SolverInterface.setsense!(m::AmplNLMathProgModel, sense::Symbol) - @assert sense == :Min || sense == :Max - return m.sense = sense +function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, x, ::Nothing) + model.x[x].start = nothing + return end -function SolverInterface.setwarmstart!( - m::AmplNLMathProgModel, - v::Vector{Float64}, +function MOI.get( + model::Optimizer, + ::MOI.VariablePrimalStart, + x::MOI.VariableIndex, ) - return m.x_0 = v -end - -function SolverInterface.optimize!(m::AmplNLMathProgModel) - m.status = :NotSolved - m.solve_exitcode = -1 - m.solve_result_num = -1 - m.solve_result = "?" - m.solve_message = "" - - # There is no non-linear binary type, only non-linear discrete, so make - # sure binary vars have bounds in [0, 1] - for i in 1:m.nvar - if m.vartypes[i] == :Bin - if m.x_l[i] < 0 - m.x_l[i] = 0 - end - if m.x_u[i] > 1 - m.x_u[i] = 1 - end - end - end - - make_var_index!(m) - make_con_index!(m) + return model.x[x].start +end - if length(m.file_basename) == 0 - # No filename specified - write to a randomly named file - file_basepath, f_prob = mktemp(solverdata_dir) - else - file_basepath = joinpath(solverdata_dir, "$(m.file_basename)") - f_prob = open(file_basepath, "w") - end - m.probfile = "$file_basepath.nl" - m.solfile = "$file_basepath.sol" +# ============================================================================== - try - write_nl_file(f_prob, m) - finally - close(f_prob) - end +struct _LinearNLPEvaluator <: MOI.AbstractNLPEvaluator end +MOI.features_available(::_LinearNLPEvaluator) = [:ExprGraph] +MOI.initialize(::_LinearNLPEvaluator, ::Vector{Symbol}) = nothing - # Rename file to have .nl extension (this is required by solvers) - # force flag added to fix issue in Windows, where temp file are not - # absolutely unique and file closing is not fast enough - # See https://github.com/jump-dev/AmplNLWriter.jl/pull/63. - mv(file_basepath, m.probfile, force = true) +MOI.Utilities.supports_default_copy_to(::Optimizer, ::Bool) = false - # Run solver and save exitcode - t = time() - try - cmd = pipeline( - `$(m.solver_command) $(m.probfile) -AMPL $(m.options)`, - stdout = stdout, - stdin = stdin, +function MOI.copy_to( + dest::Optimizer, + model::MOI.ModelLike; + copy_names::Bool = false, +) + mapping = MOI.Utilities.IndexMap() + # Initialize the NLP block. + nlp_block = try + MOI.get(model, MOI.NLPBlock()) + catch + MOI.NLPBlockData(MOI.NLPBoundsPair[], _LinearNLPEvaluator(), false) + end + if !(:ExprGraph in MOI.features_available(nlp_block.evaluator)) + error( + "Unable to use AmplNLWriter because the nonlinear evaluator " * + "does not supply expression graphs.", ) - proc = VERSION < v"0.7-" ? spawn(cmd) : run(cmd, wait = false) - wait(proc) - kill(proc) - m.solve_exitcode = proc.exitcode - catch e - m.solve_exitcode = 1 - @warn("Unable to call solver. Failed with the error: $(e)") - m.solve_result = "$(e)" - end - m.solve_time = time() - t - - if m.solve_exitcode == 0 - read_results(m) + end + MOI.initialize(nlp_block.evaluator, [:ExprGraph]) + # Objective function. + if nlp_block.has_objective + dest.f = _NLExpr(MOI.objective_expr(nlp_block.evaluator)) else - m.status = :Error - m.solution = fill(NaN, m.nvar) - m.solve_result_num = 999 + F = MOI.get(model, MOI.ObjectiveFunctionType()) + obj = MOI.get(model, MOI.ObjectiveFunction{F}()) + dest.f = _NLExpr(obj) end - - # Clean up temp files - if !debug - for temp_file in [m.probfile; m.solfile] - isfile(temp_file) && rm(temp_file) - end + # Nonlinear constraints + for (i, bound) in enumerate(nlp_block.constraint_bounds) + push!( + dest.g, + _NLConstraint(MOI.constraint_expr(nlp_block.evaluator, i), bound), + ) end -end - -function process_expression!( - nonlin_expr::Expr, - lin_expr::Dict{Int,Float64}, - varlinearities::Vector{Symbol}, -) - # Get list of all variables in the expression - extract_variables!(lin_expr, nonlin_expr) - # Extract linear and constant terms from non-linear expression - tree = LinearityExpr(nonlin_expr) - tree = pull_up_constants(tree) - _, tree, constant = prune_linear_terms!(tree, lin_expr) - # Make sure all terms remaining in the tree are .nl-compatible - nonlin_expr = convert_formula(tree) - - # Track which variables appear nonlinearly - nonlin_vars = Dict{Int,Float64}() - extract_variables!(nonlin_vars, nonlin_expr) - for j in keys(nonlin_vars) - varlinearities[j] = :Nonlin - end - - # Remove variables at coeff 0 that aren't also in the nonlinear tree - for (j, coeff) in lin_expr - if coeff == 0 && !(j in keys(nonlin_vars)) - delete!(lin_expr, j) - end + for x in MOI.get(model, MOI.ListOfVariableIndices()) + dest.x[x] = _VariableInfo(model, x) + mapping[x] = x end - - # Mark constraint as nonlinear if anything is left in the tree - linearity = nonlin_expr != 0 ? :Nonlin : :Lin - - return nonlin_expr, constant, linearity -end -function process_expression!(nonlin_expr::Real, lin_expr, varlinearities) - # Special case where body of constraint is constant - # Return empty nonlinear and linear parts, and use the body as the constant - return 0, nonlin_expr, :Lin -end - -SolverInterface.status(m::AmplNLMathProgModel) = m.status -SolverInterface.getsolution(m::AmplNLMathProgModel) = copy(m.solution) -SolverInterface.getobjval(m::AmplNLMathProgModel) = m.objval -SolverInterface.numvar(m::AmplNLMathProgModel) = m.nvar -SolverInterface.numconstr(m::AmplNLMathProgModel) = m.ncon -SolverInterface.getsolvetime(m::AmplNLMathProgModel) = m.solve_time - -# Access to AMPL solve result items -get_solve_result(m::AmplNLMathProgModel) = m.solve_result -get_solve_result_num(m::AmplNLMathProgModel) = m.solve_result_num -get_solve_message(m::AmplNLMathProgModel) = m.solve_message -get_solve_exitcode(m::AmplNLMathProgModel) = m.solve_exitcode - -# We need to track linear coeffs of all variables present in the expression tree -extract_variables!(lin_constr::Dict{Int,Float64}, c) = c -function extract_variables!(lin_constr::Dict{Int,Float64}, c::LinearityExpr) - return extract_variables!(lin_constr, c.c) -end -function extract_variables!(lin_constr::Dict{Int,Float64}, c::Expr) - if c.head == :ref - if c.args[1] == :x - @assert isa(c.args[2], Int) - lin_constr[c.args[2]] = 0 - else - error("Unrecognized reference expression $c") + dest.sense = MOI.get(model, MOI.ObjectiveSense()) + resize!(dest.order, length(dest.x)) + # Now deal with the normal MOI constraints. + for (F, S) in MOI.get(model, MOI.ListOfConstraints()) + _process_constraint(dest, model, F, S, mapping) + end + # Correct bounds of binary variables. Mainly because AMPL doesn't have the + # concept of binary nonlinear variables, but it does have binary linear + # variables! How annoying. + for (x, v) in dest.x + if v.type == _BINARY + v.lower = max(0.0, v.lower) + v.upper = min(1.0, v.upper) end - else - map(arg -> extract_variables!(lin_constr, arg), c.args) end -end - -add_constant(c, constant::Real) = c + constant -add_constant(c::Expr, constant::Real) = Expr(:call, :+, c, constant) - -function make_var_index!(m::AmplNLMathProgModel) + # Jacobian counts. The zero terms for nonlinear constraints should have + # been added when the expression was constructed. + for g in dest.g, v in keys(g.expr.linear_terms) + dest.x[v].jacobian_count += 1 + end + for h in dest.h, v in keys(h.expr.linear_terms) + dest.x[v].jacobian_count += 1 + end + # Now comes the confusing part. + # # AMPL, in all its wisdom, orders variables in a _very_ specific way. # The only hint in "Writing NL files" is the line "Variables are ordered as - # described in Tables 3 and 4 of [5]," which leads us to the following order - # - # 1) Continuous variables that appear in a nonlinear objective AND a nonlinear constraint - # 2) Discrete variables that appear in a nonlinear objective AND a nonlinear constraint - # 3) Continuous variables that appear in a nonlinear constraint, but NOT a nonlinear objective - # 4) Discrete variables that appear in a nonlinear constraint, but NOT a nonlinear objective - # 5) Continuous variables that appear in a nonlinear objective, but NOT a nonlinear constraint - # 6) Discrete variables that appear in a nonlinear objective, but NOT a nonlinear constraint - # 7) Continuous variables that DO NOT appear in a nonlinear objective or a nonlinear constraint - # 8) Binary variables that DO NOT appear in a nonlinear objective or a nonlinear constraint - # 9) Integer variables that DO NOT appear in a nonlinear objective or a nonlinear constraint + # described in Tables 3 and 4 of [5]". # - # Yes, nonlinear variables are broken into continuous/discrete, but linear - # variables are partitioned into continuous, binary, and integer. + # Reading these # # https://cfwebprod.sandia.gov/cfdocs/CompResearch/docs/nlwrite20051130.pdf # https://ampl.com/REFS/hooking2.pdf # + # leads us to the following order + # + # 1) Continuous variables that appear in a + # nonlinear objective AND a nonlinear constraint + # 2) Discrete variables that appear in a + # nonlinear objective AND a nonlinear constraint + # 3) Continuous variables that appear in a + # nonlinear constraint, but NOT a nonlinear objective + # 4) Discrete variables that appear in a + # nonlinear constraint, but NOT a nonlinear objective + # 5) Continuous variables that appear in a + # nonlinear objective, but NOT a nonlinear constraint + # 6) Discrete variables that appear in a + # nonlinear objective, but NOT a nonlinear constraint + # 7) Continuous variables that DO NOT appear in a + # nonlinear objective or a nonlinear constraint + # 8) Binary variables that DO NOT appear in a + # nonlinear objective or a nonlinear constraint + # 9) Integer variables that DO NOT appear in a + # nonlinear objective or a nonlinear constraint + # + # Yes, nonlinear variables are broken into continuous/discrete, but linear + # variables are partitioned into continuous, binary, and integer. (See also, + # the need to modify bounds for binary variables.) + if !dest.f.is_linear + for x in keys(dest.f.linear_terms) + dest.x[x].in_nonlinear_objective = true + end + for x in dest.f.nonlinear_terms + if x isa MOI.VariableIndex + dest.x[x].in_nonlinear_objective = true + end + end + end + for con in dest.g + for x in keys(con.expr.linear_terms) + dest.x[x].in_nonlinear_constraint = true + end + for x in con.expr.nonlinear_terms + if x isa MOI.VariableIndex + dest.x[x].in_nonlinear_constraint = true + end + end + end + types = dest.types + for (x, v) in dest.x + if v.in_nonlinear_constraint && v.in_nonlinear_objective + push!(v.type == _CONTINUOUS ? types[1] : types[2], x) + elseif v.in_nonlinear_constraint + push!(v.type == _CONTINUOUS ? types[3] : types[4], x) + elseif v.in_nonlinear_objective + push!(v.type == _CONTINUOUS ? types[5] : types[6], x) + elseif v.type == _CONTINUOUS + push!(types[7], x) + elseif v.type == _BINARY + push!(types[8], x) + else + @assert v.type == _INTEGER + push!(types[9], x) + end + end # However! Don't let Tables 3 and 4 fool you, because the ordering actually # depends on whether the number of nonlinear variables in the objective only # is _strictly_ greater than the number of nonlinear variables in the @@ -603,332 +480,634 @@ function make_var_index!(m::AmplNLMathProgModel) # all of the first nlvo variables appear nonlinearly in an objective. # # However, even this is slightly incorrect, because I think it should read - # "the first nlvb variables appear nonlinearly." Then, the switch on - # nlvo > nlvc determines whether the next block of variables are the ones - # that appear in the objective only, or the constraints only. + # "For all versions, the first nlvb variables appear nonlinearly." The "nlvo + # - nlvc" part is also clearly incorrect, and should probably read "nlvo - + # nlvb." # - # Here, for example, is the relevant code from Couenne: + # It's a bit confusing, so here is the relevant code from Couenne: # https://github.com/coin-or/Couenne/blob/683c5b305d78a009d59268a4bca01e0ad75ebf02/src/readnl/readnl.cpp#L76-L87 # - # Essentially, what that means is if !(nlvo > nlvc), then swap 3-4 for 5-6 in - # the variable order. - variable_orders = [Int[] for _ in 1:9] - nlvo, nlvc = 0, 0 - for i in 1:m.nvar - if m.varlinearities_obj[i] == :Nonlin && - m.varlinearities_con[i] == :Nonlin - push!(variable_orders[m.vartypes[i] == :Cont ? 1 : 2], i) - elseif m.varlinearities_obj[i] != :Nonlin && - m.varlinearities_con[i] == :Nonlin - push!(variable_orders[m.vartypes[i] == :Cont ? 3 : 4], i) - nlvc += 1 - elseif m.varlinearities_obj[i] == :Nonlin && - m.varlinearities_con[i] != :Nonlin - push!(variable_orders[m.vartypes[i] == :Cont ? 5 : 6], i) - nlvo += 1 - else - if m.vartypes[i] == :Bin - push!(variable_orders[8], i) - elseif m.vartypes[i] == :Int - push!(variable_orders[9], i) + # They interpret this paragraph to mean the switch on nlvo > nlvc determines + # whether the next block of variables are the ones that appear in the + # objective only, or the constraints only. + # + # That makes sense as a design choice, because you can read them in two + # contiguous blocks. + # + # Essentially, what all this means is if !(nlvo > nlvc), then swap 3-4 for + # 5-6 in the variable order. + nlvc = length(types[3]) + length(types[4]) + nlvo = length(types[5]) + length(types[6]) + order_i = if nlvo > nlvc + [1, 2, 3, 4, 5, 6, 7, 8, 9] + else + [1, 2, 5, 6, 3, 4, 7, 8, 9] + end + # Now we can order the variables. + n = 0 + for i in order_i + # Since variables come from a dictionary, there may be differences in + # the order depending on platform and Julia version. Sort by creation + # time for consistency. + for x in sort!(types[i]; by = y -> y.value) + dest.x[x].order = n + dest.order[n+1] = x + n += 1 + end + end + return mapping +end + +_set_to_bounds(set::MOI.Interval) = (0, set.lower, set.upper) +_set_to_bounds(set::MOI.LessThan) = (1, -Inf, set.upper) +_set_to_bounds(set::MOI.GreaterThan) = (2, set.lower, Inf) +_set_to_bounds(set::MOI.EqualTo) = (4, set.value, set.value) + +function _process_constraint(dest::Optimizer, model, F, S, mapping) + for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + f = MOI.get(model, MOI.ConstraintFunction(), ci) + s = MOI.get(model, MOI.ConstraintSet(), ci) + op, l, u = _set_to_bounds(s) + con = _NLConstraint(l, u, op, _NLExpr(f)) + if isempty(con.expr.linear_terms) && isempty(con.expr.nonlinear_terms) + if l <= con.expr.constant <= u + continue else - push!(variable_orders[7], i) + error( + "Malformed constraint. There are no variables and the " * + "bounds don't make sense.", + ) end + elseif con.expr.is_linear + push!(dest.h, con) + mapping[ci] = MOI.ConstraintIndex{F,S}(length(dest.h)) + else + push!(dest.g, con) + mapping[ci] = MOI.ConstraintIndex{F,S}(length(dest.g)) end end - if !(nlvo > nlvc) - # Swap 3 <-> 5 and 4 <-> 6 - variable_orders[3], variable_orders[5] = - variable_orders[5], variable_orders[3] - variable_orders[4], variable_orders[6] = - variable_orders[6], variable_orders[4] + return +end + +function _process_constraint( + dest::Optimizer, + model, + F::Type{MOI.SingleVariable}, + S, + mapping, +) + for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + mapping[ci] = ci + f = MOI.get(model, MOI.ConstraintFunction(), ci) + s = MOI.get(model, MOI.ConstraintSet(), ci) + _, l, u = _set_to_bounds(s) + if l > -Inf + dest.x[f.variable].lower = l + end + if u < Inf + dest.x[f.variable].upper = u + end end + return +end - # Index variables in required order - for var_list in variable_orders - add_to_index_maps!(m.v_index_map, m.v_index_map_rev, var_list) +function _process_constraint( + dest::Optimizer, + model, + F::Type{MOI.SingleVariable}, + S::Union{Type{MOI.ZeroOne},Type{MOI.Integer}}, + mapping, +) + for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + mapping[ci] = ci + f = MOI.get(model, MOI.ConstraintFunction(), ci) + dest.x[f.variable].type = S == MOI.ZeroOne ? _BINARY : _INTEGER end + return end -function make_con_index!(m::AmplNLMathProgModel) - nonlin_cons = Int[] - lin_cons = Int[] +_str(x::Float64) = isinteger(x) ? string(round(Int, x)) : string(x) - for i in 1:m.ncon - if m.conlinearities[i] == :Nonlin - push!(nonlin_cons, i) +_write_term(io, x::Float64, ::Any) = println(io, "n", _str(x)) +_write_term(io, x::Int, ::Any) = println(io, "o", x) +function _write_term(io, x::MOI.VariableIndex, nlmodel) + return println(io, "v", nlmodel.x[x].order) +end + +_is_nary(x::Int) = x in _NARY_OPCODES +_is_nary(x) = false + +function _write_nlexpr(io::IO, expr::_NLExpr, nlmodel::Optimizer) + if expr.is_linear || length(expr.nonlinear_terms) == 0 + # If the expression is linear, just write out the constant term. + _write_term(io, expr.constant, nlmodel) + return + end + # If the constant term is non-zero, we need to write it out. + skip_terms = 0 + if !iszero(expr.constant) + if expr.nonlinear_terms[1] == OPSUMLIST + # The nonlinear expression is a summation. We can write our constant + # first, but we also need to increment the number of arguments by + # one. In addition, since we're writing out the first two terms now, + # we must skip them below. + _write_term(io, OPSUMLIST, nlmodel) + println(io, expr.nonlinear_terms[2] + 1) + _write_term(io, expr.constant, nlmodel) + skip_terms = 2 else - push!(lin_cons, i) + # The nonlinear expression is something other than a summation, so + # add a new + node to the expression. + _write_term(io, OPPLUS, nlmodel) + _write_term(io, expr.constant, nlmodel) end end - for con_list in (nonlin_cons, lin_cons) - add_to_index_maps!(m.c_index_map, m.c_index_map_rev, con_list) + last_nary = false + for term in expr.nonlinear_terms + if skip_terms > 0 + skip_terms -= 1 + continue + end + if last_nary + println(io, term::Int) + last_nary = false + else + _write_term(io, term, nlmodel) + last_nary = _is_nary(term) + end end + return end -function add_to_index_maps!( - forward_map::Dict{Int,Int}, - backward_map::Dict{Int,Int}, - inds::Array{Int}, -) - for i in inds - # Indices are 0-prefixed so the next index is the current dict length - index = length(forward_map) - forward_map[i] = index - backward_map[index] = i +function _write_linear_block(io::IO, expr::_NLExpr, nlmodel::Optimizer) + elements = [(c, nlmodel.x[v].order) for (v, c) in expr.linear_terms] + for (c, x) in sort!(elements; by = i -> i[2]) + println(io, x, " ", _str(c)) end + return end -function read_results(m::AmplNLMathProgModel) - if !isfile(m.solfile) - error( - """Unable to open the solution file. The most likely cause of this - is the solver executable failing unexpectedly. Unfortunately we don't - have any other information about the solution or what went wrong.""", - ) +function Base.write(io::IO, nlmodel::Optimizer) + # ========================================================================== + # Header + # Line 1: Always the same + # Notes: + # * I think these are magic bytes used by AMPL internally for stuff. + println(io, "g3 1 1 0") + + # Line 2: vars, constraints, objectives, ranges, eqns, logical constraints + # Notes: + # * We assume there is always one objective, even if it is just `min 0`. + n_con, n_ranges, n_eqns = 0, 0, 0 + for cons in (nlmodel.g, nlmodel.h), c in cons + n_con += 1 + if c.opcode == 0 + n_ranges += 1 + elseif c.opcode == 4 + n_eqns += 1 + end end - open(m.solfile, "r") do io - return read_results(io, m) + println(io, " $(length(nlmodel.x)) $(n_con) 1 $(n_ranges) $(n_eqns) 0") + + # Line 3: nonlinear constraints, objectives + # Notes: + # * We assume there is always one objective, even if it is just `min 0`. + n_nlcon = length(nlmodel.g) + println(io, " ", n_nlcon, " ", 1) + + # Line 4: network constraints: nonlinear, linear + # Notes: + # * We don't support network constraints. I don't know how they are + # represented. + println(io, " 0 0") + + # Line 5: nonlinear vars in constraints, objectives, both + # Notes: + # * This order is confusingly different to the standard "b, c, o" order. + nlvb = length(nlmodel.types[1]) + length(nlmodel.types[2]) + nlvc = nlvb + length(nlmodel.types[3]) + length(nlmodel.types[4]) + nlvo = nlvb + length(nlmodel.types[5]) + length(nlmodel.types[6]) + println(io, " ", nlvc, " ", nlvo, " ", nlvb) + + # Line 6: linear network variables; functions; arith, flags + # Notes: + # * I don't know what this line means. It is what it is. Apparently `flags` + # is set to 1 to get suffixes in .sol file. + println(io, " 0 0 0 1") + + # Line 7: discrete variables: binary, integer, nonlinear (b,c,o) + # Notes: + # * The order is + # - binary variables in linear only + # - integer variables in linear only + # - binary or integer variables in nonlinear objective and constraint + # - binary or integer variables in nonlinear constraint + # - binary or integer variables in nonlinear objective + nbv = length(nlmodel.types[8]) + niv = length(nlmodel.types[9]) + nl_both = length(nlmodel.types[2]) + nl_cons = length(nlmodel.types[4]) + nl_obj = length(nlmodel.types[6]) + println(io, " ", nbv, " ", niv, " ", nl_both, " ", nl_cons, " ", nl_obj) + + # Line 8: nonzeros in Jacobian, gradients + # Notes: + # * Make sure to include a 0 element for every variable that appears in an + # objective or constraint, even if the linear coefficient is 0. + nnz_jacobian = 0 + for g in nlmodel.g + nnz_jacobian += length(g.expr.linear_terms) end -end - -function read_results(resultio, m::AmplNLMathProgModel) - did_read_solution = read_sol(resultio, m) - - # Convert solve_result - if 0 <= m.solve_result_num < 100 - m.status = :Optimal - m.solve_result = "solved" - elseif 100 <= m.solve_result_num < 200 - # Used to indicate solution present but likely incorrect. - m.status = :Optimal - m.solve_result = "solved?" - warn( - "The solver has returned the status :Optimal, but indicated that there might be an error in the solution. The status code returned by the solver was $(m.solve_result_num). Check the solver documentation for more info.", - ) - elseif 200 <= m.solve_result_num < 300 - m.status = :Infeasible - m.solve_result = "infeasible" - elseif 300 <= m.solve_result_num < 400 - m.status = :Unbounded - m.solve_result = "unbounded" - elseif 400 <= m.solve_result_num < 500 - m.status = :UserLimit - m.solve_result = "limit" - elseif 500 <= m.solve_result_num < 600 - m.status = :Error - m.solve_result = "failure" + for h in nlmodel.h + nnz_jacobian += length(h.expr.linear_terms) end - - # If we didn't get a valid solve_result_num, try to get the status from the - # solve_message string. - # Some solvers (e.g. SCIP) don't ever print the suffixes so we need this. - if m.status == :NotSolved - message = lowercase(m.solve_message) - if occursin("optimal", message) - m.status = :Optimal - elseif occursin("infeasible", message) - m.status = :Infeasible - elseif occursin("unbounded", message) - m.status = :Unbounded - elseif occursin("limit", message) - m.status = :UserLimit - elseif occursin("error", message) - m.status = :Error - end + nnz_gradient = length(nlmodel.f.linear_terms) + println(io, " ", nnz_jacobian, " ", nnz_gradient) + + # Line 9: max name lengths: constraints, variables + # Notes: + # * We don't add names, so this is just 0, 0. + println(io, " 0 0") + + # Line 10: common exprs: b,c,o,c1,o1 + # Notes: + # * We don't add common subexpressions (i.e., V blocks). + # * I assume the notation means + # - b = in nonlinear objective and constraint + # - c = in nonlinear constraint + # - o = in nonlinear objective + # - c1 = in linear constraint + # - o1 = in linear objective + println(io, " 0 0 0 0 0") + # ========================================================================== + # Constraints + # Notes: + # * Nonlinear constraints first, then linear. + # * For linear constraints, write out the constant term here. + for (i, g) in enumerate(nlmodel.g) + println(io, "C", i - 1) + _write_nlexpr(io, g.expr, nlmodel) end - - if did_read_solution - if m.objlinearity == :Nonlin - # Try to use NLPEvaluator if we can. - # Can fail due to unsupported functions so fallback to eval - try - m.objval = eval_f(m.d, m.solution) - return - catch - # do nothing + for (i, h) in enumerate(nlmodel.h) + println(io, "C", i - 1 + n_nlcon) + _write_nlexpr(io, h.expr, nlmodel) + end + # ========================================================================== + # Objective + # Notes: + # * NL files support multiple objectives, but we're just going to write 1, + # so it's always `O0`. + # * For linear objectives, write out the constant term here. + println(io, "O0 ", nlmodel.sense == MOI.MAX_SENSE ? "1" : "0") + _write_nlexpr(io, nlmodel.f, nlmodel) + # ========================================================================== + # VariablePrimalStart + # Notes: + # * Make sure to write out the variables in order. + println(io, "x", length(nlmodel.x)) + for (i, x) in enumerate(nlmodel.order) + start = nlmodel.x[x].start + println(io, i - 1, " ", start === nothing ? 0 : _str(start)) + end + # ========================================================================== + # Constraint bounds + # Notes: + # * Nonlinear constraints go first, then linear. + # * The constant term for linear constraints gets written out in the + # "C" block. + if n_con > 0 + println(io, "r") + # Nonlinear constraints + for g in nlmodel.g + print(io, g.opcode) + if g.opcode == 0 + println(io, " ", _str(g.lower), " ", _str(g.upper)) + elseif g.opcode == 1 + println(io, " ", _str(g.upper)) + elseif g.opcode == 2 + println(io, " ", _str(g.lower)) + else + @assert g.opcode == 4 + println(io, " ", _str(g.lower)) end end - - # Calculate objective value from nonlinear and linear parts - obj_nonlin = eval(substitute_vars!(deepcopy(m.obj), m.solution)) - obj_lin = evaluate_linear(m.lin_obj, m.solution) - m.objval = obj_nonlin + obj_lin + # Linear constraints + for h in nlmodel.h + print(io, h.opcode) + if h.opcode == 0 + println(io, " ", _str(h.lower), " ", _str(h.upper)) + elseif h.opcode == 1 + println(io, " ", _str(h.upper)) + elseif h.opcode == 2 + println(io, " ", _str(h.lower)) + else + @assert h.opcode == 4 + println(io, " ", _str(h.lower)) + end + end + end + # ========================================================================== + # Variable bounds + # Notes: + # * Not much to note, other than to make sure you iterate the variables in + # the correct order. + println(io, "b") + for x in nlmodel.order + v = nlmodel.x[x] + if v.lower == v.upper + println(io, "4 ", _str(v.lower)) + elseif -Inf < v.lower && v.upper < Inf + println(io, "0 ", _str(v.lower), " ", _str(v.upper)) + elseif -Inf == v.lower && v.upper < Inf + println(io, "1 ", _str(v.upper)) + elseif -Inf < v.lower && v.upper == Inf + println(io, "2 ", _str(v.lower)) + else + println(io, "3") + end + end + # ========================================================================== + # Jacobian block + # Notes: + # * If a variable appears in a constraint, it needs to have a corresponding + # entry in the Jacobian block, even if the linear coefficient is zero. + # AMPL uses this to determine the Jacobian sparsity. + # * As before, nonlinear constraints go first, then linear. + # * Don't write out the `k` entry for the last variable, because it can be + # inferred from the total number of elements in the Jacobian as given in + # the header. + if n_con > 0 + println(io, "k", length(nlmodel.x) - 1) + total = 0 + for i in 1:length(nlmodel.order)-1 + total += nlmodel.x[nlmodel.order[i]].jacobian_count + println(io, total) + end + for (i, g) in enumerate(nlmodel.g) + println(io, "J", i - 1, " ", length(g.expr.linear_terms)) + _write_linear_block(io, g.expr, nlmodel) + end + for (i, h) in enumerate(nlmodel.h) + println(io, "J", i - 1 + n_nlcon, " ", length(h.expr.linear_terms)) + _write_linear_block(io, h.expr, nlmodel) + end + end + # ========================================================================== + # Gradient block + # Notes: + # * You only need to write this out if there are linear terms in the + # objective. + if nnz_gradient > 0 + println(io, "G0 ", nnz_gradient) + _write_linear_block(io, nlmodel.f, nlmodel) end + return nlmodel end -function read_sol(m::AmplNLMathProgModel) - if !isfile(m.solfile) - error( - """Unable to open the solution file. The most likely cause of this - is the solver executable failing unexpectedly. Unfortunately we don't - have any other information about the solution or what went wrong.""", - ) +""" + _interpret_status(solve_result_num::Int, raw_status_string::String) + +Convert the `solve_result_num` and `raw_status_string` into MOI-type statuses. + +For the primal status, assume a solution is present. Other code is responsible +for returning `MOI.NO_SOLUTION` if no primal solution is present. +""" +function _interpret_status(solve_result_num::Int, raw_status_string::String) + if 0 <= solve_result_num < 100 + return MOI.LOCALLY_SOLVED, MOI.FEASIBLE_POINT + elseif 100 <= solve_result_num < 200 + return MOI.LOCALLY_SOLVED, MOI.UNKNOWN_RESULT_STATUS + elseif 200 <= solve_result_num < 300 + return MOI.INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS + elseif 300 <= solve_result_num < 400 + return MOI.DUAL_INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS + elseif 400 <= solve_result_num < 500 + return MOI.OTHER_LIMIT, MOI.UNKNOWN_RESULT_STATUS + elseif 500 <= solve_result_num < 600 + return MOI.OTHER_ERROR, MOI.UNKNOWN_RESULT_STATUS end - open(m.solfile, "r") do io - return readsol(io, m) + # If we didn't get a valid solve_result_num, try to get the status from the + # solve_message string. Some solvers (e.g. SCIP) don't ever print the + # suffixes so we need this. + message = lowercase(raw_status_string) + if occursin("optimal", message) + return MOI.LOCALLY_SOLVED, MOI.FEASIBLE_POINT + elseif occursin("infeasible", message) + return MOI.INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS + elseif occursin("unbounded", message) + return MOI.DUAL_INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS + elseif occursin("limit", message) + return MOI.OTHER_LIMIT, MOI.UNKNOWN_RESULT_STATUS + elseif occursin("error", message) + return MOI.OTHER_ERROR, MOI.UNKNOWN_RESULT_STATUS + else + return MOI.OTHER_ERROR, MOI.UNKNOWN_RESULT_STATUS end end -function read_sol(f::IO, m::AmplNLMathProgModel) - # Reference implementation: - # https://github.com/ampl/mp/tree/master/src/asl/solvers/readsol.c - stat = :Undefined - line = "" - - # Extract the solver message by concatenating all of the lines until - # "Options". - while true - line = readline(f) - if length(line) >= 7 && line[1:7] == "Options" - break - else - m.solve_message *= line - end +function _readline(io::IO) + if eof(io) + error("Reached end of sol file unexpectedly.") end + return strip(readline(io)) +end +_readline(io::IO, T) = parse(T, _readline(io)) +function _read_sol(filename::String, model::Optimizer) + return open(io -> _read_sol(io, model), filename, "r") +end + +""" + _read_sol(io::IO, model::Optimizer) + +This function is based on a Julia translation of readsol.c, available at +https://github.com/ampl/asl/blob/64919f75fa7a438f4b41bce892dcbe2ae38343ee/src/solvers/readsol.c +and under the following license: + +Copyright (C) 2017 AMPL Optimization, Inc.; written by David M. Gay. +Permission to use, copy, modify, and distribute this software and its +documentation for any purpose and without fee is hereby granted, +provided that the above copyright notice appear in all copies and that +both that the copyright notice and this permission notice and warranty +disclaimer appear in supporting documentation. + +The author and AMPL Optimization, Inc. disclaim all warranties with +regard to this software, including all implied warranties of +merchantability and fitness. In no event shall the author be liable +for any special, indirect or consequential damages or any damages +whatsoever resulting from loss of use, data or profits, whether in an +action of contract, negligence or other tortious action, arising out +of or in connection with the use or performance of this software. +""" +function _read_sol(io::IO, model::Optimizer) + raw_status_string = "" + line = "" + while !startswith(line, "Options") + line = _readline(io) + raw_status_string *= line + end # Read through all the options. Direct copy of reference implementation. - @assert line[1:7] == "Options" - options = [parse(Int, chomp(readline(f))) for _ in 1:3] + @assert startswith(line, "Options") + options = [_readline(io, Int), _readline(io, Int), _readline(io, Int)] num_options = options[1] - 3 <= num_options <= 9 || + if !(3 <= num_options <= 9) error("expected num_options between 3 and 9; " * "got $num_options") + end need_vbtol = false if options[3] == 3 num_options -= 2 need_vbtol = true end for j in 3:num_options - eof(f) && error() - push!(options, parse(Int, chomp(readline(f)))) + push!(options, _readline(io, Int)) end - # Read number of constraints - num_cons = parse(Int, chomp(readline(f))) - @assert(num_cons == m.ncon) - - # Read number of duals to read in - num_duals_to_read = parse(Int, chomp(readline(f))) - @assert(num_duals_to_read in [0; m.ncon]) - + num_cons = _readline(io, Int) + @assert(num_cons == length(model.g) + length(model.h)) + # Read number of dual solutions to read in + num_duals_to_read = _readline(io, Int) + @assert(num_duals_to_read == 0 || num_duals_to_read == num_cons) # Read number of variables - num_vars = parse(Int, chomp(readline(f))) - @assert(num_vars == m.nvar) - - # Read number of variables to read in - num_vars_to_read = parse(Int, chomp(readline(f))) - @assert(num_vars_to_read in [0; m.nvar]) - + num_vars = _readline(io, Int) + @assert(num_vars == length(model.x)) + # Read number of primal solutions to read in + num_vars_to_read = _readline(io, Int) + @assert(num_vars_to_read == 0 || num_vars_to_read == num_vars) # Skip over vbtol line if present - need_vbtol && readline(f) - - # Skip over duals - # TODO do something with these? - for index in 0:(num_duals_to_read-1) - eof(f) && error("End of file while reading duals.") - line = readline(f) + if need_vbtol + _readline(io) end - - # Next, read for the variable values - x = fill(NaN, m.nvar) - m.objval = NaN - for index in 0:(num_vars_to_read-1) - eof(f) && error("End of file while reading variables.") - line = readline(f) - - i = m.v_index_map_rev[index] - x[i] = parse(Float64, chomp(line)) + # Read dual solutions + # TODO(odow): read in the dual solutions! + for _ in 1:num_duals_to_read + _readline(io) + end + # Read primal solutions + primal_solution = Dict{MOI.VariableIndex,Float64}() + if num_vars_to_read > 0 + for xi in model.order + primal_solution[xi] = _readline(io, Float64) + end end - m.solution = x - # Check for status code - while !eof(f) - line = readline(f) - linevals = split(chomp(line), " ") - num_vals = length(linevals) - if num_vals > 0 && linevals[1] == "objno" - # Check for objno == 0 + solve_result_num = -1 + while !eof(io) + linevals = split(_readline(io), " ") + if length(linevals) > 0 && linevals[1] == "objno" @assert parse(Int, linevals[2]) == 0 - # Get solve_result - m.solve_result_num = parse(Int, linevals[3]) - - # We can stop looking for the 'objno' line + solve_result_num = parse(Int, linevals[3]) break end end - return num_vars_to_read > 0 + termination_status, primal_status = + _interpret_status(solve_result_num, raw_status_string) + objective_value = NaN + if length(primal_solution) > 0 + # TODO(odow): is there a better way of getting this other than + # evaluating it? + objective_value = _evaluate(model.f, primal_solution) + end + return _NLResults( + raw_status_string, + termination_status, + length(primal_solution) > 0 ? primal_status : MOI.NO_SOLUTION, + objective_value, + primal_solution, + ) end -substitute_vars!(c, x::Array{Float64}) = c -function substitute_vars!(c::Expr, x::Array{Float64}) - if c.head == :ref - if c.args[1] == :x - index = c.args[2] - @assert isa(index, Int) - c = x[index] - else - error("Unrecognized reference expression $c") - end - else - if c.head == :call - # Convert .nl unary minus (:neg) back to :- - if c.args[1] == :neg - c.args[1] = :- - # Convert .nl :sum back to :+ - elseif c.args[1] == :sum - c.args[1] = :+ +function MOI.optimize!(model::Optimizer) + temp_dir = mktempdir() + nl_file = joinpath(temp_dir, "model.nl") + open(io -> write(io, model), nl_file, "w") + try + model.optimizer() do solver_path + ret = run( + pipeline( + `$(solver_path) $(nl_file) -AMPL $(model.options)`, + stdout = stdout, + stdin = stdin, + ), + ) + if ret.exitcode != 0 + error("Nonzero exit code: $(ret.exitcode)") end end - map!(arg -> substitute_vars!(arg, x), c.args, c.args) + model.results = _read_sol(joinpath(temp_dir, "model.sol"), model) + catch err + model.results = _NLResults( + "Error calling the solver. Failed with: $(err)", + MOI.OTHER_ERROR, + MOI.NO_SOLUTION, + NaN, + Dict{MOI.VariableIndex,Float64}(), + ) end - return c + return end -function evaluate_linear(linear_coeffs::Dict{Int,Float64}, x::Array{Float64}) - total = 0.0 - for (i, coeff) in linear_coeffs - total += coeff * x[i] - end - return total +function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue) + MOI.check_result_index_bounds(model, attr) + return model.results.objective_value end -# Wrapper functions -for f in [ - :getvartype, - :getsense, - :optimize!, - :status, - :getsolution, - :getobjval, - :numvar, - :numconstr, - :getsolvetime, -] - @eval SolverInterface.$f(m::AmplNLNonlinearModel) = $f(m.inner) - @eval SolverInterface.$f(m::AmplNLLinearQuadraticModel) = $f(m.inner) +function MOI.get( + model::Optimizer, + attr::MOI.VariablePrimal, + x::MOI.VariableIndex, +) + MOI.check_result_index_bounds(model, attr) + return model.results.primal_solution[x] end -for f in [ - :get_solve_result, - :get_solve_result_num, - :get_solve_message, - :get_solve_exitcode, -] - @eval $f(m::AmplNLNonlinearModel) = $f(m.inner) - @eval $f(m::AmplNLLinearQuadraticModel) = $f(m.inner) + +function MOI.get(model::Optimizer, ::MOI.TerminationStatus) + return model.results.termination_status end -for f in [:setvartype!, :setsense!, :setwarmstart!] - @eval SolverInterface.$f(m::AmplNLNonlinearModel, x) = $f(m.inner, x) - @eval SolverInterface.$f(m::AmplNLLinearQuadraticModel, x) = $f(m.inner, x) + +function MOI.get(model::Optimizer, attr::MOI.PrimalStatus) + return attr.N == 1 ? model.results.primal_status : MOI.NO_SOLUTION end -# Utility method for deleting any leftover debug files -function clean_solverdata() - for file in readdir(solverdata_dir) - ext = splitext(file)[2] - (ext == ".nl" || ext == ".sol") && rm(joinpath(solverdata_dir, file)) - end +MOI.get(::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION + +function MOI.get(model::Optimizer, ::MOI.RawStatusString) + return model.results.raw_status_string +end + +function MOI.get(model::Optimizer, ::MOI.ResultCount) + return MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT ? 1 : 0 +end + +function MOI.get( + model::Optimizer, + attr::MOI.ConstraintPrimal, + ci::MOI.ConstraintIndex{<:MOI.SingleVariable}, +) + MOI.check_result_index_bounds(model, attr) + return model.results.primal_solution[MOI.VariableIndex(ci.value)] end -include("MOI_wrapper.jl") +function MOI.get( + model::Optimizer, + attr::MOI.ConstraintPrimal, + ci::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction}, +) + MOI.check_result_index_bounds(model, attr) + return _evaluate(model.h[ci.value].expr, model.results.primal_solution) +end + +function MOI.get( + model::Optimizer, + attr::MOI.ConstraintPrimal, + ci::MOI.ConstraintIndex{<:MOI.ScalarQuadraticFunction}, +) + MOI.check_result_index_bounds(model, attr) + return _evaluate(model.g[ci.value].expr, model.results.primal_solution) +end + +function MOI.write_to_file(model::Optimizer, filename::String) + open(io -> write(io, model), filename, "w") + return +end end diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl deleted file mode 100644 index a52cdba..0000000 --- a/src/MOI_wrapper.jl +++ /dev/null @@ -1,531 +0,0 @@ -import MathOptInterface -const MOI = MathOptInterface -const MOIU = MOI.Utilities - -import MathProgBase -const MPB = MathProgBase.SolverInterface - -MOIU.@model( - Model, - (MOI.ZeroOne, MOI.Integer), - (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval), - (), - (), - (), - (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction), - (), - (), - true # is_optimizer -) - -const MOI_SCALAR_SETS = ( - MOI.EqualTo{Float64}, - MOI.GreaterThan{Float64}, - MOI.LessThan{Float64}, - MOI.Interval{Float64}, -) - -const MOI_SCALAR_FUNCTIONS = - (MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}) - -"Struct to contain the MPB solution." -struct MPBSolution - status::Symbol - raw_status_string::String - objective_value::Float64 - primal_solution::Dict{MOI.VariableIndex,Float64} -end - -""" - Optimizer( - solver_command::String, - options::Vector{String} = String[]; - filename::String = "" - ) - -# Example - - Optimizer(Ipopt.amplexe, ["print_level=0"]) -""" -function Optimizer( - solver_command::String, - options::Vector{String} = String[]; - filename::String = "", -) - model = Model{Float64}() - model.ext[:MPBSolver] = - AmplNLSolver(solver_command, options, filename = filename) - model.ext[:VariablePrimalStart] = Dict{MOI.VariableIndex,Float64}() - return model -end - -function MOI.empty!(model::Model{Float64}) - model.name = "" - model.senseset = false - model.sense = MOI.FEASIBILITY_SENSE - model.objectiveset = false - model.objective = - MOI.ScalarAffineFunction{Float64}(MOI.ScalarAffineTerm{Float64}[], 0.0) - model.num_variables_created = 0 - model.variable_indices = nothing - empty!(model.single_variable_mask) - empty!(model.lower_bound) - empty!(model.upper_bound) - empty!(model.var_to_name) - model.name_to_var = nothing - model.nextconstraintid = 0 - empty!(model.con_to_name) - model.name_to_con = nothing - empty!(model.constrmap) - MOI.empty!(model.moi_scalaraffinefunction) - MOI.empty!(model.moi_scalarquadraticfunction) - solver = model.ext[:MPBSolver] - empty!(model.ext) - model.ext[:MPBSolver] = solver - model.ext[:VariablePrimalStart] = Dict{MOI.VariableIndex,Float64}() - return -end - -Base.show(io::IO, ::Model) = print(io, "An AmplNLWriter model") - -MOI.get(::Model, ::MOI.SolverName) = "AmplNLWriter" - -set_to_bounds(set::MOI.LessThan) = (-Inf, set.upper) -set_to_bounds(set::MOI.GreaterThan) = (set.lower, Inf) -set_to_bounds(set::MOI.EqualTo) = (set.value, set.value) -set_to_bounds(set::MOI.Interval) = (set.lower, set.upper) -set_to_cat(set::MOI.ZeroOne) = :Bin -set_to_cat(set::MOI.Integer) = :Int - -struct NLPEvaluator{T} <: MPB.AbstractNLPEvaluator - inner::T - variable_map::Dict{MOI.VariableIndex,Int} - num_inner_con::Int - has_nlobjective::Bool - objective_expr::Union{Nothing,Expr} - scalar_constraint_expr::Vector{Expr} -end - -""" -MathProgBase expects expressions with variables denoted by `x[i]` for contiguous -`i`. However, JuMP 0.19 creates expressions with `x[MOI.VariableIndex(i)]`. So -we have to recursively walk the expression replacing instances of -MOI.VariableIndex by a corresponding integer. -""" -function replace_variableindex_by_int(variable_map, expr::Expr) - for (i, arg) in enumerate(expr.args) - expr.args[i] = replace_variableindex_by_int(variable_map, arg) - end - return expr -end -function replace_variableindex_by_int(variable_map, expr::MOI.VariableIndex) - return variable_map[expr] -end -replace_variableindex_by_int(variable_map, expr) = expr - -# In the next section we match up the MPB nonlinear functions to the MOI API. -# This is basically just a rename. -function MPB.initialize(d::NLPEvaluator, requested_features::Vector{Symbol}) - if d.inner !== nothing - MOI.initialize(d.inner, requested_features) - end - return -end - -function MPB.features_available(d::NLPEvaluator) - if d.inner !== nothing - return MOI.features_available(d.inner) - else - return [:ExprGraph] - end -end - -function MPB.obj_expr(d::NLPEvaluator) - if d.has_nlobjective - # If d.objective_expr === nothing, then the objective must be nonlinear. - # Query it from the inner NLP evaluator. - expr = MOI.objective_expr(d.inner) - return replace_variableindex_by_int(d.variable_map, expr) - elseif d.objective_expr !== nothing - # d.objective_expr is a SingleVariable, ScalarAffineFunction, or a - # ScalarQuadraticFunction from MOI. If it is unset, it will be `nothing` - # (we enforce this when creating the NLPEvaluator in `optimize!`). - return d.objective_expr - else - return :(0.0) - end -end - -function MPB.constr_expr(d::NLPEvaluator, i) - # There are two types of constraints in the model: - # i = 1, 2, ..., d.num_inner_con are the nonlinear constraints. We access - # these from the inner NLP evaluator. - # i = d.num_inner_con + 1, d.num_inner_con + 2, ..., N are the linear or - # quadratic constraints added by MOI. - if i <= d.num_inner_con - expr = MOI.constraint_expr(d.inner, i) - return replace_variableindex_by_int(d.variable_map, expr) - else - return d.scalar_constraint_expr[i-d.num_inner_con] - end -end - -# In the next section, we need to convert MOI functions and sets into expression -# graphs. First, we convert functions (SingleVariable, ScalarAffine, and -# ScalarQuadratic) into expression graphs. -function func_to_expr_graph(func::MOI.SingleVariable, variable_map) - return Expr(:ref, :x, variable_map[func.variable]) -end - -function func_to_expr_graph(func::MOI.ScalarAffineFunction, variable_map) - expr = Expr(:call, :+, func.constant) - for term in func.terms - push!( - expr.args, - Expr( - :call, - :*, - term.coefficient, - Expr(:ref, :x, variable_map[term.variable_index]), - ), - ) - end - return expr -end - -function func_to_expr_graph(func::MOI.ScalarQuadraticFunction, variable_map) - expr = Expr(:call, :+, func.constant) - for term in func.affine_terms - push!( - expr.args, - Expr( - :call, - :*, - term.coefficient, - Expr(:ref, :x, variable_map[term.variable_index]), - ), - ) - end - for term in func.quadratic_terms - index_1 = variable_map[term.variable_index_1] - index_2 = variable_map[term.variable_index_2] - coef = term.coefficient - # MOI defines quadratic as 1/2 x' Q x :( - if index_1 == index_2 - coef *= 0.5 - end - push!( - expr.args, - Expr( - :call, - :*, - coef, - Expr(:ref, :x, index_1), - Expr(:ref, :x, index_2), - ), - ) - end - return expr -end - -MOI.supports(::Model, ::MOI.NLPBlock) = true - -function MOI.set(model::Model, ::MOI.NLPBlock, block) - model.ext[:NLPBlock] = block - return -end - -MOI.get(model::Model, ::MOI.NLPBlock) = get(model.ext, :NLPBlock, nothing) - -function MOI.supports( - ::Model, - ::MOI.VariablePrimalStart, - ::Type{MOI.VariableIndex}, -) - return true -end - -function MOI.set( - model::Model, - ::MOI.VariablePrimalStart, - variable::MOI.VariableIndex, - value::Union{Nothing,Real}, -) - if value === nothing - delete!(model.ext[:VariablePrimalStart], variable) - else - model.ext[:VariablePrimalStart][variable] = Float64(value) - end - return -end - -function MOI.get( - model::Model, - ::MOI.VariablePrimalStart, - variable::MOI.VariableIndex, -) - return get(model.ext[:VariablePrimalStart], variable, nothing) -end - -# Next, we need to combine a function `Expr` with a MOI set into a comparison. - -function funcset_to_expr_graph(func::Expr, set::MOI.LessThan) - return Expr(:call, :<=, func, set.upper) -end - -function funcset_to_expr_graph(func::Expr, set::MOI.GreaterThan) - return Expr(:call, :>=, func, set.lower) -end - -function funcset_to_expr_graph(func::Expr, set::MOI.EqualTo) - return Expr(:call, :(==), func, set.value) -end - -function funcset_to_expr_graph(func::Expr, set::MOI.Interval) - return Expr(:comparison, set.lower, :<=, func, :<=, set.upper) -end - -# A helper function that converts a MOI function and MOI set into a comparison -# expression. - -function moi_to_expr_graph(func, set, variable_map) - func_expr = func_to_expr_graph(func, variable_map) - return funcset_to_expr_graph(func_expr, set) -end - -# Now we're ready to optimize model. -function MOI.optimize!(model::Model) - # Extract the MathProgBase solver from the model. Recall it's a MOI - # attribute which we're careful not to delete on calls to `empty!`. - mpb_solver = model.ext[:MPBSolver] - - # Get the optimzation sense. - opt_sense = MOI.get(model, MOI.ObjectiveSense()) - sense = opt_sense == MOI.MAX_SENSE ? :Max : :Min - - # Get the NLPBlock from the model. - nlp_block = MOI.get(model, MOI.NLPBlock()) - - # Extract the nonlinear constraint bounds. We're going to append to these - # g_l and g_u vectors later. - num_con = 0 - moi_nlp_evaluator = nlp_block !== nothing ? nlp_block.evaluator : nothing - if nlp_block !== nothing - num_con += length(nlp_block.constraint_bounds) - end - g_l = fill(-Inf, num_con) - g_u = fill(Inf, num_con) - if nlp_block !== nothing - for (i, bound) in enumerate(nlp_block.constraint_bounds) - g_l[i] = bound.lower - g_u[i] = bound.upper - end - end - - # Intialize the variables. We need to form a mapping between the MOI - # VariableIndex and an Int in order to replace instances of - # `x[VariableIndex]` with `x[i]` in the expression graphs. - variables = MOI.get(model, MOI.ListOfVariableIndices()) - num_var = length(variables) - variable_map = Dict{MOI.VariableIndex,Int}() - for (i, variable) in enumerate(variables) - variable_map[variable] = i - end - - # Extract variable bounds. - x_l = fill(-Inf, num_var) - x_u = fill(Inf, num_var) - for set_type in MOI_SCALAR_SETS - for c_ref in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.SingleVariable,set_type}(), - ) - c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref) - c_set = MOI.get(model, MOI.ConstraintSet(), c_ref) - v_index = variable_map[c_func.variable] - lower, upper = set_to_bounds(c_set) - # Note that there might be multiple bounds on the same variable - # (e.g., LessThan and GreaterThan), so we should only update bounds - # if they differ from the defaults. - if lower > -Inf - x_l[v_index] = lower - end - if upper < Inf - x_u[v_index] = upper - end - end - end - - # We have to convert all ScalarAffineFunction-in-Set constraints to an - # expression graph. - scalar_constraint_expr = Expr[] - for func_type in MOI_SCALAR_FUNCTIONS - for set_type in MOI_SCALAR_SETS - for c_ref in MOI.get( - model, - MOI.ListOfConstraintIndices{func_type,set_type}(), - ) - c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref) - c_set = MOI.get(model, MOI.ConstraintSet(), c_ref) - expr = moi_to_expr_graph(c_func, c_set, variable_map) - push!(scalar_constraint_expr, expr) - lower, upper = set_to_bounds(c_set) - push!(g_l, lower) - push!(g_u, upper) - end - end - end - - # MOI objective - obj_type = MOI.get(model, MOI.ObjectiveFunctionType()) - obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_type}()) - obj_func_expr = func_to_expr_graph(obj_func, variable_map) - if obj_func_expr == :(+0.0) - obj_func_expr = nothing - end - - # Build the nlp_evaluator - nlp_evaluator = NLPEvaluator( - moi_nlp_evaluator, - variable_map, - num_con, - nlp_block === nothing ? false : nlp_block.has_objective, - obj_func_expr, - scalar_constraint_expr, - ) - - # Create the MathProgBase model. Note that we pass `num_con` and the number - # of linear constraints. - mpb_model = MPB.NonlinearModel(mpb_solver) - MPB.loadproblem!( - mpb_model, - num_var, - num_con + length(scalar_constraint_expr), - x_l, - x_u, - g_l, - g_u, - sense, - nlp_evaluator, - ) - - # Set any variables to :Bin if they are in ZeroOne and :Int if they are - # Integer. The default is just :Cont. - x_cat = fill(:Cont, num_var) - for set_type in (MOI.ZeroOne, MOI.Integer) - for c_ref in MOI.get( - model, - MOI.ListOfConstraintIndices{MOI.SingleVariable,set_type}(), - ) - c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref) - c_set = MOI.get(model, MOI.ConstraintSet(), c_ref) - v_index = variable_map[c_func.variable] - x_cat[v_index] = set_to_cat(c_set) - end - end - MPB.setvartype!(mpb_model, x_cat) - - # Set the VariablePrimalStart attributes for variables. - variable_primal_start = fill(0.0, num_var) - for (i, variable) in enumerate(variables) - start_val = MOI.get(model, MOI.VariablePrimalStart(), variable) - if start_val !== nothing - variable_primal_start[i] = start_val - end - end - MPB.setwarmstart!(mpb_model, variable_primal_start) - - # Solve the model! - MPB.optimize!(mpb_model) - - # Extract and save the MathProgBase solution. - primal_solution = Dict{MOI.VariableIndex,Float64}() - for (variable, sol) in zip(variables, MPB.getsolution(mpb_model)) - primal_solution[variable] = sol - end - model.ext[:MPBModel] = mpb_model - model.ext[:MPBSolutionAttribute] = MPBSolution( - MPB.status(mpb_model), - mpb_model.inner.solve_result, - MPB.getobjval(mpb_model), - primal_solution, - ) - return -end - -# MOI accessors for the solution info. - -function mpb_solution_attribute(model::Model)::Union{Nothing,MPBSolution} - return get(model.ext, :MPBSolutionAttribute, nothing) -end - -function MOI.get(model::Model, ::MOI.VariablePrimal, var::MOI.VariableIndex) - mpb_solution = mpb_solution_attribute(model) - if mpb_solution === nothing - return - end - return mpb_solution.primal_solution[var] -end - -function MOI.get(model::Model, ::MOI.ConstraintPrimal, idx::MOI.ConstraintIndex) - return MOIU.get_fallback(model, MOI.ConstraintPrimal(), idx) -end - -function MOI.get(model::Model, ::MOI.ObjectiveValue) - mpb_solution = mpb_solution_attribute(model) - if mpb_solution === nothing - return - end - return mpb_solution.objective_value -end - -function MOI.get(model::Model, ::MOI.TerminationStatus) - mpb_solution = mpb_solution_attribute(model) - if mpb_solution === nothing - return MOI.OPTIMIZE_NOT_CALLED - end - status = mpb_solution.status - if status == :Optimal - # TODO(odow): this is not always the case. What if Ipopt solves a - # convex problem? - return MOI.LOCALLY_SOLVED - elseif status == :Infeasible - return MOI.INFEASIBLE - elseif status == :Unbounded - return MOI.DUAL_INFEASIBLE - elseif status == :UserLimit - return MOI.OTHER_LIMIT - elseif status == :Error - return MOI.OTHER_ERROR - end - return MOI.OTHER_ERROR -end - -function MOI.get(model::Model, ::MOI.PrimalStatus) - mpb_solution = mpb_solution_attribute(model) - if mpb_solution === nothing - return MOI.NO_SOLUTION - end - status = mpb_solution.status - if status == :Optimal - return MOI.FEASIBLE_POINT - end - return MOI.NO_SOLUTION -end - -function MOI.get(::Model, ::MOI.DualStatus) - return MOI.NO_SOLUTION -end - -function MOI.get(model::Model, ::MOI.ResultCount) - if MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT - return 1 - else - return 0 - end -end - -function MOI.get(model::Model, ::MOI.RawStatusString) - mpb_solution = mpb_solution_attribute(model) - return mpb_solution.raw_status_string -end diff --git a/src/NLExpr.jl b/src/NLExpr.jl new file mode 100644 index 0000000..47a7361 --- /dev/null +++ b/src/NLExpr.jl @@ -0,0 +1,477 @@ +### ============================================================================ +### Opcodes, and their AMPL <-> Julia conversions. +### ============================================================================ + +include("opcode.jl") + +""" + _JULIA_TO_AMPL + +This dictionary is manualy curated, based on the list of opcodes in `opcode.jl`. + +The goal is to map Julia functions to their AMPL opcode equivalent. + +Sometimes, there is ambiguity, such as the `:+`, which Julia uses for unary, +binary, and n-ary addition, while AMPL doesn't support unary addition, uses +OPPLUS for binary, and OPSUMLIST for n-ary. In these cases, introduce a +different symbol to disambiguate them in the context of this dictionary, and add +logic to `_process_expr!` to rewrite the Julia expression. + +Commented out lines are opcodes supported by AMPL that don't have a clear Julia +equivalent. If you can think of one, feel free to add it. But then go and make +similar changes to `_AMPL_TO_JULIA` and `_NARY_OPCODES`. +""" +const _JULIA_TO_AMPL = Dict{Symbol,Int}( + :+ => OPPLUS, # binary-plus + :- => OPMINUS, + :* => OPMULT, + :/ => OPDIV, + :rem => OPREM, + :^ => OPPOW, + # OPLESS = 6 + :min => MINLIST, # n-ary + :max => MAXLIST, # n-ary + # FLOOR = 13 + # CEIL = 14 + :abs => ABS, + :neg => OPUMINUS, + :|| => OPOR, + :&& => OPAND, + :(<) => LT, + :(<=) => LE, + :(==) => EQ, + :(>=) => GE, + :(>) => GT, + :(!=) => NE, + :(!) => OPNOT, + :ifelse => OPIFnl, + :tanh => OP_tanh, + :tan => OP_tan, + :sqrt => OP_sqrt, + :sinh => OP_sinh, + :sin => OP_sin, + :log10 => OP_log10, + :log => OP_log, + :exp => OP_exp, + :cosh => OP_cosh, + :cos => OP_cos, + :atanh => OP_atanh, + # OP_atan2 = 48, + :atan => OP_atan, + :asinh => OP_asinh, + :asin => OP_asin, + :acosh => OP_acosh, + :acos => OP_acos, + :sum => OPSUMLIST, # n-ary plus + # OPintDIV = 55 + # OPprecision = 56 + # OPround = 57 + # OPtrunc = 58 + # OPCOUNT = 59 + # OPNUMBEROF = 60 + # OPNUMBEROFs = 61 + # OPATLEAST = 62 + # OPATMOST = 63 + # OPPLTERM = 64 + # OPIFSYM = 65 + # OPEXACTLY = 66 + # OPNOTATLEAST = 67 + # OPNOTATMOST = 68 + # OPNOTEXACTLY = 69 + # ANDLIST = 70 + # ORLIST = 71 + # OPIMPELSE = 72 + # OP_IFF = 73 + # OPALLDIFF = 74 + # OPSOMESAME = 75 + # OP1POW = 76 + # OP2POW = 77 + # OPCPOW = 78 + # OPFUNCALL = 79 + # OPNUM = 80 + # OPHOL = 81 + # OPVARVAL = 82 + # N_OPS = 83 +) + +""" + _AMPL_TO_JULIA + +This dictionary is manualy curated, based on the list of supported opcodes +`_JULIA_TO_AMPL`. + +The goal is to map AMPL opcodes to their Julia equivalents. In addition, we need +to know the arity of each opcode. + +If the opcode is an n-ary function, use `-1`. +""" +const _AMPL_TO_JULIA = Dict{Int,Tuple{Int,Function}}( + OPPLUS => (2, +), + OPMINUS => (2, -), + OPMULT => (2, *), + OPDIV => (2, /), + OPREM => (2, rem), + OPPOW => (2, ^), + # OPLESS = 6 + MINLIST => (-1, minimum), + MAXLIST => (-1, maximum), + # FLOOR = 13 + # CEIL = 14 + ABS => (1, abs), + OPUMINUS => (1, -), + OPOR => (2, |), + OPAND => (2, &), + LT => (2, <), + LE => (2, <=), + EQ => (2, ==), + GE => (2, >=), + GT => (2, >), + NE => (2, !=), + OPNOT => (1, !), + OPIFnl => (3, ifelse), + OP_tanh => (1, tanh), + OP_tan => (1, tan), + OP_sqrt => (1, sqrt), + OP_sinh => (1, sinh), + OP_sin => (1, sin), + OP_log10 => (1, log10), + OP_log => (1, log), + OP_exp => (1, exp), + OP_cosh => (1, cosh), + OP_cos => (1, cos), + OP_atanh => (1, atanh), + # OP_atan2 = 48, + OP_atan => (1, atan), + OP_asinh => (1, asinh), + OP_asin => (1, asin), + OP_acosh => (1, acosh), + OP_acos => (1, acos), + OPSUMLIST => (-1, sum), + # OPintDIV = 55 + # OPprecision = 56 + # OPround = 57 + # OPtrunc = 58 + # OPCOUNT = 59 + # OPNUMBEROF = 60 + # OPNUMBEROFs = 61 + # OPATLEAST = 62 + # OPATMOST = 63 + # OPPLTERM = 64 + # OPIFSYM = 65 + # OPEXACTLY = 66 + # OPNOTATLEAST = 67 + # OPNOTATMOST = 68 + # OPNOTEXACTLY = 69 + # ANDLIST = 70 + # ORLIST = 71 + # OPIMPELSE = 72 + # OP_IFF = 73 + # OPALLDIFF = 74 + # OPSOMESAME = 75 + # OP1POW = 76 + # OP2POW = 77 + # OPCPOW = 78 + # OPFUNCALL = 79 + # OPNUM = 80 + # OPHOL = 81 + # OPVARVAL = 82 + # N_OPS = 83 +) + +""" + _NARY_OPCODES + +A manually curated list of n-ary opcodes, taken from Table 8 of "Writing .nl +files." + +Not all of these are implemented. See `_REV_OPCODES` for more detail. +""" +const _NARY_OPCODES = Set([ + MINLIST, + MAXLIST, + OPSUMLIST, + OPCOUNT, + OPNUMBEROF, + OPNUMBEROFs, + ANDLIST, + ORLIST, + OPALLDIFF, +]) + +""" + _UNARY_SPECIAL_CASES + +This dictionary defines a set of unary functions that are special-cased. They +don't exist in the NL file format, but they may be called from Julia, and +they can easily be converted into NL-compatible expressions. + +If you have a new unary-function that you want to support, add it here. +""" +const _UNARY_SPECIAL_CASES = Dict( + :cbrt => (x) -> :($x^(1 / 3)), + :abs2 => (x) -> :($x^2), + :inv => (x) -> :(1 / $x), + :log2 => (x) -> :(log($x) / log(2)), + :log1p => (x) -> :(log(1 + $x)), + :exp2 => (x) -> :(2^$x), + :expm1 => (x) -> :(exp($x) - 1), + :sec => (x) -> :(1 / cos($x)), + :csc => (x) -> :(1 / sin($x)), + :cot => (x) -> :(1 / tan($x)), + :asec => (x) -> :(acos(1 / $x)), + :acsc => (x) -> :(asin(1 / $x)), + :acot => (x) -> :(pi / 2 - atan($x)), + :sind => (x) -> :(sin(pi / 180 * $x)), + :cosd => (x) -> :(cos(pi / 180 * $x)), + :tand => (x) -> :(tan(pi / 180 * $x)), + :secd => (x) -> :(1 / cos(pi / 180 * $x)), + :cscd => (x) -> :(1 / sin(pi / 180 * $x)), + :cotd => (x) -> :(1 / tan(pi / 180 * $x)), + :asind => (x) -> :(asin($x) * 180 / pi), + :acosd => (x) -> :(acos($x) * 180 / pi), + :atand => (x) -> :(atan($x) * 180 / pi), + :asecd => (x) -> :(acos(1 / $x) * 180 / pi), + :acscd => (x) -> :(asin(1 / $x) * 180 / pi), + :acotd => (x) -> :((pi / 2 - atan($x)) * 180 / pi), + :sech => (x) -> :(1 / cosh($x)), + :csch => (x) -> :(1 / sinh($x)), + :coth => (x) -> :(1 / tanh($x)), + :asech => (x) -> :(acosh(1 / $x)), + :acsch => (x) -> :(asinh(1 / $x)), + :acoth => (x) -> :(atanh(1 / $x)), +) + +### ============================================================================ +### Nonlinear expressions +### ============================================================================ + +# TODO(odow): This type isn't great. We should experiment with something that is +# type-stable, like +# +# @enum(_NLType, _INTEGER, _DOUBLE, _VARIABLE) +# struct _NLTerm +# type::_NLType +# data::Int64 +# end +# _NLTerm(x::Int) = _NLTerm(_INTEGER, x) +# _NLTerm(x::Float64) = _NLTerm(_DOUBLE, reinterpret(Int64, x)) +# _NLTerm(x::MOI.VariableIndex) = _NLTerm(_VARIABLE, x.value) +# function _value(x::_NLTerm) +# if x.type == _INTEGER +# return x.data +# elseif x.type == _DOUBLE +# return reinterpret(Float64, x.data) +# else +# @assert x.type == _VARIABLE +# return MOI.VariableIndex(x.data) +# end +# end + +const _NLTerm = Union{Int,Float64,MOI.VariableIndex} + +struct _NLExpr + is_linear::Bool + nonlinear_terms::Vector{_NLTerm} + linear_terms::Dict{MOI.VariableIndex,Float64} + constant::Float64 +end + +function Base.:(==)(x::_NLExpr, y::_NLExpr) + return x.is_linear == y.is_linear && + x.nonlinear_terms == y.nonlinear_terms && + x.linear_terms == y.linear_terms && + x.constant == y.constant +end + +_NLExpr(x::MOI.VariableIndex) = _NLExpr(true, _NLTerm[], Dict(x => 1.0), 0.0) + +_NLExpr(x::MOI.SingleVariable) = _NLExpr(x.variable) + +function _add_or_set(dict, key, value) + if haskey(dict, key) + dict[key] += value + else + dict[key] = value + end + return +end + +function _NLExpr(x::MOI.ScalarAffineFunction) + linear = Dict{MOI.VariableIndex,Float64}() + for (i, term) in enumerate(x.terms) + _add_or_set(linear, term.variable_index, term.coefficient) + end + return _NLExpr(true, _NLTerm[], linear, x.constant) +end + +function _NLExpr(x::MOI.ScalarQuadraticFunction) + linear = Dict{MOI.VariableIndex,Float64}() + for (i, term) in enumerate(x.affine_terms) + _add_or_set(linear, term.variable_index, term.coefficient) + end + terms = _NLTerm[] + N = length(x.quadratic_terms) + if N == 0 || N == 1 + # If there are 0 or 1 terms, no need for an addition node. + elseif N == 2 + # If there are two terms, use binary addition. + push!(terms, OPPLUS) + elseif N > 2 + # If there are more, use n-ary addition. + push!(terms, OPSUMLIST) + push!(terms, N) + end + for term in x.quadratic_terms + coefficient = term.coefficient + # MOI defines quadratic as 1/2 x' Q x :( + if term.variable_index_1 == term.variable_index_2 + coefficient *= 0.5 + end + # Optimization: no need for the OPMULT if the coefficient is 1. + if !isone(coefficient) + push!(terms, OPMULT) + push!(terms, coefficient) + end + push!(terms, OPMULT) + push!(terms, term.variable_index_1) + push!(terms, term.variable_index_2) + # For the Jacobian sparsity patterns, we need to add a linear term, even + # if the variable only appears nonlinearly. + _add_or_set(linear, term.variable_index_1, 0.0) + _add_or_set(linear, term.variable_index_2, 0.0) + end + return _NLExpr(false, terms, linear, x.constant) +end + +function _NLExpr(expr::Expr) + nlexpr = _NLExpr(false, _NLTerm[], Dict{MOI.VariableIndex,Float64}(), 0.0) + _process_expr!(nlexpr, expr) + return nlexpr +end + +function _process_expr!(expr::_NLExpr, arg::Real) + return push!(expr.nonlinear_terms, Float64(arg)) +end + +function _process_expr!(expr::_NLExpr, arg::MOI.VariableIndex) + _add_or_set(expr.linear_terms, arg, 0.0) + return push!(expr.nonlinear_terms, arg) +end + +# TODO(odow): these process_expr! functions use recursion. For large models, +# this may exceed the stack. At some point, we may have to rewrite this to not +# use recursion. +function _process_expr!(expr::_NLExpr, arg::Expr) + if arg.head == :call + f = get(_UNARY_SPECIAL_CASES, arg.args[1], nothing) + if f !== nothing + if length(arg.args) != 2 + error("Uncorrect number of arguments to $(arg.args[1]).") + end + # Some unary-functions are special cased. See the associated comment + # next to the definition of _UNARY_SPECIAL_CASES. + _process_expr!(expr, f(arg.args[2])) + else + _process_expr!(expr, arg.args) + end + elseif arg.head == :ref + _process_expr!(expr, arg.args[2]) + elseif arg == :() + return # Some evalators return a null objective of `:()`. + else + error("Unsupported expression: $(arg)") + end + return +end + +function _process_expr!(expr::_NLExpr, args::Vector{Any}) + op = first(args) + N = length(args) - 1 + # Before processing the arguments, do some re-writing. + if op == :+ + if N == 1 # +x, so we can just drop the op and process the args. + return _process_expr!(expr, args[2]) + elseif N > 2 # nary-addition! + op = :sum + end + elseif op == :- && N == 1 + op = :neg + elseif op == :* && N > 2 # nary-multiplication. + # NL doesn't define an nary multiplication operator, so we need to + # rewrite our expression as a sequence of chained binary operators. + while N > 2 + # Combine last term with previous to form a binary * expression + arg = pop!(args) + args[end] = Expr(:call, :*, args[end], arg) + N = length(args) - 1 + end + end + # Now convert the Julia expression into an _NLExpr. + opcode = get(_JULIA_TO_AMPL, op, nothing) + if opcode === nothing + error("Unsupported operation $(op)") + end + push!(expr.nonlinear_terms, opcode) + if opcode in _NARY_OPCODES + push!(expr.nonlinear_terms, N) + end + for i in 1:N + _process_expr!(expr, args[i+1]) + end + return +end + +### ============================================================================ +### Evaluate nonlinear expressions +### ============================================================================ + +function _evaluate(expr::_NLExpr, x::Dict{MOI.VariableIndex,Float64}) + y = expr.constant + for (v, c) in expr.linear_terms + y += c * x[v] + end + if length(expr.nonlinear_terms) > 0 + ret, n = _evaluate(expr.nonlinear_terms[1], expr.nonlinear_terms, x, 1) + @assert n == length(expr.nonlinear_terms) + 1 + y += ret + end + return y +end + +function _evaluate( + head::MOI.VariableIndex, + ::Vector{_NLTerm}, + x::Dict{MOI.VariableIndex,Float64}, + head_i::Int, +)::Tuple{Float64,Int} + return x[head], head_i + 1 +end + +function _evaluate( + head::Float64, + ::Vector{_NLTerm}, + ::Dict{MOI.VariableIndex,Float64}, + head_i::Int, +)::Tuple{Float64,Int} + return head, head_i + 1 +end + +function _evaluate( + head::Int, + terms::Vector{_NLTerm}, + x::Dict{MOI.VariableIndex,Float64}, + head_i::Int, +)::Tuple{Float64,Int} + N, f = _AMPL_TO_JULIA[head] + is_nary = (N == -1) + head_i += 1 + if is_nary + N = terms[head_i]::Int + head_i += 1 + end + args = Vector{Float64}(undef, N) + for n in 1:N + args[n], head_i = _evaluate(terms[head_i], terms, x, head_i) + end + return is_nary ? f(args) : f(args...), head_i +end diff --git a/src/nl_convert.jl b/src/nl_convert.jl deleted file mode 100644 index ed8e0d7..0000000 --- a/src/nl_convert.jl +++ /dev/null @@ -1,111 +0,0 @@ -# Converts terms in the expression tree to .nl supported expressions -convert_formula(c) = c -convert_formula(c::LinearityExpr) = convert_formula(c.c) - -function convert_formula(c::Expr) - map!(convert_formula, c.args, c.args) - - if c.head == :comparison - n = length(c.args) - @assert isodd(n) - - @assert(n > 3) - - # If more than binary comparison, we need to chain them together - # Get binary comparisons in sequence and chain with nested && - # It looks like we might be able to use an n-ary &&, but that's not how - # Julia parses chained && expressions. - if n > 3 - new_expr = extract_binary_comparison(c, 1) - for i in 3:2:(n-2) - new_expr = Expr(:&&, new_expr, extract_binary_comparison(c, i)) - end - c = new_expr - end - - elseif c.head == :call - op = c.args[1] - # Distinguish unary, binary and n-ary plus - if op == :+ - if length(c.args) > 3 - # n-ary plus to sum - c.args[1] = :sum - elseif length(c.args) == 2 - # Remove unary plus - c = c.args[2] - end - # Distinguish unary and binary minus - elseif op == :- - n = length(c.args) - if n == 2 - if c.args[2] == 0 - c = 0 - else - c.args[1] = :neg - end - end - # .nl has no n-ary multiplication so we need to convert to binary - elseif op == :* - # Loop for each arg after the 3rd (more than binary) - for _ in 1:(length(c.args)-3) - # Combine last term with previous to form a binary * expression - arg = pop!(c.args) - c.args[end] = :($(c.args[end]) * $arg) - end - # Handle normal conversion cases - else - if length(c.args) == 2 - c = get( - unary_special_cases, - c.args[1], - (x) -> Expr(:call, c.args[1], x), - )( - c.args[2], - ) - end - end - end - return c -end - -# Extracts `expression relation expression` from a larger comparison expression -function extract_binary_comparison(c::Expr, start::Integer) - @assert c.head == :comparison - @assert start <= length(c.args) - 2 - - return Expr(:call, c.args[start+1], c.args[start], c.args[start+2]) -end - -const unary_special_cases = Dict( - :cbrt => (x) -> :($x^(1 / 3)), - :abs2 => (x) -> :($x^2), - :inv => (x) -> :(1 / $x), - :log2 => (x) -> :(log($x) / log(2)), - :log1p => (x) -> :(log(1 + $x)), - :exp2 => (x) -> :(2^$x), - :expm1 => (x) -> :(exp($x) - 1), - :sec => (x) -> :(1 / cos($x)), - :csc => (x) -> :(1 / sin($x)), - :cot => (x) -> :(1 / tan($x)), - :asec => (x) -> :(acos(1 / $x)), - :acsc => (x) -> :(asin(1 / $x)), - :acot => (x) -> :(pi / 2 - atan($x)), - :sind => (x) -> :(sin(pi / 180 * $x)), - :cosd => (x) -> :(cos(pi / 180 * $x)), - :tand => (x) -> :(tan(pi / 180 * $x)), - :secd => (x) -> :(1 / cos(pi / 180 * $x)), - :cscd => (x) -> :(1 / sin(pi / 180 * $x)), - :cotd => (x) -> :(1 / tan(pi / 180 * $x)), - :asind => (x) -> :(asin($x) * 180 / pi), - :acosd => (x) -> :(acos($x) * 180 / pi), - :atand => (x) -> :(atan($x) * 180 / pi), - :asecd => (x) -> :(acos(1 / $x) * 180 / pi), - :acscd => (x) -> :(asin(1 / $x) * 180 / pi), - :acotd => (x) -> :((pi / 2 - atan($x)) * 180 / pi), - :sech => (x) -> :(1 / cosh($x)), - :csch => (x) -> :(1 / sinh($x)), - :coth => (x) -> :(1 / tanh($x)), - :asech => (x) -> :(acosh(1 / $x)), - :acsch => (x) -> :(asinh(1 / $x)), - :acoth => (x) -> :(atanh(1 / $x)), -) diff --git a/src/nl_linearity.jl b/src/nl_linearity.jl deleted file mode 100644 index 5e261bb..0000000 --- a/src/nl_linearity.jl +++ /dev/null @@ -1,263 +0,0 @@ -mutable struct LinearityExpr - c::Any - linearity::Symbol - coeff::Float64 -end - -LinearityExpr(c, linearity) = LinearityExpr(c, linearity, 1) - -eval(c::LinearityExpr) = eval(get_expr(c)) -Base.print(io::IO, c::LinearityExpr) = print(io::IO, "($(c.c),$(c.linearity))") -Base.show(io::IO, c::LinearityExpr) = print(io::IO, "($(c.c),$(c.linearity))") - -LinearityExpr(c) = throw("Couldn't determine linearity of expression:\n$c") -LinearityExpr(c::Real) = LinearityExpr(c, :const) -LinearityExpr(c::LinearityExpr) = c # Already processed; ignore! -function LinearityExpr(c::Expr) - if c.head == :call - for i in 2:length(c.args) - c.args[i] = LinearityExpr(c.args[i]) - end - - args = c.args[2:end] - if c.args[1] in [:+, :-] - if check_for_linearity(:nonlinear, args) - linearity = :nonlinear - elseif check_for_linearity(:linear, args) - linearity = :linear - else - linearity = :const - end - elseif c.args[1] == :* - if check_for_linearity(:nonlinear, args) || - length(filter((c) -> check_linearity(:linear, c), args)) > 1 - linearity = :nonlinear - elseif check_for_linearity(:linear, args) - linearity = :linear - else - linearity = :const - end - elseif c.args[1] == :/ - if c.args[3].linearity != :const - linearity = :nonlinear - else - linearity = c.args[2].linearity - end - elseif c.args[1] == :^ - if c.args[3].linearity != :const - linearity = :nonlinear - elseif c.args[2].linearity == :linear - if c.args[3].c == 1 - linearity = :linear - else - linearity = :nonlinear - end - else - linearity = c.args[2].linearity - end - elseif c.args[1] == :ifelse - if c.args[2].linearity == :const - # We know which branch to take already - simplify - c.args[2] = pull_up_constants(c.args[2]) - if c.args[2].c != 0 - linearity = c.args[3].linearity - c = c.args[3].c - else - linearity = c.args[4].linearity - c = c.args[4].c - end - else - linearity = :nonlinear - end - else - if check_for_linearity(:linear, args) || - check_for_linearity(:nonlinear, args) - linearity = :nonlinear - else - linearity = :const - end - end - - elseif c.head == :ref - linearity = :linear - elseif c.head == :comparison - indices = 1:2:length(c.args) - for i in indices - c.args[i] = LinearityExpr(c.args[i]) - end - if check_for_linearity(:linear, c.args[indices]) || - check_for_linearity(:nonlinear, c.args[indices]) - linearity = :nonlinear - else - linearity = :const - end - elseif c.head in [:&&, :||] - for i in 1:length(c.args) - c.args[i] = LinearityExpr(c.args[i]) - end - if check_for_linearity(:linear, c.args) || - check_for_linearity(:nonlinear, c.args) - linearity = :nonlinear - else - linearity = :const - end - end - return LinearityExpr(c, linearity) -end - -check_linearity(linearity::Symbol, c::LinearityExpr) = c.linearity == linearity - -function check_for_linearity(linearity::Symbol, args::Array) - return !isempty(filter((c) -> check_linearity(linearity, c), args)) -end - -get_expr(c) = c -get_expr(c::LinearityExpr) = get_expr(c.c) -function get_expr(c::Expr) - map!(get_expr, c.args, c.args) - return c -end - -pull_up_constants(c) = c -function pull_up_constants(c::LinearityExpr) - if c.linearity == :const - c.c = eval(c) - elseif isa(c.c, Expr) - map!(pull_up_constants, c.c.args, c.c.args) - end - return c -end - -function prune_linear_terms!( - c::LinearityExpr, - lin_constr::Dict{Int,Float64}, - constant::Float64 = 0.0, - negative_tree::Bool = false, -) - if c.linearity != :nonlinear - constant = add_linear_tree!(c, lin_constr, constant, negative_tree) - c = LinearityExpr(:(0), :const) - return true, c, constant - else - expr = c.c - if expr.head == :call - if expr.args[1] == :+ - n = length(expr.args) - pruned = falses(n - 1) - for i in 2:n - pruned[i-1], expr.args[i], constant = prune_linear_terms!( - expr.args[i], - lin_constr, - constant, - negative_tree, - ) - end - if sum(.!pruned) > 1 - inds = vcat([1], collect(2:n)[.!pruned]) - c.c.args = expr.args[inds] - else - c = expr.args[findfirst(.!pruned)+1] - end - elseif expr.args[1] == :- - if length(expr.args) == 3 - pruned_first, expr.args[2], constant = prune_linear_terms!( - expr.args[2], - lin_constr, - constant, - negative_tree, - ) - pruned_second, expr.args[3], constant = prune_linear_terms!( - expr.args[3], - lin_constr, - constant, - !negative_tree, - ) - if pruned_first - new_expr = Expr(:call, :-, expr.args[3]) - c = LinearityExpr(new_expr, :nonlinear) - elseif pruned_second - c = expr.args[2] - end - end - end - end - return false, c, constant - end -end - -function add_linear_tree!( - c::LinearityExpr, - lin_constr::Dict{Int,Float64}, - constant::Float64 = 0.0, - negative_tree::Bool = false, -) - c = collate_linear_terms(c) - negative_tree && negate(c) - constant = add_tree_to_constr!(c, lin_constr, constant) - return constant -end - -function collate_linear_terms(c::LinearityExpr) - if isa(c.c, Expr) && c.c.head == :call - for i in 2:length(c.c.args) - c.c.args[i] = collate_linear_terms(c.c.args[i]) - end - end - if c.linearity == :linear && isa(c.c, Expr) && c.c.head == :call - func = c.c.args[1] - if func == :- - if length(c.c.args) == 2 - c = negate(c.c.args[2]) - else - c.c.args[3] = negate(c.c.args[3]) - c.c.args[1] = :+ - end - elseif func == :* - # This can be n-ary multiplication, but there is only one variable - # First, find index of variable in the args - x = findfirst(ex -> ex.linearity == :linear, c.c.args[2:end]) + 1 - # Multiply variable term by each of the constants in the args - x_expr = c.c.args[x] - for i in setdiff(2:length(c.c.args), x) - x_expr = multiply(x_expr, c.c.args[i].c) - end - c = x_expr - elseif func == :/ - c = multiply(c.c.args[2], 1 / c.c.args[3].c) - elseif func == :^ - c = c.c.args[2] - end - end - return c -end - -negate(c::LinearityExpr) = multiply(c::LinearityExpr, -1) -function multiply(c::LinearityExpr, a::Real) - if isa(c.c, Expr) && c.c.head == :call - @assert c.c.args[1] == :+ - map!(arg -> multiply(arg, a), c.c.args[2:end], c.c.args[2:end]) - elseif c.linearity == :const - c.c *= a - else - c.coeff *= a - end - return c -end - -function add_tree_to_constr!( - c::LinearityExpr, - lin_constr::Dict{Int,Float64}, - constant::Float64, -) - if isa(c.c, Expr) && c.c.head == :call - @assert c.c.args[1] == :+ - for i in 2:length(c.c.args) - constant = add_tree_to_constr!(c.c.args[i], lin_constr, constant) - end - elseif c.linearity == :const - constant += c.c - else - lin_constr[c.c.args[2]] += c.coeff - end - return constant -end diff --git a/src/nl_params.jl b/src/nl_params.jl deleted file mode 100644 index 651b73f..0000000 --- a/src/nl_params.jl +++ /dev/null @@ -1,74 +0,0 @@ -const func_to_nl = Dict( - :+ => 0, - :- => 1, - :* => 2, - :/ => 3, - # :rem => 4, - :^ => 5, - # :less => 6, - :min => 11, - :max => 12, - # :floor => 13, - # :ceil => 14, - :abs => 15, - :neg => 16, - :|| => 20, - :&& => 21, - :< => 22, - :<= => 23, - :(==) => 24, - :>= => 28, - :> => 29, - :!= => 30, - :ifelse => 35, - # :not => 34, - :tanh => 37, - :tan => 38, - :sqrt => 39, - :sinh => 40, - :sin => 41, - :log10 => 42, - :log => 43, - :exp => 44, - :cosh => 45, - :cos => 46, - :atanh => 47, - # :atan2 => 48, - :atan => 49, - :asinh => 50, - :asin => 51, - :acosh => 52, - :acos => 53, - :sum => 54, - # :intdiv => 55, - # :precision => 56, - # :round => 57, - # :trunc => 58, - # :count => 59, - # :numberof => 60, - # :numberofs => 61, - # :ifs => 65, - # :and_n => 70, - # :or_n => 71, - # :implies => 72, - # :iff => 73, - # :alldiff => 74, -) - -const sense_to_nl = Dict(:Min => 0, :Max => 1) - -const relation_to_nl = Dict(:multiple => 0, :<= => 1, :>= => 2, :(==) => 4) - -const reverse_relation = Dict(:0 => 0, :1 => 2, :2 => 1, :4 => 4) - -const nary_functions = Set{Symbol}([ - :min, - :max, - :sum, - # :count, - # :numberof, - # :numberofs, - # :and_n, - # :or_n, - # :alldiff, -]) diff --git a/src/nl_write.jl b/src/nl_write.jl deleted file mode 100644 index 3ca90be..0000000 --- a/src/nl_write.jl +++ /dev/null @@ -1,266 +0,0 @@ -function write_nl_file(f::IO, m::AmplNLMathProgModel) - write_nl_header(f, m) - - if m.ncon > 0 - write_nl_c_blocks(f, m) - end - - if has_objective(m) - write_nl_o_block(f, m) - end - - if m.ncon > 0 - write_nl_d_block(f, m) - end - - write_nl_x_block(f, m) - - if m.ncon > 0 - write_nl_r_block(f, m) - end - - write_nl_b_block(f, m) - - if m.ncon > 0 - write_nl_k_block(f, m) - write_nl_j_blocks(f, m) - end - - if has_objective(m) - write_nl_g_block(f, m) - end -end - -function has_objective(m::AmplNLMathProgModel) - return !isempty(m.lin_obj) || (isa(m.obj, Expr)) -end - -function write_nl_header(f, m::AmplNLMathProgModel) - # Line 1: Always the same - println(f, "g3 1 1 0") - # Line 2: vars, constraints, objectives, ranges, eqns, logical constraints - n_ranges = sum(m.r_codes .== 0) - n_eqns = sum(m.r_codes .== 4) - nobj = has_objective(m) ? 1 : 0 - println(f, " $(m.nvar) $(m.ncon) $nobj $n_ranges $n_eqns 0") - # Line 3: nonlinear constraints, objectives - nlc = sum(m.conlinearities .== :Nonlin) - nlo = m.objlinearity == :Nonlin ? 1 : 0 - println(f, " $nlc $nlo") - # Line 4: network constraints: nonlinear, linear - println(f, " 0 0") - # Line 5: nonlinear vars in constraints, objectives, both - nonlinear_obj = m.varlinearities_obj .== :Nonlin - nonlinear_con = m.varlinearities_con .== :Nonlin - nonlinear = (nonlinear_con + nonlinear_obj) .> 0 - nonlinear_both = (nonlinear_con + nonlinear_obj) .> 1 - nlvc = sum(nonlinear_con .> 0) - nlvo = sum(nonlinear_obj .> 0) - nlvb = sum(nonlinear_both) - println(f, " $nlvc $nlvo $nlvb") - # Line 6: linear network variables; functions; arith, flags - println(f, " 0 0 0 1") # flags set to 1 to get suffixes in .sol file - # Line 7: discrete variables: binary, integer, nonlinear (b,c,o) - binary = m.vartypes .== :Bin - integer = m.vartypes .== :Int - discrete = binary + integer .> 0 - # Julia 0.6 syntax - nbv = sum(binary + .!nonlinear .> 1) - niv = sum(integer + .!nonlinear .> 1) - nlvbi = sum(nonlinear_both + discrete .> 1) - nlvci = sum(nonlinear_con - nonlinear_obj + discrete .> 1) - nlvoi = sum(nonlinear_obj - nonlinear_con + discrete .> 1) - println(f, " $nbv $niv $nlvbi $nlvci $nlvoi") - # Line 8: nonzeros in Jacobian, gradients - nzc = sum(m.j_counts) - nzo = length(m.lin_obj) - println(f, " $nzc $nzo") - # Line 9: max name lengths: constraints, variables - println(f, " 0 0") - # Line 10: common exprs: b,c,o,c1,o1 - return println(f, " 0 0 0 0 0") -end - -# Nonlinear constraint trees -function write_nl_c_blocks(f, m::AmplNLMathProgModel) - for index in 0:(m.ncon-1) - i = m.c_index_map_rev[index] - println(f, "C$index") - write_nl_expr(f, m, m.constrs[i]) - end -end - -# Nonlinear objective tree -function write_nl_o_block(f, m::AmplNLMathProgModel) - println(f, string("O0 ", sense_to_nl[m.sense])) - return write_nl_expr(f, m, m.obj) -end - -# Initial dual guesses - unused -function write_nl_d_block(f, m::AmplNLMathProgModel) - println(f, "d$(m.ncon)") - for index in 0:(m.ncon-1) - i = m.c_index_map_rev[index] - println(f, "$index 0") - end -end - -# Initial primal guesses -function write_nl_x_block(f, m::AmplNLMathProgModel) - println(f, "x$(m.nvar)") - for index in 0:(m.nvar-1) - i = m.v_index_map_rev[index] - println(f, "$index $(m.x_0[i])") - end -end - -# Constraint bounds -function write_nl_r_block(f, m::AmplNLMathProgModel) - println(f, "r") - for index in 0:(m.ncon-1) - i = m.c_index_map_rev[index] - lower = m.g_l[i] - upper = m.g_u[i] - rel = m.r_codes[i] - if rel == 0 - println(f, "$rel $lower $upper") - elseif rel == 1 - println(f, "$rel $upper") - elseif rel == 2 - println(f, "$rel $lower") - elseif rel == 3 - println(f, "$rel") - elseif rel == 4 - println(f, "$rel $lower") - end - end -end - -# Variable bounds -function write_nl_b_block(f, m::AmplNLMathProgModel) - println(f, "b") - for index in 0:(m.nvar-1) - i = m.v_index_map_rev[index] - lower = m.x_l[i] - upper = m.x_u[i] - if lower == -Inf - if upper == Inf - println(f, "3") - else - println(f, "1 $upper") - end - else - if lower == upper - println(f, "4 $lower") - elseif upper == Inf - println(f, "2 $lower") - else - println(f, "0 $lower $upper") - end - end - end -end - -# Jacobian counts -function write_nl_k_block(f, m::AmplNLMathProgModel) - println(f, "k$(m.nvar - 1)") - total = 0 - for index in 0:(m.nvar-2) - i = m.v_index_map_rev[index] - total += m.j_counts[i] - println(f, total) - end -end - -# Linear constraint expressions -function write_nl_j_blocks(f, m::AmplNLMathProgModel) - for index in 0:(m.ncon-1) - i = m.c_index_map_rev[index] - num_vars = length(m.lin_constrs[i]) - # Only print linear blocks with vars, .nl file is malformed otherwise - if num_vars > 0 - println(f, "J$index $num_vars") - - # We need to output .nl index and constraint coeff, ordered by .nl - # index. `lin_constrs` contains our variable index as the key and - # constraint coeff as the value - - # Assemble tuples of (.nl index, constraint value) - output = collect( - zip( - (m.v_index_map[j] for j in keys(m.lin_constrs[i])), - values(m.lin_constrs[i]), - ), - ) - - # Loop through output in .nl index order - for (index2, value) in sort(output) - println(f, "$index2 $value") - end - end - end -end - -# Linear objective expression -function write_nl_g_block(f, m::AmplNLMathProgModel) - println(f, string("G0 ", length(m.lin_obj))) - for index in 0:(m.nvar-1) - i = m.v_index_map_rev[index] - if i in keys(m.lin_obj) - println(f, "$index $(m.lin_obj[i])") - end - end -end - -# Convert an expression tree (with .nl formulae only) to .nl format -write_nl_expr(f, m, c) = println(f, string(c)) -# Handle numerical constants e.g. pi -write_nl_expr(f, m, c::Symbol) = write_nl_expr(f, m, float(eval(c))) -function write_nl_expr(f, m, c::Real) - return println(f, nl_number(c == round(Integer, c) ? round(Integer, c) : c)) -end -write_nl_expr(f, m, c::LinearityExpr) = write_nl_expr(f, m, c.c) -function write_nl_expr(f, m, c::Expr) - if c.head == :ref - # Output variable as `v$index` - if c.args[1] == :x - @assert isa(c.args[2], Int) - println(f, nl_variable(m.v_index_map[c.args[2]])) - else - error("Unrecognized reference expression $c") - end - elseif c.head == :call - # Output function as `o$opcode` - println(f, nl_operator(c.args[1])) - if c.args[1] in nary_functions - # Output nargs on subsequent line if n-ary function - println(f, (string(length(c.args) - 1))) - end - map(arg -> write_nl_expr(f, m, arg), c.args[2:end]) - - elseif c.head == :comparison - # .nl only handles binary comparison - @assert length(c.args) == 3 - # Output comparison type first, followed by args - println(f, nl_operator(c.args[2])) - map(arg -> write_nl_expr(f, m, arg), c.args[1:2:end]) - - elseif c.head in [:&&, :||] - # Only support binary and/or for now - @assert length(c.args) == 2 - println(f, nl_operator(c.head)) - map(arg -> write_nl_expr(f, m, arg), c.args) - else - error("Unrecognized expression $c") - end -end - -nl_variable(index::Integer) = "v$index" -nl_number(value::Real) = "n$value" - -function nl_operator(operator::Symbol) - if !haskey(func_to_nl, operator) - error("translation of the function \"$operator\" to NL is not defined") - end - return "o$(func_to_nl[operator])" -end diff --git a/src/opcode.jl b/src/opcode.jl new file mode 100644 index 0000000..2a066b1 --- /dev/null +++ b/src/opcode.jl @@ -0,0 +1,72 @@ +# Do not modify. This file is automatically created by the script in `gen.jl`. + +const OPPLUS = 0 +const OPMINUS = 1 +const OPMULT = 2 +const OPDIV = 3 +const OPREM = 4 +const OPPOW = 5 +const OPLESS = 6 +const MINLIST = 11 +const MAXLIST = 12 +const FLOOR = 13 +const CEIL = 14 +const ABS = 15 +const OPUMINUS = 16 +const OPOR = 20 +const OPAND = 21 +const LT = 22 +const LE = 23 +const EQ = 24 +const GE = 28 +const GT = 29 +const NE = 30 +const OPNOT = 34 +const OPIFnl = 35 +const OP_tanh = 37 +const OP_tan = 38 +const OP_sqrt = 39 +const OP_sinh = 40 +const OP_sin = 41 +const OP_log10 = 42 +const OP_log = 43 +const OP_exp = 44 +const OP_cosh = 45 +const OP_cos = 46 +const OP_atanh = 47 +const OP_atan2 = 48 +const OP_atan = 49 +const OP_asinh = 50 +const OP_asin = 51 +const OP_acosh = 52 +const OP_acos = 53 +const OPSUMLIST = 54 +const OPintDIV = 55 +const OPprecision = 56 +const OPround = 57 +const OPtrunc = 58 +const OPCOUNT = 59 +const OPNUMBEROF = 60 +const OPNUMBEROFs = 61 +const OPATLEAST = 62 +const OPATMOST = 63 +const OPPLTERM = 64 +const OPIFSYM = 65 +const OPEXACTLY = 66 +const OPNOTATLEAST = 67 +const OPNOTATMOST = 68 +const OPNOTEXACTLY = 69 +const ANDLIST = 70 +const ORLIST = 71 +const OPIMPELSE = 72 +const OP_IFF = 73 +const OPALLDIFF = 74 +const OPSOMESAME = 75 +const OP1POW = 76 +const OP2POW = 77 +const OPCPOW = 78 +const OPFUNCALL = 79 +const OPNUM = 80 +const OPHOL = 81 +const OPVARVAL = 82 +const N_OPS = 83 diff --git a/test/MINLPTests/run_minlptests.jl b/test/MINLPTests/run_minlptests.jl index 7ba92e0..7bbdb02 100644 --- a/test/MINLPTests/run_minlptests.jl +++ b/test/MINLPTests/run_minlptests.jl @@ -10,8 +10,14 @@ else run_with_ampl(f) = Ipopt_jll.amplexe(f) end +const MOI = AmplNLWriter.MOI + run_with_ampl() do path - OPTIMIZER = () -> AmplNLWriter.Optimizer(path, ["print_level=0"]) + OPTIMIZER = + () -> MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + AmplNLWriter.Optimizer(path, ["print_level=0"]), + ) ### ### src/nlp tests. ### @@ -35,7 +41,10 @@ run_with_ampl() do path NaN, NaN, Dict(MINLPTests.INFEASIBLE_PROBLEM => AmplNLWriter.MOI.INFEASIBLE), - Dict(MINLPTests.INFEASIBLE_PROBLEM => AmplNLWriter.MOI.NO_SOLUTION), + Dict( + MINLPTests.INFEASIBLE_PROBLEM => + AmplNLWriter.MOI.UNKNOWN_RESULT_STATUS, + ), ) end diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 633289e..4ab9911 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -15,15 +15,20 @@ const CONFIG = MOI.Test.TestConfig( ) function optimizer(path) - return MOI.Bridges.full_bridge_optimizer( - AmplNLWriter.Optimizer(path, ["print_level = 0"]), - Float64, + return MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + MOI.Bridges.full_bridge_optimizer( + MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + AmplNLWriter.Optimizer(path, ["print_level = 0"]), + ), + Float64, + ), ) end function test_name(path) - @test sprint(show, AmplNLWriter.Optimizer(path, ["print_level = 0"])) == - "An AmplNLWriter model" + @test sprint(show, AmplNLWriter.Optimizer(path)) == "An AMPL (.nl) model" end function test_unittest(path) @@ -33,7 +38,6 @@ function test_unittest(path) [ # Unsupported attributes: "number_threads", - "raw_status_string", "silent", "solve_objbound_edge_cases", "solve_time", @@ -44,21 +48,14 @@ function test_unittest(path) "solve_zero_one_with_bounds_2", "solve_zero_one_with_bounds_3", - # It seems that the AMPL NL reader declares NL files with no objective - # and no constraints as corrupt, even if they have variable bounds. Yuk. - "solve_blank_obj", - # No support for VectorOfVariables-in-SecondOrderCone "delete_soc_variables", - - # TODO(odow): fix handling of result indices. - "solve_result_index", ], ) end function test_contlinear(path) - return MOI.Test.contlineartest(optimizer(path), CONFIG, String["linear15",]) + return MOI.Test.contlineartest(optimizer(path), CONFIG) end function test_contlquadratic(path) @@ -98,13 +95,7 @@ function test_orderedindices(path) end function test_copytest(path) - return MOI.Test.copytest( - optimizer(path), - MOI.Bridges.full_bridge_optimizer( - AmplNLWriter.Optimizer(path), - Float64, - ), - ) + return MOI.Test.copytest(optimizer(path), optimizer(path)) end function test_nlptest(path) @@ -112,13 +103,25 @@ function test_nlptest(path) end function test_bad_string(::Any) - model = AmplNLWriter.Optimizer("bad_solver") + model = optimizer("bad_solver") x = MOI.add_variable(model) MOI.optimize!(model) @test MOI.get(model, MOI.TerminationStatus()) == MOI.OTHER_ERROR @test occursin("IOError", MOI.get(model, MOI.RawStatusString())) end +function test_function_constant_nonzero(path) + model = optimizer(path) + x = MOI.add_variable(model) + f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 1.0) + MOI.add_constraint(model, f, MOI.GreaterThan(3.0)) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.optimize!(model) + @test isapprox(MOI.get(model, MOI.VariablePrimal(), x), 2.0, atol = 1e-6) + @test isapprox(MOI.get(model, MOI.ObjectiveValue()), 3.0, atol = 1e-6) +end + function runtests(path) for name in names(@__MODULE__; all = true) if !startswith("$(name)", "test_") @@ -132,6 +135,10 @@ end end -run_with_ampl() do path - return TestMOIWrapper.runtests(path) +if VERSION < v"1.3" + import Ipopt + TestMOIWrapper.runtests(Ipopt.amplexe) +else + import Ipopt_jll + TestMOIWrapper.runtests(Ipopt_jll.amplexe) end diff --git a/test/NLModel.jl b/test/NLModel.jl new file mode 100644 index 0000000..ccdc8d9 --- /dev/null +++ b/test/NLModel.jl @@ -0,0 +1,692 @@ +module TestNLModel + +using AmplNLWriter +const NL = AmplNLWriter + +using MathOptInterface +const MOI = MathOptInterface + +using Test + +function _test_nlexpr(expr::NL._NLExpr, nonlinear_terms, linear_terms, constant) + @test expr.is_linear == (length(nonlinear_terms) == 0) + @test expr.nonlinear_terms == nonlinear_terms + @test expr.linear_terms == linear_terms + @test expr.constant == constant + return +end +_test_nlexpr(x, args...) = _test_nlexpr(NL._NLExpr(x), args...) + +function test_nlexpr_singlevariable() + x = MOI.VariableIndex(1) + _test_nlexpr(MOI.SingleVariable(x), NL._NLTerm[], Dict(x => 1.0), 0.0) + return +end + +function test_nlexpr_scalaraffine() + x = MOI.VariableIndex.(1:3) + f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 4.0) + return _test_nlexpr(f, NL._NLTerm[], Dict(x .=> 1), 4.0) +end + +function test_nlexpr_scalarquadratic() + x = MOI.VariableIndex(1) + f = MOI.ScalarQuadraticFunction( + [MOI.ScalarAffineTerm(1.1, x)], + [MOI.ScalarQuadraticTerm(2.0, x, x)], + 3.0, + ) + terms = [NL.OPMULT, x, x] + return _test_nlexpr(f, terms, Dict(x => 1.1), 3.0) +end + +function test_nlexpr_unary_addition() + x = MOI.VariableIndex(1) + return _test_nlexpr(:(+$x), [x], Dict(x => 0), 0.0) +end + +function test_nlexpr_binary_addition() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + return _test_nlexpr( + :($x + $y), + [NL.OPPLUS, x, y], + Dict(x => 0.0, y => 0.0), + 0.0, + ) +end + +function test_nlexpr_nary_addition() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + return _test_nlexpr( + :($x + $y + 1.0), + [NL.OPSUMLIST, 3, x, y, 1.0], + Dict(x => 0.0, y => 0.0), + 0.0, + ) +end + +function test_nlexpr_unary_subtraction() + x = MOI.VariableIndex(1) + return _test_nlexpr(:(-$x), [NL.OPUMINUS, x], Dict(x => 0.0), 0.0) +end + +function test_nlexpr_nary_multiplication() + x = MOI.VariableIndex(1) + return _test_nlexpr( + :($x * $x * 2.0), + [NL.OPMULT, x, NL.OPMULT, x, 2.0], + Dict(x => 0.0), + 0.0, + ) +end + +function test_nlexpr_unary_specialcase() + x = MOI.VariableIndex(1) + return _test_nlexpr( + :(cbrt($x)), + [NL.OPPOW, x, NL.OPDIV, 1, 3], + Dict(x => 0.0), + 0.0, + ) +end + +function test_nlexpr_unsupportedoperation() + x = MOI.VariableIndex(1) + err = ErrorException("Unsupported operation foo") + @test_throws err NL._NLExpr(:(foo($x))) + return +end + +function test_nlexpr_unsupportedexpression() + x = MOI.VariableIndex(1) + expr = :(1 <= $x <= 2) + err = ErrorException("Unsupported expression: $(expr)") + @test_throws err NL._NLExpr(expr) + return +end + +function test_nlexpr_ref() + x = MOI.VariableIndex(1) + return _test_nlexpr(:(x[$x]), [x], Dict(x => 0.0), 0.0) +end + +function test_nlconstraint_interval() + x = MOI.VariableIndex(1) + expr = :(1 <= $x <= 2) + con = NL._NLConstraint(expr, MOI.NLPBoundsPair(1.0, 2.0)) + @test con.lower == 1 + @test con.upper == 2 + @test con.opcode == 0 + @test con.expr == NL._NLExpr(expr.args[3]) +end + +function test_nlconstraint_lessthan() + x = MOI.VariableIndex(1) + expr = :($x <= 2) + con = NL._NLConstraint(expr, MOI.NLPBoundsPair(-Inf, 2.0)) + @test con.lower == -Inf + @test con.upper == 2 + @test con.opcode == 1 + @test con.expr == NL._NLExpr(expr.args[2]) +end + +function test_nlconstraint_greaterthan() + x = MOI.VariableIndex(1) + expr = :($x >= 2) + con = NL._NLConstraint(expr, MOI.NLPBoundsPair(2.0, Inf)) + @test con.lower == 2 + @test con.upper == Inf + @test con.opcode == 2 + @test con.expr == NL._NLExpr(expr.args[2]) +end + +function test_nlconstraint_equalto() + x = MOI.VariableIndex(1) + expr = :($x == 2) + con = NL._NLConstraint(expr, MOI.NLPBoundsPair(2.0, 2.0)) + @test con.lower == 2 + @test con.upper == 2 + @test con.opcode == 4 + @test con.expr == NL._NLExpr(expr.args[2]) +end + +function test_nlconstraint_interval_warn() + x = MOI.VariableIndex(1) + expr = :(2 <= $x <= 1) + @test_logs (:warn,) NL._NLConstraint(expr, MOI.NLPBoundsPair(1.0, 2.0)) +end + +function test_nlconstraint_lessthan_warn() + x = MOI.VariableIndex(1) + expr = :(1 <= $x <= 2) + @test_logs (:warn,) NL._NLConstraint(expr, MOI.NLPBoundsPair(-Inf, 2.0)) +end + +function test_nlconstraint_greaterthan_warn() + x = MOI.VariableIndex(1) + expr = :(1 <= $x <= 2) + @test_logs (:warn,) NL._NLConstraint(expr, MOI.NLPBoundsPair(1.0, Inf)) +end + +function test_nlconstraint_equalto_warn() + x = MOI.VariableIndex(1) + expr = :(1 <= $x <= 2) + @test_logs (:warn,) NL._NLConstraint(expr, MOI.NLPBoundsPair(1.0, 1.0)) +end + +function test_nlmodel_hs071() + model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + v = MOI.add_variables(model, 4) + l = [1.1, 1.2, 1.3, 1.4] + u = [5.1, 5.2, 5.3, 5.4] + start = [2.1, 2.2, 2.3, 2.4] + MOI.add_constraint.(model, MOI.SingleVariable.(v), MOI.GreaterThan.(l)) + MOI.add_constraint.(model, MOI.SingleVariable.(v), MOI.LessThan.(u)) + MOI.set.(model, MOI.VariablePrimalStart(), v, start) + lb, ub = [25.0, 40.0], [Inf, 40.0] + evaluator = MOI.Test.HS071(true) + block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), evaluator, true) + MOI.set(model, MOI.NLPBlock(), block_data) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + n = NL.Optimizer() + MOI.copy_to(n, model) + @test n.sense == MOI.MIN_SENSE + @test n.f == NL._NLExpr(MOI.objective_expr(evaluator)) + _test_nlexpr( + n.g[1].expr, + [NL.OPMULT, v[1], NL.OPMULT, v[2], NL.OPMULT, v[3], v[4]], + Dict(v .=> 0.0), + 0.0, + ) + @test n.g[1].lower == 25.0 + @test n.g[1].upper == Inf + @test n.g[1].opcode == 2 + _test_nlexpr( + n.g[2].expr, + [ + NL.OPSUMLIST, + 4, + NL.OPPOW, + v[1], + 2, + NL.OPPOW, + v[2], + 2, + NL.OPPOW, + v[3], + 2, + NL.OPPOW, + v[4], + 2, + ], + Dict(v .=> 0.0), + 0.0, + ) + @test n.g[2].lower == 40.0 + @test n.g[2].upper == 40.0 + @test n.g[2].opcode == 4 + @test length(n.h) == 0 + for i in 1:4 + @test n.x[v[i]].lower == l[i] + @test n.x[v[i]].upper == u[i] + @test n.x[v[i]].type == NL._CONTINUOUS + @test n.x[v[i]].jacobian_count == 2 + @test n.x[v[i]].in_nonlinear_constraint + @test n.x[v[i]].in_nonlinear_objective + @test 0 <= n.x[v[i]].order <= 3 + end + @test length(n.types[1]) == 4 + @test sprint(write, n) == """ + g3 1 1 0 + 4 2 1 0 1 0 + 2 1 + 0 0 + 4 4 4 + 0 0 0 1 + 0 0 0 0 0 + 8 4 + 0 0 + 0 0 0 0 0 + C0 + o2 + v0 + o2 + v1 + o2 + v2 + v3 + C1 + o54 + 4 + o5 + v0 + n2 + o5 + v1 + n2 + o5 + v2 + n2 + o5 + v3 + n2 + O0 0 + o0 + o2 + v0 + o2 + v3 + o54 + 3 + v0 + v1 + v2 + v2 + x4 + 0 2.1 + 1 2.2 + 2 2.3 + 3 2.4 + r + 2 25 + 4 40 + b + 0 1.1 5.1 + 0 1.2 5.2 + 0 1.3 5.3 + 0 1.4 5.4 + k3 + 2 + 4 + 6 + J0 4 + 0 0 + 1 0 + 2 0 + 3 0 + J1 4 + 0 0 + 1 0 + 2 0 + 3 0 + G0 4 + 0 0 + 1 0 + 2 0 + 3 0 + """ + return +end + +function test_nlmodel_hs071_linear_obj() + model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + v = MOI.add_variables(model, 4) + l = [1.1, 1.2, 1.3, 1.4] + u = [5.1, 5.2, 5.3, 5.4] + start = [2.1, 2.2, 2.3, 2.4] + MOI.add_constraint.(model, MOI.SingleVariable.(v), MOI.GreaterThan.(l)) + MOI.add_constraint.(model, MOI.SingleVariable.(v), MOI.LessThan.(u)) + MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.ZeroOne()) + MOI.add_constraint(model, MOI.SingleVariable(v[3]), MOI.Integer()) + MOI.set.(model, MOI.VariablePrimalStart(), v, start) + lb, ub = [25.0, 40.0], [Inf, 40.0] + evaluator = MOI.Test.HS071(true) + block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(lb, ub), evaluator, false) + MOI.set(model, MOI.NLPBlock(), block_data) + f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(l, v), 2.0) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + n = NL.Optimizer() + MOI.copy_to(n, model) + @test n.sense == MOI.MAX_SENSE + @test n.f == NL._NLExpr(f) + _test_nlexpr( + n.g[1].expr, + [NL.OPMULT, v[1], NL.OPMULT, v[2], NL.OPMULT, v[3], v[4]], + Dict(v .=> 0.0), + 0.0, + ) + @test n.g[1].lower == 25.0 + @test n.g[1].upper == Inf + @test n.g[1].opcode == 2 + _test_nlexpr( + n.g[2].expr, + [ + NL.OPSUMLIST, + 4, + NL.OPPOW, + v[1], + 2, + NL.OPPOW, + v[2], + 2, + NL.OPPOW, + v[3], + 2, + NL.OPPOW, + v[4], + 2, + ], + Dict(v .=> 0.0), + 0.0, + ) + @test n.g[2].lower == 40.0 + @test n.g[2].upper == 40.0 + @test n.g[2].opcode == 4 + @test length(n.h) == 0 + types = [NL._CONTINUOUS, NL._BINARY, NL._INTEGER, NL._CONTINUOUS] + u[2] = 1.0 + for i in 1:4 + @test n.x[v[i]].lower == l[i] + @test n.x[v[i]].upper == u[i] + @test n.x[v[i]].type == types[i] + @test n.x[v[i]].jacobian_count == 2 + @test n.x[v[i]].in_nonlinear_constraint + @test !n.x[v[i]].in_nonlinear_objective + @test 0 <= n.x[v[i]].order <= 3 + end + @test length(n.types[3]) == 2 + @test length(n.types[4]) == 2 + @test v[1] in n.types[3] + @test v[2] in n.types[4] + @test v[3] in n.types[4] + @test v[4] in n.types[3] + @test sprint(write, n) == """ + g3 1 1 0 + 4 2 1 0 1 0 + 2 1 + 0 0 + 4 0 0 + 0 0 0 1 + 0 0 0 2 0 + 8 4 + 0 0 + 0 0 0 0 0 + C0 + o2 + v0 + o2 + v2 + o2 + v3 + v1 + C1 + o54 + 4 + o5 + v0 + n2 + o5 + v2 + n2 + o5 + v3 + n2 + o5 + v1 + n2 + O0 1 + n2 + x4 + 0 2.1 + 1 2.4 + 2 2.2 + 3 2.3 + r + 2 25 + 4 40 + b + 0 1.1 5.1 + 0 1.4 5.4 + 0 1.2 1 + 0 1.3 5.3 + k3 + 2 + 4 + 6 + J0 4 + 0 0 + 1 0 + 2 0 + 3 0 + J1 4 + 0 0 + 1 0 + 2 0 + 3 0 + G0 4 + 0 1.1 + 1 1.4 + 2 1.2 + 3 1.3 + """ + return +end + +function test_nlmodel_linear_quadratic() + model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + x = MOI.add_variables(model, 4) + MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.GreaterThan(0.0)) + MOI.add_constraint.(model, MOI.SingleVariable.(x), MOI.LessThan(2.0)) + MOI.add_constraint(model, MOI.SingleVariable(x[2]), MOI.ZeroOne()) + MOI.add_constraint(model, MOI.SingleVariable(x[3]), MOI.Integer()) + f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x[2:4]), 2.0) + g = MOI.ScalarQuadraticFunction( + [MOI.ScalarAffineTerm(1.0, x[1])], + [MOI.ScalarQuadraticTerm(2.0, x[1], x[2])], + 3.0, + ) + h = MOI.ScalarQuadraticFunction( + [MOI.ScalarAffineTerm(1.0, x[3])], + [MOI.ScalarQuadraticTerm(1.0, x[1], x[2])], + 0.0, + ) + MOI.add_constraint(model, f, MOI.Interval(1.0, 10.0)) + MOI.add_constraint(model, g, MOI.LessThan(5.0)) + MOI.set(model, MOI.ObjectiveFunction{typeof(h)}(), h) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + n = NL.Optimizer() + MOI.copy_to(n, model) + @test n.sense == MOI.MAX_SENSE + @test n.f == NL._NLExpr(h) + terms = [NL.OPMULT, 2.0, NL.OPMULT, x[1], x[2]] + _test_nlexpr(n.g[1].expr, terms, Dict(x[1] => 1.0, x[2] => 0.0), 3.0) + @test n.g[1].opcode == 1 + @test n.g[1].lower == -Inf + @test n.g[1].upper == 5.0 + @test n.h[1].expr == NL._NLExpr(f) + @test n.h[1].opcode == 0 + @test n.h[1].lower == 1.0 + @test n.h[1].upper == 10.0 + @test n.types[1] == [x[1]] # Continuous in both + @test n.types[2] == [x[2]] # Discrete in both + @test n.types[6] == [x[3]] # Discrete in objective only + @test n.types[7] == [x[4]] # Continuous in linear + @test sprint(write, n) == """ + g3 1 1 0 + 4 2 1 1 0 0 + 1 1 + 0 0 + 2 3 2 + 0 0 0 1 + 0 0 1 0 1 + 5 3 + 0 0 + 0 0 0 0 0 + C0 + o0 + n3 + o2 + n2 + o2 + v0 + v1 + C1 + n2 + O0 1 + o2 + v0 + v1 + x4 + 0 0 + 1 0 + 2 0 + 3 0 + r + 1 5 + 0 1 10 + b + 0 0 2 + 0 0 1 + 0 0 2 + 0 0 2 + k3 + 1 + 3 + 4 + J0 2 + 0 1 + 1 0 + J1 3 + 1 1 + 2 1 + 3 1 + G0 3 + 0 0 + 1 0 + 2 1 + """ + return +end + +function test_eval_singlevariable() + x = MOI.VariableIndex(1) + f = NL._NLExpr(MOI.SingleVariable(x)) + @test NL._evaluate(f, Dict(x => 1.2)) == 1.2 +end + +function test_eval_scalaraffine() + x = MOI.VariableIndex.(1:3) + f = NL._NLExpr(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 4.0)) + @test NL._evaluate(f, Dict(x[i] => Float64(i) for i in 1:3)) == 10.0 +end + +function test_eval_scalarquadratic() + x = MOI.VariableIndex(1) + f = MOI.ScalarQuadraticFunction( + [MOI.ScalarAffineTerm(1.1, x)], + [MOI.ScalarQuadraticTerm(2.0, x, x)], + 3.0, + ) + @test NL._evaluate(NL._NLExpr(f), Dict(x => 1.1)) == 5.42 +end + +function test_eval_unary_addition() + x = MOI.VariableIndex(1) + @test NL._evaluate(NL._NLExpr(:(+$x)), Dict(x => 1.1)) == 1.1 +end + +function test_eval_binary_addition() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + @test NL._evaluate(NL._NLExpr(:($x + $y)), Dict(x => 1.1, y => 2.2)) ≈ 3.3 +end + +function test_eval_nary_addition() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + @test NL._evaluate(NL._NLExpr(:($x + $y + 1.0)), Dict(x => 1.1, y => 2.2)) ≈ + 4.3 +end + +function test_eval_unary_subtraction() + x = MOI.VariableIndex(1) + @test NL._evaluate(NL._NLExpr(:(-$x)), Dict(x => 1.1)) == -1.1 +end + +function test_eval_nary_multiplication() + x = MOI.VariableIndex(1) + @test NL._evaluate(NL._NLExpr(:($x * $x * 2.0)), Dict(x => 1.1)) ≈ 2.42 +end + +function test_eval_unary_specialcase() + x = MOI.VariableIndex(1) + @test NL._evaluate(NL._NLExpr(:(cbrt($x))), Dict(x => 1.1)) ≈ 1.1^(1 / 3) +end + +""" + test_issue_79() + +Test the problem + + min (z - 0.5)^2 = z^2 - z + 1/4 + s.t. x * z <= 0 + z ∈ {0, 1} +""" +function test_issue_79() + model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + x = MOI.add_variable(model) + z = MOI.add_variable(model) + MOI.add_constraint(model, MOI.SingleVariable(z), MOI.ZeroOne()) + f = MOI.ScalarQuadraticFunction( + [MOI.ScalarAffineTerm(-1.0, z)], + [MOI.ScalarQuadraticTerm(1.0, z, z)], + 0.25, + ) + MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.add_constraint( + model, + MOI.ScalarQuadraticFunction( + MOI.ScalarAffineTerm{Float64}[], + [MOI.ScalarQuadraticTerm(1.0, x, z)], + 0.0, + ), + MOI.LessThan(0.0), + ) + n = NL.Optimizer() + MOI.copy_to(n, model) + # x is continuous in a nonlinear constraint [Priority 3] + # z is discrete in a nonlinear constraint and objective [Priority 2] + @test n.x[x].order == 1 + @test n.x[z].order == 0 +end + +function test_malformed_constraint_error() + model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + x = MOI.add_variable(model) + MOI.add_constraint( + model, + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 1.0), + MOI.LessThan(0.0), + ) + n = NL.Optimizer() + @test_throws ErrorException MOI.copy_to(n, model) +end + +struct NoExprGraph <: MOI.AbstractNLPEvaluator end +MOI.features_available(::NoExprGraph) = Symbol[] + +function test_noexpr_graph() + model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + block_data = MOI.NLPBlockData(MOI.NLPBoundsPair[], NoExprGraph(), false) + MOI.set(model, MOI.NLPBlock(), block_data) + n = NL.Optimizer() + @test_throws ErrorException MOI.copy_to(n, model) +end + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end +end + +end + +TestNLModel.runtests() diff --git a/test/nl_convert.jl b/test/nl_convert.jl deleted file mode 100644 index 40c5847..0000000 --- a/test/nl_convert.jl +++ /dev/null @@ -1,87 +0,0 @@ -using AmplNLWriter -using Test - -@testset "[nl_convert] check special conversion cases" begin - special_cases = [ - :cbrt, - :abs2, - :inv, - :log2, - :log1p, - :exp2, - :expm1, - :sec, - :csc, - :cot, - :sind, - :cosd, - :tand, - :asind, - :acosd, - :atand, - :secd, - :cscd, - :cotd, - :sech, - :csch, - :coth, - :asech, - :acsch, - ] - for func in special_cases - x = rand() - expr = Expr(:call, func, x) - @test isapprox( - eval(AmplNLWriter.convert_formula(expr)), - eval(expr), - atol = 1e-6, - ) - end - # These functions need input >1 - for func in [:acoth, :asec, :acsc, :acot, :asecd, :acscd, :acotd] - x = rand() + 1 - expr = Expr(:call, func, x) - @test isapprox( - eval(AmplNLWriter.convert_formula(expr)), - eval(expr), - atol = 1e-6, - ) - end -end - -@testset "[nl_convert] check numeric values" begin - x = rand() - @test AmplNLWriter.convert_formula(:($x)) == :($x) - x = -rand() - @test AmplNLWriter.convert_formula(:($x)) == :($x) -end - -@testset "[nl_convert] check unary, binary and n-ary plus" begin - expr = :(+(1)) - @test AmplNLWriter.convert_formula(expr) == :(1) - expr = :(1 + 2) - @test AmplNLWriter.convert_formula(expr) == :(1 + 2) - expr = :(1 + 2 + 3) - @test AmplNLWriter.convert_formula(expr) == :(sum(1, 2, 3)) -end - -@testset "[nl_convert] check unary, binary and n-ary minus" begin - expr = :(-x) - @test AmplNLWriter.convert_formula(expr) == :(neg(x)) - expr = :(x - y) - @test AmplNLWriter.convert_formula(expr) == :(x - y) - expr = :(x - y - z) - @test AmplNLWriter.convert_formula(expr) == :((x - y) - z) -end - -@testset "[nl_convert] check n-ary multiplication" begin - expr = :(x * y * z) - @test AmplNLWriter.convert_formula(expr) == :(x * (y * z)) -end - -@testset "[nl_convert] check comparison expansion" begin - expr = :(1 < 2 < 3) - @test AmplNLWriter.convert_formula(expr) == :(1 < 2 && 2 < 3) - expr = :(1 < 2 < 3 < 4) - @test AmplNLWriter.convert_formula(expr) == :((1 < 2 && 2 < 3) && 3 < 4) -end diff --git a/test/nl_linearity.jl b/test/nl_linearity.jl deleted file mode 100644 index 2663443..0000000 --- a/test/nl_linearity.jl +++ /dev/null @@ -1,24 +0,0 @@ -using AmplNLWriter -using Test - -@testset "[nl_linearity] check simplification of formulae" begin - @testset "ifelse" begin - # First term is true, we should choose `then` - expr = :(ifelse(1 > 0, x[1]^2, x[2] + 1)) - lin_expr = AmplNLWriter.LinearityExpr(expr) - @test AmplNLWriter.convert_formula(lin_expr.c) == :(x[1]^2) - @test lin_expr.linearity == :nonlinear - - # First term is false, we should choose `else` - expr = :(ifelse(1 < 0, x[1]^2, x[2] + 1)) - lin_expr = AmplNLWriter.LinearityExpr(expr) - @test AmplNLWriter.convert_formula(lin_expr.c) == :(x[2] + 1) - @test lin_expr.linearity == :linear - - # First term isn't constant, we can't simplify - expr = :(ifelse(1 > x[1], x[1]^2, x[2] + 1)) - lin_expr = AmplNLWriter.LinearityExpr(expr) - @test AmplNLWriter.convert_formula(lin_expr.c) == expr - @test lin_expr.linearity == :nonlinear - end -end diff --git a/test/nl_write.jl b/test/nl_write.jl deleted file mode 100644 index a0bc349..0000000 --- a/test/nl_write.jl +++ /dev/null @@ -1,50 +0,0 @@ -using AmplNLWriter -using Test - -function test_temp_file_handling(path) - # Turn on debug mode so files persist - old_debug = AmplNLWriter.debug - AmplNLWriter.setdebug(true) - - filename = "test" - filepath = joinpath(AmplNLWriter.solverdata_dir, filename) - probfile = "$filepath.nl" - solfile = "$filepath.sol" - AmplNLWriter.clean_solverdata() - - @testset "all temp files deleted successfully" begin - @test length(readdir(AmplNLWriter.solverdata_dir)) == 1 - end - - MOI = AmplNLWriter.MathOptInterface - solver = AmplNLWriter.Optimizer(path, filename = filename) - x = MOI.add_variable(solver) - MOI.add_constraint( - solver, - MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0), - MOI.GreaterThan(0.0), - ) - MOI.set( - solver, - MOI.ObjectiveFunction{MOI.SingleVariable}(), - MOI.SingleVariable(x), - ) - MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE) - MOI.optimize!(solver) - - @testset "temp files present after solve in debug mode" begin - @test length(readdir(AmplNLWriter.solverdata_dir)) == 3 - end - @testset "temp files used custom name" begin - @test isfile(probfile) == true - @test isfile(solfile) == true - end - - # Reset debug mode and clean up - AmplNLWriter.setdebug(old_debug) - return AmplNLWriter.clean_solverdata() -end - -run_with_ampl() do path - return test_temp_file_handling(path) -end diff --git a/test/runtests.jl b/test/runtests.jl index cbf06a0..3cc71d5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,20 +1,9 @@ using Test -if VERSION < v"1.3" - import Ipopt - run_with_ampl(f) = f(Ipopt.amplexe) -else - import Ipopt_jll - run_with_ampl(f) = Ipopt_jll.amplexe(f) +@testset "NLModel" begin + include("NLModel.jl") end @testset "MOI" begin include("MOI_wrapper.jl") end - -@testset "Base" begin - include("nl_convert.jl") - include("nl_linearity.jl") - include("nl_write.jl") - include("sol_file_parser.jl") -end diff --git a/test/sol_file_parser.jl b/test/sol_file_parser.jl deleted file mode 100644 index 8803dfc..0000000 --- a/test/sol_file_parser.jl +++ /dev/null @@ -1,32 +0,0 @@ -using AmplNLWriter -using Test - -const FILE_80 = """ - - Couenne (/home/bgulcan/.julia/v0.6/AmplNLWriter/.solverdata/tmpMSxFvq.nl Mar 6 2018): Infeasible - -Options -3 -0 -1 -0 -140 -0 -110 -0 -objno 0 220 -""" - -@testset "Test sol file parsing" begin - @testset "Issue #80" begin - io = IOBuffer() - write(io, FILE_80) - seekstart(io) - model = AmplNLWriter.AmplNLMathProgModel("", String[], "") - model.ncon = 140 - model.nvar = 110 - @test !AmplNLWriter.read_sol(io, model) - @test model.solve_result_num == 220 - @test all(model.solution .=== NaN) - end -end