Repository: JuliaMath/Roots.jl (julia)

Issue:
Tracker works poorly for bracketing methods
See the following output:
```
julia> for solver in (Roots.Secant(), Roots.Bisection(), Roots.FalsePosition(), Roots.ITP())
           tracker = Roots.Tracks()
           find_zero(f, (0, 2), solver, tracks=tracker, atol=1e-1)
           println()
           @show solver
           @show tracker.fs
           @show length(tracker.xs)
           @show length(tracker.fs)
       end

solver = Secant()
tracker.fs = [-2.0, 2.0, -1.0, -0.22222222222222232, 0.04081632653061229]
length(tracker.xs) = 5
length(tracker.fs) = 5

solver = Bisection()
tracker.fs = [-2.0, 2.0]
length(tracker.xs) = 28
length(tracker.fs) = 2

solver = FalsePosition{:anderson_bjork}()
tracker.fs = [-2.0, 2.0]
length(tracker.xs) = 8
length(tracker.fs) = 2

solver = Roots.ITP{Float64, Int64}(0.2, 2, 1)
tracker.fs = Float64[]
length(tracker.xs) = 12
length(tracker.fs) = 0
```

The number of, and value of, the `fs` field of the tracker are plain nonsense for the bracketing methods. I have no idea what goes wrong, but a fix should not be complicated, seeing as the algorithms compute good solutions.

Pull request description:
Following #278 this reworks the `Tracks` structure:

* rather than use `xs` and `fs` to store values; this uses `xfs` and `abs` to store values depending on if the method is bracketing or not
* this needs improved documentation
* Could the interface of `ValueHistories` be useful here?
* This also addresses part of issue #279 (ITP is corrected), but the initial step for Bisection is lost
* this would close #277 too?
* Not sure about the suggestion in #276

Patch:
diff --git a/docs/src/reference.md b/docs/src/reference.md
index 36aa6e3c..ea94609d 100644
--- a/docs/src/reference.md
+++ b/docs/src/reference.md
@@ -297,3 +297,13 @@ The initial naming scheme used `fzero` instead  of `fzeros`, following the name
 ```@docs
 Roots.fzero
 ```
+
+## Tracking iterations
+
+It is possible to add the keyword arguement `verbose=true` to when calling the `find_zero` function to get detailed information about the solution, and data from each iteration. If you want to save this data instead of just printing it, you can use a `Tracks` object.
+
+----
+
+```@docs
+Roots.Tracks
+```
diff --git a/src/bracketing.jl b/src/bracketing.jl
index bdac032c..482cfda1 100644
--- a/src/bracketing.jl
+++ b/src/bracketing.jl
@@ -42,27 +42,12 @@ struct Bisection <: AbstractBisection end  # either solvable or A42
 struct BisectionExact <: AbstractBisection end
 
 ## tracks for bisection, different, we show bracketing interval
-function log_step(l::Tracks, M::AbstractBracketing, state)
-    push!(l.xs, state.xn0)
-    push!(l.xs, state.xn1) # we store [ai,bi, ai+1, bi+1, ...]
-    log_steps(l)
-end
-function show_tracks(io::IO, l::Tracks, M::AbstractBracketing)
-    xs = l.xs
-    n = length(xs)
-    for (i, j) in enumerate(1:2:(n - 1))
-        println(
-            io,
-            @sprintf(
-                "(%s, %s) = (% 18.16f, % 18.16f)",
-                "a_$(i-1)",
-                "b_$(i-1)",
-                xs[j],
-                xs[j + 1]
-            )
-        )
-    end
-    println(io, "")
+## No init here; for Bisection() [a₀, b₀] is just lost.
+function log_step(l::Tracks, M::AbstractBracketing, state; init::Bool=false)
+    a, b = state.xn0, state.xn1
+    push!(l.abₛ, a < b ? (a,b) : (b,a))
+    !init && log_iteration(l, 1)
+    nothing
 end
 
 ## helper function: floating point, sorted, finite
