diff --git a/.travis.yml b/.travis.yml index d826e155..24a3db06 100644 --- a/.travis.yml +++ b/.travis.yml @@ -7,7 +7,8 @@ os: - osx julia: - - 0.6 + - 0.7 + - 1.0 - nightly matrix: diff --git a/REQUIRE b/REQUIRE index 4e1f6eb9..1cce34e3 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,5 +1,5 @@ -julia 0.6 0.7 -IntervalArithmetic 0.14 0.15 -ForwardDiff 0.7.5 -StaticArrays 0.7.1 +julia 0.7 +IntervalArithmetic 0.15 +ForwardDiff 0.8 +StaticArrays 0.8 Polynomials diff --git a/appveyor.yml b/appveyor.yml index f7bfae29..c2588f1d 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,13 +1,18 @@ environment: matrix: - - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe" - - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe" + - julia_version: 0.7 + - julia_version: 1 + - julia_version: nightly -matrix: - allow_failures: - - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe" - - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe" +platform: + - x86 # 32-bit + - x64 # 64-bit +# # Uncomment the following lines to allow failures on nightly julia +# # (tests will run but not make your overall status red) +# matrix: +# allow_failures: +# - julia_version: nightly branches: only: @@ -21,19 +26,18 @@ notifications: on_build_status_changed: false install: - - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12" -# Download most recent Julia Windows binary - - ps: (new-object net.webclient).DownloadFile( - $env:JULIA_URL, - "C:\projects\julia-binary.exe") -# Run installer silently, output to C:\projects\julia - - C:\projects\julia-binary.exe /S /D=C:\projects\julia + - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1")) build_script: -# Need to convert from shallow to complete for Pkg.clone to work - - IF EXIST .git\shallow (git fetch --unshallow) - - C:\projects\julia\bin\julia -e "versioninfo(); - Pkg.clone(pwd(), \"IntervalRootFinding\"); Pkg.build(\"IntervalRootFinding\")" + - echo "%JL_BUILD_SCRIPT%" + - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%" test_script: - - C:\projects\julia\bin\julia -e "Pkg.test(\"IntervalRootFinding\")" + - echo "%JL_TEST_SCRIPT%" + - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%" + +# # Uncomment to support code coverage upload. Should only be enabled for packages +# # which would have coverage gaps without running on Windows +# on_success: +# - echo "%JL_CODECOV_SCRIPT%" +# - C:\julia\bin\julia -e "%JL_CODECOV_SCRIPT%" diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index d20e85c6..d70ed783 100644 --- a/src/IntervalRootFinding.jl +++ b/src/IntervalRootFinding.jl @@ -8,6 +8,8 @@ using IntervalArithmetic using ForwardDiff using StaticArrays +using LinearAlgebra: I, Diagonal + import Base: ⊆, show, big, \ @@ -120,23 +122,6 @@ function find_roots_midpoint(f::Function, a::Real, b::Real, method::Function=new end -function Base.lexcmp(a::Interval{T}, b::Interval{T}) where {T} - #@show a, b - if a.lo < b.lo - return -1 - elseif a.lo == b.lo - if a.hi < b.hi - return -1 - else - return 0 - end - end - return 1 - -end - -Base.lexcmp(a::Root{T}, b::Root{T}) where {T} = lexcmp(a.interval, b.interval) - function clean_roots(f, roots) diff --git a/src/contractors.jl b/src/contractors.jl index ca638e84..7a270d9f 100644 --- a/src/contractors.jl +++ b/src/contractors.jl @@ -1,14 +1,14 @@ export Bisection, Newton, Krawczyk -doc""" +""" Contractor{F} Abstract type for contractors. """ abstract type Contractor{F} end -doc""" +""" Bisection{F} <: Contractor{F} Contractor type for the bisection method. @@ -27,7 +27,7 @@ function (contractor::Bisection)(X, tol) return :unknown, X end -doc""" +""" newtonlike_contract(op, X, tol) Contraction operation for contractors using the first derivative of the @@ -57,7 +57,7 @@ function newtonlike_contract(op, C, X, tol) return :unknown, NX end -doc""" +""" Newton{F, FP} <: Contractor{F} Contractor type for the Newton method. @@ -76,7 +76,7 @@ function (C::Newton)(X, tol) end -doc""" +""" Single-variable Newton operator """ function 𝒩(f, X::Interval{T}) where {T} @@ -97,7 +97,7 @@ function 𝒩(f, X::Interval{T}, dX::Interval{T}) where {T} m - (f(m) / dX) end -doc""" +""" Multi-variable Newton operator. """ function 𝒩(f::Function, jacobian::Function, X::IntervalBox) # multidimensional Newton operator @@ -108,7 +108,7 @@ function 𝒩(f::Function, jacobian::Function, X::IntervalBox) # multidimension end -doc""" +""" Krawczyk{F, FP} <: Contractor{F} Contractor type for the Krawczyk method. @@ -127,7 +127,7 @@ function (C::Krawczyk)(X, tol) end -doc""" +""" Single-variable Krawczyk operator """ function 𝒦(f, f′, X::Interval{T}) where {T} @@ -137,7 +137,7 @@ function 𝒦(f, f′, X::Interval{T}) where {T} m - Y*f(m) + (1 - Y*f′(X)) * (X - m) end -doc""" +""" Multi-variable Krawczyk operator """ function 𝒦(f, jacobian, X::IntervalBox{T}) where {T} diff --git a/src/krawczyk.jl b/src/krawczyk.jl index 3bd6e64f..206fc34f 100644 --- a/src/krawczyk.jl +++ b/src/krawczyk.jl @@ -4,7 +4,7 @@ # const D = differentiate -doc"""Returns two intervals, the first being a point within the +"""Returns two intervals, the first being a point within the interval x such that the interval corresponding to the derivative of f there does not contain zero, and the second is the inverse of its derivative""" function guarded_derivative_midpoint(f::Function, f_prime::Function, x::Interval{T}) where {T} diff --git a/src/linear_eq.jl b/src/linear_eq.jl index b4242b17..1dc43b8a 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -86,7 +86,6 @@ end function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true) - n = size(A, 1) x = similar(b) x .= -1e16..1e16 x = gauss_elimination_interval!(x, A, b, precondition=precondition) @@ -94,47 +93,54 @@ function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precond return x end """ -Solves the system of linear equations using Gaussian Elimination, -with (or without) preconditioning. (kwarg - `precondition`) -Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis - Page 72 +Solves the system of linear equations using Gaussian Elimination. +Preconditioning is used when the `precondition` keyword argument is `true`. + +REF: Luc Jaulin et al., +*Applied Interval Analysis*, pg. 72 """ -function gauss_elimination_interval!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true) +function gauss_elimination_interval!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true) if precondition - ((a, b) = preconditioner(a, b)) + (A, b) = preconditioner(A, b) else - a = copy(a) + A = copy(A) b = copy(b) end - n = size(a, 1) - p = zeros(x) + n = size(A, 1) + + p = similar(b) + p .= 0 for i in 1:(n-1) - if 0 ∈ a[i, i] # diagonal matrix is not invertible + if 0 ∈ A[i, i] # diagonal matrix is not invertible p .= entireinterval(b[1]) - return p .∩ x + return p .∩ x # return x? end for j in (i+1):n - α = a[j, i] / a[i, i] + α = A[j, i] / A[i, i] b[j] -= α * b[i] for k in (i+1):n - a[j, k] -= α * a[i, k] + A[j, k] -= α * A[i, k] end end end for i in n:-1:1 - sum = 0 + + temp = zero(b[1]) + for j in (i+1):n - sum += a[i, j] * p[j] + temp += A[i, j] * p[j] end - p[i] = (b[i] - sum) / a[i, i] + + p[i] = (b[i] - temp) / A[i, i] end - p .∩ x + return p .∩ x end function gauss_elimination_interval1(A::AbstractMatrix, b::AbstractArray; precondition=true) diff --git a/src/newton.jl b/src/newton.jl index 3c29919c..be08fbca 100644 --- a/src/newton.jl +++ b/src/newton.jl @@ -2,9 +2,7 @@ # Newton method - - -doc"""If a root is known to be inside an interval, +"""If a root is known to be inside an interval, `newton_refine` iterates the interval Newton method until that root is found.""" function newton_refine(f::Function, f_prime::Function, X::Union{Interval{T}, IntervalBox{N,T}}; tolerance=eps(T), debug=false) where {N,T} @@ -30,7 +28,7 @@ end -doc"""If a root is known to be inside an interval, +"""If a root is known to be inside an interval, `newton_refine` iterates the interval Newton method until that root is found.""" function newton_refine(f::Function, f_prime::Function, X::Interval{T}; tolerance=eps(T), debug=false) where {T} @@ -60,11 +58,10 @@ end -doc"""`newton` performs the interval Newton method on the given function `f` +"""`newton` performs the interval Newton method on the given function `f` with its optional derivative `f_prime` and initial interval `x`. Optional keyword arguments give the `tolerance`, `maxlevel` at which to stop subdividing, and a `debug` boolean argument that prints out diagnostic information.""" - function newton(f::Function, f_prime::Function, x::Interval{T}, level::Int=0; tolerance=eps(T), debug=false, maxlevel=30) where {T} diff --git a/src/newton1d.jl b/src/newton1d.jl index 98d22cee..e302be2b 100644 --- a/src/newton1d.jl +++ b/src/newton1d.jl @@ -3,9 +3,8 @@ with its derivative `f′` and initial interval `x`. Optional keyword arguments give the tolerances `reltol` and `abstol`. `reltol` is the tolerance on the relative error whereas `abstol` is the tolerance on |f(X)|, and a `debug` boolean argument that prints out diagnostic information.""" - -function newton1d{T}(f::Function, f′::Function, x::Interval{T}; - reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) +function newton1d(f::Function, f′::Function, x::Interval{T}; + reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) where {T} L = Interval{T}[] # Array to hold the intervals still to be processed @@ -201,6 +200,5 @@ end Optional keyword arguments give the tolerances `reltol` and `abstol`. `reltol` is the tolerance on the relative error whereas `abstol` is the tolerance on |f(X)|, and a `debug` boolean argument that prints out diagnostic information.""" - -newton1d{T}(f::Function, x::Interval{T}; args...) = +newton1d(f::Function, x::Interval{T}; args...) where {T} = newton1d(f, x->D(f,x), x; args...) diff --git a/src/quadratic.jl b/src/quadratic.jl index 4d1f517d..c7fe6f6f 100644 --- a/src/quadratic.jl +++ b/src/quadratic.jl @@ -2,7 +2,7 @@ Helper function for quadratic_interval that computes roots of a real quadratic using interval arithmetic to bound rounding errors. """ -function quadratic_helper!{T}(a::Interval{T}, b::Interval{T}, c::Interval{T}, L::Array{Interval{T}}) +function quadratic_helper!(a::Interval{T}, b::Interval{T}, c::Interval{T}, L::Array{Interval{T}}) where {T} Δ = b^2 - 4*a*c @@ -40,7 +40,7 @@ This algorithm finds the set of points where `F.lo(x) ≥ 0` and the set of points where `F.hi(x) ≤ 0` and takes the intersection of these two sets. Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 8 """ -function quadratic_roots{T}(a::Interval{T}, b::Interval{T}, c::Interval{T}) +function quadratic_roots(a::Interval{T}, b::Interval{T}, c::Interval{T}) where {T} L = Interval{T}[] R = Interval{T}[] diff --git a/src/root_object.jl b/src/root_object.jl index cd684401..fec0374c 100644 --- a/src/root_object.jl +++ b/src/root_object.jl @@ -16,7 +16,7 @@ root_status(x::Root) = x.status show(io::IO, root::Root) = print(io, "Root($(root.interval), :$(root.status))") -isunique{T}(root::Root{T}) = (root.status == :unique) +isunique(root::Root{T}) where {T} = (root.status == :unique) ⊆(a::Interval, b::Root) = a ⊆ b.interval # the Root object has the interval in the first entry ⊆(a::Root, b::Root) = a.interval ⊆ b.interval diff --git a/src/roots.jl b/src/roots.jl index 9a8830c7..4a62e28a 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -1,9 +1,9 @@ import IntervalArithmetic: diam, isinterior -import Base: start, next, done, copy, eltype, iteratorsize +import Base: iterate, eltype, IteratorSize, copy export branch_and_prune, Bisection, Newton, RootSearch -export start, next, done, copy, step!, eltype, iteratorsize +export step! diam(x::Root) = diam(x.interval) @@ -11,7 +11,7 @@ struct RootSearchState{T <: Union{Interval,IntervalBox}} working::Vector{T} outputs::Vector{Root{T}} end -RootSearchState{T<:Union{Interval,IntervalBox}}(region::T) = +RootSearchState(region::T) where {T<:Union{Interval,IntervalBox}} = RootSearchState([region], Root{T}[]) copy(state::RootSearchState) = @@ -31,16 +31,24 @@ struct RootSearch{R <: Union{Interval,IntervalBox}, S <: Contractor, T <: Real} tolerance::T end -eltype{R, T <: RootSearch{R}}(::Type{T}) = RootSearchState{R} -iteratorsize{T <: RootSearch}(::Type{T}) = Base.SizeUnknown() +eltype(::Type{T}) where {R, T <: RootSearch{R}} = RootSearchState{R} +IteratorSize(::Type{T}) where {T <: RootSearch} = Base.SizeUnknown() -function start(iter::RootSearch) +# function start(iter::RootSearch) +# state = RootSearchState(iter.region) +# sizehint!(state.outputs, 100) +# sizehint!(state.working, 1000) +# return state +# end + +function iterate(iter::RootSearch) state = RootSearchState(iter.region) sizehint!(state.outputs, 100) sizehint!(state.working, 1000) - return state + return state, state end + """ step!(state::RootSearchState, contractor, tolerance) @@ -63,12 +71,21 @@ function step!(state::RootSearchState, contractor, tolerance) return nothing end -function next(iter::RootSearch, state::RootSearchState) +# function next(iter::RootSearch, state::RootSearchState) +# step!(state, iter.contractor, iter.tolerance) +# return state, state +# end + +function iterate(iter::RootSearch, state::RootSearchState) + + isempty(state.working) && return nothing + step!(state, iter.contractor, iter.tolerance) return state, state end -done(iter::RootSearch, state::RootSearchState) = isempty(state.working) + +# done(iter::RootSearch, state::RootSearchState) = isempty(state.working) """ branch_and_prune(X, contract, tol=1e-3) @@ -85,10 +102,12 @@ Inputs: """ function branch_and_prune(X, contractor, tol=1e-3) iter = RootSearch(X, contractor, tol) - local state + local output # complete iteration - for state in iter end - return state.outputs + for state in iter + output = state.outputs + end + return output end export recursively_branch_and_prune @@ -109,6 +128,7 @@ end const IntervalLike{T} = Union{Interval{T}, IntervalBox{T}} const NewtonLike = Union{Type{Newton}, Type{Krawczyk}} + """ roots(f, X, contractor, tol=1e-3) roots(f, deriv, X, contractor, tol=1e-3) @@ -123,9 +143,7 @@ Inputs: the status of a given box `X`. It returns the new box and a symbol indicating the status. Current possible values are `Bisection`, `Newton` and `Krawczyk` - `deriv` ; explicit derivative of `f` for `Newton` and `Krawczyk` - """ -# Contractor specific `roots` functions function roots(f, X::IntervalLike{T}, ::Type{Bisection}, tol::Float64=1e-3) where {T} branch_and_prune(X, Bisection(f), tol) end diff --git a/test/bisection.jl b/test/bisection.jl index e602e89e..a59a4b4a 100644 --- a/test/bisection.jl +++ b/test/bisection.jl @@ -1,5 +1,5 @@ using ValidatedNumerics, ValidatedNumerics.RootFinding -using Base.Test +using Test @testset "Bisection tests" begin @testset "Single variable" begin diff --git a/test/findroots.jl b/test/findroots.jl index 703750d6..602ff287 100644 --- a/test/findroots.jl +++ b/test/findroots.jl @@ -1,7 +1,7 @@ using IntervalArithmetic, IntervalRootFinding using ForwardDiff -using Base.Test +using Test const D = IntervalRootFinding.derivative diff --git a/test/linear_eq.jl b/test/linear_eq.jl index 2cc49d08..8b720998 100644 --- a/test/linear_eq.jl +++ b/test/linear_eq.jl @@ -1,6 +1,7 @@ using IntervalArithmetic, StaticArrays, IntervalRootFinding +using Test -function randVec(n::Int) +function rand_vec(n::Int) a = randn(n) A = Interval.(a) mA = MVector{n}(A) @@ -8,7 +9,7 @@ function randVec(n::Int) return A, mA, sA end -function randMat(n::Int) +function rand_mat(n::Int) a = randn(n, n) A = Interval.(a) mA = MMatrix{n, n}(A) @@ -18,23 +19,32 @@ end @testset "Linear Equations" begin - A = [[2..3 0..1;1..2 2..3], ] - b = [[0..120, 60..240], ] - x = [[-120..90, -60..240], ] + As = [[2..3 0..1; 1..2 2..3], ] + bs = [[0..120, 60..240], ] + xs = [[-120..90, -60..240], ] for i in 1:10 - rand_a = randMat(i)[1] - rand_b = randVec(i)[1] - rand_c = rand_a * rand_b - push!(A, rand_a) - push!(b, rand_c) - push!(x, rand_b) + rand_A = rand_mat(i)[1] + rand_x = rand_vec(i)[1] + rand_b = rand_A * rand_x + push!(As, rand_A) + push!(bs, rand_b) + push!(xs, rand_x) end - for i in 1:length(A) - for precondition in (false, true) - for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval, \) - @test all(x[i] .⊆ f(A[i], b[i])) - end + + n = length(As) + + + for solver in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval, \) + @testset "Solver $solver" begin + #for precondition in (false, true) + + for i in 1:n + soln = solver(As[i], bs[i]) + @test all(xs[i] .⊆ soln) + + end + #end end end end diff --git a/test/newton1d.jl b/test/newton1d.jl index 099dce00..d83229f6 100644 --- a/test/newton1d.jl +++ b/test/newton1d.jl @@ -1,7 +1,7 @@ using IntervalArithmetic, IntervalRootFinding using ForwardDiff -using Base.Test +using Test const D = IntervalRootFinding.derivative diff --git a/test/roots.jl b/test/roots.jl index 15459ee0..87ef6d5d 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -1,6 +1,6 @@ -using IntervalArithmetic, IntervalRootFinding, StaticArrays -using Base.Test +using IntervalArithmetic, IntervalRootFinding, StaticArrays, ForwardDiff +using Test function all_unique(rts) all(root_status.(rts) .== :unique) @@ -118,7 +118,7 @@ end @testset "RootSearch interface" begin contractor = Newton(sin, cos) search = RootSearch(-10..10, contractor, 1e-3) - state = start(search) + state, _ = iterate(search) # check that original and copy are independent state_copy = copy(state) @@ -127,5 +127,5 @@ end # cover optional iterator methods @test eltype(search) != Any - @test_nowarn iteratorsize(search) + # @test_nowarn iteratorsize(search) end diff --git a/test/runtests.jl b/test/runtests.jl index 0de6225b..6260591e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ using IntervalRootFinding -using Base.Test +using Test #include("bisect.jl") #include("findroots.jl") diff --git a/test/slopes.jl b/test/slopes.jl index 8513899e..5c5d51a6 100644 --- a/test/slopes.jl +++ b/test/slopes.jl @@ -1,6 +1,6 @@ using IntervalArithmetic, IntervalRootFinding using ForwardDiff -using Base.Test +using Test struct Slopes{T} f::Function diff --git a/test/smiley_examples.jl b/test/smiley_examples.jl index 92f9125f..2a26d38e 100644 --- a/test/smiley_examples.jl +++ b/test/smiley_examples.jl @@ -54,7 +54,7 @@ f(x) = SVector( (x[1]^2 + x[2]^2 + x[3]^2 - r1^2) * (x[1]^2 + x[2]^2 + x[3]^2 - r2^2), (a0 * x[1] + b0 * x[2] + c0 * x[3] - d01) * (a0 * x[1] + b0 * x[2] + c0 * x[3] - d02), - prod(as * x[1] .+ bs * x[2] .+ cs * x[3] .- ds) + prod(as[i] * x[1] + bs[i] * x[2] + cs[i] * x[3] - ds[i] for i in 1:m) ) # contains all 24 reported roots diff --git a/test/test_smiley.jl b/test/test_smiley.jl index 88c92650..70cd9299 100644 --- a/test/test_smiley.jl +++ b/test/test_smiley.jl @@ -1,8 +1,8 @@ include("smiley_examples.jl") -using Base.Test +using Test using IntervalArithmetic, IntervalRootFinding -using SmileyExample22, SmileyExample52, SmileyExample54, SmileyExample55 +using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55 function test_all_unique(xs) for x in xs @@ -14,7 +14,7 @@ end const tol = 1e-6 const method = Newton # NOTE: Bisection method performs badly in all examples -info("testing method $(method)") +@info("Testing method $(method)") @testset "$(SmileyExample22.title)" begin roots_found = roots(SmileyExample22.f, SmileyExample22.region, method, tol)