Implement max_iteration and reltol stop criterion

docs/src/roots.md

@@ -136,3 +136,5 @@ The last example shows a case where the tolerance was too large to be able to is
 !!! warning
 
     For a root `x` of some function, if the absolute tolerance is smaller than `eps(x)` i.e. if `tol + x == x`, `roots` may never be able to converge to the required tolerance and the function may get stuck in an infinite loop.
+    To avoid that, either increase `abstol` or `reltol`,
+    or set a maximum number of iterations with the `max_iteration` keyword.

src/contractors.jl

@@ -93,7 +93,7 @@ end
 Refine a root.
 """
 function refine(root_problem::RootProblem{C}, R::Root) where C
-    root_status(R) != :unique && return R
+    root_status(R) != :unique && throw(ArgumentError("trying to refine a non unique root"))
 
     X = root_region(R)
 

src/region.jl

@@ -23,6 +23,12 @@ isnai_region(X::AbstractVector) = any(isnai, X)
 diam_region(X::Interval) = diam(X)
 diam_region(X::AbstractVector) = maximum(diam, X)
 
+mag_region(X::Interval) = mag(X)
+mag_region(X::AbstractVector) = maximum(mag, X)
+
+mig_region(X::Interval) = mig(X)
+mig_region(X::AbstractVector) = minimum(mig, X)
+
 bisect_region(X::Interval, α) = bisect(X, α)
 
 function bisect_region(X::AbstractVector, α)

src/root_object.jl

@@ -47,6 +47,8 @@ show(io::IO, rt::Root) = print(io, "Root($(rt.region), :$(rt.status))")
 ⊆(a::Root, b::Root) = a.region ⊆ b.region
 
 IntervalArithmetic.diam(r::Root) = diam_region(root_region(r))
+IntervalArithmetic.mag(r::Root) = mag_region(root_region(r))
+IntervalArithmetic.mig(r::Root) = mig_region(root_region(r))
 IntervalArithmetic.isnai(r::Root) = isnai_region(root_region(r))
 
 function Base.:(==)(r1::Root, r2::Root)

src/roots.jl

@@ -1,6 +1,4 @@
 
-import IntervalArithmetic: diam, bisect, isnai
-
 export branch_and_prune, Bisection, Newton
 
 struct RootProblem{C, F, G, R, S, T}
@@ -40,14 +38,17 @@ Parameters
 - `search_order`: Order in which the sub-regions are searched.
     `BreadthFirst` (visit the largest regions first) and `DepthFirst`
     (visit the smallest regions first) are supported. Default: `BreadthFirst`.
-- `abstol`: Absolute tolerance. The search is stopped when all dimension
-    of the remaining regions are smaller than `abstol`. Default: 1e-7.
+- `abstol`: Absolute tolerance. The search is stopped when all dimensions
+    of the remaining regions are smaller than `abstol`. Default: `1e-7`.
+- `reltol`: Relative tolerance. The search is stopped when all dimensions
+    of the remaining regions are smaller than `reltol * mag(interval)`.
+    Default: `0.0``.
+- `max_iteration`: The maximum number of iteration, which also corresponds to
+    the maximum number of bisections allowed. Default: `typemax(Int)`.
 - `where_bisect`: Value used to bisect the region. It is used to avoid
     bisecting exactly on zero when starting with symmetrical regions,
     often leading to having a solution directly on the boundary of a region,
-    which prevent the contractor to prove it's unicity. Default: 127/256.
-
-`reltol` and `max_iteration` are currently ignored.
+    which prevent the contractor to prove it's unicity. Default: `127/256`.
 """
 RootProblem(f, region ; kwargs...) = RootProblem(f, Root(region, :unkown) ; kwargs...)
 
@@ -57,19 +58,10 @@ function RootProblem(
         derivative = nothing,
         search_order = BreadthFirst,
         abstol = 1e-7,
-        reltol = nothing,
-        max_iteration = nothing,
+        reltol = 0.0,
+        max_iteration = typemax(Int),
         where_bisect = 0.49609375)  # 127//256
 
-    if !isnothing(reltol) || !isnothing(max_iteration)
-        throw(
-            ArgumentError("reltol and max_iteration not yet implemented")
-        )
-    else
-        reltol = 0.0
-        max_iteration = -1
-    end
-
     N = length(root_region(root))
     if isnothing(derivative)
         if N == 1
@@ -110,23 +102,32 @@ function bisect_region(r::Root, α)
     return Root(Y1, :unknown), Root(Y2, :unknown)
 end
 
-function process(root_problem, r::Root)
-    contracted_root = contract(root_problem, r)
-    refined_root = refine(root_problem, contracted_root)
+function under_tolerance(root_problem, root::Root)
+    d = diam(root)
+    d < root_problem.abstol && return true
+    d / mag(root) < root_problem.reltol && return true
+    return false
+end
+
+function process(root_problem, root::Root)
+    contracted = contract(root_problem, root)
+    status = root_status(contracted)
 
-    status = root_status(refined_root)
+    if status == :unique 
+        refined_root = refine(root_problem, contracted)
+        return :store, refined_root
+    end
 
-    status == :unique && return :store, refined_root
-    status == :empty && return :prune, refined_root
+    status == :empty && return :prune, root
 
     if status == :unknown
         # Avoid infinite division of intervals with singularity
-        isnai(refined_root) && diam(r) < root_problem.abstol && return :store, r
-        diam(refined_root) < root_problem.abstol && return :store, refined_root
-
-        return :branch, r
+        istrivial(contracted.region) && under_tolerance(root_problem, root) && return :store, root
+        under_tolerance(root_problem, contracted) && return :store, contracted
+        
+        return :branch, root
     else
-        error("Unrecognized root status: $status")
+        throw(ArgumentError("unrecognized root status: $status"))
     end
 end
 
@@ -157,8 +158,21 @@ For information about the optional search parameters,
 see [`RootProblem`](@ref).
 """
 function roots(f, region ; kwargs...)
-    search = root_search(RootProblem(f, region ; kwargs...))
-    result = bpsearch(search)
+    problem = RootProblem(f, region ; kwargs...)
+    search = root_search(problem)
+    endstate = nothing
+
+    for (iter, state) in enumerate(search)
+        endstate = state
+        iter >= problem.max_iteration && break
+    end
+
+    result = BranchAndPrune.BranchAndPruneResult(
+        endstate.search_order,
+        search.initial_region,
+        endstate.tree
+    )
+
     return vcat(result.final_regions, result.unfinished_regions)
 end
 

test/roots.jl

@@ -286,4 +286,28 @@ end
         @test eltype(rS1) <: Root{<:SVector}
         @test eltype(rS2) <: Root{<:SVector}
     end
+end
+
+@testset "Stop criteria" begin
+    abstols = [1e-4, 1e-7, 1e-10]
+    reltols = [1e-4, 1e-7, 1e-10]
+
+    f(x) = x*sin(1/x)
+
+    for abstol in abstols
+        for reltol in reltols
+            rts = roots(f, interval(0, 1) ; abstol, reltol)
+            regions = [rt.region for rt in rts if root_status(rt) == :unknown]
+
+            @test all(
+                (diam.(regions) .< abstol) .|
+                (diam.(regions) .< (reltol .* mag.(regions)))
+            )
+        end
+    end
+
+    for max_iteration in [10, 100, 1000, 10_000]
+        rts = roots(f, interval(0, 1) ; abstol = 1e-10, max_iteration)
+        @test length(rts) <= max_iteration
+    end
 end