diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d5af33039..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" diff --git a/Project.toml b/Project.toml index da0f383e3..7b488d56a 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 ", - "Ander Gray ", -] - +authors = ["Jasper Behrensdorf ", "Ander Gray "] [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" 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 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..bc7c5ffa0 100644 --- a/src/inputs/imprecise/interval.jl +++ b/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 \ No newline at end of file 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/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 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 diff --git a/src/simulations/montecarlo.jl b/src/simulations/montecarlo.jl index 439a02a34..5c0a6ff1c 100644 --- a/src/simulations/montecarlo.jl +++ b/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 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}) && 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 54d202a10..8f458bd80 100644 --- a/test/inputs/imprecise/interval.jl +++ b/test/inputs/imprecise/interval.jl @@ -1,59 +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 +@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/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 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") 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