Simulation and Scoring

f_sim_score

Code

function [score, rst, rst_not_simd] = ...
    f_sim_score(parameters, stg, model_folders, i, j)
% This function calculates the total score of a given model by simulating
% its performance and scoring the results.
%
% Inputs:
% - parameters: A double containing the model's parameters necessary for
% the simulation.
% - stg: A structure containing the settings and configurations for the
% simulation and scoring process.
% - model_folders: A structure containing the paths to the folders where
% the model and other relevant files are stored.
%
% Outputs:
% - score: A scalar value representing the total score of the model's
% performance.
% - rst: A structure containing the results of the simulation and scoring
% process.
% - rst_not_simd: A structure containing the results of the simulation and
% scoring process with the 'simd' field removed.
%
% Used Functions:
% - f_prep_sim: Simulates the model using the provided parameters,
% settings, and model folders.
% - f_score: Calculates the score of the model based on the results of the
% simulation.
%
% Variables:
% Loaded:
% - None.
%
% Initialized:
% - None.
%
% Persistent:
% - None.
%

%Turn off Dimension analysis warning from SimBiology
warning('off','SimBiology:DimAnalysisNotDone_MatlabFcn_Dimensionless')

% Call the function that simulates the model
rst = f_prep_sim(parameters, stg, model_folders, i, j);

% Call the function that scores
rst = f_score(rst, stg, model_folders);

% Get the total score explicitly for optimization functions
score = rst.st;

rst_not_simd = rmfield( rst , 'simd');
end

This function, f_sim_score, calculates the total score of a given model by simulating the model and scoring its performance. The function takes three input arguments: parameters, stg, and model_folders. The output of the function includes the total score (score), the result of the simulation and scoring (rst), and the result of the simulation and scoring without the ‘simd’ field (rst_not_simd).

• Inputs

  • stg - Structure containing the settings and configurations for the simulation and scoring process.

  • parameters - Double containing the model’s parameters that are necessary for the simulation.

  • model_folders: Structure containing the paths to the folders where the model and other relevant files are stored.

• Outputs

  • score - rst.st Scalar value representing the total score of the model’s performance.

  • rst - rst.simd, rst.xfinal, rst.sd, rst.se, and rst.st Structure containing the results of the simulation and scoring process.

  • rst_not_simd - rst.xfinal, rst.sd, rst.se, and rst.st Structure containing the results of the simulation and scoring process with the ‘simd’ field removed.

• Calls

  • f_prep_sim This function simulates the model using the provided parameters, settings, and model folders.

  • f_score This function calculates the score of the model based on the results of the simulation.

f_prep_sim

Code

function result = f_prep_sim(parameters, settings, model_folders, i, j)
% This function prepares the parameters for a simulation by setting them to
% the default values defined in the SBTAB and then updating any parameters
% that are being tested. The parameters are also adjusted according to any
% thermodynamic constraints.
%
% Inputs:
% - parameters: Array of parameter values to be tested
% - settings: Structure containing various settings for the simulation
% - model_folders: Structure with folder paths for model data
% - i, j: Indices for tracking simulations
%
% Outputs:
% - result: Structure with simulation results and updated parameter values
%
% Used Functions:
% - update_simulation_parameters: Function to update simulation parameters
% - f_sim: Function to run the simulation
%
% Variables:
% Loaded:
% - Data: Array of structures containing experimental data
% - sbtab: SBTAB structure containing default parameters and species
% information
%
% Initialized:
% - sim_par: Array of simulation parameters
% - ssa: Array of species start amounts
%
% Persistent:
% - sbtab: SBTAB structure containing default parameters and species
% information
% - Data: Array of structures containing experimental data

% Save variables across multiple function calls
persistent sbtab Data

% Load data on the first run
if isempty(sbtab) || isempty(Data)
    if isfield(model_folders, 'model') && isfield(model_folders.model, ...
            'data') && isfield(model_folders.model.data, 'data_model')
        data_model_path = model_folders.model.data.data_model;
        load(data_model_path, 'Data', 'sbtab')
    else
        error(['Data model not found in model_folders.' ...
            ' Please check the folder paths.'])
    end
end

