update tests; bug fixes; drop pre-1.0 support (#153)

.travis.yml

@@ -3,8 +3,8 @@ os:
   - osx
   - linux
 julia:
-  - 0.6
-  - 0.7
+  - 1.0
+  - 1.1
   - nightly
 matrix:
   allow_failures:

REQUIRE

@@ -1,2 +1 @@
-julia 0.6.0
-Compat 0.59.0
+julia 1.0.0

appveyor.yml

@@ -1,7 +1,7 @@
 environment:
   matrix:
-  - julia_version: 0.6
-  - julia_version: 0.7
+  - julia_version: 1.0
+  - julia_version: 1.1
   - julia_version: latest
 
 platform:

src/Roots.jl

@@ -1,13 +1,6 @@
-VERSION < v"0.7.0-beta2.199" && __precompile__()
 module Roots
 
-
-if VERSION >= v"0.7-"
-    using Printf
-end
-
-using Compat: @nospecialize, lastindex, range, Nothing
-
+using Printf
 
 export fzero,
        fzeros,

src/find_zero.jl

@@ -178,11 +178,8 @@ function init_options(M::AbstractUnivariateZeroMethod,
                       kwargs...
                       ) where {T, S}
 
-    if VERSION < v"0.7.0"
-        d = Dict(kwargs)
-    else
-        d = kwargs
-    end
+    d = kwargs
+
     defs = default_tolerances(M, T, S)
     options = UnivariateZeroOptions(get(d, :xatol, get(d, :xabstol, defs[1])),
                                     get(d, :xrtol, get(d, :xreltol, defs[2])),

src/fzero.jl

@@ -172,5 +172,5 @@ Dispatches to `find_zeros(f, a, b; kwargs...)`.
 function fzeros(f, a::Number, b::Number; kwargs...)
     find_zeros(FnWrapper(f), float(a), float(b); kwargs...)
 end
-fzeros(f, bracket::Vector{T}; kwargs...) where {T <: Number} = fzeros(f, a, b; kwargs...)
-fzeros(f, bracket::Tuple{T,S}; kwargs...) where {T <: Number, S<:Number} = fzeros(f, a, b; kwargs...)
+fzeros(f, bracket::Vector{T}; kwargs...) where {T <: Number} = fzeros(f, bracket[1], bracket[2]; kwargs...)
+fzeros(f, bracket::Tuple{T,S}; kwargs...) where {T <: Number, S<:Number} = fzeros(f, bracket[1], bracket[2]; kwargs...)

src/newton.jl

@@ -278,6 +278,12 @@ function update_state(method::Schroder, fs, o::UnivariateZeroState{T,S}, options
     r1, r2 = o.m[1], o.m[2]
 
     delta =  r2 / (r2 - r1) * r1  # m * r1
+
+    if isissue(delta)
+        o.stopped=true
+        return
+    end
+
     xn1::T = xn - delta
 
     tmp = fΔxΔΔx(fs, xn1)

test/.gitignore

@@ -0,0 +1 @@
+*.json

test/REQUIRE

@@ -1 +1,2 @@
-SpecialFunctions 0.1.1
+SpecialFunctions 0.1.1
+JSON

test/runbenchmarks.jl

@@ -0,0 +1,126 @@
+# v1.0 only
+# Import dependency.
+using BenchmarkTools, Roots, Statistics
+
+#Create benchmark group and benchmarks
+benchmarks = BenchmarkGroup()
+
+#Put in specific benchmarks
+
+bracket_methods = (bisection=Roots.Bisection(), a42=Roots.A42(),
+                   aps=Roots.AlefeldPotraShi())
+derivative_free_methods =  (o0=Roots.Order0(), o1=Roots.Order1(), o1b=Roots.Order1B(),
+                            o2=Roots.Order2(), o2b=Roots.Order2B(),
+                            o5=Roots.Order5(), o8=Roots.Order8(), o16=Roots.Order16()
+)
+
+# collection of doable problems
+problems = Dict("f1" => (x -> sin(x), 3.0,  (3.0, 4.0)),
+                "f2" => (x -> x^5 - x - 1, 1.0, (0.5, 5.0)),
+                "f3" => (x -> exp(x) - x^4, 7.0, (5.0, 20.0)),
+                "f4" => (x -> cos(x) - x/2, pi/4, (0.0, pi/2)),
+                "f5" => (x -> x^2 - exp(x) - 3x + 2, -0.5, (-1.0, 1.0)),
+                "f6" => (x -> x^2 - exp(x) - 3x + 2, 2.0, (0.0, 3.0)),
+                "f7" => (x -> tanh(x) - tan(x), 7.6, (4.0, 8.0)),
+"f8" => (x -> exp(-x^2 + x + 2) - cos(x) + x^3 + 1, -0.5, (-2.0, 1.0)),
+"f9" => (x -> log(x) + sqrt(x) - 5, 7, (7.0, 10.0)),
+"f10" => (x -> log(x) + sqrt(x) - 5, 20, (7.0, 10.0))
+                )
+
+function run_bracket(problems, Ms)
+
+    for (nm, prob) in problems
+        fn, x0, ab = prob
+        for (mnm, M) in zip(fieldnames(typeof(Ms)), Ms)
+            find_zero(fn, ab, M)
+        end
+    end
+
+end
+
+function run_bracketing(problems, Ms)
+    rts = Float64[]
+    for (nm, prob) in problems
+        fn, x0, ab = prob
+        for M in Ms
+            rt = find_zero(fn, ab)
+            push!(rts, rt)
+        end
+    end
+    rts
+
+end
+
+function run_derivative_free(problems, Ms)
+    rts = Float64[]
+    for (nm, prob) in problems
+        fn, x0, ab = prob
+        for M in Ms
+            rt = find_zero(fn, x0, M)
+            push!(rts, rt)
+        end
+    end
+    rts
+
+end
+
+function run_simple(problems)
+    rts = Float64[]
+    for (nm, prob) in problems
+        fn, x0, ab = prob
+        push!(rts, Roots.bisection(fn, ab[1], ab[2]))
+        push!(rts, Roots.bisection(fn, ab[1], ab[2], xatol=1e-6))
+        push!(rts, Roots.secant_method(fn, x0))
+    end
+    rts
+end
+
+
+benchmarks = BenchmarkGroup()
+
+benchmarks["bracketing"] = @benchmarkable run_bracketing($problems, $bracket_methods)
+benchmarks["derivative_free"] = @benchmarkable run_derivative_free($problems, $derivative_free_methods)
+benchmarks["simple"] = @benchmarkable run_simple($problems)
+
+
+
+
+for (nm, prob) in problems
+    fn, x0, ab = prob
+    @assert fn(ab[1]) * fn(ab[2]) < 0
+
+    Ms =  bracket_methods
+    for (mnm, M) in zip(fieldnames(typeof(Ms)), Ms)
+        benchmarks[nm * "-" * string(mnm)] = @benchmarkable find_zero($fn, $ab, $M)
+    end
+
+    Ms =  derivative_free_methods
+    for (mnm, M) in zip(fieldnames(typeof(Ms)), Ms)
+        benchmarks[nm * "-" * string(mnm)] = @benchmarkable find_zero($fn, $x0, $M)
+    end
+
+    # simple methods
+    u,v = ab
+    benchmarks[nm * "-bisection"] = @benchmarkable Roots.bisection($fn, $u, $v)
+    benchmarks[nm * "-bisection-atol"] = @benchmarkable Roots.bisection($fn, $u, $v, xatol=1e-6)
+        benchmarks[nm * "-secant"] = @benchmarkable Roots.secant_method($fn, $x0)
+end
+
+results = run(benchmarks) # Get results.
+results = median(results) # Condense to median.
+
+nm = "benchmarks.json"
+fname = joinpath(@__DIR__,  nm)
+if isinteractive()
+    println("""
+To save results, manually call in the REPL: BenchmarkTools.save("benchmarks.json", results)
+""")
+end
+
+#Compare to old results
+try
+    oldresults= BenchmarkTools.load(fname)[1]
+    judge(oldresults, results)
+catch err
+    error("Couldn't load file- make sure that you've previously saved results.", err.prefix)
+end

test/runtests.jl

@@ -1,19 +1,16 @@
 using Roots
-using Compat.Test
+using Test
 import SpecialFunctions.erf
 
 include("./test_find_zero.jl")
-include("./test_fzero.jl")
+include("./test_bracketing.jl")
+include("./test_derivative_free.jl")
+include("./test_simple.jl")
 include("./test_find_zeros.jl")
+include("./test_fzero.jl")
 include("./test_newton.jl")
 include("./test_simple.jl")
-include("./RootTesting.jl")
 
 #include("./test_composable.jl")
-
-#run_benchmark_tests()
-
-#include("./test_fzero3.jl")
-#run_robustness_test()
-
-#include("./test_derivative_free.jl")
+#include("./runbenchmarks.jl")
+#include("./test_derivative_free_interactive.jl")

test/RootTesting.jl → test/test_bracketing.jl

@@ -1,9 +1,15 @@
 using Roots
+using Test
+using Printf
+
+## testing bracketing methods
+## Orignially by John Travers
+#
+#
+## This set of tests is very useful for benchmarking the number of function
+## calls, failures, and max errors for the various bracketing methods.
+## Table 1 from TOMS748 by Alefeld, Potra, Shi
 
-import Base: show
-if VERSION >= v"0.7.0-"
-    using Printf
-end
 
 mutable struct Func
     name :: Symbol
@@ -201,45 +207,91 @@ function run_tests(method; verbose=false, trace=false, name=nothing, abandon=fal
     results
 end
 
-Ms = vcat([Roots.FalsePosition(i) for i in 1:12], Roots.A42(), Roots.AlefeldPotraShi(), Roots.Brent(), Roots.Bisection())
-results = [run_tests((f,b) -> find_zero(f, b, M), name="$M") for M in Ms]
-@test all([length(result.failures) <= 10 for result in results])
+avg(x) = sum(x)/length(x)
 
-results = [run_tests((f,b) -> find_zero(f, big.(b), M), name="$M") for M in Ms]
-@test all([length(result.failures) == 0 for result in results[13:end]])
+@testset "bracketing methods" begin
 
+    ## Test for failures, ideally all of these would be 0
+    ## test for residual, ideally small
+    ## test for evaluation counts, ideally not so low for these problems
 
 
+    ## exact_bracket
+    Ms = [Roots.A42(), Roots.AlefeldPotraShi(), Roots.Bisection()]
+    results = [run_tests((f,b) -> find_zero(f, b, M), name="$M") for M in Ms]
+    maxfailures = maximum([length(result.failures) for result in results])
+    maxresidual = maximum([result.maxresidual for result in results])
+    cnts = [result.evalcount for result in results]
+    @test maxfailures == 0
+    @test maxresidual <= 5e-15
+    @test avg(cnts) <= 4700
 
-## run but not in general testing suite
-function run_benchmark_tests()
+    ## brent has some failures
+    Ms = [Roots.Brent()]
+    results = [run_tests((f,b) -> find_zero(f, b, M), name="$M") for M in Ms]
+    maxfailures = maximum([length(result.failures) for result in results])
+    maxresidual = maximum([result.maxresidual for result in results])
+    cnts = [result.evalcount for result in results]
+    @test maxfailures <= 4
+    @test maxresidual <= 1e-13
+    @test avg(cnts) <= 2600
 
-    @printf "%s\n" run_tests((f,b) -> find_zero(f, b, Roots.A42()), name="a42")
-    @printf "%s\n" run_tests((f,b) -> find_zero(f, b, Bisection()), name="Bisection")
-    @printf "%s\n" run_tests((f,b) -> find_zero(f, b, FalsePosition()), name="FalsePosition")
+    ## False position has failures, and larger residuals
+    Ms = [Roots.FalsePosition(i) for i in 1:12]
+    results = [run_tests((f,b) -> find_zero(f, b, M), name="$M") for M in Ms]
+    maxfailures = maximum([length(result.failures) for result in results])
+    maxresidual = maximum([result.maxresidual for result in results])
+    cnts = [result.evalcount for result in results]
+    @test maxfailures <= 5
+    @test maxresidual <= 5e-8
+    @test avg(cnts) <= 2500
 
 
-    for m in [Order0(), Order1(), Order2(), Roots.Order2B(), Order5(), Order8(), Order16()]
-        @printf "%s\n" run_tests((f, b) -> find_zero(f, mean(b), m), name="$m")
-    end
 
-    println("---- using BigFloat ----")
 
-    @printf "%s\n" run_tests((f,b) -> find_zero(f, big.(b), Bisection()), name="a42 (no bisection with Big values)")
-    @printf "%s\n" run_tests((f,b) -> find_zero(f, big.(b), FalsePosition()), name="FalsePosition")
+end
 
-    for m in [Order0(), Order1(), Order2(), Order5(), Order8(), Order16()]
-        @printf "%s\n" run_tests((f, b) -> find_zero(f, mean(big.(b)), m), name="$m/BigFloat")
+## Some tests for FalsePosition methods
+@testset "FalsePosition" begin
+    galadino_probs = [(x -> x^3 - 1, [.5, 1.5]),
+                      (x -> x^2 * (x^2/3 + sqrt(2) * sin(x)) - sqrt(3)/18, [.1, 1]),
+                      (x -> 11x^11 - 1, [0.1, 1]),
+                      (x ->  x^3 + 1, [-1.8, 0]),
+                      (x ->  x^3 - 2x - 5, [2.0, 3]),
+
+                      ((x,n=5)  -> 2x * exp(-n) + 1 - 2exp(-n*x) , [0,1]),
+                      ((x,n=10) -> 2x * exp(-n) + 1 - 2exp(-n*x) , [0,1]),
+    ((x,n=20) -> 2x * exp(-n) + 1 - 2exp(-n*x) , [0,1]),
+
+    ((x,n=5)  -> (1 + (1-n)^2) * x^2 - (1 - n*x)^2 , [0,1]),
+    ((x,n=10) -> (1 + (1-n)^2) * x^2 - (1 - n*x)^2 , [0,1]),
+    ((x,n=20) -> (1 + (1-n)^2) * x^2 - (1 - n*x)^2 , [0,1]),
+
+    ((x,n=5)  -> x^2 - (1-x)^n , [0,1]),
+    ((x,n=10) -> x^2 - (1-x)^n , [0,1]),
+    ((x,n=20) -> x^2 - (1-x)^n , [0,1]),
+
+    ((x,n=5)  -> (1 + (1-n)^4)*x - (1 - n*x)^4 , [0,1]),
+    ((x,n=10) -> (1 + (1-n)^4)*x - (1 - n*x)^4 , [0,1]),
+    ((x,n=20) -> (1 + (1-n)^4)*x - (1 - n*x)^4 , [0,1]),
+
+    ((x,n=5)  -> exp(-n*x)*(x-1) + x^n , [0,1]),
+    ((x,n=10) -> exp(-n*x)*(x-1) + x^n , [0,1]),
+    ((x,n=20) -> exp(-n*x)*(x-1) + x^n , [0,1]),
+
+    ((x,n=5)  -> x^2 + sin(x/n) - 1/4 , [0,1]),
+    ((x,n=10) -> x^2 + sin(x/n) - 1/4 , [0,1]),
+    ((x,n=20) -> x^2 + sin(x/n) - 1/4 , [0,1])
+    ]
+
+
+    for (fn_, ab) in galadino_probs
+        for M in (FalsePosition(i) for i in 1:12)
+            g = Cnt(fn_)
+            x0_ = find_zero(g, ab, M)
+            @test abs(fn_(x0_)) <= 1e-7
+            @test g.cnt <= 50
+        end
     end
 
-    @printf "%s\n" run_tests((f, b) -> Roots.secant_method(f, big.(b)), name="secant_method")
-    @printf "%s\n" run_tests((f, b) -> Roots.dfree(f, mean(big.(b))), name="dfree")
-#    using ForwardDiff
-#    D(f,n=1) = n > 1 ? D(D(f), n-1) : x -> ForwardDiff.derivative(f, float(x))
-#    @printf "%s\n" run_tests((f, b) -> find_zero(f, D(f), mean(big.(b)), Order5()), name="Order5/D;BigFloat")
-#    @printf "%s\n" run_tests((f, b) -> newton(f, D(f), mean(big.(b))), name="newton/BigFloat")
-#    @printf "%s\n" run_tests((f, b) -> halley(f, D(f), D(f,2), mean(big.(b))), name="halley/BigFloat")
 end
-
-
-nothing

test/test_composable.jl

@@ -3,7 +3,7 @@
 # that CI should run these
 
 using Roots
-using Compat.Test
+using Test
 using Unitful
 using SymEngine