From 544b63531475c9fe85c6ba702cf1d4e3dfce9b32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Sun, 19 Feb 2023 13:01:05 +0100 Subject: [PATCH] Use branch and prune as backend for the search (#187) * Use BranchAndPrune --- Project.toml | 2 + docs/src/internals.md | 6 +- examples/root_search_iterator.jl | 80 -------- src/IntervalRootFinding.jl | 2 +- src/branch_and_bound.jl | 326 ------------------------------- src/roots.jl | 144 ++++++-------- test/branch_and_bound.jl | 134 ------------- test/runtests.jl | 2 - test/search_interface.jl | 21 -- 9 files changed, 61 insertions(+), 656 deletions(-) delete mode 100644 examples/root_search_iterator.jl delete mode 100644 src/branch_and_bound.jl delete mode 100644 test/branch_and_bound.jl delete mode 100644 test/search_interface.jl diff --git a/Project.toml b/Project.toml index 7c873b27..3e8cbcc4 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,7 @@ uuid = "d2bf35a9-74e0-55ec-b149-d360ff49b807" version = "0.5.11" [deps] +BranchAndPrune = "d3bc4f2e-91e6-11e9-365e-cd067da536ce" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -11,6 +12,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] +BranchAndPrune = "0.2.1" ForwardDiff = "0.10" IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20" Polynomials = "0.5, 0.6, 0.7, 0.8, 1, 2, 3" diff --git a/docs/src/internals.md b/docs/src/internals.md index 95c8ef37..3a583dc9 100644 --- a/docs/src/internals.md +++ b/docs/src/internals.md @@ -47,7 +47,7 @@ Contractors play a central role in the algorithm: they are the only part of it t ## Search object -Now that we have presented the foundation of the internal algorithm, we can discuss the representation of the search. Each search strategy has a type associated, the defined types being `BreadthFirstSearch` and `DepthFirstSearch`. +Now that we have presented the foundation of the internal algorithm, we can discuss the representation of the search. Each search strategy has a type associated, the defined types being `BreadthFirst` and `DepthFirst`. A search must be given three pieces of information: 1. The region to search; @@ -64,8 +64,8 @@ df (generic function with 1 method) julia> C = Newton(f, df) Newton{typeof(f),typeof(df)}(f, df) -julia> search = DepthFirstSearch(-10..10, C, 1e-10) -DepthFirstSearch{Interval{Float64},Newton{typeof(f),typeof(df)}, +julia> search = DepthFirst(-10..10, C, 1e-10) +DepthFirst{Interval{Float64},Newton{typeof(f),typeof(df)}, Float64}(Root([-10, 10], :unknown), Newton{typeof(f),typeof(df)}(f, df), 1.0e-10) ``` diff --git a/examples/root_search_iterator.jl b/examples/root_search_iterator.jl deleted file mode 100644 index 7dbf0f71..00000000 --- a/examples/root_search_iterator.jl +++ /dev/null @@ -1,80 +0,0 @@ -using IntervalArithmetic -using IntervalRootFinding - -## -# Available strategies -## - -f(x) = sin(x) -contractor = Newton(f, cos) - -@info "Breadth first search." -# A search takes a starting region, a contractor and a tolerance as argument -search = BreadthFirstSearch(-10..10, contractor, 1e-10) - -# This is needed to avoid the data being discarded at the end of the loop use -# `local tree` instead inside a function. -global endtree -for (k, tree) in enumerate(search) - println("Tree at iteration $k") - println(tree) # The tree has custom printing - global endtree = tree -end - -rts = data(endtree) # Use `data` to only get the leaves marked as `:final` -println("Final $(length(rts)) roots: $rts") -println() - -@info "Depth first search." -# The second available type of root search is depth first -search = DepthFirstSearch(-10..10, contractor, 1e-10) -for tree in search - global endtree = tree -end # Go to the end of the iteration - -# Since the index of each node is printed and indices are attributed in the -# order nodes are processed, comparing this tree with the last iteration of the -# preceding show the difference between breadth first and depth first searches. -println(endtree) - - -## -# Custom strategy -## - -# Define custom BBSearch that store empty intervals rather than discarding it -# Functions of the interface must be explicitely imported -import IntervalRootFinding - -struct MySearch{R <: Region, C <: Contractor, T <: Real} <: BreadthFirstBBSearch{Root{R}} - initial::Root{R} - contractor::C - tol::T -end - -# This must be implemented, other functions needed are already implemented for -# BBSearch using Root as data (bisect function) or the general fallback match -# this case (root_element function) -function IntervalRootFinding.process(search::MySearch, r::Root) - contracted_root = search.contractor(r, search.tol) - status = root_status(contracted_root) - - # We only store the unrefined intervals so the final intervals cover all - # the initial region - unrefined_root = Root(r.interval, status) - status == :unique && return :store, unrefined_root - status == :empty && return :store, unrefined_root # Store largest known empty intervals - status == :unknown && diam(contracted_root) < search.tol && return :store, unrefined_root - return :bisect, unrefined_root # Always bisect the original interval to bypass [NaN, NaN] -end - -@info "Search with no interval discarded." - -search = MySearch(Root(-10..10, :unknown), contractor, 1e-10) -for tree in search - global endtree = tree -end # Go to the end of the iteration - -# In these tree we can see that the union of the resulting intervals (including -# those containing no solution) covers the starting interval. -println(endtree) diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index cce68ffe..fe414297 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -2,6 +2,7 @@ module IntervalRootFinding +using BranchAndPrune using IntervalArithmetic using ForwardDiff using StaticArrays @@ -47,7 +48,6 @@ include("krawczyk.jl") include("complex.jl") include("contractors.jl") -include("branch_and_bound.jl") include("roots.jl") include("newton1d.jl") include("quadratic.jl") diff --git a/src/branch_and_bound.jl b/src/branch_and_bound.jl deleted file mode 100644 index 39de0fcc..00000000 --- a/src/branch_and_bound.jl +++ /dev/null @@ -1,326 +0,0 @@ -import Base: copy, eltype, iterate, IteratorSize -import Base: getindex, setindex!, delete! - -export BBSearch, BreadthFirstBBSearch, DepthFirstBBSearch -export data -export copy, eltype, iterate, IteratorSize - -abstract type AbstractBBNode end - -""" - BBNode <: AbstractBBNode - -Intermediate node of a `BBTree`. Does not contain any data by itself, -only redirect toward its children. -""" -struct BBNode <: AbstractBBNode - parent::Int - children::Vector{Int} -end - -""" - BBLeaf{DATA} <: AbstractBBLeaf - -Leaf node of a `BBTree` that contains some data. Its status is either - - `:working`: the leaf will be further processed. - - `:final`: the leaf won't be touched anymore. -""" -struct BBLeaf{DATA} <: AbstractBBNode - data::DATA - parent::Int - status::Symbol -end - -function BBLeaf(data::DATA, parent::Int) where {DATA} - BBLeaf{DATA}(data, parent, :working) -end - -function BBNode(leaf::BBLeaf, child1::Int, child2::Int) - BBNode(leaf.parent, Int[child1, child2]) -end - -""" - BBTree{DATA} - -Tree storing the data used and produced by a branch and bound search in a -structured way. - -Nodes and leaves can be accessed using their index using the bracket syntax -`wt[node_id]`. However this is slow, as nodes and leaves are stored separately. - -Support the iterator interface. The element yielded by the iteration are -tuples `(node_id, lvl)` where `lvl` is the depth of the node in the tree. -""" -struct BBTree{DATA} - nodes::Dict{Int, BBNode} - leaves::Dict{Int, BBLeaf{DATA}} - working_leaves::Vector{Int} -end - -function BBTree(rootdata::DATA) where {DATA} - rootleaf = BBLeaf(rootdata, 0) - BBTree{DATA}(Dict{Int, BBNode}(), Dict(1 => rootleaf), Int[1]) -end - -show(io::IO, wn::BBNode) = print(io, "Node with children $(wn.children)") - -function show(io::IO, wl::BBLeaf) - print(io, "Leaf (:$(wl.status)) with data $(wl.data)") -end - -function show(io::IO, wt::BBTree{DATA}) where {DATA} - println(io, "Working tree with $(nnodes(wt)) elements of type $DATA") - - if nnodes(wt) > 0 - println(io, "Indices: ", vcat(collect(keys(wt.nodes)), collect(keys(wt.leaves))) |> sort) - println(io, "Structure:") - for (id, lvl) in wt - println(io, " "^lvl * "[$id] $(wt[id])") - end - end -end - -# Root node has id 1 and parent id 0 -root(wt::BBTree) = wt[1] -is_root(wt::BBTree, id::Int) = (id == 1) - -""" - nnodes(wt::BBTree) - -Number of nodes (including leaves) in a `BBTree`. -""" -nnodes(wt::BBTree) = length(wt.nodes) + length(wt.leaves) - -""" - data(leaf::BBLeaf) - -Return the data stored in the leaf. -""" -data(leaf::BBLeaf) = leaf.data - -""" - data(wt::BBTree) - -Return all the data stored in a `BBTree` as a list. The ordering of the elements -is arbitrary. -""" -data(wt::BBTree) = data.(values(wt.leaves)) - -function newid(wt::BBTree) - k1 = keys(wt.nodes) - k2 = keys(wt.leaves) - - if length(k1) > 0 - m1 = maximum(k1) - else - m1 = 0 - end - - if length(k2) > 0 - m2 = maximum(k2) - else - m2 = 0 - end - - return max(m1, m2) + 1 -end - -# Index operations (slower than manipulating the node directly in the correct -# dictionary) -function getindex(wt::BBTree, id) - haskey(wt.nodes, id) && return wt.nodes[id] - haskey(wt.leaves, id) && return wt.leaves[id] - error("getindex failed: no index $id") # TODO: make better error -end - -setindex!(wt::BBTree, val::BBNode, id) = setindex!(wt.nodes, val, id) -setindex!(wt::BBTree, val::BBLeaf, id) = setindex!(wt.leaves, val, id) - -function delete!(wt::BBTree, id) - if haskey(wt.nodes, id) - delete!(wt.nodes, id) - elseif haskey(wt.leaves, id) - delete!(wt.leaves, id) - else - error("delete! failed: no index $id") # TODO: make better error - end -end - -""" - discard_leaf!(wt::BBTree, id::Int) - -Delete the `BBLeaf` with index `id` and all its ancestors to which it is -the last descendant. -""" -function discard_leaf!(wt::BBTree, id::Int) - leaf = wt.leaves[id] - delete!(wt.leaves, id) - recursively_delete_parent!(wt, leaf.parent, id) -end - -function recursively_delete_parent!(wt, id_parent, id_child) - if !is_root(wt, id_child) - parent = wt.nodes[id_parent] - siblings = parent.children - if length(parent.children) == 1 # The child has no siblings, so delete the parent - delete!(wt.nodes, id_parent) - recursively_delete_parent!(wt, parent.parent, id_parent) - else # The child has siblings so remove it from the children list - deleteat!(parent.children, searchsortedfirst(parent.children, id_child)) - end - end -end - -function iterate(wt::BBTree, (id, lvl)=(0, 0)) - id, lvl = next_id(wt, id, lvl) - lvl == 0 && return nothing - return (id, lvl), (id, lvl) -end - -function next_id(wt::BBTree, id, lvl) - lvl == 0 && return (1, 1) - node = wt[id] - isa(node, BBNode) && return (node.children[1], lvl + 1) - return next_sibling(wt, id, lvl) -end - -function next_sibling(wt::BBTree, sibling, lvl) - parent = wt[sibling].parent - parent == 0 && return (0, 0) - children = wt[parent].children - maximum(children) == sibling && return next_sibling(wt, parent, lvl - 1) - id = minimum(filter(x -> x > sibling, children)) - return (id, lvl) -end - - -""" - BBSearch{DATA} - -Branch and bound search interface in element of type DATA. - -This interface provide an iterable that perform the search. - -There is currently three types of search supported `BreadFirstBBSearch`, -`DepthFirstBBSearch` and `KeyBBSearch`, each one processing the element of the -tree in a different order. When subtyping one of these, the following methods -must be implemented: - - `root_element(::BBSearch)`: return the element with which the search is started - - `process(::BBSearch, elem::DATA)`: return a symbol representing the action - to perform with the element `elem` and an object of type `DATA` reprensenting - the state of the element after processing (may return `elem` unchanged). - - `bisect(::BBSearch, elem::DATA)`: return two elements of type `DATA` build - by bisecting `elem` - -Subtyping `BBSearch` directly allows to have control over the order in which -the elements are process. To do this the following methods must be implemented: - - `root_element(::BBSearch)`: return the first element to be processed. Use - to build the initial tree. - - `get_leaf_id!(::BBSearch, wt::BBTree)`: return the id of the next leaf that - will be processed and remove it from the list of working leaves of `wt`. - - `insert_leaf!(::BBSearch, wt::BBTree, leaf::BBLeaf)`: insert a leaf in the - list of working leaves. - -# Valid symbols returned by the process function - - `:store`: the element is considered as final and is stored, it will not be - further processed - - `:bisect`: the element is bisected and each of the two resulting part will - be processed - - `:discard`: the element is discarded from the tree, allowing to free memory -""" -abstract type BBSearch{DATA} end - -abstract type BreadthFirstBBSearch{DATA} <: BBSearch{DATA} end -abstract type DepthFirstBBSearch{DATA} <: BBSearch{DATA} end - -""" - KeyBBSearch{DATA} <: BBSearch{DATA} - -Interface to a branch and bound search that use a key function to decide which -element to process first. The search process first the element with the largest -key as computed by `keyfunc(ks::KeyBBSearch, elem)`. - -!!! warning - Untested. -""" -abstract type KeyBBSearch{DATA} <: BBSearch{DATA} end - -""" - root_element(search::BBSearch) - -Return the initial element of the search. The `BBTree` will be build around it. - -Can be define for custom searches that are direct subtype of `BBSearch`, default -behavior is to fetch the field `initial` of the search. -""" -root_element(search::BBSearch) = search.initial - -""" - get_leaf_id!(::BBSearch, wt::BBTree) - -Return the id of the next leaf that will be processed and remove it from the -list of working leaves. - -Must be define for custom searches that are direct subtype of `BBSearch`. -""" -get_leaf_id!(::BreadthFirstBBSearch, wt::BBTree) = popfirst!(wt.working_leaves) -get_leaf_id!(::DepthFirstBBSearch, wt::BBTree) = pop!(wt.working_leaves) -get_leaf_id!(::KeyBBSearch, wt::BBTree) = popfirst!(wt.working_leaves) - -""" - insert_leaf!(::BBSearch, wt::BBTree, leaf::BBLeaf) - -Insert the id of a new leaf that has been produced by bisecting an older leaf -into the list of working leaves. - -Must be define for custom searches that are direct subtype of `BBSearch`. -""" -function insert_leaf!(::Union{BreadthFirstBBSearch{DATA}, DepthFirstBBSearch{DATA}}, - wt::BBTree{DATA}, leaf::BBLeaf{DATA}) where {DATA} - id = newid(wt) - wt.leaves[id] = leaf - push!(wt.working_leaves, id) - return id -end - -function insert_leaf!(::KS, wt::BBTree{DATA}, leaf::BBLeaf{DATA}) where {DATA, KS <: KeyBBSearch{DATA}} - id = newid(wt) - wt.leaves[id] = leaf - keys = keyfunc.(KS, wt.working_leaves) - current = keyfunc(KS, leaf) - - # Keep the working_leaves sorted - insert!(wt.working_leaves, searchsortedfirst(keys, current), id) - return id -end - -eltype(::Type{BBS}) where {DATA, BBS <: BBSearch{DATA}} = BBTree{DATA} -IteratorSize(::Type{BBS}) where {BBS <: BBSearch} = Base.SizeUnknown() - -function iterate(search::BBSearch{DATA}, - wt::BBTree=BBTree(root_element(search))) where {DATA} - - isempty(wt.working_leaves) && return nothing - - id = get_leaf_id!(search, wt) - X = wt.leaves[id] - action, newdata = process(search, data(X)) - if action == :store - wt.leaves[id] = BBLeaf(newdata, X.parent, :final) - elseif action == :bisect - child1, child2 = bisect(search, newdata) - leaf1 = BBLeaf(child1, id, :working) - leaf2 = BBLeaf(child2, id, :working) - id1 = insert_leaf!(search, wt, leaf1) - id2 = insert_leaf!(search, wt, leaf2) - wt.nodes[id] = BBNode(X, id1, id2) - delete!(wt.leaves, id) - elseif action == :discard - discard_leaf!(wt, id) - else - error("Branch and bound: process function of the search object return " * - "unknown action: $action for element $X. Valid actions are " * - ":store, :bisect and :discard.") - end - return wt, wt -end diff --git a/src/roots.jl b/src/roots.jl index 0f62a610..c5120e37 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -2,76 +2,40 @@ import IntervalArithmetic: diam, isinterior, bisect, isnan export branch_and_prune, Bisection, Newton -export BreadthFirstSearch, DepthFirstSearch diam(r::Root) = diam(interval(r)) isnan(X::IntervalBox) = any(isnan.(X)) isnan(r::Root) = isnan(interval(r)) -""" - BreadthFirstSearch <: BreadthFirstBBSearch - -Type implementing the `BreadthFirstBBSearch` interface for interval roots finding. - -# Fields: - - `initial`: region (as a `Root` object) in which roots are searched. - - `contractor`: contractor to use (`Bisection`, `Newton` or `Krawczyk`) - - `tol`: tolerance of the search -""" -struct BreadthFirstSearch{R <: Region, C <: Contractor, T <: Real} <: BreadthFirstBBSearch{Root{R}} - initial::Root{R} - contractor::C - tol::T -end - -""" - DepthFirstSearch <: DepthFirstBBSearch - -Type implementing the `DepthFirstBBSearch` interface for interval roots finding. - -# Fields: - - `initial`: region (as a `Root` object) in which roots are searched. - - `contractor`: contractor to use (`Bisection`, `Newton` or `Krawczyk`) - - `tol`: tolerance of the search -""" -struct DepthFirstSearch{R <: Region, C <: Contractor, T <: Real} <: DepthFirstBBSearch{Root{R}} - initial::Root{R} - contractor::C - tol::T +struct RootProblem{T} + abstol::T end -BreadthFirstSearch(X::Region, C::Contractor, tol::Real) = BreadthFirstSearch(Root(X, :unknown), C, tol) -DepthFirstSearch(X::Region, C::Contractor, tol::Real) = DepthFirstSearch(Root(X, :unknown), C, tol) - -root_element(search::BBSearch{Root{R}}) where {R <: Region} = search.initial - function bisect(r::Root) Y1, Y2 = bisect(interval(r)) return Root(Y1, :unknown), Root(Y2, :unknown) end -bisect(::BBSearch, r::Root) = bisect(r::Root) - -function process(search::Union{BreadthFirstSearch, DepthFirstSearch}, r::Root) - contracted_root = search.contractor(r, search.tol) +function process(contractor, root_problem, r::Root) + contracted_root = contractor(r, root_problem.abstol) status = root_status(contracted_root) status == :unique && return :store, contracted_root - status == :empty && return :discard, contracted_root + status == :empty && return :prune, contracted_root if status == :unknown # Avoid infinite division of intervals with singularity - isnan(contracted_root) && diam(r) < search.tol && return :store, r - diam(contracted_root) < search.tol && return :store, contracted_root + isnan(contracted_root) && diam(r) < root_problem.abstol && return :store, r + diam(contracted_root) < root_problem.abstol && return :store, contracted_root - return :bisect, r + return :branch, r else error("Unrecognized root status: $status") end end """ - branch_and_prune(X, contractor, strategy, tol) + branch_and_prune(X, contractor, search_order, tol) Generic branch and prune routine for finding isolated roots using the given contractor to determine the status of a given box `X`. @@ -79,24 +43,26 @@ contractor to determine the status of a given box `X`. See the documentation of the `roots` function for explanation of the other arguments. """ -function branch_and_prune(r::Root, contractor, search, tol) - iter = search(r, contractor, tol) - local endstate - # complete iteration - for state in iter - endstate = state - end - return data(endstate) +function branch_and_prune(r::Root, contractor, search_order, tol) + root_problem = RootProblem(tol) + search = BranchAndPruneSearch( + search_order, + X -> process(contractor, root_problem, X), + bisect, + r + ) + result = bpsearch(search) + return vcat(result.final_regions, result.unfinished_regions) end const NewtonLike = Union{Type{Newton}, Type{Krawczyk}} -const default_strategy = DepthFirstSearch +const default_search_order = DepthFirst const default_tolerance = 1e-7 const default_contractor = Newton """ - roots(f, X, contractor=Newton, strategy=BreadthFirstSearch, tol=1e-15) - roots(f, deriv, X, contractor=Newton, strategy=BreadthFirstSearch, tol=1e-15) + roots(f, X, contractor=Newton, search_order=BreadthFirst, tol=1e-15) + roots(f, deriv, X, contractor=Newton, search_order=BreadthFirst, tol=1e-15) roots(f, X, contractor, tol) roots(f, deriv, X, contractor, tol) @@ -110,36 +76,36 @@ Inputs: the status of a given box `X`. It returns the new box and a symbol indicating the status. Current possible values are `Bisection`, `Newton` and `Krawczyk` - `deriv`: explicit derivative of `f` for `Newton` and `Krawczyk` - - `strategy`: `SearchStrategy` determining the order in which regions are + - `search_order`: `SearchStrategy` determining the order in which regions are processed. - `tol`: Absolute tolerance. If a region has a diameter smaller than `tol`, it is returned with status `:unknown`. """ function roots(f::Function, X, contractor::Type{C}=default_contractor, - strategy::Type{S}=default_strategy, - tol::Float64=default_tolerance) where {C <: Contractor, S <: BBSearch} + search_order::Type{S}=default_search_order, + tol::Float64=default_tolerance) where {C <: Contractor, S <: SearchOrder} - _roots(f, X, contractor, strategy, tol) + _roots(f, X, contractor, search_order, tol) end function roots(f::Function, deriv::Function, X, contractor::Type{C}=default_contractor, - strategy::Type{S}=default_strategy, - tol::Float64=default_tolerance) where {C <: Contractor, S <: BBSearch} + search_order::Type{S}=default_search_order, + tol::Float64=default_tolerance) where {C <: Contractor, S <: SearchOrder} - _roots(f, deriv, X, contractor, strategy, tol) + _roots(f, deriv, X, contractor, search_order, tol) end function roots(f::Function, X, contractor::Type{C}, tol::Float64) where {C <: Contractor} - _roots(f, X, contractor, default_strategy, tol) + _roots(f, X, contractor, default_search_order, tol) end function roots(f::Function, deriv::Function, X, contractor::Type{C}, tol::Float64) where {C <: Contractor} - _roots(f, deriv, X, contractor, default_strategy, tol) + _roots(f, deriv, X, contractor, default_search_order, tol) end #=== @@ -149,99 +115,99 @@ end # For `Bisection` method function _roots(f, r::Root{T}, ::Type{Bisection}, - strategy::Type{S}, tol::Float64) where {T, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, S <: SearchOrder} - branch_and_prune(r, Bisection(f), strategy, tol) + branch_and_prune(r, Bisection(f), search_order, tol) end # For `NewtonLike` acting on `Interval` function _roots(f, r::Root{Interval{T}}, contractor::NewtonLike, - strategy::Type{S}, tol::Float64) where {T, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, S <: SearchOrder} deriv = x -> ForwardDiff.derivative(f, x) - _roots(f, deriv, r, contractor, strategy, tol) + _roots(f, deriv, r, contractor, search_order, tol) end function _roots(f, deriv, r::Root{Interval{T}}, contractor::NewtonLike, - strategy::Type{S}, tol::Float64) where {T, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, S <: SearchOrder} - branch_and_prune(r, contractor(f, deriv), strategy, tol) + branch_and_prune(r, contractor(f, deriv), search_order, tol) end # For `NewtonLike` acting on `IntervalBox` function _roots(f, r::Root{IntervalBox{N, T}}, contractor::NewtonLike, - strategy::Type{S}, tol::Float64) where {N, T, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {N, T, S <: SearchOrder} deriv = x -> ForwardDiff.jacobian(f, x) - _roots(f, deriv, r, contractor, strategy, tol) + _roots(f, deriv, r, contractor, search_order, tol) end function _roots(f, deriv, r::Root{IntervalBox{N, T}}, contractor::NewtonLike, - strategy::Type{S}, tol::Float64) where {N, T, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {N, T, S <: SearchOrder} - branch_and_prune(r, contractor(f, deriv), strategy, tol) + branch_and_prune(r, contractor(f, deriv), search_order, tol) end # Acting on `Interval` function _roots(f, X::Region, contractor::Type{C}, - strategy::Type{S}, tol::Float64) where {C <: Contractor, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {C <: Contractor, S <: SearchOrder} - _roots(f, Root(X, :unknown), contractor, strategy, tol) + _roots(f, Root(X, :unknown), contractor, search_order, tol) end function _roots(f, deriv, X::Region, contractor::Type{C}, - strategy::Type{S}, tol::Float64) where {C <: Contractor, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {C <: Contractor, S <: SearchOrder} - _roots(f, deriv, Root(X, :unknown), contractor, strategy, tol) + _roots(f, deriv, Root(X, :unknown), contractor, search_order, tol) end # Acting on `Vector` of `Root` function _roots(f, V::Vector{Root{T}}, contractor::Type{C}, - strategy::Type{S}, tol::Float64) where {T, C <: Contractor, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, C <: Contractor, S <: SearchOrder} - vcat(_roots.(f, V, contractor, strategy, tol)...) + vcat(_roots.(f, V, contractor, search_order, tol)...) end function _roots(f, deriv, V::Vector{Root{T}}, contractor::Type{C}, - strategy::Type{S}, tol::Float64) where {T, C <: Contractor, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, C <: Contractor, S <: SearchOrder} - vcat(_roots.(f, deriv, V, contractor, strategy, tol)...) + vcat(_roots.(f, deriv, V, contractor, search_order, tol)...) end # Acting on complex `Interval` function _roots(f, Xc::Complex{Interval{T}}, contractor::Type{C}, - strategy::Type{S}, tol::Float64) where {T, C <: Contractor, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, C <: Contractor, S <: SearchOrder} g = realify(f) Y = IntervalBox(reim(Xc)...) - rts = _roots(g, Root(Y, :unknown), contractor, strategy, tol) + rts = _roots(g, Root(Y, :unknown), contractor, search_order, tol) return [Root(Complex(root.interval...), root.status) for root in rts] end function _roots(f, Xc::Complex{Interval{T}}, contractor::NewtonLike, - strategy::Type{S}, tol::Float64) where {T, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, S <: SearchOrder} g = realify(f) g_prime = x -> ForwardDiff.jacobian(g, x) Y = IntervalBox(reim(Xc)...) - rts = _roots(g, g_prime, Root(Y, :unknown), contractor, strategy, tol) + rts = _roots(g, g_prime, Root(Y, :unknown), contractor, search_order, tol) return [Root(Complex(root.interval...), root.status) for root in rts] end function _roots(f, deriv, Xc::Complex{Interval{T}}, contractor::NewtonLike, - strategy::Type{S}, tol::Float64) where {T, S <: BBSearch} + search_order::Type{S}, tol::Float64) where {T, S <: SearchOrder} g = realify(f) g_prime = realify_derivative(deriv) Y = IntervalBox(reim(Xc)...) - rts = _roots(g, g_prime, Root(Y, :unknown), contractor, strategy, tol) + rts = _roots(g, g_prime, Root(Y, :unknown), contractor, search_order, tol) return [Root(Complex(root.interval...), root.status) for root in rts] end diff --git a/test/branch_and_bound.jl b/test/branch_and_bound.jl deleted file mode 100644 index 58372e9f..00000000 --- a/test/branch_and_bound.jl +++ /dev/null @@ -1,134 +0,0 @@ -import IntervalRootFinding: BBLeaf, BBNode, BBTree -import IntervalRootFinding: root, nnodes, data, newid, discard_leaf! - -import IntervalRootFinding: BreadthFirstBBSearch -import IntervalRootFinding: process, bisect - -@testset "Branch and bound tree" begin - #= Manually create a simple tree with the following structure - - (1) - | - (2) - =# - leaf = BBLeaf("I'm a leaf", 1, :final) - node = BBNode(0, [2]) - tree = BBTree(Dict(1 => node), - Dict(2 => leaf), - [2]) - - # Check that printing does not error - io = IOBuffer() - println(io, leaf) - println(io, node) - println(io, tree) - - @test root(tree) == node - @test nnodes(tree) == 2 - - d = data(tree) - @test length(d) == 1 - @test d[1] == "I'm a leaf" - - k = newid(tree) - @test k ∉ 0:2 # The new ID must be new and not 0 - @test isa(k, Integer) - - discard_leaf!(tree, 2) - # The last leaf was deleted so the tree should be empty - @test nnodes(tree) == 0 - - #= Tests on a slighty more complicated tree - - (1) - / \ - (2) (8) - / \ \ - (3) (5) (9) - | | \ \ - (4) (6) (7) (10) - =# - - n1 = BBNode(0, [2, 8]) - n2 = BBNode(1, [3, 5]) - n3 = BBNode(2, [4]) - n4 = BBLeaf("Leaf 4", 3, :working) - n5 = BBNode(2, [6, 7]) - n6 = BBLeaf("Leaf 6", 5, :final) - n7 = BBLeaf("Leaf 7", 5, :working) - n8 = BBNode(1, [9]) - n9 = BBNode(8, [10]) - n10 = BBLeaf("Leaf 10", 9, :working) - - tree = BBTree(Dict(1 => n1, 2 => n2, 3 => n3, 5 => n5, 8 => n8, 9 => n9), - Dict(4 => n4, 6 => n6, 7 => n7, 10 => n10), - [3, 4, 8]) - - # Check that printing does not error - io = IOBuffer() - println(io, tree) - - @test newid(tree) ∉ 0:10 - @test length(data(tree)) == 4 - - @test nnodes(tree) == 10 - - discard_leaf!(tree, 6) - @test nnodes(tree) == 9 - - discard_leaf!(tree, 4) - @test nnodes(tree) == 7 - - discard_leaf!(tree, 10) - @test nnodes(tree) == 4 - - @test length(data(tree)) == 1 -end - - -# Implement the interface for breadth first in a dummy way -struct DummySearch <: BreadthFirstBBSearch{Symbol} end - -process(::DummySearch, s::Symbol) = s, s -bisect(::DummySearch, s::Symbol) = :store, :store - -@testset "Interface functions" begin - #= Build the following tree for testing - (1) - / \ - (2) (5: to_bisect) - / \ - (3: to_store) (4: to_discard) - =# - rt = BBNode(0, [2, 5]) - node = BBNode(1, [3, 4]) - to_store = BBLeaf(:store, 2, :working) - to_discard = BBLeaf(:discard, 2, :working) - to_bisect = BBLeaf(:bisect, 1, :working) - - tree = BBTree(Dict(1 => rt, 2 => node), - Dict(3 => to_store, 4 => to_discard, 5 => to_bisect), - [3, 4, 5]) - - search = DummySearch() - - # First iteration processes the to_store leaf - tree, _ = iterate(search, tree) - - @test nnodes(tree) == 5 - @test tree[3].status == :final - @test length(tree.working_leaves) == 2 - - # Second iteration processes the to_discard leaf - tree, _ = iterate(search, tree) - - @test nnodes(tree) == 4 - @test length(tree.working_leaves) == 1 - - # Second iteration processes the to_bisect leaf - tree, _ = iterate(search, tree) - - @test nnodes(tree) == 6 - @test length(tree.working_leaves) == 2 - @test isa(tree[5], BBNode) -end diff --git a/test/runtests.jl b/test/runtests.jl index c49dba91..35867a48 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,6 @@ using IntervalRootFinding using Test -include("branch_and_bound.jl") -include("search_interface.jl") include("roots.jl") include("test_smiley.jl") include("newton1d.jl") diff --git a/test/search_interface.jl b/test/search_interface.jl deleted file mode 100644 index a969ee54..00000000 --- a/test/search_interface.jl +++ /dev/null @@ -1,21 +0,0 @@ -using IntervalArithmetic -using IntervalRootFinding - -@testset "Root search iterator interface" begin - contractor = Newton(sin, cos) - search = BreadthFirstSearch(-10..10, contractor, 1e-3) - - # First iteration of the search - elem, state = iterate(search) - - # Optional iterator methods - @test eltype(search) != Any -end - - -@testset "Search strategies" begin - # DepthFristSearch and BreadthFristSearch should return the same results - X = -5..5 - rts = roots(sin, cos, X, Newton, DepthFirstSearch) - @test Set(rts) == Set(roots(sin, cos, X, Newton, BreadthFirstSearch)) -end