Analysis

f_analysis

Code

function [rst,stg] = f_analysis(stg,analysis,mmf,analysis_options,rst)
% This function performs a specific analysis based on the given analysis
% options. The function calls the appropriate analysis functions
% (f_diagnostics, f_opt, f_sim_score, f_lsa, f_gsa, f_PL_m) depending on
% the selected analysis type. The analysis results are stored in an updated
% 'rst' structure, and the settings are updated as needed.
%
% Inputs:
% - stg: Settings to use in the analysis
% - analysis: Type of analysis to be performed
% - mmf: Main model files
% - analysis_options: List of possible analysis options
% - rst: Structure for storing analysis results
%
% Outputs:
% - rst: Updated structure containing the results of the analysis
% - stg: Updated settings after the analysis
%
% Used Functions:
% - f_diagnostics: Performs diagnostic analysis
% - f_opt: Performs optimization analysis
% - f_sim_score: Runs simulation for a given option
% - f_lsa: Performs Local sensititivity Analysis (LSA)
% - f_gsa: Performs Global ensititivity Analysis (GSA)
% - f_PL_m: Performs Profile Likelihood Analysisa (PLA)
%
% Variables:
% Initialized:
% - matching_option_idx: Index of the selected analysis option
% - sim_results: Cell array for storing simulation results
% - valid_options: Array of valid optimization options
%
% Persistent: None
% analysis
% analysis_options
% matching_option_idx
matching_option_idx = find(contains(analysis_options, analysis), 1);

% display start message
disp("Starting " + {analysis_options(matching_option_idx)})

% Perform the analysis based on the current analysis option index
switch matching_option_idx
    case 1
        % run diagnostic analysis
        rst.diag = f_diagnostics(stg,mmf);
    case 2
        % run optimization analysis
        [rst.opt,pa] = f_opt(stg,mmf);
        % find valid optimization options
        valid_options = find(pa(:,1) ~= 0);
        % update settings with valid options
        stg.pat = transpose(valid_options);
        % create cell array for simulation results
        sim_results = cell(length(valid_options),1);
        % parallelize simulation runs
        for j = stg.pat
            % run simulation for current valid option
            [~,sim_results{j},~] = f_sim_score(pa(j,:),stg,mmf);
        end
        % store simulation results in results structure
        rst.diag(stg.pat) = [sim_results{:}];
    case 3
        % run LSA analysis
        rst.lsa = f_lsa(stg,mmf);
    case 4
        % run GSA analysis
        rst.gsa = f_gsa(stg,mmf);
    case 5
        % run PLA analysis
        rst.PLA = f_PL_m(stg,mmf);
end

% display completion message
disp(analysis_options(matching_option_idx) + " completed successfully")
end

The f_analysis function performs a specific analysis based on the given analysis options. It takes settings, analysis type, main model files, a list of possible analysis options, and a structure for storing analysis results as input arguments. The function outputs the results of the analysis and updates the settings as needed. The function calls the appropriate analysis functions (f_diagnostics, f_opt, f_sim_score, f_lsa, f_gsa, f_PL_m) depending on the selected analysis type.

Inputs

• stg: Settings to use in the analysis.

• analysis: Type of analysis to be performed.

• mmf: Main model files.

• analysis_options: List of possible analysis options.

• rst: Structure for storing analysis results.

Outputs

• rst: Updated structure containing the results of the analysis.

• stg: Updated settings after the analysis.

Usage

[rst, stg] = f_analysis(stg, analysis, mmf, analysis_options, rst)

Example

settings = :ref:`stg<stg>`;
analysis_type = 'diagnostics';
model_files = 'path/to/model/files';
options = ["Diagnostics","Parameter Estimation",...
 "Local Sensitivity Analysis","Global Sensitivity Analysis",...
 "Profile Likelihood Analysis","Reproduce a previous analysis",...
 "Reproduce the plots of a previous analysis","Import model files"];
results = :ref:`stg<stg>`;

[results, settings] = f_analysis(settings, analysis_type, model_files, options, results);

Diagnostics

f_diagnostics

Code

function rst = f_diagnostics(stg,mmf)
% Description: This function runs a model for given parameters and
% settings, and computes the scores for fitness, both in multi-core and
% single-core mode. The function uses f_sim_score() to obtain the score for
% each parameter array. The mode (multi-core or single-core) is determined
% by the 'optmc' field in the 'stg' input.
%
% Inputs:
% - stg: A structure containing settings for the model, such as:
%        * pa: An array of parameter arrays to be tested.
%        * pat: Indices of the parameter arrays to be tested.
%        * optmc: A boolean flag, if true the function runs in multi-core
%          mode, otherwise it runs in single-core mode.
% - mmf: A custom memory-mapped file object, used for reading or writing 
% data.
% 
% Outputs:
% - rst: A vector containing the computed scores for each parameter array.
%
% Used Functions:
% - f_sim_score: A function that takes a parameter array, model settings
% (stg), and a custom memory-mapped file object (mmf) as inputs and returns
% the computed score for that parameter array.
%
% Variables:
% Loaded:
% (None)
%
% Initialized:
% - par_array_idx: The current index of the parameter array being tested.
%
% Persistent:
% (None)

