diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index dcec41d6..13a4973f 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -79,7 +79,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1'
+          version: '1.6'
       - run: |
           julia --project=docs -e '
             using Pkg
diff --git a/Project.toml b/Project.toml
index 310f0476..17b9f7a7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "Intervals"
 uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 license = "MIT"
 authors = ["Invenia Technical Computing"]
-version = "1.7.0"
+version = "1.8.0"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
diff --git a/docs/src/index.md b/docs/src/index.md
index 4fc13736..9c86ba9d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -268,6 +268,7 @@ In the plot, inclusive boundaries are marked with a vertical bar, whereas exclus
 ```@docs
 Interval
 AnchoredInterval
+IntervalSet
 HourEnding
 HourBeginning
 HE
@@ -291,4 +292,5 @@ Base.parse(::Type{Interval{T}}, ::AbstractString) where T
 union
 union!
 superset
+Intervals.find_intersections
 ```
diff --git a/src/Intervals.jl b/src/Intervals.jl
index fe00041a..7d3a939a 100644
--- a/src/Intervals.jl
+++ b/src/Intervals.jl
@@ -42,6 +42,7 @@ export Bound,
        Unbounded,
        AbstractInterval,
        Interval,
+       IntervalSet,
        AnchoredInterval,
        HourEnding,
        HourBeginning,
diff --git a/src/deprecated.jl b/src/deprecated.jl
index e91ae7c4..e6c370b5 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -195,4 +195,8 @@ function HB(anchor, inc::Inclusivity)
     return HourBeginning{L,R}(floor(anchor, Hour))
 end
 
+@deprecate union(intervals::AbstractVector{<:AbstractInterval}) convert(Vector, union(IntervalSet(intervals)))
+@deprecate union!(intervals::AbstractVector{<:AbstractInterval}) convert(Vector, union!(IntervalSet(intervals)))
+@deprecate superset(intervals::AbstractVector{<:AbstractInterval}) superset(IntervalSet(intervals))
+
 # END Intervals 1.X.Y deprecations
diff --git a/src/interval_sets.jl b/src/interval_sets.jl
index e0760a7b..16211696 100644
--- a/src/interval_sets.jl
+++ b/src/interval_sets.jl
@@ -1,7 +1,98 @@
 
 ###### Set-related Helpers #####
+"""
+    IntervalSet{T<:AbstractInterval}
+
+A set of points represented by a sequence of intervals. Set operations over interval sets
+return a new IntervalSet, with the fewest number of intervals possible. Unbounded intervals
+are not supported. The individual intervals in the set can be accessed using the iteration
+API or by passing the set to `Array`.
+
+see also: https://en.wikipedia.org/wiki/Interval_arithmetic#Interval_operators
+
+## Examples
+
+```jldoctest; setup = :(using Intervals)
+julia> union(IntervalSet(1..5), IntervalSet(3..8))
+1-interval IntervalSet{Interval{Int64, Closed, Closed}}:
+[1 .. 8]
+
+julia> intersect(IntervalSet(1..5), IntervalSet(3..8))
+1-interval IntervalSet{Interval{Int64, Closed, Closed}}:
+[3 .. 5]
+
+julia> symdiff(IntervalSet(1..5), IntervalSet(3..8))
+2-interval IntervalSet{Interval{Int64, L, R} where {L<:Bound, R<:Bound}}:
+[1 .. 3)
+(5 .. 8]
+
+julia> union(IntervalSet([1..2, 2..5]), IntervalSet(6..7))
+2-interval IntervalSet{Interval{Int64, Closed, Closed}}:
+[1 .. 5]
+[6 .. 7]
+
+julia> union(IntervalSet([1..5, 8..10]), IntervalSet([4..9, 12..14]))
+2-interval IntervalSet{Interval{Int64, Closed, Closed}}:
+[1 .. 10]
+[12 .. 14]
+
+julia> intersect(IntervalSet([1..5, 8..10]), IntervalSet([4..9, 12..14]))
+2-interval IntervalSet{Interval{Int64, Closed, Closed}}:
+[4 .. 5]
+[8 .. 9]
+
+julia> setdiff(IntervalSet([1..5, 8..10]), IntervalSet([4..9, 12..14]))
+2-interval IntervalSet{Interval{Int64, L, R} where {L<:Bound, R<:Bound}}:
+[1 .. 4)
+(9 .. 10]
+```
+"""
+struct IntervalSet{T <: AbstractInterval}
+    items::Vector{T}
+end
+
+IntervalSet(interval::T) where T <: AbstractInterval = IntervalSet{T}([interval])
+IntervalSet(interval::IntervalSet) = interval
+IntervalSet(itr) = IntervalSet{eltype(itr)}(collect(itr))
+IntervalSet() = IntervalSet(AbstractInterval[])
+
+Base.copy(intervals::IntervalSet{T}) where {T} = IntervalSet{T}(copy(intervals.items))
+Base.eltype(::IntervalSet{T}) where T = T
+Base.isempty(intervals::IntervalSet) = isempty(intervals.items) || all(isempty, intervals.items)
+Base.:(==)(a::IntervalSet, b::IntervalSet) = issetequal(a, b)
+Base.isequal(a::IntervalSet, b::IntervalSet) = isequal(a.items, b.items)
+Base.convert(::Type{T}, intervals::IntervalSet) where T <: AbstractArray = convert(T, intervals.items)
+function Base.show(io::Base.AbstractPipe, ::MIME"text/plain", x::IntervalSet)
+    intervals = union(x)
+    n = length(intervals.items)
+    iocompact = IOContext(io, :compact => true)
+    print(io, "$n-interval ")
+    show(io, MIME"text/plain"(), typeof(x))
+    println(io, ":")
+    nrows = displaysize(io)[1]
+    half = fld(nrows, 2) - 2
+    if nrows ≥ n && half > 1
+        for interval in intervals.items[1:(end-1)]
+            show(iocompact, MIME"text/plain"(), interval)
+            println(io, "")
+        end
+        isempty(intervals) || show(iocompact, MIME"text/plain"(), intervals.items[end])
+    else
+        for interval in intervals.items[1:half]
+            show(iocompact, MIME"text/plain"(), interval)
+            println(io, "")
+        end
+        println(io, "⋮")
+        for interval in intervals.items[(end-half+1):end-1]
+            show(iocompact, MIME"text/plain"(), interval)
+            println(io, "")
+        end
+        show(iocompact, MIME"text/plain"(), intervals.items[end])
+    end
+end
 
