Add multiple outputs for the evaluation of Sobol indices (#98)

Co-authored-by: andrea perin <andrea.perin@irz.uni-hannover.de>
Co-authored-by: Jasper Behrensdorf <jasper@behrensdorf.de>

Author: andreaperin

demo/sensitivity/sobolindices.jl

@@ -14,7 +14,7 @@ B = Model(
 
 mc = MonteCarlo(10000)
 
-s_mc = sobolindices([B], [X; ω], :B, mc)
+s_mc = sobolindices(B, [X; ω], :B, mc)
 
 # Compare with analytical solution
 VB = sum(σx .^ 2 .* σω .^ 2)

src/sensitivity/sobolindices.jl

@@ -1,7 +1,12 @@
+const _sobol_table_types = [Symbol[], Float64[], Float64[], Float64[], Float64[]]
+const _sobol_table_header = [
+    "Variables", "FirstOrder", "FirstOrderStdError", "TotalEffect", "TotalEffectStdError"
+]
+
 function sobolindices(
     models::Vector{<:UQModel},
     inputs::Vector{<:UQInput},
-    output::Symbol,
+    outputs::Vector{Symbol},
     sim::AbstractMonteCarlo,
 )
     sim_double_samples = @set sim.n = 2 * sim.n
@@ -12,20 +17,19 @@ function sobolindices(
     random_names = names(filter(i -> isa(i, RandomUQInput), inputs))
 
     evaluate!(models, samples)
+    indices = Dict([
+        (name, DataFrame(_sobol_table_types, _sobol_table_header)) for name in outputs
+    ])
 
     A = samples[1:(sim.n), :]
     B = samples[(sim.n + 1):end, :]
 
-    fA = A[:, output]
-    fB = B[:, output]
-
-    fA .-= mean(fA)
-    fB .-= mean(fB)
-
-    VY = var([fA; fB])
+    fA = Matrix(A[:, outputs])
+    fB = Matrix(B[:, outputs])
+    fA = fA .- mean(fA; dims=1)
+    fB = fB .- mean(fB; dims=1)
 
-    Si = zeros(length(random_names), 2)
-    STi = zeros(length(random_names), 2)
+    VY = var([fA; fB]; dims=1)
 
     for (i, name) in enumerate(random_names)
         ABi = select(A, Not(name))
@@ -35,30 +39,25 @@ function sobolindices(
             evaluate!(m, ABi)
         end
 
-        ABi[:, output] .-= mean(ABi[:, output])
+        for (j, qty) in enumerate(outputs)
+            ABi[:, qty] .-= mean(ABi[:, qty])
 
-        first_order = x -> mean(fB .* (x .- fA)) / VY # Saltelli 2009
-        total_effect = x -> (1 / (2 * sim.n)) * sum((fA .- x) .^ 2) / VY # Saltelli 2009
+            # First order effects
+            first_order = x -> mean(fB[:, j] .* (x .- fA[:, j])) / VY[j] # Saltelli 2009
+            bs = bootstrap(first_order, ABi[:, qty], BasicSampling(1000))
+            Sᵢ = first_order(ABi[:, qty])
+            σSᵢ = stderror(bs)[1]
 
-        # First order effects
-        Si[i, 1] = first_order(ABi[:, output])
-        bs = bootstrap(first_order, ABi[:, output], BasicSampling(1000))
-        Si[i, 2] = stderror(bs)[1]
+            # Total effects
+            total_effect = x -> (1 / (2 * sim.n)) * sum((fA[:, j] .- x) .^ 2) / VY[j] # Saltelli 2009
+            bs = bootstrap(total_effect, ABi[:, qty], BasicSampling(1000))
+            Sₜ = total_effect(ABi[:, qty])
+            σSₜ = stderror(bs)[1]
 
-        # Total effects
-        STi[i, 1] = total_effect(ABi[:, output])
-        bs = bootstrap(total_effect, ABi[:, output], BasicSampling(1000))
-        STi[i, 2] = stderror(bs)[1]
+            push!(indices[qty], [name, Sᵢ, σSᵢ, Sₜ, σSₜ])
+        end
     end
-
-    indices = DataFrame()
-    indices.Variables = random_names
-    indices.FirstOrder = Si[:, 1]
-    indices.FirstOrderStdError = Si[:, 2]
-    indices.TotalEffect = STi[:, 1]
-    indices.TotalEffectStdError = STi[:, 2]
-
-    return indices
+    return length(outputs) > 1 ? indices : indices[outputs[1]]
 end
 
 function sobolindices(pce::PolynomialChaosExpansion)
@@ -78,3 +77,48 @@ function sobolindices(pce::PolynomialChaosExpansion)
 
     return indices
 end
+
+function sobolindices(
+    models::Vector{<:UQModel},
+    inputs::I where {I<:UQInput},
+    outputs::Vector{Symbol},
+    sim::AbstractMonteCarlo,
+)
+    return sobolindices(models, [inputs], outputs, sim)
+end
+
+function sobolindices(
+    models::Vector{<:UQModel},
+    inputs::Vector{<:UQInput},
+    outputs::Symbol,
+    sim::AbstractMonteCarlo,
+)
+    return sobolindices(models, inputs, [outputs], sim)
+end
+
+function sobolindices(
+    models::M where {M<:UQModel},
+    inputs::Vector{<:UQInput},
+    outputs::Symbol,
+    sim::AbstractMonteCarlo,
+)
+    return sobolindices([models], inputs, [outputs], sim)
+end
+
+function sobolindices(
+    models::Vector{<:UQModel},
+    inputs::I where {I<:UQInput},
+    outputs::Symbol,
+    sim::AbstractMonteCarlo,
+)
+    return sobolindices(models, [inputs], [outputs], sim)
+end
+
+function sobolindices(
+    models::M where {M<:UQModel},
+    inputs::I where {I<:UQInput},
+    outputs::Symbol,
+    sim::AbstractMonteCarlo,
+)
+    return sobolindices([models], [inputs], [outputs], sim)
+end

test/sensitivity/sobolindices.jl

@@ -1,69 +1,88 @@
 @testset "sobolindices" begin
     x = RandomVariable.(Uniform(-π, π), [:x1, :x2, :x3])
 
-    ishigami = Model(
-        df -> sin.(df.x1) + 7 .* sin.(df.x2) .^ 2 + 0.1 .* df.x3 .^ 4 .* sin.(df.x1), :f
+    ishigami1 = Model(
+        df -> sin.(df.x1) + 7 .* sin.(df.x2) .^ 2 + 0.1 .* df.x3 .^ 4 .* sin.(df.x1), :f1
     )
+    ishigami2 = Model(
+        df -> sin.(df.x1) + 7 .* sin.(df.x2) .^ 2 + 0.05 .* df.x3 .^ 4 .* sin.(df.x1), :f2
+    )
+
+    n_mc = 4000
+    n_qmc = 2300
 
-    n = 1000
+    firstorder_analytical1 = [0.3138, 0.4424, 0.00]
+    totaleffect_analytical1 = [0.5574, 0.4424, 0.2436]
 
-    firstorder_analytical = [0.3138, 0.4424, 0.00]
-    totaleffect_analytical = [0.5574, 0.4424, 0.2436]
+    firstorder_analytical2 = [0.219, 0.687, 0.00]
+    totaleffect_analytical2 = [0.3136, 0.687, 0.0946]
 
     @testset "Monte Carlo" begin
         Random.seed!(8128)
 
-        si = sobolindices([ishigami], x, :f, MonteCarlo(n))
+        si = sobolindices(ishigami1, x, :f1, MonteCarlo(n_mc))
 
         Random.seed!()
 
-        @test all(isapprox(si.FirstOrder, firstorder_analytical; rtol=0.1))
-        @test all(isapprox(si.TotalEffect, totaleffect_analytical; rtol=0.1))
+        @test all(isapprox(si.FirstOrder, firstorder_analytical1; rtol=0.1))
+        @test all(isapprox(si.TotalEffect, totaleffect_analytical1; rtol=0.1))
     end
 
     @testset "Sobol" begin
-        Random.seed!(8128)
+        si = sobolindices(ishigami1, x, :f1, SobolSampling(n_qmc))
 
-        si = sobolindices([ishigami], x, :f, SobolSampling(n))
+        @test all(isapprox(si.FirstOrder, firstorder_analytical1; rtol=0.1))
+        @test all(isapprox(si.TotalEffect, totaleffect_analytical1; rtol=0.1))
+    end
 
-        Random.seed!()
+    @testset "Sobol Nultiple Outputs" begin
+        si = sobolindices([ishigami1, ishigami2], x, [:f1, :f2], SobolSampling(n_qmc))
 
-        @test all(isapprox(si.FirstOrder, firstorder_analytical; rtol=0.1))
-        @test all(isapprox(si.TotalEffect, totaleffect_analytical; rtol=0.1))
+        @test all(isapprox(si[:f1].FirstOrder, firstorder_analytical1; rtol=0.1))
+        @test all(isapprox(si[:f1].TotalEffect, totaleffect_analytical1; rtol=0.1))
+        @test all(isapprox(si[:f2].FirstOrder, firstorder_analytical2; rtol=0.1))
+        @test all(isapprox(si[:f2].TotalEffect, totaleffect_analytical2; rtol=0.1))
     end
 
     @testset "Halton" begin
-        Random.seed!(8128)
-
-        si = sobolindices([ishigami], x, :f, MonteCarlo(n))
+        si = sobolindices(ishigami1, x, :f1, HaltonSampling(n_qmc))
 
-        Random.seed!()
-
-        @test all(isapprox(si.FirstOrder, firstorder_analytical; rtol=0.1))
-        @test all(isapprox(si.TotalEffect, totaleffect_analytical; rtol=0.1))
+        @test all(isapprox(si.FirstOrder, firstorder_analytical1; rtol=0.1))
+        @test all(isapprox(si.TotalEffect, totaleffect_analytical1; rtol=0.1))
     end
 
     @testset "Latin Hypercube" begin
-        Random.seed!(8128)
+        si = sobolindices(ishigami1, x, :f1, LatinHypercubeSampling(n_qmc))
 
-        si = sobolindices([ishigami], x, :f, MonteCarlo(n))
-
-        Random.seed!()
+        @test all(isapprox(si.FirstOrder, firstorder_analytical1; rtol=0.1))
+        @test all(isapprox(si.TotalEffect, totaleffect_analytical1; rtol=0.1))
+    end
 
-        @test all(isapprox(si.FirstOrder, firstorder_analytical; rtol=0.1))
-        @test all(isapprox(si.TotalEffect, totaleffect_analytical; rtol=0.1))
+    @testset "Convenience Functions" begin
+        x1 = RandomVariable(Uniform(-π, +π), :x1)
+        x2 = RandomVariable(Uniform(0, +π), :x2)
+        t1 = Model(df -> sin.(df.x1), :t1)
+        t2 = Model(df -> cos.(df.x1), :t2)
+
+        @test isa(sobolindices([t1; t2], x1, :t1, SobolSampling(2)), DataFrame)
+        @test isa(sobolindices(t1, [x1; x2], :t1, SobolSampling(2)), DataFrame)
+        @test isa(sobolindices([t1; t2], [x1; x2], :t1, SobolSampling(2)), DataFrame)
+        @test isa(
+            sobolindices([t1; t2], x1, [:t1, :t2], SobolSampling(2)), Dict{Symbol,DataFrame}
+        )
+        @test isa(sobolindices(t1, x1, :t1, SobolSampling(2)), DataFrame)
     end
 
     @testset "Polynomial Chaos Expansion" begin
         p = 6
         Ψ = PolynomialChaosBasis(LegendreBasis.([p, p, p]), p)
 
         gq = GaussQuadrature()
-        pce, _ = polynomialchaos(x, ishigami, Ψ, :f, gq)
+        pce, _ = polynomialchaos(x, ishigami1, Ψ, :f1, gq)
 
         si = sobolindices(pce)
 
-        @test all(isapprox(si.FirstOrder, firstorder_analytical; rtol=0.1))
-        @test all(isapprox(si.TotalEffect, totaleffect_analytical; rtol=0.1))
+        @test all(isapprox(si.FirstOrder, firstorder_analytical1; rtol=0.1))
+        @test all(isapprox(si.TotalEffect, totaleffect_analytical1; rtol=0.1))
     end
-end
+end