% Run the model and obtain scores for fitness Multi Core
if stg.optmc
    disp("Running the model and obtaining Scores (Multicore)")

    pa = stg.pa;
    % Iterate over the parameter arrays to be tested
    parfor par_array_idx = stg.pat
        [~,rst(par_array_idx),~] =...
            f_sim_score(pa(par_array_idx,:),stg,mmf);
    end
    
    % Run the model and obtain scores for fitness single Core
else
    disp("Running the model and obtaining Scores (Singlecore)")
    
    % Iterate over the parameter arrays to be tested
    for par_array_idx = stg.pat
        [~,rst(par_array_idx),~] =...
            f_sim_score(stg.pa(par_array_idx,:),stg,mmf);
    end
end
end

Used to understand the effects of different parameters sets on model behaviour or in comparing different parameters sets.

It loads the user defined configurations, performs all the needed simulations, and calculates scores of the error functions either per experimental output, per experiment, or in total (check results).

• Inputs

  • stg - stg.optmc , stg.pat

• Outputs - rst (diagnostics results)

Optimization

f_opt

Code

function [result,parameter_array] = f_opt(settings,model_folders)
% This function performs optimization on a given objective function using
% different optimization algorithms. The available algorithms include
% fmincon, Simulated annealing, Pattern search, Genetic algorithm, Particle
% swarm, and Surrogate optimization. It iterates through the
% optimization_algorithms cell array and calls the appropriate optimization
% algorithm if its flag is set in the 'settings'.
%
% Inputs:
% - settings: A struct containing settings for the optimization algorithms.
% - model_folders: An array containing information about the models to be
% optimized.
%
% Outputs:
% - result: A vector containing the results of each optimization algorithm.
% - parameter_array: A matrix containing the parameters for each
% optimization algorithm result.
%
% Used Functions:
% - run_algorithm: Handles the execution of a specific optimization
% algorithm.
% - optimize_algorithm: Sets up the optimization algorithm and executes it.
% - run_optimization: Runs the selected optimization algorithm with
% appropriate settings.
% - run_optimizer: Executes the chosen optimization algorithm.
% - prepare_optimization_options: Prepares the optimization options based
% on the algorithm.
% - get_plot_functions: Returns the appropriate plotting function for the
% optimization algorithm.
% - f_opt_start: Returns the starting point and the starting population for
% optimization algorithms.
% 
% Loaded Variables:
% - optimization_algorithms: A cell array containing optimization algorithm
% names and flags.
% - algorithm_flag: The flag corresponding to the optimization algorithm.
% - objective_function: A function handle to the objective function for
% optimization.
% - options: A struct containing options for the optimization algorithm.
% - plot_functions: A cell array containing plot functions specific to each
% optimization algorithm.


% Set the random seed for reproducibility
rng(settings.rseed);

% Initialize the cell array containing optimization algorithm names and
% corresponding settings flags
optimization_algorithms = {
    'Genetic algorithm', 'ga';
    'fmincon', 'fmincon';
    'Simulated annealing', 'sa';
    'Pattern search', 'psearch';
   'Particle swarm', 'pswarm';
    'Surrogate optimization', 'sopt'
    };

p = gcp(); % If no pool, do not create new one.
poolsize = p.NumWorkers;
settings.pat = 1:poolsize;
for n = 1:poolsize
settings.pa(n,:) = settings.pa(1,:);
end
parfor n = 1:poolsize
[~] = f_diagnostics(settings,model_folders);
end

% Iterate through the optimization_algorithms cell array and call the
% corresponding algorithms if their flag is set in the settings
for i = 1:size(optimization_algorithms, 1)
    % algorithm_func = optimization_algorithms{i, 1};
    algorithm_flag = optimization_algorithms{i, 2};

    if settings.(algorithm_flag)
        [result(i), parameter_array(i, :)] = ...
            run_algorithm(@optimize_algorithm,...
            optimization_algorithms{i, 1}, settings, model_folders);
    end
end
end

function [algorithm_result, algorithm_parameter_array] =...
    run_algorithm(algorithm_func, optimization_algorithms_names,...
    settings, model_folders)

