Skip to content

Commit

Permalink
Working for bare intervals
Browse files Browse the repository at this point in the history
  • Loading branch information
dpsanders committed Jul 28, 2024
1 parent 18cf52d commit 86ad421
Show file tree
Hide file tree
Showing 14 changed files with 134 additions and 98 deletions.
3 changes: 1 addition & 2 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
4 changes: 2 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
26 changes: 19 additions & 7 deletions src/IntervalContractors.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
__precompile__(true)

module IntervalContractors

export plus_rev, minus_rev, inv_rev,
Expand All @@ -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")
Expand Down
54 changes: 27 additions & 27 deletions src/arithmetic.jl
Original file line number Diff line number Diff line change
@@ -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`.
Expand All @@ -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)
Expand All @@ -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`.
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -54,24 +54,24 @@ 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))

# 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)
Expand All @@ -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)
Expand All @@ -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⁻¹
Expand All @@ -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)

Expand All @@ -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.
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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`.
Expand All @@ -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)

Expand All @@ -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.
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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)
4 changes: 2 additions & 2 deletions src/decorated.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
16 changes: 8 additions & 8 deletions src/exponential.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)

Expand All @@ -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)

Expand All @@ -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)

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/extrema.jl
Original file line number Diff line number Diff line change
@@ -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);
Expand All @@ -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);
Expand Down
Loading

0 comments on commit 86ad421

Please sign in to comment.