% Set simulation parameters to default values from SBtab
if isfield(sbtab, 'defpar')
    sim_par(:, 1) = [sbtab.defpar{:, 2}]; %simulation_parameters
else
    error(['Default parameters not found in SBtab.' ...
        ' Ensure SBtab is properly formatted.'])
end

% Update parameters for Profile Likelihood if needed
if isfield(settings, 'PLind') && isfield(settings, 'PLval')
    if settings.PLind > 0 && settings.PLind <= length(parameters) + 1
        % add stg.PLval in the midle of other parameters
        parameters = ...
            [parameters(1:settings.PLind - 1), settings.PLval, ...
            parameters(settings.PLind:end)];
    else
        error(['PLind is out of bounds.' ...
            ' It should be within the range of parameter indices.'])
    end
end

% Update simulation parameters
sim_par = ...
    update_simulation_parameters(sim_par, parameters, settings, sbtab);

% Initialize the result structure
result = struct();
result.parameters = sim_par;

% Validate experiment run settings
if ~isempty(settings.exprun) && isnumeric(settings.exprun)
    result.simd = cell(1, max(settings.exprun));
else
    error(['stg.exprun is invalid.' ...
        ' Ensure exprun is a non - empty numeric array.'])
end

success_eq = zeros(1, settings.exprun(end));
error_list_E = "";

if settings.simcsl
    prep_sim_tic = tic;
end

% Run simulations for each experiment
for exp_indx = settings.exprun

    prev_success = success_eq(max(settings.exprun(1), exp_indx - 1));
    % Equilibration
    [species_start_amounts, success_eq(exp_indx), error_type] =...
        equilibrate_2(exp_indx, settings, sim_par, result, ...
        model_folders, sbtab, prev_success);

    % Main simulation
    if success_eq(exp_indx)

        tolerance_settings = ...
            create_tolerance_settings(exp_indx, settings, 0);

        [~, error_type, result] = ...
            run_with_tolerances(@(success) f_sim(exp_indx, settings, ...
            sim_par, species_start_amounts, result, ...
            model_folders, success), tolerance_settings, result);

        % if ~success_sim
        %     disp("n: " + n + " E" + (n - 1) + " fail_sim_2" + ...
        %         " last error: " + ME.identifier)
        %     result.simd{exp_indx} = 0;
        % end

        % Run detailed simulation if required
        if settings.simdetail
            [~, error_type, result] = ...
                run_with_tolerances(@(success) f_sim(exp_indx + 2 * ...
                settings.expn, settings, sim_par, species_start_amounts, ...
                result, model_folders, success), tolerance_settings, ...
                result);
        end

        % Check simulation output times
        if result.simd{exp_indx} ~= 0 && ...
                size(Data(exp_indx).Experiment.t, 1) ~= ...
                size(result.simd{exp_indx}.Data(:, end), 1)
            % disp("n: " + n + " E" + (n - 1) + " fail_sim_out_time" + ...
            %     " time_sim:  " + size(result.simd{n}.Data(:, end), 1) + ...
            %     " data_time: " + size(Data(n).Experiment.t, 1))
            error_type = "st"; %simulation time
            result.simd{exp_indx} = 0;
        end
    else
        result.simd{exp_indx} = 0;
    end
    if result.simd{exp_indx} == 0
        error_list_E = append(error_list_E, ...
            ((exp_indx - 1) + error_type + " "));
    end

end

% Display progress and errors if necessary
if settings.simcsl && ~(error_list_E == "")
    disp("i: " + i + " j: " + j + " E: " + error_list_E + ...
        " time: " + toc(prep_sim_tic))
end
end

function tolerance_settings = ...
    create_tolerance_settings(exp_indx, settings, equilibrate)
% Create a structure for tolerance settings
tolerance_settings = ...
    struct('exp_indx', exp_indx, 'reltol', settings.reltol, 'reltol_step', ...
    settings.reltol / 10, 'reltol_min', settings.reltol_min, ...
    'abstol', settings.abstol, 'abstol_step', settings.abstol / 10, ...
    'abstol_min', settings.abstol_min, 'equilibrate', equilibrate);
end

function [success, error_type, result] = ...
    run_with_tolerances(func, tolerance_settings, result)
