Determining the sign of a function over a given interval

Specific Domains Optimization (Mathematical)
11

Just wanted to build on @mforets answer specifically for the original question in case someone is interested in a similar problem in the future,

using BranchAndPrune, RangeEnclosures, IntervalArithmetic
using BenchmarkTools, Plots

# Define the search type as mutable so we can change the value of
# flags in order to terminate the search.
# Breadth first search was more applicable in this search problem.
mutable struct BoundSearch{T} <: AbstractBreadthFirstSearch{Interval{T}}
    f::Function
    initial::Interval
    algorithm::Symbol
    tol::Real
    found_positive::Bool
    found_negative::Bool
    breached_tolerance::Bool
    function BoundSearch{T}(f::Function, initial::Interval, algorithm::Symbol, tol::Real) where {T}

        if tol < 0
            throw(ArgumentError("Tolerance must be a positive number"))
        end
        new{T}(f,initial,algorithm,tol,false,false,false)
    end
end

function BranchAndPrune.process(search::BoundSearch, interval::Interval)

    f_range = enclose(search.f, interval, search.algorithm)

    if (search.found_positive && search.found_negative) || (search.breached_tolerance)
        # No further processing required in this case
        return :discard, interval
    elseif inf(f_range) > 0.0
        # We have found a sub-interval where the function is
        # uniformly positive
        search.found_positive = true
        return :store, interval
    elseif sup(f_range) < 0.0
        # We have found a sub-interval where the function is
        # uniformly negative
        search.found_negative = true
        return :store, interval
    elseif diam(interval) < search.tol
        search.breached_tolerance = true
        return :discard, interval
    else
        return :bisect, interval
    end
end

BranchAndPrune.bisect(::BoundSearch, interval) = bisect(interval,0.5)

function run_search(f, interval, algorithm, tol)

    search = BoundSearch{Float64}(f, interval, algorithm, tol)
    local endtree = nothing
    for working_tree in search
        endtree = working_tree
    end
    return endtree, search
end

"""
    sign(f::Function, int::Interval)
return
- `+1` if `f` is uniformly positive on `int`
- `-1` if `f` is uniformly negative on `int`
- `0` if `f` has at least one zero crossing in `int` (f assumed continuous)
"""
function Base.sign(f::Function, int::Interval; algorithm = :IntervalArithmetic, tol = 1e-2)
    tree, search = run_search(f,int,algorithm,tol)
    if search.found_positive && search.found_negative
        return 0
    elseif search.breached_tolerance
        error("Search tolerance not tight enough")
    elseif search.found_positive
        return 1
    else search.found_negative
        return -1
    end
end

The major change is to declare our search type as mutable so we can set the value of certain flags based on the search results – this allows us to terminate the search once certain criteria are fulfilled.
Further, we subtype AbstractBreadthFirstSearch: the status of an interval around a root of f is undetermined because the bounds will always be [+ve number, -ve number]. Thus, if there is a root, depth first search will always exceed the given tolerance which is undesirable.

f(x) = (x - 0.0)*(x - 0.1)*(x - 1.0)
int = 0 .. 1
s = @btime sign(f,int) # 32.878 μs (415 allocations: 27.25 KiB)

example_poly

2 Likes