% Define the objective function
objective_function = @(x)f_sim_score(x, settings, model_folders);

% Call the optimization algorithm
algorithm_result = algorithm_func(optimization_algorithms_names,...
    settings, objective_function);

% Find the minimum value and its corresponding parameters
[~, min_index] = min(algorithm_result.fval);
algorithm_parameter_array = algorithm_result.x(min_index, :);
end

function result_opt = optimize_algorithm(algorithm_name,...
    settings, objective_function)
p = gcp(); % If no pool, do not create new one.
poolsize = p.NumWorkers;
mst = settings.mst;
if contains(algorithm_name,{'fmincon','Simulated annealing','Pattern search','Surrogate optimization'})
settings.mst = true;
settings.msts = poolsize;
else
settings.mst = mst;
settings.msts = poolsize;
end
% Determine the starting point for the optimization
[startpoint, ~, spop] = f_opt_start(settings);

% Prepare optimization options based on settings and algorithm_name
options = prepare_optimization_options(settings, algorithm_name);

% Execute the optimization algorithm
[x, fval, exitflag, output] = run_optimization(settings, algorithm_name,...
    objective_function, startpoint, spop, settings.lb,...
    settings.ub, options);

% Save optimization results
result_opt.name = algorithm_name;
result_opt.x = x;
result_opt.fval = fval;
result_opt.exitflag = exitflag;
result_opt.output = output;
end

function [x, fval, exitflag, output] = run_optimization(settings,...
    algorithm_name, objective_function, startpoint, spop, lb, ub,...
    options)

% Run the optimization algorithm using either multiple starting points (mst)
% or a single starting point (the average of the lower and upper bounds)
if settings.mst
    parfor n = 1:settings.msts
        disp(string(n))
        [x(n,:), fval(n), exitflag(n), output(n)] = ...
            run_optimizer(algorithm_name, objective_function,...
            startpoint(n,:), spop{n}, lb, ub, options);
    end
else
    [x(1,:), fval(1), exitflag(1), output(1)] = ...
        run_optimizer(algorithm_name, objective_function,...
        (lb + ub) / 2, spop, lb, ub, options);
end
end

function [x, fval, exitflag, output] = ...
    run_optimizer(algorithm_name, objective_function, startpoint, spop,...
    lb, ub, options)

% Execute the appropriate optimization algorithm based on algorithm_name
switch algorithm_name
    case 'fmincon'
        [x, fval, exitflag, output] = ...
            fmincon(objective_function, startpoint, [], [], [], [], lb, ub, [], options);
    case 'Simulated annealing'
        [x, fval, exitflag, output] = ...
        simulannealbnd(objective_function, startpoint, lb, ub, options);
    case 'Pattern search'
        [x, fval, exitflag, output] = ...
            patternsearch(objective_function, startpoint, ...
            [], [], [], [], lb, ub, [], options);
    case 'Genetic algorithm'
        parnum = length(lb);
        [x, fval] = ga(objective_function, parnum, ...
            [], [], [], [], lb, ub, [], options);
        exitflag = [" "];
        output = [" "];
    case 'Particle swarm'
        parnum = length(lb);
        [x, fval, exitflag, output] = ...
            particleswarm(objective_function, parnum, lb, ub, options);
    case 'Surrogate optimization'
        options.InitialPoints = spop;
        [x, fval, exitflag, output] = ...
        surrogateopt(objective_function, lb, ub, options);
    otherwise
        error('Unsupported optimization algorithm.');
end
end

function options = prepare_optimization_options(settings, algorithm_name)

% Set options for each algorithm
switch algorithm_name
    case 'fmincon'
        options = settings.fm_options;
        options.UseParallel = settings.optmc;
    case 'Simulated annealing'
        options =  settings.sa_options;
        options.MaxTime = settings.optt;
        options.InitialTemperature = ones(1,settings.parnum)*2;
    case 'Pattern search'
        options = settings.psearch_options;
        options.MaxTime = settings.optt;
        options.UseParallel = settings.optmc;
    case 'Genetic algorithm'
        options = settings.ga_options;
        options.MaxTime = settings.optt;
        options.PopulationSize = settings.popsize;
        options.UseParallel = settings.optmc;
        [~,startpop,~] = f_opt_start(settings);
        options.InitialPopulationMatrix = startpop;
    case 'Particle swarm'
        options = settings.pswarm_options;
        options.MaxTime = settings.optt;
        options.SwarmSize = settings.popsize;
        options.UseParallel = settings.optmc;
        [~,startpop,~] = f_opt_start(settings);
        options.InitialSwarmMatrix = startpop;
    case 'Surrogate optimization'
        options = settings.sopt_options;
        options.MaxTime = settings.optt;
        options.UseParallel = settings.optmc;
    otherwise
        error('Unsupported optimization algorithm.');