% Run a function with varying tolerance settings until it succeeds
success = true;
error_type = '';

for reltol = tolerance_settings.reltol: - tolerance_settings.reltol_step:...
        min(tolerance_settings.reltol_min,tolerance_settings.reltol)

    % disp("reltol: " + reltol)
    for abstol = ...
            tolerance_settings.abstol: - tolerance_settings.abstol_step:...
            min(tolerance_settings.abstol_min,tolerance_settings.abstol)
        % disp("abstol: " + abstol)
        tolerance_settings.stg.reltol = reltol;
        tolerance_settings.stg.abstol = abstol;
        if tolerance_settings.equilibrate
            try
            [result, error_type] = func();
            success = true;
            catch
                error_type = "e"; %equilibration
                success = false;
            end
        else
            try
                result = func(success);
                success = true;
            catch
                % tolerance_settings.exp_indx
                result.simd{tolerance_settings.exp_indx} = 0;
                error_type = "s"; %simulation
                success = false;
            end
        end
        if success
            break;
        end
    end
    if success
        break;
    end
end
end

function [species_start_amounts_1, success_eq, error_type] = ...
    equilibrate_2(exp_indx, settings, sim_par, ...
    result, model_folders, sbtab, success_eq)
% Equilibrate the model for the specified experiment
persistent species_start_amounts %species start amount
error_type = "";

% Check if this is not the first experiment, the start values are equal to
% the previous experiment, and the previous simulation was valid
is_not_first_experiment = exp_indx ~= settings.exprun(1);

x = find(settings.exprun == exp_indx) - 1;

start_values_equal = min([sbtab.datasets(exp_indx).start_amount{:, 2}] == ...
    [sbtab.datasets(settings.exprun(max(x, 1))).start_amount{:, 2}]);
previous_simulation_valid = result.simd{settings.exprun(max(x, 1))} ~= 0;

if is_not_first_experiment && start_values_equal && ...
        previous_simulation_valid && success_eq
    % Set the start amounts based on the previous experiment
    species_start_amounts(:, exp_indx) =...
        species_start_amounts(:, settings.exprun(x));
    if settings.simdetail
        species_start_amounts(:, exp_indx + 2*settings.expn) = ...
            species_start_amounts(:, settings.exprun(x));
    end
    success_eq = true;
else
    % Set start amounts for species with non - zero values
    for j = 1:size(sbtab.datasets(exp_indx).start_amount, 1)
        % Set the start amount of the species to the number defined in
        % the sbtab for each experiment
        species_start_amounts(sbtab.datasets(exp_indx).start_amount{j, 3}, ...
            exp_indx + settings.expn) =...
            sbtab.datasets(exp_indx).start_amount{j, 2};
    end

    tolerance_settings = create_tolerance_settings(exp_indx, settings, 1);

    [success_eq, error_type, species_start_amounts] = ...
        run_with_tolerances(@() equilibrate(exp_indx, settings, sim_par, ...
        result, model_folders, sbtab, species_start_amounts, success_eq), ...
        tolerance_settings, species_start_amounts);
end
species_start_amounts_1 = species_start_amounts;
end

function sim_par = ...
    update_simulation_parameters(sim_par, parameters, settings, sbtab)
% Update simulation parameters based on input and thermodynamic constraints

% Iterate over all the parameters of the model
for n = 1:length(sim_par)
    % Update tested parameters
    if settings.partest(n) > 0
        sim_par(n) = 10 .^ parameters(settings.partest(n, 1));
    end
    % Update thermodynamic constrained parameters
    if isfield(settings, 'tci') && ~isempty(settings.tci) &&...
            ismember(n, settings.tci) && settings.partest(n) > 0
        sim_par = ...
            update_thermo_constrained(sim_par, parameters, settings, ...
            n, sbtab);
    end
end
end

function sim_par = ...
    update_thermo_constrained(sim_par, parameters, settings, n, sbtab)
% Update parameters constrained by thermodynamic laws

key_1 = ["tcm", "tcd"];
key_2 = ["*", "/"];

