From 1aefcbf43fb7ce9ebde2f8b01ce97c8c481ab249 Mon Sep 17 00:00:00 2001 From: schillic Date: Sun, 27 Aug 2023 11:56:09 +0200 Subject: [PATCH 1/2] partial update to IntervalArithmetic v0.21 --- .gitignore | 3 ++- benchmark/benchmarks.jl | 6 +++--- perf/linear_eq.jl | 4 ++-- src/IntervalRootFinding.jl | 9 ++++++--- src/contractors.jl | 8 ++++---- src/krawczyk.jl | 10 +++++----- src/linear_eq.jl | 4 ++-- src/newton.jl | 10 +++++----- src/quadratic.jl | 14 +++++++------- src/roots.jl | 2 +- test/bisection.jl | 20 ++++++++++---------- test/findroots.jl | 4 ++-- test/linear_eq.jl | 4 ++-- test/newton1d.jl | 1 - test/roots.jl | 1 - 15 files changed, 51 insertions(+), 49 deletions(-) diff --git a/.gitignore b/.gitignore index 4feed2c0..4cffedf8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ docs/build/ docs/site/ benchmark/tune.json -Manifest.toml \ No newline at end of file +Manifest.toml +.vscode diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index d0d7ff26..a3d959e9 100644 --- a/benchmark/benchmarks.jl +++ b/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() diff --git a/perf/linear_eq.jl b/perf/linear_eq.jl index a418bfba..ae25d51e 100644 --- a/perf/linear_eq.jl +++ b/perf/linear_eq.jl @@ -2,7 +2,7 @@ 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 @@ -10,7 +10,7 @@ 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 diff --git a/src/IntervalRootFinding.jl b/src/IntervalRootFinding.jl index fe414297..7012bd5b 100644 --- a/src/IntervalRootFinding.jl +++ b/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 diff --git a/src/contractors.jl b/src/contractors.jl index 5fb6dd95..fbc4e467 100644 --- a/src/contractors.jl +++ b/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 diff --git a/src/krawczyk.jl b/src/krawczyk.jl index 206fc34f..64068659 100644 --- a/src/krawczyk.jl +++ b/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...) diff --git a/src/linear_eq.jl b/src/linear_eq.jl index 22e0f94a..a8928759 100644 --- a/src/linear_eq.jl +++ b/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) diff --git a/src/newton.jl b/src/newton.jl index 83f81805..47d7fe02 100644 --- a/src/newton.jl +++ b/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...) diff --git a/src/quadratic.jl b/src/quadratic.jl index c7fe6f6f..4c9e75f3 100644 --- a/src/quadratic.jl +++ b/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 diff --git a/src/roots.jl b/src/roots.jl index 4e87dd53..d9e8c442 100644 --- a/src/roots.jl +++ b/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) diff --git a/test/bisection.jl b/test/bisection.jl index a59a4b4a..9654f395 100644 --- a/test/bisection.jl +++ b/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 diff --git a/test/findroots.jl b/test/findroots.jl index 602ff287..4a7192c0 100644 --- a/test/findroots.jl +++ b/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 diff --git a/test/linear_eq.jl b/test/linear_eq.jl index 8b720998..fd392e46 100644 --- a/test/linear_eq.jl +++ b/test/linear_eq.jl @@ -3,7 +3,7 @@ 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 @@ -11,7 +11,7 @@ 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 diff --git a/test/newton1d.jl b/test/newton1d.jl index 2c11c0bb..e4b95fbe 100644 --- a/test/newton1d.jl +++ b/test/newton1d.jl @@ -1,4 +1,3 @@ - using IntervalArithmetic, IntervalRootFinding using ForwardDiff using Test diff --git a/test/roots.jl b/test/roots.jl index 3d7fb6bf..1f20a466 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -1,4 +1,3 @@ - using IntervalArithmetic, IntervalRootFinding, StaticArrays, ForwardDiff using Test From ddddd2de967c3a2e020802a7f1fbefe081ca0fcb Mon Sep 17 00:00:00 2001 From: schillic Date: Mon, 28 Aug 2023 22:23:43 +0200 Subject: [PATCH 2/2] more fixes --- Project.toml | 2 +- src/complex.jl | 4 +++- src/linear_eq.jl | 3 +-- src/slopes.jl | 2 +- test/slopes.jl | 20 ++++++++++---------- 5 files changed, 16 insertions(+), 15 deletions(-) diff --git a/Project.toml b/Project.toml index 3e8cbcc4..55d991aa 100644 --- a/Project.toml +++ b/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" diff --git a/src/complex.jl b/src/complex.jl index 567345d7..8360eed9 100644 --- a/src/complex.jl +++ b/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") diff --git a/src/linear_eq.jl b/src/linear_eq.jl index a8928759..bc6e042a 100644 --- a/src/linear_eq.jl +++ b/src/linear_eq.jl @@ -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 diff --git a/src/slopes.jl b/src/slopes.jl index 82119666..319e9eb3 100644 --- a/src/slopes.jl +++ b/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) diff --git a/test/slopes.jl b/test/slopes.jl index 5c5d51a6..4feb2203 100644 --- a/test/slopes.jl +++ b/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