Use branch and prune as backend for the search

[Kolaru]

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]

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

[lucaferranti]

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

[lucaferranti]

out of curiosity, do you have any benchmarks / insights / literature references on how the two search strategies compare for interval root finding? I guess they are pretty much equivalent and there's not a strong advantage of one over the other?

[Kolaru]

As far as I know the main difference is that is you you have a limited number of iteration, going breadth first give you evenly sized intervals over the search domain, while going depth first gives you smallest possible intervals for the subset of the search domain that it had time to go through.

[dpsanders]

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]

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.

[lbenet]

LGTM. Go ahead and merge!

[Kolaru]

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.