for k=1:2
    count = 0;
    for m = 1:size(settings.(key_1(k)), 2)
        if settings.(key_1(k))(n, m) ~= 0
            count = count + 1;
        end
    end
    % Iterate over the parameters that need to be multiplied for
    % calculating the parameter that depends on the thermodynamic
    % constraints
    for m = 1:count
        % Check that the parameter that is going to be used to calculate
        % the parameter dependent on thermodynamic constraints is not the
        % default
        if settings.partest(settings.(key_1(k))(n, m), 1) > 0
            % Check if the parameter needs to be set to the value relevant
            % for Profile Likelihood
            if isfield(settings, "PLind") && ...
                    settings.partest(settings.(key_1(k))(n, m), 1) == ...
                    settings.PLind
                parameters(settings.partest(settings.(key_1(k))(n, m), 1)) = ...
                    settings.PLval;
            end
            % Make the appropriate multiplications to get the
            % thermodynamically constrained parameter
            eval("sim_par(n) = sim_par(n) " + key_2(k) + ...
                " (10 ^ (parameters(settings.partest(" + ...
                "settings.(key_1(k))(n, m), 1))));");
        else
            % Make the appropriate multiplications to get the
            % thermodynamically constrained parameter
            eval("sim_par(n) = sim_par(n) " + key_2(k) + ...
                "(sbtab.defpar{settings.(key_1(k))(n, m), 2});");
        end
    end
end
end

function [species_start_amounts, error_type] = ...
    equilibrate(n, settings, sim_par, result, model_folders, sbtab, ...
    species_start_amounts, success)
% Equilibrate the model

result = f_sim(n + settings.expn, settings, sim_par, ...
    species_start_amounts, result, model_folders, success);
error_type = "";

if result.simd{n + settings.expn}.Time(end, 1) ~= settings.eqt
    % disp("n: " + n + " E" + (n - 1) + " time_eq: " + ...
    %     result.simd{n + settings.expn}.Time(end, 1) + ...
    %     " settings.eqt: " + settings.eqt)
    % disp("pecado")
    error("E" + (n - 1) + " fail_eq_out_time")
    error_type = "et"; %equilibration time
end

% Update the start amounts based on equilibrium results
for j = 1:size(sbtab.species, 1)
    final_amount = result.simd{n + settings.expn}.Data(end, j);
    if final_amount < 10 ^ -10
        species_start_amounts(j, n) = 0;
        if settings.simdetail
            species_start_amounts(j, n + 2 * settings.expn) = 0;
        end
    else
        species_start_amounts(j, n) = final_amount;
        if settings.simdetail
            species_start_amounts(j, n + 2 * settings.expn) = final_amount;
        end
    end
end
end

The code defines a main function f_prep_sim and three additional helper functions update_thermo_constrained_multiplications, update_thermo_constrained_divisions, and f_sim. The main function f_prep_sim prepares the parameters for a simulation, setting default values and updating any parameters being tested. It also adjusts the parameters according to thermodynamic constraints.

The main function f_prep_sim takes the following inputs:

• parameters: parameters for the simulation

• :ref:`stg<stg>`: settings for the simulation

• model_folders: folder paths for the models

And it outputs :

• rst: results of the simulation

The function initializes several persistent variables, imports data on the first run, and sets the default parameters for the simulation.

The function checks if the parameters need to be updated for Profile Likelihood.

It iterates through all model parameters, updating tested parameters and thermodynamic constrained parameters accordingly.

The function initializes the start amount for the species in the model to 0 and sets up a loop for each experiment being run.

Within the loop, the function tries to simulate the model, performing several checks and updates. If an error occurs during the simulation, the function catches the error and sets the simulation output to 0, indicating the simulation did not work properly.

The helper functions update_thermo_constrained_multiplications and update_thermo_constrained_divisions update the parameters according to the thermodynamic constraints. They iterate through parameters that need to be multiplied or divided, respectively, and make the appropriate adjustments.

The helper function f_sim runs simulations using SimBiology models for a set of experiments. It takes the following inputs:

• experiment_idx: indices of experiments to run

• settings: simulation settings

• simulation_parameters: parameter values for simulations

• species_start_amount: start amounts for species in simulations

• results: output variable to save simulation results

• main_model_folders: paths for model files

It outputs:

• results: simulation results

