From b99ce9ca880b859e2fc3a60cc740cdab3ecedceb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Fri, 20 Oct 2023 09:39:35 +0200 Subject: [PATCH 01/15] Non-linear equality constraints in cobyla --- NEWS.md | 7 + README.md | 19 ++- src/PRIMA.jl | 369 ++++++++++++++++++++++++++++++++++------------- test/runtests.jl | 19 +-- 4 files changed, 298 insertions(+), 116 deletions(-) diff --git a/NEWS.md b/NEWS.md index f8b7c17..245ecf0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,12 @@ # User visible changes in `PRIMA.jl` +## Version 0.1.2 + +- 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)` as for other algorithms. + - The functions implementing the non-linear constraints are passed by + keywords `nonlinear_eq` and `nonlinear_ineq`. ## Version 0.1.1 diff --git a/README.md b/README.md index daed0bb..7706ce6 100644 --- a/README.md +++ b/README.md @@ -140,6 +140,7 @@ 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` | @@ -174,11 +175,19 @@ Assuming `n = length(x)` is the number of variables, then: elementwise lower and upper bounds for the variables. Feasible variables are such that `xl ≤ x ≤ xu` (elementwise). -- `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_eq` (default `nothing`) may be specified with a function, say + `c_eq`, implementing `n_eq` non-linear equality constraints defined by + `c_eq(x) = 0`. If the caller is interested in the values of `c_eq(x)` at the + returned solution, the keyword may be set with a 2-tuple `(v_eq, c_eq)` or + `(c_eq, v_eq)` with `v_eq` a vector of `n_eq` floating-point values to store + `c_eq(x)`. + +- `nonlinear_ineq` (default `nothing`) may be specified with a function, say + `c_ineq`, implementing `n_ineq` non-linear inequality constraints defined by + `c_ineq(x) ≤ 0`. If the caller is interested in the values of `c_ineq(x)` at + the returned solution, the keyword may be set with a 2-tuple `(v_ineq, + c_ineq)` or `(c_ineq, v_ineq)` with `v_ineq` a vector of `n_ineq` + floating-point values to store `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 diff --git a/src/PRIMA.jl b/src/PRIMA.jl index 22fa002..aec31b3 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -81,11 +81,20 @@ 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 `n_eq` non-linear equality constraints defined by + `c_eq(x) = 0`. If the caller is interested in the values of `c_eq(x)` at the + returned solution, the keyword may be set with a 2-tuple `(v_eq, c_eq)` or + `(c_eq, v_eq)` with `v_eq` a vector of `n_eq` floating-point values to store + `c_eq(x)`. + +- `nonlinear_ineq` (default `nothing`) may be specified with a function, say + `c_ineq`, implementing `n_ineq` non-linear inequality constraints defined by + `c_ineq(x) ≤ 0`. If the caller is interested in the values of `c_ineq(x)` at + the returned solution, the keyword may be set with a 2-tuple `(v_ineq, + c_ineq)` or `(c_ineq, v_ineq)` with `v_ineq` a vector of `n_ineq` + floating-point values to store `c_ineq(x)`. + """ const _doc_linear_constraints = """ @@ -203,14 +212,10 @@ 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 +The objective function takes a single argument, the variables `x`, and returns +the function value, it shall implement the following signature: -where the e,tries of `cx` are to be overwritten by the non-linear consttaints -`c(x)`. + f(x::Vector{Cdouble})::Real Allowed keywords are (`n = length(x)` is the number of variables): @@ -283,16 +288,20 @@ function bobyqa!(f, x::DenseVector{Cdouble}; xl = _get_lower_bound(xl, n) xu = _get_upper_bound(xu, 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) # 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) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end @@ -316,16 +325,20 @@ function newuoa!(f, x::DenseVector{Cdouble}; _check_rho(rhobeg, rhoend) _check_npt(npt, 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) # 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) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end @@ -347,16 +360,20 @@ function uobyqa!(f, x::DenseVector{Cdouble}; n = length(x) # number of variables _check_rho(rhobeg, rhoend) + # 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) # 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) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end @@ -379,7 +396,8 @@ variables; on return, `x` is overwritten by an approximate solution. """ function cobyla!(f, x::DenseVector{Cdouble}; - nonlinear_ineq::Union{AbstractVector{<:Real},Integer,Nothing} = nothing, + nonlinear_ineq = nothing, + nonlinear_eq = nothing, linear_ineq::Union{LinearConstraints,Nothing} = nothing, linear_eq::Union{LinearConstraints,Nothing} = nothing, xl::Union{AbstractVector{<:Real},Nothing} = nothing, @@ -389,29 +407,52 @@ function cobyla!(f, x::DenseVector{Cdouble}; ftarget::Real = -Inf, maxfun::Integer = 100*length(x), iprint::Union{Integer,Message} = MSG_NONE) - # Check arguments and get constaints. + # 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) + n_nl_eq, v_nl_eq, c_nl_eq = _get_nonlinear_constraints(nonlinear_eq, x, "equality") + n_nl_ineq, v_nl_ineq, c_nl_ineq = _get_nonlinear_constraints(nonlinear_ineq, x, "inequality") + n_lin_eq, A_eq, b_eq = _get_linear_constraints(linear_eq, n) + n_lin_ineq, A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) + # Total number of non-linear inequalities and vector to store them. + n_nl = 2*n_nl_eq + n_nl_ineq + v_nl = Vector{Cdouble}(undef, n_nl) + + # 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, c_nl_eq, n_nl_eq, c_nl_ineq, n_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, + rc = prima_cobyla(n_nl, fp, n, x, fx, cstrv, v_nl, + n_lin_ineq, A_ineq, b_ineq, n_lin_eq, A_eq, b_eq, xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, iprint) + # Unpack constraints. + if length(v_nl_eq) > 0 + i = firstindex(v_nl) + for j in eachindex(v_nl_eq) + v_nl_eq[j] = v_nl[i] + i += 2 + end + end + if length(v_nl_ineq) > 0 + i = firstindex(v_nl) + 2*n_nl_eq + for j in eachindex(v_nl_ineq) + v_nl_ineq[j] = v_nl[i] + i += 1 + end + end return (fx[], Int(nf[]), rc, cstrv[]) finally - _pop_wrapper(fw) + _pop_objfun(fw) end end @@ -434,104 +475,223 @@ function lincoa!(f, x::DenseVector{Cdouble}; maxfun::Integer = 100*length(x), npt::Integer = 2*length(x) + 1, iprint::Union{Integer,Message} = MSG_NONE) - # Check arguments and get constaints. + # 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) + n_lin_eq, A_eq, b_eq = _get_linear_constraints(linear_eq, n) + n_lin_ineq, A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) + # 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) # 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, + n_lin_ineq, A_ineq, b_ineq, n_lin_eq, A_eq, b_eq, xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) return (fx[], Int(nf[]), rc, 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: + + 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 +end + +unconstrained(x::AbstractVector{T}) where {T} = NullVector{T}() + +ObjFun(f::F, n::Integer) where {F} = ObjFun(f, n, unconstrained, 0, unconstrained, 0) + +""" + 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. -# 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 +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 -# 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 +""" + 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)` coresponding 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 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 = unsafe_wrap(Array, x_ptr, f.n) + return as(Cdouble, call(f, x)) +end + +function unsafe_call(f::ObjFun, x_ptr::Ptr{Cdouble}, c_ptr::Ptr{Cdouble}) + x = unsafe_wrap(Array, x_ptr, f.n) + c = unsafe_wrap(Array, c_ptr, 2*f.n_eq + f.n_ineq) + return as(Cdouble, call!(f, x, c)) 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. +# 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 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} +# 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] -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 + return _objfun_stack[i] 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) @@ -601,12 +761,27 @@ function _get_linear_constraints(Ab::LinearConstraints, n::Integer) return m, 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 (n,v,c) with n the number of constraints, v a vector to store the +# constraints on output (possibly a NullVector if the user is not interested in +# that), and c the callable object implementing the constraints. +_get_nonlinear_constraints(::Nothing, x::AbstractArray, str::AbstractString) = + 0, NullVector{Cdouble}(), nothing +_get_nonlinear_constraints(c::Tuple{Integer,Any}, x::AbstractArray, str::AbstractString) = + as(Int, c[1]), NullVector{Cdouble}(), c[2] +_get_nonlinear_constraints(c::Tuple{AbstractVector,Any}, x::AbstractArray, str::AbstractString) = + length(c[1]), c[1], c[2] +_get_nonlinear_constraints(c::Tuple{Any,Union{Integer,AbstractVector}}, x::AbstractArray, str::AbstractString) = + _get_nonlinear_constraints(reverse(c), x, str) +function _get_nonlinear_constraints(c::Any, x::AbstractArray, str::AbstractString) + # Assume c is callable. + try + v = c(x) + return length(v), NullVector{Cdouble}(), c + catch ex + ex isa MethodError || throw(ex) + throw(ArgumentError("method `c(x)` not implemented for non-linear $str constraints")) + end +end for (uplo, def) in ((:lower, typemin), (:upper, typemax)) diff --git a/test/runtests.jl b/test/runtests.jl index 38088c9..2bc0efe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,20 +40,11 @@ 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. @@ -132,7 +123,7 @@ optimizer(algo::Symbol) = println("\nCOBYLA:") # 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), + nonlinear_ineq = c_ineq, linear_ineq = (A_ineq, b_ineq), xl, xu, rhobeg = 1.0, rhoend = 1e-3, @@ -147,7 +138,7 @@ optimizer(algo::Symbol) = @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, + nonlinear_ineq = (c, c_ineq), linear_ineq = (A_ineq, b_ineq), xl, xu, rhobeg = 1.0, rhoend = 1e-3, From 5933fc01f246f65cc6ee58b9ebc429f1aa330128 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Fri, 20 Oct 2023 09:52:37 +0200 Subject: [PATCH 02/15] Update doc. --- README.md | 6 ++---- src/PRIMA.jl | 10 +++++----- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 7706ce6..496084d 100644 --- a/README.md +++ b/README.md @@ -190,12 +190,10 @@ Assuming `n = length(x)` is the number of variables, then: floating-point values to store `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/src/PRIMA.jl b/src/PRIMA.jl index aec31b3..6f035c1 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -112,7 +112,7 @@ const _doc_linear_constraints = """ 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 @@ -138,7 +138,7 @@ $(_doc_common_keywords) 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 @@ -165,7 +165,7 @@ $(_doc_npt) approximately solves the bound constrained optimization problem: - min f(x) s.t. x ∈ Ω ⊆ ℝⁿ + min f(x) subject to x ∈ Ω ⊆ ℝⁿ with @@ -199,7 +199,7 @@ $(_doc_bound_constraints) approximately solves the constrained optimization problem: - min f(x) s.t. x ∈ Ω ⊆ ℝⁿ + min f(x) subject to x ∈ Ω ⊆ ℝⁿ with @@ -234,7 +234,7 @@ $(_doc_nonlinear_constraints) approximately solves the constrained optimization problem: - min f(x) s.t. x ∈ Ω ⊆ ℝⁿ + min f(x) subject to x ∈ Ω ⊆ ℝⁿ with From 80ed3d904e7c137f1c4c7adecaaaf78dff06e827 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Fri, 20 Oct 2023 10:00:26 +0200 Subject: [PATCH 03/15] Mention BlackBoxOptim package --- README.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/README.md b/README.md index 496084d..1913aee 100644 --- a/README.md +++ b/README.md @@ -9,6 +9,13 @@ Powell](https://en.wikipedia.org/wiki/Michael_J._D._Powell) for minimizing a multi-variate objective function possibly under constraints and without derivatives. +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. + Formally, these algorithms are designed to solve problems of the form: ``` julia From a73f854a193217f0893a1f7893d42c868c3d7af7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Fri, 20 Oct 2023 10:33:26 +0200 Subject: [PATCH 04/15] Mention NOMAD package --- README.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/README.md b/README.md index 1913aee..082b734 100644 --- a/README.md +++ b/README.md @@ -16,6 +16,10 @@ objectives may be of interest: 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 blackbox + optimization problems. + Formally, these algorithms are designed to solve problems of the form: ``` julia From 762a1a1fabb67ec9aafe0c036111bc07bce73b6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Fri, 20 Oct 2023 18:58:38 +0200 Subject: [PATCH 05/15] Algorithms return a Status --- gen/wrapper.jl | 8 +++++--- src/wrappers.jl | 14 +++++++------- 2 files changed, 12 insertions(+), 10 deletions(-) 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/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 From f09c8dfa647c479dc0fe3da42dd68940fa2eda78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Fri, 20 Oct 2023 19:21:15 +0200 Subject: [PATCH 06/15] All algorithms yield (x,info) --- NEWS.md | 18 +++- Project.toml | 1 + README.md | 105 +++++++++--------- src/PRIMA.jl | 276 +++++++++++++++++++++++------------------------ test/runtests.jl | 213 +++++++++++++++++------------------- 5 files changed, 303 insertions(+), 310 deletions(-) diff --git a/NEWS.md b/NEWS.md index 245ecf0..48652b8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,10 +4,26 @@ - 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)` as for other algorithms. + - 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. + ## Version 0.1.1 - Keywords for other constraints than bounds have been renamed as diff --git a/Project.toml b/Project.toml index a3dbf56..6f39fc6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Éric Thiébaut and contributors"] version = "0.1.1" [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 082b734..9a9fbb6 100644 --- a/README.md +++ b/README.md @@ -82,62 +82,67 @@ 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...) +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...) ``` -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. +where `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). -For example, `rc` can be: +The objective function is called as `f(x)` with `x` the variables, it must +implement the following signature: -- `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 +f(x::Vector{Cdouble})::Real +``` -- `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. +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. -There are other possibilities which all indicate an error. Calling: +The output `info` has the following properties: ``` julia -PRIMA.reason(rc) +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 ``` -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: +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. @@ -169,8 +174,8 @@ 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`. - `maxfun` (default `100n`) is the maximum number of function evaluations allowed for the algorithm. If the number of calls to `f(x)` exceeds this @@ -183,22 +188,18 @@ 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 `n_eq` non-linear equality constraints defined by - `c_eq(x) = 0`. If the caller is interested in the values of `c_eq(x)` at the - returned solution, the keyword may be set with a 2-tuple `(v_eq, c_eq)` or - `(c_eq, v_eq)` with `v_eq` a vector of `n_eq` floating-point values to store - `c_eq(x)`. + `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 `n_ineq` non-linear inequality constraints defined by - `c_ineq(x) ≤ 0`. If the caller is interested in the values of `c_ineq(x)` at - the returned solution, the keyword may be set with a 2-tuple `(v_ineq, - c_ineq)` or `(c_ineq, v_ineq)` with `v_ineq` a vector of `n_ineq` - floating-point values to store `c_ineq(x)`. + `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 impose the linear equality constraints `Aₑ⋅x = bₑ`. diff --git a/src/PRIMA.jl b/src/PRIMA.jl index 6f035c1..68019d9 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -4,20 +4,65 @@ using PRIMA_jll const libprimac = PRIMA_jll.libprimac include("wrappers.jl") -export bobyqa, cobyla, lincoa, newuoa, uobyqa +export bobyqa, cobyla, lincoa, newuoa, uobyqa, issuccess using TypeUtils +using LinearAlgebra #------------------------------------------------------------------------------ # PUBLIC INTERFACE """ - PRIMA.reason(rc) -> str + PRIMA.Info + +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 sucessful. 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 + +LinearAlgebra.issuccess(info::Info) = issuccess(info.status) +LinearAlgebra.issuccess(status::Status) = -yields a textual message explaining `rc`, the code returned by one of the PRIMA -optimizers. +""" + PRIMA.reason(info::PRIMA.Info) -> str + PRIMA.reason(status::PRIMA.Status) -> str + +yield a textual message explaining `info.statu` 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)) # The high level wrappers. First the methods, then their documentation. @@ -26,28 +71,25 @@ for func in (:bobyqa, :newuoa, :uobyqa, :lincoa, :cobyla) @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 result of the algorithm 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): -const _doc_common_keywords = """ - `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 @@ -82,29 +124,22 @@ const _doc_bound_constraints = """ const _doc_nonlinear_constraints = """ - `nonlinear_eq` (default `nothing`) may be specified with a function, say - `c_eq`, implementing `n_eq` non-linear equality constraints defined by - `c_eq(x) = 0`. If the caller is interested in the values of `c_eq(x)` at the - returned solution, the keyword may be set with a 2-tuple `(v_eq, c_eq)` or - `(c_eq, v_eq)` with `v_eq` a vector of `n_eq` floating-point values to store - `c_eq(x)`. + `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 `n_ineq` non-linear inequality constraints defined by - `c_ineq(x) ≤ 0`. If the caller is interested in the values of `c_ineq(x)` at - the returned solution, the keyword may be set with a 2-tuple `(v_ineq, - c_ineq)` or `(c_ineq, v_ineq)` with `v_ineq` a vector of `n_ineq` - floating-point values to store `c_ineq(x)`. - + `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ᵢ`. """ """ @@ -120,21 +155,12 @@ 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: @@ -145,23 +171,14 @@ 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: @@ -177,16 +194,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) @@ -210,16 +218,9 @@ Approximations\") method. This algorithm is based on a trust region methood 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 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) -$(_doc_common_keywords) +$(_doc_npt) $(_doc_bound_constraints) @@ -245,16 +246,7 @@ This algorithm is based on a trust region methood 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) @@ -265,7 +257,7 @@ $(_doc_linear_constraints) """ lincoa """ - PRIMA.bobyqa!(f, x; kwds...) -> (fx, nf, rc) + 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 @@ -298,15 +290,15 @@ function bobyqa!(f, x::DenseVector{Cdouble}; 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) + status = prima_bobyqa(fp, n, x, fx, xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) + return Info(; fx = fx[], nf = nf[], status) finally _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 @@ -335,15 +327,15 @@ function newuoa!(f, x::DenseVector{Cdouble}; 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) + status = prima_newuoa(fp, n, x, fx, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) + return Info(; fx = fx[], nf = nf[], status) finally _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 @@ -370,15 +362,15 @@ function uobyqa!(f, x::DenseVector{Cdouble}; 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) + status = prima_uobyqa(fp, n, x, fx, nf, rhobeg, rhoend, ftarget, maxfun, iprint) + return Info(; fx = fx[], nf = nf[], status) finally _pop_objfun(fw) end end """ - PRIME.LinearConstraints + 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` @@ -388,7 +380,7 @@ 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 @@ -412,14 +404,14 @@ function cobyla!(f, x::DenseVector{Cdouble}; _check_rho(rhobeg, rhoend) xl = _get_lower_bound(xl, n) xu = _get_upper_bound(xu, n) - n_nl_eq, v_nl_eq, c_nl_eq = _get_nonlinear_constraints(nonlinear_eq, x, "equality") - n_nl_ineq, v_nl_ineq, c_nl_ineq = _get_nonlinear_constraints(nonlinear_ineq, x, "inequality") - n_lin_eq, A_eq, b_eq = _get_linear_constraints(linear_eq, n) - n_lin_ineq, A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) + 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) + A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) - # Total number of non-linear inequalities and vector to store them. - n_nl = 2*n_nl_eq + n_nl_ineq - v_nl = Vector{Cdouble}(undef, n_nl) + # Alocate 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 @@ -428,36 +420,37 @@ function cobyla!(f, x::DenseVector{Cdouble}; # Create wrapper to objective function and push it on top of per-thread # stack before calling optimizer. - fw = ObjFun(f, n, c_nl_eq, n_nl_eq, c_nl_ineq, n_nl_ineq) + fw = ObjFun(f, n, 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(n_nl, fp, n, x, fx, cstrv, v_nl, - n_lin_ineq, A_ineq, b_ineq, n_lin_eq, A_eq, b_eq, - xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, iprint) + 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) # Unpack constraints. - if length(v_nl_eq) > 0 - i = firstindex(v_nl) - for j in eachindex(v_nl_eq) - v_nl_eq[j] = v_nl[i] + 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(v_nl_ineq) > 0 - i = firstindex(v_nl) + 2*n_nl_eq - for j in eachindex(v_nl_ineq) - v_nl_ineq[j] = v_nl[i] + 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 (fx[], Int(nf[]), rc, cstrv[]) + return Info(; fx = fx[], nf = nf[], status, cstrv = cstrv[], nl_eq, nl_ineq) finally _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 @@ -481,8 +474,8 @@ function lincoa!(f, x::DenseVector{Cdouble}; _check_npt(npt, n) xl = _get_lower_bound(xl, n) xu = _get_upper_bound(xu, n) - n_lin_eq, A_eq, b_eq = _get_linear_constraints(linear_eq, n) - n_lin_ineq, A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) + A_eq, b_eq = _get_linear_constraints(linear_eq, n) + A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) # References for output values. cstrv = Ref{Cdouble}(NaN) # to store constraint violation @@ -495,10 +488,12 @@ function lincoa!(f, x::DenseVector{Cdouble}; fp = _push_objfun(lincoa, fw) # pointer to C-callable function try # Call low-level optimizer. - rc = prima_lincoa(fp, n, x, fx, cstrv, - n_lin_ineq, A_ineq, b_ineq, n_lin_eq, A_eq, b_eq, xl, xu, - nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint) - return (fx[], Int(nf[]), rc, cstrv[]) + 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) + return Info(; fx = fx[], nf = nf[], status, cstrv = cstrv[]) finally _pop_objfun(fw) end @@ -734,7 +729,7 @@ 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}() + NullMatrix{Cdouble}(), NullVector{Cdouble}() function _get_linear_constraints(Ab::LinearConstraints, n::Integer) A, b = Ab Base.has_offset_axes(A) && error( @@ -758,29 +753,26 @@ function _get_linear_constraints(Ab::LinearConstraints, n::Integer) end end b_ = _dense_array(T, b) - return m, A_, b_ + return A_, b_ end -# Yield (n,v,c) with n the number of constraints, v a vector to store the -# constraints on output (possibly a NullVector if the user is not interested in -# that), and c the callable object implementing the constraints. +# 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) = - 0, NullVector{Cdouble}(), nothing + Cdouble[], unconstrained _get_nonlinear_constraints(c::Tuple{Integer,Any}, x::AbstractArray, str::AbstractString) = - as(Int, c[1]), NullVector{Cdouble}(), c[2] -_get_nonlinear_constraints(c::Tuple{AbstractVector,Any}, x::AbstractArray, str::AbstractString) = - length(c[1]), c[1], c[2] -_get_nonlinear_constraints(c::Tuple{Any,Union{Integer,AbstractVector}}, 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) - # Assume c is callable. - try - v = c(x) - return length(v), NullVector{Cdouble}(), c - catch ex - ex isa MethodError || throw(ex) - throw(ArgumentError("method `c(x)` not implemented for non-linear $str constraints")) - end + # 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 for (uplo, def) in ((:lower, typemin), diff --git a/test/runtests.jl b/test/runtests.jl index 2bc0efe..082df58 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -47,9 +47,6 @@ optimizer(algo::Symbol) = 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) @@ -75,46 +72,46 @@ 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") + x, info = @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(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx + @test f(x) ≈ info.fx @test x0 == x0_sav 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") + x, info = @inferred PRIMA.uobyqa(f, x0; + rhobeg = 1.0, rhoend = 1e-3, + ftarget = -Inf, + maxfun = 200n, + iprint = PRIMA.MSG_EXIT) + msg = PRIMA.reason(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ fx + @test f(x) ≈ info.fx @test x0 == x0_sav 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") + x, info = @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(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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) end @@ -122,51 +119,51 @@ optimizer(algo::Symbol) = @testset "COBYLA" begin println("\nCOBYLA:") # First call with just the number of non-linear inequality constraints. - x, fx, nf, rc, cstrv = @inferred PRIMA.cobyla(f, x0; - nonlinear_ineq = c_ineq, - 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") + x, info = @inferred PRIMA.cobyla(f, x0; + nonlinear_ineq = c_ineq, + 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(info) + println("x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), c(x) = $(info.nl_ineq), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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) # 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, c_ineq), - 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; + nonlinear_ineq = (length(c_ineq(x0)), c_ineq), + 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(info) + println("x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), c(x) = $(info.nl_ineq), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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) 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") + x, info = @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(info) + println("x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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) end @@ -233,27 +230,22 @@ 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") + msg = PRIMA.reason(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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) @@ -262,24 +254,21 @@ 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") + msg = PRIMA.reason(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @test x ≈ [1.095247,1.2] rtol=0 atol=2e-2 - @test f(x) ≈ fx + @test f(x) ≈ info.fx end # Linearly constrained optimization. @@ -288,40 +277,36 @@ 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") + msg = PRIMA.reason(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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") + msg = PRIMA.reason(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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:") @@ -329,20 +314,18 @@ 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") + msg = PRIMA.reason(info) + println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") @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 From c20caf73782dfdea556fa8d091a65b05b1b3caa4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Fri, 20 Oct 2023 19:56:29 +0200 Subject: [PATCH 07/15] Update doc. --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 9a9fbb6..a6cd055 100644 --- a/README.md +++ b/README.md @@ -17,10 +17,11 @@ objectives may be of interest: 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 blackbox + Adaptive Direct Search algorithm* (MADS) designed for difficult black-box optimization problems. -Formally, these algorithms are designed to solve problems of the form: +Formally, the algorithms provided by `PRIMA` are designed to solve problems of +the form: ``` julia min f(x) subject to x ∈ Ω ⊆ ℝⁿ From 757249dc069368ada17de52a789750e9280f8f71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Sun, 22 Oct 2023 14:44:37 +0200 Subject: [PATCH 08/15] Use functions to have same keyword defaults --- src/PRIMA.jl | 86 +++++++++++++++++++++++++++------------------------- 1 file changed, 45 insertions(+), 41 deletions(-) diff --git a/src/PRIMA.jl b/src/PRIMA.jl index 68019d9..d0dba65 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -65,6 +65,22 @@ 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) +default_rhobeg() = 1.0 +default_rhoend(rhobeg::Real) = 1e-4*rhobeg + # The high level wrappers. First the methods, then their documentation. for func in (:bobyqa, :newuoa, :uobyqa, :lincoa, :cobyla) func! = Symbol(func, "!") @@ -220,8 +236,6 @@ objective function. No derivatives of the objective function are needed. $(_doc_common) -$(_doc_npt) - $(_doc_bound_constraints) $(_doc_linear_constraints) @@ -267,12 +281,12 @@ variables; on return, `x` is overwritten by an approximate solution. function bobyqa!(f, x::DenseVector{Cdouble}; 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, + 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 and get constraints. n = length(x) # number of variables _check_rho(rhobeg, rhoend) @@ -306,11 +320,11 @@ 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, + 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 @@ -343,10 +357,10 @@ 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, + 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 @@ -369,16 +383,6 @@ function uobyqa!(f, x::DenseVector{Cdouble}; end end -""" - 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}} - """ PRIMA.cobyla!(f, x; kwds...) -> info::PRIMA.Info @@ -388,17 +392,17 @@ variables; on return, `x` is overwritten by an approximate solution. """ function cobyla!(f, x::DenseVector{Cdouble}; - nonlinear_ineq = nothing, - nonlinear_eq = nothing, - linear_ineq::Union{LinearConstraints,Nothing} = nothing, - linear_eq::Union{LinearConstraints,Nothing} = nothing, + 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) @@ -458,16 +462,16 @@ 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, + 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) From 48177b1c82e446cd84937fc0b626e544216155ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Sun, 22 Oct 2023 14:56:37 +0200 Subject: [PATCH 09/15] Fix doc. --- src/PRIMA.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/PRIMA.jl b/src/PRIMA.jl index d0dba65..5035359 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -99,10 +99,10 @@ with `x` the variables, it must implement the following signature: f(x::Vector{Cdouble})::Real -The result of the algorithm 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. +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): @@ -159,7 +159,7 @@ const _doc_linear_constraints = """ """ """ - uobyqa(f, x0; kwds...) -> (x, fx, nf, rc) + uobyqa(f, x0; kwds...) -> x, info approximately solves the unconstrained optimization problem: @@ -219,7 +219,7 @@ $(_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: @@ -230,7 +230,7 @@ with Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, Aᵢ⋅x ≤ bᵢ, 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 methood +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. @@ -245,7 +245,7 @@ $(_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: @@ -256,7 +256,7 @@ with Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, and Aᵢ⋅x ≤ bᵢ } by M.J.D. Powell's LINCOA (for \"LINearly Constrained Optimization\") method. -This algorithm is based on a trust region methood where variables are updated +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. From 476aad1610c3078e76cedc65142cbd00071009a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Mon, 23 Oct 2023 09:11:57 +0200 Subject: [PATCH 10/15] General PRIMA.prima driver --- NEWS.md | 4 ++ README.md | 16 ++++-- src/PRIMA.jl | 81 ++++++++++++++++++++++++++- test/runtests.jl | 140 ++++++++++++++++++++++++----------------------- 4 files changed, 165 insertions(+), 76 deletions(-) diff --git a/NEWS.md b/NEWS.md index 48652b8..aeacd1f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,10 @@ ## Version 0.1.2 +- 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. diff --git a/README.md b/README.md index a6cd055..b8d8be7 100644 --- a/README.md +++ b/README.md @@ -79,10 +79,16 @@ 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 +x, info = prima(f, x0; kwds...) +``` + +To explicitly use one the specific Powell's algorithms, call one of: ``` julia -using PRIMA x, info = uobyqa(f, x0; kwds...) x, info = newuoa(f, x0; kwds...) x, info = bobyqa(f, x0; kwds...) @@ -90,9 +96,9 @@ x, info = cobyla(f, x0; kwds...) x, info = lincoa(f, x0; kwds...) ``` -where `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). +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). The objective function is called as `f(x)` with `x` the variables, it must implement the following signature: diff --git a/src/PRIMA.jl b/src/PRIMA.jl index 5035359..e521cb0 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -4,7 +4,7 @@ using PRIMA_jll const libprimac = PRIMA_jll.libprimac include("wrappers.jl") -export bobyqa, cobyla, lincoa, newuoa, uobyqa, issuccess +export bobyqa, cobyla, lincoa, newuoa, prima, uobyqa, issuccess using TypeUtils using LinearAlgebra @@ -51,6 +51,10 @@ struct Info 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) = @@ -82,7 +86,7 @@ default_rhobeg() = 1.0 default_rhoend(rhobeg::Real) = 1e-4*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...) @@ -158,6 +162,39 @@ const _doc_linear_constraints = """ impose the linear inequality constraints `Aᵢ⋅x ≤ bᵢ`. """ +""" + 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ᵢ, 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 @@ -270,6 +307,46 @@ $(_doc_linear_constraints) """ lincoa +""" + 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}; + 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; 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; 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; rhobeg, rhoend, iprint, ftarget, maxfun, npt, + xl, xu) + else + # NEWUOA is more efficient than UOBYQA for unconstrained problems. + return newuoa!(f, x; rhobeg, rhoend, iprint, ftarget, maxfun, npt) + end +end + """ PRIMA.bobyqa!(f, x; kwds...) -> info::PRIMA.Info diff --git a/test/runtests.jl b/test/runtests.jl index 082df58..c6c429c 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`") +prt_1(x::AbstractVector, info::PRIMA.Info) = prt_1(stdout, x, info) +function prt_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 + +prt_2(x::AbstractVector, info::PRIMA.Info) = prt_2(stdout, x, info) +function prt_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 + +prt_3(x::AbstractVector, info::PRIMA.Info) = prt_3(stdout, x, info) +function prt_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) @@ -72,28 +92,25 @@ optimizer(algo::Symbol) = @testset "NEWUOA" begin println("\nNEWUOA:") - x, info = @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(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.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...) + prt_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @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 end @testset "UOBYQA" begin println("\nUOBYQA:") - x, info = @inferred PRIMA.uobyqa(f, x0; - rhobeg = 1.0, rhoend = 1e-3, - ftarget = -Inf, - maxfun = 200n, - iprint = PRIMA.MSG_EXIT) - msg = PRIMA.reason(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") + kwds = (rhobeg = 1.0, rhoend = 1e-3, ftarget = -Inf, + maxfun = 200n, iprint = PRIMA.MSG_EXIT) + x, info = @inferred PRIMA.uobyqa(f, x0; kwds...) + prt_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @test f(x) ≈ info.fx @test x0 == x0_sav @@ -101,71 +118,61 @@ optimizer(algo::Symbol) = @testset "BOBYQA" begin println("\nBOBYQA:") - x, info = @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(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.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...) + prt_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @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 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, info = @inferred PRIMA.cobyla(f, x0; - nonlinear_ineq = c_ineq, - 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(info) - println("x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), c(x) = $(info.nl_ineq), status = $(info.status), msg = '$msg', evals = $(info.nf)") - @test x ≈ [3,2] atol=2e-2 rtol=0 - @test f(x) ≈ info.fx - @test x0 == x0_sav - @test check_bounds(xl, x, xu) - # Second call with an array to store the non-linear inequality constraints. - x, info = @inferred PRIMA.cobyla(f, x0; - nonlinear_ineq = (length(c_ineq(x0)), c_ineq), - 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(info) - println("x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), c(x) = $(info.nl_ineq), status = $(info.status), msg = '$msg', evals = $(info.nf)") + x, info = @inferred PRIMA.cobyla(f, x0; kwds..., + nonlinear_ineq = c_ineq) + prt_3(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @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 end @testset "LINCOA" begin println("\nLINCOA:") - x, info = @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(info) - println("x = $x, f(x) = $(info.fx), cstrv = $(info.cstrv), status = $(info.status), msg = '$msg', evals = $(info.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...) + prt_2(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @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 end end @@ -242,8 +249,7 @@ optimizer(algo::Symbol) = else continue end - msg = PRIMA.reason(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") + prt_1(x, info) @test x ≈ [1,1] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f(x) ≈ info.fx @@ -265,8 +271,7 @@ optimizer(algo::Symbol) = else continue end - msg = PRIMA.reason(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") + prt_1(x, info) @test x ≈ [1.095247,1.2] rtol=0 atol=2e-2 @test f(x) ≈ info.fx end @@ -285,8 +290,7 @@ optimizer(algo::Symbol) = else continue end - msg = PRIMA.reason(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") + prt_1(x, info) @test x ≈ [1.0,1.0] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f(x) ≈ info.fx @@ -303,8 +307,7 @@ optimizer(algo::Symbol) = else continue end - msg = PRIMA.reason(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") + prt_1(x, info) @test x ≈ [1.0,1.0] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f(x) ≈ info.fx @@ -322,8 +325,7 @@ optimizer(algo::Symbol) = else continue end - msg = PRIMA.reason(info) - println("x = $x, f(x) = $(info.fx), status = $(info.status), msg = '$msg', evals = $(info.nf)") + prt_1(x, info) @test x ≈ [1.441832,2.077557] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f2(x) ≈ info.fx end From a93c13564c1f10e8c83882109d071b119daeaea5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Mon, 23 Oct 2023 09:16:45 +0200 Subject: [PATCH 11/15] Fix doc. about non-linear constraints --- README.md | 6 ++++-- src/PRIMA.jl | 4 ++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index b8d8be7..63531d8 100644 --- a/README.md +++ b/README.md @@ -32,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` diff --git a/src/PRIMA.jl b/src/PRIMA.jl index e521cb0..4663f13 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -171,7 +171,7 @@ approximately solves the constrained optimization problem: 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 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 @@ -264,7 +264,7 @@ approximately solves the constrained optimization problem: 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 From 3fe03a50dc2cfffbdd46bfce1bc48f418e9dcfa6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Mon, 23 Oct 2023 11:31:58 +0200 Subject: [PATCH 12/15] Bump version --- NEWS.md | 2 +- Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index aeacd1f..c4be776 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,6 @@ # User visible changes in `PRIMA.jl` -## Version 0.1.2 +## 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) diff --git a/Project.toml b/Project.toml index 6f39fc6..d3c3131 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PRIMA" uuid = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed" authors = ["Éric Thiébaut and contributors"] -version = "0.1.1" +version = "0.1.2" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From 07b17f698751d60e8eef392c656ae0d3e4e9cbcb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Mon, 23 Oct 2023 14:27:53 +0200 Subject: [PATCH 13/15] Scaling of the variables --- NEWS.md | 9 ++ README.md | 20 ++++- src/PRIMA.jl | 219 +++++++++++++++++++++++++++++++++++------------ test/runtests.jl | 30 +++++++ 4 files changed, 222 insertions(+), 56 deletions(-) diff --git a/NEWS.md b/NEWS.md index c4be776..e6aeca9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -28,6 +28,15 @@ - `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 - Keywords for other constraints than bounds have been renamed as diff --git a/README.md b/README.md index 63531d8..497da4c 100644 --- a/README.md +++ b/README.md @@ -172,6 +172,23 @@ following table. 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 @@ -184,7 +201,8 @@ Assuming `n = length(x)` is the number of variables, then: - `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). - `maxfun` (default `100n`) is the maximum number of function evaluations allowed for the algorithm. If the number of calls to `f(x)` exceeds this diff --git a/src/PRIMA.jl b/src/PRIMA.jl index 4663f13..c52efd1 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -82,8 +82,10 @@ 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) -default_rhobeg() = 1.0 -default_rhoend(rhobeg::Real) = 1e-4*rhobeg +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, :cobyla, :lincoa, :newuoa, :prima, :uobyqa) @@ -110,10 +112,23 @@ successful and `x` is an approximate solution of the problem. Allowed keywords are (`n = length(x)` is the number of variables): -- `rhobeg` (default value `1.0`) is the initial radius of the trust region. +- `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. -- `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. +- `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 `$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 @@ -126,7 +141,8 @@ Allowed keywords are (`n = length(x)` is the number of variables): - `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 = """ @@ -316,7 +332,8 @@ variables; on return, `x` is overwritten by an approximate solution. """ function prima!(f, x::DenseVector{Cdouble}; - rhobeg::Real = default_rhobeg(), + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, rhoend::Real = default_rhoend(rhobeg), ftarget::Real = -Inf, maxfun::Integer = default_maxfun(x), @@ -330,20 +347,20 @@ function prima!(f, x::DenseVector{Cdouble}; nonlinear_eq = nothing) if nonlinear_eq !== nothing || nonlinear_ineq !== nothing # Only COBYLA can cope with non-linear constraints. - return cobyla!(f, x; rhobeg, rhoend, iprint, ftarget, maxfun, + 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; rhobeg, rhoend, iprint, ftarget, maxfun, npt, + 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; rhobeg, rhoend, iprint, ftarget, maxfun, npt, + 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; rhobeg, rhoend, iprint, ftarget, maxfun, npt) + return newuoa!(f, x; scale, rhobeg, rhoend, iprint, ftarget, maxfun, npt) end end @@ -356,20 +373,22 @@ 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 = 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) # 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 @@ -377,11 +396,13 @@ function bobyqa!(f, x::DenseVector{Cdouble}; # Create wrapper to objective function and push it on top of per-thread # stack before calling optimizer. - fw = ObjFun(f, n) # wrapper to objective function + fw = ObjFun(f, n, scl) # wrapper to objective function fp = _push_objfun(bobyqa, fw) # pointer to C-callable function try - # Call low-level optimizer. + # 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_objfun(fw) @@ -397,7 +418,8 @@ variables; on return, `x` is overwritten by an approximate solution. """ function newuoa!(f, x::DenseVector{Cdouble}; - rhobeg::Real = default_rhobeg(), + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, rhoend::Real = default_rhoend(rhobeg), ftarget::Real = -Inf, maxfun::Integer = default_maxfun(x), @@ -405,8 +427,9 @@ function newuoa!(f, x::DenseVector{Cdouble}; 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 @@ -414,11 +437,13 @@ function newuoa!(f, x::DenseVector{Cdouble}; # Create wrapper to objective function and push it on top of per-thread # stack before calling optimizer. - fw = ObjFun(f, n) # wrapper to objective function + fw = ObjFun(f, n, scl) # wrapper to objective function fp = _push_objfun(newuoa, fw) # pointer to C-callable function try - # Call low-level optimizer. + # 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_objfun(fw) @@ -434,7 +459,8 @@ variables; on return, `x` is overwritten by an approximate solution. """ function uobyqa!(f, x::DenseVector{Cdouble}; - rhobeg::Real = default_rhobeg(), + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, rhoend::Real = default_rhoend(rhobeg), ftarget::Real = -Inf, maxfun::Integer = default_maxfun(x), @@ -442,6 +468,7 @@ function uobyqa!(f, x::DenseVector{Cdouble}; # 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 @@ -449,11 +476,13 @@ function uobyqa!(f, x::DenseVector{Cdouble}; # Create wrapper to objective function and push it on top of per-thread # stack before calling optimizer. - fw = ObjFun(f, n) # wrapper to objective function + fw = ObjFun(f, n, scl) # wrapper to objective function fp = _push_objfun(uobyqa, fw) # pointer to C-callable function try - # Call low-level optimizer. + # 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_objfun(fw) @@ -469,7 +498,8 @@ variables; on return, `x` is overwritten by an approximate solution. """ function cobyla!(f, x::DenseVector{Cdouble}; - rhobeg::Real = default_rhobeg(), + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, rhoend::Real = default_rhoend(rhobeg), ftarget::Real = -Inf, maxfun::Integer = default_maxfun(x), @@ -483,12 +513,13 @@ function cobyla!(f, x::DenseVector{Cdouble}; # 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) + 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) - A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n) + A_eq, b_eq = _get_linear_constraints(linear_eq, n, scl) + A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n, scl) # Alocate vector to store all non-linear constraints (the non-linear # equalities being implemented as 2 inequalities each). @@ -501,14 +532,16 @@ function cobyla!(f, x::DenseVector{Cdouble}; # Create wrapper to objective function and push it on top of per-thread # stack before calling optimizer. - fw = ObjFun(f, n, c_nl_eq, length(nl_eq), c_nl_ineq, length(nl_ineq)) + 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. + # 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) @@ -539,7 +572,8 @@ variables; on return, `x` is overwritten by an approximate solution. """ function lincoa!(f, x::DenseVector{Cdouble}; - rhobeg::Real = default_rhobeg(), + scale::Union{AbstractVector{<:Real},Nothing} = default_scale, + rhobeg::Real = default_rhobeg, rhoend::Real = default_rhoend(rhobeg), ftarget::Real = -Inf, maxfun::Integer = default_maxfun(x), @@ -551,12 +585,13 @@ function lincoa!(f, x::DenseVector{Cdouble}; 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) - A_eq, b_eq = _get_linear_constraints(linear_eq, n) - 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 @@ -565,15 +600,17 @@ function lincoa!(f, x::DenseVector{Cdouble}; # Create wrapper to objective function and push it on top of per-thread # stack before calling optimizer. - fw = ObjFun(f, n) # wrapper to objective function + fw = ObjFun(f, n, scl) # wrapper to objective function fp = _push_objfun(lincoa, fw) # pointer to C-callable function try - # Call low-level optimizer. + # 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_objfun(fw) @@ -617,11 +654,17 @@ struct ObjFun{F,E,I} 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 + wrk::Vector{Cdouble} # workspace for variables end unconstrained(x::AbstractVector{T}) where {T} = NullVector{T}() -ObjFun(f::F, n::Integer) where {F} = ObjFun(f, n, unconstrained, 0, unconstrained, 0) +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 @@ -727,16 +770,38 @@ function _objfun(x_ptr::Ptr{Cdouble}, # (input) variables end function unsafe_call(f::ObjFun, x_ptr::Ptr{Cdouble}) - x = unsafe_wrap(Array, x_ptr, f.n) + 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 = unsafe_wrap(Array, x_ptr, f.n) + 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.wrk, 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 @@ -809,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) = +_get_linear_constraints(::Nothing, n::Integer, scl::AbstractVector) = NullMatrix{Cdouble}(), NullVector{Cdouble}() -function _get_linear_constraints(Ab::LinearConstraints, n::Integer) +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")) @@ -828,9 +893,20 @@ 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) @@ -856,17 +932,30 @@ function _get_nonlinear_constraints(c::Any, x::AbstractArray, str::AbstractStrin " 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 @@ -876,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/test/runtests.jl b/test/runtests.jl index c6c429c..388a968 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -71,6 +71,8 @@ end 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; @@ -103,6 +105,11 @@ end 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 @@ -114,6 +121,11 @@ end @test x ≈ [3,2] atol=2e-2 rtol=0 @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 @@ -131,6 +143,11 @@ end 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 @@ -156,6 +173,13 @@ end 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 @@ -173,6 +197,12 @@ end 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 From 58999bbf12c785707aea52f1b22fe3dbbe638f4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Mon, 23 Oct 2023 17:04:38 +0200 Subject: [PATCH 14/15] Fix typos detected by spell-checking GitHub Action --- .github/actions/spelling/allow.txt | 1 + src/PRIMA.jl | 16 +++++++-------- test/runtests.jl | 32 +++++++++++++++--------------- 3 files changed, 25 insertions(+), 24 deletions(-) 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/src/PRIMA.jl b/src/PRIMA.jl index cdb67a9..feae3ad 100644 --- a/src/PRIMA.jl +++ b/src/PRIMA.jl @@ -28,7 +28,7 @@ An object, say `info`, of this type has the following properties: info.nl_ineq # non-linear inequality constraints, empty vector if none Call `issuccess(info)` or `issuccess(info.status)` to check whether algorithm -was sucessful. Call `PRIMA.reason(info)` or `PRIMA.reason(info.status)` to +was successful. Call `PRIMA.reason(info)` or `PRIMA.reason(info.status)` to retrieve a textual description of the algorithm termination. """ @@ -62,7 +62,7 @@ LinearAlgebra.issuccess(status::Status) = PRIMA.reason(info::PRIMA.Info) -> str PRIMA.reason(status::PRIMA.Status) -> str -yield a textual message explaining `info.statu` or `status`, the status code +yield a textual message explaining `info.status` or `status`, the status code returned by one of the PRIMA optimizers. """ @@ -521,7 +521,7 @@ function cobyla!(f, x::DenseVector{Cdouble}; A_eq, b_eq = _get_linear_constraints(linear_eq, n, scl) A_ineq, b_ineq = _get_linear_constraints(linear_ineq, n, scl) - # Alocate vector to store all non-linear constraints (the non-linear + # 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)) @@ -655,7 +655,7 @@ struct ObjFun{F,E,I} 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 - wrk::Vector{Cdouble} # workspace for variables + x::Vector{Cdouble} # workspace for variables end unconstrained(x::AbstractVector{T}) where {T} = NullVector{T}() @@ -691,7 +691,7 @@ 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)` coresponding to the non-linear inequality constraints +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 @@ -757,8 +757,8 @@ function _objfun(x_ptr::Ptr{Cdouble}, # (input) variables return nothing end -# C-callable objective function for for problems with non-linear constraints -# (for COBYLA algorithm). +# 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 @@ -781,7 +781,7 @@ function unsafe_call(f::ObjFun, x_ptr::Ptr{Cdouble}, c_ptr::Ptr{Cdouble}) end function get_variables(f::ObjFun, x_ptr::Ptr{Cdouble}) - n, x, scl = f.n, f.wrk, f.scl + 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. diff --git a/test/runtests.jl b/test/runtests.jl index 388a968..65d0ff1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,20 +18,20 @@ optimizer(algo::Symbol) = algo === :prima ? PRIMA.prima : error("unknown optimizer `:$algo`") -prt_1(x::AbstractVector, info::PRIMA.Info) = prt_1(stdout, x, info) -function prt_1(io::IO, x::AbstractVector, info::PRIMA.Info) +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 -prt_2(x::AbstractVector, info::PRIMA.Info) = prt_2(stdout, x, info) -function prt_2(io::IO, x::AbstractVector, info::PRIMA.Info) +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 -prt_3(x::AbstractVector, info::PRIMA.Info) = prt_3(stdout, x, info) -function prt_3(io::IO, x::AbstractVector, info::PRIMA.Info) +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 @@ -97,7 +97,7 @@ end 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...) - prt_1(x, info) + print_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @test f(x) ≈ info.fx @test x0 == x0_sav @@ -117,7 +117,7 @@ end kwds = (rhobeg = 1.0, rhoend = 1e-3, ftarget = -Inf, maxfun = 200n, iprint = PRIMA.MSG_EXIT) x, info = @inferred PRIMA.uobyqa(f, x0; kwds...) - prt_1(x, info) + print_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @test f(x) ≈ info.fx @test x0 == x0_sav @@ -134,7 +134,7 @@ end 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...) - prt_1(x, info) + print_1(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @test f(x) ≈ info.fx @test x0 == x0_sav @@ -158,7 +158,7 @@ end # First call with just the number of non-linear inequality constraints. x, info = @inferred PRIMA.cobyla(f, x0; kwds..., nonlinear_ineq = c_ineq) - prt_3(x, info) + print_3(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @test f(x) ≈ info.fx @test x0 == x0_sav @@ -188,7 +188,7 @@ end 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...) - prt_2(x, info) + print_2(x, info) @test x ≈ [3,2] atol=2e-2 rtol=0 @test f(x) ≈ info.fx @test x0 == x0_sav @@ -279,7 +279,7 @@ end else continue end - prt_1(x, info) + print_1(x, info) @test x ≈ [1,1] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f(x) ≈ info.fx @@ -301,7 +301,7 @@ end else continue end - prt_1(x, info) + print_1(x, info) @test x ≈ [1.095247,1.2] rtol=0 atol=2e-2 @test f(x) ≈ info.fx end @@ -320,7 +320,7 @@ end else continue end - prt_1(x, info) + print_1(x, info) @test x ≈ [1.0,1.0] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f(x) ≈ info.fx @@ -337,7 +337,7 @@ end else continue end - prt_1(x, info) + print_1(x, info) @test x ≈ [1.0,1.0] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f(x) ≈ info.fx @@ -355,7 +355,7 @@ end else continue end - prt_1(x, info) + print_1(x, info) @test x ≈ [1.441832,2.077557] rtol=0 atol=(optim == :cobyla ? 3e-2 : 2e-2) @test f2(x) ≈ info.fx end From b18232ee77f03bbc4cf8aa03093076b486033733 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Thi=C3=A9baut?= Date: Mon, 23 Oct 2023 17:14:55 +0200 Subject: [PATCH 15/15] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d3c3131..98547d5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PRIMA" uuid = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed" authors = ["Éric Thiébaut and contributors"] -version = "0.1.2" +version = "0.2.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"