Roots change callable_function (#227)

* replace callable_function
* clean up interface (find_zero, solve(ZeroProblem), solve!(ZeroProblemIterator)...)

Project.toml

@@ -1,6 +1,6 @@
 name = "Roots"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "1.0.11"
+version = "1.0.12"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"

docs/src/roots.md

@@ -15,6 +15,7 @@ julia> using Roots
 
 julia> using Plots, ForwardDiff
 
+
 ```
 
 ## Bracketing
@@ -622,7 +623,7 @@ julia> u=1/3; for i in 1:10 (global  u=prevfloat(u);push!(ns, u)) end
 julia> sort!(ns);
 
 julia> maximum(abs.(f.(ns) - f1.(ns)))
-1.5543122344752192e-15
+1.887379141862766e-15
 
 ```
 
@@ -632,16 +633,16 @@ is like $10^{15}$. This is roughly as expected, where even one
 addition may introduce a relative error as big as $2\cdot 10^{-16}$ and here
 there are several such.
 
-!!! Note:
+!!! note
     (These values are subject to the vagaries of floating point evaluation, so may differ depending on the underlying computer architecture.)
 
 Generally this variation is not even a thought, as the differences are generally
 negligible, but when we want to identify if a value is zero, these
 small differences might matter. Here we look at the signs of the
-function values:
+function values for a run of the above:
 
 
-```jldoctest roots
+```
 julia> fs = sign.(f.(ns));
 
 julia> f1s = sign.(f1.(ns));
@@ -677,7 +678,7 @@ Parsing this shows a few surprises. First, there are two zeros of
 floating point value of `1/3` and the next largest floating point
 number.
 
-```jldoctest roots
+```
 julia> findall(iszero, fs)
 2-element Vector{Int64}:
  11
@@ -691,7 +692,7 @@ point value for `1/3` but rather 10 floating point numbers
 away.
 
 
-```jldoctest roots
+```
 julia> findall(iszero, f1s)
 1-element Vector{Int64}:
  21
@@ -701,7 +702,7 @@ julia> findall(iszero, f1s)
 
 Further, there are several sign changes of the function values for `f1s`:
 
-```jldoctest roots
+```
 julia> findall(!iszero,diff(sign.(f1s)))
 9-element Vector{Int64}:
   2

src/bracketing.jl

@@ -8,13 +8,9 @@ Consider a different bracket or try fzero(f, c) with an initial guess c.
 """
 
 
-# pushed this into both __middle and bisection(f,a,b)
-@deprecate bisection64(f, a, b) Roots.bisection(f, a, b)
-
-
-
 
 abstract type AbstractBisection <: AbstractBracketing end
+fn_argout(::AbstractBracketing) = 1
 
 """
 
@@ -76,7 +72,6 @@ end
 # In init_state the bracketing interval is left as (a,b) with
 # a,b both finite and of the same sign
 function init_state(method::AbstractBisection, fs, x)
-    length(x) > 1 || throw(ArgumentError(bracketing_error))
 
     x0, x1 = adjust_bracket(x) # now finite, right order
     fx0, fx1 = promote(fs(x0), fs(x1))
@@ -325,32 +320,35 @@ end
 ## for zero tolerance, we have either BisectionExact or A42 methods
 ## for non-zero tolerances, we have use thegeneral Bisection method
 function find_zero(fs, x0, method::M;
+                   p=nothing,
                    tracks = NullTracks(),
                    verbose=false,
                    kwargs...) where {M <: Union{Bisection}}
 
-    x = adjust_bracket(x0)
-    T = eltype(x[1])
-    F = callable_function(fs)
-    state = init_state(method, F, x)
+#    x = adjust_bracket(x0)
+#    T = eltype(x[1])
+
+    F = Callable_Function(method, fs, p)
+    state = init_state(method, F, x0)
+    x = (state.xn0, state.xn1)
     options = init_options(method, state; kwargs...)
+    l = Tracks(verbose, tracks, state)
 
     # check if tolerances are exactly 0
     iszero_tol = iszero(options.xabstol) && iszero(options.xreltol) && iszero(options.abstol) && iszero(options.reltol)
 
-    l = (verbose && isa(tracks, NullTracks)) ? Tracks(eltype(state.xn1)[], eltype(state.fxn1)[]) : tracks
 
     if iszero_tol
-        if T <: FloatNN
+        if eltype(state.xn1) <: FloatNN
             return find_zero(F, x, BisectionExact(); tracks=l, verbose=verbose, kwargs...)
         else
             return find_zero(F, x, A42(); tracks=l, verbose=verbose, kwargs...)
         end
     end
 
-    find_zero(method, F, state, options, l)
-
-    verbose && show_trace(method, nothing, state, l)
+    ZPI = ZeroProblemIterator(method, F, state, options, l)
+    solve!(ZPI)
+    verbose && tracks(ZPI)
 
     state.xstar
 
@@ -1010,14 +1008,15 @@ With the default tolerances, one  of these should be the case: `exact`  is `true
 function find_bracket(fs, x0, method::M=A42(); kwargs...) where {M <: Union{AbstractAlefeldPotraShi, BisectionExact}}
     x = adjust_bracket(x0)
     T = eltype(x[1])
-    F = callable_function(fs)
+    F = Callable_Function(method, fs) #callable_function(fs)
     state = init_state(method, F, x)
     options = init_options(method, state; kwargs...)
 
     # check if tolerances are exactly 0
     iszero_tol = iszero(options.xabstol) && iszero(options.xreltol) && iszero(options.abstol) && iszero(options.reltol)
 
-    find_zero(method, F, state, options, NullTracks())
+    ZPI = ZeroProblemIterator(method, F, state, options, NullTracks())
+    solve!(ZPI)
 
     (xstar=state.xstar, bracket=(state.xn0, state.xn1), exact=iszero(state.fxstar))
 

src/derivative_free.jl

@@ -43,10 +43,8 @@ end
 # check if we should guard against step for method M; call N if yes, P if not
 function update_state_guarded(M::AbstractSecant,N::AbstractUnivariateZeroMethod, P::AbstractUnivariateZeroMethod, fs, o, options)
     if do_guarded_step(M, o)
-        #@debug "do secant step"
         return update_state(N, fs, o, options)
     else
-        #@debug "do $N step"
         update_state(P, fs, o, options)
     end
 end
@@ -60,8 +58,8 @@ function init_state(method::AbstractSecant, fs, x::Number)
 end
 
 function init_state(method::AbstractSecant, fs, x)
-    x0, x1 = x₀x₁(x)
-    fx0, fx1 = fs(x0), fs(x1)
+    x0, x1 = x₀x₁(x);
+    fx0, fx1 = first(fs(x0)), first(fs(x1))
     state = UnivariateZeroState(x1, x0, zero(x1)/zero(x1)*oneunit(x1), typeof(x1)[],
                                 fx1, fx0, fx1, typeof(fx1)[],
                                 0, 2,
@@ -114,13 +112,14 @@ superlinear, and relatively robust to non-reasonable starting points.
 struct Order0 <: AbstractSecant end
 
 function find_zero(fs, x0, method::Order0;
+                   p=nothing,
                    tracks::AbstractTracks=NullTracks(),
                    verbose=false,
                    kwargs...)
     M = Order1()
     N = AlefeldPotraShi()
-
-    find_zero(callable_function(fs), x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
+    F = Callable_Function(M, fs, p)
+    find_zero(F, x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
 end
 
 ##################################################
@@ -460,15 +459,14 @@ function update_state(M::Union{Order5, KumarSinghAkanksha}, fs, o::UnivariateZer
     nothing
 end
 
-
-
-## If we have a derivative, we have this. (Deprecate?)
-function update_state(method::Order5, fs::Union{FirstDerivative,SecondDerivative},
+struct Order5Derivative <: AbstractSecant end
+fn_argout(::Order5Derivative) = 2
+function update_state(method::Order5Derivative, f,
                       o::UnivariateZeroState{T,S}, options)  where {T, S}
 
 
     xn, fxn = o.xn1, o.fxn1
-    a::T, b::S = fΔx(fs, xn)
+    a::T, b::S = f(xn)
     fpxn = a/b
     incfn(o)
 
@@ -478,7 +476,7 @@ function update_state(method::Order5, fs::Union{FirstDerivative,SecondDerivative
     end
 
     yn::T = xn - fxn / fpxn
-    fyn::S, Δyn::T = fΔx(fs, yn)
+    fyn::S, Δyn::T = f(yn)
     fpyn = fyn / Δyn
     incfn(o, 2)
 
@@ -491,12 +489,12 @@ function update_state(method::Order5, fs::Union{FirstDerivative,SecondDerivative
 
 
     zn::T = xn  - (fxn + fyn) / fpxn
-    fzn::S = fs(zn)
-    incfn(o, 1)
+    fzn::S, _ = f(zn)
+    incfn(o, 2)
 
     xn1 = zn - fzn / fpyn
-    fxn1 = fs(xn1)
-    incfn(o, 1)
+    fxn1,_ = f(xn1)
+    incfn(o, 2)
 
     o.xn0, o.xn1 = xn, xn1
     o.fxn0, o.fxn1 = fxn, fxn1
@@ -610,7 +608,7 @@ this method generally isn't faster (fewer function calls/steps) over
 other methods when using `Float64` values, but may be useful for
 solving over `BigFloat`.
 
-The error, `eᵢ = xᵢ - α`, is expressed as `eᵢ₊₁ = K⋅eᵢ¹⁶` for an explicit `K` 
+The error, `eᵢ = xᵢ - α`, is expressed as `eᵢ₊₁ = K⋅eᵢ¹⁶` for an explicit `K`
 in equation (50) of the paper.
 
 """