The function f_sim maintains the state of the loaded models between calls using persistent variables, loads the appropriate models, compiles the code for the simulation run, substitutes the start amounts of species and parameter values based on real-time results, and runs the simulation.

The simulation results are saved in the output variable, and the function can be called multiple times for different experiments. The function checks if the times of the simulation output and the simulation data from SBTAB match. If they do not match, it sets the simulation output to 0, indicating that the simulation did not work properly.

In summary, the main function f_prep_sim prepares the parameters for a simulation by setting them to the default values and then updating any parameters being tested. It adjusts the parameters according to any thermodynamic constraints and iterates through all the experiments to be run. The function then calls the helper function f_sim to run the simulation using SimBiology models for the set of experiments. The simulation results are saved in the output variable, and any errors encountered during the simulation are caught and handled appropriately.

f_sim

Code

function results = f_sim(experiment_idx,settings,simulation_parameters,...
    species_start_amount,results,main_model_folders,success)
% This function runs simulations using SimBiology models for a set of
% experiments. It loads the appropriate models and compiles the code for
% simulation run, then substitutes the start amounts of species and
% parameter values based on real-time results and runs the simulation. The
% results of the simulation are saved in the output variable. The function
% can be called multiple times for different experiments, and it maintains
% the state of the loaded models between calls using persistent variables.
%
% Inputs:
% - experiment_idx: Indices of experiments to run
% - settings: Simulation settings (e.g., sbioaccelerate, simdetail, expn,
% exprun)
% - simulation_parameters: Parameter values for simulations
% - species_start_amount: Start amounts for species in simulations
% - results: Output variable to save simulation results
% - main_model_folders: Paths for model files (e.g., model_exp_default,
% model_exp_eq, model_exp_detail)
%
% Outputs:
% - results: Simulation results (e.g., simd)
%
% Functions called:
% - sbioaccelerate: Compile model code for faster simulation run
% - sbiosimulate: Run simulation of SimBiology model with specified
% configuration
%
% Variables:
% Loaded:
% - model_exp: SimBiology model for each experiment
% - config_exp: SimBiology model configuration for each experiment
%
% Initialized:
% None
%
% Persistent:
% - models: Cell array containing the SimBiology models for each experiment
% - configs: Cell array containing the configurations for each SimBiology
% model

% Save variables that need to be mantained over multiple function calls
persistent models
persistent configs

rel_tol_ind = abs(log10(settings.reltol_min)-log10(settings.reltol))+1;
abs_tol_ind = abs(log10(settings.abstol_min)-log10(settings.abstol))+1;
% If the function is called for the first time, load the appropriate model
% and compile the code for simulation run
if isempty(models)

    % Turn off warning messages
    warning('off','all')

    % Generate an empty array to be populated with the model suited for
    % each equilibration and experiment
    models = cell(1, 1, settings.expn * (2 + settings.simdetail));
    configs = cell(1, 1, settings.expn * (2 + settings.simdetail));

    % Set the file paths for the different models
    model_exp_default = main_model_folders.model.data.model_exp.default;
    model_exp_eq = main_model_folders.model.data.model_exp.equilibration;
    model_exp_detail = main_model_folders.model.data.model_exp.detail;

    % Iterate over the experiments thar are being run
    for n = settings.exprun
        for j = 1:rel_tol_ind
            for k = 1:abs_tol_ind
                % Load and compile SimBiology models for each experiment

                % Load models for main simulation
                load(model_exp_default + n + ".mat", 'model_exp', ...
                    'config_exp')
                models{j, k, n} = model_exp;
                configs{j, k, n} = config_exp;

                % Load models for equilibrium
                load(model_exp_eq + n + ".mat", 'model_exp', 'config_exp')
                models{j, k, n+settings.expn} = model_exp;
                configs{j, k, n+settings.expn} = config_exp;

                % Load models for detailed simulation
                if settings.simdetail
                    load(model_exp_detail + n + ".mat", 'model_exp', ...
                        'config_exp')
                    models{j, k, n+2*settings.expn} = model_exp;
                    configs{j, k, n+2*settings.expn} = config_exp;
                end
                % Compile the model code if the option is selected in
                % settings
                if settings.sbioacc
                    sbioaccelerate(models{j, k, n}, configs{j,k,n});
                    sbioaccelerate(models{j, k, n+settings.expn}, ...
                        configs{j, k, n+settings.expn});
                    if settings.simdetail
                        sbioaccelerate(models{j, k, n+2*settings.expn}, ...
                            configs{j, k, n+2*settings.expn});
                    end
                end
            end
        end
    end