end

% Enable plots if chosen in settings
if settings.optplots
    % Set the appropriate PlotFcn based on the optimization algorithm
    options.PlotFcn = get_plot_functions(algorithm_name);
end

% Enable console messages if chosen in settings
if settings.optcsl
    options.Display = 'iter';
end
end

function plot_functions = get_plot_functions(algorithm_name)

% Set the appropriate PlotFcn based on the optimization algorithm
switch algorithm_name
    case 'fmincon'
        plot_functions = ...
            {@optimplotx, @optimplotfunccount, @optimplotfval, ...
            @optimplotstepsize, @optimplotfirstorderopt};
    case 'Simulated annealing'
        plot_functions = ...
            {@saplotbestf, @saplottemperature, @saplotf, ...
            @saplotstopping, @saplotx};
    case 'Pattern search'
        plot_functions = ...
            {@psplotbestf, @psplotfuncount, @psplotmeshsize, @psplotbestx};
    case 'Genetic algorithm'
        plot_functions = ...
            {@gaplotbestf, @gaplotexpectation, @gaplotrange, ...
            @gaplotscorediversity, @gaplotstopping, ...
            @gaplotscores, @gaplotdistance, @gaplotselection};
    case 'Particle swarm'
        plot_functions = {@pswplotbestf};
    case 'Surrogate optimization'
        plot_functions = {@surrogateoptplot};
    otherwise
        error('Unsupported optimization algorithm.');
end
end

function [spoint,spop_mc,spop_sc] = f_opt_start(stg)
% Set the randomm seed for reproducibility
rng(stg.rseed);

% Optimization Start method 1
if stg.osm == 1
    
    % Get a random starting point or group of starting points, if using
    % multistart, inside the bounds
    spoint = lhsdesign(stg.msts,stg.parnum).*(stg.ub-stg.lb)+stg.lb;

    % Get a group of ramdom starting points inside the bounds
    spop_mc = lhsdesign(stg.popsize,stg.parnum).*(stg.ub-stg.lb)+stg.lb;

    % Get a group of ramdom starting points inside the bounds
    for n = 1:stg.msts
    spop_sc{n} = lhsdesign(ceil(stg.popsize/36),stg.parnum).*(stg.ub-stg.lb)+stg.lb;
    end
    % Optimization Start method 2
elseif stg.osm == 2
    
    % Get a random starting point or group of starting points, if using
    % multistart, near the best point
    spoint = stg.bestpa - stg.dbs +...
        (stg.dbs*2*lhsdesign(stg.msts,stg.parnum));
    
    % Get a group of ramdom starting points near the best point
    spop_mc = stg.bestpa - stg.dbs +...
        (stg.dbs*2*lhsdesign(stg.popsize,stg.parnum));

    for n = 1:stg.msts
    spop_sc{n} = stg.bestpa - stg.dbs +...
        (stg.dbs*2*lhsdesign(ceil(stg.popsize/36),stg.parnum));
    end
end
end

Calls the correct optmizer or optimizers that have been chosen in the settings file.

• Inputs

  • stg - stg.fmincon, stg.sa, stg.psearch, stg.ga, stg.pswarm, stg.sopt

• Outputs - rst (optimization results)

f_opt_start

Code

Creates the starting parameter set or sets of the optimizations, if single or multistart selected in settings file.

It supports two different random distributions for the starting points.

• Inputs

  • stg - stg.rseed, stg.osm, stg.msts, stg.parnum, stg.ub, stg.lb, stg.popsize, stg.bestpa, stg.dbpa

• Outputs

  • spoint - (double) starting parameter set for the optimization

  • spop - (double) Starting parameter sets for multiple start optimizations

f_opt_fmincon/sa/psearch/ga/pswarm/sopt

Code

f_opt_fmincon

f_opt_sa

f_opt_psearch

f_opt_ga

f_opt_pswarm

f_opt_sopt

These functions call built in MATLAB® functions that perform parameter optimization . For furher information relating to how these optimizers work please follow the links to the MATLAB® documentation. Optimizers used:

• f_opt_fmincon - fmincon

• f_opt_sa - Simmulated annealing

• f_opt_psearch - Pattern search

• f_opt_ga - Genetic algorihtm

• f_opt_pswarm - Particle swarm

• f_opt_sopt - Surrogate optmization

• Inputs - stg

• Outputs - Optimization results

Global Sensitivity Analysis

f_gsa

Code