-const IntervalSet = AbstractVector{<:AbstractInterval}
+# currently (to avoid breaking changes) new methods for `Base`
+# accept `IntervalSet` objects and Interval singletons.
 const AbstractIntervals = Union{AbstractInterval, IntervalSet}
 
 # During merge operations used to compute unions, intersections etc...,
@@ -39,8 +130,10 @@ function endpoint_tracking(
 )
     return TrackEachEndpoint()
 end
-endpoint_tracking(a::AbstractVector, b::AbstractVector) = endpoint_tracking(eltype(a), eltype(b))
+
+endpoint_tracking(a::IntervalSet, b::IntervalSet) = endpoint_tracking(eltype(a), eltype(b))
 endpoint_tracking(a::AbstractInterval, b::AbstractInterval) = endpoint_tracking(typeof(a), typeof(b))
+endpoint_tracking(a::AbstractVector, b::AbstractVector) = endpoint_tracking(eltype(a), eltype(b))
 
 # track: run a thunk, but only if we are tracking endpoints dynamically
 track(fn::Function, ::TrackEachEndpoint, args...) = fn(args...)
@@ -58,12 +151,21 @@ interval_type(::TrackRightOpen{T}) where T = Interval{T, Closed, Open}
 # all vectors of sets be represented by their endpoints. The functions unbunch
 # and bunch convert between an interval and an endpoint representation
 
-function unbunch(interval::AbstractInterval, tracking::EndpointTracking; lt=isless) 
+function unbunch(interval::AbstractInterval, tracking::EndpointTracking; lt=isless)
     return endpoint_type(tracking)[LeftEndpoint(interval), RightEndpoint(interval)]
 end
+function unbunch(intervals::IntervalSet, tracking::EndpointTracking; kwargs...)
+    return unbunch(convert(Vector, intervals), tracking; kwargs...)
+end
 unbunch_by_fn(_) = identity
-function unbunch(intervals::Union{AbstractIntervals, Base.Iterators.Enumerate{<:AbstractIntervals}}, 
-                 tracking::EndpointTracking; lt=isless)
+function unbunch(
+    intervals::Union{
+        AbstractVector{<:AbstractInterval},
+        Base.Iterators.Enumerate{<:Union{AbstractIntervals, AbstractVector{<:AbstractInterval}}}
+    },
+    tracking::EndpointTracking;
+    lt=isless,
+)
     by = unbunch_by_fn(intervals)
     filtered = Iterators.filter(!isempty ∘ by, intervals)
     isempty(filtered) && return endpoint_type(tracking)[]
