GSA using Morris and Sobol's methods gives output on different forms

[TorkelE]

I realise there might be a reason for this given that Sobol's S2 field is a matrix, but this is till a bit odd:

using Catalyst
seir_model = @reaction_network begin
    10^β, S + I --> E + I
    10^a, E --> I
    10^γ, I --> R
end

using OrdinaryDiffEq
u0 = [:S => 999.0, :I => 1.0, :E => 0.0, :R => 0.0]
function peak_cases(p)
    oprob = ODEProblem(seir_model, u0, (0.0, 10000.0), p)
    sol = solve(oprob, Tsit5(); maxiters=100000, verbose=false)
    SciMLBase.successful_retcode(sol) || return Inf
    return maximum(sol[:I])
end

using GlobalSensitivity
g_sens_1 = gsa(peak_cases, Morris(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)])
g_sens_2 = gsa(peak_cases, Sobol(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)]; samples=500)

Here, g_sens_1.means_star is a 1x3 Matrix, while g_sens_2.S1 is a length-3 vector. If I want to e.g. plot these using bar charts I have to use different notations. Is this intended?

Next, if I want to consider a vector output:

function peak_cases(p)
    oprob = ODEProblem(seir_model, u0, (0.0, 10000.0), p)
    sol = solve(oprob, Tsit5(); maxiters=100000, verbose=false)
    SciMLBase.successful_retcode(sol) || return Inf
    return [maximum(sol[:I]), maximum(sol[:E])]
end

using GlobalSensitivity
g_sens_1 = gsa(peak_cases, Morris(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)])
g_sens_2 = gsa(peak_cases, Sobol(), [(-3.0,-1.0), (-2.0,0.0), (-2.0,0.0)]; samples=500)

Suddenly we have both g_sens_1.means_starand g_sens_2.S1 being 1x3 matrices. Also, g_sens_2.S2 is entirely absent. This all feels inconsistent.

[Spinachboul]

@TorkelE I realize that for uniform plotting, we would need to normalize the input values for a consistent output. For this we can write a simle function that transforms both the output formats 1x3 Matrix and length-3 vector into a unform vector format.

#Function to convert different output formats to uniform Vector format
function convert_to_vector(output)
    if isa(output, Matrix)
        return vec(output)
        
    else if isa(output, AbstractVector)
        return output
    
    else
        return Vector{FLoat64}()

Now use this function to convert the outputs:

using Plots

# Assuming g_sens_1.means_star and g_sens_2.S1 are your outputs
output_1 = g_sens_1.means_star
output_2 = g_sens_2.S1

# Convert outputs to Vector
output_1_vec = convert_to_vector(output_1)
output_2_vec = convert_to_vector(output_2)

# Plotting using Plots.jl
bar(["Parameter 1", "Parameter 2", "Parameter 3"], [output_1_vec, output_2_vec], 
    label = ["Means Star 1" "S1 2"], title = "Sensitivity Analysis Results", 
    ylabel = "Sensitivity Index")

Please test this code and check if this solution works!