diff --git a/Project.toml b/Project.toml index 8ab743c4..8f53ca03 100644 --- a/Project.toml +++ b/Project.toml @@ -18,11 +18,4 @@ InteractiveUtils = "1.10.0" IntervalArithmetic = "1.0.3" Reexport = "1" StaticArrays = "1" -julia = "1" - -[extras] -Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["Test", "Polynomials"] +julia = "1" \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 00000000..5cae9d08 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,22 @@ +[deps] +BranchAndPrune = "d3bc4f2e-91e6-11e9-365e-cd067da536ce" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" +IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +BranchAndPrune = "0.2.1" +ForwardDiff = "0.10, 1" +InteractiveUtils = "1.10.0" +IntervalArithmetic = "1.0.3" +Reexport = "1" +StaticArrays = "1" +Suppressor = "0.2.8" +julia = "1" diff --git a/test/findroots.jl b/test/findroots.jl deleted file mode 100644 index 602ff287..00000000 --- a/test/findroots.jl +++ /dev/null @@ -1,171 +0,0 @@ - -using IntervalArithmetic, IntervalRootFinding -using ForwardDiff -using Test - -const D = IntervalRootFinding.derivative - -include("wilkinson.jl") - - -function generate_wilkinson(n)#, T=BigFloat) # SLOW - - p = poly(collect(1:n)) - p_prime = polyder(p) - - coeffs = p.a #[convert(T, coeff) for coeff in p.a] - # coeffs_prime = [convert(T, coeff) for coeff in p_prime.a] - - function f(x) - total = coeffs[1] - for i in 2:length(coeffs) - total += coeffs[i]*x^(i-1) - end - total - end - - f -end - - -setprecision(Interval, Float64) -float_pi = @interval(pi) - -setprecision(Interval, 10000) -big_pi = @interval(pi) -# Using precision "only" 256 leads to overestimation of the true roots for `cos` -# i.e the Newton method gives more accurate results! - -half_pi = big_pi / 2 -three_halves_pi = 3*big_pi/2 - - -# Format: (function, derivative, lower_bound, upper_bound, [true_roots]) -function_list = [ - (sin, cos, -5, 5, [ -big_pi, @interval(0), big_pi ] ) , - (cos, x->D(cos, x), -7.5, 7.5, [ -three_halves_pi, -half_pi, half_pi, three_halves_pi ] ), - (W₃, x->D(W₃, x), -10, 10, [ @interval(1), @interval(2), @interval(3) ] ), - (W₇, x->D(W₇, x), -10, 10, [ @interval(i) for i in 1:7 ] ), - (x->exp(x)-2, y->D(x->exp(x),y), -20, 20, [log(@biginterval(2))] ), - (x->asin(sin(x)) - 0.1, y->D(x->asin(sin(x)),y), 0, 1, [@biginterval(0.1)]) - ] - - -@testset "Testing root finding" begin - - for rounding_type in (:wide, :narrow) - @testset "Interval rounding: $rounding_type" begin - # setrounding(Interval, rounding_type) - - for prec in ( (BigFloat,53), (BigFloat,256), (Float64,64) ) - @testset "Precision: $prec" begin - setprecision(Interval, prec) - - for method in (IntervalRootFinding.newton, IntervalRootFinding.krawczyk) - @testset "Method $method" begin - - for func in function_list - - f, f_prime, a_lower, a_upper, true_roots = func - a = @interval(a_lower, a_upper) - - @testset "Function $f; interval $a" begin - - for autodiff in (false, true) - - @testset "With autodiff=$autodiff" begin - - if autodiff - rts = find_roots(f, a, method) - else - rts = find_roots(f, f_prime, a, method) - end - - @test length(rts) == length(true_roots) - - for i in 1:length(rts) - root = rts[i] - - @test isa(root, Root) - @test is_unique(root) - @test true_roots[i] ⊆ root.interval - end - end - end - end - end - end - end - end - end - end - end -end - - - -setprecision(Interval, Float64) - -@testset "find_roots tests" begin - f(x) = x^2 - 2 - - rts = IntervalRootFinding.newton(f, @interval(10, 11)) - @test length(rts) == 0 - - - rts = IntervalRootFinding.find_roots_midpoint(f, -5, 5) - @test length(rts) == 3 - @test length(rts[1]) == 2 - - rts = find_roots(f, -5, 5) - @test length(rts) == 2 - - setprecision(Interval, 256) - - for method in (IntervalRootFinding.newton, IntervalRootFinding.krawczyk) - new_roots = method(f, rts) - @test length(new_roots) == length(rts) - - for (root1, root2) in zip(new_roots, rts) - @test root1 ⊆ root2 - end - end -end - -@testset "Multiple roots" begin - setprecision(Interval, Float64) - let - f(x) = (x-1) * (x^2 - 2)^3 * (x^3 - 2)^4 - - rts = IntervalRootFinding.newton(f, -5..5.1, maxlevel=1000) - - @test length(rts) == 4 - @test rts[1].status == :unknown - @test rts[1].interval == Interval(-1.4142135623730954, -1.414213562373095) - @test rts[3].interval == Interval(1.259921049894873, 1.2599210498948734) - - end -end - -@testset "Iterated logistic function fixed points" begin - ∘(f, g) = x -> f(g(x)) - iterate(f, n) = n == 1 ? f : f ∘ iterate(f, n-1) - - f(x) = 4x*(1-x) # logistic map - - for n in 1:10 - g = x -> iterate(f, n)(x) - x # look for fixed points of f^n - - rts = IntervalRootFinding.newton(g, 0..1) - @test length(rts) == 2^n - - rts = IntervalRootFinding.krawczyk(g, 0..1) - @test length(rts) == 2^n - - end -end - - - -# Example of a function with a double root at 0 from Burden & Faires, 9th ed, p.84 -# exp(x) - x - 1 diff --git a/test/multidimensional.jl b/test/multidimensional.jl new file mode 100644 index 00000000..be0bf16c --- /dev/null +++ b/test/multidimensional.jl @@ -0,0 +1,74 @@ +@testset "2D roots" begin + f(x, y) = [x^2 + y^2 - 1, y - 2x] + f(X) = f(X...) + X = [interval(-6, 6), interval(-6, 6)] + + # Bisection + rts = roots(f, X ; contractor = Bisection, abstol = 1e-3) + exact_sol = [sqrt(1/5), 2sqrt(1/5)] + @test in_solution_set(exact_sol, root_region.(rts)) + @test in_solution_set(-exact_sol, root_region.(rts)) + + for contractor in newtonlike_methods + deriv = xx -> ForwardDiff.jacobian(f, xx) + test_newtonlike(f, deriv, X, contractor, 2) + rts = roots(f, X ; contractor) + @test in_solution_set(exact_sol, root_region.(rts)) + @test in_solution_set(-exact_sol, root_region.(rts)) + end + + # Infinite interval + X = [interval(-Inf, Inf), interval(-Inf, Inf)] + @suppress rts = roots(f, X ; contractor = Newton) + @test length(rts) == 2 +end + + +# From R docs: +# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden + +@testset "3D roots" begin + function g(x) + (x1, x2, x3) = x + SVector( x1^2 + x2^2 + x3^2 - 1, + x1^2 + x3^2 - 0.25, + x1^2 + x2^2 - 4x3 + ) + end + dg = xx -> ForwardDiff.jacobian(g, xx) + + X = interval(-5, 5) + XX = [X, X, X] + + for contractor in newtonlike_methods + rts = roots(g, XX ; contractor) + @test length(rts) == 4 + @test all(isunique, rts) + @test all(rts .== roots(g, XX ; contractor, derivative = dg)) + end +end + +@testset "10D roots" begin + include("../examples/10_dimensional.jl") + + for contractor in newtonlike_methods + rts = roots(f10d, X10d_large ; contractor, max_iteration = 50_000) + @test any(X -> all(in_interval.(sol10d, X.region)), rts) + + rts = roots(f10d, X10d_close ; contractor, max_iteration = 1_000_000) + @test length(rts) == 1 + @test isunique(rts[1]) + end +end + +@testset "Dimension mismatch" begin + f21(xy) = [xy[1]^2 - 2] + f23(xy) = [xy[1]^2 - 2, xy[2]^2 - 3, xy[1] + xy[2]] + + X = [interval(0, 5), interval(0, 5)] + + for contractor in newtonlike_methods + @test_throws DimensionMismatch roots(f21, X ; contractor) + @test_throws DimensionMismatch roots(f23, X ; contractor) + end +end diff --git a/test/newton1d.jl b/test/newton1d.jl deleted file mode 100644 index 2c11c0bb..00000000 --- a/test/newton1d.jl +++ /dev/null @@ -1,86 +0,0 @@ - -using IntervalArithmetic, IntervalRootFinding -using ForwardDiff -using Test - -const D = IntervalRootFinding.derivative - -setprecision(Interval, Float64) -float_pi = @interval(pi) - -setprecision(Interval, 10000) -big_pi = @interval(pi) -# Using precision "only" 256 leads to overestimation of the true roots for `cos` -# i.e the Newton method gives more accurate results! - -half_pi = big_pi / 2 -three_halves_pi = 3*big_pi/2 - -@testset "Testing newton1d" begin - - f(x) = exp(x^2) - cos(x) #double root - f′(x) = 2*x*exp(x^2) + sin(x) - f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24 #four unique roots - f1′(x) = 4x^3 - 30x^2 + 70x - 50 - f2(x) = 4567x^2 - 9134x + 4567 #double root - f2′(x) = 9134x - 9134 - f3(x) = (x^2 - 2)^2 #two double roots - f3′(x) = 4x * (x^2 - 2) - f4(x) = sin(x) - x #triple root at 0 - f4′(x) = cos(x) - 1 - f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4 #two quadruple roots - f5′(x) = 8x * (-3 + 2x^2) * (2 - 3x^2 + x^4)^3 - for autodiff in (false, true) - if autodiff - rts1 = newton1d(sin, -5..5) - rts2 = newton1d(f, -∞..∞) - rts3 = newton1d(f1, -10..10) - rts4 = newton1d(f2, -10..11) - rts5 = newton1d(f3, -10..10) - rts6 = newton1d(f4, -10..10) - rts7 = newton1d(f5, -10..10) - - else - rts1 = newton1d(sin, cos, -5..5) - rts2 = newton1d(f, f′, -∞..∞) - rts3 = newton1d(f1, f1′, -10..10) - rts4 = newton1d(f2, f2′, -10..11) - rts5 = newton1d(f3, f3′, -10..10) - rts6 = newton1d(f4, f4′, -10..10) - rts7 = newton1d(f5, f5′, -10..10, reltol=0) - end - - @test length(rts1) == 3 - L = [-pi, 0, pi] - for i = 1:length(rts1) - @test L[i] in rts1[i].interval && :unique == rts1[i].status - end - - @test length(rts2) == 1 - @test (0..0) == rts2[1].interval && :unknown == rts2[1].status - - @test length(rts3) == 4 - L = [1, 2, 3, 4] - for i = 1:length(rts3) - @test L[i] in rts3[i].interval && :unique == rts3[i].status - end - - @test length(rts4) == 1 - @test 1 in rts4[1].interval && :unknown == rts4[1].status - - L1 = [-sqrt(2), sqrt(2)] - for i = 1:length(rts5) - @test L1[i] in rts5[i].interval && :unknown == rts5[i].status - end - - @test length(rts6) == 1 - @test 0 in rts6[1].interval && :unknown == rts6[1].status - - @test length(rts7) == 4 - L = [-sqrt(2), -1, 1, sqrt(2)] - for i = 1:length(rts7) - @test L[i] in rts7[i].interval && :unknown == rts7[i].status - end - - end -end diff --git a/test/quadratic.jl b/test/quadratic.jl deleted file mode 100644 index 5e4dda02..00000000 --- a/test/quadratic.jl +++ /dev/null @@ -1,21 +0,0 @@ -using IntervalArithmetic, IntervalRootFinding - -struct Quad{T} - a::Interval{T} - b::Interval{T} - c::Interval{T} - x::Vector{Interval{T}} -end - - -@testset "Quadratic Interval Equations" begin - rts = Quad{Float64}[] - push!(rts, Quad(interval(1, 5), interval(2, 4), interval(0, 2), [interval(-4, -2), interval(-0, -0)])) - push!(rts, Quad(interval(-1, 1), interval(-2, 2), interval(-1, 1), [interval(-∞, -1), interval(-1, -1), interval(-1, ∞)])) - push!(rts, Quad(interval(1, 2), interval(1, 3), interval(2, 6), [interval(-2, -1)])) - push!(rts, Quad(interval(1, 1), interval(2, 2), interval(1, 1), [interval(-1, -1), interval(-1, -1)])) - - for i in 1:length(rts) - @test quadratic_roots(rts[i].a, rts[i].b, rts[i].c) == rts[i].x - end -end diff --git a/test/roots.jl b/test/roots.jl index a5a2ac8a..c0b9c931 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -1,38 +1,8 @@ -using IntervalArithmetic, IntervalRootFinding, StaticArrays, ForwardDiff -using Test - -function all_unique(rts) - all(root_status.(rts) .== :unique) -end - -function roots_dist(rt1::Root{<:Interval}, rt2::Root{<:Interval}) - d = dist(root_region(rt1), root_region(rt2)) - return sum(d) -end - -function roots_dist(rt1::Root{<:AbstractVector}, rt2::Root{<:AbstractVector}) - return sum(dist.(root_region(rt1), root_region(rt2))) -end - -function roots_dist(rt1::Root{Complex{<:Interval}}, rt2::Root{Complex{<:Interval}}) - dreal = dist(real(root_region(rt1)), real(root_region(rt2))) - dimag = dist(imag(root_region(rt1)), imag(root_region(rt2))) - - return sum(dreal) + sum(dimag) -end - -function test_newtonlike(f, derivative, X, contractor, nsol, tol=1e-10) - rts = roots(f, X ; contractor) - @test length(rts) == nsol - @test all_unique(rts) - @test sum(roots_dist.(rts, roots(f, X ; contractor, derivative))) < tol -end - -@testset "1D roots" begin +@testset "Simple tests" begin # Default rts = roots(sin, interval(-5, 5)) @test length(rts) == 3 - @test all_unique(rts) + @test all(isunique, rts) # Bisection rts = roots(sin, interval(-5, 6) ; contractor = Bisection, abstol = 1e-3) @@ -49,10 +19,6 @@ end test_newtonlike(sin, cos, interval(-5, 5), contractor, 3) - # Infinite interval - rts = roots(x -> x^2 - 2, interval(-Inf, Inf) ; contractor) - @test length(rts) == 2 - # abs g(p) = abs(1 / ( (1+p)^30 ) * 10_000 - 100) g2(p) = (1 / ( (1+p)^30 ) * 10_000 - 100)^2 @@ -75,107 +41,95 @@ end @test length(rts) == 1 @test root_status(only(rts)) == :unknown @test in_interval(0, root_region(only(rts))) - end -end - -function in_solution_set(point, solution_intervals) - return any(map(Y -> in_region(point, Y), solution_intervals)) + # Double reciprocal from #131 + rr(x) = 1/(1/(1 + x)*(1 - x)) + rts = roots(rr, interval(-10, 10) ; contractor) + @test length(rts) == 2 + @test all(root_status(rt) == :unknown for rt in rts) + @test in_interval(-1, root_region(rts[1])) + @test in_interval(1, root_region(rts[2])) + end end -@testset "2D roots" begin - f(x, y) = [x^2 + y^2 - 1, y - 2x] - f(X) = f(X...) - X = [interval(-6, 6), interval(-6, 6)] - - # Bisection - rts = roots(f, X ; contractor = Bisection, abstol = 1e-3) - exact_sol = [sqrt(1/5), 2sqrt(1/5)] - @test in_solution_set(exact_sol, root_region.(rts)) - @test in_solution_set(-exact_sol, root_region.(rts)) +@testset "Multiple roots" begin + f(x) = (x-1) * (x^2 - 2)^3 * (x^3 - 2)^4 + g(x) = exp(x) - x - 1 # Double root at 0 from Burden & Faires, 9th ed, p.84 for contractor in newtonlike_methods - deriv = xx -> ForwardDiff.jacobian(f, xx) - test_newtonlike(f, deriv, X, contractor, 2) - rts = roots(f, X ; contractor) - @test in_solution_set(exact_sol, root_region.(rts)) - @test in_solution_set(-exact_sol, root_region.(rts)) - end - - # Infinite interval - X = [interval(-Inf, Inf), interval(-Inf, Inf)] - rts = roots(f, X ; contractor = Newton) - @test length(rts) == 2 -end + rts = roots(f, interval(-5, 5) ; contractor) + @test length(rts) == 4 + @test rts[1].status == :unknown + @test in_interval(-sqrt(2), rts[1].region) -# From R docs: -# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden + @test rts[2].status == :unique + @test in_interval(1, rts[2].region) -@testset "3D roots" begin - function g(x) - (x1, x2, x3) = x - SVector( x1^2 + x2^2 + x3^2 - 1, - x1^2 + x3^2 - 0.25, - x1^2 + x2^2 - 4x3 - ) - end - dg = xx -> ForwardDiff.jacobian(g, xx) + @test rts[3].status == :unknown + @test in_interval(cbrt(2), rts[3].region) - X = interval(-5, 5) - XX = [X, X, X] + @test rts[4].status == :unknown + @test in_interval(sqrt(2), rts[4].region) - for contractor in newtonlike_methods - rts = roots(g, XX ; contractor) - @test length(rts) == 4 - @test all_unique(rts) - @test all(rts .== roots(g, XX ; contractor, derivative = dg)) + rts = roots(g, interval(-5, 5) ; contractor) + @test all(rt.status == :unknown for rt in rts) + @test in_interval(0, hull([rt.region for rt in rts]...)) end end - -@testset "10D roots" begin - include("../examples/10_dimensional.jl") +@testset "Out of domain" begin + s(x) = sqrt(x^2 - 4) - 2 # root: -6 and 6 + a(xy) = [asin(xy[1]), acos(xy[2]) - π/2] # root: (0, 0) for contractor in newtonlike_methods - rts = roots(f10d, X10d_large ; contractor, max_iteration = 50_000) - @test any(X -> all(in_interval.(sol10d, X.region)), rts) + @suppress @test length(roots(log, interval(-100, 2) ; contractor)) == 1 + @suppress @test length(roots(log, interval(-100, -1) ; contractor)) == 0 - rts = roots(f10d, X10d_close ; contractor, max_iteration = 1_000_000) - @test length(rts) == 1 - @test isunique(rts[1]) + @suppress @test length(roots(s, interval(-20, 20) ; contractor)) == 2 + @suppress @test length(roots(s, interval(0, 20) ; contractor)) == 1 + @suppress @test length(roots(s, interval(-1, 1) ; contractor)) == 0 + + @suppress @test length(roots(a, fill(interval(-20, 20), 2) ; contractor)) == 1 + @suppress @test length(roots(a, fill(interval(0.5, 20), 2) ; contractor)) == 0 end end -@testset "Dimension mismatch" begin - f21(xy) = [xy[1]^2 - 2] - f23(xy) = [xy[1]^2 - 2, xy[2]^2 - 3, xy[1] + xy[2]] - - X = [interval(0, 5), interval(0, 5)] +@testset "Infinite domain" begin + # Infinite interval + f(xy) = [sin(1/xy[1]), xy[1]*xy[2]] # largest roots: [1/π, 0] and [infinity, 0] for contractor in newtonlike_methods - @test_throws DimensionMismatch roots(f21, X ; contractor) - @test_throws DimensionMismatch roots(f23, X ; contractor) - end -end + rts = @suppress roots(x -> x^2 - 2, interval(-Inf, Inf) ; contractor) + @test length(rts) == 2 + @test all(isbounded, [rt.region for rt in rts]) + @test all(rt.region.decoration == dac for rt in rts) + @test all(isunique, rts) -@testset "Out of domain" begin - for contractor in newtonlike_methods - @test length(roots(log, interval(-100, 2) ; contractor)) == 1 - @test length(roots(log, interval(-100, -1) ; contractor)) == 0 - end -end + rts = @suppress roots(f, [interval(0.9/π, Inf), interval(-Inf, Inf)] ; contractor) + @test length(rts) == 2 + sort!(rts, by = mig) + @test isunique(rts[1]) + @test all(isbounded, rts[1].region) + @test minimum(X.decoration for X in rts[1].region) == dac -@testset "Infinite domain" begin - for contractor in newtonlike_methods - rts = roots(x -> x^2 - 2, interval(-Inf, Inf) ; contractor) - @test length(filter(isunique, rts)) == 2 + @test !isunique(rts[2]) + @test !all(isbounded, rts[2].region) + @test minimum(X.decoration for X in rts[2].region) == dac end end @testset "NaN return value" begin - f(xx) = ( (x, y) = xx; [log(y/x) + 3x, y - 2x] ) + f1(x) = x^2/x - 1 + f(xy) = [log(xy[2]/xy[1]) + 3xy[1], xy[2] - 2xy[1]] X = [interval(-100, 100), interval(-100, 100)] + for contractor in newtonlike_methods + rts = roots(f1, interval(-2, 2) ; contractor) + @test length(rts) == 2 + @test !isunique(rts[1]) + @test in_interval(0, root_region(rts[1])) + @test isunique(rts[2]) + rts = roots(f, X ; contractor) @test length(filter(isunique, rts)) == 1 @test length(filter(x -> all(in_interval.(0, x)), root_region.(rts))) == 1 @@ -412,4 +366,59 @@ end g(x) = (x < 1 ? x : error()) rts = roots(g, interval(-10, 10) ; ignored_errors = [ErrorException, IntervalArithmetic.InconclusiveBooleanOperation]) @test any(isunique, rts) +end + +# Wilkinson-type polynomial defined by its roots: +# Wn(x) = (x-1)⋅(x-2)⋅ ⋯ ⋅ (x-n) +W3(x) = prod(x .- (1:3)) +W7(x) = prod(x .- (1:7)) + +# Format: (function, derivative, lower_bound, upper_bound, [true_roots]) +function_list = [ + (sin, cos, -5, 5, [-big(π), 0, big(π)]) , + (cos, x -> -sin(x), -7.5, 7.5, [-3big(π)/2, -big(π)/2, big(π)/2, 3big(π)/2]), + (W3, nothing, -10, 10, [1, 2, 3]), + (W7, nothing, -10, 10, collect(1:7)), + (x -> exp(x) - 2, exp, -20, 20, [log(big(2))]), + (x -> asin(sin(x)) - big"0.1", Returns(1.0), 0, 1, [big"0.1"]), +] + +@testset "Various precisions" begin + @testset "$T" for T in [Float64, BigFloat] + @testset "Function $(func[1])" for func in function_list + @testset "$contractor" for contractor in newtonlike_methods + f, f_prime, a_lower, a_upper, true_roots = func + a = interval(T, a_lower, a_upper) + + @testset "autodiff = $autodiff" for autodiff in (false, true) + isnothing(f_prime) && continue + + if autodiff + rts = roots(f, a ; contractor) + else + rts = roots(f, a ; contractor, derivative = f_prime) + end + + @test length(rts) == length(true_roots) + + for (root, true_root) in zip(rts, true_roots) + @test isa(root, Root) + @test isunique(root) + @test in_interval(true_root, root.region) + end + end + end + end + end +end + + +@testset "Interval parameters" begin + a = interval(1.9, 2.1) + for contractor in newtonlike_methods + rts = roots(x -> x^2 - a^2, interval(0, 10) ; contractor) + + @test length(rts) == 1 + @test issubset_interval(a, rts[1].region) + end end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 7367f5dd..f5e6282d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,16 +1,57 @@ +using BranchAndPrune using IntervalRootFinding using IntervalArithmetic.Symbols -using BranchAndPrune +using ForwardDiff +using StaticArrays +using Suppressor using Test newtonlike_methods = [Newton, Krawczyk] +# Compute the distance between set of roots +function roots_dist(rt1::Root{<:Interval}, rt2::Root{<:Interval}) + d = dist(root_region(rt1), root_region(rt2)) + return sum(d) +end + +function roots_dist(rt1::Root{<:AbstractVector}, rt2::Root{<:AbstractVector}) + return sum(dist.(root_region(rt1), root_region(rt2))) +end + +function roots_dist(rt1::Root{Complex{<:Interval}}, rt2::Root{Complex{<:Interval}}) + dreal = dist(real(root_region(rt1)), real(root_region(rt2))) + dimag = dist(imag(root_region(rt1)), imag(root_region(rt2))) + + return sum(dreal) + sum(dimag) +end + +# Test that a contractor that use derivatives +# Check that all return roots are unique, and that giving the derivative explictly +# makes no difference +function test_newtonlike(f, derivative, X, contractor, nsol, tol=1e-10) + rts = roots(f, X ; contractor) + @test length(rts) == nsol + @test all(isunique, rts) + @test sum(roots_dist.(rts, roots(f, X ; contractor, derivative))) < tol +end + +# Check that a given point is in one of the given intervals +function in_solution_set(point, solution_intervals) + return any(map(Y -> in_region(point, Y), solution_intervals)) +end + include("roots.jl") include("stationary_points.jl") + +include("multidimensional.jl") include("svectors.jl") + include("test_smiley.jl") include("linear_eq.jl") -# include("newton1d.jl") -# include("quadratic.jl") -# include("slopes.jl") +@testset "Example files" begin + for file in readdir("../examples/", join = true) + include(file) + @test true + end +end \ No newline at end of file diff --git a/test/stationary_points.jl b/test/stationary_points.jl index ce742659..e46301c1 100644 --- a/test/stationary_points.jl +++ b/test/stationary_points.jl @@ -1,7 +1,3 @@ -using Test - -import ForwardDiff - @testset "Stationary points" begin f(xx) = ( (x, y) = xx; sin(x) * sin(y) ) gradf = xx -> ForwardDiff.gradient(f, xx) diff --git a/test/test_smiley.jl b/test/test_smiley.jl index 749d08cc..bd1cc456 100644 --- a/test/test_smiley.jl +++ b/test/test_smiley.jl @@ -1,33 +1,23 @@ include("../examples/smiley_examples.jl") -using Test -using IntervalArithmetic, IntervalRootFinding using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55 -function test_all_unique(xs) - for x in xs - @test x.status == :unique - end - return nothing -end - const abstol = 1e-6 -for contractor in (Newton, Krawczyk) # NOTE: Bisection method performs badly in all examples - -@info("Testing method $contractor") @testset "$(SmileyExample22.title)" begin - roots_found = roots(SmileyExample22.f, SmileyExample22.region ; contractor, abstol) - @test length(roots_found) == 8 - test_all_unique(roots_found) - # no reference data for roots given + for contractor in (Newton, Krawczyk) + roots_found = roots(SmileyExample22.f, SmileyExample22.region ; contractor, abstol) + @test length(roots_found) == 8 + @test all(isunique, roots_found) + # no reference data for roots given + end end -for example in (SmileyExample52, SmileyExample54)#, SmileyExample55) - @testset "$(example.title)" begin +@testset "$(example.title)" for example in (SmileyExample52, SmileyExample54) + for contractor in (Newton, Krawczyk) roots_found = roots(example.f, example.region ; contractor, abstol) @test length(roots_found) == length(example.known_roots) - test_all_unique(roots_found) + @test all(isunique, roots_found) for rf in roots_found # check there is exactly one known root for each found root @test 1 == sum(example.known_roots) do rk @@ -36,5 +26,3 @@ for example in (SmileyExample52, SmileyExample54)#, SmileyExample55) end end end - -end # method loop diff --git a/test/wilkinson.jl b/test/wilkinson.jl deleted file mode 100644 index 0ef6b8d2..00000000 --- a/test/wilkinson.jl +++ /dev/null @@ -1,54 +0,0 @@ -# Wilkinson-type polynomial defined by its roots: -# wₙ(x) = (x-1)⋅(x-2)⋅ ⋯ ⋅ (x-n) - -using Polynomials - - -# function generate_wilkinson(n)#, T=BigFloat) # SLOW - -# p = poly(collect(1:n)) -# p_prime = polyder(p) - -# coeffs = p.a #[convert(T, coeff) for coeff in p.a] -# # coeffs_prime = [convert(T, coeff) for coeff in p_prime.a] - -# function f(x) -# total = coeffs[1] -# for i in 2:length(coeffs) -# total += coeffs[i]*x^(i-1) -# end -# total -# end - -# f -# end - -function generate_wilkinson_horner(n) - p = poly(collect(1:n)) # make polynomial with roots 1 to n - - coeffs = p.a - - expr = :($(coeffs[end])) - - for i in length(coeffs)-1:-1:1 - expr = :($expr*x + $(coeffs[i])) - end - - expr -end - -function subscriptify(n::Int) - subscript_digits = [c for c in "₀₁₂₃₄₅₆₇₈₉"] - dig = reverse(digits(n)) - join([subscript_digits[i+1] for i in dig]) -end - -# Generate Wilkinson functions W₃ etc.: -for n in 1:15 - - fn_name = Symbol("W", subscriptify(n)) - - expr = generate_wilkinson_horner(n) - - @eval $(fn_name)(x) = $(expr) -end