end

% substitute the start amount of the species in the model with the correct
% ones for  simulations
set(models{rel_tol_ind,abs_tol_ind,experiment_idx}.species(...
    1:size(species_start_amount(:,experiment_idx),1)),...
    {'InitialAmount'}, num2cell(species_start_amount(:,experiment_idx)));

% Substitute the values of the parameters in the model for the correct one
% for simultaions
set(models{rel_tol_ind,abs_tol_ind,experiment_idx}.parameters(...
    1:size(simulation_parameters,1)),...
    {'Value'}, num2cell(simulation_parameters));

if ~success
    configs_fail = configs{rel_tol_ind,abs_tol_ind,experiment_idx};
    configs_fail.SolverOptions.AbsoluteTolerance = settings.abstol;
    configs_fail.SolverOptions.RelativeTolerance = settings.reltol;
    results.simd{experiment_idx} = ...
    sbiosimulate(models{rel_tol_ind,abs_tol_ind,experiment_idx},...
        configs_fail);
else
    % simulate the model using matlab built in function
    results.simd{experiment_idx} = ...
        sbiosimulate(models{rel_tol_ind,abs_tol_ind,experiment_idx},...
        configs{rel_tol_ind,abs_tol_ind,experiment_idx});
end
end

Description

Simulates the model with the provided configurations. The first time it is run it loads a representation of the model and the simulation, and compiles this information to C code.

Input Arguments

• exp_n - (double) Unique number to identify the model for each experiment or equilibrium reaction (it needs a new model object for each one)

• stg - stg.expn, stg.name, stg.sbioacc

• rt

  • rt.ssa - (double) steady state amounts

  • rt.par - (double) All parameters of the model, takes the default ones from SBtab and then replaces the ones being worked on.

• rst - rst.simd

Output Arguments

• rst - rst.simd

Functions called - Sbioaccelerate, Sbiosimulate Loaded variables - Ready to run model, Ready to run model equilibration

f_score

Code

function results = f_score(results, settings, model_folders)
% This function calculates the score for a set of simulated results by
% comparing them to experimental data. It computes the score for each
% dataset and experiment and then calculates the total score based on the
% selected scoring strategy. The score serves as a metric for comparing the
% accuracy of different simulations or models.
%
% Inputs:
% - rst: Structure containing the simulation results and scores.
% - stg: Structure containing the settings for the scoring strategy, such
% as the option to use log10 scaling, error score, and other options.
% - mmf: Structure containing the model directory information.
%
% Outputs:
% - rst: Updated structure containing the calculated scores for each
% dataset, experiment, and the total score.
%
% Used Functions:
% - calculate_score_per_experiment: Calculates the score for each
% experiment.
% - calculate_score_per_dataset: Calculates the score for each dataset in
% an experiment.
% - f_normalize: Normalizes the simulation results for a dataset.
%
% Variables:
% Loaded:
% - Data: Experimental data.
% - sbtab: sbtab structure containing dataset information.
%
% Persistent:
% - sbtab: sbtab structure containing dataset information (persistent).
% - Data: Experimental data (persistent).

persistent sbtab
persistent Data

data_model = model_folders.model.data.data_model;

% Import the data on the first run
if isempty(Data)
    load(data_model,'Data','sbtab')
end

% Iterate over the number of experiments
for n = settings.exprun
    results = calculate_score_per_experiment(results, settings, n, ...
        model_folders, sbtab, Data);
end

% Calculate the total score
results.st = sum(results.se);

% Calculate the log10 of total score if option selected
if settings.useLog == 3
    results.st = log10(results.st);
end
end

function results = ...
    calculate_score_per_experiment(results, settings, n, model_folders, ...
    sbtab, Data)

% Iterate over the number of datasets per experiment
for j = 1:size(sbtab.datasets(n).output, 2)
    results = ...
        calculate_score_per_dataset(results, settings, n, j, ...
        model_folders, Data);
