eFAST yields non-zero indices when parameter bounds are identical

[Ge0rges]

Hello,

I have a model which I seek to conduct a sensitivity analysis of. To do so efficiently I am trying to use eFAST. I have normalized the model outputs to be between 0 and 1 in all cases by simple scaling. However, eFAST yields non-zero sobol indices even when the parameter range passed has an identical upper and lower bound. I have confirmed that these are equal using the same loop and if statement used in the Julia code.

I would like to understand how each parameter affects the final output of the model, which is a vector of 4 states. I would also like to understand what the confidence level is for each of those results.

MWE:

using DifferentialEquations
using DiffEqSensitivity
using ForwardDiff
using GlobalSensitivity
using Statistics
using Printf

# Runs the sensitivity analysis using Sobol method.
function run_sensitivity_analysis()
    p_bounds = [[1.0e-6, 100.0], [1.0e-5, 500], [1, 500], [1000000000, 1000000000], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 100], [1000.0, 1.0e8], [-1, -1]]
    u0 = [1.6382099999999998e13, 3.4061512877350002e10, 0, 100000, 1.461e7]

    # Define a function that remakes the problem and gets its result. Called for each sample.
    f1 = function (p)
        sol = solve_model(p, u0)

        scaled1 = sol[1, end]/1e14
        scaled2 = sol[2, end]/1e14
        scaled3 = sol[3, end]/1e19
        scaled4 = sol[4, end]/1e9

        if scaled1 < 0 || scaled1 > 1 || scaled2 < 0 || scaled2 > 1 || scaled3 < 0 || scaled3 > 1 || scaled4 < 0 || scaled4 > 1
            @show [scaled1, scaled2, scaled3, scaled4]
        end
        [scaled1, scaled2, scaled3, scaled4]
    end

    # Run sobol - This takes over 30GBs of memory and gets killed
#     sobol_result = GlobalSensitivity.gsa(f1, Sobol(order=[0, 1], nboot=2 conf_level=0.95), p_bounds, samples=2^13)

    sobol_result = GlobalSensitivity.gsa(f1, eFAST(), p_bounds, samples=2^13)

    @show sobol_result.ST[1,:]
    return (sobol_result.ST[1,:], sobol_result.S1[1,:])
end


# Defines a problem object and solves it.
function solve_model(p, u0)
    # Build a callback to introduce a carbon addition as it is a discontinuity
    stops = []
    carbon = []
    punctual_organic_carbon_addition = p[end] == -1 ? [] : p[end]
    for (time_to_add, (pOC_to_add, dOC_to_add)) in punctual_organic_carbon_addition
        push!(stops, time_to_add)
        push!(carbon, (convert(Float64, pOC_to_add), convert(Float64, dOC_to_add)))
    end

    # Convert p and u0, and remove puncutal_punctual_organic_carbon_addition from p
    p = convert(Array{Float64}, p[1:end-1])
    u0 = convert(Array{Float64}, u0)

    # Callback for punctual addition - causes a discontinuity
    additions = Dict(Pair.(stops, carbon))
    function addition!(integrator)
        integrator.u[1] += additions[integrator.t][1]
        integrator.u[2] += additions[integrator.t][2]

        if integrator.u[4] < 1
            integrator.u[4] = 1  # Add a viable cell if none exist
        end
    end
    carbon_add_cb = PresetTimeCallback(stops, addition!)

    # Callback for the max() function in the model - causes a discontinuity.
    # If max equation changes, this condition will have to change.
    max_condition(u, t, integrator) = integrator.p[2] * u[4] - u[2]
    do_nothing(integrator) = nothing
    max_cb = ContinuousCallback(max_condition, do_nothing)

    # Min condition for EEA rate
    min_condition(u, t, integrator) = integrator.p[10] * u[4] - u[1]
    min_cb = ContinuousCallback(min_condition, do_nothing)

    # Set things to 0 if they are less than 1e-100
    zero_condition(u, t, integrator) = u[1] < 1e-20 || u[2] < 1e-20 || u[3] < 1e-20 || u[4] < 1e-20
    function zero_out!(integrator)
        for i in 1:4
            if integrator.u[i] < 1e-20
                integrator.u[i] = 0
            end
        end
    end
    zero_cb = DiscreteCallback(zero_condition, zero_out!)

    # Callback list
    cbs = CallbackSet(carbon_add_cb, max_cb, min_cb, zero_cb)

    # Out of domain function
    is_invalid_domain(u,p,t) = u[1] < 0 || u[2] < 0 || u[3] < 0 || u[4] < 0

    # Build the ODE Problem and Solve
    prob = ODEProblem(model, u0[1:end-1], (0.0, last(u0)), p)
    sol = solve(prob, Rosenbrock23(), callback=cbs, isoutofdomain=is_invalid_domain, maxiters=1e6)

    return sol
end


# Implements the differential equations that define the model
function model(du, u, p, t)
    # Paramaters
    mu_max = p[1]
    maintenance_per_cell = p[2]
    dOC_per_cell = p[3]

    carrying_capacity = p[4]
    pOC_input_rate = p[5]
    dOC_input_rate = p[6]
    inorganic_carbon_input_rate = p[7]
    inorganic_carbon_fixing_rate = p[8]
    inorganic_carbon_per_cell = p[9]
    eea_rate = p[10]
    Kd = p[11]

    # Load state conditions
    pOC_content = u[1]
    dOC_content = u[2]
    inorganic_carbon_content = u[3]
    cell_count = u[4]

    ## CELL COUNT
    # Growth
    growth = mu_max * (dOC_content / (Kd + dOC_content)) * cell_count * (1 - cell_count / carrying_capacity)

    # Organic carbon requirement
    required_dOC_per_cell = maintenance_per_cell
    required_dOC = required_dOC_per_cell * cell_count

    # Starvation Deaths
    dOC_missing = max(required_dOC - dOC_content, 0)
    starvation_deaths = dOC_missing / required_dOC_per_cell

    # Total Deaths
    deaths = starvation_deaths

    ## CARBON
    dOC_consumption = required_dOC_per_cell * (cell_count - deaths)
    fixed_carbon = inorganic_carbon_fixing_rate * inorganic_carbon_content

    # EEA rate
    eea_removal = min(eea_rate*cell_count, pOC_content)

    # Particulate Organic carbon
    du[1] = pOC_input_rate - eea_removal

    # Dissolved Organic carbon
    du[2] = dOC_input_rate + dOC_per_cell * (deaths - growth) - dOC_consumption + eea_removal + fixed_carbon

    # Inorganic carbon
    du[3] = inorganic_carbon_input_rate + inorganic_carbon_per_cell * (deaths - growth) + required_dOC - fixed_carbon

    # Net cell count change
    du[4] = growth - deaths
end