function results = f_gsa(settings,model_folders)
% f_gsa - Performs global sensitivity analysis (GSA) using the
% Sobol-Saltelli method on a given model
%
% This function conducts a GSA based on the Sobol-Saltelli method for a
% given model, as inspired by Geir Halnes et al.'s 2009 paper (J. comp.
% neuroscience 27.3 (2009): 471). The function takes two inputs:
% 'settings', a structure containing various settings for the analysis, and
% 'model_folders', the model folders structure. The function returns
% 'results', a structure containing the GSA results, including sensitivity
% indices (Si) and total sensitivity indices (SiT) for each score type.
%
% Inputs:
% - settings: A structure containing settings for the GSA, such as
% parameter bounds, sampling mode, and random seed.
% - model_folders: A structure containing the necessary model folder paths.
%
% Outputs:
% - results: A structure containing the GSA results, including sensitivity
% indices (Si) and total sensitivity indices (SiT) for each score type.
%
% Functions called:
% - f_make_par_samples(settings): Generates sample matrices of model
% parameters.
% - f_make_output_sample(results, settings, model_folders): Calculates the
% outputs for three different matrices (GSA M1, GSA M2, and GSA N) based on
% the input parameters, settings, and mathematical model.
% - f_calc_sensitivities(results, settings): Removes simulation errors
% from the data and calculates Si and SiT using the bootstrap method.

results = f_make_par_samples(settings);

results = f_make_output_sample(results,settings,model_folders);

results = f_calc_sensitivities(results,settings);
end

Calls the global sensitivity analysis functions in the correct order.

f_make_par_samples

Code

function results= f_make_par_samples(settings)
% This function creates sample matrices of model parameters according to
% specified distributions and settings, based on Geir Halnes et al.'s 2009
% paper (J. comp. neuroscience 27.3 (2009): 471). It generates two sample
% matrices M1 and M2, and a 3D matrix N, where N(:,:,i) is M2 with its i:th
% column replaced by M1(:,i).
%
% Inputs:
% - settings: A structure containing fields for generating samples,
% including sansamples, parnum, rseed, sasamplemode, lb, ub, bestpa, and
% sasamplesigma.
%
% Outputs:
% - results: A structure containing the following fields:
% - M1: Sample matrix 1.
% - M2: Sample matrix 2.
% - N: A 3D matrix where N(:,:,i) is M2 with its i:th column replaced by
% M1(:,i).
%
% Functions called:
% - makedist: Create a probability distribution object.
% - truncate: Truncate a probability distribution.
% - random: Generate random numbers from a probability distribution.
%
% Variables:
% Initialized:
% - M1: Sample matrix 1.
% - M2: Sample matrix 2.
% - N: A 3D matrix where N(:,:,i) is M2 with its i:th column replaced by
% M1(:,i).
% - pd: An array of probability distribution objects for each parameter.
%
% Loaded (from settings structure):
% - sansamples: Number of samples to generate for each parameter.
% - parnum: Number of model parameters.
% - rseed: Seed for the random number generator.
% - sasamplemode: Sampling mode.
% - lb: Vector containing lower bounds for each parameter.
% - ub: Vector containing upper bounds for each parameter.
% - bestpa: Vector containing the best parameter values.
% - sasamplesigma: Standard deviation for normal distribution-based
% sampling modes.

% Pre-allocate memory for sample matrices
M1 = zeros(settings.sansamples, settings.parnum);
M2 = zeros(settings.sansamples, settings.parnum);
N = zeros(settings.sansamples, settings.parnum, settings.parnum);
rng(settings.rseed)

% Loop through each parameter and create a distribution based on the
% specified settings(Note that the sampling is being performed in log space
% as the parameters and its bounds are in log space)
for i = 1:settings.parnum
    % Initialize distribution parameters depending on the sample mode
    switch settings.sasamplemode
        case 0
            % Uniform distribution with bounds defined by parameter limits
            lb = settings.lb(i);
            ub = settings.ub(i);
        case {1, 2}
            % Normal distribution centered at the best parameter value with
            % specified standard deviation
            mu = settings.bestpa(i);
            sigma = settings.sasamplesigma;
        case {3, 4}
            % Normal distribution centered at the mean of parameter bounds
            % with specified standard deviation
            mu = settings.lb(i) + (settings.ub(i) - settings.lb(i)) / 2;
            sigma = settings.sasamplesigma;
    end
    % Generate samples based on the distribution parameters
    if settings.sasamplemode == 0
        % Uniform distribution
        M1(:, i) = lb + (ub - lb) .* rand(1, settings.sansamples);
        M2(:, i) = lb + (ub - lb) .* rand(1, settings.sansamples);
    else
        % Normal distribution
        pd(i) = makedist('Normal', 'mu', mu, 'sigma', sigma);

        % Truncate the distribution if specified
        if ismember(settings.sasamplemode, [1, 3])
            pd(i) = truncate(pd(i), settings.lb(i), settings.ub(i));
        end
        M1(:, i) = random(pd(i), settings.sansamples, 1);
        M2(:, i) = random(pd(i), settings.sansamples, 1);
    end
