From dc7bdf0bb8ffa6967a30d4b9a9369caf04913c76 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:16 +0100
Subject: [PATCH 01/20] Prototype

---
 src/piecewise.jl | 54 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)
 create mode 100644 src/piecewise.jl

diff --git a/src/piecewise.jl b/src/piecewise.jl
new file mode 100644
index 00000000..04a02822
--- /dev/null
+++ b/src/piecewise.jl
@@ -0,0 +1,54 @@
+using IntervalArithmetic
+
+# represents a piecewise 1D function
+struct Piecewise{T <: Tuple}
+    pieces::T 
+    continuous::Bool
+end
+
+Piecewise(pairs... ; continuous = false) = Piecewise(pairs, continuous)
+
+domain(piecewise::Piecewise) = reduce(hull, subdomains(piecewise))
+subdomains(piecewise::Piecewise) = first.(piecewise.pieces)
+
+intersecting(X::Interval, Y::Interval) = !isempty_interval(intersect_interval(X, Y))
+
+function (piecewise::Piecewise)(X::Interval)
+    intersections = intersecting.(X, subdomains(piecewise))
+    dom = domain(piecewise)
+    if !issubset_interval(X, dom)
+        dec = trv
+    elseif count(intersections) > 1
+        if piecewise.continuous
+            dec = dac
+        else
+            dec = def
+        end
+    else
+        dec = com
+    end
+
+    used_pieces = filter(piece -> intersecting(X, piece[1]), piecewise.pieces)
+    outputs = map(used_pieces) do (region, f)
+        S = IntervalArithmetic.setdecoration(intersect_interval(X, region), decoration(X))
+        return f(S)
+    end
+
+    dec = min(dec, minimum(decoration.(outputs)))
+    return IntervalArithmetic.setdecoration(reduce(hull, outputs), dec)
+end
+
+function (piecewise::Piecewise)(x::Real)
+    for (region, f) in piecewise.pieces
+        in_interval(x, region) && return f(x)
+    end
+    throw(DomainError("piecewise function was called with $x which is outside of its domain $(domain(piecewise))"))
+end
+
+f = Piecewise(
+        0..3 => x -> x + 1,
+        3..6 => identity,
+        6..Inf => x -> x + 2
+)
+
+f(2..5)
\ No newline at end of file

From b503f57c89cd287060a9dd1dfb2206f970f71c7a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:16 +0100
Subject: [PATCH 02/20] Add basic tests and Constant

---
 src/IntervalArithmetic.jl        |  5 +++++
 src/piecewise.jl                 | 21 ++++++++++-----------
 test/interval_tests/piecewise.jl | 32 ++++++++++++++++++++++++++++++++
 3 files changed, 47 insertions(+), 11 deletions(-)
 create mode 100644 test/interval_tests/piecewise.jl

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 19e1b445..a1f13138 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -23,6 +23,11 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 
 #
 
+include("piecewise.jl")
+    export Constant, Piecewise
+
+#
+
 include("display.jl")
     export setdisplay
 
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 04a02822..edc144c1 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -1,6 +1,10 @@
-using IntervalArithmetic
+struct Constant{T}
+    value::T
+end
+
+(constant::Constant)(::Any) = constant.value
+(constant::Constant)(::Interval) = interval(constant.value)
 
-# represents a piecewise 1D function
 struct Piecewise{T <: Tuple}
     pieces::T 
     continuous::Bool
@@ -31,9 +35,12 @@ function (piecewise::Piecewise)(X::Interval)
     used_pieces = filter(piece -> intersecting(X, piece[1]), piecewise.pieces)
     outputs = map(used_pieces) do (region, f)
         S = IntervalArithmetic.setdecoration(intersect_interval(X, region), decoration(X))
+        @show S
         return f(S)
     end
 
+    @show outputs
+
     dec = min(dec, minimum(decoration.(outputs)))
     return IntervalArithmetic.setdecoration(reduce(hull, outputs), dec)
 end
@@ -43,12 +50,4 @@ function (piecewise::Piecewise)(x::Real)
         in_interval(x, region) && return f(x)
     end
     throw(DomainError("piecewise function was called with $x which is outside of its domain $(domain(piecewise))"))
-end
-
-f = Piecewise(
-        0..3 => x -> x + 1,
-        3..6 => identity,
-        6..Inf => x -> x + 2
-)
-
-f(2..5)
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
new file mode 100644
index 00000000..976db0af
--- /dev/null
+++ b/test/interval_tests/piecewise.jl
@@ -0,0 +1,32 @@
+@testset "Step function" begin
+    step = Piecewise(
+        interval(-Inf, 0) => Constant(0),
+        interval(0, Inf) => Constant(1)
+    )
+
+    @test step(-1) == 0
+    @test step(100) == 1
+    @test isequal_interval(step(interval(-3.2, -2.1)), interval(0))
+    @test decoration(step(interval(-3.33))) == com
+    @test isequal_interval(step(interval(2.3, 3.4)), interval(1))
+    @test decoration(step(interval(4.44))) == com
+    @test isequal_interval(step(interval(-22.2, 33.3)), interval(0, 1))
+    @test decoration(step(interval(-11, 11))) == def
+end
+
+@testset "abs" begin
+    myabs = Piecewise(
+        interval(-Inf, 0) => x -> -x,
+        interval(0, Inf) => identity ;
+        continuous = true
+    )
+
+    @test myabs(-1) == 1
+    @test myabs(100) == 100
+    @test isequal_interval(myabs(interval(-3.2, -2.1)), interval(2.1, 3.2))
+    @test decoration(myabs(interval(-3.33))) == com
+    @test isequal_interval(myabs(interval(2.3, 3.4)), interval(2.3, 3.4))
+    @test decoration(myabs(interval(4.444))) == com
+    @test isequal_interval(myabs(interval(-22.2, 33.3)), interval(0, 33.3))
+    @test decoration(myabs(interval(-11, 11))) == dac
+end
\ No newline at end of file

From 67436967c71a372b4550f6f9a5b468d4126e4f19 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:16 +0100
Subject: [PATCH 03/20] Add tests for and correct out-of-domain behavior

---
 src/piecewise.jl                 |  4 +---
 test/interval_tests/piecewise.jl | 13 +++++++++++++
 2 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/piecewise.jl b/src/piecewise.jl
