Differentiate find_zero and fzero (#148)

* find_zero(s) specializes on fn call f; fzero(s) does not; remove _find_zero; update README, docstrings

NEWS.md

@@ -1,3 +1,19 @@
+CHANGES in v0.7.4
+
+* set find_zero(s) to specialize on the function, fzero(s) to not. (#148)
+* adjust Steffensen method logic to take secant step or steffensen
+  step, rather than modified steffensen step. Seems to improve
+  robustness. (#147)
+* add Schroder method (order 2 for multiplicity with derivative), King (1B)
+  (superlinear for multiplicity, no derivative), Esser (2B) (order 2
+  for multipicity, no derivative) (#143, #147)
+* close issue #143 by allowing fns to Newton, Halley to compute f, f/fp, fp/fpp
+* add `newton` function to simple.jl
+* change find_zeros to identify zeros on [a,b], not (a,b). Closes #141.
+* bug fix: issue with quad step after a truncated M-step in find_zero(M,N,...)
+* bug fix: verbose argument for Bisection method (#139)
+* bug fix: unintentional widening of types in intial secant step (#139)
+
 CHANGES in v0.7.3
 
 * fix bug with find_zeros and Float32

README.md

@@ -1,5 +1,5 @@
 [![Roots](http://pkg.julialang.org/badges/Roots_0.6.svg)](http://pkg.julialang.org/?pkg=Roots&ver=0.6)
-[![Roots](http://pkg.julialang.org/badges/Roots_0.7.svg)](http://pkg.julialang.org/?pkg=Roots&ver=0.7)  
+[![Roots](http://pkg.julialang.org/badges/Roots_0.7.svg)](http://pkg.julialang.org/?pkg=Roots&ver=0.7)
 Linux: [![Build Status](https://travis-ci.org/JuliaMath/Roots.jl.svg?branch=master)](https://travis-ci.org/JuliaMath/Roots.jl)
 Windows: [![Build status](https://ci.appveyor.com/api/projects/status/goteuptn5kypafyl?svg=true)](https://ci.appveyor.com/project/jverzani/roots-jl)
 
@@ -13,15 +13,11 @@ specification of a method. These include:
 
 * Bisection-like algorithms. For functions where a bracketing interval
   is known (one where `f(a)` and `f(b)` have alternate signs), the
-  `Bisection` method can be specified. For most floating point number types, bisection occurs
-  in a manner exploiting floating point storage conventions. For
-  others, an algorithm of Alefeld, Potra, and Shi is used. These
-  methods are guaranteed to converge.
+  `Bisection` method can be specified. For most floating point number
+  types, bisection occurs in a manner exploiting floating point
+  storage conventions. For others, an algorithm of Alefeld, Potra, and
+  Shi is used. These methods are guaranteed to converge.
 
-  For typically faster convergence -- though not guaranteed -- the
-  `FalsePosition` method can be specified. This method has one of 12
-  implementations for a modified secant method to potentially
-  accelerate convergence.
 
 * Several derivative-free methods are implemented. These are specified
   through the methods `Order0`, `Order1` (the secant method), `Order2`
@@ -30,15 +26,19 @@ specification of a method. These include:
   method is the default, and the most robust, but may take many more
   function calls to converge. The higher order methods promise higher
   order (faster) convergence, though don't always yield results with
-  fewer function calls than `Order1` or `Order2`.
+  fewer function calls than `Order1` or `Order2`. The methods
+  `Roots.Order1B` and `Roots.Order2B` are superlinear and quadratically converging
+  methods independent of the multiplicity of the zero.
 
 
-* There are two historic methods that require a derivative or two:
-  `Roots.Newton` and `Roots.Halley`. (Neither is exported.)
+* There are historic methods that require a derivative or two:
+  `Roots.Newton` and `Roots.Halley`.  `Roots.Schroder` provides a
+  quadratic method, like Newton's method, which is independent of the
+  multiplicity of the zero.
 
 Each method's documentation has additional detail.
 
-Some examples: 
+Some examples:
 
 
 ```julia
@@ -116,26 +116,27 @@ Now we have:
 ```
 f(x) = x^3 - 2x - 5
 x0 = 2
-find_zero((f,D(f)) x0, Roots.Newton())   # 2.0945514815423265
+find_zero((f,D(f)), x0, Roots.Newton())   # 2.0945514815423265
 ```
 
 Automatic derivatives allow for easy solutions to finding critical
 points of a function.
 
 ```julia
 ## mean
+using Statistics
 as = rand(5)
-function M(x) 
+function M(x)
   sum([(x-a)^2 for a in as])
 end
-fzero(D(M), .5) - mean(as)	  # 0.0
+find_zero(D(M), .5) - mean(as)	  # 0.0
 
 ## median
-function m(x) 
+function m(x)
   sum([abs(x-a) for a in as])
 
 end
-fzero(D(m), 0, 1)  - median(as)	# 0.0
+find_zero(D(m), (0, 1)) - median(as)	# 0.0
 ```
 
 ### Multiple zeros
@@ -179,23 +180,22 @@ accepts keyword arguments `atol`, `rtol`, `xatol`, and `xrtol`.
 The `Bisection` and `Roots.A42` methods are guaranteed to converge
 even if the tolerances are set to zero, so these are the
 defaults. Non-zero values for `xatol` and `xrtol` can be specified to
-reduce the number of function calls when lower precision is required. 
+reduce the number of function calls when lower precision is required.
 
 
 ## An alternate interface
 
-For MATLAB users, this functionality is provided by the `fzero`
-function. `Roots` also provides this alternative interface:
+This functionality is provided by the `fzero` function, familiar to
+MATLAB users. `Roots` also provides this alternative interface:
 
 
-* `fzero(f, a::Real, b::Real)` and `fzero(f,
-  bracket::Vector)` call the `find_zero` algorithm with the
-  `Bisection` method.
-
-* `fzero(f, x0::Real; order::Int=0)` calls a
-  derivative-free method. with the order specified matching one of
+* `fzero(f, x0::Real; order=0)` calls a
+  derivative-free method. with the order specifying one of
   `Order0`, `Order1`, etc.
-
+
+* `fzero(f, a::Real, b::Real)` calls the `find_zero` algorithm with the
+  `Bisection` method.
+
 * `fzeros(f, a::Real, b::Real)` will call `find_zeros`.
 
 
@@ -207,16 +207,22 @@ f(x) = exp(x) - x^4
 ## bracketing
 fzero(f, 8, 9)		          # 8.613169456441398
 fzero(f, -10, 0)		      # -0.8155534188089606
-fzeros(f, -10, 10)            # -0.815553, 1.42961  and 8.61317 
+fzeros(f, -10, 10)            # -0.815553, 1.42961  and 8.61317
 
 ## use a derivative free method
 fzero(f, 3)			          # 1.4296118247255558
 
 ## use a different order
-fzero(sin, 3, order=16)		  # 3.141592653589793
+fzero(sin, big(3), order=16)  # 3.141592653589793...
 ```
 
+### Technical difference between find_zero and fzero
 
+The `fzero` function is not identical to `find_zero`. When a function, `f`,
+is passed to `find_zero` the code is specialized to the function `f`
+which means the first use of `f` will be slower due to compilation,
+but subsequent uses will be faster. For `fzero`, the code is not
+specialized to the function `f`, so the story is reversed.
 
 
 

src/derivative_free.jl

@@ -71,7 +71,7 @@ function find_zero(fs, x0, method::Order0;
     M = Order1()
     N = AlefeldPotraShi()
 
-    _find_zero(callable_function(fs), x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
+    find_zero(callable_function(fs), x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
 end
 
 ##################################################

src/find_zero.jl

@@ -243,22 +243,26 @@ end
 ### Functions
 abstract type CallableFunction end
 
-struct DerivativeFree <: CallableFunction
-    f
+struct FnWrapper <: CallableFunction
+   f
 end
 
-struct FirstDerivative <: CallableFunction
-    f
-    fp
+struct DerivativeFree{F} <: CallableFunction
+    f::F
 end
 
-struct SecondDerivative <: CallableFunction
-    f
-    fp
-    fpp
+struct FirstDerivative{F,FP} <: CallableFunction
+    f::F
+    fp::FP
 end
 
+struct SecondDerivative{F,FP,FPP} <: CallableFunction
+    f::F
+    fp::FP
+    fpp::FPP
+end
 
+(F::FnWrapper)(x::Number) = first(F.f(x))
 (F::DerivativeFree)(x::Number) = first(F.f(x))
 (F::FirstDerivative)(x::Number) = first(F.f(x))
 (F::SecondDerivative)(x::Number) = first(F.f(x))
@@ -286,24 +290,19 @@ end
 fΔxΔΔx(F, x) = F(x)
 
 
-## It is faster the first time a function is used if we do not
-## parameterize this. (As this requires less compilation) It is slower
-## the second time a function is used. This seems like the proper
-## tradeoff.
-##
-## We can override this with, say:
-## import FunctionWrappers
-## import FunctionWrappers: FunctionWrapper
-## struct CallbackF64F64
-## f::FunctionWrapper{Float64,Tuple{Float64}}
-## end
-## (cb::CallbackF64F64)(v) = cb.f(v)
-## Roots.callable_function(f::CallbackF64F64) = f
-
+# allows override for function, if desired
+# the default for this specializes on the function passed
+# in. When specialization occurs there is overhead due to compilation
+# costs which can be amortized over subsequent calls to `find_zero`.
+# However, if a function is only used once, then using `fzero` will be
+# faster. (The difference is clear between `@time` and `@btime` to measure
+# execution time)
+#                 first call    subsequent calls  default for
+# specialize       slower         faster          find_zero
+# no specialize    faster         slower          fzero
+#
 callable_function(fs::Any) = _callable_function(fs)
-
-## Default does not specialize on function
-function _callable_function(@nospecialize(fs))
+function _callable_function(fs)
     if isa(fs, Tuple)
         length(fs)==1 && return DerivativeFree(fs[1])
         length(fs)==2 && return FirstDerivative(fs[1],fs[2])
@@ -524,23 +523,6 @@ function find_zero(fs, x0, method::AbstractUnivariateZeroMethod,
                    kwargs...)
 
     F = callable_function(fs)
-
-    _find_zero(F, x0, method, N; tracks=tracks, verbose=verbose, kwargs...)
-end
-
-## This allows for bypassing of `callable_function`.  With
-## callable_function there is no specialization on the function. This
-## makes the first call faster, but subsequent calls slower, due to
-## type instability. Without callable_function the first usage is
-## slower, though this can be effectively reduced using
-## `FunctionWrappers`.
-function _find_zero(F, x0, method::AbstractUnivariateZeroMethod,
-                    N::Union{Nothing, AbstractBracketing}=nothing;
-                    tracks::AbstractTracks=NullTracks(),
-                    verbose=false,
-                    kwargs...)
-
-
     state = init_state(method, F, x0)
     options = init_options(method, state; kwargs...)