@@ -72,12 +174,13 @@ function unbunch(intervals::Union{AbstractIntervals, Base.Iterators.Enumerate{<:
 end
 # support for `unbunch(enumerate(vcat(x)))` (transforming [(i, interval)] -> [(i, endpoint), (i,endpoint)])
 unbunch_by_fn(::Base.Iterators.Enumerate) = last
-function unbunch((i, interval)::Tuple, tracking; lt=isless) 
+function unbunch((i, interval)::Tuple, tracking; lt=isless)
     eltype = Tuple{Int, endpoint_type(tracking)}
     return eltype[(i, LeftEndpoint(interval)), (i, RightEndpoint(interval))]
 end
 
-function unbunch(a::AbstractIntervals, b::AbstractIntervals; kwargs...)
+function unbunch(a::Union{AbstractVector{<:AbstractInterval}, AbstractIntervals},
+                 b::Union{AbstractVector{<:AbstractInterval}, AbstractIntervals}; kwargs...)
     tracking = endpoint_tracking(a, b)
     a_ = unbunch(a, tracking; kwargs...)
     b_ = unbunch(b, tracking; kwargs...)
@@ -87,16 +190,17 @@ end
 # represent a sequence of endpoints as a sequence of one or more intervals
 function bunch(endpoints, tracking)
     @assert iseven(length(endpoints))
-    isempty(endpoints) && return interval_type(tracking)[]
-    return map(Iterators.partition(endpoints, 2)) do pair
+    isempty(endpoints) && return IntervalSet(interval_type(tracking)[])
+    res = map(Iterators.partition(endpoints, 2)) do pair
         return Interval(pair..., tracking)
     end
+    return IntervalSet(res)
 end
 Interval(a::Endpoint, b::Endpoint, ::TrackEachEndpoint) = Interval(a, b)
 Interval(a::Endpoint, b::Endpoint, ::TrackLeftOpen{T}) where T = Interval{T,Open,Closed}(a.endpoint, b.endpoint)
 Interval(a::Endpoint, b::Endpoint, ::TrackRightOpen{T}) where T = Interval{T,Closed,Open}(a.endpoint, b.endpoint)
 
-# the sentinel endpoint reduces the number of edgecases 
+# the sentinel endpoint reduces the number of edgecases
 # we have to deal with when comparing endpoints during a merge
 # NOTE: it's tempting to replace this with an unbounded endpoint
 # but if we ever want to support unbounded endpoints in mergesets
@@ -108,7 +212,7 @@ function first_endpoint(x)
     return eltype(x) <: Tuple ? last(first(x)) : first(x)
 end
 function last_endpoint(x)
-    isempty(x) && return SentinelEndpoint() 
+    isempty(x) && return SentinelEndpoint()
     # if the endpoints are enumerated, eltype will be a tuple
     return eltype(x) <: Tuple ? last(last(x)) : last(x)
 end
@@ -139,7 +243,7 @@ isleft(::RightEndpoint) = false
 # will remain unchanged moving left to right along the real-number line until we encounter a
 # new endpoint.
 #
-# For each endpoint, we determine two things: 
+# For each endpoint, we determine two things:
 #   1. whether subsequent points should be included in the merge operation or not (based on
 #        its membership in both `x` and `y`) by using `op`
 #   2. whether the next step will define a left (start including) or right endpoint (stop
@@ -264,35 +368,37 @@ right_endpoint(t, ::TrackRightOpen{T}) where T = RightEndpoint{T,Open}(endpoint(
 
 # There is power in a union.
 """
-    union(intervals::AbstractVector{<:AbstractInterval})
+    union(intervals::IntervalSets)
 
-Flattens a vector of overlapping intervals into a new, smaller vector containing only
-non-overlapping intervals.
+Flattens any overlapping intervals within the `IntervalSet` into a new, smaller set
+containing only non-overlapping intervals.
 """
-function Base.union(intervals::AbstractVector{<:AbstractInterval})
-    return union!(convert(Vector{AbstractInterval}, intervals))
-end
-
-# allow a concretely typed array for `Interval` objects (as opposed to e.g. anchored intervals
-# which may change type during the union process)
-function Base.union(intervals::AbstractVector{T}) where T <: Interval
-    return union!(copy(intervals))
+Base.union(intervals::IntervalSet{<:Interval}) = union!(copy(intervals))
+
+# In the case where we're dealing with a non-concrete interval type like AnchoredIntervals then simply
+# allocate a AbstractInterval vector
+function Base.union(intervals::IntervalSet{<:AbstractInterval})
+    T = AbstractInterval
+    dest = Vector{T}(undef, length(intervals.items))
+    copyto!(dest, intervals.items)
+    return union!(IntervalSet{T}(dest))
 end
 
 """
-    union!(intervals::AbstractVector{<:AbstractInterval})
+    union!(intervals::IntervalSet)
 
 Flattens a vector of overlapping intervals in-place to be a smaller vector containing only
 non-overlapping intervals.
 """
-function Base.union!(intervals::AbstractVector{<:AbstractInterval})
-    sort!(intervals)
+function Base.union!(intervals::IntervalSet)
+    items = intervals.items
+    sort!(items)
 
     i = 2
-    n = length(intervals)
+    n = length(items)
     while i <= n
-        prev = intervals[i - 1]
-        curr = intervals[i]
+        prev = items[i - 1]
+        curr = items[i]
 
         # If the current and previous intervals don't meet then move along
         if !overlaps(prev, curr) && !contiguous(prev, curr)
@@ -301,8 +407,8 @@ function Base.union!(intervals::AbstractVector{<:AbstractInterval})
         # If the two intervals meet then we absorb the current interval into
         # the previous one.
         else
-            intervals[i - 1] = merge(prev, curr)
-            deleteat!(intervals, i)
+            items[i - 1] = merge(prev, curr)
+            deleteat!(items, i)
             n -= 1
         end
     end
@@ -311,13 +417,13 @@ function Base.union!(intervals::AbstractVector{<:AbstractInterval})
 end
 
 """
-    superset(intervals::AbstractArray{<:AbstractInterval}) -> Interval
+    superset(intervals::IntervalSet) -> Interval
 
 Create the smallest single interval which encompasses all of the provided intervals.
 """
-function superset(intervals::AbstractArray{<:AbstractInterval})
-    left = minimum(LeftEndpoint.(intervals))
-    right = maximum(RightEndpoint.(intervals))
+function superset(intervals::IntervalSet)
+    left = minimum(LeftEndpoint.(intervals.items))
+    right = maximum(RightEndpoint.(intervals.items))
 
     return Interval(left, right)
 end
@@ -328,18 +434,31 @@ Base.intersect(x::IntervalSet, y::IntervalSet) = mergesets((inx, iny) -> inx &&
 Base.union(x::IntervalSet, y::IntervalSet) = mergesets((inx, iny) -> inx || iny, x, y)
 Base.setdiff(x::IntervalSet, y::IntervalSet) = mergesets((inx, iny) -> inx && !iny, x, y)
 Base.symdiff(x::IntervalSet, y::IntervalSet) = mergesets((inx, iny) -> inx ⊻ iny, x, y)
-Base.issubset(x::AbstractIntervals, y::AbstractIntervals) = isempty(setdiff(x, y))
 Base.isdisjoint(x::AbstractIntervals, y::AbstractIntervals) = isempty(intersect(x, y))
 
+Base.issubset(x, y::IntervalSet) = x in y
+Base.issubset(x::AbstractInterval, y::IntervalSet) = any(Base.Fix1(issubset, x), y.items)
+Base.issubset(x::IntervalSet, y::AbstractInterval) = all(Base.Fix2(issubset, y), x.items)
+Base.issubset(x::IntervalSet, y::IntervalSet) = isempty(setdiff(x, y))
+
+# Add methods where just 1 argument is an Interval
+for f in (:intersect, :union, :setdiff, :symdiff)
+    @eval Base.$f(x::AbstractInterval, y::IntervalSet) = $f(IntervalSet([x]), y)
+    @eval Base.$f(x::IntervalSet, y::AbstractInterval) = $f(x, IntervalSet([y]))
+end
+
 function Base.issetequal(x::AbstractIntervals, y::AbstractIntervals)
-    x, y, tracking = unbunch(union(vcat(x)), union(vcat(y)))
-    return x == y || all(isempty, bunch(x, tracking)) && all(isempty, bunch(y, tracking))
+    x, y, tracking = unbunch(union(IntervalSet(x)), union(IntervalSet(y)))
+    return x == y || isempty(bunch(x, tracking)) && isempty(bunch(y, tracking))
 end
 
+# when `x` is number-like object (Number, Date, Time, etc...):
+Base.in(x, y::IntervalSet) = any(Base.Fix1(in, x), y.items)
+
 # order edges so that closed boundaries are on the outside: e.g. [( )]
 intersection_order(x::Endpoint) = isleft(x) ? !isclosed(x) : isclosed(x)
 intersection_isless_fn(::TrackStatically) = isless
-function intersection_isless_fn(::TrackEachEndpoint) 
+function intersection_isless_fn(::TrackEachEndpoint)
     function (x,y)
         if isequal(x, y)
             return isless(intersection_order(x), intersection_order(y))
@@ -351,22 +470,22 @@ end
 
 """
     find_intersections(
-        x::Union{AbstractInterval, AbstractVector{<:AbstractInterval}},
-        y::Union{AbstractInterval, AbstractVector{<:AbstractInterval}},
+        x::AbstractVector{<:AbstractInterval},
+        y::AbstractVector{<:AbstractInterval}
     )
 
 Returns a `Vector{Vector{Int}}` where the value at index `i` gives the indices to all
 intervals in `y` that intersect with `x[i]`.
 """
-function find_intersections(x_::AbstractIntervals, y_::AbstractIntervals)
-    xa, ya = vcat(x_), vcat(y_)
-    tracking = endpoint_tracking(xa, ya)
+find_intersections(x, y) = find_intersections(vcat(x), vcat(y))
+function find_intersections(x::AbstractVector{<:AbstractInterval}, y::AbstractVector{<:AbstractInterval})
+    tracking = endpoint_tracking(x, y)
     lt = intersection_isless_fn(tracking)
-    x = unbunch(enumerate(xa), tracking; lt)
-    y = unbunch(enumerate(ya), tracking; lt)
-    result = [Vector{Int}() for _ in 1:length(xa)]
+    x_endpoints = unbunch(enumerate(x), tracking; lt)
+    y_endpoints = unbunch(enumerate(y), tracking; lt)
+    result = [Vector{Int}() for _ in 1:length(x)]
 
-    return find_intersections_helper!(result, x, y, lt)
+    return find_intersections_helper!(result, x_endpoints, y_endpoints, lt)
 end
 
 function find_intersections_helper!(result, x, y, lt)
@@ -402,4 +521,3 @@ function find_intersections_helper!(result, x, y, lt)
 
     return unique!.(result)
 end
-
diff --git a/test/comparisons.jl b/test/comparisons.jl
index 1ddcd7ca..2c8b131d 100644
--- a/test/comparisons.jl
+++ b/test/comparisons.jl
@@ -89,16 +89,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([earlier, later]) == [earlier, later]
+            @test union(IntervalSet([earlier, later])) == IntervalSet([earlier, later])
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(earlier) && isbounded(later)
                 @test union([earlier], [later]) == [earlier, later]
+                @test union(IntervalSet(earlier), IntervalSet(later)) == IntervalSet([earlier, later])
+
                 @test intersect([earlier], [later]) == []
+                @test intersect(IntervalSet(earlier), IntervalSet(later)) == IntervalSet()
 
                 @test setdiff([earlier], [later]) == expected_xor[1:1]
+                @test setdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor[1:1])
+
                 @test setdiff([later], [earlier]) == expected_xor[2:2]
+                @test setdiff(IntervalSet(later), IntervalSet(earlier)) == IntervalSet(expected_xor[2:2])
 
                 @test symdiff([earlier], [later]) == expected_xor
+                @test symdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -168,16 +176,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([earlier, later]) == [earlier, later]
+            @test union(IntervalSet([earlier, later])) == IntervalSet([earlier, later])
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(earlier) && isbounded(later)
                 @test union([earlier], [later]) == [earlier, later]
+                @test union(IntervalSet(earlier), IntervalSet(later)) == IntervalSet([earlier, later])
+
                 @test intersect([earlier], [later]) == []
+                @test intersect(IntervalSet(earlier), IntervalSet(later)) == IntervalSet()
 
                 @test setdiff([earlier], [later]) == expected_xor[1:1]
+                @test setdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor[1:1])
+
                 @test setdiff([later], [earlier]) == expected_xor[2:2]
+                @test setdiff(IntervalSet(later), IntervalSet(earlier)) == IntervalSet(expected_xor[2:2])
 
                 @test symdiff([earlier], [later]) == expected_xor
+                @test symdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -247,20 +263,30 @@ end
 
             # Using a vector of intervals as sets
             @test union([earlier, later]) == [expected_superset]
+            @test union(IntervalSet([earlier, later])) == IntervalSet([expected_superset])
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(earlier) && isbounded(later)
-                @test union([earlier], [later]) == [expected_superset]
+                @test union([earlier], [later]) == expected_xor != [expected_superset]
+                @test union(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_superset)
+
                 @test intersect([earlier], [later]) == []
+                @test intersect(IntervalSet(earlier), IntervalSet(later)) == IntervalSet()
 
                 @test setdiff([earlier], [later]) == expected_xor[1:1]
+                @test setdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor[1:1])
+
                 @test setdiff([later], [earlier]) == expected_xor[2:2]
+                @test setdiff(IntervalSet(later), IntervalSet(earlier)) == IntervalSet(expected_xor[2:2])
 
-                @test symdiff([earlier], [later]) == union(expected_xor)
+                # TODO: Sometimes expected_xor would get mutated in this call
+                @test symdiff([earlier], [later]) == expected_xor != union(expected_xor)
+                @test symdiff(IntervalSet(earlier), IntervalSet(later)) == union(IntervalSet(expected_xor))
             end
         end
     end
 
+
     # Compare two intervals which "touch" and the earlier interval includes that point:
     # Visualization of the finite case:
     #
@@ -326,16 +352,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([earlier, later]) == [expected_superset]
+            @test union(IntervalSet([earlier, later])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(earlier) && isbounded(later)
-                @test union([earlier], [later]) == [expected_superset]
+                @test union([earlier], [later]) == expected_xor != [expected_superset]
+                @test union(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_superset)
+
                 @test intersect([earlier], [later]) == []
+                @test intersect(IntervalSet(earlier), IntervalSet(later)) == IntervalSet()
 
                 @test setdiff([earlier], [later]) == expected_xor[1:1]
+                @test setdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor[1:1])
+
                 @test setdiff([later], [earlier]) == expected_xor[2:2]
+                @test setdiff(IntervalSet(later), IntervalSet(earlier)) == IntervalSet(expected_xor[2:2])
 
-                @test symdiff([earlier], [later]) == union(expected_xor)
+                @test symdiff([earlier], [later]) == expected_xor != union(expected_xor)
+                @test symdiff(IntervalSet(earlier), IntervalSet(later)) == union(IntervalSet(expected_xor))
             end
         end
     end
@@ -410,16 +444,25 @@ end
 
             # Using a vector of intervals as sets
             @test union([earlier, later]) == [expected_superset]
+            @test union(IntervalSet([earlier, later])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(earlier) && isbounded(later)
-                @test union([earlier], [later]) == [expected_superset]
-                @test intersect([earlier], [later]) == [expected_overlap]
+                # NOTE: expected_xor may have different bounds than [earlier, later]
+                @test union([earlier], [later]) == [earlier, later] != [expected_superset]
+                @test union(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_superset)
 
-                @test setdiff([earlier], [later]) == expected_xor[1:1]
-                @test setdiff([later], [earlier]) == expected_xor[2:2]
+                @test intersect([earlier], [later]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_overlap)
 
-                @test symdiff([earlier], [later]) == expected_xor
+                @test setdiff([earlier], [later]) == [earlier] != expected_xor[1:1]
+                @test setdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor[1:1])
+
+                @test setdiff([later], [earlier]) == [later] != expected_xor[2:2]
+                @test setdiff(IntervalSet(later), IntervalSet(earlier)) == IntervalSet(expected_xor[2:2])
+
+                @test symdiff([earlier], [later]) == [earlier, later]
+                @test symdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -494,16 +537,25 @@ end
 
             # Using a vector of intervals as sets
             @test union([earlier, later]) == [expected_superset]
+            @test union(IntervalSet([earlier, later])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(earlier) && isbounded(later)
-                @test union([earlier], [later]) == [expected_superset]
-                @test intersect([earlier], [later]) == [expected_overlap]
+                # NOTE: expected_xor may have different bounds than [earlier, later]
+                @test union([earlier], [later]) == [earlier, later] != [expected_superset]
+                @test union(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_superset)
 
-                @test setdiff([earlier], [later]) == expected_xor[1:1]
-                @test setdiff([later], [earlier]) == expected_xor[2:2]
+                @test intersect([earlier], [later]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_overlap)
 
-                @test symdiff([earlier], [later]) == expected_xor
+                @test setdiff([earlier], [later]) == [earlier] != expected_xor[1:1]
+                @test setdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor[1:1])
+
+                @test setdiff([later], [earlier]) == [later] != expected_xor[2:2]
+                @test setdiff(IntervalSet(later), IntervalSet(earlier)) == IntervalSet(expected_xor[2:2])
+
+                @test symdiff([earlier], [later]) == [earlier, later] != expected_xor
+                @test symdiff(IntervalSet(earlier), IntervalSet(later)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -560,16 +612,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
                 @test union([a], [b]) == [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
+
                 @test intersect([a], [b]) == [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
 
                 @test setdiff([a], [b]) == []
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
+
                 @test setdiff([b], [a]) == []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
 
                 @test symdiff([a], [b]) == []
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
             end
         end
     end
@@ -627,18 +687,27 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             # TODO: will have to think carefully about the `expected_` variables
             # when we allow for unbounded values
             if isbounded(a) && isbounded(b)
-                @test union([a], [b]) == [expected_superset]
-                @test intersect([a], [b]) == [expected_overlap]
+                # NOTE: expected_xor may have different bounds than [a, b]
+                @test union([a], [b]) == [a, b] != [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
 
-                @test setdiff([a], [b]) == expected_xor[1:1]
-                @test setdiff([b], [a]) == []
+                @test intersect([a], [b]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
+
+                @test setdiff([a], [b]) == [a] != expected_xor[1:1]
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor[1:1])
 
-                @test symdiff([a], [b]) == expected_xor
+                @test setdiff([b], [a]) == [b] != []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
+
+                @test symdiff([a], [b]) == [a, b] != expected_xor
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -696,18 +765,27 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             # TODO: will have to think carefully about the `expected_` variables
             # when we allow for unbounded values
             if isbounded(a) && isbounded(b)
-                @test union([a], [b]) == [expected_superset]
-                @test intersect([a], [b]) == [expected_overlap]
+                # NOTE: expected_xor may have different bounds than [a, b]
+                @test union([a], [b]) == [a, b] != [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
 
-                @test setdiff([a], [b]) == expected_xor[1:1]
-                @test setdiff([b], [a]) == []
+                @test intersect([a], [b]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
+
+                @test setdiff([a], [b]) == [a] != expected_xor[1:1]
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor[1:1])
 
-                @test symdiff([a], [b]) == expected_xor
+                @test setdiff([b], [a]) == [b] != []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
+
+                @test symdiff([a], [b]) == [a, b] != expected_xor
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -769,18 +847,27 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             # TODO: will have to think carefully about the `expected_` variables
             # when we allow for unbounded values
             if isbounded(a) && isbounded(b)
-                @test union([a], [b]) == [expected_superset]
-                @test intersect([a], [b]) == [expected_overlap]
+                # NOTE: expected_xor may have different bounds than [a, b]
+                @test union([a], [b]) == [a, b] != [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
 
-                @test setdiff([a], [b]) == expected_xor
-                @test setdiff([b], [a]) == []
+                @test intersect([a], [b]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
 
-                @test symdiff([a], [b]) == expected_xor
+                @test setdiff([a], [b]) == [a] != expected_xor
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor)
+
+                @test setdiff([b], [a]) == [b] != []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
+
+                @test symdiff([a], [b]) == [a, b] != expected_xor
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -840,16 +927,25 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
-                @test union([a], [b]) == [expected_superset]
-                @test intersect([a], [b]) == [expected_overlap]
+                # NOTE: expected_xor may have different bounds than [a, b]
+                @test union([a], [b]) == [a, b] != [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
 
-                @test setdiff([a], [b]) == []
-                @test setdiff([b], [a]) == expected_xor[1:1]
+                @test intersect([a], [b]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
+
+                @test setdiff([a], [b]) == [a] != []
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
 
-                @test symdiff([a], [b]) == expected_xor
+                @test setdiff([b], [a]) == [b] != expected_xor[1:1]
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet(expected_xor[1:1])
+
+                @test symdiff([a], [b]) == [a, b] != expected_xor
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -909,16 +1005,27 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
-            @test intersect([a], [b]) == [expected_overlap]
+            @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
+            @test intersect([a], [b]) != [expected_overlap]
+            @test_broken intersect(IntervalSet([a, b])) == IntervalSet(expected_overlap) # Internal type issue
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
-                @test union([a], [b]) == [expected_superset]
+                # NOTE: expected_xor may have different bounds than [a, b]
+                @test union([a], [b]) == [a, b] != [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
 
-                @test setdiff([a], [b]) == []
-                @test setdiff([b], [a]) == expected_xor[1:1]
+                @test intersect([a], [b]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
+
+                @test setdiff([a], [b]) == [a] != []
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
+
+                @test setdiff([b], [a]) == [b] != expected_xor[1:1]
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet(expected_xor[1:1])
 
-                @test symdiff([a], [b]) == expected_xor
+                @test symdiff([a], [b]) == [a, b] != expected_xor
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor)
             end
         end
     end
@@ -975,16 +1082,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
                 @test union([a], [b]) == [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
+
                 @test intersect([a], [b]) == [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
 
                 @test setdiff([a], [b]) == []
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
+
                 @test setdiff([b], [a]) == []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
 
                 @test symdiff([a], [b]) == []
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
             end
         end
     end
@@ -1042,16 +1157,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
                 @test union([a], [b]) == [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
+
                 @test intersect([a], [b]) == [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
+
+                @test setdiff([a], [b]) == [a] != expected_xor[1:1]
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor[1:1])
 
-                @test setdiff([a], [b]) == expected_xor[1:1]
                 @test setdiff([b], [a]) == []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
 
                 @test symdiff([a], [b]) == []
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
             end
         end
     end
@@ -1109,16 +1232,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
                 @test union([a], [b]) == [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
+
                 @test intersect([a], [b]) == [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
+
+                @test setdiff([a], [b]) == [a] != expected_xor[1:1]
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor[1:1])
 
-                @test setdiff([a], [b]) == expected_xor[1:1]
                 @test setdiff([b], [a]) == []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
 
                 @test symdiff([a], [b]) == []
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
             end
         end
     end
@@ -1174,16 +1305,24 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
                 @test union([a], [b]) == [expected_superset]
+                @test union(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_superset)
+
                 @test intersect([a], [b]) == [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
 
                 @test setdiff([a], [b]) == []
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
+
                 @test setdiff([b], [a]) == []
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
 
                 @test symdiff([a], [b]) == []
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
             end
         end
     end
@@ -1231,14 +1370,20 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
-            @test union([a], [b]) == [expected_superset]
+            @test union([a], [b]) == [a, b] != [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
+
+            @test intersect([a], [b]) == [] != [expected_overlap]
+            @test intersect(IntervalSet([a, b])) == IntervalSet(expected_overlap)
 
-            @test intersect([a], [b]) == [expected_overlap]
+            @test setdiff([a], [b]) == [a] != []
+            @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
 
-            @test setdiff([a], [b]) == []
-            @test setdiff([b], [a]) == []
+            @test setdiff([b], [a]) == [b] != []
+            @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet()
 
-            @test symdiff([a], [b]) == []
+            @test symdiff([a], [b]) == [a, b] != []
+            @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
         end
     end
 
@@ -1319,16 +1464,25 @@ end
 
             # Using a vector of intervals as sets
             @test union([a, b]) == [expected_superset]
+            @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
             # TODO: These functions should be compatible with unbounded intervals
             if isbounded(a) && isbounded(b)
-                @test union([a], [b]) == [expected_superset]
-                @test intersect([a], [b]) == [expected_overlap]
+                # NOTE: expected_xor may have different bounds than [a, b]
+                @test union([a], [b]) == [a, b] != [expected_superset]
+                @test union(IntervalSet([a, b])) == IntervalSet(expected_superset)
 
-                @test setdiff([a], [b]) == []
-                @test setdiff([b], [a]) == expected_xor[1:2]
+                @test intersect([a], [b]) == [] != [expected_overlap]
+                @test intersect(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_overlap)
+
+                @test setdiff([a], [b]) == [a] != []
+                @test setdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet()
+
+                @test setdiff([b], [a]) == [b] != expected_xor[1:2]
+                @test setdiff(IntervalSet(b), IntervalSet(a)) == IntervalSet(expected_xor[1:2])
 
-                @test symdiff([a], [b]) == expected_xor
+                @test symdiff([a], [b]) == [a, b] != expected_xor
+                @test symdiff(IntervalSet(a), IntervalSet(b)) == IntervalSet(expected_xor)
             end
         end
     end
diff --git a/test/interval.jl b/test/interval.jl
index 197a06a9..b85814d6 100644
--- a/test/interval.jl
+++ b/test/interval.jl
@@ -765,6 +765,7 @@
             Interval{Open, Open}(-10, -1),
             Interval{Open, Open}(13, 20),
         ]
+
         @test union!(intervals) == expected
         @test intervals == expected
 
diff --git a/test/sets.jl b/test/sets.jl
index 839a38c4..e28e24fc 100644
--- a/test/sets.jl
+++ b/test/sets.jl
@@ -28,12 +28,29 @@ end
     area(x::Interval) = last(x) - first(x)
     # note: `mapreduce` fails here for empty vectors
     area(x::AbstractVector{<:AbstractInterval{T}}) where T = mapreduce(area, +, x, init=zero(T))
+    area(x::IntervalSet) = area(x.items)
     area(x) = isempty(x) ? 0 : error("Undefined area for object of type $(typeof(x))")
     myunion(x::Interval) = x
     myunion(x::AbstractVector{<:Interval}) = union(x)
+    myunion(x::IntervalSet) = union(x)
 
     rand_bound_type(rng) = rand(rng, (Closed, Open))
 
+    # verify case where we interpret array as a set of elements (rather than an
+    # interval-bound point set)
+    @test intersect([1..2, 2..3, 3..4, 4..5], [2..3, 3..4]) == [2..3, 3..4]
+
+    # verify that elements are in / subsets of interval sets
+    @test 2 ∈ IntervalSet([1..3, 5..10])
+    @test 0 ∉ IntervalSet([1..3, 5..10])
+    @test 4 ∉ IntervalSet([1..3, 5..10])
+    @test 11 ∉ IntervalSet([1..3, 5..10])
+    @test issubset(2, IntervalSet([1..3, 5..10]))
+    @test !issubset(0, IntervalSet([1..3, 5..10]))
+    @test !issubset(4, IntervalSet([1..3, 5..10]))
+    @test !issubset(11, IntervalSet([1..3, 5..10]))
+    @test issubset(2, IntervalSet([1.0 .. 3.0, 5.0 .. 10.0]))
+
     function testsets(a, b)
         @test area(a ∪ b) ≤ area(myunion(a)) + area(myunion(b))
         @test area(setdiff(a, b)) ≤ area(myunion(a))
@@ -43,36 +60,53 @@ end
         @test issetequal(a, a)
         @test isdisjoint(setdiff(a, b), b)
         @test !isdisjoint(a, a)
-        
-        intersections = find_intersections(a, b)
-        a, b = vcat(a), vcat(b)  # Always make `a` / `b` vectors
+
+        intersections = find_intersections(convert(Array, a), convert(Array, b))
 
         # verify that all indices returned in `find_intersections` correspond to sets
         # in b that overlap with the given set in a
         @test all(enumerate(intersections)) do (i, x)
-            isempty(x) || !isempty(intersect(a[i:i], b[x]))
+            isempty(x) || !isempty(intersect(IntervalSet(a.items[i]), IntervalSet(b.items[x])))
         end
 
         # verify that all indices not returned in `find_intersections` correspond to
         # sets in b that do not overlap with the given set in akk
         @test all(enumerate(intersections)) do (i, x)
-            isempty(intersect(a[i:i], b[Not(x)]))
+            isempty(intersect(IntervalSet(a.items[i]), IntervalSet(b.items[Not(x)])))
+        end
+
+        # Test f(interval, intervalset) and f(intervalset, interval) methods work the same
+        # as the regular f(intervalset, intervalset) methods
+        if a isa IntervalSet && b isa IntervalSet
+            for f in (intersect, union, setdiff, symdiff, isdisjoint, issubset)
+                # We'll just use the first and last entries as sample intervals for these
+                # tests rather than an exhausting search
+                x = first(a.items)
+                @test f(x, b) == f(IntervalSet([x]), b)
+                x = last(a.items)
+                @test f(x, b) == f(IntervalSet([x]), b)
+
+                y = first(b.items)
+                @test f(a, y) == f(a, IntervalSet([y]))
+                y = last(b.items)
+                @test f(a, y) == f(a, IntervalSet([y]))
+            end
         end
     end
 
     # verify empty interval set
     @test isempty(union(Interval[]))
-    
+
     # a few taylored interval sets
-    a = [Interval(i, i + 3) for i in 1:5:15]
-    b = a .+ (1:2:5)
-    @test all(first.(a) .∈ a)
+    a = IntervalSet([Interval(i, i + 3) for i in 1:5:15])
+    b = IntervalSet(a.items .+ (1:2:5))
+    @test all(x -> first(x) ∈ a, a.items)
     testsets(a, b)
-    testsets(a[1:1], b)
-    testsets(a, b[1:1])
+    testsets(IntervalSet(a.items[1]), b)
+    testsets(a, IntervalSet(b.items[1]))
 
     # verify that `last` need not be ordered
-    intervals = [Interval(0, 5), Interval(0, 3)]
+    intervals = IntervalSet([Interval(0, 5), Interval(0, 3)])
     @test superset(union(intervals)) == Interval(0, 5)
 
     # try out some more involved (random) intervals
@@ -82,47 +116,47 @@ end
     ends = starts .+ rand(rng, 1:10_000, n)
     offsets = round.(Int, (ends .- starts) .* (2 .* rand(rng, n) .- 1))
 
-    a = Interval.(starts, ends)
-    b = Interval.(starts .+ offsets, ends .+ offsets)
-    @test all(first.(a) .∈ a)
+    a = IntervalSet(Interval.(starts, ends))
+    b = IntervalSet(Interval.(starts .+ offsets, ends .+ offsets))
+    @test all(x -> first(x) ∈ a, a.items)
     testsets(a, b)
-    testsets(a[1:1], b)
-    testsets(a, b[1:1])
+    testsets(IntervalSet(first(a.items)), b)
+    testsets(a, IntervalSet(first(b.items)))
 
-    a = Interval{rand_bound_type(rng), rand_bound_type(rng)}.(starts, ends)
-    b = Interval{rand_bound_type(rng), rand_bound_type(rng)}.(starts .+ offsets, ends .+ offsets)
+    a = IntervalSet(Interval{rand_bound_type(rng), rand_bound_type(rng)}.(starts, ends))
+    b = IntervalSet(Interval{rand_bound_type(rng), rand_bound_type(rng)}.(starts .+ offsets, ends .+ offsets))
     testsets(a, b)
-    testsets(a[1:1], b)
-    testsets(a, b[1:1])
+    testsets(IntervalSet(first(a.items)), b)
+    testsets(a, IntervalSet(first(b.items)))
 
-    a = Interval{Closed, Open}.(starts, ends)
-    b = Interval{Closed, Open}.(starts .+ offsets, ends .+ offsets)
+    a = IntervalSet(Interval{Closed, Open}.(starts, ends))
+    b = IntervalSet(Interval{Closed, Open}.(starts .+ offsets, ends .+ offsets))
     @test Intervals.endpoint_tracking(a, b) isa Intervals.TrackStatically
-    @test Intervals.endpoint_tracking(a[1:1], b) isa Intervals.TrackStatically
-    @test Intervals.endpoint_tracking(a, b[1:1]) isa Intervals.TrackStatically
+    @test Intervals.endpoint_tracking(IntervalSet(first(a.items)), b) isa Intervals.TrackStatically
+    @test Intervals.endpoint_tracking(a, IntervalSet(first(b.items))) isa Intervals.TrackStatically
     testsets(a, b)
-    testsets(a[1:1], b)
-    testsets(a, b[1:1])
+    testsets(IntervalSet(first(a.items)), b)
+    testsets(a, IntervalSet(first(b.items)))
 
-    a = Interval{Open, Closed}.(starts, ends)
-    b = Interval{Open, Closed}.(starts .+ offsets, ends .+ offsets)
+    a = IntervalSet(Interval{Open, Closed}.(starts, ends))
+    b = IntervalSet(Interval{Open, Closed}.(starts .+ offsets, ends .+ offsets))
     @test Intervals.endpoint_tracking(a, b) isa Intervals.TrackStatically
-    @test Intervals.endpoint_tracking(a[1:1], b) isa Intervals.TrackStatically
-    @test Intervals.endpoint_tracking(a, b[1:1]) isa Intervals.TrackStatically
+    @test Intervals.endpoint_tracking(IntervalSet(first(a.items)), b) isa Intervals.TrackStatically
+    @test Intervals.endpoint_tracking(a, IntervalSet(first(b.items))) isa Intervals.TrackStatically
     testsets(a, b)
-    testsets(a[1:1], b)
-    testsets(a, b[1:1])
+    testsets(IntervalSet(first(a.items)), b)
+    testsets(a, IntervalSet(first(b.items)))
 
     randint(x::Interval) = Interval{rand_bound_type(rng), rand_bound_type(rng)}(first(x), last(x))
     leftint(x::Interval) = Interval{Closed, Open}(first(x), last(x))
     rightint(x::Interval) = Interval{Open, Closed}(first(x), last(x))
 
-    a = Interval{Closed, Open}.(starts, ends)
-    b = Interval{Closed, Open}.(starts .+ offsets, ends .+ offsets)
-    testsets(a, randint.(b))
-    testsets(a[1:1], randint.(b))
-    testsets(a, leftint.(b[1:1]))
-    testsets(a, rightint.(b))
-    testsets(a[1:1], rightint.(b))
-    testsets(a, rightint.(b[1:1]))
+    a = IntervalSet(Interval{Closed, Open}.(starts, ends))
+    b = IntervalSet(Interval{Closed, Open}.(starts .+ offsets, ends .+ offsets))
+    testsets(a, IntervalSet(randint.(b.items)))
+    testsets(IntervalSet(first(a.items)), IntervalSet(randint.(b.items)))
+    testsets(a, IntervalSet(leftint.(first(b.items))))
+    testsets(a, IntervalSet(rightint.(b.items)))
+    testsets(IntervalSet(first(a.items)), IntervalSet(rightint.(b.items)))
+    testsets(a, IntervalSet(rightint.(first(b.items))))
 end