index edc144c1..8e839836 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -20,6 +20,7 @@ intersecting(X::Interval, Y::Interval) = !isempty_interval(intersect_interval(X,
 function (piecewise::Piecewise)(X::Interval)
     intersections = intersecting.(X, subdomains(piecewise))
     dom = domain(piecewise)
+    !intersecting(X, dom) && return emptyinterval()
     if !issubset_interval(X, dom)
         dec = trv
     elseif count(intersections) > 1
@@ -35,12 +36,9 @@ function (piecewise::Piecewise)(X::Interval)
     used_pieces = filter(piece -> intersecting(X, piece[1]), piecewise.pieces)
     outputs = map(used_pieces) do (region, f)
         S = IntervalArithmetic.setdecoration(intersect_interval(X, region), decoration(X))
-        @show S
         return f(S)
     end
 
-    @show outputs
-
     dec = min(dec, minimum(decoration.(outputs)))
     return IntervalArithmetic.setdecoration(reduce(hull, outputs), dec)
 end
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index 976db0af..5552d7f2 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -29,4 +29,17 @@ end
     @test decoration(myabs(interval(4.444))) == com
     @test isequal_interval(myabs(interval(-22.2, 33.3)), interval(0, 33.3))
     @test decoration(myabs(interval(-11, 11))) == dac
+end
+
+@testset "Out of domain" begin
+    window = Piecewise(
+        interval(-π, π) => x -> 1/2 * (cos(x) + 1)
+    )
+
+    @test_throws DomainError window(123)
+    @test isequal_interval(window(interval(0, π)), interval(0, 1))
+    @test decoration(window(interval(-π, 0))) == com
+    @test isequal_interval(window(interval(-10, 10)), interval(0, 1))
+    @test decoration(window(interval(-10, 10))) == trv
+    @test isempty_interval(window(interval(100, 1000)))
 end
\ No newline at end of file

From b934dad93743e97d5379b3d119fdca4c138de587 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:16 +0100
Subject: [PATCH 04/20] New strategy defining pieces

---
 src/IntervalArithmetic.jl                     |  2 +-
 .../interval_operations/set_operations.jl     | 14 ++++
 src/piecewise.jl                              | 73 +++++++++++++------
 3 files changed, 66 insertions(+), 23 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index a1f13138..69a4e1eb 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -24,7 +24,7 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 #
 
 include("piecewise.jl")
-    export Constant, Piecewise
+    export Constant, Piece, Piecewise, intersecting, domain
 
 #
 
diff --git a/src/intervals/interval_operations/set_operations.jl b/src/intervals/interval_operations/set_operations.jl
index b95750f6..792d36c5 100644
--- a/src/intervals/interval_operations/set_operations.jl
+++ b/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]
+    isequal_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
\ No newline at end of file
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 8e839836..1acb4129 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -5,38 +5,67 @@ end
 (constant::Constant)(::Any) = constant.value
 (constant::Constant)(::Interval) = interval(constant.value)
 
-struct Piecewise{T <: Tuple}
-    pieces::T 
-    continuous::Bool
+struct Piece{T}
+    f
+    domain::Interval{T}
+    left_continuity::Int
 end
 
-Piecewise(pairs... ; continuous = false) = Piecewise(pairs, continuous)
+Piece(f, domain::Interval) = Piece(f, domain, -1)
+Piece(x::Real, domain::Interval) = Piece(Constant(x), domain)
 
-domain(piecewise::Piecewise) = reduce(hull, subdomains(piecewise))
-subdomains(piecewise::Piecewise) = first.(piecewise.pieces)
+struct Piecewise{T}
+    pieces::Vector{Piece{T}}
+    holes::Vector{Interval{T}}
+end
+
+function Piecewise(pieces...)
+    pieces = sort(collect(pieces) ; lt = (p1, p2) -> precedes(p1.domain, p2.domain))
+    subdomains = domain.(pieces)
+
+    for (S1, S2) in zip(subdomains[1:end-1], subdomains[2:end])
+        sup(S1) > inf(S2) && throw(ArgumentError("pieces are not disjoint"))
+    end
+
+    holes = [entireinterval()]
+
+    for S in subdomains
+        holes = mapreduce(vcat, holes) do hole
+            return interval_diff(hole, S)
+        end
+    end
+
+    return Piecewise(pieces, holes)
+end
+
+domain(piece::Piece) = piece.domain
+domain(piecewise::Piecewise) = domain.(piecewise.pieces)
 
 intersecting(X::Interval, Y::Interval) = !isempty_interval(intersect_interval(X, Y))
 
+function interior_intersecting(X::Interval, Y::Interval)
+    inter = intersect_interval(X, Y)
+    (inter == inf(X) || inter == sup(X)) && return false
+    return !isempty_interval(inter)
+end
+
+# TODO This is not true for interval bordering a hole.
 function (piecewise::Piecewise)(X::Interval)
-    intersections = intersecting.(X, subdomains(piecewise))
-    dom = domain(piecewise)
-    !intersecting(X, dom) && return emptyinterval()
-    if !issubset_interval(X, dom)
+    if any(interior_intersecting.(X, piecewise.holes))
         dec = trv
-    elseif count(intersections) > 1
-        if piecewise.continuous
-            dec = dac
-        else
-            dec = def
-        end
     else
         dec = com
+        for piece in piecewise.pieces
+            if in_interval(inf(domain(piece)), X) && piece.left_continuity < 0
+                dec = def
+            end
+        end
     end
 
-    used_pieces = filter(piece -> intersecting(X, piece[1]), piecewise.pieces)
-    outputs = map(used_pieces) do (region, f)
-        S = IntervalArithmetic.setdecoration(intersect_interval(X, region), decoration(X))
-        return f(S)
+    used_pieces = filter(piece -> intersecting(X, domain(piece)), piecewise.pieces)
+    outputs = map(used_pieces) do piece
+        S = IntervalArithmetic.setdecoration(intersect_interval(X, domain(piece)), decoration(X))
+        return piece.f(S)
     end
 
     dec = min(dec, minimum(decoration.(outputs)))
@@ -44,8 +73,8 @@ function (piecewise::Piecewise)(X::Interval)
 end
 
 function (piecewise::Piecewise)(x::Real)
-    for (region, f) in piecewise.pieces
-        in_interval(x, region) && return f(x)
+    for piece in piecewise.pieces
+        in_interval(x, domain(piece)) && return piece.f(x)
     end
     throw(DomainError("piecewise function was called with $x which is outside of its domain $(domain(piecewise))"))
 end
\ No newline at end of file

From 6966ddd38ed3178f52223e8e3b7797554c34f17c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:16 +0100
Subject: [PATCH 05/20] Use Intervals for the domain of the pieces

---
 Project.toml              |   1 +
 src/IntervalArithmetic.jl |   2 +-
 src/piecewise.jl          | 115 +++++++++++++++++++++++---------------
 3 files changed, 71 insertions(+), 47 deletions(-)

diff --git a/Project.toml b/Project.toml
index 75dd7071..0fcf1bc0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.22.22"
 
 [deps]
 CRlibm_jll = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
+Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 69a4e1eb..84668ff1 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -24,7 +24,7 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 #
 
 include("piecewise.jl")
-    export Constant, Piece, Piecewise, intersecting, domain
+    export Constant, Piecewise, Domain, domain, subdomains, discontinuities
 
 #
 
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 1acb4129..9ca3b9f2 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -1,3 +1,11 @@
+import Intervals
+import Intervals: IntervalSet, less_than_disjoint
+
+const Domain = Intervals.Interval
+
+inf(x::Domain) = first(x)
+sup(x::Domain) = last(x)
+
 struct Constant{T}
     value::T
 end
@@ -5,76 +13,91 @@ end
 (constant::Constant)(::Any) = constant.value
 (constant::Constant)(::Interval) = interval(constant.value)
 
-struct Piece{T}
-    f
-    domain::Interval{T}
-    left_continuity::Int
+struct Piecewise
+    domain::IntervalSet
+    fs
+    discontinuities::Vector
 end
 
-Piece(f, domain::Interval) = Piece(f, domain, -1)
-Piece(x::Real, domain::Interval) = Piece(Constant(x), domain)
+function Piecewise(
+        subdomains::Vector{<:Domain},
+        fs,
+        continuous::Vector{Bool} = fill(false, length(subdomains) - 1))
 
-struct Piecewise{T}
-    pieces::Vector{Piece{T}}
-    holes::Vector{Interval{T}}
-end
+    if length(subdomains) != length(fs)
+        throw(ArgumentError("the number of domains and the number of functions don't match"))
+    end
+    
+    if length(subdomains) - 1 != length(continuous)
+        n = length(subdomains)
+        throw(ArgumentError("$(length(sub)) junction points but $(n - 1) are expected based on the number of domains ($n)"))
+    end
 
-function Piecewise(pieces...)
-    pieces = sort(collect(pieces) ; lt = (p1, p2) -> precedes(p1.domain, p2.domain))
-    subdomains = domain.(pieces)
+    for k in 1:length(subdomains) - 1
+        s1 = subdomains[k]
+        s2 = subdomains[k + 1]
 
-    for (S1, S2) in zip(subdomains[1:end-1], subdomains[2:end])
-        sup(S1) > inf(S2) && throw(ArgumentError("pieces are not disjoint"))
+        !less_than_disjoint(s1, s2) && throw(ArgumentError("domains are either not ordered or not disjoint"))
     end
 
-    holes = [entireinterval()]
+    singularities = sup.(subdomains[1:end-1])
 
-    for S in subdomains
-        holes = mapreduce(vcat, holes) do hole
-            return interval_diff(hole, S)
-        end
-    end
+    return Piecewise(IntervalSet(subdomains), fs, singularities[.!continuous])
+end
 
-    return Piecewise(pieces, holes)
+function Piecewise(pairs::Vararg{<:Pair} ; continuous = fill(false, length(pairs) - 1))
+    pairs = collect(pairs)
+    return Piecewise(first.(pairs), last.(pairs), continuous)
 end
 
-domain(piece::Piece) = piece.domain
-domain(piecewise::Piecewise) = domain.(piecewise.pieces)
+subdomains(piecewise::Piecewise) = convert(Vector, piecewise.domain)
+domain(piecewise::Piecewise) = piecewise.domain
+pieces(piecewise::Piecewise) = zip(subdomains(piecewise), piecewise.fs)
+discontinuities(piecewise::Piecewise) = piecewise.discontinuities
+
+domain_string(x::Domain) = repr(x ; context = IOContext(stdout, :compact => true))
 
-intersecting(X::Interval, Y::Interval) = !isempty_interval(intersect_interval(X, Y))
+function domain_string(S::IntervalSet)
+    r = repr("text/plain", S ; context = IOContext(stdout, :compact => true))
+    return join(split(r, "\n")[2:end], " ∪ ")
+end
+
+domain_string(piecewise::Piecewise) = domain_string(domain(piecewise))
 
-function interior_intersecting(X::Interval, Y::Interval)
-    inter = intersect_interval(X, Y)
-    (inter == inf(X) || inter == sup(X)) && return false
-    return !isempty_interval(inter)
+function Base.show(io::IO, ::MIME"text/plain", piecewise::Piecewise)
+    n = length(pieces(piecewise))
+    print(io, "Piecewise function with $n pieces:")
+    
+    for (subdomain, f) in pieces(piecewise)
+        println(io)
+        print(io, "  $(domain_string(subdomain)) -> $(repr(f))")
+    end
 end
 
-# TODO This is not true for interval bordering a hole.
-function (piecewise::Piecewise)(X::Interval)
-    if any(interior_intersecting.(X, piecewise.holes))
+function (piecewise::Piecewise)(X::Interval{T}) where T
+    set = Domain(inf(X), sup(X), true, true)
+    if !isempty(setdiff(set, domain(piecewise)))
         dec = trv
+    elseif any(in(set), discontinuities(piecewise))
+        dec = def 
     else
         dec = com
-        for piece in piecewise.pieces
-            if in_interval(inf(domain(piece)), X) && piece.left_continuity < 0
-                dec = def
-            end
-        end
     end
 
-    used_pieces = filter(piece -> intersecting(X, domain(piece)), piecewise.pieces)
-    outputs = map(used_pieces) do piece
-        S = IntervalArithmetic.setdecoration(intersect_interval(X, domain(piece)), decoration(X))
-        return piece.f(S)
+    results = Interval{T}[]
+    for (subdomain, f) in pieces(piecewise)
+        subset = intersect(set, subdomain)
+        isempty(subset) && continue
+        push!(results, f(interval(inf(subset), sup(subset), decoration(X))))
     end
 
-    dec = min(dec, minimum(decoration.(outputs)))
-    return IntervalArithmetic.setdecoration(reduce(hull, outputs), dec)
+    dec = min(dec, minimum(decoration.(results)))
+    return IntervalArithmetic.setdecoration(reduce(hull, results), dec)
 end
 
 function (piecewise::Piecewise)(x::Real)
-    for piece in piecewise.pieces
-        in_interval(x, domain(piece)) && return piece.f(x)
+    for (subdomain, f) in pieces(piecewise)
+        (x in subdomain) && return f(x)
     end
-    throw(DomainError("piecewise function was called with $x which is outside of its domain $(domain(piecewise))"))
+    throw(DomainError(x, "piecewise function was called outside of its domain $(domain_string(piecewise))"))
 end
\ No newline at end of file

From 50574153296e8531a1a62e766d68adbc5385a64d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:55 +0100
Subject: [PATCH 06/20] Compat with forwarddiff

---
 ext/IntervalArithmeticForwardDiffExt.jl | 58 +++++++++++++++++++++++++
 src/IntervalArithmetic.jl               |  2 +-
 src/piecewise.jl                        | 17 ++++----
 3 files changed, 68 insertions(+), 9 deletions(-)

diff --git a/ext/IntervalArithmeticForwardDiffExt.jl b/ext/IntervalArithmeticForwardDiffExt.jl
index 80f1b818..c64cbd9d 100644
--- a/ext/IntervalArithmeticForwardDiffExt.jl
+++ b/ext/IntervalArithmeticForwardDiffExt.jl
@@ -90,4 +90,62 @@ function Base.:(^)(x::ExactReal, y::Dual{<:Ty}) where {Ty}
     end
 end
 
+# resolve ambiguity
+
+Base.convert(::Type{Dual{T,V,N}}, x::ExactReal) where {T,V,N} = Dual{T}(V(x), zero(Partials{N,V}))
+
+
+# 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 ForwardDiff.derivative(piecewise::Piecewise)
+    dfs = map(piecewise.fs)  do f
+        x -> ForwardDiff.derivative(f, x)
+    end
+
+    return Piecewise(
+        piecewise.domain,
+        dfs,
+        continuity .- 1,
+        piecewise.singularities
+    )
+end
+
+function (piecewise::Piecewise)(dual::Dual{T, <:Interval}) where {T}
+    X = value(dual)
+    set = Domain(inf(X), sup(X), true, true)
+    if !isempty(setdiff(set, domain(piecewise)))
+        dec = trv
+    elseif any(in(set), discontinuities(piecewise))
+        dec = def 
+    else
+        dec = com
+    end
+
+    dual_results = []
+    for (subdomain, f) in pieces(piecewise)
+        subset = intersect(set, subdomain)
+        isempty(subset) && continue
+        sub_X = interval(inf(subset), sup(subset), decoration(X))
+        sub_dual = Dual{T}(sub_X, partials(dual))
+        push!(dual_results, f(sub_dual))
+    end
+
+    results = value.(dual_results)
+    dec = min(dec, minimum(decoration.(results)))
+    primal = IntervalArithmetic.setdecoration(reduce(hull, results), dec)
+
+    dresults = partials.(dual_results)
+    partial = map(zip(dresults...)) do pp
+        pdec = min(dec, minimum(decoration.(pp)))
+        return IntervalArithmetic.setdecoration(reduce(hull, pp), pdec)
+    end
+
+    return Dual{T}(primal, tuple(partial...))
+end
+
+
+end
\ No newline at end of file
diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 84668ff1..926f72da 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -24,7 +24,7 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 #
 
 include("piecewise.jl")
-    export Constant, Piecewise, Domain, domain, subdomains, discontinuities
+    export Constant, Piecewise, Domain, domain, subdomains, discontinuities, pieces
 
 #
 
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 9ca3b9f2..47297f47 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -16,19 +16,20 @@ end
 struct Piecewise
     domain::IntervalSet
     fs
-    discontinuities::Vector
+    continuity::Vector{Int}
+    singularities::Vector
 end
 
 function Piecewise(
         subdomains::Vector{<:Domain},
         fs,
-        continuous::Vector{Bool} = fill(false, length(subdomains) - 1))
+        continuity::Vector{Int} = fill(-1, length(subdomains) - 1))
 
     if length(subdomains) != length(fs)
         throw(ArgumentError("the number of domains and the number of functions don't match"))
     end
     
