From 86ad421f81bf2f47f65ff9361d69ccc5cd09722e Mon Sep 17 00:00:00 2001 From: David Sanders Date: Sun, 28 Jul 2024 18:46:19 -0500 Subject: [PATCH] Working for bare intervals --- NEWS.md | 3 +-- README.md | 4 +-- src/IntervalContractors.jl | 26 +++++++++++++----- src/arithmetic.jl | 54 +++++++++++++++++++------------------- src/decorated.jl | 4 +-- src/exponential.jl | 16 +++++------ src/extrema.jl | 4 +-- src/hyperbolic.jl | 12 ++++----- src/inverse_hyperbolic.jl | 6 ++--- src/inverse_trig.jl | 26 +++++++++--------- src/powers.jl | 2 +- src/transformations.jl | 25 +++++++++--------- src/trig.jl | 31 ++++++++++++++-------- test/runtests.jl | 19 ++++++++++++-- 14 files changed, 134 insertions(+), 98 deletions(-) diff --git a/NEWS.md b/NEWS.md index 9e8d53a..d83f07a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,5 +8,4 @@ - Adds a dependency on the separate `IntervalBoxes.jl` package (to replace the functionality removed from `IntervalArithmetic.jl`) - - Although the user-facing API of `IntervalContractors.jl` has *not changed* - in this release, the extra dependency warrants a minor-version bump + - `Interval` now refers to a *decorated* interval \ No newline at end of file diff --git a/README.md b/README.md index 1832b86..f5173ab 100644 --- a/README.md +++ b/README.md @@ -18,10 +18,10 @@ The reverse function of a function `f` calculates the (interval hull of) its inverse function, `f⁻¹`. -For example, `sin_rev(Y::Interval, X::Interval)` calculates the (interval hull of) those `x ∈ X` such that `sin(x) ∈ Y`. This can also be thought of as an inverse function, calculating `X_new := sin⁻¹(Y) ⊓ X`. +For example, `sin_rev(Y::IntervalType, X::IntervalType)` calculates the (interval hull of) those `x ∈ X` such that `sin(x) ∈ Y`. This can also be thought of as an inverse function, calculating `X_new := sin⁻¹(Y) ⊓ X`. The return value is `(Y, X_new)`. -Functions such as `mul_rev(C::Interval, A::Interval, B::Interval)` take three arguments, and correspond to `C = A * B`; they return `(C, A_new, B_new)`, with `A_new` and `B_new` similarly defined to be the corresponding inverse images of the multiplication operator in each component. +Functions such as `mul_rev(C::IntervalType, A::IntervalType, B::IntervalType)` take three arguments, and correspond to `C = A * B`; they return `(C, A_new, B_new)`, with `A_new` and `B_new` similarly defined to be the corresponding inverse images of the multiplication operator in each component. ### Contractors diff --git a/src/IntervalContractors.jl b/src/IntervalContractors.jl index 2f9acca..58fe3e3 100644 --- a/src/IntervalContractors.jl +++ b/src/IntervalContractors.jl @@ -1,5 +1,3 @@ -__precompile__(true) - module IntervalContractors export plus_rev, minus_rev, inv_rev, @@ -18,11 +16,25 @@ export plus_rev, minus_rev, inv_rev, using IntervalArithmetic, IntervalArithmetic.Symbols using IntervalBoxes -const half_pi = ExactReal(0.5) * interval(pi) -const two_pi = ExactReal(2.0) * interval(pi) -# -# Base.:⊔(f::Function, g::Function) = X -> ( f(X) ⊔ g(X) ) -# Base.:⊓(f::Function, g::Function) = X -> ( f(X) ⊓ g(X) ) # or f(g(X)) for contractors +const IntervalType{T} = Union{Interval{T}, BareInterval{T}} + +# @generated +# half_pi(::Type{T}) where {T <: IntervalType} = :(ExactReal(0.5) * convert($T, ExactReal(pi))) +# @generated +# two_pi(::Type{T}) where {T <: IntervalType} = :(ExactReal(2.0) * convert($T, ExactReal(pi))) + +@generated function half_pi(x::T) where {T <: IntervalType} + return ExactReal(0.5) * convert(T, ExactReal(pi)) +end + +@generated function two_pi(x::T) where {T <: IntervalType} + return ExactReal(2.0) * convert(T, ExactReal(pi)) +end + +@generated function pi_interval(x::T) where {T <: IntervalType} + return convert(T, ExactReal(pi)) +end + include("arithmetic.jl") include("transformations.jl") diff --git a/src/arithmetic.jl b/src/arithmetic.jl index bf51b45..85a2ffb 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -1,6 +1,6 @@ """ - plus_rev(a::Interval, b::Interval[, c::Interval]) + plus_rev(a::IntervalType, b::IntervalType[, c::IntervalType]) Reverse addition. Calculates the preimage of `a = b + c` for `b` and `c`. @@ -12,7 +12,7 @@ The triplet `(a, b_new, c_new)` where - `b_new` is the interval hull of the set ``{x ∈ b : ∃ y ∈ c, x + y ∈ a}`` - `c_new` is the interval hull of the set ``{y ∈ c : ∃ x ∈ b, x + y ∈ a}`` """ -function plus_rev(a::Interval, b::Interval, c::Interval) # a = b + c +function plus_rev(a::IntervalType, b::IntervalType, c::IntervalType) # a = b + c # a = a ⊓ (b + c) # add this line for plus contractor (as opposed to reverse function) b_new = b ⊓ (a - c) c_new = c ⊓ (a - b) @@ -23,7 +23,7 @@ end plus_rev(a,b,c) = plus_rev(promote(a,b,c)...) """ - minus_rev(a::Interval, b::Interval[, c::Interval]) + minus_rev(a::IntervalType, b::IntervalType[, c::IntervalType]) Reverse subtraction. Calculates the preimage of `a = b - c` for `b` and `c`. @@ -35,7 +35,7 @@ The triplet `(a, b_new, c_new)` where - `b_new` is the interval hull of the set ``{x ∈ b : ∃ y ∈ c, x - y ∈ a}`` - `c_new` is the interval hull of the set ``{y ∈ c : ∃ x ∈ b, x - y ∈ a}`` """ -function minus_rev(a::Interval, b::Interval, c::Interval) # a = b - c +function minus_rev(a::IntervalType, b::IntervalType, c::IntervalType) # a = b - c b_new = b ⊓ (a + c) c_new = c ⊓ (b - a) @@ -45,7 +45,7 @@ end minus_rev(a,b,c) = minus_rev(promote(a,b,c)...) -function minus_rev(a::Interval, b::Interval) # a = -b +function minus_rev(a::IntervalType, b::IntervalType) # a = -b b_new = b ⊓ (-a) return (a, b_new) end @@ -54,7 +54,7 @@ end """ Reverse multiplication """ -function mul_rev(a::Interval, b::Interval, c::Interval) # a = b * c +function mul_rev(a::IntervalType, b::IntervalType, c::IntervalType) # a = b * c # ((0.0 ∉ a) || (0.0 ∉ b)) && (c = c ⊓ (a / b)) # ((0.0 ∉ a) || (0.0 ∉ c)) && (b = b ⊓ (a / c)) @@ -62,16 +62,16 @@ function mul_rev(a::Interval, b::Interval, c::Interval) # a = b * c # a = a ⊓ (b * c) # ? if in_interval(0.0, b) - temp = c .⊓ extended_div(a, b) - c′ = hull(temp[1], temp[2]) + temp = extended_div(a, b) + c′ = hull(c ⊓ temp[1], c ⊓ temp[2]) else c′ = c ⊓ (a / b) end if in_interval(0.0, c) - temp = b .⊓ extended_div(a, c) - b′ = hull(temp[1], temp[2]) + temp = extended_div(a, c) + b′ = hull(b ⊓ temp[1], b ⊓ temp[2]) else b′ = b ⊓ (a / c) @@ -85,7 +85,7 @@ mul_rev(a,b,c) = mul_rev(promote(a,b,c)...) """ Reverse division """ -function div_rev(a::Interval, b::Interval, c::Interval) # a = b / c +function div_rev(a::IntervalType, b::IntervalType, c::IntervalType) # a = b / c b = b ⊓ (a * c) c = c ⊓ (b / a) @@ -96,7 +96,7 @@ end div_rev(a,b,c) = div_rev(promote(a,b,c)...) """ - inv_rev(a::Interval, b::Interval) + inv_rev(a::IntervalType, b::IntervalType) Reverse inverse. Calculates the interval hull of the preimage of a = b⁻¹ @@ -107,7 +107,7 @@ Pair `(a, b_new)` where - `a` is unchanged - `b_new` is the interval hull of the set ``{x ∈ b : x⁻¹ ∈ a}`` """ -function inv_rev(a::Interval, b::Interval) # a = inv(b) +function inv_rev(a::IntervalType, b::IntervalType) # a = inv(b) b_new = b ⊓ inv(a) @@ -117,7 +117,7 @@ end inv_rev(a,b) = inv_rev(promote(a,b)...) """ - power_rev(a::Interval, b::Interval, n::Integer) + power_rev(a::IntervalType, b::IntervalType, n::Integer) Reverse power. Calculates the preimage of `a = bⁿ`. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -129,7 +129,7 @@ The triplet `(a, b_new, n)` where - `a` and `n` are unchanged - `b_new` is the interval hull of the set ``{x ∈ b : xⁿ ∈ a}`` """ -function power_rev(a::Interval{T}, b::Interval{T}, n::Integer) where T # a = b^n, log(a) = inf(n)g(b), b = a^(1/n) +function power_rev(a::IntervalType{T}, b::IntervalType{T}, n::Integer) where T # a = b^n, log(a) = inf(n)g(b), b = a^(1/n) if iszero(n) in_interval(1.0, a) && return (a, entireinterval(T) ⊓ b, n) @@ -161,9 +161,9 @@ function power_rev(a::Interval{T}, b::Interval{T}, n::Integer) where T # a = b^ return (a, b, n) end -power_rev(a::Interval{T}, n::Integer) where {T} = power_rev(a, entireinterval(T), n) +power_rev(a::IntervalType{T}, n::Integer) where {T} = power_rev(a, entireinterval(T), n) -function power_rev(a::Interval, b::Interval, c::Interval) # a = b^c +function power_rev(a::IntervalType, b::IntervalType, c::IntervalType) # a = b^c if isinteger(c) temp = power_rev(a, b, Int(inf(c))) # use version with integer @@ -180,7 +180,7 @@ power_rev(a, b, c) = power_rev(promote(a, b, c)...) """ - sqrt_rev(a::Interval, b::Interval) + sqrt_rev(a::IntervalType, b::IntervalType) Reverse square root. Calculates the preimage of `a = √b`. @@ -191,7 +191,7 @@ The pair `(a, b_new)` where - `a` is unchanged - `b_new` is the interval hull of the set ``{x ∈ b : √x ∈ a}`` """ -function sqrt_rev(a::Interval, b::Interval) # a = sqrt(b) +function sqrt_rev(a::IntervalType, b::IntervalType) # a = sqrt(b) b_new = b ⊓ (a^2) @@ -204,7 +204,7 @@ sqrt_rev(a,b) = sqrt_rev(promote(a,b)...) # IEEE-1788 style """ - sqr_rev(c::Interval[, x::Interval]) + sqr_rev(c::IntervalType[, x::IntervalType]) Reverse square. Calculates the preimage of `a = x²`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -227,7 +227,7 @@ function sqr_rev(c, x = entireinterval(c)) # c = x^2; refine x end """ - abs_rev(c::Interval[, x::Interval]) + abs_rev(c::IntervalType[, x::IntervalType]) Reverse absolute value. Calculates the preimage of `a = |x|`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -252,7 +252,7 @@ end """ Reverse sign """ -function sign_rev(a::Interval, b::Interval) # a = sqrt(b) +function sign_rev(a::IntervalType, b::IntervalType) # a = sqrt(b) (a == 1.0) && b = b ⊓ (interval(0, Inf)) (a == 0.0) && b = b ⊓ (0.interval(0, 0).0) @@ -266,7 +266,7 @@ sign_rev(a,b) = sign_rev(promote(a,b)...) ## IEEE-1788 versions: """ - mul_rev_IEEE1788(b::Interval, c::Interval[, x::Interval]) + mul_rev_IEEE1788(b::IntervalType, c::IntervalType[, x::IntervalType]) Reverse multiplication. Computes the preimage of ``c = x * b`` with respect to `x`. If `x` is not provided, then by default ``[-Inf, Inf]`` is used. @@ -279,7 +279,7 @@ See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. mul_rev_IEEE1788(b, c, x = entireinterval(b)) = mul_rev(c, x, b)[2] """ - pow_rev1(b::Interval, c::Interval[, x::Interval]) + pow_rev1(b::IntervalType, c::IntervalType[, x::IntervalType]) Reverse power 1. Computes the preimage of ``c=xᵇ`` with respect to `x`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used.. See section 10.5.4 of the @@ -294,7 +294,7 @@ function pow_rev1(b, c, x) # c = x^b end """ - pow_rev2(b::Interval, c::Interval[, x::Interval]) + pow_rev2(b::IntervalType, c::IntervalType[, x::IntervalType]) Reverse power 2. Computes the preimage of ``c = aˣ`` with respect to `x`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -308,7 +308,7 @@ function pow_rev2(a, c, x) # c = a^x end """ - mul_rev_to_pair(b::Interval, c::Interval) + mul_rev_to_pair(b::IntervalType, c::IntervalType) Computes the division c / b, but returns a pair of intervals instead of a single interval. If the set corresponding to c / b is composed of two disjoint intervals, then it returns the @@ -326,4 +326,4 @@ julia> mul_rev_to_pair(interval(1, 2), interval(3, 4)) ([1.5, 4], emptyinterval()) """ -mul_rev_to_pair(b::Interval, c::Interval) = extended_div(c, b) +mul_rev_to_pair(b::IntervalType, c::IntervalType) = extended_div(c, b) diff --git a/src/decorated.jl b/src/decorated.jl index a920015..121b1d9 100644 --- a/src/decorated.jl +++ b/src/decorated.jl @@ -8,7 +8,7 @@ for op in (:sqr_rev, :abs_rev, :sin_rev, :cos_rev, :tan_rev, :cosh_rev, :sinh_re return (Decoratedinterval(bare[1], trv), Decoratedinterval(bare[2], trv)) end end - @eval $op(a::Interval{T}) where T = $op(a, entireinterval(T)) + @eval $op(a::IntervalType{T}) where T = $op(a, entireinterval(T)) @eval $op(a::DecoratedInterval{T}) where T = $op(a, entiredecorated(T)) end @@ -29,7 +29,7 @@ for op in (:mul_rev_IEEE1788, :pow_rev1, :pow_rev2) end end - @eval $op(a::Interval{T}, b::Interval{T}) where T = $op(a, b, entireinterval(T)) + @eval $op(a::IntervalType{T}, b::IntervalType{T}) where T = $op(a, b, entireinterval(T)) @eval $op(a::DecoratedInterval{T}, b::DecoratedInterval{T}) where T = $op(a, b, entiredecorated(T)) end \ No newline at end of file diff --git a/src/exponential.jl b/src/exponential.jl index 2c373f9..483b1fe 100644 --- a/src/exponential.jl +++ b/src/exponential.jl @@ -10,7 +10,7 @@ end """ Reverse function for `exp`. """ -function exp_rev(y::Interval, x::Interval) +function exp_rev(y::IntervalType, x::IntervalType) y_new = y ⊓ (interval(0, Inf)) x_new = x ⊓ log(y) return y_new, x_new @@ -19,7 +19,7 @@ end """ Reverse function for `exp2`. """ -function exp2_rev(y::Interval, x::Interval) +function exp2_rev(y::IntervalType, x::IntervalType) y_new = y ⊓ (interval(0, Inf)) x_new = x ⊓ log2(y) @@ -29,7 +29,7 @@ end """ Reverse function for `exp10`. """ -function exp10_rev(y::Interval, x::Interval) +function exp10_rev(y::IntervalType, x::IntervalType) y_new = y ⊓ (interval(0, Inf)) x_new = x ⊓ log10(y) @@ -39,7 +39,7 @@ end """ Reverse function for `expm1`. """ -function expm1_rev(y::Interval, x::Interval) +function expm1_rev(y::IntervalType, x::IntervalType) y_new = y ⊓ (interval(-1, Inf)) x_new = x ⊓ log1p(y) @@ -59,7 +59,7 @@ end """ Reverse function for `log`: ``y = \\log(x)`` """ -function log_rev(y::Interval, x::Interval) +function log_rev(y::IntervalType, x::IntervalType) x_new = x ⊓ exp(y) return y, x_new @@ -68,7 +68,7 @@ end """ Reverse function for `log2`: ``y = \\log2(x)`` """ -function log2_rev(y::Interval, x::Interval) +function log2_rev(y::IntervalType, x::IntervalType) x_new = x ⊓ exp2(y) return y, x_new @@ -78,7 +78,7 @@ end """ Reverse function for `log10`: ``y = \\log10(x)`` """ -function log10_rev(y::Interval, x::Interval) +function log10_rev(y::IntervalType, x::IntervalType) x_new = x ⊓ exp10(y) return y, x_new @@ -87,7 +87,7 @@ end """ Reverse function for `log1p`: ``y = \\log1p(x)`` """ -function log1p_rev(y::Interval, x::Interval) +function log1p_rev(y::IntervalType, x::IntervalType) x_new = x ⊓ expm1(y) return y, x_new diff --git a/src/extrema.jl b/src/extrema.jl index f08684e..30263fc 100644 --- a/src/extrema.jl +++ b/src/extrema.jl @@ -1,7 +1,7 @@ """ Reverse max """ -function max_rev(a::Interval, b::Interval, c::Interval) # a = max(b,c) +function max_rev(a::IntervalType, b::IntervalType, c::IntervalType) # a = max(b,c) B_lo = inf(b); B_hi = sup(b); C_lo = inf(c); C_hi = sup(c); @@ -24,7 +24,7 @@ max_rev(a,b,c) = max_rev(promote(a,b,c)...) """ Reverse min """ -function min_rev(a::Interval, b::Interval, c::Interval) +function min_rev(a::IntervalType, b::IntervalType, c::IntervalType) B_lo = inf(b); B_hi = sup(b); C_lo = inf(c); C_hi = sup(c); diff --git a/src/hyperbolic.jl b/src/hyperbolic.jl index 53dcfcf..b795f4d 100644 --- a/src/hyperbolic.jl +++ b/src/hyperbolic.jl @@ -1,5 +1,5 @@ """ - sinh_rev(c::Interval[, x::Interval]) + sinh_rev(c::IntervalType[, x::IntervalType]) Reverse hyperbolic sine. Calculates the preimage of `a = sinh(x)`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -11,14 +11,14 @@ The pair `(c, x_new)` where - `c` is unchanged - `x_new` is the interval hull of the set ``{x ∈ b : sinh(x) ∈ a}`` """ -function sinh_rev(y::Interval, x::Interval = entireinterval(y)) +function sinh_rev(y::IntervalType, x::IntervalType = entireinterval(y)) x = x ⊓ asinh(y) return y, x end """ - cosh_rev(c::Interval[, x::Interval]) + cosh_rev(c::IntervalType[, x::IntervalType]) Reverse square root. Calculates the preimage of `a = cosh(x)`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -30,7 +30,7 @@ The pair `(c, x_new)` where - `c` is unchanged - `x_new` is the interval hull of the set ``{x ∈ b : cosh(x) ∈ a}`` """ -function cosh_rev(y::Interval, x::Interval = entireinterval(y)) +function cosh_rev(y::IntervalType, x::IntervalType = entireinterval(y)) y_new = y ⊓ interval(1.,Inf) x = (x ⊓ acosh(y)) ⊔ (x ⊓ -acosh(y)) @@ -38,7 +38,7 @@ function cosh_rev(y::Interval, x::Interval = entireinterval(y)) end """ - tanh_rev(c::Interval[, x::Interval]) + tanh_rev(c::IntervalType[, x::IntervalType]) Reverse square root. Calculates the preimage of `a = tanh(x)`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -50,7 +50,7 @@ The pair `(c, x_new)` where - `c` is unchanged - `x_new` is the interval hull of the set ``{x ∈ b : tanh(x) ∈ a}`` """ -function tanh_rev(y::Interval, x::Interval = entireinterval(y)) +function tanh_rev(y::IntervalType, x::IntervalType = entireinterval(y)) y_new = y ⊓ interval(-1.,1.) x = x ⊓ atanh(y) diff --git a/src/inverse_hyperbolic.jl b/src/inverse_hyperbolic.jl index c8ae4d9..e1c182d 100644 --- a/src/inverse_hyperbolic.jl +++ b/src/inverse_hyperbolic.jl @@ -1,7 +1,7 @@ """ Reverse function for `asinh`. """ -function asinh_rev(y::Interval, x::Interval = entireinterval(y)) +function asinh_rev(y::IntervalType, x::IntervalType = entireinterval(y)) x = x ⊓ sinh(y) return y, x @@ -10,7 +10,7 @@ end """ Reverse function for `acosh`. """ -function acosh_rev(y::Interval, x::Interval = entireinterval(y)) +function acosh_rev(y::IntervalType, x::IntervalType = entireinterval(y)) y_new = y ⊓ interval(0.0,Inf) x = x ⊓ cosh(y) @@ -20,7 +20,7 @@ end """ Reverse function for `atanh`. """ -function atanh_rev(y::Interval, x::Interval = entireinterval(y)) +function atanh_rev(y::IntervalType, x::IntervalType = entireinterval(y)) x = x ⊓ tanh(y) return y, x diff --git a/src/inverse_trig.jl b/src/inverse_trig.jl index 0b9680e..7ea40d1 100644 --- a/src/inverse_trig.jl +++ b/src/inverse_trig.jl @@ -1,7 +1,7 @@ function asin!(X::IntervalBox) x, y = X - h = inf(half_pi) + h = inf(half_pi(y)) y_new = y ⊓ interval(-h, h) # range of asin x_new = sin(y_new) @@ -11,9 +11,8 @@ end """ Reverse `asin`. """ -function asin_rev(y::Interval, x::Interval) # y = asin(x) - - h = inf(half_pi) +function asin_rev(y::IntervalType, x::IntervalType = entireinterval(y)) # y = asin(x) + h = inf(half_pi(y)) y_new = y ⊓ interval(-h, h) # range of asin x_new = sin(y_new) @@ -24,22 +23,23 @@ end """ Reverse `acos`. """ -function acos_rev(y::Interval, x::Interval) - y_new = y ⊓ interval(0.0, sup(two_pi)) - x_new = x ⊓ cos(y_new) +function acos_rev(y::IntervalType, x::IntervalType = entireinterval(y)) + y_new = y ⊓ interval(0.0, sup(two_pi(y))) + x_new = x ⊓ cos(y_new) - return y_new, x_new + return y_new, x_new end """ - atan_rev(y::Interval, x::Interval) + atan_rev(y::IntervalType, x::IntervalType) Inverse of `y = atan(x)`. Returns the new `y` and `x`. """ -function atan_rev(y::Interval, x::Interval) - y_new = y ⊓ interval(-sup(half_pi), sup(half_pi)) - x_new = x ⊓ tan(y_new) +function atan_rev(y::IntervalType, x::IntervalType = entireinterval(y)) + h = sup(half_pi(y)) + y_new = y ⊓ interval(-h, h) + x_new = x ⊓ tan(y_new) - return y_new, x_new + return y_new, x_new end diff --git a/src/powers.jl b/src/powers.jl index 1384e9b..20db535 100644 --- a/src/powers.jl +++ b/src/powers.jl @@ -19,7 +19,7 @@ function square_pos(X::IntervalBox) return IntervalBox(x, y) end -square_neg = symmetrise(square_pos, reflect_x(0.0)) +square_neg = symmetrise(square_pos, reflect_x(zero)) square!(X::IntervalBox) = square_pos(X) ⊔ square_neg(X) diff --git a/src/transformations.jl b/src/transformations.jl index fb35b0c..c5aab2e 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -1,27 +1,28 @@ """ - integer_contractor(x::Interval) + integer_contractor(x::IntervalType) Return the integers enclosed in the interval `x`. """ -function integer_contractor(x::Interval) - a = floor(inf(x))+1 +function integer_contractor(x::IntervalType) + a = floor(inf(x)) + 1 b = floor(sup(x)) a > b && return emptyinterval(x) - return interval(a, b) + return _build_interval(x, a, b) end ### Transformations on IntervalBoxes -"""Reflect in mirror at position x_mirror""" +"""Reflect in mirror at position x_mirror +x_mirror is a function that returns an interval giving the position of the mirror.""" function reflect_x(x_mirror) X -> begin x, y = X - x = 2*x_mirror - x + x = (ExactReal(2) * x_mirror(x)) - x return IntervalBox(x, y) end @@ -62,23 +63,23 @@ symmetrise(C, op) = op ∘ C ∘ op -"Periodize the contractor C" +"Periodize the contractor C. period is a function that returns an interval giving the period" function periodise(C, period) X -> begin x, y = X - x2 = entireinterval() + x2 = entireinterval(x) x2, y = C(IntervalBox(x2, y)) - isempty(IntervalBox(x2, y)) && return(IntervalBox(emptyinterval(), emptyinterval())) + isempty(IntervalBox(x2, y)) && return(IntervalBox(emptyinterval(x), emptyinterval(x))) # periods where the periodization of x intersects with x2: - periods = integer_contractor((x - x2) / period) + periods = integer_contractor((x - x2) / period(x)) - isempty_interval(periods) && return(IntervalBox(emptyinterval(), emptyinterval())) + isempty_interval(periods) && return(IntervalBox(emptyinterval(x), emptyinterval(x))) - x3 = x2 + periods*period + x3 = x2 + periods * period(x) x = x ⊓ x3 return IntervalBox(x, y) diff --git a/src/trig.jl b/src/trig.jl index 5a528ec..356d616 100644 --- a/src/trig.jl +++ b/src/trig.jl @@ -8,8 +8,10 @@ function sin_main(X::IntervalBox) x, y = X - x_range = interval(-sup(half_pi), sup(half_pi)) - y_range = interval(-1, 1) + h = sup(half_pi(x)) + + x_range = _build_interval(x, -h, h) + y_range = _build_interval(x, -1, 1) x = x ⊓ x_range y = y ⊓ y_range @@ -37,7 +39,7 @@ sin!(X::IntervalBox) = periodise(sin_main, two_pi)(X) ⊔ periodise(sin_reverse, # Reverse function for sin; does not alter y """ - sin_rev(c::Interval[, x::Interval]) + sin_rev(c::IntervalType[, x::IntervalType]) Reverse sine. Calculates the preimage of `a = sin(x)`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -49,7 +51,7 @@ The pair `(c, x_new)` where - `c` is unchanged - `x_new` is the interval hull of the set ``{x ∈ b : sin(x) ∈ a}`` """ -function sin_rev(y::Interval, x::Interval = entireinterval(y)) +function sin_rev(y::IntervalType, x::IntervalType = entireinterval(y)) X = IntervalBox(x, y) @@ -67,8 +69,8 @@ function cos_main(X::IntervalBox) x, y = X - x_range = interval(0, inf(interval(π))) - y_range = interval(-1, 1) + x_range = _build_interval(x, 0, inf(pi_interval(x))) + y_range = _build_interval(x, -1, 1) x = x ⊓ x_range y = y ⊓ y_range @@ -97,7 +99,7 @@ cos!(X::IntervalBox) = periodise(cos_main, two_pi)(X) ⊔ periodise(cos_reverse, # Reverse function for cos; does not alter y """ - cos_rev(c::Interval[, x::Interval]) + cos_rev(c::IntervalType[, x::IntervalType]) Reverse cosine. Calculates the preimage of `a = cos(x)`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -109,7 +111,7 @@ The pair `(c, x_new)` where - `c` is unchanged - `x_new` is the interval hull of the set ``{x ∈ b : cos(x) ∈ a}`` """ -function cos_rev(y::Interval, x::Interval = entireinterval(y)) +function cos_rev(y::IntervalType, x::IntervalType = entireinterval(y)) X = IntervalBox(x, y) @@ -126,7 +128,9 @@ function tan_main(X::IntervalBox) x, y = X - x_range = interval(-sup(half_pi), sup(half_pi)) + h = sup(half_pi) + + x_range = _build_interval(x, -h, h) x = x ⊓ x_range @@ -142,7 +146,7 @@ end tan!(X::IntervalBox) = periodise(tan_main, interval(π))(X) """ - tan_rev(c::Interval[, x::Interval]) + tan_rev(c::IntervalType[, x::IntervalType]) Reverse tangent. Calculates the preimage of `a = tan(x)`. If `x` is not provided, then byt default ``[-Inf, Inf]`` is used. See section 10.5.4 of the IEEE 1788-2015 standard for interval arithmetic. @@ -154,7 +158,7 @@ The pair `(c, x_new)` where - `c` is unchanged - `x_new` is the interval hull of the set ``{x ∈ b : tan(x) ∈ a}`` """ -function tan_rev(y::Interval, x::Interval = entireinterval(y)) +function tan_rev(y::IntervalType, x::IntervalType = entireinterval(y)) X = IntervalBox(x, y) @@ -162,3 +166,8 @@ function tan_rev(y::Interval, x::Interval = entireinterval(y)) return X_new[2], X_new[1] # return in order y, x end + + +# build an interval of the corresponding type: +_build_interval(x::Interval, a, b) = interval(a, b) +_build_interval(x::BareInterval, a, b) = bareinterval(a, b) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index edc9a05..90e3932 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,8 +8,11 @@ orig_power_mode = IntervalArithmetic.power_mode() # set power mode to "slow" (using MPFR for correct rounding): IntervalArithmetic.power_mode() = IntervalArithmetic.PowerMode{:slow}() -eq(a, b) = isequal_interval(bareinterval(a), bareinterval(b)) -approx_eq(x::Interval, y::Interval) = +eq(a::IntervalType, b::IntervalType) = isequal_interval(bareinterval(a), bareinterval(b)) + +eq(a::Tuple, b::Tuple) = all(eq.(a, b)) + +approx_eq(x::IntervalType, y::IntervalType) = isapprox(inf(x), inf(y), atol=1e-4) && isapprox(sup(x), sup(y), atol=1e-4) @testset "IntervalContractors tests" begin @@ -38,5 +41,17 @@ approx_eq(x::Interval, y::Interval) = end end +@testset "Bare intervals" begin + x = bareinterval(0, 1) + y = bareinterval(0, 2) + z = bareinterval(0, 3) + + @test eq(plus_rev(x, y, z), + (bareinterval(0, 1), bareinterval(0, 1), bareinterval(0, 1))) + + @test eq(sin_rev(bareinterval(0..1), bareinterval(3..4)), + (bareinterval(0, 1), bareinterval(3, 3.1415926535897936))) +end + # reset power mode: IntervalArithmetic.power_mode() = orig_power_mode