end
% Calculate score per experiment
results.se(n, 1) = sum(results.sd(:, n));

% Calculate the log10 of experiment score if option selected
if settings.useLog == 2
    results.se(n, 1) = log10(results.se(n, 1));
end

end

function results = ...
    calculate_score_per_dataset(results, settings, n, j, model_folders, ...
    Data)

% Calculate score per dataset if there are no errors
if results.simd{n} ~= 0
    data = Data(n).Experiment.x(:, j);
    data_sd = Data(n).Experiment.x_SD(:, j);

    number_points = size(Data(n).Experiment.x(:, j), 1);
    [sim_results_norm,~,sim_results] = ...
        f_normalize(results, settings, n, j, model_folders);

    if ~isempty(sim_results_norm)
        sim_results = sim_results_norm;
    end
    results.xfinal{n, 1}(j) = sim_results(end);

    % Calculate score using formula that accounts for normalization
    % with the starting point of the result
    switch settings.useLog
        case {0, 2, 3}
            results.sd(j, n) = ...
                sum(((data - sim_results) ./ (data_sd)).^2) / ...
                number_points;
        case 1
            results.sd(j, n) = ...
                max(0, log10(sum(((data - sim_results) ./ ...
                (data_sd)).^2) / number_points));
        case 4
            results.sd(j, n) = ...
                sum(((data - sim_results) ./ ...
                (data_sd * sqrt(number_points))).^2);
        otherwise
            error('Invalid value for stage.useLog: %d', settings.useLog);
    end

    if results.sd(j, n) == inf || isnan(results.sd(j, n))
        results.sd(j, n) = settings.errorscore;
        results.xfinal{n, 1}(j) = 0;
    end

    % If there are errors output a very high score value (10^10)
elseif results.simd{n} == 0 %|| rst.sd(j, n) == inf || isnan(rst.sd(j, n))
    switch settings.useLog
        case {1}
            results.sd(j, n) = log10(settings.errorscore);
        case {0, 2, 3, 4}
            results.sd(j, n) = settings.errorscore;
    end
    results.xfinal{n, 1}(j) = 0;
end
end

Description

The f_score function computes the score for a given set of simulated results by comparing them with the experimental data. The function calculates the score for each dataset and experiment, and then computes the total score based on the selected scoring strategy. The score serves as a metric for comparing the accuracy of different simulations or models.

Input Arguments

• rst: Structure containing the simulation results and scores.

• stg: Structure containing the settings for the scoring strategy, such as the option to use log10 scaling, error score, and other options.

• mmf: Structure containing the model information, including the data model.

Output Arguments

• rst.st: Updated structure containing the calculated scores for each dataset, experiment, and the total score.

Example

% Define the input structures and settings
stg.useLog = 1;
stg.errorscore = 1e10;
stg.exprun = 1:3;
mmf.model.data.data_model = matlab model file
rst.simd = matlab output from f_prep_sim

% Call the f_score function
rst = f_score(rst, stg, mmf);

This example demonstrates how to call the f_score function with the input structures and settings. The function calculates the scores for the given simulation results based on the scoring strategy defined in the stg structure.

f_normalize

Code

function [sim_results_norm,sim_results_detailed,sim_results] = ...
f_normalize(results, settings, exp_number, output_number, model_folders)
% This function processes and normalizes simulation results based on a
% specified normalization method. It accepts a set of inputs, including the
% simulation results, settings, experiment and output numbers, and a model
% metafile structure. The function returns normalized simulation results,
% along with detailed normalized simulation results if the 'simdetail'
% setting is enabled.
%
% Inputs:
% - rst: Structure containing the simulation results and scores
% - stg: Structure containing the settings for the simulation
% - exp_number: An integer representing the experiment number
% - output_number: An integer representing the output number
% - mmf: Structure containing the model directory information
%
% Outputs:
% - sim_results: A matrix containing the normalized simulation results
% - sim_results_detailed: A matrix containing the detailed normalized
% simulation results (if stg.simdetail is true)
%
% Used Functions:
% - extract_data: Helper function to extract data from the simulation
% results
%
% Variables
% Persistent:
% - sbtab: SBtab structure
% - Data: Loaded data from the model
% - sb: Structure containing experiment IDs

