From 2179beacec39b746bfcf6e9215d6093e8321ec3d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 16 Apr 2026 08:03:53 +0200 Subject: [PATCH] Remove ExaModels wrapper --- Project.toml | 2 - src/GenOpt.jl | 1 - src/examodels.jl | 1045 ---------------------------------------------- 3 files changed, 1048 deletions(-) delete mode 100644 src/examodels.jl diff --git a/Project.toml b/Project.toml index 57646b7..953e2b0 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,10 @@ version = "0.1.1" authors = ["Benoît Legat "] [deps] -ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" [compat] -ExaModels = "0.9" JuMP = "1.29" MathOptInterface = "1.46" julia = "1.10" diff --git a/src/GenOpt.jl b/src/GenOpt.jl index bb8e003..0d69dd0 100644 --- a/src/GenOpt.jl +++ b/src/GenOpt.jl @@ -8,7 +8,6 @@ module GenOpt include("MOI_wrapper.jl") include("bridge.jl") include("JuMP_wrapper.jl") -include("examodels.jl") # Copied from JuMP.jl: const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include] diff --git a/src/examodels.jl b/src/examodels.jl deleted file mode 100644 index 0deacc9..0000000 --- a/src/examodels.jl +++ /dev/null @@ -1,1045 +0,0 @@ -# Copyright (c) 2023 Sungho Shin -# -# See licence file in ExaModels.jl -# Inspired from ExaModels/ext/ExaModelsMOI.jl - -# Copyright (c) 2025: Benoît Legat and contributors -# -# Use of this source code is governed by an MIT-style license that can be found -# in the LICENSE.md file or at https://opensource.org/licenses/MIT. - -import ExaModels - -import MathOptInterface - -const SUPPORTED_FUNC_TYPE{T} = Union{ - MOI.ScalarAffineFunction{T}, - MOI.ScalarQuadraticFunction{T}, - MOI.ScalarNonlinearFunction, -} -const SUPPORTED_FUNC_TYPE_WITH_VAR{T} = - Union{SUPPORTED_FUNC_TYPE{T},MOI.VariableIndex} -const SUPPORTED_FUNC_SET_TYPE{T} = - Union{MOI.GreaterThan{T},MOI.LessThan{T},MOI.EqualTo{T},MOI.Interval{T}} -const SUPPORTED_VAR_SET_TYPE{T} = - Union{MOI.GreaterThan{T},MOI.LessThan{T},MOI.EqualTo{T},MOI.Parameter{T}} -const PARAMETER_INDEX_THRESHOLD = Int64(4_611_686_018_427_387_904) # div(typemax(Int64),2)+1 - -""" - Abstract data structure for storing expression tree and data arrays -""" -abstract type AbstractBin end - -# The linked list with `inner` represents a sum of expressions `∑_{i in data} head(i)` -# It is a linked list and not a `Vector` as each `head` may have a different type. -# We append new ones at the beginning because `Bin` is non-mutable and -# its fields are concretely typed. -struct Bin{E,P,I} <: AbstractBin - head::E - data::P - inner::I -end - -struct BinNull <: AbstractBin end - -function update_bin!(bin, e, p) - if _update_bin!(bin, e, p) # if update succeeded, return the original bin - return bin - else # if update has failed, return a new bin - if p isa Tuple - p = [p] - end - return Bin(e, p, bin) - end -end -function _update_bin!(bin::Bin{E,P,I}, e, p) where {E,P,I} - if e == bin.head && p isa eltype(bin.data) - push!(bin.data, p) - return true - else - return _update_bin!(bin.inner, e, p) - end -end -function _update_bin!(::BinNull, e, p) - return false -end - -function check_supported(T, moim) - con_types = MOI.get(moim, MOI.ListOfConstraintTypesPresent()) - for (F, S) in con_types - if F <: FunctionGenerator - continue - end - !(F <: SUPPORTED_FUNC_TYPE_WITH_VAR) && - error("Unsupported function type $F.") - if F <: MOI.VariableIndex - !(S <: SUPPORTED_VAR_SET_TYPE) && - error("Unsupported variable index constraint $F in $S") - else - !(S <: SUPPORTED_FUNC_SET_TYPE) && error("Unsupported set type $S") - end - end - - obj_type = MOI.get(moim, MOI.ObjectiveFunctionType()) - !(obj_type <: SUPPORTED_FUNC_TYPE_WITH_VAR) && - error("Unsupported objective function type $obj_type.") - - obj_sense = MOI.get(moim, MOI.ObjectiveSense()) - !(obj_sense in (MOI.MIN_SENSE, MOI.MAX_SENSE)) && - error("Unsupported objective sense $obj_sense.") - return obj_sense === MOI.MIN_SENSE -end - -function to_exacore(moim::MOI.ModelLike; backend = nothing, T = Float64) - minimize = check_supported(T, moim) - - c = ExaModels.ExaCore(T; backend = backend, minimize = minimize) - - var_to_idx = copy_variables!(c, moim, T) - con_to_idx = copy_constraints!(c, moim, var_to_idx, T) - copy_objective!(c, moim, var_to_idx) - - return c, (var_to_idx, con_to_idx) -end - -function fill_variable_bounds!(moim, lvar, uvar, var_to_idx, T) - for ci in MOI.get( - moim, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.GreaterThan{T}}(), - ) - vi = MOI.get(moim, MOI.ConstraintFunction(), ci) - lvar[var_to_idx[vi]] = MOI.get(moim, MOI.ConstraintSet(), ci).lower - end - for ci in MOI.get( - moim, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.LessThan{T}}(), - ) - vi = MOI.get(moim, MOI.ConstraintFunction(), ci) - uvar[var_to_idx[vi]] = MOI.get(moim, MOI.ConstraintSet(), ci).upper - end - for ci in MOI.get( - moim, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.EqualTo{T}}(), - ) - vi = MOI.get(moim, MOI.ConstraintFunction(), ci) - fixed_val = MOI.get(moim, MOI.ConstraintSet(), ci).value - lvar[var_to_idx[vi]] = fixed_val - uvar[var_to_idx[vi]] = fixed_val - end -end - -function fill_variable_start!(moim, x0, param_vis) - var_to_idx = Dict{MOI.VariableIndex,Int}() - i = 0 - for vi in MOI.get(moim, MOI.ListOfVariableIndices()) - vi ∈ param_vis && continue - i += 1 - var_to_idx[vi] = i - start = if MOI.supports(moim, MOI.VariablePrimalStart(), typeof(vi)) - MOI.get(moim, MOI.VariablePrimalStart(), vi) - else - nothing - end - isnothing(start) && continue - x0[i] = start - end - return var_to_idx -end - -function _get_parameters(moim::MOI.ModelLike, T) - cis = MOI.get( - moim, - MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Parameter{T}}(), - ) - parameters = Vector{Tuple{MOI.VariableIndex,MOI.Parameter{T}}}() - for ci in cis - vi = MOI.get(moim, MOI.ConstraintFunction(), ci) - set = MOI.get(moim, MOI.ConstraintSet(), ci) - push!(parameters, (vi, set)) - end - sort!(parameters, by = x -> x[1].value) - return parameters -end - -function copy_variables!(c, moim, T) - nvarpar = MOI.get(moim, MOI.NumberOfVariables()) - parameters = _get_parameters(moim, T) - npar = length(parameters) - nvar = nvarpar - npar - - x0 = zeros(T, nvar) - var_to_idx = fill_variable_start!(moim, x0, first.(parameters)) - - lvar = fill(T(-Inf), nvar) - uvar = fill(T(Inf), nvar) - fill_variable_bounds!(moim, lvar, uvar, var_to_idx, T) - - ExaModels.variable(c, nvar; start = x0, lvar = lvar, uvar = uvar) - - varpar_to_idx = Dict() - for (vi, i) in var_to_idx - varpar_to_idx[vi] = (type = :variable, idx = i) - end - - if npar > 0 - p0 = zeros(T, npar) - for (i, (vi, set)) in enumerate(parameters) - p0[i] = T(set.value) - varpar_to_idx[vi] = (type = :parameter, idx = i) - end - ExaModels.parameter(c, p0) - end - - return varpar_to_idx -end - -function copy_objective!(c, moim, var_to_idx) - obj_type = MOI.get(moim, MOI.ObjectiveFunctionType()) - - bin = BinNull() - bin = exafy_obj( - MOI.get(moim, MOI.ObjectiveFunction{obj_type}()), - bin, - var_to_idx, - ) - - return build_objective!(c, bin) -end - -function copy_constraints!(c, moim, var_to_idx, T) - bin = BinNull() - offset = 0 - lcon = zeros(T, 0) - ucon = zeros(T, 0) - y0 = zeros(T, 0) - con_to_idx = Dict{MOI.ConstraintIndex,Int}() - - con_types = MOI.get(moim, MOI.ListOfConstraintTypesPresent()) - for (F, S) in con_types - F <: MOI.VariableIndex && continue - F <: FunctionGenerator && continue - cis = MOI.get(moim, MOI.ListOfConstraintIndices{F,S}()) - bin, offset = exafy_con( - moim, - cis, - bin, - offset, - lcon, - ucon, - y0, - var_to_idx, - con_to_idx, - ) - end - cons = ExaModels.constraint(c, offset; start = y0, lcon = lcon, ucon = ucon) - build_constraint!(c, cons, bin) - - for (F, S) in con_types - F <: FunctionGenerator || continue - cis = MOI.get(moim, MOI.ListOfConstraintIndices{F,S}()) - copy_generator_constraints!(c, moim, cis, var_to_idx, con_to_idx, T) - end - - return con_to_idx -end - -exagen(α::Number, _) = α - -#exagen(i::IteratorIndex) = ExaModels.ParIndexed(ExaModels.ParSource(), i.value) - -function exagen(f::MOI.ScalarNonlinearFunction, offsets) - if f.head == :getindex - v = f.args[1] - if v isa ContiguousArrayOfVariables - idx = exagen(f.args[2], offsets) - if !iszero(v.offset) - idx = v.offset + idx - end - cp = cumprod(v.size) - for i in 3:length(f.args) - idx += cp[i-2] * (exagen(f.args[i], offsets) - 1) - end - return ExaModels.Var(idx) - elseif v isa IteratorIndex - @assert length(f.args) == 2 - @assert f.args[2] isa Integer - if isnothing(offsets) - @assert isone(f.args[2]) - return ExaModels.ParSource() - else - return ExaModels.ParIndexed( - ExaModels.ParSource(), - offsets[v.value] + f.args[2], - ) - end - else - error( - "Unexpected the first operatand of `getindex` to be of type `$(typeof(v))`", - ) - end - else - return op(f.head)((exagen(e, offsets) for e in f.args)...) - end -end - -_lower_bounds(::Union{MOI.Zeros,MOI.Nonnegatives}, T) = zero(T) -_lower_bounds(::MOI.Nonpositives, T) = typemin(T) -_upper_bounds(::Union{MOI.Zeros,MOI.Nonpositives}, T) = zero(T) -_upper_bounds(::MOI.Nonnegatives, T) = typemax(T) - -function _exagen(func::MOI.ScalarNonlinearFunction, iterators) - lengths = map(it -> length(first(it.values)), iterators) - if length(lengths) == 1 && lengths[] == 1 - cs = nothing - pars = only.(iterators[].values) - else - cs = [0; cumsum(lengths)[1:(end-1)]] - pars = vec( - map( - Base.Iterators.ProductIterator( - ntuple(i -> iterators[i].values, length(iterators)), - ), - ) do I - return reduce((i, j) -> tuple(i..., j...), I) - end, - ) - end - expr = exagen(func, cs) - return expr, pars -end - -function _exafy(func::SumGenerator, _) - return _exagen(func.func, func.iterators) -end - -function copy_generator_constraints!(c, moim, cis, var_to_idx, con_to_idx, T) - # FIXME we assume that `var_to_idx` is the identity - for ci in cis - func = MOI.get(moim, MOI.ConstraintFunction(), ci) - set = MOI.get(moim, MOI.ConstraintSet(), ci) - con_to_idx[ci] = c.ncon - expr, pars = _exagen(func.func, func.iterators) - ExaModels.constraint( - c, - expr, - pars; - lcon = _lower_bounds(set, T), - ucon = _upper_bounds(set, T), - ) - end -end - -function _exafy_con( - i, - c::C, - bin, - var_to_idx, - con_to_idx; - pos = true, -) where {C<:MOI.ScalarAffineFunction} - for mm in c.terms - e, p = _exafy(mm, var_to_idx) - e = pos ? e : -e - bin = update_bin!( - bin, - ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) => e, - (p..., con_to_idx[i]), - ) # augment data with constraint index - end - bin = update_bin!(bin, ExaModels.Null(c.constant), (1,)) - return bin -end -function _exafy_con( - i, - c::C, - bin, - var_to_idx, - con_to_idx; - pos = true, -) where {C<:MOI.ScalarQuadraticFunction} - for mm in c.affine_terms - e, p = _exafy(mm, var_to_idx) - e = pos ? e : -e - bin = update_bin!( - bin, - ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) => e, - (p..., con_to_idx[i]), - ) # augment data with constraint index - end - for mm in c.quadratic_terms - e, p = _exafy(mm, var_to_idx) - e = pos ? e : -e - bin = update_bin!( - bin, - ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) => e, - (p..., con_to_idx[i]), - ) # augment data with constraint index - end - bin = update_bin!(bin, ExaModels.Null(c.constant), (1,)) - return bin -end -function _exafy_con( - i, - c::C, - bin, - var_to_idx, - con_to_idx; - pos = true, -) where {C<:MOI.ScalarNonlinearFunction} - if c.head == :+ - for mm in c.args - bin = _exafy_con(i, mm, bin, var_to_idx, con_to_idx) - end - # elseif c.head == :- - # bin, offset = _exafy_con(i, c.args[1], bin, offset) - # bin, offset = _exafy_con(i, c.args[2], bin, offset; pos = false) - else - e, p = _exafy(c, var_to_idx) - e = pos ? e : -e - bin = update_bin!( - bin, - ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) => e, - (p..., con_to_idx[i]), - ) # augment data with constraint index - end - return bin -end -function _exafy_con( - i, - c::C, - bin, - var_to_idx, - con_to_idx; - pos = true, -) where {C<:Real} - e = - pos ? ExaModels.ParIndexed(ExaModels.ParSource(), 1) : - -ExaModels.ParIndexed(ExaModels.ParSource(), 1) - bin = update_bin!( - bin, - ExaModels.ParIndexed(ExaModels.ParSource(), 2) => - 0 * ExaModels.Var(1) + e, - (c, con_to_idx[i]), - ) - - return bin -end - -function exafy_con( - moim, - cons::V, - bin, - offset, - lcon, - ucon, - y0, - var_to_idx, - con_to_idx, -) where {V<:Vector{<:MOI.ConstraintIndex}} - l = sum(cons) do ci - return MOI.dimension(MOI.get(moim, MOI.ConstraintSet(), ci)) - end - - resize!(lcon, offset + l) - resize!(ucon, offset + l) - resize!(y0, offset + l) - for (i, ci) in enumerate(cons) - func = MOI.get(moim, MOI.ConstraintFunction(), ci) - set = MOI.get(moim, MOI.ConstraintSet(), ci) - con_to_idx[ci] = offset + 1 - start = if MOI.supports(moim, MOI.ConstraintPrimalStart(), typeof(ci)) - MOI.get(moim, MOI.ConstraintPrimalStart(), ci) - else - nothing - end - _exafy_con_update_start(ci, start, y0, con_to_idx) - _exafy_con_update_vector(ci, set, lcon, ucon, con_to_idx) - bin = _exafy_con(ci, func, bin, var_to_idx, con_to_idx) - offset += MOI.dimension(set) - end - return bin, offest -end - -function _exafy_con_update_start(i, start, y0, con_to_idx) - return y0[con_to_idx[i]] = start -end - -function _exafy_con_update_start(i, ::Nothing, y0, con_to_idx) - return y0[con_to_idx[i]] = zero(eltype(y0)) -end - -function _exafy_con_update_vector( - i, - e::MOI.Interval{T}, - lcon, - ucon, - con_to_idx, -) where {T} - lcon[con_to_idx[i]] = e.lower - return ucon[con_to_idx[i]] = e.upper -end - -function _exafy_con_update_vector( - i, - e::MOI.LessThan{T}, - lcon, - ucon, - con_to_idx, -) where {T} - lcon[con_to_idx[i]] = -Inf - return ucon[con_to_idx[i]] = e.upper -end - -function _exafy_con_update_vector( - i, - e::MOI.GreaterThan{T}, - lcon, - ucon, - con_to_idx, -) where {T} - ucon[con_to_idx[i]] = Inf - return lcon[con_to_idx[i]] = e.lower -end - -function _exafy_con_update_vector( - i, - e::MOI.EqualTo{T}, - lcon, - ucon, - con_to_idx, -) where {T} - lcon[con_to_idx[i]] = e.value - return ucon[con_to_idx[i]] = e.value -end - -function build_constraint!(c, cons, bin) - build_constraint!(c, cons, bin.inner) - return ExaModels.constraint!(c, cons, bin.head, bin.data) -end - -function build_constraint!(c, cons, ::BinNull) end - -function build_objective!(c, bin) - build_objective!(c, bin.inner) - return ExaModels.objective(c, bin.head, bin.data) -end - -function build_objective!(c, ::BinNull) end - -function exafy_obj(o::Nothing, bin, var_to_idx) - return bin -end - -function exafy_obj(o::MOI.VariableIndex, bin, var_to_idx) - e, p = _exafy(o, var_to_idx) - return update_bin!(bin, e, p) -end - -function exafy_obj(o::MOI.ScalarQuadraticFunction{T}, bin, var_to_idx) where {T} - for m in o.affine_terms - e, p = _exafy(m, var_to_idx) - bin = update_bin!(bin, e, p) - end - for m in o.quadratic_terms - e, p = _exafy(m, var_to_idx) - bin = update_bin!(bin, e, p) - end - - return update_bin!(bin, ExaModels.Null(o.constant), (1,)) -end - -function exafy_obj(o::MOI.ScalarAffineFunction{T}, bin, var_to_idx) where {T} - for m in o.terms - e, p = _exafy(m, var_to_idx) - bin = update_bin!(bin, e, p) - end - - return update_bin!(bin, ExaModels.Null(o.constant), (1,)) -end - -function exafy_obj(o::MOI.ScalarNonlinearFunction, bin, var_to_idx) - constant = 0.0 - if o.head == :+ - for m in o.args - if m isa MOI.ScalarAffineFunction - for mm in m.terms - e, p = _exafy(mm, var_to_idx) - bin = update_bin!(bin, e, p) - end - elseif m isa MOI.ScalarQuadraticFunction - for mm in m.affine_terms - e, p = _exafy(mm, var_to_idx) - bin = update_bin!(bin, e, p) - end - for mm in m.quadratic_terms - e, p = _exafy(mm, var_to_idx) - bin = update_bin!(bin, e, p) - end - constant += m.constant - else - e, p = _exafy(m, var_to_idx) - bin = update_bin!(bin, e, p) - end - end - else - e, p = _exafy(o, var_to_idx) - bin = update_bin!(bin, e, p) - end - - return update_bin!(bin, ExaModels.Null(constant), (1,)) # TODO see if this can be empty tuple -end - -function _exafy(v::MOI.VariableIndex, var_to_idx, p = ()) - i = ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) - vartype, idx = var_to_idx[v] - if vartype === :variable - return ExaModels.Var(i), (p..., idx) - elseif vartype === :parameter - return ExaModels.ParameterNode(i), (p..., idx) - else - error("Unknown variable type: $vartype") - end -end - -function _exafy(i::R, var_to_idx, p) where {R<:Real} - return ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1), (p..., i) -end - -function _exafy(e::MOI.ScalarNonlinearFunction, var_to_idx, p = ()) - args = map(e.args) do expr - c, p = _exafy(expr, var_to_idx, p) - return c - end - return op(e.head)(args...), p -end - -function _exafy(e::MOI.ScalarAffineFunction{T}, var_to_idx, p = ()) where {T} - ec = if !isempty(e.terms) - sum(begin - c1, p = _exafy(term, var_to_idx, p) - c1 - end for term in e.terms) + - ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) - else - ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) - end - - return ec, (p..., e.constant) -end - -function _exafy(e::MOI.ScalarAffineTerm{T}, var_to_idx, p = ()) where {T} - c1, p = _exafy(e.variable, var_to_idx, p) - return *(c1, ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1)), - (p..., e.coefficient) -end - -function _exafy(e::MOI.ScalarQuadraticFunction{T}, var_to_idx, p = ()) where {T} - t = ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) - p = (p..., e.constant) - - if !isempty(e.affine_terms) - t += sum( - begin - c1, p = _exafy(term, var_to_idx, p) - c1 - end for term in e.affine_terms - ) - end - - if !isempty(e.quadratic_terms) - t += sum( - begin - c1, p = _exafy(term, var_to_idx, p) - c1 - end for term in e.quadratic_terms - ) - end - - return t, p -end - -function _exafy(e::MOI.ScalarQuadraticTerm{T}, var_to_idx, p = ()) where {T} - if e.variable_1 == e.variable_2 - v, p = _exafy(e.variable_1, var_to_idx, p) - return ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) * - abs2(v), - (p..., e.coefficient / 2) # it seems that MOI assumes this by default - else - v1, p = _exafy(e.variable_1, var_to_idx, p) - v2, p = _exafy(e.variable_2, var_to_idx, p) - return ExaModels.ParIndexed(ExaModels.ParSource(), length(p) + 1) * - v1 * - v2, - (p..., e.coefficient) - end -end - -# eval can be a performance killer -- we want to explicitly include symbols for frequently used operations. -function op(s::Symbol) - # uni/multi - if s === :+ - return + - elseif s === :- - return - - # multi - elseif s === :* - return * - elseif s === :^ - return ^ - elseif s === :/ - return / - # uni - elseif s === :abs - return abs - elseif s === :sign - error("sign not supported") - elseif s === :sqrt - return sqrt - elseif s === :cbrt - return cbrt - elseif s === :abs2 - return abs2 - elseif s === :inv - return inv - elseif s === :log - return log - elseif s === :log10 - return log10 - elseif s === :log2 - return log2 - elseif s === :log1p - return log1p - elseif s === :exp - return exp - elseif s === :exp2 - return exp2 - elseif s === :expm1 - error("expm1 not supported") - # trig - elseif s === :sin - return sin - elseif s === :cos - return cos - elseif s === :tan - return tan - elseif s === :sec - return sec - elseif s === :csc - return csc - elseif s === :cot - return cot - elseif s === :sind - return sind - elseif s === :cosd - return cosd - elseif s === :tand - return tand - elseif s === :secd - return secd - elseif s === :cscd - return cscd - elseif s === :cotd - return cotd - elseif s === :asin - return asin - elseif s === :acos - return acos - elseif s === :atan - return atan - elseif s === :asec - error("asec not supported") - elseif s === :acsc - error("acsc not supported") - elseif s === :acot - return acot - elseif s === :asind - error("asind not supported") - elseif s === :acosd - error("acosd not supported") - elseif s === :atand - return atand - elseif s === :asecd - error("aced not supported") - elseif s === :acscd - error("acscd not supported") - elseif s === :acotd - return acotd - elseif s === :sinh - return sinh - elseif s === :cosh - return cosh - elseif s === :tanh - return tanh - elseif s === :sech - return sech - elseif s === :csch - return csch - elseif s === :coth - return coth - elseif s === :asinh - return asinh - elseif s === :acosh - return acosh - elseif s === :atanh - return atanh - elseif s === :asech - error("asech not supported") - elseif s === :acsch - error("acsch not supported") - elseif s === :acoth - return acoth - # special (commented will use `eval` which would succeed if SpecialFunctions is loaded) - elseif s === :deg2rad - error("deg2rad not supported") - elseif s === :rad2deg - error("rad2deg not supported") - # elseif s === :erf error("erf not supported") - # elseif s === :erfinv error("erfinv not supported") - # elseif s === :erfc error("erfc not supported") - # elseif s === :erfcinv error("erfcinv not supported") - # elseif s === :erfi error("erfi not supported") - # elseif s === :gamma error("gamma not supported") - elseif s === :lgamma - error("lgamma not supported") - # elseif s === :digamma error("digamma not supported") - # elseif s === :invdigamma error("invdigamma not supported") - # elseif s === :trigamma error("trigamma not supported") - # elseif s === :airyai error("airyai not supported") - # elseif s === :airybi error("airybi not supported") - # elseif s === :airyaiprime error("airyaiprime not supported") - # elseif s === :airybiprime error("airybiprime not supported") - # elseif s === :besselj0 error("besselj0 not supported") - # elseif s === :besselj1 error("besselj1 not supported") - # elseif s === :bessely0 error("bessely0 not supported") - # elseif s === :bessely1 error("bessely1 not supported") - # elseif s === :erfcx error("erfcx not supported") - # elseif s === :dawson error("dawson not supported") - - # not in MOI - elseif s === :exp10 - return exp10 - elseif s === :beta - return beta - elseif s === :logbeta - return logbeta - else - return eval(s) - end -end - -mutable struct ExaOptimizer{B,S} <: MOI.AbstractOptimizer - solver::S - backend::B - model::Union{Nothing,ExaModels.ExaModel} - result::Any - solve_time::Float64 - options::Dict{Symbol,Any} -end - -MOI.is_empty(model::ExaOptimizer) = isnothing(model.model) - -function MOI.supports_constraint( - ::ExaOptimizer, - ::Type{<:FunctionGenerator}, - ::Type{<:Union{MOI.Zeros,MOI.Nonnegatives,MOI.Nonpositives}}, -) - return true -end - -function MOI.supports_constraint( - ::ExaOptimizer, - ::Type{<:SUPPORTED_FUNC_TYPE}, - ::Type{<:SUPPORTED_FUNC_SET_TYPE}, -) - return true -end -function MOI.supports_constraint( - ::ExaOptimizer, - ::Type{MOI.VariableIndex}, - ::Type{<:SUPPORTED_VAR_SET_TYPE}, -) - return true -end -function MOI.supports(::ExaOptimizer, ::MOI.ObjectiveSense) - return true -end -function MOI.supports( - ::ExaOptimizer, - ::MOI.ObjectiveFunction{<:SUPPORTED_FUNC_TYPE_WITH_VAR}, -) - return true -end -function MOI.supports( - ::ExaOptimizer, - ::MOI.VariablePrimalStart, - ::Type{MOI.VariableIndex}, -) - return true -end - -function ExaOptimizer(solver, backend = nothing; kwargs...) - return ExaOptimizer( - solver, - backend, - nothing, - nothing, - 0.0, - Dict{Symbol,Any}(kwargs...), - ) -end - -function MOI.empty!(model::ExaOptimizer) - return model.model = nothing -end - -function MOI.copy_to(dest::ExaOptimizer, src::MOI.ModelLike) - core, maps = to_exacore(src; backend = dest.backend) - dest.model = ExaModels.ExaModel(core; prod = true) - return _make_index_map(src, maps) -end - -function MOI.optimize!(optimizer::ExaOptimizer) - optimizer.solve_time = @elapsed begin - result = optimizer.solver(optimizer.model; optimizer.options...) - optimizer.result = ( - objective = result.objective, - solution = Array(result.solution), - multipliers = Array(result.multipliers), - multipliers_L = Array(result.multipliers_L), - multipliers_U = Array(result.multipliers_U), - status = result.status, - ) - end - - return optimizer -end - -function MOI.get(optimizer::ExaOptimizer, ::MOI.TerminationStatus) - return ExaModels.termination_status_translator( - optimizer.solver, - optimizer.result.status, - ) -end -function MOI.get( - model::ExaOptimizer, - attr::Union{MOI.PrimalStatus,MOI.DualStatus}, -) - return ExaModels.result_status_translator(model.solver, model.result.status) -end - -function MOI.get( - model::ExaOptimizer, - attr::MOI.VariablePrimal, - vi::MOI.VariableIndex, -) - MOI.check_result_index_bounds(model, attr) - if vi.value > PARAMETER_INDEX_THRESHOLD - return model.model.θ[vi.value-PARAMETER_INDEX_THRESHOLD] - else - return model.result.solution[vi.value] - end -end - -function MOI.get( - model::ExaOptimizer, - attr::MOI.ConstraintDual, - ci::MOI.ConstraintIndex{<:SUPPORTED_FUNC_TYPE,<:SUPPORTED_FUNC_SET_TYPE}, -) - MOI.check_result_index_bounds(model, attr) - # MOI.throw_if_not_valid(model, ci) - s = -1.0 - return s * model.result.multipliers[ci.value] -end - -function MOI.get( - model::ExaOptimizer, - attr::MOI.ConstraintDual, - ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}, -) - MOI.check_result_index_bounds(model, attr) - # MOI.throw_if_not_valid(model, ci) - rc = - model.result.multipliers_L[ci.value] - - model.result.multipliers_U[ci.value] - return min(0.0, rc) -end - -function MOI.get( - model::ExaOptimizer, - attr::MOI.ConstraintDual, - ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}, -) - MOI.check_result_index_bounds(model, attr) - # MOI.throw_if_not_valid(model, ci) - rc = - model.result.multipliers_L[ci.value] - - model.result.multipliers_U[ci.value] - return max(0.0, rc) -end - -function MOI.get( - model::ExaOptimizer, - attr::MOI.ConstraintDual, - ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}, -) - MOI.check_result_index_bounds(model, attr) - # MOI.throw_if_not_valid(model, ci) - rc = - model.result.multipliers_L[ci.value] - - model.result.multipliers_U[ci.value] - return rc -end - -function MOI.get(model::ExaOptimizer, ::MOI.ResultCount) - return (model.result !== nothing) ? 1 : 0 -end - -function MOI.get(model::ExaOptimizer, attr::MOI.ObjectiveValue) - MOI.check_result_index_bounds(model, attr) - # scale = (model.sense == MOI.MAX_SENSE) ? -1 : 1 - # return scale * model.result.objective - return model.result.objective -end - -MOI.get(model::ExaOptimizer, ::MOI.SolveTimeSec) = model.solve_time -function MOI.get(model::ExaOptimizer, ::MOI.SolverName) - return "$(string(model.solver)) running with ExaModels" -end - -function MOI.set(model::ExaOptimizer, p::MOI.RawOptimizerAttribute, value) - model.options[Symbol(p.name)] = value - # No need to reset model.solver because this gets handled in optimize!. - return -end - -function _make_index_map(model::MOI.ModelLike, maps) - return _make_index_map(model, maps[1], maps[2]) -end -function _make_index_map(model::MOI.ModelLike, var_to_idx, con_to_idx) - variables = MOI.get(model, MOI.ListOfVariableIndices()) - map = MOI.Utilities.IndexMap() - for x in variables - vartype, rawidx = var_to_idx[x] - if vartype === :variable - map[x] = typeof(x)(rawidx) - elseif vartype === :parameter - map[x] = typeof(x)(rawidx + PARAMETER_INDEX_THRESHOLD) - else - error("Unknown variable type $vartype") - end - end - for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent()) - _make_constraints_map(model, map.con_map[F, S], con_to_idx) - end - return map -end -function _make_constraints_map( - model, - map::MOI.Utilities.DoubleDicts.IndexDoubleDictInner{F,S}, - con_to_idx, -) where {F,S} - for c in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) - map[c] = typeof(c)(con_to_idx[c]) - end - return -end