end
% Replace the i:th column in M2 by the i:th column from M1 to create N
% matrices
for i=1:settings.parnum
    N(:,:,i) = M2;
    N(:,i,i) = M1(:,i);
end

% Store resulting matrices in the output structure
results.M1=M1;
results.M2=M2;
results.N=N;
end

Creates parameter sets samples with specific parameter distributions that are used to perform the global sensitivity analysis.

• Inputs

  • stg - stg.sansamples, stg.parnum, stg.sasamplemode, stg.ub, stg.lb

• Outputs - M1, M2, N

Code inspired by Geir Halnes et al. 2009 paper.

f_make_output_sample

Code

function results = f_make_output_sample(results,settings,model_folders)
% This function calculates the outputs for three different matrices (GSA
% M1, GSA M2, and GSA N) based on the input parameters, settings, and a
% mathematical model. This code is inspired by Geir Halnes et al. 2009
% paper.
%
% Inputs:
% - results: A structure containing M1, M2, and N matrices
% - settings: A structure with settings, including sansamples and parnum
% - model_folders: A variable containing folders for the models
%
% Outputs:
% - results: A structure containing fM1, fM2, and fN fields, each
% containing computed outputs for corresponding GSA methods
%
% Used Functions :
% - compute_f: Compute the results for GSA M1, GSA M2, and GSA N methods
% - f_sim_score: Simulation score calculation function
% - progress_track: Track progress of the calculations
%
% Variables
% Loaded:
% - nSamples: Number of samples (loaded from settings.sansamples)
% - nPars: Number of parameters (loaded from settings.parnum)
%
% Initialized:
% - time_begin: The start time of the function
% - f_out: Structure holding the computed values for fM1, fM2, and fN
%
% Persistent:
% - current_sample: Keeps track of the current sample number for progress
% tracking
% - last_time: Keeps track of the last time progress was printed

% Number of samples
nSamples = settings.sansamples;

% Number of parameters
nPars = settings.parnum;

% Start time of the function
time_begin = datetime;

% Compute the results for GSA M1, GSA M2, and GSA N methods
results.fM1 = compute_f(results.M1, nSamples, settings, model_folders, nPars, time_begin,"GSA M1","fM");
results.fM2 = compute_f(results.M2, nSamples, settings, model_folders, nPars, time_begin,"GSA M2","fM");
results.fN = compute_f(results.N, nSamples, settings, model_folders, nPars, time_begin,"GSA N","fN");
end

function [f_out] =...
    compute_f(parameter_array, nSamples, settings, model_folders, nPars, time_begin, task_name, matrix_type)

% Create a DataQueue for parallel processing
D = parallel.pool.DataQueue;

% Set up the progress tracker for the DataQueue
afterEach(D, @progress_track);

clear R RN

% Choose the matrix type and perform the calculations accordingly
switch matrix_type
    case 'fM'
        R = cell(nSamples,1);
      
        % Parallel loop for computing the output
        parfor i = 1:nSamples
            [~,~,R{i}] = f_sim_score(parameter_array(i,:), settings, model_folders);
            send(D, {task_name, 1, time_begin, nSamples, nPars});
        end
        % Extract the computed values

        sd = zeros(nSamples,size(R{1}.sd,1)*size(R{1}.sd,2));
        se = zeros(nSamples,size(R{1}.se,1));
        st = zeros(nSamples,size(R{1}.st,1));
        xfinal = zeros(nSamples,size(R{1}.sd,1)*(max(settings.exprun)-min(settings.exprun)+1));

        for i = 1:nSamples
            sd(i,:) = reshape(R{i}.sd(:,:), 1, []);
            se(i,:) = R{i}.se(:);
            st(i,:) = R{i}.st;
            xfinal(i,:) = [R{i}.xfinal{:}];
        end

    case 'fN'
        % Nested parallel loop for computing the output
        RN = cell(nSamples,nPars);
        
        parfor i = 1:nSamples
            for j = 1:nPars
                [~,~,RN{i,j}] = f_sim_score(parameter_array(i,:,j), settings, model_folders);
                
            end
            send(D, {task_name, 1, time_begin, nSamples, nPars});
          
        end

        sd = zeros(nSamples,size(RN{1,1}.sd,1)*size(RN{1,1}.sd,2),nPars);
        se = zeros(nSamples,size(RN{1,1}.se,1),nPars);
        st = zeros(nSamples,size(RN{1,1}.st,1),nPars);
        xfinal = zeros(nSamples,size(RN{1,1}.sd,1)*(max(settings.exprun)-min(settings.exprun)+1),nPars);

        for i = 1:nSamples
            for j = 1:nPars
                sd(i,:,j) = reshape(RN{i,j}.sd, 1, []);
                se(i,:,j) = RN{i,j}.se;
                st(i,:,j) = RN{i,j}.st;
                xfinal(i,:,j) = [RN{i,j}.xfinal{:}];
            end
        end
