From 105e0a2144b2ba2b173816b83c3174f36d7585cc Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 11 May 2026 12:18:22 +0200
Subject: [PATCH 01/15] Add JointInterval

---
 Project.toml                     | 10 +++++-----
 src/UncertaintyQuantification.jl |  3 +++
 src/inputs/imprecise/interval.jl | 20 ++++++++++++++++++++
 src/inputs/transportmaps.jl      |  4 ++--
 4 files changed, 30 insertions(+), 7 deletions(-)

diff --git a/Project.toml b/Project.toml
index da0f383e3..fecf04012 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,11 +1,7 @@
 name = "UncertaintyQuantification"
 uuid = "7183a548-a887-11e9-15ce-a56ab60bad7a"
 version = "0.14.0"
-authors = [
-    "Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>",
-    "Ander Gray <ander.gray@polytechnique.edu>",
-]
-
+authors = ["Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>", "Ander Gray <ander.gray@polytechnique.edu>"]
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -22,12 +18,14 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
+LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MeshAdaptiveDirectSearch = "f4d74008-4565-11e9-04bd-4fe404e6a92a"
 Monomials = "272bfe72-f66c-432f-a94d-600f29493792"
 Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
 Mustache = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
@@ -51,11 +49,13 @@ Distributions = "0.24, 0.25"
 FastGaussQuadrature = "0.4, 0.5, 1"
 FiniteDifferences = "0.12"
 Format = "1.3"
+LazySets = "6.0.3"
 MeshAdaptiveDirectSearch = "0.1.0"
 Monomials = "1.0"
 Mooncake = "0.4.196, 0.5"
 Mustache = "1.0"
 Optim = "1.9.4, 2.0"
+Polyhedra = "0.8.1"
 Primes = "0.5"
 QuadGK = "2.11.1"
 QuasiMonteCarlo = "0.3"
diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl
index 4902f4c5d..067356097 100644
--- a/src/UncertaintyQuantification.jl
+++ b/src/UncertaintyQuantification.jl
@@ -13,12 +13,14 @@ using Distributed
 using FastGaussQuadrature
 using FiniteDifferences
 using Format
+using LazySets
 using LinearAlgebra
 using MeshAdaptiveDirectSearch
 using Monomials
 using Mooncake: Mooncake
 using Mustache
 using Optim
+using Polyhedra
 using Primes
 using QuadGK
 using QuasiMonteCarlo
@@ -126,6 +128,7 @@ export ImportanceSampling
 export Interval
 export IntervalVariable
 export JointDistribution
+export JointInterval
 export KanaiTajimi
 export LaplaceEstimateBayesian
 export LatinHypercubeSampling
diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index 10ec31c34..2e5c561d0 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -107,3 +107,23 @@ Base.in(u, i::IntervalVariable) = i.lb <= u <= i.ub
 
 mean(i::IntervalVariable) = Interval(i)
 var(i::IntervalVariable) = Interval(0.0, (i.ub - i.lb)^2 / 4)
