diff --git a/ext/IntervalArithmeticForwardDiffExt.jl b/ext/IntervalArithmeticForwardDiffExt.jl index a789cc21..05529beb 100644 --- a/ext/IntervalArithmeticForwardDiffExt.jl +++ b/ext/IntervalArithmeticForwardDiffExt.jl @@ -28,6 +28,32 @@ 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 tests partials for structural zero with `iszero`, whose generic +# fallback is `x == zero(x)` and therefore throws `InconclusiveBooleanOperation` +# for a non-thin interval. A partial carries no perturbation exactly when it is +# (recursively) the thin interval [0, 0], so `isthinzero` is a safe and correct +# answer. +# +# `_iszero` is kept separate from `Base.iszero(::Interval)`, so this change is +# restricted to this extension. +# +# Higher-order derivatives nest `Dual`s around intervals to a depth equal to the +# derivative order, and dispatch cannot express "an `Interval` under any number +# of `Dual` layers"; `NestedInterval` spells out up to 4 layers, which seems +# sufficient for practical use. +const NestedInterval = let + U = Interval + for _ in 1:4 + U = Union{U, Dual{T,<:U} where {T}} + end + U +end +_iszero(x::Interval) = isthinzero(x) +_iszero(d::Dual) = _iszero(value(d)) & _iszero(partials(d)) +_iszero(p::Partials) = all(_iszero, p.values) +Base.iszero(d::Dual{T,<:NestedInterval}) where {T} = _iszero(d) +Base.iszero(p::Partials{N,<:NestedInterval}) where {N} = _iszero(p) + function Base.:(^)(x::Dual{Txy,<:Interval}, y::Dual{Txy,<:Interval}) where {Txy} vx, vy = value(x), value(y) expv = vx^vy diff --git a/test/interval_tests/forwarddiff.jl b/test/interval_tests/forwarddiff.jl index 906a973e..55cf7bea 100644 --- a/test/interval_tests/forwarddiff.jl +++ b/test/interval_tests/forwarddiff.jl @@ -105,4 +105,38 @@ end @exact g(x) = 2^x + 6sin(x^3) - 33 @test isguaranteed(ForwardDiff.derivative(f, interval(1))) end + + @testset "thick partials" begin + # differentiating with respect to a real leaves interval-valued constants + # in the partials, where a partial is a structural zero only when thin + @test !iszero(ForwardDiff.Partials((interval(-0.5, 0.5),))) + @test iszero(ForwardDiff.Partials((interval(0),))) + + x, w = 2.0, interval(-0.5, 0.5) + + # the exponent reaches `^` as a `Real`, an `Interval` or an `ExactReal` + for n ∈ (4, 4.0, interval(4), exact(4)) + @test isequal_interval(ForwardDiff.derivative(t -> (x + t*w)^n, 0), interval(4x^3) * w) + end + + # each derivative order wraps another `Dual` around the interval + ϕ(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) + + # ψ(v) = (v₁w + v₂)^5 has ∂ψ/∂v₁ = 5(v₁w + v₂)^4 w and ∂²ψ/∂v₁² = 20(v₁w + v₂)^3 w², + # evaluated here at v = [0, 1], where v₁w + v₂ is thin + ψ(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