diff --git a/Project.toml b/Project.toml index 2b2deb79..8c3b93d9 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ version = "0.1.1" [deps] DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" @@ -11,6 +12,7 @@ TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" [compat] DynamicPolynomials = "0.3, 0.4" +ForwardDiff = "0.10" IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20" IntervalOptimisation = "0.4.1" Reexport = "0.2, 1" @@ -21,8 +23,8 @@ julia = "1.6" [extras] IntervalOptimisation = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "IntervalOptimisation", "SumOfSquares"] diff --git a/docs/src/lib/types.md b/docs/src/lib/types.md index 3fac1041..e1073e58 100644 --- a/docs/src/lib/types.md +++ b/docs/src/lib/types.md @@ -13,6 +13,7 @@ CurrentModule = RangeEnclosures ```@docs AbstractEnclosureAlgorithm +BranchAndBoundEnclosure NaturalEnclosure MooreSkelboeEnclosure SumOfSquaresEnclosure diff --git a/src/RangeEnclosures.jl b/src/RangeEnclosures.jl index 42e20de8..828b6316 100644 --- a/src/RangeEnclosures.jl +++ b/src/RangeEnclosures.jl @@ -1,6 +1,7 @@ module RangeEnclosures using DynamicPolynomials, Requires, Reexport +@reexport using ForwardDiff @reexport using IntervalArithmetic const Interval_or_IntervalBox = Union{Interval, IntervalBox} @@ -25,7 +26,8 @@ API include("enclose.jl") export enclose, relative_precision, - NaturalEnclosure, MooreSkelboeEnclosure, SumOfSquaresEnclosure, TaylorModelsEnclosure, + NaturalEnclosure, MooreSkelboeEnclosure, SumOfSquaresEnclosure, + TaylorModelsEnclosure, BranchAndBoundEnclosure, AbstractEnclosureAlgorithm end # module diff --git a/src/algorithms.jl b/src/algorithms.jl index 7d0d1ad4..1081dc04 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -114,3 +114,37 @@ Base.@kwdef struct SumOfSquaresEnclosure{T} <: AbstractEnclosureAlgorithm backend::T order::Int = 5 end + +""" + BranchAndBoundEnclosure + +Data type to bound the range of `f` over `X` using the branch and bound algorithm. + +### Fields + +- `maxdepth` (default `10`): maximum depth of the search tree +- `tol` (default `1e-3`): tolerance to compute the range of the function + +### Algorithm + +The algorithm evaluates a function `f` over an interval `X`. If the maximum depth is reached or the +width of `f(X)` is below the tolerance, the algorithm returns the computed range; otherwise it bisects +the interval/interval box. + +The algorithm also looks at the sign of the derivative / gradient to see if the range can be +computed directly. By default, the derivative / gradient is computed using `ForwardDiff.jl`, +but a custom value can be passed via the `df` keyword argument to [`enclose`](@ref). + +### Examples + +```jldoctest +julia> enclose(x -> -x^3/6 + 5x, 1..4, BranchAndBoundEnclosure()) +[4.83333, 10.5709] + +julia> enclose(x -> -x^3/6 + 5x, 1..4, BranchAndBoundEnclosure(tol=1e-2); df=x->-x^2/2+5) +[4.83333, 10.5709] +""" +Base.@kwdef struct BranchAndBoundEnclosure <: AbstractEnclosureAlgorithm + maxdepth = 10 + tol = 1e-3 +end diff --git a/src/branchandbound.jl b/src/branchandbound.jl index d63dd1b0..6eb93e29 100644 --- a/src/branchandbound.jl +++ b/src/branchandbound.jl @@ -1,101 +1,59 @@ # univariate -function enclose_BranchandBound(f::Function, dom::Interval; order=10, tol=0.6) - x0 = Interval(mid(dom)) - x = TaylorModel1(order, x0, dom) - return branchandbound(f(x-x0).pol, dom, tol) -end -# multivariate -function enclose_BranchandBound(f::Function, dom::IntervalBox{D}; - order=10, tol=0.6) where {D} - x0 = mid(dom) - set_variables(Float64, "x", order=2order, numvars=D) - x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i=1:D] - return branchandbound(f(x...).pol, dom - x0, tol) +@inline function enclose(f::Function, X::Interval, ba::BranchAndBoundEnclosure; + df=Base.Fix1(ForwardDiff.derivative, f)) + return _branch_bound(ba, f, X, df) end -@inline _Rnext(R::Vector{<:Interval}) = Interval(minimum(inf.(R)), maximum(sup.(R))) +@inline function enclose(f::Function, X::IntervalBox, ba::BranchAndBoundEnclosure; + df=t->ForwardDiff.gradient(w->f(w...), t.v)) + return _branch_bound(ba, f, X, df) +end -const Kmax = 1000 +function _branch_bound(ba::BranchAndBoundEnclosure, f::Function, X::Interval_or_IntervalBox, df; + initial=emptyinterval(first(X)), + cnt=1) -function branchandbound(p::Union{Taylor1, TaylorN}, - dom::Interval_or_IntervalBox, - tol::Number) - K = 1 - Rprev = evaluate(p, dom) - D = bisect(dom) - R = [evaluate(p, D[1]), evaluate(p, D[2])] - Rnext = R[1] ∩ R[2] + dfX = df(X) + range_extrema, flag = _monotonicity_check(f, X, dfX) + flag && return hull(range_extrema, initial) - while (Rprev.hi - Rnext.hi) >= tol * diam(Rnext) && - (Rprev.lo - Rnext.lo) >= tol * diam(Rnext) + fX = f(X...) # TODO: allow user to choose how to evaluate this (mean value, natural enclosure) + # if tolerance or maximum number of iteration is met, return current enclosure + if diam(fX) <= ba.tol || cnt == ba.maxdepth + return hull(fX, initial) + end - Rprev = _Rnext(R) - R_x = sup.(R) - R_n = inf.(R) - max_range = maximum(R_x) - max_index = findall(x->x == max_range, R_x)[1] - min_range = minimum(R_n) - min_index = findall(x->x == min_range, R_n)[1] - l_D = length(D[1]) - K += 1 - if K == Kmax - error("convergence not achieved in branch and bound after $Kmax steps") - end + X1, X2 = bisect(X) + y1 = _branch_bound(ba, f, X1, df; initial=initial, cnt=cnt+1) + return _branch_bound(ba, f, X2, df; initial=y1, cnt=cnt+1) +end - D, R = divide_dom!(p, D, R, max_index) - if max_index < min_index - min_index += 1 - D, R = divide_dom!(p, D, R, min_index) - end - Rnext = _Rnext(R) +function _monotonicity_check(f::Function, X::Interval, dfX::Interval) + if inf(dfX) >= 0 || sup(dfX) <= 0 # monotone function, evaluate at extrema + lo = Interval(X.lo) + hi = Interval(X.hi) + return hull(f(lo), f(hi)), true end - return _Rnext(R) -end -function divide_dom!(p::Union{TaylorN{T}, Taylor1{T}}, - D::Union{Array{IntervalBox{N, W}, M}, Array{Interval{W}, M}}, - R::Array{Interval{W}, M}, index::Number) where {N, T, M, W} - BA = [ ((D[index][i]).hi - (D[index][i]).lo) for i = 1:length(D[1])] - Beta1 = maximum(BA) - β = findall(x->x == Beta1, BA)[1] - D1, D2 = bisect(D[index], β) - D[index] = D1 - DD = push!(D[1:index], D2) - D = append!(DD, D[(index + 1):length(D)]) - R[index] = evaluate(p, D[index]) - RR = append!(R[1:index], evaluate(p, D[index + 1])) - R = append!(RR, R[(index + 1):length(R)]) - return D, R + return zero(eltype(dfX)), false end -function enclose_binary(f, dom::Interval; kmax=3, tol=1e-3, algorithm=:IntervalArithmetic) - y = enclose(f, dom, algorithm) - yinf, ysup = inf(y), sup(y) - kmax == 0 && return Interval(yinf, ysup) - x = bisect(dom) - fx1 = enclose(f, x[1], algorithm) - fx2 = enclose(f, x[2], algorithm) - ynew = hull(fx1, fx2) - ynew_inf, ynew_sup = inf(ynew), sup(ynew) - inf_close = abs(yinf - ynew_inf) <= tol - sup_close = abs(ysup - ynew_sup) <= tol - both_close = inf_close && sup_close - inf_improves = ynew_inf > yinf - sup_improves = ynew_sup < ysup - both_improve = inf_improves && sup_improves - if both_close || !both_improve - return Interval(yinf, ysup) - end - yinf = max(yinf, ynew_inf) - ysup = min(ysup, ynew_sup) - if inf_improves - yinf = ynew_inf - end - if sup_improves - ysup = ynew_sup +function _monotonicity_check(f::Function, X::IntervalBox{N}, ∇fX::AbstractVector) where {N} + low = zeros(eltype(X), N) + high = zeros(eltype(X), N) + + @inbounds for (i, di) in enumerate(∇fX) + if di >=0 # increasing + high[i] = sup(X[i]) + low[i] = inf(X[i]) + elseif di <= 0 # decreasing + high[i] = inf(X[i]) + low[i] = sup(X[i]) + else + return di, false + end end - e1 = enclose_binary(f, x[1], kmax=kmax-1, algorithm=algorithm) - e2 = enclose_binary(f, x[2], kmax=kmax-1, algorithm=algorithm) - return Interval(hull(e1, e2)) + + return hull(f(low...), f(high...)), true end diff --git a/src/enclose.jl b/src/enclose.jl index 9ec3d4c7..272d39c8 100644 --- a/src/enclose.jl +++ b/src/enclose.jl @@ -43,22 +43,6 @@ function enclose(f::Function, dom::Interval_or_IntervalBox, return _enclose(solver, f, dom; kwargs...) end -function enclose(f::Function, dom::Interval_or_IntervalBox, solver::Symbol; kwargs...) - - 𝑂 = Dict(kwargs) - - if solver == :BranchandBound - tol = :tol ∈ keys(𝑂) ? 𝑂[:tol] : 0.6 - order = :order ∈ keys(𝑂) ? 𝑂[:order] : 10 - #solver - R = enclose_BranchandBound(f, dom, order=order, tol=tol) - else - error("algorithm $algorithm unknown") - end - - return R -end - function enclose(p::AbstractPolynomialLike, dom::Interval_or_IntervalBox, solver::AbstractEnclosureAlgorithm=NaturalEnclosure(); kwargs...) f(x...) = p(variables(p) => x) diff --git a/test/runtests.jl b/test/runtests.jl index 91cfc560..a66f971c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,8 @@ using IntervalOptimisation, Test, RangeEnclosures using DynamicPolynomials: @polyvar -available_solvers = (NaturalEnclosure(), TaylorModelsEnclosure(), :BranchandBound, MooreSkelboeEnclosure()) +available_solvers = (NaturalEnclosure(), TaylorModelsEnclosure(), + BranchAndBoundEnclosure(), MooreSkelboeEnclosure()) include("univariate.jl") include("multivariate.jl")