diff --git a/Project.toml b/Project.toml
index 9f83f9a..37ad6e9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ DiscreteMeasures = "7766d772-2108-41ee-a4bd-11c51440a39b"
 OUQBase = "01930cae-99d2-7439-8f4f-ace2ece9f1b9"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [sources]
 CanonicalMoments = { path = "lib/CanonicalMoments" }
@@ -24,6 +25,7 @@ Pkg = "1.10"
 PrecompileTools = "1.2"
 Reexport = "1.2.2"
 SafeTestsets = "0.1"
+Symbolics = "7.37"
 Test = "1.10"
 SciMLTesting = "2"
 
diff --git a/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl b/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl
index 62ae7ca..f42ef52 100644
--- a/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl
+++ b/lib/DiscreteMeasures/ext/IntervalArithmeticExt.jl
@@ -2,10 +2,9 @@ module IntervalArithmeticExt
 
 using DiscreteMeasures, IntervalArithmetic
 
-println("Load Ext")
 function DiscreteMeasures.clamp_domain(x::Interval, lb, ub)
     bounds = interval(lb, ub)
-    return intersect(x, bounds)
+    return intersect_interval(x, bounds)
 end
 
 function DiscreteMeasures.clamp_weight(w::Interval{T}, maxw::T = one(T)) where {T}
diff --git a/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl b/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl
index 14a48ab..ac42b98 100644
--- a/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl
+++ b/lib/DiscreteMeasures/test/Core/discrete_measures_tests.jl
@@ -91,11 +91,12 @@ end
     end
 
     @testset "Interval Ext" begin
-        @test clamp_domain(interval(-2, -1), 0, 2) == emptyinterval()
-        @test clamp_domain(interval(-2, 1), 0, 2) == interval(0, 1)
-        @test clamp_domain(interval(0.1, 0.5), 0, 2) == interval(0.1, 0.5)
-        @test clamp_domain(interval(0.5, 5), 0, 2) == interval(0.5, 2)
-        @test clamp_domain(interval(3, 4), 0, 2) == emptyinterval()
+        # IntervalArithmetic 1.x: empty intervals compare false under `==` even when equal.
+        @test isequal_interval(clamp_domain(interval(-2, -1), 0, 2), emptyinterval())
+        @test isequal_interval(clamp_domain(interval(-2, 1), 0, 2), interval(0, 1))
+        @test isequal_interval(clamp_domain(interval(0.1, 0.5), 0, 2), interval(0.1, 0.5))
+        @test isequal_interval(clamp_domain(interval(0.5, 5), 0, 2), interval(0.5, 2))
+        @test isequal_interval(clamp_domain(interval(3, 4), 0, 2), emptyinterval())
     end
 
 end
diff --git a/lib/OUQBase/src/interface.jl b/lib/OUQBase/src/interface.jl
index 831dd9c..58b9ba2 100644
--- a/lib/OUQBase/src/interface.jl
+++ b/lib/OUQBase/src/interface.jl
@@ -203,6 +203,18 @@ function get_group_name(var, ouq_sys::OUQSystem)
     return get_group_name(var, ouq_sys.admissible_set)
 end
 
