From 6788a5d62912e41c34e69acf7e87f86b9295085d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Thu, 21 Sep 2023 23:09:48 +0200 Subject: [PATCH] First version passing basic tests --- Project.toml | 8 +- README.md | 159 ++++++++++++++++++ src/PRIMA.jl | 429 ++++++++++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 136 ++++++++++++++- 4 files changed, 728 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 45d39cd..3d18c99 100644 --- a/Project.toml +++ b/Project.toml @@ -3,8 +3,14 @@ uuid = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed" authors = ["Éric Thiébaut and contributors"] version = "0.1.0" +[deps] +PRIMA_jll = "eead6e0c-2d5b-5641-a95c-b722de96d551" +TypeUtils = "c3b1956e-8857-4d84-9b79-890df85b1e67" + [compat] -julia = "1" +julia = "1.6" +PRIMA_jll = "v0.7.0+1" +TypeUtils = "0.3" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index ec53396..3e6df2f 100644 --- a/README.md +++ b/README.md @@ -1 +1,160 @@ # PRIMA [![Build Status](https://github.com/emmt/PRIMA.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/emmt/PRIMA.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Build Status](https://ci.appveyor.com/api/projects/status/github/emmt/PRIMA.jl?svg=true)](https://ci.appveyor.com/project/emmt/PRIMA-jl) [![Coverage](https://codecov.io/gh/emmt/PRIMA.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/emmt/PRIMA.jl) + +This package is a Julia interface to [PRIMA](https://github.com/libprima/prima) +a **R**eference **I**mplementation for **P**owell's **M**ethods with +**M**odernization and **A**melioration which implement algorithms by M.J.D. +Powell for minimizing a multi-variate objective function possibly under +constraints and without derivatives. + +Formally, these algorithms are designed to solve problems of the form: + +``` julia +min f(x) s.t. x ∈ Ω ⊆ ℝⁿ +``` + +where `f(x)` is the function to minimize and `Ω ⊆ ℝⁿ` is the set of feasible +variables. + +Five algorithms are provided: + +- `uobyqa` for *Unconstrained Optimization BY Quadratic Approximations* is for + unconstrained optimization, that is `Ω = ℝⁿ`; + +- `newuoa` for *Unconstrained Optimization BY Quadratic Approximations* is also + for unconstrained optimization (according to M.J.D. Powell, `newuoa` is + superior to `uobyqa`); + +- `bobyqa` for *Bounded Optimization BY Quadratic Approximations* is for + simple bound constrained problems; + +- `cobyla` for *Constrained Optimization BY Linear Approximations* is for + general constrained problems; + +- `lincoa` is also for general constrained problems but, compared to `cobyla`, + linear equality and inequality constraints can be explicitly specified for + efficiency. + +These methods are called as follows: + +``` julia +using PRIMA +x, fx, nf, rc = uobyqa(f, x0; kwds...) +x, fx, nf, rc = newuoa(f, x0; kwds...) +x, fx, nf, rc = bobyqa(f, x0; kwds...) +x, fx, nf, rc, cstrv = cobyla(f, x0; kwds...) +x, fx, nf, rc, cstrv = lincoa(f, x0; kwds...) +``` + +where `$optim` is one of `uobyqa`, `newuoa`, `bobyqa`, `cobyla`, or `lincoa`, +`f` is the objective function, and `x0` is the initial solution. Constraints +and options may be specified by keywords `kwds...` (see below). The result is +the 4-tuple `(x, fx, nf,rc)` or the 5-tuple `(x, fx, nf, rc, cstrv)` where `x` +is the (approximate) solution found by the algorithm, `fx` is the value of +`f(x)`, `nf` is the number of calls to `f`, `rc` is a status code (an +enumeration value of type `PRIMA.Status`), and `cstrv` is the amount of +constraint violation. + +For example, `rc` can be: + +- `PRIMA.SMALL_TR_RADIUS` if the radius of the trust region becomes smaller or + equal the value of keyword `rhobeg`, in other words, the algorithm has + converged in terms of variable precision; + +- `PRIMA.FTARGET_ACHIEVED` if the objective function is smaller of equal the + value of keyword `ftarget`, in other words, the algorithm has converged in + terms of function value. + +There are other possibilities which all indicate an error. Calling: + +``` julia +PRIMA.reason(rc) +``` + +yields a textual explanation about the reason that leads the algorithm to stop. + +For all algorithms, except `cobyla`, the user-defined function takes a single +argument, the variables of the problem, and returns the value of the objective +function. It has the following signature: + +``` julia +function objfun(x::Vector{Cdouble}) + return f(x) +end +``` + +The `cobyla` algorithm aims at solving the problem: + +``` julia +min f(x) s.t. c(x) ≤ 0 +``` + +where `f(x)` is the objective function while `c(x)` implements `m` inequality +constraints. For this algorithm, the user-defined function takes two arguments, +the `x` variables `x` and the values taken by the constraints `cx`, it shall +overwrite the array of constraints with `cx = c(x)` and return the value of the +objective function. It has the following signature: + +``` julia +function objfuncon(x::Vector{Cdouble}, cx::Vector{Cdouble}) + copyto!(cx, c(x)) # set constraints + return f(x) # return value of objective function +end +``` + +The keywords allowed by the different algorithms are summarized by the +following table. + +| Keyword | Description | Algorithms | +|:-------------|:---------------------------------------|:-----------------------------| +| `rhobeg` | Initial trust region radius | all | +| `rhoend` | Final trust region radius | all | +| `ftarget` | Target objective function value | all | +| `maxfun` | Maximum number of function evaluations | all | +| `iprint` | Verbosity level | all | +| `npt` | Number of points in local model | `bobyqa`, `lincoa`, `newuoa` | +| `xl` | Lower bound | `bobyqa`, `cobyla`, `lincoa` | +| `xu` | Upper bound | `bobyqa`, `cobyla`, `lincoa` | +| `nlconstr` | Non-linear constraints | `cobyla` | +| `eqconstr` | Linear equality constraints | `cobyla`, `lincoa` | +| `ineqconstr` | Linear inequality constraints | `cobyla`, `lincoa` | + +Assuming `n = length(x)` is the number of variables, then: + +- `rhobeg` (default value `1.0`) is the initial radius of the trust region. + +- `rhoend` (default value `rhobeg*1e-4`) is the final radius of the trust + region used to decide whether the algorithm has converged in the variables. + +- `ftarget` (default value `-Inf`) is another convergence setting. The + algorithm is stopped as soon as `f(x) ≤ ftarget` and the status + `PRIMA.FTARGET_ACHIEVED` is returned. + +- `iprint` (default value `PRIMA.MSG_NONE`) sets the level of verbosity of the + algorithm. Possible values are `PRIMA.MSG_EXIT`, `PRIMA.MSG_RHO`, or + `PRIMA.MSG_FEVL`. + +- `maxfun` (default `100n`) is the maximum number of function evaluations + allowed for the algorithm. If the number of calls to `f(x)` exceeds this + value, the algorithm is stopped and the status `PRIMA.MAXFUN_REACHED` is + returned. + +- `npt` (default value `2n + 1`) is the number of points used to approximate + the local behavior of the objective function and such that `n + 1 ≤ npt ≤ + (n + 1)*(n + 2)/2`. The default value corresponds to the one recommended by + M.J.D. Powell. + +- `xl` and `xu` (default `fill(+Inf, n)` and `fill(-Inf, n)`) are the + elementwise lower and upper bounds for the variables. Feasible variables are + such that `xl ≤ x ≤ xu` (elementwise). + +- `nlconstr` (default `nothing`) may be specified as a vector of `m` double + precision floating-point values which are passed by `cobyla` to the + user-defined function to store `c(x)` the non-linear constraints in `x`. + +- `eqconstr` (default `nothing`) may be specified as a tuple `(A,b)` + to represent linear equality constraints. Feasible variables are + such that `A'⋅x = b` holds elementwise. + +- `neqconstr` (default `nothing`) may be specified as a tuple `(A,b)` to + represent linear inequality constraints. Feasible variables are such that + `A'⋅x ≤ b` holds elementwise. diff --git a/src/PRIMA.jl b/src/PRIMA.jl index 3d53ed5..1482a76 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -1,5 +1,432 @@ module PRIMA -# Write your package code here. +export bobyqa, cobyla, lincoa, newuoa, uobyqa + +using TypeUtils +using PRIMA_jll +const libprimac = PRIMA_jll.libprimac + +#------------------------------------------------------------------------------ +# PUBLIC INTERFACE + +# Verbosity level +@enum Message::Cint begin + MSG_NONE = 0 # No messages + MSG_EXIT = 1 # Exit reasons + MSG_RHO = 2 # Rho changes + MSG_FEVL = 3 # The object/constraint functions get evaluated +end + +# Possible return values +@enum Status::Cint begin + SMALL_TR_RADIUS = 0 + FTARGET_ACHIEVED = 1 + TRSUBP_FAILED = 2 + MAXFUN_REACHED = 3 + MAXTR_REACHED = 20 + NAN_INF_X = -1 + NAN_INF_F = -2 + NAN_INF_MODEL = -3 + NO_SPACE_BETWEEN_BOUNDS = 6 + DAMAGING_ROUNDING = 7 + ZERO_LINEAR_CONSTRAINT = 8 + INVALID_INPUT = 100 + ASSERTION_FAILS = 101 + VALIDATION_FAILS = 102 + MEMORY_ALLOCATION_FAILS = 103 +end + +""" + PRIMA.reason(rc) -> str + +yields a textual message explaining `rc`, the code returned by one of the PRIMA +optimizers. + +""" +reason(status::Union{Integer,Status}) = + unsafe_string(@ccall libprimac.prima_get_rc_string(Integer(status)::Cint)::Cstring) + +# The high level wrappers. +for func in (:bobyqa, :newuoa, :uobyqa, :lincoa, :cobyla) + func! = Symbol(func, "!") + @eval begin + @doc """ + $($func)(f, x0; kwds...) -> (x, fx, nf, rc) + + attempts to minimize objective function `f` by the + [`PRIMA.$($func)!`](@ref) method starting with variables `x0` (which + are left unchanged on output). The result is the 4-tuple `(x, fx, nf, + rc)` with `x` the (approximate) solution found, `fx` the value of + `f(x)`, `nf` the number of calls to `f`, and `rc` a status code (see + [`PRIMA.Status`](@ref) and [`PRIMA.reason`](@ref)). + + See the documention of [`PRIMA.$($func)!`](@ref) for a description of + the method and of the allowed keywords `kwds...`. + + """ + function $func(f, x0::AbstractVector{<:Real}; kwds...) + x = copyto!(Vector{Cdouble}(undef, length(x0)), x0) + return x, $func!(f, x; kwds...)... + end + end +end + +""" + PRIMA.bobyqa!(f, x; kwds...) -> (fx, nf, rc) + +attempts to minimize objective function `f` by M.J.D. Powell's BOBYQA method +for "Bounded Optimization BY Quadratic Approximations" and without derivatives. +On entry, argument `x` is a dense vector of double precision value with the +initial variables; on return, `x` is overwritten by an approximate solution and +the 3-tuple `(fx, nf, rc)` is returned with, `fx` the value of `f(x)`, `nf` the +number of calls to `f`, and `rc` a status code (see [`PRIMA.Status`](@ref) and +[`PRIMA.reason`](@ref)). + +Call [`PRIMA.bobyqa`](@ref) for a version that preserves the initial variables. + +""" +function bobyqa!(f, x::DenseVector{Cdouble}; + xl::DenseVector{Cdouble} = fill(typemin(Cdouble), length(x)), + xu::DenseVector{Cdouble} = fill(typemax(Cdouble), length(x)), + rhobeg::Real = 1.0, + rhoend::Real = rhobeg*1e-4, + iprint::Union{Integer,Message} = MSG_NONE, + ftarget::Real = -Inf, + maxfun::Integer = 100*length(x), + npt::Integer = 2*length(x) + 1) + n = length(x) # number of variables + @assert length(xl) == n + @assert length(xu) == n + @assert isfinite(rhobeg) && rhobeg > zero(rhobeg) + @assert isfinite(rhoend) && rhoend ≥ zero(rhoend) && rhoend ≤ rhobeg + @assert n + 1 ≤ npt ≤ (n + 1)*(n + 2)/2 + + fx = Ref{Cdouble}(NaN) # to store f(x) on return + nf = Ref{Cint}(0) # to store number of function calls + fw = ObjFun(f, n) # wrapper to objective function + fp = _push_wrapper(fw) # pointer to C-callable function + try + # Call low-level optimizer. + rc = @ccall libprimac.prima_bobyqa( + fp ::Ptr{Cvoid}, # C-type: prima_obj + n ::Cint, + x ::Ptr{Cdouble}, # n elements + fx ::Ref{Cdouble}, + xl ::Ptr{Cdouble}, # const, n elements + xu ::Ptr{Cdouble}, # const, n elements + nf ::Ref{Cint}, + rhobeg ::Cdouble, + rhoend ::Cdouble, + ftarget ::Cdouble, + maxfun ::Cint, # maxfun + npt ::Cint, + Integer(iprint) ::Cint)::Cint + return (fx[], Int(nf[]), rc) + finally + _pop_wrapper(fw) + end +end + +""" + PRIMA.newuoa!(f, x; kwds...) -> (fx, nf, rc) + +attempts to minimize objective function `f` by M.J.D. Powell's NEWUOA method +for unconstrained optimization without derivatives. On entry, argument `x` is +a dense vector of double precision value with the initial variables; on return, +`x` is overwritten by an approximate solution and the 3-tuple `(fx, nf, rc)` is +returned with, `fx` the value of `f(x)`, `nf` the number of calls to `f`, and +`rc` a status code (see [`PRIMA.Status`](@ref) and [`PRIMA.reason`](@ref)). + +Call [`PRIMA.newuoa`](@ref) for a version that preserves the +initial variables. + +""" +function newuoa!(f, x::DenseVector{Cdouble}; + rhobeg::Real = 1.0, + rhoend::Real = rhobeg*1e-4, + ftarget::Real = -Inf, + maxfun::Integer = 100*length(x), + npt::Integer = 2*length(x) + 1, + iprint::Union{Integer,Message} = MSG_NONE) + n = length(x) # number of variables + @assert isfinite(rhobeg) && rhobeg > zero(rhobeg) + @assert isfinite(rhoend) && rhoend ≥ zero(rhoend) && rhoend ≤ rhobeg + @assert n + 1 ≤ npt ≤ (n + 1)*(n + 2)/2 + + fx = Ref{Cdouble}(NaN) # to store f(x) on return + nf = Ref{Cint}(0) # to store number of function calls + fw = ObjFun(f, n) # wrapper to objective function + fp = _push_wrapper(fw) # pointer to C-callable function + try + # Call low-level optimizer. + rc = @ccall libprimac.prima_newuoa( + fp ::Ptr{Cvoid}, # C-type: prima_obj + n ::Cint, + x ::Ptr{Cdouble}, # n elements + fx ::Ref{Cdouble}, + nf ::Ref{Cint}, + rhobeg ::Cdouble, + rhoend ::Cdouble, + ftarget ::Cdouble, + maxfun ::Cint, + npt ::Cint, + Integer(iprint) ::Cint)::Cint + return (fx[], Int(nf[]), rc) + finally + _pop_wrapper(fw) + end +end + +""" + PRIMA.uobyqa!(f, x; kwds...) -> (fx, nf, rc) + +attempts to minimize objective function `f` by M.J.D. Powell's UOBYQA method +for "Unconstrained Optimization BY Quadratic Approximation" and without +derivatives. On entry, argument `x` is a dense vector of double precision value +with the initial variables; on return, `x` is overwritten by an approximate +solution and the 3-tuple `(fx, nf, rc)` is returned with, `fx` the value of +`f(x)`, `nf` the number of calls to `f`, and `rc` a status code (see +[`PRIMA.Status`](@ref) and [`PRIMA.reason`](@ref)). + +Call [`PRIMA.uobyqa`](@ref) for a version that preserves the initial variables. + +""" +function uobyqa!(f, x::DenseVector{Cdouble}; + rhobeg::Real = 1.0, + rhoend::Real = rhobeg*1e-4, + ftarget::Real = -Inf, + maxfun::Integer = 100*length(x), + iprint::Union{Integer,Message} = MSG_NONE) + n = length(x) # number of variables + @assert isfinite(rhobeg) && rhobeg > zero(rhobeg) + @assert isfinite(rhoend) && rhoend ≥ zero(rhoend) && rhoend ≤ rhobeg + + fx = Ref{Cdouble}(NaN) # to store f(x) on return + nf = Ref{Cint}(0) # to store number of function calls + fw = ObjFun(f, n) # wrapper to objective function + fp = _push_wrapper(fw) # pointer to C-callable function + try + # Call low-level optimizer. + rc = @ccall libprimac.prima_uobyqa( + fp ::Ptr{Cvoid}, # C-type: prima_obj + n ::Cint, + x ::Ptr{Cdouble}, # n elements + fx ::Ref{Cdouble}, + nf ::Ref{Cint}, + rhobeg ::Cdouble, + rhoend ::Cdouble, + ftarget ::Cdouble, + maxfun ::Cint, + Integer(iprint) ::Cint)::Cint + return (fx[], Int(nf[]), rc) + finally + _pop_wrapper(fw) + end +end + +const LinearConstraint = Tuple{DenseMatrix{Cdouble},DenseVector{Cdouble}} + +function cobyla!(f, x::DenseVector{Cdouble}; + nlconstr::Union{DenseVector{Cdouble},Nothing} = nothing, + ineqconstr::Union{LinearConstraint,Nothing} = nothing, + eqconstr::Union{LinearConstraint,Nothing} = nothing, + xl::DenseVector{Cdouble} = fill(typemin(Cdouble), length(x)), + xu::DenseVector{Cdouble} = fill(typemax(Cdouble), length(x)), + rhobeg::Real = 1.0, + rhoend::Real = rhobeg*1e-4, + ftarget::Real = -Inf, + maxfun::Integer = 100*length(x), + iprint::Union{Integer,Message} = MSG_NONE) + n = length(x) # number of variables + @assert length(xl) == n + @assert length(xu) == n + @assert isfinite(rhobeg) && rhobeg > zero(rhobeg) + @assert isfinite(rhoend) && rhoend ≥ zero(rhoend) && rhoend ≤ rhobeg + + m_nlcon, nlcon_ptr = _get_nonlinear_constraints(nlconstr) + m_eq, A_eq, b_eq = _get_linear_constraints(eqconstr, n) + m_ineq, A_ineq, b_ineq = _get_linear_constraints(ineqconstr, n) + + cstrv = Ref{Cdouble}(NaN) # to store constraint violation + fx = Ref{Cdouble}(NaN) # to store f(x) on return + nf = Ref{Cint}(0) # to store number of function calls + fw = ObjFunCon(f, n, m_nlcon) # wrapper to objective function + fp = _push_wrapper(fw) # pointer to C-callable function + + try + # Call low-level optimizer. + rc = GC.@preserve nlconstr ineqconstr eqconstr @ccall libprimac.prima_cobyla( + m_nlcon ::Cint, + fp ::Ptr{Cvoid}, # C-type: prima_objcon + n ::Cint, + x ::Ptr{Cdouble}, # n elements + fx ::Ref{Cdouble}, + cstrv ::Ref{Cdouble}, + nlcon_ptr ::Ptr{Cdouble}, # m_nlcon elements + m_ineq ::Cint, + A_ineq ::Ptr{Cdouble}, # const, m_ineq*n elements + b_ineq ::Ptr{Cdouble}, # const, m_ineq elements + m_eq ::Cint, + A_eq ::Ptr{Cdouble}, # const, m_eq*n elements + b_eq ::Ptr{Cdouble}, # const, m_eq elements + xl ::Ptr{Cdouble}, # const, n elements + xu ::Ptr{Cdouble}, # const, n elements + nf ::Ref{Cint}, + rhobeg ::Cdouble, + rhoend ::Cdouble, + ftarget ::Cdouble, + maxfun ::Cint, + Integer(iprint) ::Cint)::Cint + return (fx[], Int(nf[]), rc, cstrv[]) + finally + _pop_wrapper(fw) + end +end + +function lincoa!(f, x::DenseVector{Cdouble}; + ineqconstr::Union{LinearConstraint,Nothing} = nothing, + eqconstr::Union{LinearConstraint,Nothing} = nothing, + xl::DenseVector{Cdouble} = fill(typemin(Cdouble), length(x)), + xu::DenseVector{Cdouble} = fill(typemax(Cdouble), length(x)), + rhobeg::Real = 1.0, + rhoend::Real = rhobeg*1e-4, + ftarget::Real = -Inf, + maxfun::Integer = 100*length(x), + npt::Integer = 2*length(x) + 1, + iprint::Union{Integer,Message} = MSG_NONE) + n = length(x) # number of variables + @assert length(xl) == n + @assert length(xu) == n + @assert isfinite(rhobeg) && rhobeg > zero(rhobeg) + @assert isfinite(rhoend) && rhoend ≥ zero(rhoend) && rhoend ≤ rhobeg + @assert n + 1 ≤ npt ≤ (n + 1)*(n + 2)/2 + + m_eq, A_eq, b_eq = _get_linear_constraints(eqconstr, n) + m_ineq, A_ineq, b_ineq = _get_linear_constraints(ineqconstr, n) + + cstrv = Ref{Cdouble}(NaN) # to store constraint violation + fx = Ref{Cdouble}(NaN) # to store f(x) on return + nf = Ref{Cint}(0) # to store number of function calls + fw = ObjFun(f, n) # wrapper to objective function + fp = _push_wrapper(fw) # pointer to C-callable function + try + # Call low-level optimizer. + rc = GC.@preserve ineqconstr eqconstr @ccall libprimac.prima_lincoa( + fp ::Ptr{Cvoid}, # C-type: prima_objcon + n ::Cint, + x ::Ptr{Cdouble}, # n elements + fx ::Ref{Cdouble}, + cstrv ::Ref{Cdouble}, + m_ineq ::Cint, + A_ineq ::Ptr{Cdouble}, # const, m_ineq*n elements + b_ineq ::Ptr{Cdouble}, # const, m_ineq elements + m_eq ::Cint, + A_eq ::Ptr{Cdouble}, # const, m_eq*n elements + b_eq ::Ptr{Cdouble}, # const, m_eq elements + xl ::Ptr{Cdouble}, # const, n elements + xu ::Ptr{Cdouble}, # const, n elements + nf ::Ref{Cint}, + rhobeg ::Cdouble, + rhoend ::Cdouble, + ftarget ::Cdouble, + maxfun ::Cint, + npt ::Cint, + Integer(iprint) ::Cint)::Cint + return (fx[], Int(nf[]), rc, cstrv[]) + finally + _pop_wrapper(fw) + end +end + +#------------------------------------------------------------------------------ +# PRIVATE METHODS AND TYPES + +abstract type AbstractObjFun end + +# Small structure to wrap simple objective functions for BOBYQA, NEWUOA, +# UOBYQA, and LINCOA at low level. +struct ObjFun <: AbstractObjFun + f::Any # user-defined objective function + nx::Int # number of variables +end + +# Small structure to wrap objective functions with constraints for COBYLA at +# low level. +struct ObjFunCon <: AbstractObjFun + f::Any # user-defined onjective function + nx::Int # number of variables + nc::Int # number of non-linear constraints +end + +# Global variable storing the per-thread stacks of objective functions indexed +# by thread identifier and then by execution order. On start of an +# optimization, an object linked to the user-defined objective function is +# pushed. This object is popped out of the stack on return of the optimization +# call, whatever happens. It is therefore to wrap th call to the optimization +# method in a `try-finally` clause. +const _objfun_stack = Vector{Vector{ObjFun}}(undef, 0) +const _objfuncon_stack = Vector{Vector{ObjFunCon}}(undef, 0) + +# Private function `_get_stack` yields the stack for the caller thread and for +# a given type of objectve function. +_get_stack(fw::ObjFun) = _get_stack(_objfun_stack) +_get_stack(fw::ObjFunCon) = _get_stack(_objfuncon_stack) +function _get_stack(stack::Vector{Vector{T}}) where {T} + i = Threads.threadid() + while length(stack) < i + push!(stack, Vector{T}(undef, 0)) + end + return stack[i] +end + +# The following wrappers are intended to be C-callable functions. Their +# respective signatures must follow the prototypes `prima_obj` and +# `prima_objcon` in `prima.h` header. +function _objfun_wrapper(x_ptr::Ptr{Cdouble}, # (input) variables + f_ptr::Ptr{Cdouble}) # (output) function value + fw = last(_objfun_stack[Threads.threadid()]) # retrieve the objective function + x = unsafe_wrap(Array, x_ptr, fw.nx) + unsafe_store!(f_ptr, as(Cdouble, fw.f(x))) + return nothing +end +function _objfuncon_wrapper(x_ptr::Ptr{Cdouble}, # (input) variables + f_ptr::Ptr{Cdouble}, # (output) function value + c_ptr::Ptr{Cdouble}) # (output) constraints + fw = last(_objfuncon_stack[Threads.threadid()]) # retrieve the objective function + x = unsafe_wrap(Array, x_ptr, fw.nx) + c = unsafe_wrap(Array, c_ptr, fw.nc) + unsafe_store!(f_ptr, as(Cdouble, fw.f(x, c))) + return nothing +end + +# Private functions `_push_wrapper` and `_pop_wrapper` are to be used in a +# `try-finally` clause as explained above. +_c_wrapper(::ObjFun) = @cfunction(_objfun_wrapper, Cvoid, (Ptr{Cdouble}, Ptr{Cdouble},)) +_c_wrapper(::ObjFunCon) = @cfunction(_objfuncon_wrapper, Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},)) +function _push_wrapper(fw::AbstractObjFun) + push!(_get_stack(fw), fw) + return _c_wrapper(fw) +end +function _pop_wrapper(fw::AbstractObjFun) + stack = _get_stack(fw) + last(stack) === fw || error + ("objective function is not the last one in the caller tread stask") + resize!(stack, length(stack) - 1) + return nothing +end + +# FIXME: Matrix A of linear constraints is in row-major order (this is imposed +# by the C interface which transposes the matrix A). +_get_linear_constraints(::Nothing, n::Integer) = 0, Ptr{Cdouble}(0), Ptr{Cdouble}(0) +function _get_linear_constraints(Ab::LinearConstraint, n::Integer) + A, b = Ab + m = length(b) # number of constraints + @assert size(A) == (n, m) + return m, pointer(A), pointer(b) +end + +_get_nonlinear_constraints(::Nothing, n::Integer) = 0, Ptr{Cdouble}(0) +_get_nonlinear_constraints(c::DenseVector{Cdouble}) = + length(c), pointer(c) end diff --git a/test/runtests.jl b/test/runtests.jl index ab6b2c6..f848bed 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,138 @@ using PRIMA -using Test +using Test, TypeUtils @testset "PRIMA.jl" begin - # Write your tests here. + @testset "Utils " begin + for status in instances(PRIMA.Status) + @test PRIMA.reason(status) isa String + end + end + + # User defined objective function: + function f(x::AbstractVector{T}) where {T<:AbstractFloat} + @assert length(x) == 2 + x1, x2 = x + return as(T, 5*(x1 - 3)*(x1 - 3) + 7*(x2 - 2)*(x2 - 2) + as(T, 1//10)*(x1 + x2) - 10) + end + + # Objective function for COBYLA (same name but different signature to implement non-linear constraint). + function f(x::AbstractVector{T}, c::AbstractVector{T}) where {T<:AbstractFloat} + nlconstr!(x, c) + return f(x) + end + + # Non-linear constraint. + function nlconstr!(x::AbstractVector{T}, c::AbstractVector{T}) where {T<:AbstractFloat} + @assert length(x) == 2 + @assert length(c) == 1 + x1, x2 = x + c[firstindex(c)] = x1*x1 + x2*x2 - 13 # ‖x‖² ≤ 13 + return nothing + end + + # Array to store non-linear constraints. + nlconstr = Array{Cdouble}(undef, 1) + + # Initial solution. + x0 = [0.0, 0.0] + n = length(x0) + x0_sav = copy(x0) # for testing that x0 does not get overwritten + + # Inequality constraints: x₁ ≤ 4, x₂ ≤ 3, x₁ + x₂ ≤ 10 + A_ineq = [1.0 0.0 1.0; + 0.0 1.0 1.0] + b_ineq = [4.0, 3.0, 10.0] + + # Bounds. + xl = [-6.0, -6.0] + xu = [ 6.0, 6.0] + + function check_bounds(xl, x, xu) + flag = true + for (lᵢ, xᵢ, uᵢ) in zip(xl, x, xu) + flag &= lᵢ ≤ xᵢ ≤ uᵢ + end + return flag + end + + @testset "NEWUOA" begin + x, fx, nf, rc = @inferred PRIMA.newuoa(f, x0; + rhobeg = 1.0, rhoend = 1e-3, + ftarget = -Inf, + maxfun = 200n, + npt = 2n + 1, + iprint = PRIMA.MSG_EXIT) + msg = PRIMA.reason(rc) + println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + @test abs(x[1] - 3) ≤ 2e-2 && abs(x[2] - 2) ≤ 2e-2 + @test f(x) ≈ fx + @test x0 == x0_sav + end + + @testset "UOBYQA" begin + x, fx, nf, rc = @inferred PRIMA.uobyqa(f, x0; + rhobeg = 1.0, rhoend = 1e-3, + ftarget = -Inf, + maxfun = 200n, + iprint = PRIMA.MSG_EXIT) + msg = PRIMA.reason(rc) + println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + @test abs(x[1] - 3) ≤ 2e-2 && abs(x[2] - 2) ≤ 2e-2 + @test f(x) ≈ fx + @test x0 == x0_sav + end + + @testset "BOBYQA" begin + x, fx, nf, rc = @inferred PRIMA.bobyqa(f, x0; + xl, xu, + rhobeg = 1.0, rhoend = 1e-3, + ftarget = -Inf, + maxfun = 200n, + npt = 2n + 1, + iprint = PRIMA.MSG_EXIT) + msg = PRIMA.reason(rc) + println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + @test abs(x[1] - 3) ≤ 2e-2 && abs(x[2] - 2) ≤ 2e-2 + @test f(x) ≈ fx + @test x0 == x0_sav + @test check_bounds(xl, x, xu) + end + + @testset "COBYLA" begin + x, fx, nf, rc, cstrv = @inferred PRIMA.cobyla(f, x0; + nlconstr, + ineqconstr = (A_ineq, b_ineq), + xl, xu, + rhobeg = 1.0, rhoend = 1e-3, + ftarget = -Inf, + maxfun = 200*n, + iprint = PRIMA.MSG_EXIT) + msg = PRIMA.reason(rc) + println("x = $x, f(x) = $fx, cstrv = $cstrv, nlconstr = $nlconstr, rc = $rc, msg = '$msg', evals = $nf") + @test abs(x[1] - 3) ≤ 2e-2 && abs(x[2] - 2) ≤ 2e-2 + @test f(x) ≈ fx + @test x0 == x0_sav + @test check_bounds(xl, x, xu) + end + + + @testset "LINCOA" begin + x, fx, nf, rc, cstrv = @inferred PRIMA.lincoa(f, x0; + ineqconstr = (A_ineq, b_ineq), + xl, xu, + rhobeg = 1.0, rhoend = 1e-3, + ftarget = -Inf, + maxfun = 200*n, + npt = 2n + 1, + iprint = PRIMA.MSG_EXIT) + msg = PRIMA.reason(rc) + println("x = $x, f(x) = $fx, cstrv = $cstrv, nlconstr = $nlconstr, rc = $rc, msg = '$msg', evals = $nf") + @test abs(x[1] - 3) ≤ 2e-2 && abs(x[2] - 2) ≤ 2e-2 + @test f(x) ≈ fx + @test x0 == x0_sav + @test check_bounds(xl, x, xu) + end + end + +nothing