-    if length(subdomains) - 1 != length(continuous)
+    if length(subdomains) - 1 != length(continuity)
         n = length(subdomains)
         throw(ArgumentError("$(length(sub)) junction points but $(n - 1) are expected based on the number of domains ($n)"))
     end
@@ -42,18 +43,18 @@ function Piecewise(
 
     singularities = sup.(subdomains[1:end-1])
 
-    return Piecewise(IntervalSet(subdomains), fs, singularities[.!continuous])
+    return Piecewise(IntervalSet(subdomains), fs, continuity, singularities)
 end
 
-function Piecewise(pairs::Vararg{<:Pair} ; continuous = fill(false, length(pairs) - 1))
+function Piecewise(pairs::Vararg{<:Pair} ; continuity = fill(-1, length(pairs) - 1))
     pairs = collect(pairs)
-    return Piecewise(first.(pairs), last.(pairs), continuous)
+    return Piecewise(first.(pairs), last.(pairs), continuity)
 end
 
 subdomains(piecewise::Piecewise) = convert(Vector, piecewise.domain)
 domain(piecewise::Piecewise) = piecewise.domain
 pieces(piecewise::Piecewise) = zip(subdomains(piecewise), piecewise.fs)
-discontinuities(piecewise::Piecewise) = piecewise.discontinuities
+discontinuities(piecewise::Piecewise) = piecewise.singularities[piecewise.continuity .< 0]
 
 domain_string(x::Domain) = repr(x ; context = IOContext(stdout, :compact => true))
 
@@ -74,7 +75,7 @@ function Base.show(io::IO, ::MIME"text/plain", piecewise::Piecewise)
     end
 end
 
-function (piecewise::Piecewise)(X::Interval{T}) where T
+function (piecewise::Piecewise)(X::Interval{T}) where {T}
     set = Domain(inf(X), sup(X), true, true)
     if !isempty(setdiff(set, domain(piecewise)))
         dec = trv

From 72f08f113e5501bddb06cee6e66369fbbc50c15c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:55 +0100
Subject: [PATCH 07/20] Add tests

---
 ext/IntervalArithmeticForwardDiffExt.jl | 17 ++----
 src/piecewise.jl                        |  2 +
 test/Project.toml                       |  1 +
 test/interval_tests/piecewise.jl        | 75 ++++++++++++++++++++++---
 4 files changed, 74 insertions(+), 21 deletions(-)

diff --git a/ext/IntervalArithmeticForwardDiffExt.jl b/ext/IntervalArithmeticForwardDiffExt.jl
index c64cbd9d..5017912f 100644
--- a/ext/IntervalArithmeticForwardDiffExt.jl
+++ b/ext/IntervalArithmeticForwardDiffExt.jl
@@ -101,22 +101,13 @@ function (constant::Constant)(::Dual{T, Interval{S}}) where {T, S}
     return Dual{T}(interval(S, constant.value), interval(S, 0.0))
 end
 
-function ForwardDiff.derivative(piecewise::Piecewise)
-    dfs = map(piecewise.fs)  do f
-        x -> ForwardDiff.derivative(f, x)
-    end
-
-    return Piecewise(
-        piecewise.domain,
-        dfs,
-        continuity .- 1,
-        piecewise.singularities
-    )
-end
-
 function (piecewise::Piecewise)(dual::Dual{T, <:Interval}) where {T}
     X = value(dual)
     set = Domain(inf(X), sup(X), true, true)
+    if isdisjoint(set, domain(piecewise))
+        return Dual{T}(emptyinterval(X), emptyinterval(X) .* partials(dual))
+    end
+
     if !isempty(setdiff(set, domain(piecewise)))
         dec = trv
     elseif any(in(set), discontinuities(piecewise))
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 47297f47..d068c8b2 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -77,6 +77,8 @@ end
 
 function (piecewise::Piecewise)(X::Interval{T}) where {T}
     set = Domain(inf(X), sup(X), true, true)
+    isdisjoint(set, domain(piecewise)) && return emptyinterval(T)
+
     if !isempty(setdiff(set, domain(piecewise)))
         dec = trv
     elseif any(in(set), discontinuities(piecewise))
diff --git a/test/Project.toml b/test/Project.toml
index 2dd97b19..9b7863b7 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,6 +2,7 @@
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index 5552d7f2..adfa3d4c 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -1,7 +1,7 @@
 @testset "Step function" begin
     step = Piecewise(
-        interval(-Inf, 0) => Constant(0),
-        interval(0, Inf) => Constant(1)
+        Domain(-Inf, 0, false, true) => Constant(0),
+        Domain(0, Inf, false, false) => Constant(1)
     )
 
     @test step(-1) == 0
@@ -16,9 +16,9 @@ end
 
 @testset "abs" begin
     myabs = Piecewise(
-        interval(-Inf, 0) => x -> -x,
-        interval(0, Inf) => identity ;
-        continuous = true
+        Domain(-Inf, 0, false, true) => x -> -x,
+        Domain(0, Inf, false, false) => identity ;
+        continuity = [0]
     )
 
     @test myabs(-1) == 1
@@ -28,12 +28,12 @@ end
     @test isequal_interval(myabs(interval(2.3, 3.4)), interval(2.3, 3.4))
     @test decoration(myabs(interval(4.444))) == com
     @test isequal_interval(myabs(interval(-22.2, 33.3)), interval(0, 33.3))
-    @test decoration(myabs(interval(-11, 11))) == dac
+    @test decoration(myabs(interval(-11, 11))) == com
 end
 
 @testset "Out of domain" begin
     window = Piecewise(
-        interval(-π, π) => x -> 1/2 * (cos(x) + 1)
+        Domain(-π, π) => x -> 1/2 * (cos(x) + 1)
     )
 
     @test_throws DomainError window(123)
@@ -42,4 +42,63 @@ end
     @test isequal_interval(window(interval(-10, 10)), interval(0, 1))
     @test decoration(window(interval(-10, 10))) == trv
     @test isempty_interval(window(interval(100, 1000)))
-end
\ No newline at end of file
+end
+
+@testset "Derivatives" begin
+    slide = Piecewise(
+        Domain(-Inf, -1, false, true) => x -> -2x - 1,
+        Domain(-1, 0, false, true) => x -> x^2,
+        Domain(0, Inf, false, false) => Constant(0) ;
+        continuity = [1, 1]
+    )
+
+    @test ForwardDiff.derivative(slide, -5.5) == -2
+    @test ForwardDiff.derivative(slide, -0.5) == -1
+    @test ForwardDiff.derivative(slide, 1.2) == 0 
+
+    @test isequal_interval(
+        ForwardDiff.derivative(slide, interval(-7, -3)),
+        interval(-2)
+    )
+    @test isequal_interval(
+        ForwardDiff.derivative(slide, interval(-0.7, -0.3)),
+        interval(-1.4, -0.6)
+    )
+    @test isequal_interval(
+        ForwardDiff.derivative(slide, interval(0.7, 1.3)),
+        interval(0)
+    )
+
+    @test isequal_interval(
+        ForwardDiff.derivative(slide, interval(-1.7, -0.3)),
+        interval(-2, -0.6)
+    )
+
+    @test isequal_interval(
+        ForwardDiff.derivative(slide, interval(-0.7, 1.3)),
+        interval(-1.4, 0)
+    )
+
+    @test isequal_interval(
+        ForwardDiff.derivative(slide, interval(-1.7, 1.3)),
+        interval(-2, 0)
+    )
+
+    x1 = interval(-0.5, 0)
+    x2 = interval(-3, -2)
+
+    grad1 = ForwardDiff.gradient(xx -> slide(-xx[1]^2), [x1, x2])
+    grad2 = ForwardDiff.gradient(xx -> slide(0.7xx[2]), [x1, x2])
+
+    g1 = -2x1 * ForwardDiff.derivative(slide, -x1^2)
+    g2 = 0.7 * ForwardDiff.derivative(slide, x2)
+
+    @test isequal_interval(grad1[1], g1)
+    @test isequal_interval(grad1[2], interval(0))
+    @test isequal_interval(grad2[1], interval(0))
+    @test isequal_interval(grad2[2], g2)
+
+    grad = ForwardDiff.gradient(xx -> slide(-xx[1]^2 + 0.7xx[2]), [x1, x2])
+    @test isequal_interval(grad[1], g1)
+    @test isequal_interval(grad[2], g2)
+end

From fe77b77393b5234545707c68fb5b5ca67880840c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:08:55 +0100
Subject: [PATCH 08/20] Finish and fix tests

---
 ext/IntervalArithmeticForwardDiffExt.jl |  2 +-
 src/IntervalArithmetic.jl               |  1 +
 src/piecewise.jl                        |  4 +--
 test/interval_tests/piecewise.jl        | 35 ++++++++++++++++++++-----
 4 files changed, 32 insertions(+), 10 deletions(-)

diff --git a/ext/IntervalArithmeticForwardDiffExt.jl b/ext/IntervalArithmeticForwardDiffExt.jl
index 5017912f..b22977a8 100644
--- a/ext/IntervalArithmeticForwardDiffExt.jl
+++ b/ext/IntervalArithmeticForwardDiffExt.jl
@@ -110,7 +110,7 @@ function (piecewise::Piecewise)(dual::Dual{T, <:Interval}) where {T}
 
     if !isempty(setdiff(set, domain(piecewise)))
         dec = trv
-    elseif any(in(set), discontinuities(piecewise))
+    elseif any(in(set), discontinuities(piecewise, 1))
         dec = def 
     else
         dec = com
diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 926f72da..65f80a6b 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -24,6 +24,7 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 #
 
 include("piecewise.jl")
+    export Open, Closed, Domain
     export Constant, Piecewise, Domain, domain, subdomains, discontinuities, pieces
 
 #
diff --git a/src/piecewise.jl b/src/piecewise.jl
index d068c8b2..e7a1dc8e 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -1,5 +1,5 @@
 import Intervals
-import Intervals: IntervalSet, less_than_disjoint
+import Intervals: IntervalSet, less_than_disjoint, Open, Closed
 
 const Domain = Intervals.Interval
 
@@ -54,7 +54,7 @@ end
 subdomains(piecewise::Piecewise) = convert(Vector, piecewise.domain)
 domain(piecewise::Piecewise) = piecewise.domain
 pieces(piecewise::Piecewise) = zip(subdomains(piecewise), piecewise.fs)
-discontinuities(piecewise::Piecewise) = piecewise.singularities[piecewise.continuity .< 0]
+discontinuities(piecewise::Piecewise, order = 0) = piecewise.singularities[piecewise.continuity .< order]
 
 domain_string(x::Domain) = repr(x ; context = IOContext(stdout, :compact => true))
 
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index adfa3d4c..9917f382 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -1,7 +1,7 @@
 @testset "Step function" begin
     step = Piecewise(
-        Domain(-Inf, 0, false, true) => Constant(0),
-        Domain(0, Inf, false, false) => Constant(1)
+        Domain{Open, Closed}(-Inf, 0) => Constant(0),
+        Domain{Open, Open}(0, Inf) => Constant(1)
     )
 
     @test step(-1) == 0
@@ -16,8 +16,8 @@ end
 
 @testset "abs" begin
     myabs = Piecewise(
-        Domain(-Inf, 0, false, true) => x -> -x,
-        Domain(0, Inf, false, false) => identity ;
+        Domain{Open, Closed}(-Inf, 0) => x -> -x,
+        Domain{Open, Open}(0, Inf) => identity ;
         continuity = [0]
     )
 
@@ -46,9 +46,9 @@ end
 
 @testset "Derivatives" begin
     slide = Piecewise(
-        Domain(-Inf, -1, false, true) => x -> -2x - 1,
-        Domain(-1, 0, false, true) => x -> x^2,
-        Domain(0, Inf, false, false) => Constant(0) ;
+        Domain{Open, Closed}(-Inf, -1) => x -> -2x - 1,
+        Domain{Open, Closed}(-1, 0) => x -> x^2,
+        Domain{Open, Open}(0, Inf) => Constant(0) ;
         continuity = [1, 1]
     )
 
@@ -99,6 +99,27 @@ end
     @test isequal_interval(grad2[2], g2)
 
     grad = ForwardDiff.gradient(xx -> slide(-xx[1]^2 + 0.7xx[2]), [x1, x2])
+    g1 = -2x1 * ForwardDiff.derivative(slide, -x1^2 + 0.7x2)
+    g2 = 0.7 * ForwardDiff.derivative(slide, -x1^2 + 0.7x2)
     @test isequal_interval(grad[1], g1)
     @test isequal_interval(grad[2], g2)
 end
+
+@testset "Singularities" begin
+    f = Piecewise(
+        Domain{Open, Closed}(0, 1) => Constant(0),
+        Domain{Open, Closed}(1, 2) => x -> 0.5x,
+        Domain{Open, Closed}(2, 3) => Constant(1),
+        Domain{Open, Open}(3, 4) => x -> (x-3)^2 + 1 ;
+        continuity = [-1, 0, 1]
+    )
+
+    @test decoration(f(interval(0.5, 1.5))) == def
+    @test decoration(f(interval(1.5, 2.5))) == com
+    @test decoration(f(interval(2.5, 3.5))) == com
+
+    df = x -> ForwardDiff.derivative(f, x)
+    @test decoration(df(interval(0.5, 1.5))) == def
+    @test decoration(df(interval(1.5, 2.5))) == def
+    @test decoration(df(interval(2.5, 3.5))) == com
+end
\ No newline at end of file

From 7cec2d67467f2c9467fb1710a8177883076abc9f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 16:29:18 +0100
Subject: [PATCH 09/20] Add deps

---
 Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Project.toml b/Project.toml
index 0fcf1bc0..07df01e9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,6 +27,7 @@ CRlibm_jll = "1"
 DiffRules = "1"
 ForwardDiff = "0.10"
 IntervalSets = "0.7"
+Intervals = "1.10"
 LinearAlgebra = "1.9"
 MacroTools = "0.5"
 RecipesBase = "1"

From aa069af7e3d111e2e6ed352f585dc57c4cdcaceb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Thu, 6 Feb 2025 18:25:04 +0100
Subject: [PATCH 10/20] Remove Intervals.jl dep (not tested, need to go)

---
 src/piecewise.jl | 86 +++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 77 insertions(+), 9 deletions(-)

diff --git a/src/piecewise.jl b/src/piecewise.jl
index e7a1dc8e..88b643cb 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -1,10 +1,54 @@
-import Intervals
-import Intervals: IntervalSet, less_than_disjoint, Open, Closed
+struct Domain{T <: NumTypes, L, R}
+    lo::T
+    hi::T
+end
+
+Domain(lo, hi) = Domain{:open, :closed}(lo, hi)
+Domain(lo::Tuple, hi::Tuple) = Domain{lo[2], hi[2]}(lo[1], lo[2])
+Domain() = Domain{:open, :open}(Inf, -Inf)
+Domain(X::Interval) = Domain{:closed, :closed}(inf(X), sup(X))
+
+lowerbound(x::Domain{L, R}) where {L, R} = (x.lo, L)
+upperbound(x::Domain{L, R}) where {L, R} = (x.hi, R)
 
-const Domain = Intervals.Interval
+inf(x::Domain) = x.lo
+sup(x::Domain) = x.hi
 
-inf(x::Domain) = first(x)
-sup(x::Domain) = last(x)
+function rightof(x, (val, bound))
+    bound == :closed && return val <= x
+    return val < x
+end
+
+function rightof((val1, bound1), (val2, bound2))
+    if val1 < val2 || (val1 == val2 && bound1 == bound2 == :closed)
+        return false
+    end
+
+    return true
+end
+
+function leftof(x, (val, bound))
+    bound == :closed && return x <= val
+    return x < val
+end
+
+function leftof((val1, bound1), (val2, bound2))
+    if val2 < val1 || (val1 == val2 && bound1 == bound2 == :closed)
+        return false
+    end
+
+    return true
+end
+
+Base.in(x::Real, domain::Domain) = rightof(x, lowerbound(domain)) && leftof(x, upperbound(domain))
+
+function Base.intersect(d1::Domain, d2::Domain)
+    left = max(lowerbound(d1), lowerbound(d2))
+    right = min(upperbound(d1), upperbound(d2))
+
+    left > right && return Domain()
+    return Domain(left, right)
+end
 
 struct Constant{T}
     value::T
@@ -38,7 +82,9 @@ function Piecewise(
         s1 = subdomains[k]
         s2 = subdomains[k + 1]
 
-        !less_than_disjoint(s1, s2) && throw(ArgumentError("domains are either not ordered or not disjoint"))
+        if rightof(upperbound(s1), lowerbound(s2))
+           throw(ArgumentError("domains are either not ordered or not disjoint"))
+        end
     end
 
     singularities = sup.(subdomains[1:end-1])
@@ -76,10 +122,32 @@ function Base.show(io::IO, ::MIME"text/plain", piecewise::Piecewise)
 end
 
 function (piecewise::Piecewise)(X::Interval{T}) where {T}
-    set = Domain(inf(X), sup(X), true, true)
-    isdisjoint(set, domain(piecewise)) && return emptyinterval(T)
+    set = Domain(X)
+    all(isempty, intersect.(Ref(set), subdomains(piecewise))) && return emptyinterval(T)
+
+    # This logic exploits the fact that subdomains are ordered
+    lo = lowerbound(set)
+    istrivial = false
+    for subdomain in subdomains(piecewise)
+        if leftof(lo, lowerbound(subdomain))
+            istrivial = true
+            break
+        end
+
+        val, bound = lowerbound(subdomain)
+
+        if bound == :closed
+            lo = (val, :open)
+        else
+            lo = (val, :closed)
+        end
+    end
+
+    if rightof(upperbound(set), upperbound(subdomains(piecewise)[end]))
+        istrivial = true
+    end
 
-    if !isempty(setdiff(set, domain(piecewise)))
+    if istrivial
         dec = trv
     elseif any(in(set), discontinuities(piecewise))
         dec = def 

From fde8b0a7875b12e164743a0583c0d8dc709193b9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Fri, 7 Feb 2025 16:01:36 +0100
Subject: [PATCH 11/20] Fix everything (except the forwarddiff extension)

---
 ext/IntervalArithmeticForwardDiffExt.jl |  4 --
 src/IntervalArithmetic.jl               |  3 +-
 src/piecewise.jl                        | 56 ++++++++++++++-----
 test/interval_tests/piecewise.jl        | 71 +++++++++++++++++++++----
 4 files changed, 105 insertions(+), 29 deletions(-)

diff --git a/ext/IntervalArithmeticForwardDiffExt.jl b/ext/IntervalArithmeticForwardDiffExt.jl
index b22977a8..9072fc25 100644
--- a/ext/IntervalArithmeticForwardDiffExt.jl
+++ b/ext/IntervalArithmeticForwardDiffExt.jl
@@ -90,10 +90,6 @@ function Base.:(^)(x::ExactReal, y::Dual{<:Ty}) where {Ty}
     end
 end
 
-# resolve ambiguity
-
-Base.convert(::Type{Dual{T,V,N}}, x::ExactReal) where {T,V,N} = Dual{T}(V(x), zero(Partials{N,V}))
-
 
 # Piecewise functions
 
diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 65f80a6b..4eefd41f 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -24,8 +24,7 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 #
 
 include("piecewise.jl")
-    export Open, Closed, Domain
-    export Constant, Piecewise, Domain, domain, subdomains, discontinuities, pieces
+    export Domain, Constant, Piecewise, Domain, domain, subdomains, discontinuities, pieces
 
 #
 
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 88b643cb..37d2eb39 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -1,12 +1,13 @@
-struct Domain{T <: NumTypes, L, R}
+struct Domain{L, R, T, S}
     lo::T
-    hi::T
+    hi::S
 end
 
+Domain{L, R}(lo::T, hi::S) where {T, S, L, R} = Domain{L, R, T, S}(lo, hi)
 Domain(lo, hi) = Domain{:open, :closed}(lo, hi)
-Domain(lo::Tuple, hi::Tuple) = Domain{lo[2], hi[2]}(lo[1], lo[2])
-Domain() = Domain{:open, :open}(Inf, -Inf)
+Domain(lo::Tuple, hi::Tuple) = Domain{lo[2], hi[2]}(lo[1], hi[1])
 Domain(X::Interval) = Domain{:closed, :closed}(inf(X), sup(X))
+Domain() = Domain{:open, :open}(Inf, -Inf)
 
 lowerbound(x::Domain{L, R}) where {L, R} = (x.lo, L)
 upperbound(x::Domain{L, R}) where {L, R} = (x.hi, R)
@@ -14,7 +15,12 @@ upperbound(x::Domain{L, R}) where {L, R} = (x.hi, R)
 inf(x::Domain) = x.lo
 sup(x::Domain) = x.hi
 
-function rightof(x, (val, bound))
+"""
+    rightof(val, lowerbound)
+
+Determine if a value is on the right of a lower bound.
+"""
+function rightof(x::Real, (val, bound))
     bound == :closed && return val <= x
     return val < x
 end
@@ -27,7 +33,12 @@ function rightof((val1, bound1), (val2, bound2))
     return true
 end
 
-function leftof(x, (val, bound))
+"""
+    leftof(val, upperbound)
+
+Determine if a value is on the left of an upper bound.
+"""
+function leftof(x::Real, (val, bound))
     bound == :closed && return x <= val
     return x < val
 end
@@ -40,6 +51,14 @@ function leftof((val1, bound1), (val2, bound2))
     return true
 end
 
+function leftof(d1::Domain, d2::Domain)
+    val1, bound1 = upperbound(d1)
+    val2, bound2 = lowerbound(d2)
+
+    val1 == val2 && return !(bound1 == bound2 == :closed)
+    return val1 < val2
+end
+
 Base.in(x::Real, domain::Domain) = rightof(x, lowerbound(domain)) && leftof(x, upperbound(domain))
 
 function Base.intersect(d1::Domain, d2::Domain)
@@ -50,6 +69,14 @@ function Base.intersect(d1::Domain, d2::Domain)
     return Domain(left, right)
 end
 
+function Base.isempty(domain::Domain)
+    lo, lobound = lowerbound(domain)
+    hi, hibound = upperbound(domain)
+
+    lo == hi && return !(lobound == hibound == :closed)
+    return lo > hi
+end
+
 struct Constant{T}
     value::T
 end
@@ -57,8 +84,9 @@ end
 (constant::Constant)(::Any) = constant.value
 (constant::Constant)(::Interval) = interval(constant.value)
 
+# TODO Proper type parameters
 struct Piecewise
-    domain::IntervalSet
+    domain::Vector
     fs
     continuity::Vector{Int}
     singularities::Vector
@@ -82,14 +110,14 @@ function Piecewise(
         s1 = subdomains[k]
         s2 = subdomains[k + 1]
 
-        if rightof(upperbound(s1), lowerbound(s2))
+        if !leftof(s1, s2)
            throw(ArgumentError("domains are either not ordered or not disjoint"))
         end
     end
 
     singularities = sup.(subdomains[1:end-1])
 
-    return Piecewise(IntervalSet(subdomains), fs, continuity, singularities)
+    return Piecewise(subdomains, fs, continuity, singularities)
 end
 
 function Piecewise(pairs::Vararg{<:Pair} ; continuity = fill(-1, length(pairs) - 1))
@@ -104,7 +132,7 @@ discontinuities(piecewise::Piecewise, order = 0) = piecewise.singularities[piece
 
 domain_string(x::Domain) = repr(x ; context = IOContext(stdout, :compact => true))
 
-function domain_string(S::IntervalSet)
+function domain_string(S::Vector{<:Domain})
     r = repr("text/plain", S ; context = IOContext(stdout, :compact => true))
     return join(split(r, "\n")[2:end], " ∪ ")
 end
@@ -128,13 +156,17 @@ function (piecewise::Piecewise)(X::Interval{T}) where {T}
     # This logic exploits the fact that subdomains are ordered
     lo = lowerbound(set)
     istrivial = false
+
     for subdomain in subdomains(piecewise)
-        if leftof(lo, lowerbound(subdomain))
+
+        if !rightof(lo, lowerbound(subdomain))
             istrivial = true
             break
         end
 
-        val, bound = lowerbound(subdomain)
+        val, bound = upperbound(subdomain)
+
+        val > upperbound(set)[1] && break
 
         if bound == :closed
             lo = (val, :open)
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index 9917f382..5f0d3b2c 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -1,7 +1,55 @@
+using IntervalArithmetic: leftof
+
+@testset "Domain" begin
+    @testset "leftof" begin
+        d1 = Domain{:open, :closed}(0, 1)
+        d2 = Domain{:open, :open}(0, 1)
+        d3 = Domain{:open, :closed}(1, 2)
+        d4 = Domain{:closed, :closed}(1, 2)
+
+        @test !leftof(d1, d2)
+        @test leftof(d1, d3)
+        @test !leftof(d1, d4)
+
+        @test leftof(d2, d3)
+        @test leftof(d2, d4)
+    end
+
+    @testset "intersect" begin
+        d1 = Domain{:closed, :open}(0, 10)
+        d2 = Domain{:closed, :open}(2, 15)
+        d3 = Domain{:open, :closed}(4, 7)
+        d4 = Domain{:open, :closed}(-20, 3)
+
+        @test intersect(d1, d2) == Domain{:closed, :open}(2, 10)
+        @test intersect(d1, d3) == Domain{:open, :closed}(4, 7)
+        @test intersect(d1, d4) == Domain{:closed, :closed}(0, 3)
+        @test intersect(d2, d3) == Domain{:open, :closed}(4, 7)
+        @test intersect(d2, d4) == Domain{:closed, :closed}(2, 3)
+        @test intersect(d3, d4) == Domain()
+    end
+
+    @testset "isempty" begin
+        @test isempty(Domain{:open, :open}(1, 1))
+        @test !isempty(Domain{:closed, :closed}(1, 1))
+        @test !isempty(Domain{:open, :open}(1, 2))
+        @test isempty(Domain())
+    end
+end
+
+@testset "Piecewise constructor" begin
+    d1 = Domain{:open, :closed}(0, 1)
+    d2 = Domain{:open, :closed}(1, 2)
+    d3 = Domain{:closed, :closed}(1, 2)
+
+    @test_throws ArgumentError Piecewise(d2 => Constant(0), d1 => Constant(1))
+    @test_throws ArgumentError Piecewise(d1 => Constant(0), d3 => Constant(1))
+end
+
 @testset "Step function" begin
     step = Piecewise(
-        Domain{Open, Closed}(-Inf, 0) => Constant(0),
-        Domain{Open, Open}(0, Inf) => Constant(1)
+        Domain{:open, :closed}(-Inf, 0) => Constant(0),
+        Domain{:open, :open}(0, 1000) => Constant(1)
     )
 
     @test step(-1) == 0
@@ -12,12 +60,13 @@
     @test decoration(step(interval(4.44))) == com
     @test isequal_interval(step(interval(-22.2, 33.3)), interval(0, 1))
     @test decoration(step(interval(-11, 11))) == def
+    @test decoration(step(interval(500, 2000))) == trv
 end
 
 @testset "abs" begin
     myabs = Piecewise(
-        Domain{Open, Closed}(-Inf, 0) => x -> -x,
-        Domain{Open, Open}(0, Inf) => identity ;
+        Domain{:open, :closed}(-Inf, 0) => x -> -x,
+        Domain{:open, :open}(0, Inf) => identity ;
         continuity = [0]
     )
 
@@ -46,9 +95,9 @@ end
 
 @testset "Derivatives" begin
     slide = Piecewise(
-        Domain{Open, Closed}(-Inf, -1) => x -> -2x - 1,
-        Domain{Open, Closed}(-1, 0) => x -> x^2,
-        Domain{Open, Open}(0, Inf) => Constant(0) ;
+        Domain{:open, :closed}(-Inf, -1) => x -> -2x - 1,
+        Domain{:open, :closed}(-1, 0) => x -> x^2,
+        Domain{:open, :open}(0, Inf) => Constant(0) ;
         continuity = [1, 1]
     )
 
@@ -107,10 +156,10 @@ end
 
 @testset "Singularities" begin
     f = Piecewise(
-        Domain{Open, Closed}(0, 1) => Constant(0),
-        Domain{Open, Closed}(1, 2) => x -> 0.5x,
-        Domain{Open, Closed}(2, 3) => Constant(1),
-        Domain{Open, Open}(3, 4) => x -> (x-3)^2 + 1 ;
+        Domain{:open, :closed}(0, 1) => Constant(0),
+        Domain{:open, :closed}(1, 2) => x -> 0.5x,
+        Domain{:open, :closed}(2, 3) => Constant(1),
+        Domain{:open, :open}(3, 4) => x -> (x-3)^2 + 1 ;
         continuity = [-1, 0, 1]
     )
 

From 81492c7de188783940fc9751d422a35f1e7a6edc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Sat, 8 Feb 2025 01:02:31 +0100
Subject: [PATCH 12/20] Fix the derivative too

---
 Project.toml                            |  2 --
 ext/IntervalArithmeticForwardDiffExt.jl |  6 ++---
 src/piecewise.jl                        | 32 +++++++++++++------------
 3 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/Project.toml b/Project.toml
index 07df01e9..75dd7071 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ version = "0.22.22"
 
 [deps]
 CRlibm_jll = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
-Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
@@ -27,7 +26,6 @@ CRlibm_jll = "1"
 DiffRules = "1"
 ForwardDiff = "0.10"
 IntervalSets = "0.7"
-Intervals = "1.10"
 LinearAlgebra = "1.9"
 MacroTools = "0.5"
 RecipesBase = "1"
diff --git a/ext/IntervalArithmeticForwardDiffExt.jl b/ext/IntervalArithmeticForwardDiffExt.jl
index 9072fc25..a6dc5940 100644
--- a/ext/IntervalArithmeticForwardDiffExt.jl
+++ b/ext/IntervalArithmeticForwardDiffExt.jl
@@ -99,12 +99,12 @@ end
 
 function (piecewise::Piecewise)(dual::Dual{T, <:Interval}) where {T}
     X = value(dual)
-    set = Domain(inf(X), sup(X), true, true)
-    if isdisjoint(set, domain(piecewise))
+    set = Domain(X)
+    if !IntervalArithmetic.overlapdomain(set, piecewise) 
         return Dual{T}(emptyinterval(X), emptyinterval(X) .* partials(dual))
     end
 
-    if !isempty(setdiff(set, domain(piecewise)))
+    if !IntervalArithmetic.indomain(set, piecewise)
         dec = trv
     elseif any(in(set), discontinuities(piecewise, 1))
         dec = def 
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 37d2eb39..e3ad454d 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -77,6 +77,7 @@ function Base.isempty(domain::Domain)
     return lo > hi
 end
 
+
 struct Constant{T}
     value::T
 end
@@ -120,7 +121,7 @@ function Piecewise(
     return Piecewise(subdomains, fs, continuity, singularities)
 end
 
-function Piecewise(pairs::Vararg{<:Pair} ; continuity = fill(-1, length(pairs) - 1))
+function Piecewise(pairs::Vararg{Pair} ; continuity = fill(-1, length(pairs) - 1))
     pairs = collect(pairs)
     return Piecewise(first.(pairs), last.(pairs), continuity)
 end
@@ -149,24 +150,20 @@ function Base.show(io::IO, ::MIME"text/plain", piecewise::Piecewise)
     end
 end
 
-function (piecewise::Piecewise)(X::Interval{T}) where {T}
-    set = Domain(X)
-    all(isempty, intersect.(Ref(set), subdomains(piecewise))) && return emptyinterval(T)
+function indomain(domain, piecewise)
+    rightof(upperbound(domain), upperbound(subdomains(piecewise)[end])) && return false
 
-    # This logic exploits the fact that subdomains are ordered
-    lo = lowerbound(set)
-    istrivial = false
+    # This relies on the fact that subdomains are ordered
+    lo = lowerbound(domain)
 
     for subdomain in subdomains(piecewise)
-
         if !rightof(lo, lowerbound(subdomain))
-            istrivial = true
-            break
+            return false
         end
 
         val, bound = upperbound(subdomain)
 
-        val > upperbound(set)[1] && break
+        val > upperbound(domain)[1] && break
 
         if bound == :closed
             lo = (val, :open)
@@ -175,11 +172,16 @@ function (piecewise::Piecewise)(X::Interval{T}) where {T}
         end
     end
 
-    if rightof(upperbound(set), upperbound(subdomains(piecewise)[end]))
-        istrivial = true
-    end
+    return true
+end
+
+overlapdomain(domain, piecewise) = any(!isempty, intersect.(Ref(domain), subdomains(piecewise)))
+
+function (piecewise::Piecewise)(X::Interval{T}) where {T}
+    set = Domain(X)
+    !overlapdomain(set, piecewise) && return emptyinterval(T)
 
-    if istrivial
+    if !indomain(set, piecewise)
         dec = trv
     elseif any(in(set), discontinuities(piecewise))
         dec = def 

From 51db09e697d9edd52ff1f4b60082b07f4234ad1c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Sun, 9 Feb 2025 19:13:37 +0100
Subject: [PATCH 13/20] Fix interval_diff

---
 .../interval_operations/set_operations.jl     |  2 +-
 src/intervals/intervals.jl                    |  2 +-
 test/interval_tests/set_operations.jl         | 19 +++++++++++++++++++
 3 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/src/intervals/interval_operations/set_operations.jl b/src/intervals/interval_operations/set_operations.jl
index 792d36c5..abffa339 100644
--- a/src/intervals/interval_operations/set_operations.jl
+++ b/src/intervals/interval_operations/set_operations.jl
@@ -204,7 +204,7 @@ end
 
 function interval_diff(x::Interval{T}, y::Interval) where {T<:NumTypes}
     isdisjoint_interval(x, y) && return [x]
-    isequal_interval(x, y) && return Interval{T}[]
+    issubset_interval(x, y) && return Interval{T}[]
 
     intersection = intersect_interval(x, y)
     inf(x) == inf(intersection) && return [interval(sup(intersection), sup(x))]
diff --git a/src/intervals/intervals.jl b/src/intervals/intervals.jl
index cb4685e9..bb5d2cd0 100644
--- a/src/intervals/intervals.jl
+++ b/src/intervals/intervals.jl
@@ -41,7 +41,7 @@ include("interval_operations/overlap.jl")
 include("interval_operations/numeric.jl")
     export inf, sup, bounds, mid, diam, radius, midradius, mag, mig, dist
 include("interval_operations/set_operations.jl")
-    export intersect_interval, hull, interiordiff
+    export intersect_interval, hull, interiordiff, interval_diff
 include("interval_operations/bisect.jl")
     export bisect, mince, mince!
 
diff --git a/test/interval_tests/set_operations.jl b/test/interval_tests/set_operations.jl
index 4b5c28bd..d415b12b 100644
--- a/test/interval_tests/set_operations.jl
+++ b/test/interval_tests/set_operations.jl
@@ -37,3 +37,22 @@ end
     @test interiordiff(x, z) == Interval{Float64}[]
     @test all(isequal_interval.(interiordiff(z, x), [interval(0, 1), interval(3, 5)]))
 end
+
+@testset "interval_diff" begin
+    A, B = interval_diff(interval(1, 10), interval(2, 5))
+    @test isequal_interval(A, interval(1, 2))
+    @test isequal_interval(B, interval(5, 10))
+
+    @test isequal_interval(
+        only(interval_diff(interval(1, 10), interval(1, 5))),
+        interval(5, 10)
+    )
+
+    @test isequal_interval(
+        only(interval_diff(interval(1, 10), interval(7, 12))),
+        interval(1, 7)
+    )
+
+    @test interval_diff(interval(1, 10), interval(-1, 14)) == []
+    @test interval_diff(interval(1, 10), interval(1, 10)) == []
+end
\ No newline at end of file

From af320698ff5ebf7ae6aba3e7f807b72b61db7576 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Sun, 9 Feb 2025 20:02:12 +0100
Subject: [PATCH 14/20] Check bounds type at construction

---
 src/piecewise.jl                 | 7 +++++++
 test/interval_tests/piecewise.jl | 5 +++++
 2 files changed, 12 insertions(+)

diff --git a/src/piecewise.jl b/src/piecewise.jl
index e3ad454d..01e7d1a8 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -1,6 +1,13 @@
 struct Domain{L, R, T, S}
     lo::T
     hi::S
+
+    function Domain{L, R, T, S}(lo::T, hi::S) where {L, R, T, S}
+        (!(L in (:open, :closed)) || !(R in (:open, :closed))) && throw(ErrorException(
+            "Domain bound must be either :open or :closed, got $L and $R instead"
+        ))
+        return new{L, R, T, S}(lo, hi)
+    end
 end
 
 Domain{L, R}(lo::T, hi::S) where {T, S, L, R} = Domain{L, R, T, S}(lo, hi)
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index 5f0d3b2c..0cafa398 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -1,6 +1,11 @@
 using IntervalArithmetic: leftof
 
 @testset "Domain" begin
+    @testset "Construction" begin
+        @test_throws ErrorException Domain{:oopen, :closed}(0, 1)
+        @test_throws ErrorException Domain{:open, :close}(0, 1)
+    end
+
     @testset "leftof" begin
         d1 = Domain{:open, :closed}(0, 1)
         d2 = Domain{:open, :open}(0, 1)

From dbe04631f5813f4294d85996a8fe90ac7d12b162 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Sun, 9 Feb 2025 21:49:04 +0100
Subject: [PATCH 15/20] Fix exports

---
 src/IntervalArithmetic.jl             | 2 +-
 src/intervals/intervals.jl            | 2 +-
 test/interval_tests/set_operations.jl | 2 ++
 3 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 4eefd41f..d60049c1 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -24,7 +24,7 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 #
 
 include("piecewise.jl")
-    export Domain, Constant, Piecewise, Domain, domain, subdomains, discontinuities, pieces
+    export Domain, Constant, Piecewise, domain, subdomains, discontinuities, pieces
 
 #
 
diff --git a/src/intervals/intervals.jl b/src/intervals/intervals.jl
index bb5d2cd0..cb4685e9 100644
--- a/src/intervals/intervals.jl
+++ b/src/intervals/intervals.jl
@@ -41,7 +41,7 @@ include("interval_operations/overlap.jl")
 include("interval_operations/numeric.jl")
     export inf, sup, bounds, mid, diam, radius, midradius, mag, mig, dist
 include("interval_operations/set_operations.jl")
-    export intersect_interval, hull, interiordiff, interval_diff
+    export intersect_interval, hull, interiordiff
 include("interval_operations/bisect.jl")
     export bisect, mince, mince!
 
diff --git a/test/interval_tests/set_operations.jl b/test/interval_tests/set_operations.jl
index d415b12b..eed5e45a 100644
--- a/test/interval_tests/set_operations.jl
+++ b/test/interval_tests/set_operations.jl
@@ -1,3 +1,5 @@
+using IntervalArithmetic: interval_diff
+
 @testset "removed interval" begin
     @test_throws ArgumentError intersect(interval(1))
     @test_throws ArgumentError intersect(interval(1), 2, [1], 4., 5)

From 86564e3141f45ab156a3fd3cf5e25ba768671a62 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Sun, 9 Feb 2025 22:40:58 +0100
Subject: [PATCH 16/20] A bit of code cleaning

---
 src/IntervalArithmetic.jl        |  2 +-
 src/piecewise.jl                 | 90 ++++++++++++++++++--------------
 test/interval_tests/piecewise.jl |  2 +-
 3 files changed, 54 insertions(+), 40 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index d60049c1..b00b6734 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -24,7 +24,7 @@ const RealIntervalType{T} = Union{BareInterval{T},Interval{T}}
 #
 
 include("piecewise.jl")
-    export Domain, Constant, Piecewise, domain, subdomains, discontinuities, pieces
+    export Domain, Constant, Piecewise, domains, discontinuities, pieces
 
 #
 
diff --git a/src/piecewise.jl b/src/piecewise.jl
index 01e7d1a8..ec817bfd 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -11,7 +11,6 @@ struct Domain{L, R, T, S}
 end
 
 Domain{L, R}(lo::T, hi::S) where {T, S, L, R} = Domain{L, R, T, S}(lo, hi)
-Domain(lo, hi) = Domain{:open, :closed}(lo, hi)
 Domain(lo::Tuple, hi::Tuple) = Domain{lo[2], hi[2]}(lo[1], hi[1])
 Domain(X::Interval) = Domain{:closed, :closed}(inf(X), sup(X))
 Domain() = Domain{:open, :open}(Inf, -Inf)
@@ -92,40 +91,51 @@ end
 (constant::Constant)(::Any) = constant.value
 (constant::Constant)(::Interval) = interval(constant.value)
 
-# TODO Proper type parameters
-struct Piecewise
-    domain::Vector
-    fs
-    continuity::Vector{Int}
-    singularities::Vector
+struct Piecewise{N, D<:Tuple, F<:Tuple, S<:Tuple}
+    domains::D
+    fs::F
+    continuity::NTuple{N, Int}
+    singularities::S
+
+    function Piecewise(domains::D, fs::F, continuity::NTuple{N, Int}, singularities::S) where {D, F, N, S}
+        if !(N + 1 == length(domains) == length(fs))
+            throw(ArgumentError(
+                "a Piecewise function must have as many domains as functions, " *
+                "and one less continuity point. " *
+                "Given: $(length(domains)) domains, $(length(fs)) functions, " *
+                "$N continuity points."))
+        end
+
+        return new{N, D, F, S}(domains, fs, continuity, singularities)
+    end
 end
 
 function Piecewise(
-        subdomains::Vector{<:Domain},
+        domains::Vector{<:Domain},
         fs,
-        continuity::Vector{Int} = fill(-1, length(subdomains) - 1))
+        continuity::Vector{Int} = fill(-1, length(domains) - 1))
 
-    if length(subdomains) != length(fs)
+    if length(domains) != length(fs)
         throw(ArgumentError("the number of domains and the number of functions don't match"))
     end
     
-    if length(subdomains) - 1 != length(continuity)
-        n = length(subdomains)
+    if length(domains) - 1 != length(continuity)
+        n = length(domains)
         throw(ArgumentError("$(length(sub)) junction points but $(n - 1) are expected based on the number of domains ($n)"))
     end
 
-    for k in 1:length(subdomains) - 1
-        s1 = subdomains[k]
-        s2 = subdomains[k + 1]
+    for k in 1:length(domains) - 1
+        s1 = domains[k]
+        s2 = domains[k + 1]
 
         if !leftof(s1, s2)
            throw(ArgumentError("domains are either not ordered or not disjoint"))
         end
     end
 
-    singularities = sup.(subdomains[1:end-1])
+    singularities = sup.(domains[1:end-1])
 
-    return Piecewise(subdomains, fs, continuity, singularities)
+    return Piecewise(Tuple(domains), Tuple(fs), Tuple(continuity), Tuple(singularities))
 end
 
 function Piecewise(pairs::Vararg{Pair} ; continuity = fill(-1, length(pairs) - 1))
@@ -133,42 +143,46 @@ function Piecewise(pairs::Vararg{Pair} ; continuity = fill(-1, length(pairs) - 1
     return Piecewise(first.(pairs), last.(pairs), continuity)
 end
 
-subdomains(piecewise::Piecewise) = convert(Vector, piecewise.domain)
-domain(piecewise::Piecewise) = piecewise.domain
-pieces(piecewise::Piecewise) = zip(subdomains(piecewise), piecewise.fs)
-discontinuities(piecewise::Piecewise, order = 0) = piecewise.singularities[piecewise.continuity .< order]
+domains(piecewise::Piecewise) = piecewise.domains
+pieces(piecewise::Piecewise) = zip(domains(piecewise), piecewise.fs)
 
-domain_string(x::Domain) = repr(x ; context = IOContext(stdout, :compact => true))
+function discontinuities(piecewise::Piecewise, order = 0)
+    return [s for (s, C) in zip(piecewise.singularities, piecewise.continuity) if C .< order]
+end
+
+function domain_string(domain::Domain{L, R}) where {L, R}
+    left = (L == :closed) ? "[" : "("
+    right = (R == :closed) ? "]" : ")"
 
-function domain_string(S::Vector{<:Domain})
-    r = repr("text/plain", S ; context = IOContext(stdout, :compact => true))
-    return join(split(r, "\n")[2:end], " ∪ ")
+    return "$left$(domain.lo), $(domain.hi)$right"
 end
 
-domain_string(piecewise::Piecewise) = domain_string(domain(piecewise))
+function domain_string(piecewise::Piecewise)
+    join(domain_string.(domains(piecewise)), " ∪ ")
+end
 
 function Base.show(io::IO, ::MIME"text/plain", piecewise::Piecewise)
     n = length(pieces(piecewise))
     print(io, "Piecewise function with $n pieces:")
     
-    for (subdomain, f) in pieces(piecewise)
+    for (domain, f) in pieces(piecewise)
         println(io)
-        print(io, "  $(domain_string(subdomain)) -> $(repr(f))")
+        print(io, "  $(domain_string(domain)) -> $(repr(f))")
     end
 end
 
 function indomain(domain, piecewise)
-    rightof(upperbound(domain), upperbound(subdomains(piecewise)[end])) && return false
+    rightof(upperbound(domain), upperbound(domains(piecewise)[end])) && return false
 
-    # This relies on the fact that subdomains are ordered
+    # This relies on the fact that domains are ordered
     lo = lowerbound(domain)
 
-    for subdomain in subdomains(piecewise)
-        if !rightof(lo, lowerbound(subdomain))
+    for domain in domains(piecewise)
+        if !rightof(lo, lowerbound(domain))
             return false
         end
 
-        val, bound = upperbound(subdomain)
+        val, bound = upperbound(domain)
 
         val > upperbound(domain)[1] && break
 
@@ -182,7 +196,7 @@ function indomain(domain, piecewise)
     return true
 end
 
-overlapdomain(domain, piecewise) = any(!isempty, intersect.(Ref(domain), subdomains(piecewise)))
+overlapdomain(domain, piecewise) = any(!isempty, intersect.(Ref(domain), domains(piecewise)))
 
 function (piecewise::Piecewise)(X::Interval{T}) where {T}
     set = Domain(X)
@@ -197,8 +211,8 @@ function (piecewise::Piecewise)(X::Interval{T}) where {T}
     end
 
     results = Interval{T}[]
-    for (subdomain, f) in pieces(piecewise)
-        subset = intersect(set, subdomain)
+    for (domain, f) in pieces(piecewise)
+        subset = intersect(set, domain)
         isempty(subset) && continue
         push!(results, f(interval(inf(subset), sup(subset), decoration(X))))
     end
@@ -208,8 +222,8 @@ function (piecewise::Piecewise)(X::Interval{T}) where {T}
 end
 
 function (piecewise::Piecewise)(x::Real)
-    for (subdomain, f) in pieces(piecewise)
-        (x in subdomain) && return f(x)
+    for (domain, f) in pieces(piecewise)
+        (x in domain) && return f(x)
     end
     throw(DomainError(x, "piecewise function was called outside of its domain $(domain_string(piecewise))"))
 end
\ No newline at end of file
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index 0cafa398..b98df54f 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -87,7 +87,7 @@ end
 
 @testset "Out of domain" begin
     window = Piecewise(
-        Domain(-π, π) => x -> 1/2 * (cos(x) + 1)
+        Domain{:open, :closed}(-π, π) => x -> 1/2 * (cos(x) + 1)
     )
 
     @test_throws DomainError window(123)

From 5078432b3eaa38aa0f174a2af4eec582da23a7fb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Sun, 9 Feb 2025 23:16:17 +0100
Subject: [PATCH 17/20] Doc doc doc

---
 docs/src/manual/usage.md | 17 ++++++++
 src/piecewise.jl         | 83 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 100 insertions(+)

diff --git a/docs/src/manual/usage.md b/docs/src/manual/usage.md
index f0bb645c..94332d83 100644
--- a/docs/src/manual/usage.md
+++ b/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
 
diff --git a/src/piecewise.jl b/src/piecewise.jl
index ec817bfd..d722167a 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -1,3 +1,14 @@
+"""
+    Domain{LeftBound, RightBound}(lo, hi)
+
+The domain of a function.
+
+`LeftBound` and `RightBound` must be symbols and are either `:closed` or
+`:open` determining if the corresponding endpoint is (respectively)
+included or not in the domain.
+
+If `hi > lo`, the domain is considered to be empty.
+"""
 struct Domain{L, R, T, S}
     lo::T
     hi::S
@@ -83,7 +94,30 @@ function Base.isempty(domain::Domain)
     return lo > hi
 end
 
+"""
+    Constant(value)
+
+A constant function compatible with interval arithmetic.
 
+Return an interval containing only the value for an interval input,
+and the value directly otherwise.
+
+```jldoctest
+julia> c = Constant(1.2)
+Constant{Float64}(1.2)
+
+julia> c(22.2)
+1.2
+
+julia> c(interval(0, 1.3))
+[1.19999, 1.20001]_com
+```
+
+Note that this is not equivalent to `Returns(value)` from base,
+which always outputs `value`, even for an interval input.
+This can shortcircuit the propagation of intervals in the computation
+and lose the associated guarantee of correctness.
+"""
 struct Constant{T}
     value::T
 end
@@ -91,6 +125,55 @@ end
 (constant::Constant)(::Any) = constant.value
 (constant::Constant)(::Interval) = interval(constant.value)
 
+"""
+    Piecewise(pairs... ; continuity = fill(-1, length(pairs) - 1))
+
+A function defined by pieces (each associating a domain to a function).
+Support both intervals and standard numbers.
+
+```jldoctest
+julia> myabs = Piecewise(
+          Domain{:open, :closed}(-Inf, 0) => x -> -x,
+          Domain{:open, :open}(0, Inf) => identity
+       )
+Piecewise function with 2 pieces:
+  (-Inf, 0] -> var"#119#120"()
+  (0, Inf) -> identity
+
+julia> myabs(-22.3)
+22.3
+
+julia> myabs(interval(-5, 5))
+[0.0, 5.0]_def
+```
+
+For constant pieces, it is recommended to use `Constant`
+for full compatibility with intervals.
+
+The domains must be specified in increasing order and must not overlap.
+
+The `continuity` optional argument takes a vector of `N - 1` integers
+(where  `N` is the number of domains)
+determining how the piecewise function behaves at the endpoints between
+the subdomains.
+The possibility are:
+- `-1` : the function is discontinuous between the domains.
+- `0` : the function is continuous but not differentiable between the domains.
+- `n > 0` : the function is `n` times continuously differentiable between the
+    domains. This only matter when using `ForwardDiff` to compute derivative
+    of the function.
+
+This information is used to determine the decoration of intervals that
+covers the endpoint of several domains.
+
+If an input interval goes outside the domain of definition of the piecewise
+function, the output will always have the trivial (`trv`) decoration.
+For standard number, it throws a `DomainError`.
+
+The piecewise function can have a gap between two pieces.
+In this case, the `continuity` optional argument is ignored,
+and interval spanning over the gap always as the `trv` decoration.
+"""
 struct Piecewise{N, D<:Tuple, F<:Tuple, S<:Tuple}
     domains::D
     fs::F

From d29aec39f572ee3640396b13d92d832be78ff919 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Sun, 9 Feb 2025 23:35:33 +0100
Subject: [PATCH 18/20] Fix jldoctest (when possible)

---
 src/piecewise.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/piecewise.jl b/src/piecewise.jl
index d722167a..e8d6ee5f 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -110,7 +110,7 @@ julia> c(22.2)
 1.2
 
 julia> c(interval(0, 1.3))
-[1.19999, 1.20001]_com
+Interval{Float64}(1.2, 1.2, com)
 ```
 
 Note that this is not equivalent to `Returns(value)` from base,
@@ -131,7 +131,7 @@ end
 A function defined by pieces (each associating a domain to a function).
 Support both intervals and standard numbers.
 
-```jldoctest
+```jl
 julia> myabs = Piecewise(
           Domain{:open, :closed}(-Inf, 0) => x -> -x,
           Domain{:open, :open}(0, Inf) => identity
@@ -144,7 +144,7 @@ julia> myabs(-22.3)
 22.3
 
 julia> myabs(interval(-5, 5))
-[0.0, 5.0]_def
+Interval{Float64}(0.0, 5.0, def)
 ```
 
 For constant pieces, it is recommended to use `Constant`

From aa65113dc207f254e70c0325ceba6b3cd6adac23 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Mon, 10 Feb 2025 11:40:16 +0100
Subject: [PATCH 19/20] Exploit NTuple better

---
 src/piecewise.jl                 | 38 ++++++++++++++++++--------------
 test/interval_tests/piecewise.jl |  9 ++++++++
 2 files changed, 30 insertions(+), 17 deletions(-)

diff --git a/src/piecewise.jl b/src/piecewise.jl
index e8d6ee5f..e53470dc 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -174,29 +174,31 @@ The piecewise function can have a gap between two pieces.
 In this case, the `continuity` optional argument is ignored,
 and interval spanning over the gap always as the `trv` decoration.
 """
-struct Piecewise{N, D<:Tuple, F<:Tuple, S<:Tuple}
+struct Piecewise{N, M, D<:NTuple{N, Domain}, F<:NTuple{N, Any}, S<:NTuple{M, Real}}
     domains::D
     fs::F
-    continuity::NTuple{N, Int}
+    continuity::NTuple{M, Int}
     singularities::S
 
-    function Piecewise(domains::D, fs::F, continuity::NTuple{N, Int}, singularities::S) where {D, F, N, S}
-        if !(N + 1 == length(domains) == length(fs))
-            throw(ArgumentError(
-                "a Piecewise function must have as many domains as functions, " *
-                "and one less continuity point. " *
-                "Given: $(length(domains)) domains, $(length(fs)) functions, " *
-                "$N continuity points."))
-        end
+    function Piecewise(
+            domains::NTuple{N, Domain},
+            fs::NTuple{N, Any},
+            continuity::NTuple{M, Int},
+            singularities::NTuple{M, Real}) where {N, M}
+
+        N != M + 1 && throw(ArgumentError(
+                "a Piecewise function with N pieces must have N - 1 singularities. " *
+                "Given: $N pieces and $M singularities."
+        ))
 
-        return new{N, D, F, S}(domains, fs, continuity, singularities)
+        return new{N, M, typeof(domains), typeof(fs), typeof(singularities)}(domains, fs, continuity, singularities)
     end
 end
 
 function Piecewise(
-        domains::Vector{<:Domain},
-        fs,
-        continuity::Vector{Int} = fill(-1, length(domains) - 1))
+        domains::NTuple{N, Domain},
+        fs::NTuple{N, Any},
+        continuity::NTuple{M, Int} = ntuple(i -> -1, Vla(N-1))) where {N, M}
 
     if length(domains) != length(fs)
         throw(ArgumentError("the number of domains and the number of functions don't match"))
@@ -221,9 +223,11 @@ function Piecewise(
     return Piecewise(Tuple(domains), Tuple(fs), Tuple(continuity), Tuple(singularities))
 end
 
-function Piecewise(pairs::Vararg{Pair} ; continuity = fill(-1, length(pairs) - 1))
-    pairs = collect(pairs)
-    return Piecewise(first.(pairs), last.(pairs), continuity)
+function Piecewise(
+        pairs::Vararg{Pair, N} ;
+        continuity = ntuple(i -> -1, Val(N - 1))) where N
+
+    return Piecewise(first.(pairs), last.(pairs), Tuple(continuity))
 end
 
 domains(piecewise::Piecewise) = piecewise.domains
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index b98df54f..8e9c7658 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -4,6 +4,15 @@ using IntervalArithmetic: leftof
     @testset "Construction" begin
         @test_throws ErrorException Domain{:oopen, :closed}(0, 1)
         @test_throws ErrorException Domain{:open, :close}(0, 1)
+        @test_throws MethodError Piecewise(
+            (Domain{:closed, :closed}(1, 2), Domain{:open, :open}(2, 3)),
+            (sin, cos), (-1,), (2, 3))
+        @test_throws ArgumentError Piecewise(
+            (Domain{:closed, :closed}(1, 2), Domain{:open, :open}(2, 3)),
+            (sin, cos), (-1,-1), (2, 3))
+        @test_throws MethodError Piecewise(
+            (Domain{:closed, :closed}(1, 2), Domain{:open, :open}(2, 3)),
+            (sin, cos, log), (-1,), (2))
     end
 
     @testset "leftof" begin

From 8982c716df137d1347c7b13cac964470096316c2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Richard?= <kolaru@hotmail.com>
Date: Mon, 10 Feb 2025 11:57:02 +0100
Subject: [PATCH 20/20] Don't export set operations on Domain (for now)

---
 ext/IntervalArithmeticForwardDiffExt.jl | 31 +++++++++++-----------
 src/piecewise.jl                        | 34 ++++++++++++-------------
 test/interval_tests/piecewise.jl        | 26 +++++++++----------
 3 files changed, 46 insertions(+), 45 deletions(-)

diff --git a/ext/IntervalArithmeticForwardDiffExt.jl b/ext/IntervalArithmeticForwardDiffExt.jl
index a6dc5940..bcb22f6c 100644
--- a/ext/IntervalArithmeticForwardDiffExt.jl
+++ b/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
@@ -99,34 +100,34 @@ end
 
 function (piecewise::Piecewise)(dual::Dual{T, <:Interval}) where {T}
     X = value(dual)
-    set = Domain(X)
-    if !IntervalArithmetic.overlapdomain(set, piecewise) 
+    input_domain = Domain(X)
+    if !overlap_domain(input_domain, piecewise) 
         return Dual{T}(emptyinterval(X), emptyinterval(X) .* partials(dual))
     end
 
-    if !IntervalArithmetic.indomain(set, piecewise)
+    if !in_domain(input_domain, piecewise)
         dec = trv
-    elseif any(in(set), discontinuities(piecewise, 1))
+    elseif any(x -> in_domain(x, input_domain), discontinuities(piecewise, 1))
         dec = def 
     else
         dec = com
     end
 
-    dual_results = []
-    for (subdomain, f) in pieces(piecewise)
-        subset = intersect(set, subdomain)
-        isempty(subset) && continue
-        sub_X = interval(inf(subset), sup(subset), decoration(X))
+    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_results, f(sub_dual))
+        push!(dual_piece_outputs, f(sub_dual))
     end
 
-    results = value.(dual_results)
-    dec = min(dec, minimum(decoration.(results)))
-    primal = IntervalArithmetic.setdecoration(reduce(hull, results), dec)
+    piece_outputs = value.(dual_piece_outputs)
+    dec = min(dec, minimum(decoration.(piece_outputs)))
+    primal = IntervalArithmetic.setdecoration(reduce(hull, piece_outputs), dec)
 
-    dresults = partials.(dual_results)
-    partial = map(zip(dresults...)) do pp
+    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
diff --git a/src/piecewise.jl b/src/piecewise.jl
index e53470dc..6aa39b32 100644
--- a/src/piecewise.jl
+++ b/src/piecewise.jl
@@ -76,9 +76,9 @@ function leftof(d1::Domain, d2::Domain)
     return val1 < val2
 end
 
-Base.in(x::Real, domain::Domain) = rightof(x, lowerbound(domain)) && leftof(x, upperbound(domain))
+in_domain(x::Real, domain::Domain) = rightof(x, lowerbound(domain)) && leftof(x, upperbound(domain))
 
-function Base.intersect(d1::Domain, d2::Domain)
+function intersect_domain(d1::Domain, d2::Domain)
     left = max(lowerbound(d1), lowerbound(d2))
     right = min(upperbound(d1), upperbound(d2))
 
@@ -86,7 +86,7 @@ function Base.intersect(d1::Domain, d2::Domain)
     return Domain(left, right)
 end
 
-function Base.isempty(domain::Domain)
+function isempty_domain(domain::Domain)
     lo, lobound = lowerbound(domain)
     hi, hibound = upperbound(domain)
 
@@ -258,7 +258,7 @@ function Base.show(io::IO, ::MIME"text/plain", piecewise::Piecewise)
     end
 end
 
-function indomain(domain, piecewise)
+function in_domain(domain, piecewise)
     rightof(upperbound(domain), upperbound(domains(piecewise)[end])) && return false
 
     # This relies on the fact that domains are ordered
@@ -283,34 +283,34 @@ function indomain(domain, piecewise)
     return true
 end
 
-overlapdomain(domain, piecewise) = any(!isempty, intersect.(Ref(domain), domains(piecewise)))
+overlap_domain(domain, piecewise) = any(!isempty_domain, intersect_domain.(Ref(domain), domains(piecewise)))
 
 function (piecewise::Piecewise)(X::Interval{T}) where {T}
-    set = Domain(X)
-    !overlapdomain(set, piecewise) && return emptyinterval(T)
+    input_domain = Domain(X)
+    !overlap_domain(input_domain, piecewise) && return emptyinterval(T)
 
-    if !indomain(set, piecewise)
+    if !in_domain(input_domain, piecewise)
         dec = trv
-    elseif any(in(set), discontinuities(piecewise))
+    elseif any(x -> in_domain(x, input_domain), discontinuities(piecewise))
         dec = def 
     else
         dec = com
     end
 
-    results = Interval{T}[]
-    for (domain, f) in pieces(piecewise)
-        subset = intersect(set, domain)
-        isempty(subset) && continue
-        push!(results, f(interval(inf(subset), sup(subset), decoration(X))))
+    piece_outputs = Interval{T}[]
+    for (piece_domain, f) in pieces(piecewise)
+        piece_input = intersect_domain(input_domain, piece_domain)
+        isempty_domain(piece_input) && continue
+        push!(piece_outputs, f(interval(inf(piece_input), sup(piece_input), decoration(X))))
     end
 
-    dec = min(dec, minimum(decoration.(results)))
-    return IntervalArithmetic.setdecoration(reduce(hull, results), dec)
+    dec = min(dec, minimum(decoration.(piece_outputs)))
+    return IntervalArithmetic.setdecoration(reduce(hull, piece_outputs), dec)
 end
 
 function (piecewise::Piecewise)(x::Real)
     for (domain, f) in pieces(piecewise)
-        (x in domain) && return f(x)
+        (in_domain(x, domain)) && return f(x)
     end
     throw(DomainError(x, "piecewise function was called outside of its domain $(domain_string(piecewise))"))
 end
\ No newline at end of file
diff --git a/test/interval_tests/piecewise.jl b/test/interval_tests/piecewise.jl
index 8e9c7658..7cd18353 100644
--- a/test/interval_tests/piecewise.jl
+++ b/test/interval_tests/piecewise.jl
@@ -1,4 +1,4 @@
-using IntervalArithmetic: leftof
+using IntervalArithmetic: leftof, intersect_domain, isempty_domain
 
 @testset "Domain" begin
     @testset "Construction" begin
@@ -29,25 +29,25 @@ using IntervalArithmetic: leftof
         @test leftof(d2, d4)
     end
 
-    @testset "intersect" begin
+    @testset "intersect_domain" begin
         d1 = Domain{:closed, :open}(0, 10)
         d2 = Domain{:closed, :open}(2, 15)
         d3 = Domain{:open, :closed}(4, 7)
         d4 = Domain{:open, :closed}(-20, 3)
 
-        @test intersect(d1, d2) == Domain{:closed, :open}(2, 10)
-        @test intersect(d1, d3) == Domain{:open, :closed}(4, 7)
-        @test intersect(d1, d4) == Domain{:closed, :closed}(0, 3)
-        @test intersect(d2, d3) == Domain{:open, :closed}(4, 7)
-        @test intersect(d2, d4) == Domain{:closed, :closed}(2, 3)
-        @test intersect(d3, d4) == Domain()
+        @test intersect_domain(d1, d2) == Domain{:closed, :open}(2, 10)
+        @test intersect_domain(d1, d3) == Domain{:open, :closed}(4, 7)
+        @test intersect_domain(d1, d4) == Domain{:closed, :closed}(0, 3)
+        @test intersect_domain(d2, d3) == Domain{:open, :closed}(4, 7)
+        @test intersect_domain(d2, d4) == Domain{:closed, :closed}(2, 3)
+        @test intersect_domain(d3, d4) == Domain()
     end
 
-    @testset "isempty" begin
-        @test isempty(Domain{:open, :open}(1, 1))
-        @test !isempty(Domain{:closed, :closed}(1, 1))
-        @test !isempty(Domain{:open, :open}(1, 2))
-        @test isempty(Domain())
+    @testset "isempty_domain" begin
+        @test isempty_domain(Domain{:open, :open}(1, 1))
+        @test !isempty_domain(Domain{:closed, :closed}(1, 1))
+        @test !isempty_domain(Domain{:open, :open}(1, 2))
+        @test isempty_domain(Domain())
     end
 end