Partial update to IntervalArithmetic v0.21

.gitignore

@@ -1,4 +1,5 @@
 docs/build/
 docs/site/
 benchmark/tune.json
-Manifest.toml
+Manifest.toml
+.vscode

Project.toml

@@ -14,7 +14,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 [compat]
 BranchAndPrune = "0.2.1"
 ForwardDiff = "0.10"
-IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
+IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21"
 Polynomials = "0.5, 0.6, 0.7, 0.8, 1, 2, 3"
 Reexport = "1"
 StaticArrays = "0.11, 0.12, 1.0"

benchmark/benchmarks.jl

@@ -49,8 +49,8 @@ sizes = (2, 5, 10)
 
 for n in sizes
     s = S["n = $n"] = BenchmarkGroup()
-    A = Interval.(randn(n, n))
-    b = Interval.(randn(n))
+    A = interval.(randn(n, n))
+    b = interval.(randn(n))
 
     s["Gauss seidel"] = @benchmarkable gauss_seidel_interval($A, $b)
     s["Gauss seidel contractor"] = @benchmarkable gauss_seidel_contractor($A, $b)
@@ -59,7 +59,7 @@ end
 
 
 S = SUITE["Dietmar-Ratz"] = BenchmarkGroup()
-X = Interval(0.75, 1.75)
+X = interval(0.75, 1.75)
 
 for (k, dr) in enumerate(dr_functions)
     s = S["Dietmar-Ratz $k"] = BenchmarkGroup()

perf/linear_eq.jl

@@ -2,15 +2,15 @@ using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmeti
 
 function randVec(n::Int)
     a = randn(n)
-    A = Interval.(a)
+    A = interval.(a)
     mA = MVector{n}(A)
     sA = SVector{n}(A)
     return A, mA, sA
 end
 
 function randMat(n::Int)
     a = randn(n, n)
-    A = Interval.(a)
+    A = interval.(a)
     mA = MMatrix{n, n}(A)
     sA = SMatrix{n, n}(A)
     return A, mA, sA

src/IntervalRootFinding.jl

@@ -31,7 +31,10 @@ export isunique, root_status
 using Reexport
 @reexport using IntervalArithmetic
 
-import IntervalArithmetic: interval, wideinterval
+import IntervalArithmetic: interval
+
+wideinterval(x) = interval(prevfloat(x), nextfloat(x))
+IntervalArithmetic.mid(a, α) = scaled_mid(a, α)
 
 
 