+
+
+struct JointInterval <: UQInput
+    intervals::Vector{IntervalVariable}
+    hull::VPolytope
+    function JointInterval(x::AbstractMatrix, names::AbstractVector{Symbol})
+        if size(x,1) != length(names)
+            error("Number of names must match dimensionality of points")
+        end
+        lb = vec(minimum(x,dims=2))
+        ub = vec(maximum(x, dims=2))
+
+        intervals = IntervalVariable.(lb,ub, names)
+        hull = VPolytope(convex_hull(collect(eachcol(x))))
+        return new(intervals, hull)
+    end
+end
+
+Base.in(x::AbstractVector, ji::JointInterval) = x ∈ ji.hull
+bounds(ji::JointInterval) = bounds.(ji.intervals)
\ No newline at end of file
diff --git a/src/inputs/transportmaps.jl b/src/inputs/transportmaps.jl
index 5d96b9651..9e9477b56 100644
--- a/src/inputs/transportmaps.jl
+++ b/src/inputs/transportmaps.jl
@@ -231,7 +231,7 @@ end
 A transport map constructed from samples in the physical space X which are mapped to the standard normal space Z. The `map` is the optimized composed map, `samples` is a DataFrame of samples in the physical space X used to fit the map, and `names` is a vector of variable names.
 """
 struct TransportMapFromSamples <: AbstractTransportMap
-    map::ComposedMap{LinearMap}
+    map::ComposedMap{TransportMaps.LinearMap}
     samples::DataFrame
     names::Vector{Symbol}
 end
@@ -274,7 +274,7 @@ function mapfromsamples(
     target_samples = Matrix(X)
 
     # First, fit a linear map
-    linear_map = LinearMap(target_samples)
+    linear_map = TransportMaps.LinearMap(target_samples)
     # Optimize transportmap
     optimize!(
         transportmap, target_samples, linear_map; optimizer=optimizer, options=options

From 1a8c4dbaa6e4a4067cabaa7123b3ce69c6d8a6d3 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 11 May 2026 12:57:35 +0200
Subject: [PATCH 02/15] Add JointInterval to isimprecise

---
 src/inputs/imprecise/interval.jl | 9 ++++-----
 src/util/imprecise.jl            | 1 +
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index 2e5c561d0..18c8e174d 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -108,18 +108,17 @@ Base.in(u, i::IntervalVariable) = i.lb <= u <= i.ub
 mean(i::IntervalVariable) = Interval(i)
 var(i::IntervalVariable) = Interval(0.0, (i.ub - i.lb)^2 / 4)
 
-
 struct JointInterval <: UQInput
     intervals::Vector{IntervalVariable}
     hull::VPolytope
     function JointInterval(x::AbstractMatrix, names::AbstractVector{Symbol})
-        if size(x,1) != length(names)
+        if size(x, 1) != length(names)
             error("Number of names must match dimensionality of points")
         end
-        lb = vec(minimum(x,dims=2))
-        ub = vec(maximum(x, dims=2))
+        lb = vec(minimum(x; dims=2))
+        ub = vec(maximum(x; dims=2))
 
-        intervals = IntervalVariable.(lb,ub, names)
+        intervals = IntervalVariable.(lb, ub, names)
         hull = VPolytope(convex_hull(collect(eachcol(x))))
         return new(intervals, hull)
     end
diff --git a/src/util/imprecise.jl b/src/util/imprecise.jl
index 372500831..2eabfa7eb 100644
--- a/src/util/imprecise.jl
+++ b/src/util/imprecise.jl
@@ -4,6 +4,7 @@ end
 
 function isimprecise(input::UQInput)
     return isa(input, IntervalVariable) ||
+           isa(input, JointInterval) ||
            isa(input, RandomVariable{<:ProbabilityBox}) ||
            (
                isa(input, JointDistribution{<:Copulas.Copula,<:RandomVariable}) &&

From 7ca73efd8f59db42a6870deabf6c96dede829ea2 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 11 May 2026 13:03:10 +0200
Subject: [PATCH 03/15] Sample JointInterval

---
 src/inputs/imprecise/interval.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index 18c8e174d..7f0f0c76f 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -125,4 +125,5 @@ struct JointInterval <: UQInput
 end
 
 Base.in(x::AbstractVector, ji::JointInterval) = x ∈ ji.hull
-bounds(ji::JointInterval) = bounds.(ji.intervals)
\ No newline at end of file
+bounds(ji::JointInterval) = bounds.(ji.intervals)
+sample(ji::JointInterval, n::Integer=1) = sample(ji.intervals, n)
\ No newline at end of file

From d159f7cd4511e31dbedcd62cbea7f44f810edd00 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 11 May 2026 13:05:32 +0200
Subject: [PATCH 04/15] Add empty SNS mappings

---
 src/inputs/imprecise/interval.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index 7f0f0c76f..a0e5270c1 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -126,4 +126,6 @@ end
 
 Base.in(x::AbstractVector, ji::JointInterval) = x ∈ ji.hull
 bounds(ji::JointInterval) = bounds.(ji.intervals)
-sample(ji::JointInterval, n::Integer=1) = sample(ji.intervals, n)
\ No newline at end of file
+sample(ji::JointInterval, n::Integer=1) = sample(ji.intervals, n)
+to_standard_normal_space!(_::JointInterval, _::DataFrame) = nothing
+to_physical_space!(_::JointInterval, _::DataFrame) = nothing
\ No newline at end of file

From 92222c898830979f1c9fc21a652203571cba85af Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 11 May 2026 13:06:00 +0200
Subject: [PATCH 05/15] Add JointInterval to QMC sampling

---
 src/simulations/montecarlo.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/simulations/montecarlo.jl b/src/simulations/montecarlo.jl
index 439a02a34..8adb1ef9d 100644
--- a/src/simulations/montecarlo.jl
+++ b/src/simulations/montecarlo.jl
@@ -76,7 +76,8 @@ end
 
 function sample(inputs::Vector{<:UQInput}, sim::AbstractQuasiMonteCarlo)
     random_inputs = filter(i -> isa(i, RandomUQInput) || isa(i, ProbabilityBox), inputs)
-    deterministic_inputs = filter(i -> isa(i, Parameter) || isa(i, Interval), inputs)
+    deterministic_inputs =
+        filter(i -> isa(i, Parameter) || isa(i, Interval), inputs) || isa(i, JointInterval)
 
     n_rv = count_rvs(random_inputs)
 

From de1c8c7e73bb5d1a27f0dda1dedac393c286a751 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Tue, 12 May 2026 13:59:30 +0200
Subject: [PATCH 06/15] Add sampling of dependent intervals

---
 src/inputs/imprecise/interval.jl |  3 ++
 src/inputs/inputs.jl             |  2 +-
 src/simulations/montecarlo.jl    | 70 +++++++++++++++++++++++++++-----
 3 files changed, 64 insertions(+), 11 deletions(-)

diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index a0e5270c1..708c7f4f4 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -107,6 +107,7 @@ Base.in(u, i::IntervalVariable) = i.lb <= u <= i.ub
 
 mean(i::IntervalVariable) = Interval(i)
 var(i::IntervalVariable) = Interval(0.0, (i.ub - i.lb)^2 / 4)
+dimensions(i::IntervalVariable) = 1
 
 struct JointInterval <: UQInput
     intervals::Vector{IntervalVariable}
@@ -125,6 +126,8 @@ struct JointInterval <: UQInput
 end
 
 Base.in(x::AbstractVector, ji::JointInterval) = x ∈ ji.hull
+names(ji::JointInterval) = names(ji.intervals)
+dimensions(ji::JointInterval) = length(ji.intervals)
 bounds(ji::JointInterval) = bounds.(ji.intervals)
 sample(ji::JointInterval, n::Integer=1) = sample(ji.intervals, n)
 to_standard_normal_space!(_::JointInterval, _::DataFrame) = nothing
diff --git a/src/inputs/inputs.jl b/src/inputs/inputs.jl
index 6a169b1f7..96a0677b6 100644
--- a/src/inputs/inputs.jl
+++ b/src/inputs/inputs.jl
@@ -27,7 +27,7 @@ function names(inputs::Vector{<:UQInput})
     _names = Symbol[]
 
     for i in inputs
-        if i isa JointDistribution || i isa SpectralRepresentation
+        if i isa JointDistribution || i isa SpectralRepresentation || i isa JointInterval
             append!(_names, names(i))
         else
             push!(_names, getproperty(i, :name))
diff --git a/src/simulations/montecarlo.jl b/src/simulations/montecarlo.jl
index 8adb1ef9d..652c25393 100644
--- a/src/simulations/montecarlo.jl
+++ b/src/simulations/montecarlo.jl
@@ -74,22 +74,72 @@ function sample(inputs::Vector{<:UQInput}, sim::MonteCarlo)
     return sample(inputs, sim.n)
 end
 
-function sample(inputs::Vector{<:UQInput}, sim::AbstractQuasiMonteCarlo)
-    random_inputs = filter(i -> isa(i, RandomUQInput) || isa(i, ProbabilityBox), inputs)
-    deterministic_inputs =
-        filter(i -> isa(i, Parameter) || isa(i, Interval), inputs) || isa(i, JointInterval)
+function sample(
+    inputs::Vector{<:UQInput}, sim::AbstractQuasiMonteCarlo; intervals::Bool=true
+)
+    rvs = filter(i -> isa(i, RandomUQInput), inputs)
+    ivs = filter(i -> isa(i, IntervalVariable) || isa(i, JointInterval), inputs)
+    parameters = filter(i -> isa(i, Parameter), inputs)
+
+    dependent = any(isa.(ivs, JointInterval))
+    # verify randomized QMC for dependent intervals
+    if dependent &
+        !intervals &
+        !(isa(sim, LatinHypercubeSampling) || isa(sim, RandomizedHaltonSample))
+        if sim.randomization == :none
+            error("QMC sampling must be randomized for joint intervals")
+        end
+    end
 
-    n_rv = count_rvs(random_inputs)
+    n_rv = count_rvs(rvs)
+    n_int = !isempty(ivs) ? mapreduce(dimensions, +, ivs) : 0
 
-    u = qmc_samples(sim, n_rv)
+    # if intervals is true only rvs need to be sampled. If not we also obtain qmc samples for intervals
+    # if dependent intervals are involved sample much more than requested
+    u = qmc_samples(
+        dependent & !intervals ? typeof(sim)(10^6, sim.randomization) : sim,
+        intervals ? n_rv : n_rv + n_int,
+    )
+
+    samples = if intervals
+        DataFrame(names(rvs) .=> eachrow(u))
+    else
+        DataFrame(vcat(names(rvs), names(ivs)) .=> eachrow(u))
+    end
 
-    samples = quantile.(Normal(), u)
-    samples = DataFrame(names(random_inputs) .=> eachrow(samples))
+    # map rvs into standard normal space
+    samples[:, names(rvs)] = quantile.(Normal(), samples[:, names(rvs)])
+
+    if !isempty(ivs)
+        if intervals
+            # append intervals if not sampled
+            DataFrames.hcat!(samples, sample(ivs, size(samples, 1)))
+        else
+            # translate qmc samples to interval ranges
+            for i in mapreduce(i -> isa(i, IntervalVariable) ? i : i.intervals, vcat, ivs)
+                samples[:, i.name] = samples[:, i.name] .* (i.ub - i.lb) .+ i.lb
+            end
+        end
+    end
+
+    # discard samples outside the permissible set
+    if dependent & !intervals
+        for ji in filter(i -> isa(i, JointInterval), ivs)
+            idx = findall(.!in.(eachrow(Matrix(samples[:, names(ji)])), ji))
+            deleteat!(samples, idx)
+        end
+        # discard excess samples
+        if size(samples)[1] > sim.n
+            deleteat!(samples, (sim.n + 1):size(samples)[1])
+        end
+    end
 
-    if !isempty(deterministic_inputs)
-        DataFrames.hcat!(samples, sample(deterministic_inputs, size(samples, 1)))
+    # finally append any parameters
+    if !isempty(parameters)
+        DataFrames.hcat!(samples, sample(parameters, size(samples, 1)))
     end
 
+    # map rvs to physical space before returning the samples
     to_physical_space!(inputs, samples)
 
     return samples

From 192689a3b6c046f26cd913b24143547ffa8b5a25 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Tue, 12 May 2026 14:48:49 +0200
Subject: [PATCH 07/15] Add docstring

---
 Project.toml                     |  2 ++
 src/UncertaintyQuantification.jl |  1 +
 src/simulations/montecarlo.jl    | 28 ++++++++++++++++++++--------
 3 files changed, 23 insertions(+), 8 deletions(-)

diff --git a/Project.toml b/Project.toml
index fecf04012..9aca2a68c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -34,6 +34,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TransportMaps = "6fd49bf0-03ac-4dd1-b9b3-42e511274d6e"
 
 [compat]
@@ -63,5 +64,6 @@ Reexport = "0.2, 1.0"
 Roots = "2.2.2, 3.0"
 Statistics = "1"
 StatsBase = "0.33, 0.34"
+Suppressor = "0.2.8"
 TransportMaps = "0.1.5"
 julia = "1.10"
diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl
index 067356097..614895fbc 100644
--- a/src/UncertaintyQuantification.jl
+++ b/src/UncertaintyQuantification.jl
@@ -28,6 +28,7 @@ using Random
 using Reexport
 using Roots
 using StatsBase
+using Suppressor
 using TransportMaps
 
 @reexport using TransportMaps
diff --git a/src/simulations/montecarlo.jl b/src/simulations/montecarlo.jl
index 652c25393..da530e2bf 100644
--- a/src/simulations/montecarlo.jl
+++ b/src/simulations/montecarlo.jl
@@ -74,8 +74,16 @@ function sample(inputs::Vector{<:UQInput}, sim::MonteCarlo)
     return sample(inputs, sim.n)
 end
 
+"""
+    sample(inputs::Vector{UQInput}, sim::AbstractQuasiMonteCarlo; intervals::Bool=true, n_internal::Integer=sim.n * 10)
+
+    Generate Quasi-Monte Carlo samples of the `inputs` using the QMC sampling method `sim`. By default any [`IntervalVariable`](@ref) or [`'JointInterval`]('ref) will be included as the original intervals. To apply the QMC sampling also to the intervals pass the keyword `;intervals=true`. For [`JointInterval`]('ref)s this will internally sample 10 times the desired samples and discard the samples outside the permissible set and any excess samples. If the permissible set is small more samples might be required to generated sufficient samples. In this case, the keyword `n_internal` can be used to increase the number of samples used.
+"""
 function sample(
-    inputs::Vector{<:UQInput}, sim::AbstractQuasiMonteCarlo; intervals::Bool=true
+    inputs::Vector{<:UQInput},
+    sim::AbstractQuasiMonteCarlo;
+    intervals::Bool=true,
+    n_internal::Integer=sim.n * 10,
 )
     rvs = filter(i -> isa(i, RandomUQInput), inputs)
     ivs = filter(i -> isa(i, IntervalVariable) || isa(i, JointInterval), inputs)
@@ -87,19 +95,21 @@ function sample(
         !intervals &
         !(isa(sim, LatinHypercubeSampling) || isa(sim, RandomizedHaltonSample))
         if sim.randomization == :none
-            error("QMC sampling must be randomized for joint intervals")
+            error("QMC sampling must be randomized to be applied to joint intervals")
         end
     end
 
     n_rv = count_rvs(rvs)
     n_int = !isempty(ivs) ? mapreduce(dimensions, +, ivs) : 0
 
-    # if intervals is true only rvs need to be sampled. If not we also obtain qmc samples for intervals
-    # if dependent intervals are involved sample much more than requested
-    u = qmc_samples(
-        dependent & !intervals ? typeof(sim)(10^6, sim.randomization) : sim,
-        intervals ? n_rv : n_rv + n_int,
-    )
+    u = @suppress_err begin
+        # if intervals is true only rvs need to be sampled. If not we also obtain qmc samples for intervals
+        # if dependent intervals are involved sample much more than requested
+        qmc_samples(
+            dependent & !intervals ? typeof(sim)(n_internal, sim.randomization) : sim,
+            intervals ? n_rv : n_rv + n_int,
+        )
+    end
 
     samples = if intervals
         DataFrame(names(rvs) .=> eachrow(u))
@@ -131,6 +141,8 @@ function sample(
         # discard excess samples
         if size(samples)[1] > sim.n
             deleteat!(samples, (sim.n + 1):size(samples)[1])
+        else
+            @warn "Only $(size(samples)[1]) of $(sim.n) samples generated. Try increasing 'n_internal'"
         end
     end
 

From 926ee9064c6cef76e15238caa692a88bed49bfc7 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Tue, 12 May 2026 16:36:24 +0200
Subject: [PATCH 08/15] Fix n-dimensional JointInterval

---
 Project.toml                      |  4 ++--
 src/UncertaintyQuantification.jl  |  1 -
 src/inputs/imprecise/interval.jl  |  2 +-
 src/simulations/montecarlo.jl     | 14 ++++++--------
 test/inputs/imprecise/interval.jl | 13 +++++++++++++
 5 files changed, 22 insertions(+), 12 deletions(-)

diff --git a/Project.toml b/Project.toml
index 9aca2a68c..7b488d56a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,6 +18,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MeshAdaptiveDirectSearch = "f4d74008-4565-11e9-04bd-4fe404e6a92a"
@@ -34,7 +35,6 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TransportMaps = "6fd49bf0-03ac-4dd1-b9b3-42e511274d6e"
 
 [compat]
@@ -50,6 +50,7 @@ Distributions = "0.24, 0.25"
 FastGaussQuadrature = "0.4, 0.5, 1"
 FiniteDifferences = "0.12"
 Format = "1.3"
+GLPK = "1.2.1"
 LazySets = "6.0.3"
 MeshAdaptiveDirectSearch = "0.1.0"
 Monomials = "1.0"
@@ -64,6 +65,5 @@ Reexport = "0.2, 1.0"
 Roots = "2.2.2, 3.0"
 Statistics = "1"
 StatsBase = "0.33, 0.34"
-Suppressor = "0.2.8"
 TransportMaps = "0.1.5"
 julia = "1.10"
diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl
index 614895fbc..067356097 100644
--- a/src/UncertaintyQuantification.jl
+++ b/src/UncertaintyQuantification.jl
@@ -28,7 +28,6 @@ using Random
 using Reexport
 using Roots
 using StatsBase
-using Suppressor
 using TransportMaps
 
 @reexport using TransportMaps
diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index 708c7f4f4..6d0a4d3bf 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -120,7 +120,7 @@ struct JointInterval <: UQInput
         ub = vec(maximum(x; dims=2))
 
         intervals = IntervalVariable.(lb, ub, names)
-        hull = VPolytope(convex_hull(collect(eachcol(x))))
+        hull = VPolytope(convex_hull([x[:, i] for i in axes(x, 2)]))
         return new(intervals, hull)
     end
 end
diff --git a/src/simulations/montecarlo.jl b/src/simulations/montecarlo.jl
index da530e2bf..5c0a6ff1c 100644
--- a/src/simulations/montecarlo.jl
+++ b/src/simulations/montecarlo.jl
@@ -102,14 +102,12 @@ function sample(
     n_rv = count_rvs(rvs)
     n_int = !isempty(ivs) ? mapreduce(dimensions, +, ivs) : 0
 
-    u = @suppress_err begin
-        # if intervals is true only rvs need to be sampled. If not we also obtain qmc samples for intervals
-        # if dependent intervals are involved sample much more than requested
-        qmc_samples(
-            dependent & !intervals ? typeof(sim)(n_internal, sim.randomization) : sim,
-            intervals ? n_rv : n_rv + n_int,
-        )
-    end
+    # if intervals is true only rvs need to be sampled. If not we also obtain qmc samples for intervals
+    # if dependent intervals are involved sample much more than requested
+    u = qmc_samples(
+        dependent & !intervals ? typeof(sim)(n_internal, sim.randomization) : sim,
+        intervals ? n_rv : n_rv + n_int,
+    )
 
     samples = if intervals
         DataFrame(names(rvs) .=> eachrow(u))
diff --git a/test/inputs/imprecise/interval.jl b/test/inputs/imprecise/interval.jl
index 54d202a10..cc7d90d7d 100644
--- a/test/inputs/imprecise/interval.jl
+++ b/test/inputs/imprecise/interval.jl
@@ -57,3 +57,16 @@ end
 
     @test UncertaintyQuantification.sample(interval) == DataFrame(; x=Interval(interval))
 end
+
+@testset "JointInterval" begin
+    x = rand(3, 10)
+    ji = JointInterval(x, [:x1, :x2, :x3])
+
+    @test names(ji) == [:x1, :x2, :x3]
+    for i in 1:3
+        @test ji.intervals[i].lb == minimum(x[i, :])
+        @test ji.intervals[i].ub == maximum(x[i, :])
+    end
+
+    @test all(in.(eachcol(x), ji))
+end

From d931a7be0d11d10e6e8bb908258d63577ba0b6eb Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Tue, 12 May 2026 16:37:32 +0200
Subject: [PATCH 09/15] Add JointInterval to API docs

---
 docs/src/api/inputs.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/api/inputs.md b/docs/src/api/inputs.md
index 8cc1ea200..969dd5d1c 100644
--- a/docs/src/api/inputs.md
+++ b/docs/src/api/inputs.md
@@ -16,6 +16,7 @@ RandomVariable
 IntervalVariable
 EmpiricalDistribution
 JointDistribution
+JointInterval
 ```
 
 ## Functions

From 7454017af4ba561523b32ebee4e3543ae3bb37b1 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Wed, 13 May 2026 11:04:48 +0200
Subject: [PATCH 10/15] Use TestItems

---
 src/inputs/imprecise/interval.jl  |   2 +-
 test/Project.toml                 |   2 +
 test/inputs/imprecise/interval.jl | 148 +++++++++++++++---------------
 test/runtests.jl                  |   4 +
 4 files changed, 83 insertions(+), 73 deletions(-)

diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index 6d0a4d3bf..bc7c5ffa0 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -128,7 +128,7 @@ end
 Base.in(x::AbstractVector, ji::JointInterval) = x ∈ ji.hull
 names(ji::JointInterval) = names(ji.intervals)
 dimensions(ji::JointInterval) = length(ji.intervals)
-bounds(ji::JointInterval) = bounds.(ji.intervals)
+bounds(ji::JointInterval) = bounds(ji.intervals)
 sample(ji::JointInterval, n::Integer=1) = sample(ji.intervals, n)
 to_standard_normal_space!(_::JointInterval, _::DataFrame) = nothing
 to_physical_space!(_::JointInterval, _::DataFrame) = nothing
\ No newline at end of file
diff --git a/test/Project.toml b/test/Project.toml
index 37cc67ced..0f1b82e14 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -9,3 +9,5 @@ QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
+TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
diff --git a/test/inputs/imprecise/interval.jl b/test/inputs/imprecise/interval.jl
index cc7d90d7d..8f458bd80 100644
--- a/test/inputs/imprecise/interval.jl
+++ b/test/inputs/imprecise/interval.jl
@@ -1,72 +1,76 @@
-@testset "Interval" begin
-    name = :l
-    lb = 0.14
-    ub = 0.16
-    @test_throws ErrorException(
-        "Lower bound parameter must be smaller than upper bound parameter for Intervals."
-    ) Interval(ub, lb)
-    interval = Interval(lb, ub)
-    @test interval.lb == lb
-    @test interval.ub == ub
-
-    @test !(0.13 ∈ interval)
-    @test 0.14 ∈ interval
-    @test 0.15 ∈ interval
-    @test 0.16 ∈ interval
-    @test !(0.17 ∈ interval)
-
-    @test sprint(show, interval) == "[0.14, 0.16]"
-end
-
-@testset "IntervalVariable" begin
-    name = :l
-    lb = 0.14
-    ub = 0.16
-    @test_throws ErrorException(
-        "Lower bound parameter must be smaller than upper bound parameter for Interval $name.",
-    ) IntervalVariable(ub, lb, name)
-    interval = IntervalVariable(lb, ub, name)
-    @test interval.lb == lb
-    @test interval.ub == ub
-    @test interval.name == name
-    @test !(0.13 ∈ interval)
-    @test 0.14 ∈ interval
-    @test 0.15 ∈ interval
-    @test 0.16 ∈ interval
-    @test !(0.17 ∈ interval)
-
-    @test sprint(show, interval) == "l ∈ [0.14, 0.16]"
-
-    par = 0.13
-    @test_throws ErrorException("0.13 not in [0.14, 0.16] for Interval l.") UncertaintyQuantification.map_to_precise(
-        par, interval
-    )
-
-    par = 0.17
-    @test_throws ErrorException("0.17 not in [0.14, 0.16] for Interval l.") UncertaintyQuantification.map_to_precise(
-        par, interval
-    )
-
-    @test UncertaintyQuantification.map_to_precise(0.15, interval) ==
-        Parameter(0.15, interval.name)
-
-    interval = IntervalVariable(0, 1, :x)
-
-    @test mean(interval) == Interval(0, 1)
-    @test var(interval) == Interval(0.0, 0.25)
-
-    @test UncertaintyQuantification.sample(interval) == DataFrame(; x=Interval(interval))
-end
-
-@testset "JointInterval" begin
-    x = rand(3, 10)
-    ji = JointInterval(x, [:x1, :x2, :x3])
-
-    @test names(ji) == [:x1, :x2, :x3]
-    for i in 1:3
-        @test ji.intervals[i].lb == minimum(x[i, :])
-        @test ji.intervals[i].ub == maximum(x[i, :])
-    end
-
-    @test all(in.(eachcol(x), ji))
-end
+@testitem "Interval" begin
+    name = :l
+    lb = 0.14
+    ub = 0.16
+    @test_throws ErrorException(
+        "Lower bound parameter must be smaller than upper bound parameter for Intervals."
+    ) Interval(ub, lb)
+    interval = Interval(lb, ub)
+    @test interval.lb == lb
+    @test interval.ub == ub
+
+    @test !(0.13 ∈ interval)
+    @test 0.14 ∈ interval
+    @test 0.15 ∈ interval
+    @test 0.16 ∈ interval
+    @test !(0.17 ∈ interval)
+
+    @test sprint(show, interval) == "[0.14, 0.16]"
+end
+
+@testitem "IntervalVariable" begin
+    using DataFrames
+    name = :l
+    lb = 0.14
+    ub = 0.16
+    @test_throws ErrorException(
+        "Lower bound parameter must be smaller than upper bound parameter for Interval $name.",
+    ) IntervalVariable(ub, lb, name)
+    interval = IntervalVariable(lb, ub, name)
+    @test interval.lb == lb
+    @test interval.ub == ub
+    @test interval.name == name
+    @test !(0.13 ∈ interval)
+    @test 0.14 ∈ interval
+    @test 0.15 ∈ interval
+    @test 0.16 ∈ interval
+    @test !(0.17 ∈ interval)
+
+    @test sprint(show, interval) == "l ∈ [0.14, 0.16]"
+
+    par = 0.13
+    @test_throws ErrorException("0.13 not in [0.14, 0.16] for Interval l.") UncertaintyQuantification.map_to_precise(
+        par, interval
+    )
+
+    par = 0.17
+    @test_throws ErrorException("0.17 not in [0.14, 0.16] for Interval l.") UncertaintyQuantification.map_to_precise(
+        par, interval
+    )
+
+    @test UncertaintyQuantification.map_to_precise(0.15, interval) ==
+        Parameter(0.15, interval.name)
+
+    interval = IntervalVariable(0, 1, :x)
+
+    @test mean(interval) == Interval(0, 1)
+    @test var(interval) == Interval(0.0, 0.25)
+
+    @test UncertaintyQuantification.sample(interval) == DataFrame(; x=Interval(interval))
+end
+
+@testitem "JointInterval" begin
+    x = rand(3, 10)
+    ji = JointInterval(x, [:x1, :x2, :x3])
+
+    @test names(ji) == [:x1, :x2, :x3]
+    @test UncertaintyQuantification.bounds(ji) ==
+        UncertaintyQuantification.bounds(ji.intervals)
+    @test dimensions(ji) == 3
+    for i in 1:3
+        @test ji.intervals[i].lb == minimum(x[i, :])
+        @test ji.intervals[i].ub == maximum(x[i, :])
+    end
+
+    @test all(in.(eachcol(x), ji))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index d85395aa7..b1c3b52d8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -8,8 +8,12 @@ using QuasiMonteCarlo
 using Random
 using StatsBase: fit, Histogram, corkendall
 using Test
+using TestItems
+using TestItemRunner
 using UncertaintyQuantification
 
+@run_package_tests
+
 include("inputs/empiricaldistribution.jl")
 include("dynamics/psd.jl")
 include("inputs/parameter.jl")

From 2af8e2f46c4d826654c2374e298602bdb9ade011 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Wed, 13 May 2026 11:34:38 +0200
Subject: [PATCH 11/15] Test JointInterval QMC sampling

---
 test/simulations/montecarlo.jl | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/test/simulations/montecarlo.jl b/test/simulations/montecarlo.jl
index 1894e5433..39ecac289 100644
--- a/test/simulations/montecarlo.jl
+++ b/test/simulations/montecarlo.jl
@@ -111,3 +111,27 @@ end
         @test sample(inputs, lattice) != sample(inputs, lattice)
     end
 end
+
+@testitem "JointInterval sampling" begin
+    x = rand(3, 10)
+    ji = JointInterval(x, [:x1, :x2, :x3])
+    rv = RandomVariable(Normal(), :y)
+    i = IntervalVariable(-1, 5, :z)
+
+    inputs = [i, ji, rv]
+
+    samples = sample(inputs, SobolSampling(256, :matousek); intervals=false)
+
+    @test size(samples)[1] == 256
+
+    @test all(insupport.(rv, samples.y))
+    @test all(in.(samples.z, i))
+    @test all(in.(samples.x1, ji.intervals[1]))
+    @test all(in.(samples.x2, ji.intervals[2]))
+    @test all(in.(samples.x3, ji.intervals[3]))
+    @test all(in.(eachrow(Matrix(samples[:, [:x1, :x2, :x3]])), ji))
+
+    @test_throws ErrorException(
+        "QMC sampling must be randomized to be applied to joint intervals"
+    ) sample(inputs, SobolSampling(256, :none); intervals=false)
+end

From acd896f3c56f83dd3dd7af38bd8904a7f9f53536 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 18 May 2026 10:13:31 +0200
Subject: [PATCH 12/15] Add constraints to DoubleLoop

---
 .../probabilityoffailure_imprecise.jl         | 23 +++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/src/reliability/probabilityoffailure_imprecise.jl b/src/reliability/probabilityoffailure_imprecise.jl
index 2ce200dba..1b204ca93 100644
--- a/src/reliability/probabilityoffailure_imprecise.jl
+++ b/src/reliability/probabilityoffailure_imprecise.jl
@@ -92,6 +92,23 @@ function probability_of_failure(
     imprecise_inputs = filter(x -> isimprecise(x), inputs)
     precise_inputs = filter(x -> !isimprecise(x), inputs)
 
+    imprecise_names = names(imprecise_inputs)
+    joint_int = filter(x -> isa(x, JointInterval), imprecise_inputs)
+    joint_int_idx = [findfirst(x -> x == n, imprecise_names) for n in names(joint_int)]
+
+    function constraints(x)
+        i = 1
+        for ji in joint_int
+            idx = joint_int_idx[i:(i + dimensions(ji) - 1)]
+            if !(x[idx] ∈ ji)
+                x, x[idx] ∈ ji
+                return false
+            end
+            i += dimensions(ji)
+        end
+        return true
+    end
+
     function pf_low(x)
         imprecise_inputs_x = map_to_precise_inputs(x, imprecise_inputs)
         mc_inputs = [precise_inputs..., imprecise_inputs_x...]
@@ -116,6 +133,7 @@ function probability_of_failure(
         lowerbound=lb,
         upperbound=ub,
         min_mesh_size=1e-13,
+        constraints=isempty(joint_int) ? [] : [constraints],
     )
 
     result_ub = minimize(
@@ -125,6 +143,7 @@ function probability_of_failure(
         lowerbound=lb,
         upperbound=ub,
         min_mesh_size=1e-13,
+        constraints=isempty(joint_int) ? [] : [constraints],
     )
 
     pf_lb = result_lb.f
@@ -167,6 +186,10 @@ function map_to_precise_inputs(x::AbstractVector, inputs::AbstractVector{<:UQInp
                 end
             end
             push!(precise_inputs, JointDistribution(i.d, precise_marginals))
+        elseif isa(i, JointInterval)
+            for iv in i.intervals
+                push!(precise_inputs, map_to_precise(popfirst!(params), iv))
+            end
         end
     end
     return precise_inputs

From ae78199599b4eba4ec7034c6387e19cbf5353885 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 18 May 2026 16:46:39 +0200
Subject: [PATCH 13/15] Set fail-fast to false

---
 .github/workflows/ci.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index d5af33039..0c8e30e09 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -22,6 +22,7 @@ jobs:
           - "ubuntu-latest"
           - "windows-latest"
           - "macos-latest"
+        fail-fast: false
     steps:
       - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v3

From 38b54e8bacfdf0f11561036b832850ef71a325e2 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 18 May 2026 16:54:06 +0200
Subject: [PATCH 14/15] Fix workflow file

---
 .github/workflows/ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0c8e30e09..fa6e6b220 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -14,6 +14,7 @@ jobs:
     name: Julia ${{ matrix.julia-version }} - ${{ matrix.os }}
     runs-on: ${{ matrix.os }}
     strategy:
+      fail-fast: false
       matrix:
         julia-version:
           - "lts"
@@ -22,7 +23,6 @@ jobs:
           - "ubuntu-latest"
           - "windows-latest"
           - "macos-latest"
-        fail-fast: false
     steps:
       - uses: actions/checkout@v6
       - uses: julia-actions/setup-julia@v3

From 6f613f0eeec392a0629961407f0d61b98ae7b173 Mon Sep 17 00:00:00 2001
From: Jasper Behrensdorf <behrensdorf@irz.uni-hannover.de>
Date: Mon, 18 May 2026 17:56:03 +0200
Subject: [PATCH 15/15] Disable show TM show test below julia 1.11

---
 test/inputs/transportmaps.jl | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/test/inputs/transportmaps.jl b/test/inputs/transportmaps.jl
index 754992787..1f1e83464 100644
--- a/test/inputs/transportmaps.jl
+++ b/test/inputs/transportmaps.jl
@@ -93,9 +93,12 @@
         @test insupport(tm, samp)
 
         # Test show methods
-        @test_nowarn sprint(show, tm.d)
-        @test_nowarn sprint(print, tm.d)
-        @test_nowarn display(tm.d)
+        # problematic on 1.10
+        if VERSION >= v"1.11"
+            @test_nowarn sprint(show, tm.d)
+            @test_nowarn sprint(print, tm.d)
+            @test_nowarn display(tm.d)
+        end
     end
 
     @testset "TransportMap from samples" begin
@@ -163,12 +166,12 @@
         tm_var = var(tm_samples.d)
         @test length(tm_var) == 2
         @test all(tm_var .> 0)
-        @test isapprox(tm_var, [1.5, 1.]; atol=3e-1)  # Approximation from samples
+        @test isapprox(tm_var, [1.5, 1.0]; atol=3e-1)  # Approximation from samples
 
         tm_std = std(tm_samples.d)
         @test length(tm_std) == 2
         @test all(tm_std .> 0)
-        @test isapprox(tm_std, sqrt.([1.5, 1.]); atol=2e-1)  # Approximation from samples
+        @test isapprox(tm_std, sqrt.([1.5, 1.0]); atol=2e-1)  # Approximation from samples
         @test isapprox(tm_std, sqrt.(tm_var); atol=1e-10)
 
         # Test with custom quadrature
@@ -186,9 +189,12 @@
         @test insupport(tm_samples, samp)
 
         # Test show methods
-        @test_nowarn sprint(show, tm_samples.d)
-        @test_nowarn sprint(print, tm_samples.d)
-        @test_nowarn display(tm_samples.d)
+        # problematic on 1.10
+        if VERSION >= v"1.11"
+            @test_nowarn sprint(show, tm_samples.d)
+            @test_nowarn sprint(print, tm_samples.d)
+            @test_nowarn display(tm_samples.d)
+        end
     end
 
     @testset "TransportMap with transform_density" begin
@@ -245,7 +251,7 @@
         @test isapprox(mean(df.x2), tm_mean[2]; atol=3e-1)
 
         # Test pdf
-        x_test = [-10., 0]
+        x_test = [-10.0, 0]
         @test pdf(tm, x_test) ≈ 0
         @test pdf(tm, [2, 0]) > 0