Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 32 additions & 3 deletions ext/IntervalArithmeticForwardDiffExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,20 @@ Base.:(==)(x::Dual, y::Interval) = value(x) == y
Base.:<(x::Interval, y::Dual) = x < value(y)
Base.:<(x::Dual, y::Interval) = value(x) < y

# ForwardDiff's `^` methods use `iszero(partials(x))`, which is undecidable for
# non-thin intervals. Use a recursive thin-zero test instead. `NestedInterval`
# supports intervals nested in up to four `Dual` layers.
const NestedInterval = let
U = Interval
for _ in 1:4
U = Union{U, Dual{T,<:U} where {T}}
end
U
end
_isthinzero(x::Interval) = isthinzero(x)
_isthinzero(x::Real) = iszero(x)
_isthinzero(d::Dual) = _isthinzero(value(d)) && all(_isthinzero, partials(d))

function Base.:(^)(x::Dual{Txy,<:Interval}, y::Dual{Txy,<:Interval}) where {Txy}
vx, vy = value(x), value(y)
expv = vx^vy
Expand All @@ -49,16 +63,31 @@ function Base.:(^)(x::Dual{Tx,<:Interval}, y::Dual{Ty,<:Interval}) where {Tx,Ty}
end
end

function Base.:(^)(x::Dual{Tx,<:Interval}, y::Interval) where {Tx}
function Base.:(^)(x::Dual{Tx,<:NestedInterval}, y::Interval) where {Tx}
v = value(x)
expv = v^y
if isthinzero(y) || all(isthinzero, values(partials(x)))
if isthinzero(y) || all(_isthinzero, values(partials(x)))
return Dual{Tx}(expv, zero(partials(x)))
else
return Dual{Tx}(expv, partials(x) * y * v^(y - interval(1)))
end
end

# A `y::Real` method would be ambiguous with ForwardDiff's generated `^`
# methods and its `Dual^Dual` tag methods, so define the relevant concrete
# exponent types instead.
for R in (:Integer, :Rational, :AbstractFloat, :Irrational)
@eval function Base.:(^)(x::Dual{Tx,<:NestedInterval}, y::$R) where {Tx}
v = value(x)
expv = v^y
if iszero(y) || all(_isthinzero, values(partials(x)))
return Dual{Tx}(expv, zero(partials(x)))
else
return Dual{Tx}(expv, partials(x) * y * v^(y - 1))
end
end
end

function Base.:(^)(x::Interval, y::Dual{Ty,<:Interval}) where {Ty}
v = value(y)
expv = x^v
Expand All @@ -76,7 +105,7 @@ Base.:(^)(x::ExactReal, y::Dual{<:Any,I}) where {I<:Interval} = convert(I, x)^y
function Base.:(^)(x::Dual{Tx}, y::ExactReal) where {Tx}
v = value(x)
expv = v^y
if iszero(y.value) || all(iszero, values(partials(x)))
if iszero(y.value) || all(_isthinzero, values(partials(x)))
return Dual{Tx}(expv, zero(partials(x)))
else
return Dual{Tx}(expv, partials(x) * y * v^(y - 1))
Expand Down
31 changes: 31 additions & 0 deletions test/interval_tests/forwarddiff.jl
Original file line number Diff line number Diff line change
Expand Up @@ -105,4 +105,35 @@ end
@exact g(x) = 2^x + 6sin(x^3) - 33
@test isguaranteed(ForwardDiff.derivative(f, interval(1)))
end

@testset "thick partials" begin
# Differentiating by a real leaves interval-valued constants in the partials.
x, w = 2.0, interval(-0.5, 0.5)

# Test real, interval, and exact exponents.
for n ∈ (4, 4.0, interval(4), exact(4))
@test isequal_interval(ForwardDiff.derivative(t -> (x + t*w)^n, 0), interval(4x^3) * w)
@test isequal_interval(ForwardDiff.derivative(t -> ForwardDiff.derivative(s -> (x + s*w)^n, t), 0),
interval(12x^2) * w * w)
end

# Each derivative order adds a `Dual` layer.
ϕ(t) = (x + t*w)^4
dϕ(t) = ForwardDiff.derivative(ϕ, t)
ddϕ(t) = ForwardDiff.derivative(dϕ, t)
dddϕ(t) = ForwardDiff.derivative(ddϕ, t)
ddddϕ(t) = ForwardDiff.derivative(dddϕ, t)

@test isequal_interval(dϕ(0) , interval(4x^3) * w)
@test isequal_interval(ddϕ(0) , interval(12x^2) * w * w)
@test isequal_interval(dddϕ(0) , interval(24x) * w * w * w)
@test isequal_interval(ddddϕ(0), interval(24) * w * w * w * w)

# Test multivariate first- and second-order derivatives.
ψ(v) = (v[1]*w + v[2])^5
@test all(isequal_interval.(ForwardDiff.gradient(ψ, [0, 1]), [interval(5) * w, interval(5)]))
@test all(isequal_interval.(ForwardDiff.hessian(ψ, [0, 1]),
[interval(20) * w * w interval(20) * w
interval(20) * w interval(20) ]))
end
end
Loading