I am running a model (written in Julia) with 42 parameters and 6 outputs (cost values between modelled variables and observed variables). If I run the gsa function with N=10000 (or larger number), it returns ST, S1 and S2 values with very large values like 1e+16...
Here are the S1 values I got...42 is the num of parameter and 6 is the num of outputs.
You could see that the 13th row (parameter epsilon_max) and the last two rows have very large S1 values...But if I narrow down the parameter group to a small group (still including 13th row), then there is no 1e+16 values for that parameter (epsilon_max)...
For MWE, the model is large and not public yet...but I can show you the codes about GSA as below...for cost function I am using inversed normalised Nash Efficiency like below:
function metric(y::AbstractArray, yσ::AbstractArray, ŷ::AbstractArray, ::NNSE)
NSE_v = metric(y, yσ, ŷ, NSE())
NNSE = one(eltype(ŷ)) / (one(eltype(ŷ)) + one(eltype(ŷ)) - NSE_v)
return NNSE
end
function metric(y::AbstractArray, yσ::AbstractArray, ŷ::AbstractArray, ::NNSEInv)
NNSEInv = one(eltype(ŷ)) - metric(y, yσ, ŷ, NNSE())
return NNSEInv
end
function metric(y::AbstractArray, yσ::AbstractArray, ŷ::AbstractArray, ::NSE)
NSE = one(eltype(ŷ)) .- sum(abs2.((y .- ŷ))) / sum(abs2.((y .- mean(y))))
return NSE
end
function globalSensitivity(cost_function, method_options, p_bounds, ::GlobalSensitivitySobolDM; batch=true)
sampler = getproperty(SindbadOptimization.GlobalSensitivity, Symbol(method_options.sampler))(; method_options.sampler_options...);
samples = method_options.samples;
lb = first.(p_bounds);
ub = last.(p_bounds);
@debug samples
A, B = QuasiMonteCarlo.generate_design_matrices(samples, lb, ub, sampler);
@debug size(A)
@debug method_options
results = gsa(cost_function, Sobol(; method_options.method_options...), A, B; method_options..., batch=batch);
return results
end
function costSensitivityMatrix(param_matrix, _, selected_models, space_forcing, space_spinup_forcing, loc_forcing_t,
output_array, space_output, space_land, tem_info, observations, param_updater, cost_options,
multi_constraint_method, parameter_scaling_type, cost_out::Matrix, ::CostModelObsMT)
@info size(param_matrix)
param_set_size = size(param_matrix, 2)
done_params=1
Threads.@threads for param_index in eachindex(1:param_set_size)
idx = Threads.threadid()
param_vector = param_matrix[:, param_index]
@debug param_vector
updated_models = updateModels(param_vector, param_updater, parameter_scaling_type, selected_models)
coreTEM!(updated_models, space_forcing, space_spinup_forcing, loc_forcing_t, space_output[idx], space_land, tem_info)
cost_vector = metricVector(space_output[idx], observations, cost_options)
cost_metric = combineMetric(cost_vector, multi_constraint_method)
# cost_out[param_index] = cost_metric
@debug length(cost_vector)
cost_out[:, param_index] = cost_vector
@debug cost_out
@info idx, round(100 * done_params/param_set_size,digits=2), param_set_size, cost_metric
@debug idx, round(100 * done_params/param_set_size,digits=2), param_set_size, cost_metric, cost_vector
done_params += 1
end
return cost_out
end
function prepOptiSensitivityMatrix(forcing, observations, info, ::CostModelObsMT; algorithm_info_field=:algorithm_optimization)
algorithm_info = getproperty(info.optimization, algorithm_info_field);
opti_helpers = prepOpti(forcing, observations, info, CostModelObs());
run_helpers = opti_helpers.run_helpers;
cost_vector_size = getCostVectorSize(getproperty(algorithm_info, :options), opti_helpers.default_values, getproperty(algorithm_info, :method));
cost_vector = Vector{eltype(opti_helpers.default_values)}(undef, cost_vector_size);
num_observations = length(info.optimization.observational_constraints);
cost_matrix = Matrix{eltype(opti_helpers.default_values)}(undef, num_observations, cost_vector_size);
space_index = 1 # the parallelization of cost computation only runs in single pixel runs
cost_function = x -> costSensitivityMatrix(x, opti_helpers.default_values, info.models.forward,
run_helpers.space_forcing[space_index],
run_helpers.space_spinup_forcing[space_index],
run_helpers.loc_forcing_t, run_helpers.output_array,
run_helpers.space_output_mt,
deepcopy(run_helpers.space_land[space_index]),
run_helpers.tem_info, observations, opti_helpers.tbl_params,
opti_helpers.cost_options, info.optimization.multi_constraint_method,
info.optimization.optimization_parameter_scaling, cost_matrix,
info.optimization.optimization_cost_method)
@debug cost_function(opti_helpers.default_values)
opti_helpers = (; opti_helpers..., cost_function=cost_function, cost_vector=cost_matrix)
return opti_helpers
end
function getCostVectorSize(algo_options, param_vector, ::GlobalSensitivitySobol)
default_opt = sindbad_default_options(GlobalSensitivitySobol())
samples = default_opt.samples
nparam = length(param_vector)
norder = length(algo_options.method_options.order) - 1
if hasproperty(algo_options, :samples)
samples = algo_options.samples
end
cost_vector_size = samples * (norder * nparam + 2)
return cost_vector_size
end
**Environment (please complete the following information):**
- Output of `using Pkg; Pkg.status()`
```Julia
julia> using Pkg; Pkg.status()
Status `/Net/Groups/BGI/work_4/scratch/eebiomass/sindbad.jl.github/examples/exp_WROASTED/Project.toml`
[af5da776] GlobalSensitivity v2.7.0
[033835bb] JLD2 v0.5.12
[91a5bcdd] Plots v1.40.11
[8a4e6c94] QuasiMonteCarlo v0.3.3
[295af30f] Revise v3.7.3
[6686e6de] Sindbad v0.1.0 `../..`
[772b809f] SindbadData v0.1.0 `../../lib/SindbadData`
[bbd31fc6] SindbadExperiment v0.1.0 `../../lib/SindbadExperiment`
[eb1d8004] SindbadML v0.1.0 `../../lib/SindbadML`
[7a8077d0] SindbadMetrics v0.1.0 `../../lib/SindbadMetrics`
[64fff8bd] SindbadOptimization v0.1.0 `../../lib/SindbadOptimization`
[2b7f2987] SindbadSetup v0.1.0 `../../lib/SindbadSetup`
[f6108451] SindbadTEM v0.1.0 `../../lib/SindbadTEM`
[a2ad09c8] SindbadUtils v0.1.0 `../../lib/SindbadUtils`
[860ef19b] StableRNGs v1.0.2
I am running a model (written in Julia) with 42 parameters and 6 outputs (cost values between modelled variables and observed variables). If I run the
gsafunction withN=10000(or larger number), it returns ST, S1 and S2 values with very large values like 1e+16...Here are the S1 values I got...42 is the num of parameter and 6 is the num of outputs.
You could see that the 13th row (parameter epsilon_max) and the last two rows have very large S1 values...But if I narrow down the parameter group to a small group (still including 13th row), then there is no 1e+16 values for that parameter (epsilon_max)...
For MWE, the model is large and not public yet...but I can show you the codes about GSA as below...for cost function I am using inversed normalised Nash Efficiency like below:
The sensitivity code is
using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()