end

f_out.sd = sd;
f_out.se = se;
f_out.st = st;
f_out.xfinal = xfinal;

% Display the runtime information for the task
disp(task_name + " Runtime: " + string(datetime - time_begin) +...
    "  All " + nSamples + " samples executed")
end

% Function to track progress of the LSA
function progress_track(arg)
persistent current_sample
persistent last_time

% Initialize or update the current sample counter
if isempty(current_sample) || current_sample == arg{4}
    current_sample = arg{2};
end

% Retrieve other required arguments
task_name = arg{1};
start_time = arg{3};
num_samples = arg{4};
par_n = arg{5};
current_sample = current_sample + 1;

% Print progress information at each step
if mod(current_sample,ceil(num_samples/10)) == 0 &&...
        current_sample ~= num_samples
    
    if ((num_samples-current_sample)/num_samples)*10 <...
            ((num_samples-(num_samples/10))/num_samples)*10
        % Calculate the remaining time
        dt = (datetime-last_time);
       
        remaining_time = seconds(dt);
        remaining_time =...
            remaining_time*(num_samples-current_sample)/num_samples*10;
        remaining_time = seconds(remaining_time);
        remaining_time.Format = 'hh:mm:ss';
     
        % Print the progress and remaining time
        fprintf('%s Runtime: %s  Time to finish: %s  Samples: %d/%d\n', ...
            task_name, string(datetime - start_time), string(remaining_time), ...
            current_sample, num_samples);
    else

        % Print the progress without remaining time
        fprintf('%s Runtime: %s  Samples: %d/%d\n', ...
            task_name, string(datetime - start_time), ...
            current_sample, num_samples);
    end

    last_time = datetime;
end
end

For each parameter set given in the matrices M1, M2, and N it runs the function f_sim_score generating new matrices fM1, fM2, and fN respectively.

• Inputs - M1, M2, N, stg.sansamples, stg.parnum,

• Outputs - fM1, fM2, fN

Code inspired by Geir Halnes et al. 2009 paper.

f_calc_sensitivities

Code

function results = f_calc_sensitivities(results,settings)
% Computes the first-order (Si) and total-order (SiT) sensitivity indices
% for each score type present in the input data using the bootstrap method.
% This function takes in the raw data and various settings for the
% calculation process, and returns the results structure with calculated Si
% and SiT for each score type.
%
% Inputs:
% - results: a structure containing raw data (fM1, fM2, and fN) for
% different score types
% - settings: a structure containing settings for the calculation, such as
% bootstrap size, random seed, and error score
%
% Outputs:
% - results: a modified version of the input results structure containing
% the computed Si and SiT for each score type
%
% Functions called:
% - remove_sim_error: removes simulation errors from the input data
% - bootstrap: calculates Si and SiT using the bootstrap method
% - bootstrap_q: creates bootstrapped samples and calculates Si and SiT for
% the bootstrapped samples
% - calcSobolSaltelli: calculates Si and SiT using the Sobol-Saltelli
% method
% 
% Loaded variables:
% - None
%
% Initialized variables:
% - fM1: a structure containing the first-order score data for each score type
% - fM2: a structure containing the total-order score data for each score type
% - fN: a structure containing the score data for each score type and parameter
% - settings: a structure containing settings for the calculation
%
% Persistent variables:
% - None

% Remove simulation errors from the input data
results = remove_sim_error(results,settings);

% Calculate Si and SiT using the bootstrap method and store them in the results
% structure
[results.SiQ,results.SiTQ,results.Si,results.SiT] = bootstrap(results,settings);
end

function [SiQ,SiTQ,Si,SiT]=bootstrap(results,settings)
% calculates confidence intervals for the sensitivity indices using the
% bootstrap method.

% Initialize variables
fM1 = results.fM1;
fM2 = results.fM2;
fN = results.fN;

% Set the number of resamples for bootstrapping if not provided
if (isempty(settings.gsabootstrapsize))
    settings.gsabootstrapsize=ceil(sqrt(size(fM1.sd,2)));
end%if

% Define the list of scores to be used in the calculation
scores_list = ["sd","se","st","xfinal"];

