JuliaIntervals / IntervalRootFinding.jl

Library for finding the roots of a function using interval arithmetic
https://juliaintervals.github.io/IntervalRootFinding.jl/
Other
129 stars 26 forks source link

Use branch and prune as backend for the search #187

Closed Kolaru closed 1 year ago

Kolaru commented 1 year ago

Replaces the internal branch and prune search by BranchAndPrune.jl v0.2.

This is the first of 3 planned PR

  1. Use BranchAndPrune package
  2. Refactor contractors to dissociate them totally from the iteration configuration (only the absolute tolerance for now)
  3. Refactor the whole interface, using a RootProblem struct for better handling of iteration config, and add options like max iteration and relative tolerance
Kolaru commented 1 year ago

@lucaferranti @lbenet @dpsanders If there is no strong feeling against this, I'd like to proceed in the coming days.

lucaferranti commented 1 year ago

Sorry! overall looks good. I'll give a more detailed review during the weekend. If I faik to review by monday, feel free to merge then

dpsanders commented 1 year ago

In general I think this is a great idea, but it would be good to benchmark to make sure that performance is OK - at some point I remember comparing with a simple implementation which was considerably faster.

I seem to remember that there is a separate implementation of Newton in 1D that it would be useful to keep.

lbenet commented 1 year ago

I seem to remember that there is a separate implementation of Newton in 1D that it would be useful to keep.

It is still there: src/newton1d.jl.

Kolaru commented 1 year ago

In general I think this is a great idea, but it would be good to benchmark to make sure that performance is OK - at some point I remember comparing with a simple implementation which was considerably faster.

To get a bit of peace of mind I quickly benchmarked the general method versus newton1d and it is seems to be faring pretty well (the previous implementation not using branch and prune is on a quite old version). newton1d do some additional checks, so it is not totally surprising it is slower.

julia> f(x) = x^2 - 2
f (generic function with 1 method)

julia> df(x) = 2x
df (generic function with 1 method)

julia> abstol = 1e-7
1.0e-7

julia> X = -5..5
[-5, 5]

julia> @benchmark newton1d(f, df, X ; abstol)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  166.700 μs …  15.693 ms  ┊ GC (min … max): 0.00% … 73.56%
 Time  (median):     175.000 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   214.903 μs ± 588.151 μs  ┊ GC (mean ± σ):  8.09% ±  2.92%

  ▇█▇▅▄▄▄▂▄▃▃▃▂▁▁▁▁                                             ▂
  █████████████████████▇█▇██▇▇█▇▇▇▇▇▇▇▇▇▇▆▆▆▆▆▆▆▆▅▄▆▅▆▆▄▆▆▄▄▄▄▄ █
  167 μs        Histogram: log(frequency) by time        402 μs <

 Memory estimate: 86.91 KiB, allocs estimate: 2086.

julia> @benchmark roots(f, df, X, Newton, abstol)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  137.800 μs …  17.256 ms  ┊ GC (min … max): 0.00% … 78.47%
 Time  (median):     146.900 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   179.884 μs ± 568.193 μs  ┊ GC (mean ± σ):  8.74% ±  2.74%

  ██▇▆▆▅▅▃▃▃▃▂▂▁  ▁                                             ▂
  ████████████████████▇▇▇▇▇▇▇▆▇▆▇▇▇▇▇▆▆▆▅▅▅▆▅▅▅▆▅▆▆▅▅▅▅▅▅▆▅▄▄▅▄ █
  138 μs        Histogram: log(frequency) by time        359 μs <

 Memory estimate: 73.88 KiB, allocs estimate: 1774.