@@ -83,13 +86,13 @@ function find_roots(f::Function, a::Real, b::Real, method::Function=newton;
 
     if precision >= 0
         setprecision(Interval, precision) do
-            find_roots(f, @interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
+            find_roots(f, interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
         end
 
 
     else  # use current precision
 
-        find_roots(f, @interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
+        find_roots(f, interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
 
     end
 end

src/complex.jl

@@ -7,7 +7,9 @@ struct Compl{T} <: Number
 end
 
 Base.promote_rule(::Type{Compl{T}}, ::Type{S}) where {T, S<:Real} = Compl{T}
-Base.convert(::Type{Compl{T}}, x::Real) where {T} = Compl{T}(x, 0)
+Base.convert(::Type{Compl{T}}, x::Real) where {T} = Compl{T}(x, zero(T))
+Base.convert(::Type{Compl{T}}, x::Real) where {S,T<:Interval{S}} = Compl{T}(interval(S, x), zero(T))
+
 
 Base.show(io::IO, c::Compl) = println(io, c.re, " + ", c.im, "im")
 

src/contractors.jl

@@ -42,12 +42,12 @@ struct Newton{F, FP} <: AbstractContractor{F}
 end
 
 function (N::Newton)(X::Interval ; α=where_bisect)
-    m = Interval(mid(X, α))
+    m = interval(mid(X, α))
     return m - (N.f(m) / N.f′(X))
 end
 
 function (N::Newton)(X::IntervalBox ; α=where_bisect)
-    m = Interval.(mid(X, α))
+    m = interval.(mid(X, α))
     J = N.f′(X)
     y = gauss_elimination_interval(J, N.f(m))  # J \ f(m)
     return IntervalBox(m .- y)
@@ -79,7 +79,7 @@ struct Krawczyk{F, FP} <: AbstractContractor{F}
 end
 
 function (K::Krawczyk)(X::Interval ; α=where_bisect)
-    m = Interval(mid(X, α))
+    m = interval(mid(X, α))
     Y = 1 / K.f′(m)
 
     return m - Y*K.f(m) + (1 - Y*K.f′(X)) * (X - m)
@@ -135,7 +135,7 @@ function contract(C::Union{Newton, Krawczyk}, R::Root)
 
     NX = contracted_X ∩ X
 
-    isinf(X) && return Root(NX, :unknown)  # force bisection
+    isentire(X) && return Root(NX, :unknown)  # force bisection
     safe_isempty(NX) && return Root(X, :empty)
 
     if R.status == :unique || NX ⪽ X  # isinterior, we know there's a unique root inside

src/krawczyk.jl

@@ -11,14 +11,14 @@ function guarded_derivative_midpoint(f::Function, f_prime::Function, x::Interval
 
     α = convert(T, 0.46875)   # close to 0.5, but exactly representable as a floating point
 
-    m = Interval(mid(x))
+    m = interval(mid(x))
 
     C = inv(f_prime(m))
 
     # Check that 0 is not in C; if so, consider another point rather than m
     i = 0
     while zero(T) ∈ C || isempty(C)
-        m = Interval( α*x.lo + (one(T)-α)*x.hi )
+        m = interval( α*x.lo + (one(T)-α)*x.hi )
         C = inv(f_prime(m))
 
         i += 1
@@ -84,9 +84,9 @@ function krawczyk(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
 
     # bisecting
     roots = vcat(
-        krawczyk(f, f_prime, Interval(x.lo, m), level+1,
+        krawczyk(f, f_prime, interval(x.lo, m), level+1,
                  tolerance=tolerance, debug=debug, maxlevel=maxlevel),
-        krawczyk(f, f_prime, Interval(m, x.hi), level+1,
+        krawczyk(f, f_prime, interval(m, x.hi), level+1,
                  tolerance=tolerance, debug=debug, maxlevel=maxlevel)
         )
 
@@ -103,7 +103,7 @@ krawczyk(f::Function,x::Interval{T}; args...) where {T} =
 
 # krawczyk for vector of intervals:
 krawczyk(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) where {T} =
-    vcat([krawczyk(f, f_prime, @interval(x); args...) for x in xx]...)
+    vcat([krawczyk(f, f_prime, x; args...) for x in xx]...)
 
 krawczyk(f::Function,  xx::Vector{Interval{T}}, level; args...) where {T} =
     krawczyk(f, x->D(f,x), xx; args...)

src/linear_eq.jl

@@ -67,9 +67,9 @@ function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::Abstra
     extdiagA = copy(A)
     for i in 1:n
         if (typeof(b) <: SArray)
-            extdiagA = setindex(extdiagA, Interval(0), i, i)
+            extdiagA = setindex(extdiagA, interval(0), i, i)
         else
-            extdiagA[i, i] = Interval(0)
+            extdiagA[i, i] = interval(0)
         end
     end
     inv_diagA = inv(diagA)
@@ -113,8 +113,7 @@ function gauss_elimination_interval!(x::AbstractArray, A::AbstractMatrix, b::Abs
 
     n = size(A, 1)
 
-    p = similar(b)
-    p .= 0
+    p = zeros(eltype(b), length(b))
 
     for i in 1:(n-1)
         if 0 ∈ A[i, i] # diagonal matrix is not invertible

src/newton.jl

@@ -65,9 +65,9 @@ function newton(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
 
         # bisect:
         roots = vcat(
-            newton(f, f_prime, Interval(x.lo, m), level+1,
+            newton(f, f_prime, interval(x.lo, m), level+1,
                    tolerance=tolerance, debug=debug, maxlevel=maxlevel),
-            newton(f, f_prime, Interval(m, x.hi), level+1,
+            newton(f, f_prime, interval(m, x.hi), level+1,
                    tolerance=tolerance, debug=debug, maxlevel=maxlevel)
             # note the nextfloat here to prevent repetition
             )
@@ -76,8 +76,8 @@ function newton(f::Function, f_prime::Function, x::Interval{T}, level::Int=0;
 
     else  # 0 in deriv; this does extended interval division by hand
 
-        y1 = Interval(deriv.lo, -z)
-        y2 = Interval(z, deriv.hi)
+        y1 = interval(deriv.lo, -z)
+        y2 = interval(z, deriv.hi)
 
         if debug
             @show (y1, y2)
@@ -114,7 +114,7 @@ newton(f::Function, x::Interval{T};  args...) where {T} =
 
 # newton for vector of intervals:
 newton(f::Function, f_prime::Function, xx::Vector{Interval{T}}; args...) where {T} =
-    vcat([newton(f, f_prime, @interval(x); args...) for x in xx]...)
+    vcat([newton(f, f_prime, x; args...) for x in xx]...)
 
 newton(f::Function,  xx::Vector{Interval{T}}, level; args...) where {T} =
     newton(f, x->D(f,x), xx, 0, args...)

src/quadratic.jl

@@ -45,27 +45,27 @@ function quadratic_roots(a::Interval{T}, b::Interval{T}, c::Interval{T}) where {
     L = Interval{T}[]
     R = Interval{T}[]
 
-    quadratic_helper!(Interval(a.lo), Interval(b.lo), Interval(c.lo), L)
-    quadratic_helper!(Interval(a.hi), Interval(b.hi), Interval(c.hi), L)
-    quadratic_helper!(Interval(a.lo), Interval(b.hi), Interval(c.lo), L)
-    quadratic_helper!(Interval(a.hi), Interval(b.lo), Interval(c.hi), L)
+    quadratic_helper!(interval(a.lo), interval(b.lo), interval(c.lo), L)
+    quadratic_helper!(interval(a.hi), interval(b.hi), interval(c.hi), L)
+    quadratic_helper!(interval(a.lo), interval(b.hi), interval(c.lo), L)
+    quadratic_helper!(interval(a.hi), interval(b.lo), interval(c.hi), L)
 
     if (length(L) == 8)
         resize!(L, 4)
     end
 
     if (a.lo < 0 || (a.lo == 0 && b.hi == 0) || (a.lo == 0 && b.hi == 0 && c.lo ≤ 0))
-        push!(L, Interval(-∞))
+        push!(L, -∞..nextfloat(-∞))
     end
 
     if (a.lo < 0 || (a.lo == 0 && b.lo == 0) || (a.lo == 0 && b.lo == 0 && c.lo ≤ 0))
-        push!(L, Interval(∞))
+        push!(L, prevfloat(∞)..∞)
     end
 
     sort!(L, by = x -> x.lo)
 
     for i in 1:2:length(L)
-        push!(R, Interval(L[i].lo, L[i+1].hi))
+        push!(R, interval(L[i].lo, L[i+1].hi))
     end
 
     R

src/roots.jl

@@ -10,7 +10,7 @@ isnan(r::Root) = isnan(interval(r))
 struct RootProblem{T}
     abstol::T
 end
-   
+
 function bisect(r::Root)
     Y1, Y2 = bisect(interval(r))
     return Root(Y1, :unknown), Root(Y2, :unknown)

src/slopes.jl

@@ -14,7 +14,7 @@ struct Slope{T}
 end
 
 Slope(c) = Slope(c, c, 0)
-Slope(a, b, c) = Slope(promote(convert(Interval, a), b, c)...)
+Slope(a, b, c) = Slope(promote(interval(a), b, c)...)
 
 function slope_var(v::Number)
     Slope(v, v, 1)

test/bisection.jl

@@ -8,9 +8,9 @@ using Test
         roots = bisection(f, -∞..∞, tolerance=1e-3)
         roots2 = [root.interval for root in roots]
 
-        @test roots2 ==  [  Interval(-0.8408203124999999, -0.8398437499999999),
-                            Interval(3.6425781249999996, 3.6435546874999996),
-                            Interval(6.062499999999999, 6.063476562499999)
+        @test roots2 ==  [  interval(-0.8408203124999999, -0.8398437499999999),
+                            interval(3.6425781249999996, 3.6435546874999996),
+                            interval(6.062499999999999, 6.063476562499999)
                          ]
 
     end
@@ -21,13 +21,13 @@ using Test
         X = (-∞..∞) × (-∞..∞)
         roots = bisection(g, X, tolerance=1e-3)
         roots2 = [root.interval for root in roots]
-        
-        @test roots2 == [IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7080078124999999, -0.7070312499999999)),
- IntervalBox(Interval(-0.7080078124999999, -0.7070312499999999), Interval(-0.7070312499999999, -0.7060546874999999)),
- IntervalBox(Interval(-0.7070312499999999, -0.7060546874999999), Interval(-0.7080078124999999, -0.7070312499999999)),
- IntervalBox(Interval(0.7060546874999999, 0.7070312499999999), Interval(0.7070312499999999, 0.7080078124999999)),
- IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7060546874999999, 0.7070312499999999)),
- IntervalBox(Interval(0.7070312499999999, 0.7080078124999999), Interval(0.7070312499999999, 0.7080078124999999))]
+
+        @test roots2 == [IntervalBox(interval(-0.7080078124999999, -0.7070312499999999), interval(-0.7080078124999999, -0.7070312499999999)),
+ IntervalBox(interval(-0.7080078124999999, -0.7070312499999999), interval(-0.7070312499999999, -0.7060546874999999)),
+ IntervalBox(interval(-0.7070312499999999, -0.7060546874999999), interval(-0.7080078124999999, -0.7070312499999999)),
+ IntervalBox(interval(0.7060546874999999, 0.7070312499999999), interval(0.7070312499999999, 0.7080078124999999)),
+ IntervalBox(interval(0.7070312499999999, 0.7080078124999999), interval(0.7060546874999999, 0.7070312499999999)),
+ IntervalBox(interval(0.7070312499999999, 0.7080078124999999), interval(0.7070312499999999, 0.7080078124999999))]
 
     end
 

test/findroots.jl

@@ -141,8 +141,8 @@ end
 
         @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)
+        @test rts[1].interval == interval(-1.4142135623730954, -1.414213562373095)
+        @test rts[3].interval == interval(1.259921049894873, 1.2599210498948734)
 
     end
 end

test/linear_eq.jl

@@ -3,15 +3,15 @@ using Test
 
 function rand_vec(n::Int)
     a = randn(n)
-    A = Interval.(a)
+    A = interval.(a)
     mA = MVector{n}(A)
     sA = SVector{n}(A)
     return A, mA, sA
 end
 
 function rand_mat(n::Int)
     a = randn(n, n)
-    A = Interval.(a)
+    A = interval.(a)
     mA = MMatrix{n, n}(A)
     sA = SMatrix{n, n}(A)
     return A, mA, sA

test/newton1d.jl

@@ -1,4 +1,3 @@
-
 using IntervalArithmetic, IntervalRootFinding
 using ForwardDiff
 using Test

test/roots.jl

@@ -1,4 +1,3 @@
-
 using IntervalArithmetic, IntervalRootFinding, StaticArrays, ForwardDiff
 using Test
 

test/slopes.jl

@@ -11,17 +11,17 @@ end
 
 @testset "Automatic slope expansion" begin
     for T in [Float64, BigFloat]
-        s = interval(T(0.75), T(1.75))
+        s = interval(T, 0.75, 1.75)
         example = Slopes{T}[]
-        push!(example, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(T(-2.8), T(0.1))))
-        push!(example, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(T(-44), T(38.5))))
-        push!(example, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(T(-0.16), T(0.44))))
-        push!(example, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(T(0.03), T(0.33))))
-        push!(example, Slopes(x->(exp(x^2)), s, mid(s), interval(T(6.03), T(33.23))))
-        push!(example, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(T(-39), T(65.56))))
-        push!(example, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(T(-146.9), T(67.1))))
-        push!(example, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(T(1), T(2))))
-        push!(example, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(T(1.36), T(∞))))
+        push!(example, Slopes(x->((x + sin(x)) * exp(-x^2)), s, mid(s), interval(T, -2.8, 0.1)))
+        push!(example, Slopes(x->(x^4 - 10x^3 + 35x^2 - 50x + 24), s, mid(s), interval(T, -44, 38.5)))
+        push!(example, Slopes(x->((log(x + 1.25) - 0.84x) ^ 2), s, mid(s), interval(T, -0.16, 0.44)))
+        push!(example, Slopes(x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)), s, mid(s), interval(T, 0.03, 0.33)))
+        push!(example, Slopes(x->(exp(x^2)), s, mid(s), interval(T, 6.03, 33.23)))
+        push!(example, Slopes(x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)), s, mid(s), interval(T, -39, 65.56)))
+        push!(example, Slopes(x->(x^6 - 15x^4 + 27x^2 + 250), s, mid(s), interval(T, -146.9, 67.1)))
+        push!(example, Slopes(x->(atan(cos(tan(x)))), s, mid(s), interval(T, 1, 2)))
+        push!(example, Slopes(x->(asin(cos(acos(sin(x))))), s, mid(s), interval(T, 1.36, ∞)))
 
         for i in 1:length(example)
             @test slope(example[i].f, example[i].x, example[i].c) ⊆ example[i].sol