% Generate random indices for bootstrapping
for j=1:settings.gsabootstrapsize
    rng(j*settings.rseed)
    I(j,:)=ceil(rand(size(fM1.sd,1),1)*size(fM1.sd,1));
end

% Calculate Si and SiT for each score type and perform bootstrapping
for n = 1:size(scores_list,2)
    fM1h = fM1.(scores_list(n));
    fM2h = fM2.(scores_list(n));
    fNh = fN.(scores_list(n));

    % Calculate Si and SiT without bootstrapping
    [Sih,SiTh] = calcSobolSaltelli(fM1h,fM2h,fNh,settings);
    Si.(scores_list(n)) = Sih;
    SiT.(scores_list(n)) = SiTh;

    % Initialize variables for bootstrapped Si and SiT
    SiQh = cell(1, settings.gsabootstrapsize);
    SiTQh = cell(1, settings.gsabootstrapsize);

    % Perform bootstrapping for the current score type
    parfor j=1:settings.gsabootstrapsize
        [SiQh{j},SiTQh{j}] = bootstrap_q(fM1h,fM2h,fNh,settings,j,I);
    end%parfor

    % Store bootstrapped Si and SiT values in the output structures
    for j=1:settings.gsabootstrapsize
        SiQ.(scores_list(n))(j,:,:) = SiQh{j};
        SiTQ.(scores_list(n))(j,:,:) = SiTQh{j};
    end
end
end%function

function [Si,SiT] = bootstrap_q(fM1,fM2,fN,settings,j,I)
    % Create bootstrapped samples using the generated indices
    fM1q = fM1(I(j,:),:);
    fM2q = fM2(I(j,:),:);
    fNq = fN(I(j,:),:,:);

    % Calculate Si and SiT for the bootstrapped samples
    [Si,SiT] = calcSobolSaltelli(fM1q,fM2q,fNq,settings);
end

function [Si,SiT] = calcSobolSaltelli(fM1,fM2,fN,settings)
% Calculates Si and SiT using the Sobol-Saltelli method.
% Code inspired by Geir Halnes et al. 2009 paper. (Halnes, Geir, et al. J.
% comp. neuroscience 27.3 (2009): 471.)

% Get the dimensions of the input data
[Nsamples,Nvars,Npars]=size(fN);

% Optionally subtract the mean to stabilize the model
if(settings.sasubmean)
    fM1 = fM1 - mean(fM1,1);
    fM2 = fM2 - mean(fM2,1);
    for i=1:Npars
        fN(:,:,i) =  fN(:,:,i) - mean(fN(:,:,i),1);
    end
end

% Calculate EY2 and variances
EY2 = mean(fM1.*fM2); % Valid definition (see Halnes et. al. Appendix)
VY = sum(fM1.^2)/(Nsamples-1) - EY2;
VYT = sum(fM2.^2)/(Nsamples-1) - EY2;

% Initialize Si and SiT arrays
Si = zeros(Nvars,Npars);
SiT= zeros(Nvars,Npars);

% Calculate Si and SiT for each parameter
for i=1:Npars
    Si(:,i) = (sum(fM1.*fN(:,:,i))/(Nsamples-1) - EY2)./VY;
    SiT(:,i) = 1 - (sum(fM2.*fN(:,:,i))/(Nsamples-1) - EY2)./VYT;
end
end

function results = remove_sim_error(results,settings)
% Removes any simulation errors present in the input data.

error_indices = any(results.fM1.sd == settings.errorscore, 2) | ...
                any(results.fM2.sd == settings.errorscore, 2) | ...
                any(any(results.fN.sd == settings.errorscore, 3), 2);

% Remove error rows from all fields of the input data
fields = fieldnames(results.fM1);
for i = 1:numel(fields)
    results.fM1.(fields{i})(error_indices, :) = [];
    results.fM2.(fields{i})(error_indices, :) = [];
    results.fN.(fields{i})(error_indices, :, :) = [];
end
end

Takes the matrices fM1, fM2, and fN and calculates sensitivity indexes. It calculates indexes based on the following oputputs of the f_sim_score function:

• The scores of each experimental output

• The scores of each experiment

• The total score

• The value of each experimental outputs at the end of the simulation

• Inputs - fM1, fM2, fN, stg.sasubmean

• Outputs - SI, STI

Code modified from the Geir Halnes et al. 2009 paper.

References

Halnes, G., Ulfhielm, E., Ljunggren, E.E., Kotaleski, J.H. and Rospars, J.P., 2009. Modelling and sensitivity analysis of the reactions involving receptor, G-protein and effector in vertebrate olfactory receptor neurons. Journal of Computational Neuroscience, 27(3), p.471.