diff --git a/.github/actions/spelling/allow.txt b/.github/actions/spelling/allow.txt index 7ab0963..5165ba3 100644 --- a/.github/actions/spelling/allow.txt +++ b/.github/actions/spelling/allow.txt @@ -1103,6 +1103,7 @@ isrvub isscalar isspace isstruct +issuccess issymmetric istp istr diff --git a/NEWS.md b/NEWS.md index f8b7c17..e6aeca9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,41 @@ # User visible changes in `PRIMA.jl` +## Version 0.2.0 + +- General method `prima(f, x0; kwds...)` which solves the problem with the most + suitable Powell's algorithm (among COBYLA, LINCOA, BOBYQA, or NEWUOA) + depending on the constraints imposed via the keywords `kwds...`. + +- Objective function and non-linear constraints for `cobyla`: + - Non-linear equality constraints can be specified by keyword `nonlinear_eq`. + - The objective function is called as `f(x)` like in other algorithms. + - The functions implementing the non-linear constraints are passed by + keywords `nonlinear_eq` and `nonlinear_ineq`. + +- Algorithms have a more similar interface: + - All algorithms have the same positional input, only the available keywords + may change. + - All algorithms have the same convention for the objective function. + Non-linear constraints, if any, are provided by other user-defined + functions. + - All algorithms yield a 2-tuple `(x,info)` with `x` an approximate solution + (provided algorithm was successful) and `info` a structured object + collecting other outputs from the algorithm. The following properties are + available for all algorithms: `info.fx` is the objective function value + `f(x)`, `info.nf` is the number of evaluations of the objective function, + and `info.status` is the termination status of the algorithm. + - `issuccess(info)` yields whether algorithm was successful. + - `PRIMA.reason(info)` yields a textual description of the termination status + of the algorithm. + +- Scaling of the variables: the variables `x ∈ ℝⁿ` may be scaled by the scaling + factors provided by the caller via keywords `scale` to re-express the problem + in the scaled variables `u ∈ ℝⁿ` such that `u[i] = x[i]/scale[i]`. Note that + the objective function, the constraints (linear and non-linear), and the + bounds remain specified in the variables. Scaling the variables is useful to + improve the conditioning of the problem, to make the scaled variables `u` + having approximately the same magnitude, and to adapt to heterogeneous + variables or with different units. ## Version 0.1.1 diff --git a/Project.toml b/Project.toml index a3dbf56..98547d5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "PRIMA" uuid = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed" authors = ["Éric Thiébaut and contributors"] -version = "0.1.1" +version = "0.2.0" [deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PRIMA_jll = "eead6e0c-2d5b-5641-a95c-b722de96d551" TypeUtils = "c3b1956e-8857-4d84-9b79-890df85b1e67" diff --git a/README.md b/README.md index daed0bb..497da4c 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,19 @@ Powell](https://en.wikipedia.org/wiki/Michael_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: +Depending on the problem to solve, other Julia package(s) with similar +objectives may be of interest: + +- [BlackBoxOptim](https://github.com/robertfeldt/BlackBoxOptim.jl) is a global + optimization package for single- and multi-variate problems that does not + require the objective function to be differentiable. + +- [NOMAD](https://github.com/bbopt/NOMAD.jl) is an interface to the *Mesh + Adaptive Direct Search algorithm* (MADS) designed for difficult black-box + optimization problems. + +Formally, the algorithms provided by `PRIMA` are designed to solve problems of +the form: ``` julia min f(x) subject to x ∈ Ω ⊆ ℝⁿ @@ -20,12 +32,14 @@ variables, and `n ≥ 1` is the number of variables. The most general feasible set is: ``` julia -Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, Aᵢ⋅x ≤ bᵢ, and c(x) ≤ 0 } +Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, Aᵢ⋅x ≤ bᵢ, cₑ(x) = 0, and cᵢ(x) ≤ 0 } ``` where `xl ∈ ℝⁿ` and `xu ∈ ℝⁿ` are lower and upper bounds, `Aₑ` and `bₑ` implement linear equality constraints, `Aᵢ` and `bᵢ` implement linear -inequality constraints, and `c: ℝⁿ → ℝᵐ` implements `m` non-linear constraints. +inequality constraints, `cₑ: ℝⁿ → ℝᵐ` implements `m` non-linear equality +constraints, and `cᵢ: ℝⁿ → ℝʳ` implements `r` non-linear inequality +constraints. The five Powell's algorithms of the [PRIMA library](https://github.com/libprima/prima) are provided by the `PRIMA` @@ -67,66 +81,77 @@ constraints). | `lincoa` | quadratic | bounds, linear | | `cobyla` | affine | bounds, linear, non-linear | -These methods are called as follows: +To use the most suitable Powell's algorithm depending on the constraints, +simply do: ``` 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...) +x, info = prima(f, x0; kwds...) ``` -where `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: +To explicitly use one the specific Powell's algorithms, call one of: -- `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; +``` julia +x, info = uobyqa(f, x0; kwds...) +x, info = newuoa(f, x0; kwds...) +x, info = bobyqa(f, x0; kwds...) +x, info = cobyla(f, x0; kwds...) +x, info = lincoa(f, x0; kwds...) +``` -- `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. +In any case, `f` is the objective function and `x0` specifies the initial +values of the variables (and is left unchanged). Constraints and options may be +specified by keywords `kwds...` (see below). -There are other possibilities which all indicate an error. Calling: +The objective function is called as `f(x)` with `x` the variables, it must +implement the following signature: ``` julia -PRIMA.reason(rc) +f(x::Vector{Cdouble})::Real ``` -yields a textual explanation about the reason that leads the algorithm to stop. +All the algorithms return a 2-tuple `(x, info)` with `x` the variables and +`info` a structured object collecting all other information. If +`issuccess(info)` is true, then the algorithm was successful and `x` is an +approximate solution of the problem. + +The output `info` has the following properties: + +``` julia +info.fx # value of the objective function f(x) on return +info.nf # number of calls to the objective function +info.status # final status code +info.cstrv # amount of constraints violation, 0.0 if unconstrained +info.nl_eq # non-linear equality constraints, empty vector if none +info.nl_ineq # non-linear inequality constraints, empty vector if none +``` -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: +Calling one of: ``` julia -function objfun(x::Vector{Cdouble}) - return f(x) -end +issuccess(info) +issuccess(info.status) ``` -The `cobyla` algorithm can account for `m` non-linear constraints expressed by -`c(x) ≤ 0`. 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: +yield whether the algorithm has converged. If this is the case, `info.status` +can be one of: + +- `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 a failure. Calling one of: ``` julia -function objfuncon(x::Vector{Cdouble}, cx::Vector{Cdouble}) - copyto!(cx, c(x)) # set constraints - return f(x) # return value of objective function -end +PRIMA.reason(info) +PRIMA.reason(info.status) ``` +yield a textual explanation about the reason that leads the algorithm to stop. + The keywords allowed by the different algorithms are summarized by the following table. @@ -140,12 +165,30 @@ following table. | `npt` | Number of points in local model | `bobyqa`, `lincoa`, `newuoa` | | `xl` | Lower bound | `bobyqa`, `cobyla`, `lincoa` | | `xu` | Upper bound | `bobyqa`, `cobyla`, `lincoa` | +| `nonlinear_eq` | Non-linear equality constraints | `cobyla` | | `nonlinear_ineq` | Non-linear inequality constraints | `cobyla` | | `linear_eq` | Linear equality constraints | `cobyla`, `lincoa` | | `linear_ineq` | Linear inequality constraints | `cobyla`, `lincoa` | Assuming `n = length(x)` is the number of variables, then: +- `scale` (default value `nothing`) may be set with a vector of `n` positive + scaling factors. If specified, the problem is solved in the scaled variables + `u ∈ ℝⁿ` such that `u[i] = x[i]/scale[i]`. If unspecified, it is assumed that + `scale[i] = 1` for all variables. Note that the objective function, the + constraints (linear and non-linear), and the bounds remain specified in the + variables. Scaling the variables is useful to improve the conditioning of the + problem, to make the scaled variables `u` having approximately the same + magnitude, and to adapt to heterogeneous variables or with different units. + +- `rhobeg` (default value `1.0`) is the initial radius of the trust region. The + radius of the trust region is given by the Euclidean norm of the scaled + variables (see keyword `scale` above). + +- `rhoend` (default value `1e-$*rhobeg`) is the final radius of the trust + region used to decide whether the algorithm has converged in the scaled + variables. + - `rhobeg` (default value `1.0`) is the initial radius of the trust region. - `rhoend` (default value `1e-4*rhobeg`) is the final radius of the trust @@ -157,8 +200,9 @@ Assuming `n = length(x)` is the number of variables, then: `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`. + algorithm. Possible values are `PRIMA.MSG_EXIT`, `PRIMA.MSG_RHO`, or + `PRIMA.MSG_FEVL`. Note that the values that are printed by the software are + those of the scaled variables (see keyword `scale` above). - `maxfun` (default `100n`) is the maximum number of function evaluations allowed for the algorithm. If the number of calls to `f(x)` exceeds this @@ -171,22 +215,24 @@ Assuming `n = length(x)` is the number of variables, then: 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). + element-wise lower and upper bounds for the variables. Feasible variables are + such that `xl ≤ x ≤ xu`. + +- `nonlinear_eq` (default `nothing`) may be specified with a function, say + `c_eq`, implementing non-linear equality constraints defined by `c_eq(x) = + 0`. On return, the values of the non-linear equality constraints are given by + `info.nl_eq` to save calling `c_eq(x)`. -- `nonlinear_ineq` (default `nothing`) may be specified with the number `m` of - non-linear inequality constraints expressed `c(x) ≤ 0`. If the caller is - interested in the values of `c(x)` at the returned solution, the keyword may - be set with a vector of `m` double precision floating-point values to store - `c(x)`. This keyword only exists for `cobyla`. +- `nonlinear_ineq` (default `nothing`) may be specified with a function, say + `c_ineq`, implementing non-linear inequality constraints defined by + `c_ineq(x) ≤ 0`. On return, the values of the non-linear inequality + constraints are given by `info.nl_ineq` to save calling `c_ineq(x)`. - `linear_eq` (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. + impose the linear equality constraints `Aₑ⋅x = bₑ`. - `linear_ineq` (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. + impose the linear inequality constraints `Aᵢ⋅x ≤ bᵢ`. ## References diff --git a/gen/wrapper.jl b/gen/wrapper.jl index 626fff5..0a63a25 100644 --- a/gen/wrapper.jl +++ b/gen/wrapper.jl @@ -27,9 +27,11 @@ function main() r"\bprima_message\b" => "Message", r"\bprima_rc\b" => "Status", r"\bPRIMA_" => "", - # Make enum signed (see PRIMA library doc. about negative `iprint` - # values). - r"^(\s*@enum\s+\w+)\s*::\s*U(Int[0-9]*)\s+begin\s*$" => s"\1::\2 begin", + # Force enums to have signed value of type Cint (see PRIMA library doc. + # about negative `iprint` values). + r"^ *@enum +(\w+) *:: *U?Int\d+ +begin *$" => s"@enum \1::Cint begin", + # All algorithms return a Status. + r"\) *:: *Cint *$" => s")::Status", # Remove some useless code. r"^\s*const\s+PRIMAC_API\s*=.*$" => "", ] diff --git a/src/PRIMA.jl b/src/PRIMA.jl index aa1e5f5..feae3ad 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -4,54 +4,131 @@ using PRIMA_jll const libprimac = PRIMA_jll.libprimac include("wrappers.jl") -export bobyqa, cobyla, lincoa, newuoa, uobyqa +export bobyqa, cobyla, lincoa, newuoa, prima, uobyqa, issuccess using TypeUtils +using LinearAlgebra #------------------------------------------------------------------------------ # PUBLIC INTERFACE """ - PRIMA.reason(rc) -> str + PRIMA.Info -yields a textual message explaining `rc`, the code returned by one of the PRIMA -optimizers. +is the type of the structured object used to collect various information +provided by the optimization algorithms of the `PRIMA` package. + +An object, say `info`, of this type has the following properties: + + info.fx # value of the objective function f(x) on return + info.nf # number of calls to the objective function + info.status # final status code + info.cstrv # amount of constraints violation, 0.0 if unconstrained + info.nl_eq # non-linear equality constraints, empty vector if none + info.nl_ineq # non-linear inequality constraints, empty vector if none + +Call `issuccess(info)` or `issuccess(info.status)` to check whether algorithm +was successful. Call `PRIMA.reason(info)` or `PRIMA.reason(info.status)` to +retrieve a textual description of the algorithm termination. + +""" +struct Info + fx::Cdouble # f(x) on return + nf::Int # number of function calls + status::Status # returned status code + cstrv::Cdouble # amount of constraints violation + nl_eq::Vector{Cdouble} # non-linear equality constraints + nl_ineq::Vector{Cdouble} # non-linear inequality constraints + function Info(; # Mandatory keywords. + fx::Real, + nf::Integer, + status::Status, + # Optional keywords. + cstrv::Real = 0.0, + nl_eq::AbstractVector{<:Real} = Cdouble[], + nl_ineq::AbstractVector{<:Real} = Cdouble[]) + return new(fx, nf, status, cstrv, nl_eq, nl_ineq) + end +end + +Base.:(==)(a::Info, b::Info) = + ((a.fx == b.fx) & (a.nf == b.nf) & (a.status == b.status) & (a.cstrv == b.cstrv)) && + (a.nl_eq == b.nl_eq) && (a.nl_ineq == b.nl_ineq) + +LinearAlgebra.issuccess(info::Info) = issuccess(info.status) +LinearAlgebra.issuccess(status::Status) = + +""" + PRIMA.reason(info::PRIMA.Info) -> str + PRIMA.reason(status::PRIMA.Status) -> str + +yield a textual message explaining `info.status` or `status`, the status code +returned by one of the PRIMA optimizers. """ +reason(info::Info) = reason(info.status) reason(status::Union{Integer,Status}) = unsafe_string(prima_get_rc_string(status)) +""" + PRIMA.LinearConstraints + +is the type of `(A,b)`, the 2-tuple representing linear equality constraints +`A⋅x = b` or linear inequality constraints `A⋅x ≤ b` where `A` is a matrix, `x` +is the vector of variables, and `b` is a vector. + +""" +const LinearConstraints = Tuple{AbstractMatrix{<:Real},AbstractVector{<:Real}} + +# Default settings. +default_npt(x::AbstractVector{<:Real}) = 2*length(x) + 1 +default_maxfun(x::AbstractVector{<:Real}) = 100*length(x) +const default_scale = nothing +const default_rhobeg = 1.0 +const default_rhoend_relative = 1e-4 +default_rhoend(rhobeg::Real) = default_rhoend_relative*rhobeg + # The high level wrappers. First the methods, then their documentation. -for func in (:bobyqa, :newuoa, :uobyqa, :lincoa, :cobyla) +for func in (:bobyqa, :cobyla, :lincoa, :newuoa, :prima, :uobyqa) func! = Symbol(func, "!") @eval begin function $func(f, x0::AbstractVector{<:Real}; kwds...) x = copyto!(Vector{Cdouble}(undef, length(x0)), x0) - return x, $func!(f, x; kwds...)... + return x, $func!(f, x; kwds...) end end end -const _doc_2_inputs_4_outputs = """ +const _doc_common = """ The arguments are the objective function `f` and the initial variables `x0` -which are left unchanged on output. The returned value is the 4-tuple `(x, fx, -nf, rc)` with `x` the (approximate) solution, `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)). -""" +which are left unchanged on output. The objective function is called as `f(x)` +with `x` the variables, it must implement the following signature: -const _doc_2_inputs_5_outputs = """ -The arguments are the objective function `f` and the initial variables `x0` -which are left unchanged on output. The returned value is the 5-tuple `(x, fx, -nf, rc, cstrv)` with `x` the (approximate) solution, `fx` the value of `f(x)`, -`nf` the number of calls to `f`, `rc` a status code (see [`PRIMA.Status`](@ref) -and [`PRIMA.reason`](@ref)), and `cstrv` is the amount of constraint violation. -""" + f(x::Vector{Cdouble})::Real + +The returned result is a 2-tuple `(x, info)` with `x` the variables and `info` +a structured object collecting other information provided by the algorithm (see +[`PRIMA.Info`](@ref)). If `issuccess(info)` is true, then the algorithm was +successful and `x` is an approximate solution of the problem. + +Allowed keywords are (`n = length(x)` is the number of variables): + +- `scale` (default value `$default_scale`) may be set with a vector of `n` + positive scaling factors. If specified, the problem is solved in the scaled + variables `u ∈ ℝⁿ` such that `u[i] = x[i]/scale[i]`. If unspecified, it is + assumed that `scale[i] = 1` for all variables. Note that the objective + function, the constraints (linear and non-linear), and the bounds remain + specified in the variables. Scaling the variables is useful to improve the + conditioning of the problem, to make the scaled variables `u` having + approximately the same magnitude, and to adapt to heterogeneous variables or + with different units. -const _doc_common_keywords = """ -- `rhobeg` (default value `1.0`) is the initial radius of the trust region. +- `rhobeg` (default value `$default_rhobeg`) is the initial radius of the trust + region. The radius of the trust region is given by the Euclidean norm of the + scaled variables (see keyword `scale` above). -- `rhoend` (default value `1e-4*rhobeg`) is the final radius of the trust - region used to decide whether the algorithm has converged in the variables. +- `rhoend` (default value `$default_rhoend_relative*rhobeg`) is the final + radius of the trust region used to decide whether the algorithm has converged + in the scaled variables. - `ftarget` (default value `-Inf`) is another convergence setting. The algorithm is stopped as soon as `f(x) ≤ ftarget` and the status @@ -64,7 +141,8 @@ const _doc_common_keywords = """ - `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`. + `PRIMA.MSG_FEVL`. Note that the values that are printed by the software + are those of the scaled variables (see keyword `scale` above). """ const _doc_npt = """ @@ -81,29 +159,64 @@ const _doc_bound_constraints = """ """ const _doc_nonlinear_constraints = """ -- `nonlinear_ineq` (default `nothing`) may be specified with the number `m` of - non-linear inequality constraints expressed `c(x) ≤ 0`. If the caller is - interested in the values of `c(x)` at the returned solution, the keyword may - be set with a vector of `m` double precision floating-point values to store - `c(x)`. +- `nonlinear_eq` (default `nothing`) may be specified with a function, say + `c_eq`, implementing non-linear equality constraints defined by `c_eq(x) = + 0`. On return, the values of the non-linear equality constraints are given by + `info.nl_eq` to save calling `c_eq(x)`. + +- `nonlinear_ineq` (default `nothing`) may be specified with a function, say + `c_ineq`, implementing non-linear inequality constraints defined by + `c_ineq(x) ≤ 0`. On return, the values of the non-linear inequality + constraints are given by `info.nl_ineq` to save calling `c_ineq(x)`. """ const _doc_linear_constraints = """ -- `linear_eq` (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. +- `linear_eq` (default `nothing`) may be specified as a tuple `(Aₑ,bₑ)` to + impose the linear equality constraints `Aₑ⋅x = bₑ`. - `linear_ineq` (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. + impose the linear inequality constraints `Aᵢ⋅x ≤ bᵢ`. """ """ - uobyqa(f, x0; kwds...) -> (x, fx, nf, rc) + prima(f, x0; kwds...) -> x, info + +approximately solves the constrained optimization problem: + + min f(x) subject to x ∈ Ω ⊆ ℝⁿ + +with + + Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, Aᵢ⋅x ≤ bᵢ, cₑ(x) = 0, and cᵢ(x) ≤ 0 } + +by one of the M.J.D. Powell's algorithms COBYLA, LINCOA, BOBYQA, or NEWUOA +depending on the constraints set by `Ω`. These algorithms are based on a trust +region method where variables are updated according to a linear or a quadratic +local approximation of the objective function. No derivatives of the objective +function are needed. + +$(_doc_common) + +$(_doc_npt) + +$(_doc_bound_constraints) + +$(_doc_linear_constraints) + +$(_doc_nonlinear_constraints) + +See also [`PRIMA.cobyla`](@ref), [`PRIMA.lincoa`](@ref), +[`PRIMA.bobyqa`](@ref), [`PRIMA.newuoa`](@ref), or [`PRIMA.uobyqa`](@ref) to +use one of the specific Powell's algorithms. + +""" prima + +""" + uobyqa(f, x0; kwds...) -> x, info approximately solves the unconstrained optimization problem: - min f(x) s.t. x ∈ ℝⁿ + min f(x) subject to x ∈ ℝⁿ by M.J.D. Powell's UOBYQA (for \"Unconstrained Optimization BY Quadratic Approximations\") method. This algorithm is based on a trust region method @@ -111,52 +224,34 @@ where variables are updated according to a quadratic local approximation interpolating the objective function. No derivatives of the objective function are needed. -$(_doc_2_inputs_4_outputs) - -The objective function takes a single argument, the variables `x`, and returns -the function value, it shall implement the following signature: - - f(x::Vector{Cdouble})::Real - -Allowed keywords are (`n = length(x)` is the number of variables): - -$(_doc_common_keywords) +$(_doc_common) """ uobyqa """ - newuoa(f, x0; kwds...) -> (x, fx, nf, rc) + newuoa(f, x0; kwds...) -> x, info approximately solves the unconstrained optimization problem: - min f(x) s.t. x ∈ ℝⁿ + min f(x) subject to x ∈ ℝⁿ by M.J.D. Powell's NEWUOA method. This algorithm is based on a trust region method where variables are updated according to a quadratic local approximation interpolating the objective function at a number of `npt` points. No derivatives of the objective function are needed. -$(_doc_2_inputs_4_outputs) - -The objective function takes a single argument, the variables `x`, and returns -the function value, it shall implement the following signature: - - f(x::Vector{Cdouble})::Real - -Allowed keywords are (`n = length(x)` is the number of variables): - -$(_doc_common_keywords) +$(_doc_common) $(_doc_npt) """ newuoa """ - bobyqa(f, x0; kwds...) -> (x, fx, nf, rc) + bobyqa(f, x0; kwds...) -> x, info approximately solves the bound constrained optimization problem: - min f(x) s.t. x ∈ Ω ⊆ ℝⁿ + min f(x) subject to x ∈ Ω ⊆ ℝⁿ with @@ -168,16 +263,7 @@ where variables are updated according to a quadratic local approximation interpolating the objective function at a number of `npt` points. No derivatives of the objective function are needed. -$(_doc_2_inputs_4_outputs) - -The objective function takes a single argument, the variables `x`, and returns -the function value, it shall implement the following signature: - - f(x::Vector{Cdouble})::Real - -Allowed keywords are (`n = length(x)` is the number of variables): - -$(_doc_common_keywords) +$(_doc_common) $(_doc_npt) @@ -186,35 +272,22 @@ $(_doc_bound_constraints) """ bobyqa """ - cobyla(f, x0; kwds...) -> (x, fx, nf, rc, cstrv) + cobyla(f, x0; kwds...) -> x, info approximately solves the constrained optimization problem: - min f(x) s.t. x ∈ Ω ⊆ ℝⁿ + min f(x) subject to x ∈ Ω ⊆ ℝⁿ with - Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, Aᵢ⋅x ≤ bᵢ, and c(x) ≤ 0 } + Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, Aᵢ⋅x ≤ bᵢ, cₑ(x) = 0, and cᵢ(x) ≤ 0 } by M.J.D. Powell's COBYLA (for \"Constrained Optimization BY Linear Approximations\") method. This algorithm is based on a trust region method where variables are updated according to a linear local approximation of the objective function. No derivatives of the objective function are needed. -$(_doc_2_inputs_5_outputs) - -The objective function takes two arguments, the variables `x` and a vector `cx` -to store the non-linear constraints, and returns the function value, it shall -implement the following signature: - - f(x::Vector{Cdouble}, cx::Vector{Cdouble})::Real - -where the e,tries of `cx` are to be overwritten by the non-linear constraints -`c(x)`. - -Allowed keywords are (`n = length(x)` is the number of variables): - -$(_doc_common_keywords) +$(_doc_common) $(_doc_bound_constraints) @@ -225,11 +298,11 @@ $(_doc_nonlinear_constraints) """ cobyla """ - lincoa(f, x0; kwds...) -> (x, fx, nf, rc, cstrv) + lincoa(f, x0; kwds...) -> x, info approximately solves the constrained optimization problem: - min f(x) s.t. x ∈ Ω ⊆ ℝⁿ + min f(x) subject to x ∈ Ω ⊆ ℝⁿ with @@ -240,16 +313,7 @@ This algorithm is based on a trust region method where variables are updated according to a quadratic local approximation of the objective function. No derivatives of the objective function are needed. -$(_doc_2_inputs_5_outputs) - -The objective function takes a single argument, the variables `x`, and returns -the function value, it shall implement the following signature: - - f(x::Vector{Cdouble})::Real - -Allowed keywords are (`n = length(x)` is the number of variables): - -$(_doc_common_keywords) +$(_doc_common) $(_doc_npt) @@ -260,7 +324,48 @@ $(_doc_linear_constraints) """ lincoa """ - PRIMA.bobyqa!(f, x; kwds...) -> (fx, nf, rc) + PRIMA.prima!(f, x; kwds...) -> info::PRIMA.Info + +in-place version of [`PRIMA.prima`](@ref) which to see for details. 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. + +""" +function prima!(f, x::DenseVector{Cdouble}; + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, + rhoend::Real = default_rhoend(rhobeg), + ftarget::Real = -Inf, + maxfun::Integer = default_maxfun(x), + npt::Integer = default_npt(x), + iprint::Union{Integer,Message} = MSG_NONE, + xl::Union{AbstractVector{<:Real},Nothing} = nothing, + xu::Union{AbstractVector{<:Real},Nothing} = nothing, + linear_ineq::Union{LinearConstraints,Nothing} = nothing, + linear_eq::Union{LinearConstraints,Nothing} = nothing, + nonlinear_ineq = nothing, + nonlinear_eq = nothing) + if nonlinear_eq !== nothing || nonlinear_ineq !== nothing + # Only COBYLA can cope with non-linear constraints. + return cobyla!(f, x; scale, rhobeg, rhoend, iprint, ftarget, maxfun, + xl, xu, linear_ineq, linear_eq, + nonlinear_ineq, nonlinear_eq) + elseif linear_eq !== nothing || linear_ineq !== nothing + # LINCOA is the most efficient for linearly constrained problems. + return lincoa!(f, x; scale, rhobeg, rhoend, iprint, ftarget, maxfun, npt, + xl, xu, linear_ineq, linear_eq) + elseif xu !== nothing || xl !== nothing + # BOBYQA is designed for bound constrained problems. + return bobyqa!(f, x; scale, rhobeg, rhoend, iprint, ftarget, maxfun, npt, + xl, xu) + else + # NEWUOA is more efficient than UOBYQA for unconstrained problems. + return newuoa!(f, x; scale, rhobeg, rhoend, iprint, ftarget, maxfun, npt) + end +end + +""" + PRIMA.bobyqa!(f, x; kwds...) -> info::PRIMA.Info in-place version of [`PRIMA.bobyqa`](@ref) which to see for details. On entry, argument `x` is a dense vector of double precision value with the initial @@ -268,36 +373,44 @@ variables; on return, `x` is overwritten by an approximate solution. """ function bobyqa!(f, x::DenseVector{Cdouble}; + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, + rhoend::Real = default_rhoend(rhobeg), xl::Union{AbstractVector{<:Real},Nothing} = nothing, xu::Union{AbstractVector{<:Real},Nothing} = nothing, - rhobeg::Real = 1.0, - rhoend::Real = 1e-4*rhobeg, - iprint::Union{Integer,Message} = MSG_NONE, ftarget::Real = -Inf, - maxfun::Integer = 100*length(x), - npt::Integer = 2*length(x) + 1) + maxfun::Integer = default_maxfun(x), + npt::Integer = default_npt(x), + iprint::Union{Integer,Message} = MSG_NONE) # Check arguments and get constraints. n = length(x) # number of variables - _check_rho(rhobeg, rhoend) _check_npt(npt, n) - xl = _get_lower_bound(xl, n) - xu = _get_upper_bound(xu, n) + _check_rho(rhobeg, rhoend) + scl = _get_scaling(scale, n) + xl = _get_lower_bound(xl, n, scl) + xu = _get_upper_bound(xu, n, scl) + # References for output values. 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 + + # Create wrapper to objective function and push it on top of per-thread + # stack before calling optimizer. + fw = ObjFun(f, n, scl) # wrapper to objective function + fp = _push_objfun(bobyqa, fw) # pointer to C-callable function try - # Call low-level optimizer. - rc = prima_bobyqa(fp, n, x, fx, xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) - return (fx[], Int(nf[]), rc) + # Call low-level optimizer on the (scaled) variables. + isempty(scl) || _scale!(x, /, scl) + status = prima_bobyqa(fp, n, x, fx, xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) + isempty(scl) || _scale!(x, *, scl) + return Info(; fx = fx[], nf = nf[], status) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end """ - PRIMA.newuoa!(f, x; kwds...) -> (fx, nf, rc) + PRIMA.newuoa!(f, x; kwds...) -> info::PRIMA.Info in-place version of [`PRIMA.newuoa`](@ref) which to see for details. On entry, argument `x` is a dense vector of double precision value with the initial @@ -305,32 +418,40 @@ variables; on return, `x` is overwritten by an approximate solution. """ function newuoa!(f, x::DenseVector{Cdouble}; - rhobeg::Real = 1.0, - rhoend::Real = 1e-4*rhobeg, + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, + rhoend::Real = default_rhoend(rhobeg), ftarget::Real = -Inf, - maxfun::Integer = 100*length(x), - npt::Integer = 2*length(x) + 1, + maxfun::Integer = default_maxfun(x), + npt::Integer = default_npt(x), iprint::Union{Integer,Message} = MSG_NONE) # Check arguments. n = length(x) # number of variables - _check_rho(rhobeg, rhoend) _check_npt(npt, n) + _check_rho(rhobeg, rhoend) + scl = _get_scaling(scale, n) + # References for output values. 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 + + # Create wrapper to objective function and push it on top of per-thread + # stack before calling optimizer. + fw = ObjFun(f, n, scl) # wrapper to objective function + fp = _push_objfun(newuoa, fw) # pointer to C-callable function try - # Call low-level optimizer. - rc = prima_newuoa(fp, n, x, fx, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) - return (fx[], Int(nf[]), rc) + # Call low-level optimizer on the (scaled) variables. + isempty(scl) || _scale!(x, /, scl) + status = prima_newuoa(fp, n, x, fx, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) + isempty(scl) || _scale!(x, *, scl) + return Info(; fx = fx[], nf = nf[], status) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end """ - PRIMA.uobyqa!(f, x; kwds...) -> (fx, nf, rc) + PRIMA.uobyqa!(f, x; kwds...) -> info::PRIMA.Info in-place version of [`PRIMA.uobyqa`](@ref) which to see for details. On entry, argument `x` is a dense vector of double precision value with the initial @@ -338,40 +459,38 @@ variables; on return, `x` is overwritten by an approximate solution. """ function uobyqa!(f, x::DenseVector{Cdouble}; - rhobeg::Real = 1.0, - rhoend::Real = 1e-4*rhobeg, + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, + rhoend::Real = default_rhoend(rhobeg), ftarget::Real = -Inf, - maxfun::Integer = 100*length(x), + maxfun::Integer = default_maxfun(x), iprint::Union{Integer,Message} = MSG_NONE) # Check arguments. n = length(x) # number of variables _check_rho(rhobeg, rhoend) + scl = _get_scaling(scale, n) + # References for output values. 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 + + # Create wrapper to objective function and push it on top of per-thread + # stack before calling optimizer. + fw = ObjFun(f, n, scl) # wrapper to objective function + fp = _push_objfun(uobyqa, fw) # pointer to C-callable function try - # Call low-level optimizer. - rc = prima_uobyqa(fp, n, x, fx, nf, rhobeg, rhoend, ftarget, maxfun, iprint) - return (fx[], Int(nf[]), rc) + # Call low-level optimizer on the (scaled) variables. + isempty(scl) || _scale!(x, /, scl) + status = prima_uobyqa(fp, n, x, fx, nf, rhobeg, rhoend, ftarget, maxfun, iprint) + isempty(scl) || _scale!(x, *, scl) + return Info(; fx = fx[], nf = nf[], status) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end """ - PRIME.LinearConstraints - -is the type of `(A,b)`, the 2-tuple representing linear equality constraints -`A⋅x = b` or linear inequality constraints `A⋅x ≤ b` where `A` is a matrix, `x` -is the vector of variables, and `b` is a vector. - -""" -const LinearConstraints = Tuple{AbstractMatrix{<:Real},AbstractVector{<:Real}} - -""" - PRIMA.cobyla!(f, x; kwds...) -> (fx, nf, rc) + PRIMA.cobyla!(f, x; kwds...) -> info::PRIMA.Info in-place version of [`PRIMA.cobyla`](@ref) which to see for details. On entry, argument `x` is a dense vector of double precision value with the initial @@ -379,44 +498,73 @@ variables; on return, `x` is overwritten by an approximate solution. """ function cobyla!(f, x::DenseVector{Cdouble}; - nonlinear_ineq::Union{AbstractVector{<:Real},Integer,Nothing} = nothing, - linear_ineq::Union{LinearConstraints,Nothing} = nothing, - linear_eq::Union{LinearConstraints,Nothing} = nothing, + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, + rhoend::Real = default_rhoend(rhobeg), + ftarget::Real = -Inf, + maxfun::Integer = default_maxfun(x), + iprint::Union{Integer,Message} = MSG_NONE, xl::Union{AbstractVector{<:Real},Nothing} = nothing, xu::Union{AbstractVector{<:Real},Nothing} = nothing, - rhobeg::Real = 1.0, - rhoend::Real = 1e-4*rhobeg, - ftarget::Real = -Inf, - maxfun::Integer = 100*length(x), - iprint::Union{Integer,Message} = MSG_NONE) + linear_ineq::Union{LinearConstraints,Nothing} = nothing, + linear_eq::Union{LinearConstraints,Nothing} = nothing, + nonlinear_ineq = nothing, + nonlinear_eq = nothing) # Check arguments and get constraints. n = length(x) # number of variables _check_rho(rhobeg, rhoend) - xl = _get_lower_bound(xl, n) - xu = _get_upper_bound(xu, n) - m_nlcon, nlcon = _get_nonlinear_constraints(nonlinear_ineq) - m_eq, A_eq, b_eq = _get_linear_constraints(linear_eq, n) - m_ineq, A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) - + scl = _get_scaling(scale, n) + xl = _get_lower_bound(xl, n, scl) + xu = _get_upper_bound(xu, n, scl) + nl_eq, c_nl_eq = _get_nonlinear_constraints(nonlinear_eq, x, "equality") + nl_ineq, c_nl_ineq = _get_nonlinear_constraints(nonlinear_ineq, x, "inequality") + A_eq, b_eq = _get_linear_constraints(linear_eq, n, scl) + A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n, scl) + + # Allocate vector to store all non-linear constraints (the non-linear + # equalities being implemented as 2 inequalities each). + nl_all = Vector{Cdouble}(undef, 2*length(nl_eq) + length(nl_ineq)) + + # References for output values. 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 + # Create wrapper to objective function and push it on top of per-thread + # stack before calling optimizer. + fw = ObjFun(f, n, scl, c_nl_eq, length(nl_eq), c_nl_ineq, length(nl_ineq)) + fp = _push_objfun(cobyla, fw) # pointer to C-callable function try - # Call low-level optimizer. - rc = prima_cobyla(m_nlcon, fp, n, x, fx, cstrv, nlcon, - m_ineq, A_ineq, b_ineq, m_eq, A_eq, b_eq, - xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, iprint) - return (fx[], Int(nf[]), rc, cstrv[]) + # Call low-level optimizer on the (scaled) variables. + isempty(scl) || _scale!(x, /, scl) + status = prima_cobyla(length(nl_all), fp, n, x, fx, cstrv, nl_all, + length(b_ineq), A_ineq, b_ineq, + length(b_eq), A_eq, b_eq, + xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, iprint) + isempty(scl) || _scale!(x, *, scl) + # Unpack constraints. + if length(nl_eq) > 0 + i = firstindex(nl_all) + for j in eachindex(nl_eq) + nl_eq[j] = nl_all[i] + i += 2 + end + end + if length(nl_ineq) > 0 + i = firstindex(nl_all) + 2*length(nl_eq) + for j in eachindex(nl_ineq) + nl_ineq[j] = nl_all[i] + i += 1 + end + end + return Info(; fx = fx[], nf = nf[], status, cstrv = cstrv[], nl_eq, nl_ineq) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end """ - PRIMA.lincoa!(f, x; kwds...) -> (fx, nf, rc) + PRIMA.lincoa!(f, x; kwds...) -> info::PRIMA.Info in-place version of [`PRIMA.lincoa`](@ref) which to see for details. On entry, argument `x` is a dense vector of double precision value with the initial @@ -424,114 +572,267 @@ variables; on return, `x` is overwritten by an approximate solution. """ function lincoa!(f, x::DenseVector{Cdouble}; - linear_ineq::Union{LinearConstraints,Nothing} = nothing, - linear_eq::Union{LinearConstraints,Nothing} = nothing, + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, + rhoend::Real = default_rhoend(rhobeg), + ftarget::Real = -Inf, + maxfun::Integer = default_maxfun(x), + npt::Integer = default_npt(x), + iprint::Union{Integer,Message} = MSG_NONE, xl::Union{AbstractVector{<:Real},Nothing} = nothing, xu::Union{AbstractVector{<:Real},Nothing} = nothing, - rhobeg::Real = 1.0, - rhoend::Real = 1e-4*rhobeg, - ftarget::Real = -Inf, - maxfun::Integer = 100*length(x), - npt::Integer = 2*length(x) + 1, - iprint::Union{Integer,Message} = MSG_NONE) + linear_ineq::Union{LinearConstraints,Nothing} = nothing, + linear_eq::Union{LinearConstraints,Nothing} = nothing) # Check arguments and get constraints. n = length(x) # number of variables - _check_rho(rhobeg, rhoend) _check_npt(npt, n) - xl = _get_lower_bound(xl, n) - xu = _get_upper_bound(xu, n) - m_eq, A_eq, b_eq = _get_linear_constraints(linear_eq, n) - m_ineq, A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) + _check_rho(rhobeg, rhoend) + scl = _get_scaling(scale, n) + xl = _get_lower_bound(xl, n, scl) + xu = _get_upper_bound(xu, n, scl) + A_eq, b_eq = _get_linear_constraints(linear_eq, n, scl) + A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n, scl) + # References for output values. 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 + + # Create wrapper to objective function and push it on top of per-thread + # stack before calling optimizer. + fw = ObjFun(f, n, scl) # wrapper to objective function + fp = _push_objfun(lincoa, fw) # pointer to C-callable function try - # Call low-level optimizer. - rc = prima_lincoa(fp, n, x, fx, cstrv, - m_ineq, A_ineq, b_ineq, m_eq, A_eq, b_eq, xl, xu, - nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) - return (fx[], Int(nf[]), rc, cstrv[]) + # Call low-level optimizer on the (scaled) variables. + isempty(scl) || _scale!(x, /, scl) + status = prima_lincoa(fp, n, x, fx, cstrv, + length(b_ineq), A_ineq, b_ineq, + length(b_eq), A_eq, b_eq, + xl, xu, + nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) + isempty(scl) || _scale!(x, *, scl) + return Info(; fx = fx[], nf = nf[], status, cstrv = cstrv[]) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end #------------------------------------------------------------------------------ # PRIVATE METHODS AND TYPES -abstract type AbstractObjFun end +""" + PRIMA.ObjFun(f, n) + +builds an object to wrap a multi-variate user-defined objective function +`f: ℝⁿ → ℝ` with `n` the number of variables. + + PRIMA.ObjFun(n, f, c_eq, n_eq, c_ineq, n_ineq) + +builds an object to wrap a multi-variate user-defined objective function `f: ℝⁿ +→ ℝ` with `n` the number of variables and the `n_eq` and `n_ineq` non-linear +equality and inequality constraints: -# 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 + c_eq(x) = 0 + c_ineq(x) ≤ 0 + +Objective function object `objfun` has the following fields: + + objfun.f # callable implementing user-defined objective function + objfun.n # number of variables + objfun.c_eq # callable implementing non-linear equalities as `c_eq(x) = 0` + objfun.n_eq # number of non-linear equalities + objfun.c_ineq # callable implementing non-linear inequalities as `c_ineq(x) ≤ 0` + objfun.n_ineq # number of non-linear inequalities + +See also [`PRIMA.call`](@ref) and [`PRIMA.call!`](@ref). + +""" +struct ObjFun{F,E,I} + f::F # callable implementing user-defined objective function + n::Int # number of variables + c_eq::E # callable implementing non-linear equalities as `c_eq(x) = 0` + n_eq::Int # number of non-linear equalities + c_ineq::I # callable implementing non-linear inequalities as `c_ineq(x) ≤ 0` + n_ineq::Int # number of non-linear inequalities + scl::Vector{Cdouble} # scaling factors or empty + x::Vector{Cdouble} # workspace for variables end -# Small structure to wrap objective functions with constraints for COBYLA at -# low level. -struct ObjFunCon <: AbstractObjFun - f::Any # user-defined objective function - nx::Int # number of variables - nc::Int # number of non-linear constraints +unconstrained(x::AbstractVector{T}) where {T} = NullVector{T}() + +ObjFun(f, n::Integer, scl::AbstractVector) = + ObjFun(f, n, scl, unconstrained, 0, unconstrained, 0) + +ObjFun(f, n::Integer, scl::AbstractVector, c_eq, n_eq::Integer, c_ineq, n_ineq::Integer) = + ObjFun(f, n, c_eq, n_eq, c_ineq, n_ineq, scl, Vector{Cdouble}(undef, n)) + +""" + PRIMA.call(f::ObjFun, x::DenseVector{Cdouble}) -> fx + +yields value of objective function `f(x)` for variables `x`. + +The default method may be extended by specializing on the type of `f`. It is +guaranteed that: + + length(x) = f.n + +holds. + +See also [`PRIMA.ObjFun`](@ref) and [`PRIMA.call!`](@ref). + +""" +function call(f::ObjFun, x::DenseVector{Cdouble}) + length(x) == f.n || throw(DimensionMismatch( + "invalid number of variables")) + return f.f(x) +end + +""" + PRIMA.call!(f::ObjFun, x::DenseVector{Cdouble}, cx::DenseVector{Cdouble}) -> fx + +yields value of objective function `f(x)` for variables `x` and overwrites `cx` +with the values `c(x)` corresponding to the non-linear inequality constraints +`c(x) ≤ 0`. + +The default method may be extended by specializing on the type of `f`. It is +guaranteed that: + + length(x) = f.n + length(cx) = 2*f.n_eq + f.n_ineq + +both hold. The second equality is because non-linear equality constraints are +equivalent to two non-linear equality constraints: `c_eq(x) = 0` is rewritten +as `c_eq(x) ≤ 0` and `-c_eq(x) ≤ 0`. + +See also [`PRIMA.ObjFun`](@ref) and [`PRIMA.call`](@ref). + +""" +function call!(f::ObjFun, x::DenseVector{Cdouble}, cx::DenseVector{Cdouble}) + length(x) == f.n || throw(DimensionMismatch( + "invalid number of variables")) + i = firstindex(cx) - 1 + if f.n_eq != 0 + c_eq = f.c_eq(x)::Union{Real,AbstractVector{<:Real}} + length(c_eq) == f.n_eq || throw(DimensionMismatch( + "invalid number of equalities")) + for j in eachindex(c_eq) + cx[i += 1] = c_eq[j] + cx[i += 1] = -c_eq[j] + end + end + if f.n_ineq != 0 + c_ineq = f.c_ineq(x)::Union{Real,AbstractVector{<:Real}} + length(c_ineq) == f.n_ineq || throw(DimensionMismatch( + "invalid number of inequalities")) + for j in eachindex(c_ineq) + cx[i += 1] = c_ineq[j] + end + end + return f.f(x) +end + +# Computation of objective function and non-linear constraints, if any, is done +# in several stages: +# +# 1. `_objfun` retrieves the objective function; +# +# 2. `unsafe_call` wraps pointers to the variables and the constraints, if any, +# into Julia arrays; +# +# 3. `call` or `call!` execute the user-defined functions implementing the +# objective function and the non-linear constraints. +# +# The motivations are (i) to dispatch as soon as possible on code depending on +# the type of the user-defined Julia functions and (ii) to make possible to +# extend `call` or `call!` at the last stage . + +# C-callable objective function for problems with no non-linear constraints +# (for other algorithms than COBYLA). +function _objfun(x_ptr::Ptr{Cdouble}, # (input) variables + f_ptr::Ptr{Cdouble}) # (output) function value + # Retrieve objective function object and dispatch on its type to compute + # f(x). + f = last(_objfun_stack[Threads.threadid()]) + unsafe_store!(f_ptr, unsafe_call(f, x_ptr)) + return nothing +end + +# C-callable objective function for problems with non-linear constraints (for +# COBYLA algorithm). +function _objfun(x_ptr::Ptr{Cdouble}, # (input) variables + f_ptr::Ptr{Cdouble}, # (output) function value + c_ptr::Ptr{Cdouble}) # (output) constraints + # Retrieve objective function object and dispatch on its type to compute + # f(x) and the non-linear constraints. + f = last(_objfun_stack[Threads.threadid()]) + unsafe_store!(f_ptr, unsafe_call(f, x_ptr, c_ptr)) + return nothing end +function unsafe_call(f::ObjFun, x_ptr::Ptr{Cdouble}) + x = get_variables(f, x_ptr) + return as(Cdouble, call(f, x)) +end + +function unsafe_call(f::ObjFun, x_ptr::Ptr{Cdouble}, c_ptr::Ptr{Cdouble}) + x = get_variables(f, x_ptr) + c = unsafe_wrap(Array, c_ptr, 2*f.n_eq + f.n_ineq) + return as(Cdouble, call!(f, x, c)) +end + +function get_variables(f::ObjFun, x_ptr::Ptr{Cdouble}) + n, x, scl = f.n, f.x, f.scl + length(x) == n || corrupted_structure("invalid number of variables") + if isempty(scl) + # No scaling of variables. + @inbounds for i in 1:n + x[i] = unsafe_load(x_ptr, i) + end + elseif length(scl) == n + # Scale variables. + @inbounds for i in 1:n + x[i] = scl[i]*unsafe_load(x_ptr, i) + end + else + corrupted_structure("invalid number of scaling factors") + end + return x +end + +@noinline corrupted_structure(msg::AbstractString) = + throw(AssertionError("corrupted structure ($msg)")) + # 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. +# call, whatever happens. It is therefore necessary to wrap the 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 objective 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} +# Private function `_get_objfun_stack` yields the stack of objective functions +# for the caller thread. +function _get_objfun_stack() i = Threads.threadid() - while length(stack) < i - push!(stack, Vector{T}(undef, 0)) + while length(_objfun_stack) < i + push!(_objfun_stack, Vector{ObjFun}(undef, 0)) end - return stack[i] + return _objfun_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 +# Private functions `_push_objfun` and `_pop_objfun` 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) +function _push_objfun(algorithm, fw::ObjFun) + _objfun_ptr(::Any) = + @cfunction(_objfun, Cvoid, (Ptr{Cdouble}, Ptr{Cdouble})) + _objfun_ptr(::typeof(cobyla)) = + @cfunction(_objfun, Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble})) + push!(_get_objfun_stack(), fw) + return _objfun_ptr(algorithm) end -function _pop_wrapper(fw::AbstractObjFun) - stack = _get_stack(fw) + +function _pop_objfun(fw::ObjFun) + stack = _get_objfun_stack() last(stack) === fw || error( "objective function is not the last one in the caller thread stask") resize!(stack, length(stack) - 1) @@ -573,14 +874,14 @@ Base.unsafe_convert(::Type{Ptr{S}}, ::NullArray{T,N}) where {T,N,S<:Union{Cvoid, # 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, NullMatrix{Cdouble}(), NullVector{Cdouble}() -function _get_linear_constraints(Ab::LinearConstraints, n::Integer) +_get_linear_constraints(::Nothing, n::Integer, scl::AbstractVector) = + NullMatrix{Cdouble}(), NullVector{Cdouble}() +function _get_linear_constraints(Ab::LinearConstraints, n::Integer, scl::AbstractVector) A, b = Ab - Base.has_offset_axes(A) && error( - "matrix `A` of linear constraints must have 1-based indices") - Base.has_offset_axes(b) && error( - "vector `b` of linear constraints must have 1-based indices") + Base.has_offset_axes(A) && throw(ArgumentError( + "matrix `A` of linear constraints must have 1-based indices")) + Base.has_offset_axes(b) && throw(ArgumentError( + "vector `b` of linear constraints must have 1-based indices")) m = length(b) # number of constraints size(A) == (m,n) || throw(DimensionMismatch( "matrix `A` of linear constraints has incompatible dimensions")) @@ -592,33 +893,69 @@ function _get_linear_constraints(Ab::LinearConstraints, n::Integer) # twice. This isn't a big issue for a small number of variables and # constraints, but it's not completely satisfactory either. A_ = Matrix{T}(undef, n, m) - @inbounds for i ∈ 1:m - for j ∈ 1:n - A_[j,i] = A[i,j] + if isempty(scl) + # No scaling of the variables. + @inbounds for i ∈ 1:m + for j ∈ 1:n + A_[j,i] = A[i,j] + end + end + else + # Scaling matrix of linear constraints to account for the scaling of + # the variables. + @inbounds for i ∈ 1:m + for j ∈ 1:n + A_[j,i] = A[i,j]*scl[j] + end end end b_ = _dense_array(T, b) - return m, A_, b_ + return A_, b_ end -_get_nonlinear_constraints(::Nothing) = - 0, NullVector{Cdouble}() -_get_nonlinear_constraints(m::Integer) = - _get_nonlinear_constraints(Vector{Cdouble}(undef, m)) -_get_nonlinear_constraints(c::AbstractVector{<:Real}) = - length(c), _dense_array(Cdouble, c) +# Yield (v,c) with v a vector to store the constraints on output (possibly an +# empty vector if there are no constraints), and c the callable object +# implementing the constraints. +_get_nonlinear_constraints(::Nothing, x::AbstractArray, str::AbstractString) = + Cdouble[], unconstrained +_get_nonlinear_constraints(c::Tuple{Integer,Any}, x::AbstractArray, str::AbstractString) = + Vector{Cdouble}(undef, c[1]), c[2] +_get_nonlinear_constraints(c::Tuple{Any,Integer}, x::AbstractArray, str::AbstractString) = + _get_nonlinear_constraints(reverse(c), x, str) +function _get_nonlinear_constraints(c::Any, x::AbstractArray, str::AbstractString) + # If c(x) does not yield a result convertible into a vector of double precisions + # values, an error will be thrown to the caller. + cx = c(x) + cx isa Real && return [as(Cdouble, cx)], c + cx isa AbstractVector{<:Real} && return as(Vector{Cdouble}, cx), c + throw(ArgumentError(string("method implementing non-linear ", str, + " constraints must return a real or a vector of reals"))) +end +# Encode _get_lower_bound and _get_upper_bound. These methods assume that scl +# and n have been checked. for (uplo, def) in ((:lower, typemin), (:upper, typemax)) func = Symbol("_get_$(uplo)_bound") @eval begin - $func(::Nothing, n::Integer) = fill($def(Cdouble), n) - function $func(b::AbstractVector, n::Integer) + $func(::Nothing, n::Integer, scl::AbstractVector) = fill($def(Cdouble), n) + function $func(b::AbstractVector, n::Integer, scl::AbstractVector) Base.has_offset_axes(b) && error( $("$uplo bound must have 1-based indices")) length(b) == n || throw(DimensionMismatch(string( $("$uplo bound must have "), n, " elements"))) - return _dense_array(Cdouble, b) + if length(scl) == 0 + # No scaling of variables. For type-stability, ensure that an + # ordinary vector is returned. + return convert(Vector{Cdouble}, b) + else + # Modify the bounds to account for the scaling of the variables. + b_ = Vector{Cdouble}(undef, n) + @inbounds for i in 1:n + b_[i] = b[i]/scl[i] + end + return b_ + end end end end @@ -628,4 +965,24 @@ end _dense_array(::Type{T}, A::DenseArray{T,N}) where {T,N} = A _dense_array(::Type{T}, A::AbstractArray{<:Any,N}) where {N,T} = convert(Array{T,N}, A) +# Scale variables. +function _scale!(x::AbstractVector, op, scl::AbstractVector) + @inbounds for i in eachindex(x, scl) + x[i] = op(x[i], scl[i]) + end + return nothing end + +# Yield scl, rhobeg, and rhoend given the arguments. +_get_scaling(scl::Nothing, n::Int) = Cdouble[] +function _get_scaling(scl::AbstractVector{<:Real}, n::Int) + Base.has_offset_axes(scl) && throw(ArgumentError( + "vector of scaling factors must have 1-based indices")) + length(scl) == n || throw(ArgumentError( + "vector of scaling factors have incompatible number of elements")) + all(x -> isfinite(x) & (x > zero(x)), scl) || throw(ArgumentError( + "scaling factor(s) must be finite and positive")) + return convert(Vector{Cdouble}, scl) +end + +end # module diff --git a/src/wrappers.jl b/src/wrappers.jl index b35acb2..ae95c06 100644 --- a/src/wrappers.jl +++ b/src/wrappers.jl @@ -1,11 +1,11 @@ -@enum Message::Int32 begin +@enum Message::Cint begin MSG_NONE = 0 MSG_EXIT = 1 MSG_RHO = 2 MSG_FEVL = 3 end -@enum Status::Int32 begin +@enum Status::Cint begin SMALL_TR_RADIUS = 0 FTARGET_ACHIEVED = 1 TRSUBP_FAILED = 2 @@ -39,21 +39,21 @@ function prima_bobyqa(calfun, n, x, f, xl, xu, nf, rhobeg, rhoend, ftarget, maxf f::Ptr{Cdouble}, xl::Ptr{Cdouble}, xu::Ptr{Cdouble}, nf::Ptr{Cint}, rhobeg::Cdouble, rhoend::Cdouble, ftarget::Cdouble, maxfun::Cint, npt::Cint, - iprint::Cint)::Cint + iprint::Cint)::Status end function prima_newuoa(calfun, n, x, f, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) @ccall libprimac.prima_newuoa(calfun::prima_obj, n::Cint, x::Ptr{Cdouble}, f::Ptr{Cdouble}, nf::Ptr{Cint}, rhobeg::Cdouble, rhoend::Cdouble, ftarget::Cdouble, maxfun::Cint, - npt::Cint, iprint::Cint)::Cint + npt::Cint, iprint::Cint)::Status end function prima_uobyqa(calfun, n, x, f, nf, rhobeg, rhoend, ftarget, maxfun, iprint) @ccall libprimac.prima_uobyqa(calfun::prima_obj, n::Cint, x::Ptr{Cdouble}, f::Ptr{Cdouble}, nf::Ptr{Cint}, rhobeg::Cdouble, rhoend::Cdouble, ftarget::Cdouble, maxfun::Cint, - iprint::Cint)::Cint + iprint::Cint)::Status end function prima_cobyla(m_nlcon, calcfc, n, x, f, cstrv, nlconstr, m_ineq, Aineq, bineq, m_eq, @@ -64,7 +64,7 @@ function prima_cobyla(m_nlcon, calcfc, n, x, f, cstrv, nlconstr, m_ineq, Aineq, bineq::Ptr{Cdouble}, m_eq::Cint, Aeq::Ptr{Cdouble}, beq::Ptr{Cdouble}, xl::Ptr{Cdouble}, xu::Ptr{Cdouble}, nf::Ptr{Cint}, rhobeg::Cdouble, rhoend::Cdouble, - ftarget::Cdouble, maxfun::Cint, iprint::Cint)::Cint + ftarget::Cdouble, maxfun::Cint, iprint::Cint)::Status end function prima_lincoa(calfun, n, x, f, cstrv, m_ineq, Aineq, bineq, m_eq, Aeq, beq, xl, xu, @@ -75,5 +75,5 @@ function prima_lincoa(calfun, n, x, f, cstrv, m_ineq, Aineq, bineq, m_eq, Aeq, b Aeq::Ptr{Cdouble}, beq::Ptr{Cdouble}, xl::Ptr{Cdouble}, xu::Ptr{Cdouble}, nf::Ptr{Cint}, rhobeg::Cdouble, rhoend::Cdouble, ftarget::Cdouble, maxfun::Cint, - npt::Cint, iprint::Cint)::Cint + npt::Cint, iprint::Cint)::Status end diff --git a/test/runtests.jl b/test/runtests.jl index 38088c9..65d0ff1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ optimizer_name(::typeof(PRIMA.newuoa)) = "NEWUOA" optimizer_name(::typeof(PRIMA.bobyqa)) = "BOBYQA" optimizer_name(::typeof(PRIMA.cobyla)) = "COBYLA" optimizer_name(::typeof(PRIMA.lincoa)) = "LINCOA" +optimizer_name(::typeof(PRIMA.prima)) = "PRIMA" optimizer_name(algo::Symbol) = optimizer_name(optimizer(algo)) optimizer(algo::Symbol) = @@ -14,8 +15,27 @@ optimizer(algo::Symbol) = algo === :bobyqa ? PRIMA.bobyqa : algo === :cobyla ? PRIMA.cobyla : algo === :lincoa ? PRIMA.lincoa : + algo === :prima ? PRIMA.prima : error("unknown optimizer `:$algo`") +print_1(x::AbstractVector, info::PRIMA.Info) = print_1(stdout, x, info) +function print_1(io::IO, x::AbstractVector, info::PRIMA.Info) + msg = PRIMA.reason(info) + println(io, "x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") +end + +print_2(x::AbstractVector, info::PRIMA.Info) = print_2(stdout, x, info) +function print_2(io::IO, x::AbstractVector, info::PRIMA.Info) + msg = PRIMA.reason(info) + println(io, "x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), status = $(info.status), msg = '$msg', evals = $(info.nf)") +end + +print_3(x::AbstractVector, info::PRIMA.Info) = print_3(stdout, x, info) +function print_3(io::IO, x::AbstractVector, info::PRIMA.Info) + msg = PRIMA.reason(info) + println("x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), c(x) = $(info.nl_ineq), status = $(info.status), msg = '$msg', evals = $(info.nf)") +end + @testset "PRIMA.jl" begin @testset "Utils " begin for status in instances(PRIMA.Status) @@ -40,29 +60,19 @@ optimizer(algo::Symbol) = 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} - c!(x, c) - return f(x) - end - - # Non-linear constraint. - function c!(x::AbstractVector{T}, c::AbstractVector{T}) where {T<:AbstractFloat} + # Non-linear inequality constraints. + function c_ineq(x::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 + return x1*x1 + x2*x2 - 13 # ‖x‖² ≤ 13 end - # Array to store non-linear constraints. - c = 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 + scl = 2.0 # scaling factor, a power of two should not change results + scale = [scl, scl] # Inequality constraints: x₁ ≤ 4, x₂ ≤ 3, x₁ + x₂ ≤ 10 A_ineq = [1.0 0.0; @@ -84,100 +94,115 @@ optimizer(algo::Symbol) = @testset "NEWUOA" begin println("\nNEWUOA:") - 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") + kwds = (rhobeg = 1.0, rhoend = 1e-3, ftarget = -Inf, + maxfun = 200n, npt = 2n + 1, iprint = PRIMA.MSG_EXIT) + x, info = @inferred PRIMA.newuoa(f, x0; kwds...) + print_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx + @test f(x) ≈ info.fx @test x0 == x0_sav + # Solve problem with general driver. + x1, info1 = @inferred PRIMA.prima(f, x0; kwds...) + @test x1 == x + @test info1 == info + # Solve problem with scaling factors. + kwds = (scale, rhobeg = 1.0/scl, rhoend = 1e-3/scl, ftarget = -Inf, + maxfun = 200n, npt = 2n + 1, iprint = PRIMA.MSG_EXIT) + x1, info1 = @inferred PRIMA.newuoa(f, x0; kwds...) + @test x1 ≈ x end @testset "UOBYQA" begin println("\nUOBYQA:") - 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") + kwds = (rhobeg = 1.0, rhoend = 1e-3, ftarget = -Inf, + maxfun = 200n, iprint = PRIMA.MSG_EXIT) + x, info = @inferred PRIMA.uobyqa(f, x0; kwds...) + print_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx + @test f(x) ≈ info.fx @test x0 == x0_sav + # Solve problem with scaling factors. + kwds = (scale, rhobeg = 1.0/scl, rhoend = 1e-3/scl, ftarget = -Inf, + maxfun = 200n, iprint = PRIMA.MSG_EXIT) + x1, info1 = @inferred PRIMA.uobyqa(f, x0; kwds...) + @test x1 ≈ x end @testset "BOBYQA" begin println("\nBOBYQA:") - 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") + kwds = (xl = xl, xu = xu, + rhobeg = 1.0, rhoend = 1e-3, ftarget = -Inf, + maxfun = 200n, npt = 2n + 1, iprint = PRIMA.MSG_EXIT) + x, info = @inferred PRIMA.bobyqa(f, x0; kwds...) + print_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx + @test f(x) ≈ info.fx @test x0 == x0_sav @test check_bounds(xl, x, xu) + # Solve problem with general driver. + x1, info1 = @inferred PRIMA.prima(f, x0; kwds...) + @test x1 == x + @test info1 == info + # Solve problem with scaling factors. + kwds = (scale, rhobeg = 1.0/scl, rhoend = 1e-3/scl, ftarget = -Inf, + maxfun = 200n, npt = 2n + 1, iprint = PRIMA.MSG_EXIT) + x1, info1 = @inferred PRIMA.bobyqa(f, x0; kwds...) + @test x1 ≈ x end @testset "COBYLA" begin println("\nCOBYLA:") + kwds = (xl = xl, xu = xu, linear_ineq = (A_ineq, b_ineq), + rhobeg = 1.0, rhoend = 1e-3, ftarget = -Inf, + maxfun = 200*n, iprint = PRIMA.MSG_EXIT) # First call with just the number of non-linear inequality constraints. - x, fx, nf, rc, cstrv = @inferred PRIMA.cobyla(f, x0; - nonlinear_ineq = length(c), - linear_ineq = (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, rc = $rc, msg = '$msg', evals = $nf") - @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx - @test x0 == x0_sav - @test check_bounds(xl, x, xu) - # Second call with an array to store the non-linear inequality constraints. - x, fx, nf, rc, cstrv = @inferred PRIMA.cobyla(f, x0; - nonlinear_ineq = c, - linear_ineq = (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, c(x) = $c, rc = $rc, msg = '$msg', evals = $nf") + x, info = @inferred PRIMA.cobyla(f, x0; kwds..., + nonlinear_ineq = c_ineq) + print_3(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx + @test f(x) ≈ info.fx @test x0 == x0_sav @test check_bounds(xl, x, xu) + # Call with given number of non-linear inequality constraints. + x1, info1 = @inferred PRIMA.cobyla(f, x0; kwds..., + nonlinear_ineq = (length(c_ineq(x0)), c_ineq)) + @test x1 == x + @test info1 == info + # Solve problem with general driver. + x1, info1 = @inferred PRIMA.prima(f, x0; kwds..., + nonlinear_ineq = c_ineq) + @test x1 == x + @test info1 == info + # Solve problem with scaling factors. + kwds = (xl = xl, xu = xu, linear_ineq = (A_ineq, b_ineq), + scale, rhobeg = 1.0/scl, rhoend = 1e-3/scl, ftarget = -Inf, + maxfun = 200n, iprint = PRIMA.MSG_EXIT) + x1, info1 = @inferred PRIMA.cobyla(f, x0; kwds..., + nonlinear_ineq = c_ineq) + @test x1 ≈ x end @testset "LINCOA" begin println("\nLINCOA:") - x, fx, nf, rc, cstrv = @inferred PRIMA.lincoa(f, x0; - linear_ineq = (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, rc = $rc, msg = '$msg', evals = $nf") + kwds = (xl = xl, xu = xu, linear_ineq = (A_ineq, b_ineq), + rhobeg = 1.0, rhoend = 1e-3, ftarget = -Inf, + maxfun = 200*n, npt = 2n + 1, iprint = PRIMA.MSG_EXIT) + x, info = @inferred PRIMA.lincoa(f, x0; kwds...) + print_2(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx + @test f(x) ≈ info.fx @test x0 == x0_sav @test check_bounds(xl, x, xu) + # Solve problem with general driver. + x1, info1 = @inferred PRIMA.prima(f, x0; kwds...) + @test x1 == x + @test info1 == info + # Solve problem with scaling factors. + kwds = (xl = xl, xu = xu, linear_ineq = (A_ineq, b_ineq), + scale, rhobeg = 1.0/scl, rhoend = 1e-3/scl, ftarget = -Inf, + maxfun = 200n, npt = 2n + 1, iprint = PRIMA.MSG_EXIT) + x1, info1 = @inferred PRIMA.lincoa(f, x0; kwds...) + @test x1 ≈ x end end @@ -242,27 +267,21 @@ optimizer(algo::Symbol) = println("\nUnconstrained minimization of Rosenbrock function by $(optimizer_name(optim)):") x0 = [-1, 2] if optim === :uobyqa - x, fx, nf, rc = @inferred uobyqa(f, x0; - rhobeg, rhoend, ftarget, maxfun, iprint) + x, info = @inferred uobyqa(f, x0; rhobeg, rhoend, ftarget, maxfun, iprint) elseif optim === :newuoa - x, fx, nf, rc = @inferred newuoa(f, x0; - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred newuoa(f, x0; rhobeg, rhoend, ftarget, maxfun, iprint, npt) elseif optim === :bobyqa - x, fx, nf, rc = @inferred bobyqa(f, x0; - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred bobyqa(f, x0; rhobeg, rhoend, ftarget, maxfun, iprint, npt) elseif optim === :cobyla - x, fx, nf, rc = @inferred cobyla(f, x0; - rhobeg, rhoend, ftarget, maxfun, iprint) + x, info = @inferred cobyla(f, x0; rhobeg, rhoend, ftarget, maxfun, iprint) elseif optim === :lincoa - x, fx, nf, rc = @inferred lincoa(f, x0; - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred lincoa(f, x0; rhobeg, rhoend, ftarget, maxfun, iprint, npt) else continue end - msg = PRIMA.reason(rc) - println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + print_1(x, info) @test x ≈ [1,1] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) - @test f(x) ≈ fx + @test f(x) ≈ info.fx # Bound constrained optimization. if optim ∈ (:bobyqa, :cobyla, :lincoa) @@ -271,24 +290,20 @@ optimizer(algo::Symbol) = xl = [-Inf, 1.2] xu = [+Inf, +Inf] if optim === :bobyqa - x, fx, nf, rc = @inferred bobyqa(f, x0; - xl, xu, - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred bobyqa(f, x0; xl, xu, + rhobeg, rhoend, ftarget, maxfun, iprint, npt) elseif optim === :cobyla - x, fx, nf, rc = @inferred cobyla(f, x0; - xl, xu, - rhobeg, rhoend, ftarget, maxfun, iprint) + x, info = @inferred cobyla(f, x0; xl, xu, + rhobeg, rhoend, ftarget, maxfun, iprint) elseif optim === :lincoa - x, fx, nf, rc = @inferred lincoa(f, x0; - xl, xu, - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred lincoa(f, x0; xl, xu, + rhobeg, rhoend, ftarget, maxfun, iprint, npt) else continue end - msg = PRIMA.reason(rc) - println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + print_1(x, info) @test x ≈ [1.095247,1.2] rtol=0 atol=2e-2 - @test f(x) ≈ fx + @test f(x) ≈ info.fx end # Linearly constrained optimization. @@ -297,40 +312,34 @@ optimizer(algo::Symbol) = x0 = [-1, 2] # starting point linear_ineq = (A_ineq, b_ineq) if optim === :cobyla - x, fx, nf, rc = @inferred cobyla(f, x0; - linear_ineq, - rhobeg, rhoend, ftarget, maxfun, iprint) + x, info = @inferred cobyla(f, x0; linear_ineq, + rhobeg, rhoend, ftarget, maxfun, iprint) elseif optim === :lincoa - x, fx, nf, rc = @inferred lincoa(f, x0; - linear_ineq, - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred lincoa(f, x0; linear_ineq, + rhobeg, rhoend, ftarget, maxfun, iprint, npt) else continue end - msg = PRIMA.reason(rc) - println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + print_1(x, info) @test x ≈ [1.0,1.0] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) - @test f(x) ≈ fx + @test f(x) ≈ info.fx println("\nIdem but with one linear inequality constraint replaced by a bound constraint:") x0 = [1, 2] # starting point linear_ineq = (A_ineq[1:2,:], b_ineq[1:2]) xl = [-1,-Inf] if optim === :cobyla - x, fx, nf, rc = @inferred cobyla(f, x0; - xl, linear_ineq, - rhobeg, rhoend, ftarget, maxfun, iprint) + x, info = @inferred cobyla(f, x0; xl, linear_ineq, + rhobeg, rhoend, ftarget, maxfun, iprint) elseif optim === :lincoa - x, fx, nf, rc = @inferred lincoa(f, x0; - xl, linear_ineq, - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred lincoa(f, x0; xl, linear_ineq, + rhobeg, rhoend, ftarget, maxfun, iprint, npt) else continue end - msg = PRIMA.reason(rc) - println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + print_1(x, info) @test x ≈ [1.0,1.0] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) - @test f(x) ≈ fx + @test f(x) ≈ info.fx # The solution is on the line: 4x + 3y = 12 (the boundary of the first constraint). println("\nIdem but one linear constraint is active at the solution:") @@ -338,20 +347,17 @@ optimizer(algo::Symbol) = linear_ineq = (A_ineq, b_ineq) xl = [-1,-Inf] if optim === :cobyla - x, fx, nf, rc = @inferred cobyla(f2, x0; - linear_ineq, - rhobeg, rhoend, ftarget, maxfun, iprint) + x, info = @inferred cobyla(f2, x0; linear_ineq, + rhobeg, rhoend, ftarget, maxfun, iprint) elseif optim === :lincoa - x, fx, nf, rc = @inferred lincoa(f2, x0; - linear_ineq, - rhobeg, rhoend, ftarget, maxfun, iprint, npt) + x, info = @inferred lincoa(f2, x0; linear_ineq, + rhobeg, rhoend, ftarget, maxfun, iprint, npt) else continue end - msg = PRIMA.reason(rc) - println("x = $x, f(x) = $fx, rc = $rc, msg = '$msg', evals = $nf") + print_1(x, info) @test x ≈ [1.441832,2.077557] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) - @test f2(x) ≈ fx + @test f2(x) ≈ info.fx end end end