+# `𝔼`/`ℙ` are modeled as `Symbolics.Operator`s (the same abstraction used for
+# e.g. `Differential`). `SymbolicUtils.default_is_atomic` therefore treats any
+# `𝔼(Q)`/`ℙ(Q)` application as an atomic leaf (correct for `D(x)`, but not for
+# us): we need to recurse through it to discover the underlying random
+# variable `Q` that it wraps.
+function _random_variable_is_atomic(ex)
+    if SymbolicUtils.iscall(ex) && SymbolicUtils.operation(ex) isa Operator
+        return false
+    end
+    return SymbolicUtils.default_is_atomic(ex)
+end
+
 function get_ordered_group_names(
         expression,
         admissible_set::AdmissibleSet;
@@ -210,11 +222,11 @@ function get_ordered_group_names(
         ensure_all = false,
     )
     if ensure_singleton
-        vars = get_variables(expression)
+        vars = get_variables(expression; is_atomic = _random_variable_is_atomic)
         @assert length(vars) == 1 "Expression $expression has multiple random variables $vars. Disable `ensure_singleton` if this is intended."
         return unique([get_group_name(only(vars), admissible_set)])
     else
-        vars = get_variables(expression)
+        vars = get_variables(expression; is_atomic = _random_variable_is_atomic)
         @debug "Expression $expression has multiple random variables $vars"
         unordered_group_names = Set([get_group_name(var, admissible_set) for var in vars])
         if ensure_all
diff --git a/lib/OUQBase/src/operators.jl b/lib/OUQBase/src/operators.jl
index adb9562..4afb743 100644
--- a/lib/OUQBase/src/operators.jl
+++ b/lib/OUQBase/src/operators.jl
@@ -1,14 +1,20 @@
 import Symbolics: Operator, value
-import SymbolicUtils: Term, symtype
+import SymbolicUtils: symtype, term
 
 struct 𝔼_ <: Operator end
 
 const 𝔼 = 𝔼_() # To get same object
 
-(::𝔼_)(x) = Term{symtype(x)}(𝔼, Any[x])
+# `term` (rather than constructing `Term{T}` directly) picks the SymVariant
+# type parameter itself and infers the symbolic type via `promote_symtype`.
+(::𝔼_)(x) = term(𝔼, x)
 (::𝔼_)(x::Num) = Num(𝔼(value(x)))
 
 SymbolicUtils.promote_symtype(::𝔼_, x) = x
+# Disambiguates against SymbolicUtils's generic `promote_symtype(::Operator,
+# ::Type{T}) where T`, which is equally specific to the untyped method above
+# when called with a `Type` argument (e.g. `symtype(x) === Real`).
+SymbolicUtils.promote_symtype(::𝔼_, ::Type{T}) where {T} = T
 Base.nameof(::𝔼_) = :𝔼
 SymbolicUtils.isbinop(::𝔼_) = false
 
@@ -28,10 +34,11 @@ struct ℙ_ <: Operator end
 
 const ℙ = ℙ_() # To get same object
 
-(::ℙ_)(x) = Term{symtype(x)}(ℙ, Any[x])
+(::ℙ_)(x) = term(ℙ, x)
 (::ℙ_)(x::Num) = Num(ℙ(value(x)))
 
 SymbolicUtils.promote_symtype(::ℙ_, x) = x
+SymbolicUtils.promote_symtype(::ℙ_, ::Type{T}) where {T} = T
 Base.nameof(::ℙ_) = :ℙ
 SymbolicUtils.isbinop(::ℙ_) = false
 
diff --git a/src/OptimalUncertaintyQuantification.jl b/src/OptimalUncertaintyQuantification.jl
index bc5b3ad..58a8d61 100644
--- a/src/OptimalUncertaintyQuantification.jl
+++ b/src/OptimalUncertaintyQuantification.jl
@@ -3,6 +3,11 @@ module OptimalUncertaintyQuantification
 using Reexport: @reexport
 @reexport using OUQBase
 
+# `@random_variables` expands to code that references `Symbolics.@variables`
+# in the caller's scope (via `esc`), so `Symbolics` must be bound here for the
+# precompile workload below to expand successfully.
+using Symbolics: Symbolics
+
 using PrecompileTools: @compile_workload, @setup_workload
 
 @setup_workload begin
@@ -11,7 +16,7 @@ using PrecompileTools: @compile_workload, @setup_workload
             Independent(Q, bounds = (0.0, 1.0))
         end
         admissible_set = AdmissibleSet(random_vars, [𝔼(Q) ~ 0.5])
-        OUQSystem(
+        OUQSystem(;
             objective = 𝔼(Q),
             admissible_set,
             reduction_alg = WinklerExtremalMeasures(),