persistent sbtab
persistent Data
persistent sb

sim_results_norm = [];

% Load Data and sbtab if empty
if isempty(Data)
    load(model_folders.model.data.data_model, 'Data', 'sbtab','sb')
end

% Extract the data from the simulation results
sim_results = extract_data(results, exp_number, output_number, sbtab);
sim_results_detailed = [];

% Check if detailed simulation results are requested
if settings.simdetail
    sim_results_detailed =...
        extract_data(results, exp_number + 2 * settings.expn, ...
        output_number, sbtab);
end

% Get the normalization method
normalize = sbtab.datasets(exp_number).Normalize;

% Perform normalization if a method is specified
if ~isempty(normalize)
    output_ID = sbtab.datasets(exp_number).output_ID{output_number}{:};
    norm_factor = extract_data(results, exp_number, output_number, sbtab);
    sim_results_norm = sim_results;

    % Normalize by the maximum value
    if contains(normalize, 'Max') && contains(normalize, output_ID)
        max_norm = max(norm_factor);
        if max_norm == 0
            max_norm = 1;
        end
        sim_results_norm = sim_results / max_norm;
        if settings.simdetail
            sim_results_detailed = sim_results_detailed / max_norm;
        end
    end

    % Normalize to a value between 0 and 1
    if contains(normalize, 'Norm') && contains(normalize, output_ID)
        max_norm = max(norm_factor(2:end));% Here because sometimes the first value presents strange results
        min_norm = min(norm_factor(2:end));% Here because sometimes the first value presents strange results
        if max_norm-min_norm == 0
            max_norm = 1;
            min_norm = 0;
        end
        sim_results_norm = (sim_results-min_norm) / (max_norm-min_norm);
        sim_results_norm(1) = sim_results_norm(2);% Here because sometimes the first value presents strange results
        if settings.simdetail
            sim_results_detailed = (sim_results_detailed-min_norm) / (max_norm-min_norm);
        end
    end


    % Normalize by the minimum value
    % if contains(normalize, 'Min') && contains(normalize, output_ID)
    % 
    %     min_norm = min(norm_factor);
    %     sim_results = sim_results / min_norm;
    %     min_norm
    %     if settings.simdetail
    %         sim_results_detailed = sim_results_detailed / min_norm;
    %     end
    % end

    % Normalize by time-dependent factors
    if contains(normalize, 'Time') && contains(normalize, output_ID)
        t_size = size(Data(exp_number).Experiment.t, 1);
        for n = 1:t_size
            exp_ID = getfield(sb, ['E' num2str(exp_number - 1) '.ID'], n);
            if contains(normalize, exp_ID)
                sim_results_norm = sim_results / norm_factor(n);
                if settings.simdetail
                    sim_results_detailed = ...
                    sim_results_detailed / norm_factor(n);
                end
            end
        end
    end
end
end

% Helper function to extract data
function data = extract_data(results, exp_number, output_number, sbtab)
% Helper function to extract data from the simulation results
% Input arguments:
% rst  - Structure containing the simulation results and scores.
% stg  - Structure containing the settings for the simulation.
% output_number - An integer representing the output number.
% sbtab - SBtab structure
% Output argument:
% data - extracted data from the simulation results
    data = results.simd{1, exp_number}.Data(:, end - ...
        size(sbtab.datasets(exp_number).output, 2) + output_number);
end

Description

The f_normalize function processes and normalizes simulation results based on a specified normalization method. It accepts a set of inputs, including the simulation results, settings, experiment and output numbers, and a model metafile structure. The function returns normalized simulation results, along with detailed normalized simulation results if the ‘simdetail’ setting is enabled.

Input Arguments

• rst: A structure containing the simulation results.

• stg: A structure containing the settings for the simulation.

• exp_number: An integer representing the experiment number.

• output_number: An integer representing the output number.

• mmf: A structure containing the model metafile information.

Output Arguments

• sim_results: A matrix containing the normalized simulation results.

• sim_results_detailed: A matrix containing the detailed normalized simulation results (if stg.simdetail is true).