Add JointInterval

.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"

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,15 @@ 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"
 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 +50,14 @@ 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"
 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"

docs/src/api/inputs.md

@@ -16,6 +16,7 @@ RandomVariable
 IntervalVariable
 EmpiricalDistribution
 JointDistribution
+JointInterval
 ```
 
 ## Functions

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

src/inputs/imprecise/interval.jl

@@ -107,3 +107,28 @@ 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}
+    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([x[:, i] for i in axes(x, 2)]))
+        return new(intervals, hull)
+    end
+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
+to_physical_space!(_::JointInterval, _::DataFrame) = nothing

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))

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

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

src/simulations/montecarlo.jl

@@ -74,21 +74,82 @@ 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)
+"""
+    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,
+    n_internal::Integer=sim.n * 10,
+)
+    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 to be applied to 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)(n_internal, 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
+
+    # 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
 
-    samples = quantile.(Normal(), u)
-    samples = DataFrame(names(random_inputs) .=> eachrow(samples))
+    # 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])
+        else
+            @warn "Only $(size(samples)[1]) of $(sim.n) samples generated. Try increasing 'n_internal'"
+        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

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}) &&

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"