@@ -130,6 +115,14 @@ function default_tolerances(::AbstractBisection, ::Type{T}, ::Type{S}) where {T,
     (xatol, xrtol, atol, rtol, maxevals, maxfnevals, strict)
 end
 
+function log_step(l::Tracks, M::Bisection, state; init::Bool=false)
+    a, b = state.xn0, state.xn1
+    push!(l.abₛ, (a,b))
+    init && log_iteration(l, 1) # c is computed
+    !init && log_iteration(l, 1)
+    nothing
+end
+
 # find middle of (a,b) with convention that
 # * if a, b finite, they are made non-finite
 # if a,b of different signs, middle is 0
@@ -468,12 +461,13 @@ function default_tolerances(::AbstractAlefeldPotraShi, ::Type{T}, ::Type{S}) whe
 end
 
 ## initial step, needs to log a,b,d
-function log_step(l::Tracks, M::AbstractAlefeldPotraShi, state, ::Any)
+function log_step(l::Tracks, M::AbstractAlefeldPotraShi, state; init::Bool=false)
     a, b, c = state.xn0, state.xn1, state.d
-    append!(l.xs, extrema((a, b, c)))
-    push!(l.xs, a)
-    push!(l.xs, b) # we store [ai,bi, ai+1, bi+1, ...] for brackecting methods
-    log_steps(l, 1)
+    init && push!(l.abₛ, extrema((a, b, c)))
+    init && log_iteration(l, 1) # take an initial step
+    push!(l.abₛ, (a,b))
+    !init && log_iteration(l, 1)
+    nothing
 end
 
 struct A42State{T,S} <: AbstractUnivariateZeroState{T,S}
@@ -762,6 +756,7 @@ struct BrentState{T,S} <: AbstractUnivariateZeroState{T,S}
     mflag::Bool
 end
 
+
 # # we store mflag as -1, or +1 in state.mflag
 function init_state(::Brent, F, x₀, x₁, fx₀, fx₁)
     u, v, fu, fv = x₀, x₁, fx₀, fx₁
@@ -850,14 +845,6 @@ function update_state(
 end
 
 
-function log_step(l::Tracks, M::Brent, state)
-    a, b = state.xn0, state.xn1
-    u, v = a < b ? (a, b) : (b, a)
-    push!(l.xs, a)
-    push!(l.xs, b) # we store [ai,bi, ai+1, bi+1, ...]
-    log_steps(l)
-end
-
 
 ## --------------------------------------------------
 
@@ -872,14 +859,14 @@ Example:
 ```jldoctest
 julia> using Roots
 
-julia> find_zero(x -> exp(x) - x^4, (5, 15), Roots.Ridders())
-8.6131694564414
+julia> find_zero(x -> exp(x) - x^4, (5, 15), Roots.Ridders()) ≈ 8.61316945644
+true
 
-julia> find_zero(x -> x*exp(x) - 10, (-100, 100), Roots.Ridders())
-1.7455280027406994
+julia> find_zero(x -> x*exp(x) - 10, (-100, 100), Roots.Ridders()) ≈ 1.74552800274
+true
 
-julia> find_zero(x -> tan(x)^tan(x) - 1e3, (0, 1.5), Roots.Ridders())
-1.3547104419635592
+julia> find_zero(x -> tan(x)^tan(x) - 1e3, (0, 1.5), Roots.Ridders()) ≈ 1.3547104419
+true
 ```
 
 [Ridders](https://cs.fit.edu/~dmitra/SciComp/Resources/RidderMethod.pdf) showed the error satisfies `eₙ₊₁ ≈ 1/2 eₙeₙ₋₁eₙ₋₂ ⋅ (g^2-2fh)/f` for
@@ -965,6 +952,7 @@ struct ITP{T,S} <: AbstractAcceleratedBisection
 end
 ITP(;κ₁ = 0.2, κ₂ = 2, n₀=1) = ITP(κ₁, κ₂, n₀)
 
+
 struct ITPState{T,S,R} <: AbstractUnivariateZeroState{T,S}
     xn1::T
     xn0::T
diff --git a/src/find_zero.jl b/src/find_zero.jl
index b4766325..6b5394de 100644
--- a/src/find_zero.jl
+++ b/src/find_zero.jl
@@ -149,113 +149,159 @@ end
 # # reset options to default values
 @deprecate init_options!(options, M) init_options(M)
 
+### --------------------------------------------------
+
 ## Tracks (for logging actual steps)
 ## when no logging this should get optimized out
+## when logging, this allocates
 abstract type AbstractTracks end
 struct NullTracks <: AbstractTracks end
-# api
-log_step(s::NullTracks, M, x, init=false) = nothing
-log_steps(::NullTracks, n=1) = nothing
-incfn(::NullTracks, i=1) = nothing
-log_message(::NullTracks, msg) = nothing
-log_convergence(::NullTracks, msg) = nothing
-log_state(::NullTracks, state) = nothing
-log_method(::NullTracks, method) = nothing
-log_nmethod(::NullTracks, method) = nothing
+# loggin api
+log_step(s::NullTracks, M, x; init=false) = nothing  # log a step (x,f(x)) or (a,b)
+log_iteration(::NullTracks, n=1) = nothing           # log an iteration
+incfn(::NullTracks, i=1) = nothing                   # add a function call (should rename log_fncall); only one called in update_step
+log_message(::NullTracks, msg) = nothing             # append to a message
+log_convergence(::NullTracks, msg) = nothing         # flag for convergence
+log_last(::NullTracks, state) = nothing              # log α
+log_method(::NullTracks, method) = nothing           # record M
+log_nmethod(::NullTracks, method) = nothing          # record N (if hybrid)
 
 # a tracks object to record tracks
 """
     Tracks(T, S)
+    Tracks()
+
 
 Construct a `Tracks` object used to record the progress of the algorithm.
+`T` is the type of function inputs, and `S` is the type of function outputs. They
+both default to `Float64`. Note that because this type is not exported, you have to
+write `Roots.Tracks()` to construct a `Tracks` object.
 
-By default, a null tracks object used, but if `verbose=true` is specified of a `Tracks` object passed to the keyword `tracks` then the steps of the algorithm are recorded in the `Tracks` object. If `verbose=true`, the `Tracks` object will be shown.
+By default, no tracking is done while finding a root (with `find_zero`, `solve`, or `fzero`).
+To change this, construct a `Tracks` object, and pass it to the keyword argument `tracks`.
+This will modify the `Tracks` object, storing the input and function values at each iteration,
+along with additional information about the root-finding process.
 
-Alternatively, a `Tracks` object can be constructed and then shown; `find_zero` modifies a `Tracks` object.
+`Tracker` objects are printed in a nice, easy-to-read
+format. Internally either a tuple of `(x,f(x))` pairs or `(aₙ, bₙ)`
+pairs are stored, the latter for bracketing methods, as illustrated
+in the examples below. (These implementation details may change without notice.)
 
-## Example
+If you only want to print the information, but you don't need it later, this can conveniently be
+done by passing `verbose=true` to the root-finding function. This will not
+effect the return value, which will still be the root of the function.
 
-(This example uses a default type for `Tracks` which may not be suitable for some uses.)
+## Examples
 
-```
-julia> tracks = Roots.Tracks()  # default type is Float64 for both x, fx values
-julia> find_zero(sin, 3, Secant(); tracks=tracks)
-julia> tracks.steps # 4
-julia> tracks  # show summary of algorithm
+```jldoctest Tracks
+julia> using Roots
+
+julia> f(x) = x^2-2
+f (generic function with 1 method)
+
+julia> tracker = Roots.Tracks()
+Algorithm has not been run
+
+julia> find_zero(f, (0, 2), Roots.Secant(), tracks=tracker) ≈ √2
+true
+
+julia> tracker
+Results of univariate zero finding:
+
+* Converged to: 1.4142135623730947
+* Algorithm: Secant()
+* iterations: 7
+* function evaluations ≈ 9
+* stopped as |f(x_n)| ≤ max(δ, |x|⋅ϵ) using δ = atol, ϵ = rtol
+
+Trace:
+x_0 =  0.0000000000000000,	 fx_0 = -2.0000000000000000
+x_1 =  2.0000000000000000,	 fx_1 =  2.0000000000000000
+x_2 =  1.0000000000000000,	 fx_2 = -1.0000000000000000
+x_3 =  1.3333333333333333,	 fx_3 = -0.2222222222222223
+x_4 =  1.4285714285714286,	 fx_4 =  0.0408163265306123
+x_5 =  1.4137931034482758,	 fx_5 = -0.0011890606420930
+x_6 =  1.4142114384748701,	 fx_6 = -0.0000060072868389
+x_7 =  1.4142135626888697,	 fx_7 =  0.0000000008931456
+x_8 =  1.4142135623730947,	 fx_8 = -0.0000000000000009
+
+julia> tracker.xfₛ  # stored as (x, f(x)) pairs
+9-element Vector{Tuple{Float64, Float64}}:
+ (0.0, -2.0)
+ (2.0, 2.0)
+ (1.0, -1.0)
+ (1.3333333333333333, -0.22222222222222232)
+ (1.4285714285714286, 0.04081632653061229)
+ (1.4137931034482758, -0.0011890606420930094)
+ (1.41421143847487, -6.007286838860537e-6)
+ (1.4142135626888697, 8.931455575122982e-10)
+ (1.4142135623730947, -8.881784197001252e-16)
+
+julia> tracker = Roots.Tracks()
+Algorithm has not been run
+
+julia> find_zero(sin, (3, 4), Roots. A42(), tracks=tracker) ≈ π
+true
+
+julia> tracker.abₛ # stored as (aₙ, bₙ) pairs
+4-element Vector{Tuple{Float64, Float64}}:
+ (3.0, 4.0)
+ (3.0, 3.5)
+ (3.14156188905231, 3.1416247553172956)
+ (3.141592653589793, 3.1415926535897936)
 ```
 
 """
 mutable struct Tracks{T,S} <: AbstractTracks
-    xs::Vector{T}
-    fs::Vector{S}
+    xfₛ::Vector{Tuple{T,S}} # (x,f(x))
+    abₛ::Vector{Tuple{T,T}} # (aᵢ, bᵢ)
     steps::Int
     fncalls::Int
     convergence_flag::Symbol
     message::String
-    state
+    alpha::T
     method
     nmethod
 end
-Tracks(T, S) = Tracks(T[], S[], 0, 0, :not_converged, "", nothing, nothing, nothing)
+Tracks(T, S) = Tracks(Tuple{T,S}[], Tuple{T,T}[], 0, 0, :not_converged, "",  NaN*zero(T), nothing, nothing)
 Tracks(s::AbstractUnivariateZeroState{T,S}) where {T,S} = Tracks(T, S)
-Tracks(verbose, tracks, state) =
-    (verbose && isa(tracks, NullTracks)) ? Tracks(state) : tracks
+Tracks(verbose, tracks, state::AbstractUnivariateZeroState{T,S}) where {T,S} =
+    (verbose && isa(tracks, NullTracks)) ? Tracks(T,S) : tracks
 Tracks() = Tracks(Float64, Float64) # give default
 
-function log_step(l::Tracks, M::Any, state, init=nothing)
-    if init !== nothing
-        x₀, fx₀ = state.xn0, state.fxn0
-        push!(l.xs, x₀)
-        push!(l.fs, fx₀)
-    end
-    x₁, fx₁ = state.xn1, state.fxn1
-    push!(l.xs, x₁)
-    push!(l.fs, fx₁)
-
-    init === nothing && log_steps(l, 1)
+function log_step(l::Tracks, M::AbstractNonBracketing, state; init=false)
+    init && push!(l.xfₛ, (state.xn0, state.fxn0))
+    push!(l.xfₛ, (state.xn1, state.fxn1))
+    !init  && log_iteration(l, 1)
     nothing
 end
 
 incfn(l::Tracks, i=1) = (l.fncalls += i; nothing)
-log_steps(l::Tracks, n=1) = (l.steps += n; nothing)
+log_iteration(l::Tracks, n=1) = (l.steps += n; nothing)
 log_message(l::Tracks, msg) = (l.message *= msg; nothing)
 log_convergence(l::Tracks, msg) = (l.convergence_flag = msg; nothing)
-log_state(l::Tracks, state) = (l.state = state; nothing)
+log_last(l::Tracks, α) = (l.alpha = α; nothing)
 log_method(l::Tracks, method) = (l.method = method; nothing)
 log_nmethod(l::Tracks, method) = (l.nmethod = method; nothing)
 
-Base.show(io::IO, l::Tracks) = show_trace(io, l.method, l.nmethod, l.state, l)
+Base.show(io::IO, l::Tracks) = show_trace(io, l.method, l.nmethod, l)
 
-function show_tracks(io::IO, s::Tracks, M::AbstractUnivariateZeroMethod)
-    for (i, (xi, fxi)) in enumerate(zip(s.xs, s.fs))
-        println(
-            io,
-            @sprintf(
-                "%s = % 18.16f,\t %s = % 18.16f",
-                "x_$(i-1)",
-                float(xi),
-                "fx_$(i-1)",
-                float(fxi)
-            )
-        )
-    end
-    println(io, "")
-end
 
-function show_trace(io::IO, method, N, state, tracks)
-    if state === nothing
+function show_trace(io::IO, method, N, tracks)
+
+    if length(tracks.xfₛ) == 0 && length(tracks.abₛ) == 0
         print(io, "Algorithm has not been run")
         return nothing
     end
 
-    converged = !isnan(state.xn1)
+    converged = !isnan(tracks.alpha)
     println(io, "Results of univariate zero finding:\n")
     if converged
-        println(io, "* Converged to: $(state.xn1)")
-        if N === nothing || method isa AbstractBracketing
+        println(io, "* Converged to: $(tracks.alpha)")
+        if N === nothing || length(tracks.abₛ) == 0
             println(io, "* Algorithm: $(method)")
         else
-            println(io, "* Algorithm: $(method), with possible bracketing with $N")
+            println(io, "* Algorithm: $(method); finished with bracketing method $N")
         end
         println(io, "* iterations: $(tracks.steps)")
         println(io, "* function evaluations ≈ $(tracks.fncalls)")
@@ -274,6 +320,42 @@ function show_trace(io::IO, method, N, state, tracks)
     show_tracks(io, tracks, method)
 end
 
+function show_tracks(io::IO, s::Tracks, M::AbstractUnivariateZeroMethod)
+
+    # show (x,f(x))
+    for (i, (xi, fxi)) in enumerate(s.xfₛ)
+        println(
+            io,
+            @sprintf(
+                "%s = % 18.16f,\t %s = % 18.16f",
+                "x_$(i-1)",
+                float(xi),
+                "fx_$(i-1)",
+                float(fxi)
+            )
+        )
+    end
+
+    # show bracketing
+    i₀ = length(s.xfₛ)
+    for (i, (a,b)) in enumerate(s.abₛ)
+        j = i₀ + i
+        println(
+            io,
+            @sprintf(
+                "(%s, %s) = (% 18.16f, % 18.16f )",
+                "a_$(j-1)",
+                "b_$(j-1)",
+                a, b
+            )
+        )
+    end
+    println(io, "")
+
+end
+
+## --------------------------------------------------
+
 ### Functions
 # A hacky means to call a function so that parameters can be passed as desired
 # and the correct number of outputs are computed
@@ -963,7 +1045,7 @@ function solve!(P::ZeroProblemIterator; verbose=false)
 
     val, stopped = :not_converged, false
     ctr = 1
-    log_step(l, M, state, :init)
+    log_step(l, M, state; init=true)
 
     while !stopped
         val, stopped = assess_convergence(M, state, options)
@@ -976,15 +1058,17 @@ function solve!(P::ZeroProblemIterator; verbose=false)
         ctr += 1
     end
 
+    val, stopped = assess_convergence(M, state, options) # udpate val flag
+    α = decide_convergence(M, F, state, options, val)
+
     if !(l isa NullTracks)
         log_convergence(l, val)
-        log_state(l, state)
         log_method(l, M)
+        log_last(l, α)
         verbose && display(l)
     end
+    α
 
-    val, stopped = assess_convergence(M, state, options) # udpate val flag
-    decide_convergence(M, F, state, options, val)
 end
 
 # thread verbose through
diff --git a/src/find_zeros.jl b/src/find_zeros.jl
index 0beb5c38..db217844 100644
--- a/src/find_zeros.jl
+++ b/src/find_zeros.jl
@@ -262,7 +262,7 @@ The algorithm is derived from one in a
 
 For example, this function (due to `@truculentmath`) is particularly tricky, as it is positive at every floating point number, but has two zeros (the vertical asymptote at `15//11` is only negative within adjacent floating point values):
 
-```jldoctest interval_root_finding
+```
 julia> using IntervalArithmetic, IntervalRootFinding, Roots
 
 julia> g(x) = x^2 + 1 +log(abs( 11*x-15 ))/99
@@ -278,7 +278,7 @@ julia> IntervalRootFinding.roots(g, -3..3, IntervalRootFinding.Bisection)
 
 A less extreme usage might be the following, where `unique` indicates Bisection could be useful and indeed `find_zeros` will identify these values:
 
-```jldoctest interval_root_finding
+```
 julia> g(x) = exp(x) - x^5
 g (generic function with 1 method)
 
diff --git a/src/order0.jl b/src/order0.jl
index 32d2a3c8..4133d978 100644
--- a/src/order0.jl
+++ b/src/order0.jl
@@ -67,9 +67,7 @@ function solve!(𝐙::ZeroProblemIterator{𝐌,𝐍}; verbose=false) where {𝐌
     M, N, F, state, options, l = 𝐙.M, 𝐙.N, 𝐙.F, 𝐙.state, 𝐙.options, 𝐙.logger
 
     incfn(l, 2)
-    log_step(l, M, state, :init)
-    log_method(l, M)
-    log_nmethod(l, N)
+    log_step(l, M, state; init=true)
 
     quad_ctr = 0
     flag = :not_converged
@@ -90,15 +88,16 @@ function solve!(𝐙::ZeroProblemIterator{𝐌,𝐍}; verbose=false) where {𝐌
             state = state0
             break
         elseif sign(state0.fxn0) * sign(state0.fxn1) < 0
+            log_step(l, M, state0)
             !isa(l, NullTracks) && log_message(
                 l,
-                "Used bracketing method $N on  [$(state0.xn0),$(state0.xn1)], those steps not recorded",
+                "Used bracketing method $N on  [$(min(state0.xn0,state0.xn1)),$(max(state0.xn0,state0.xn1))]",
             )
 
             Fₙ = Callable_Function(N, F)
             stateₙ = init_state(N, state0, Fₙ) # save function calls by using state0 values
             optionsₙ = init_options(N, stateₙ)
-            α = solve!(init(N, Fₙ, stateₙ, optionsₙ))
+            α = solve!(init(N, Fₙ, stateₙ, optionsₙ, l))
             break
         end
 
@@ -122,18 +121,18 @@ function solve!(𝐙::ZeroProblemIterator{𝐌,𝐍}; verbose=false) where {𝐌
 
         # a sign change after shortening?
         if sign(state.fxn1) * sign(state0.fxn1) < 0
+            log_step(l, M, state)
             a, b = state.xn1, state0.xn1
             fa, fb = state.fxn1, state0.fxn1
-
             !isa(l, NullTracks) && log_message(
                 l,
-                "Used bracketing method $N on  [$a,$b], those steps not recorded",
+                "Used bracketing method $N on  [$a,$b]",
             )
 
             Fₙ = Callable_Function(N, F)
             stateₙ = init_state(N, Fₙ, a, b, fa, fb)
             optionsₙ = init_options(N, stateₙ)
-            α = solve!(init(N, Fₙ, stateₙ, optionsₙ))
+            α = solve!(init(N, Fₙ, stateₙ, optionsₙ, l))
             break
         end
 
@@ -180,10 +179,14 @@ function solve!(𝐙::ZeroProblemIterator{𝐌,𝐍}; verbose=false) where {𝐌
         log_step(l, M, state)
     end
 
-    log_state(l, state)
+
+    val, stopped = assess_convergence(M, state, options)
+    log_convergence(l, val)
+    log_method(l, M)
+    log_nmethod(l, N)
+    log_last(l, α)
     verbose && display(l)
 
-    flag, converged = assess_convergence(M, state, options)
     isnan(α) ? decide_convergence(M, F, state, options, flag) : α
 end
 


Test patch:
diff --git a/test/test_find_zero.jl b/test/test_find_zero.jl
index fe9b070f..6f40c8e5 100644
--- a/test/test_find_zero.jl
+++ b/test/test_find_zero.jl
@@ -68,7 +68,7 @@ struct Order3_Test <: Roots.AbstractSecant end
     M = Roots.Bisection()
     tracks = Roots.Tracks(Float64, Float64)
     @inferred(find_zero(fn, (a, b), M, tracks=tracks, xatol=h, xrtol=0.0))
-    u, v = tracks.xs[(end - 1):end]
+    u, v = tracks.abₛ[end]
     @test h >= abs(u - v) >= h / 2
 
     ## test of strict