Implement piecewise function

docs/src/manual/usage.md

@@ -123,6 +123,23 @@ One can refer to the following:
 - `setdiff`: cannot be used with intervals. See instead [`interiordiff`](@ref).
 
 
+## Piecewise functions
+
+Since intervals don't play well with `if .. else .. end` statement,
+we provide a utility to define function by pieces:
+
+```@repl usage
+myabs = Piecewise(
+    Domain{:open, :closed}(-Inf, 0) => x -> -x,
+    Domain{:open, :open}(0, Inf) => identity
+)
+myabs(-1.23)
+myabs(interval(-1, 23))
+```
+
+The resulting function work with both standard numbers and intervals,
+and deal properly with the decorations of intervals.
+
 
 ## Custom interval bounds type
 

ext/IntervalArithmeticForwardDiffExt.jl

@@ -1,6 +1,7 @@
 module IntervalArithmeticForwardDiffExt
 
 using IntervalArithmetic, ForwardDiff
+using IntervalArithmetic: isempty_domain, overlap_domain, intersect_domain, in_domain, leftof
 using ForwardDiff: Dual, Partials, ≺, value, partials
 
 ForwardDiff.can_dual(::Type{ExactReal}) = true
@@ -90,4 +91,49 @@ function Base.:(^)(x::ExactReal, y::Dual{<:Ty}) where {Ty}
     end
 end
 
+
+# Piecewise functions
+
+function (constant::Constant)(::Dual{T, Interval{S}}) where {T, S}
+    return Dual{T}(interval(S, constant.value), interval(S, 0.0))
+end
+
+function (piecewise::Piecewise)(dual::Dual{T, <:Interval}) where {T}
+    X = value(dual)
+    input_domain = Domain(X)
+    if !overlap_domain(input_domain, piecewise) 
+        return Dual{T}(emptyinterval(X), emptyinterval(X) .* partials(dual))
+    end
+
+    if !in_domain(input_domain, piecewise)
+        dec = trv
+    elseif any(x -> in_domain(x, input_domain), discontinuities(piecewise, 1))
+        dec = def 
+    else
+        dec = com
+    end
+
+    dual_piece_outputs = []
+    for (piece_domain, f) in pieces(piecewise)
+        piece_input = intersect_domain(input_domain, piece_domain)
+        isempty_domain(piece_input) && continue
+        sub_X = interval(inf(piece_input), sup(piece_input), decoration(X))
+        sub_dual = Dual{T}(sub_X, partials(dual))
+        push!(dual_piece_outputs, f(sub_dual))
+    end
+
+    piece_outputs = value.(dual_piece_outputs)
+    dec = min(dec, minimum(decoration.(piece_outputs)))
+    primal = IntervalArithmetic.setdecoration(reduce(hull, piece_outputs), dec)
+
+    doutputs = partials.(dual_piece_outputs)
+    partial = map(zip(doutputs...)) do pp
+        pdec = min(dec, minimum(decoration.(pp)))
+        return IntervalArithmetic.setdecoration(reduce(hull, pp), pdec)
+    end
+
+    return Dual{T}(primal, tuple(partial...))
 end
+
+
+end

src/IntervalArithmetic.jl

@@ -23,6 +23,11 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 
 #
 
+include("piecewise.jl")
+    export Domain, Constant, Piecewise, domains, discontinuities, pieces
+
+#
+
 include("display.jl")
     export setdisplay
 

src/intervals/interval_operations/set_operations.jl

@@ -201,3 +201,17 @@ function _interiordiff(x::Interval, y::Interval)
     t = isguaranteed(x) & isguaranteed(y)
     return (_unsafe_interval(h₁, d, t), _unsafe_interval(h₂, d, t), _unsafe_interval(inter, d, t))
 end
+
+function interval_diff(x::Interval{T}, y::Interval) where {T<:NumTypes}
+    isdisjoint_interval(x, y) && return [x]
+    issubset_interval(x, y) && return Interval{T}[]
+
+    intersection = intersect_interval(x, y)
+    inf(x) == inf(intersection) && return [interval(sup(intersection), sup(x))]
+    sup(x) == sup(intersection) && return [interval(inf(x), inf(intersection))]
+
+    return [
+        interval(inf(x), inf(intersection)),
+        interval(sup(intersection), sup(x))
+    ]
+end