main.m

clear; clc; tic;

% Load experimental data
expDataPrep  % Produces: Rmatrix, Rdotmatrix, noisyR, noisyRdot, R_stats, Rdot_stats, tmatrix

% Parameters
variability = 1:10;                     % noise realizations
betaVals = 0.1:0.1:1.0;                  % std. deviation scaling factors
material = 'UM1_Prob.mat';               % output filename
addpath('UM1')

% Parallel setup
if isempty(gcp('nocreate'))
    parpool('local');
end

% Model setup
modelFiles = {'newtonian.mat', 'NH.mat', 'NHKV.mat', 'qNH.mat', 'linmax.mat', 'qKV.mat', 'SLS.mat'};
modelNames = {'newtonianProb', 'NHProb', 'NHKVProb', 'qNHProb', 'linmaxProb', 'qKVProb', 'SLSProb'};
nModels = numel(modelFiles);

% Preallocate futures
futures = parallel.FevalFuture.empty(nModels, 0);

% Launch all returnProb jobs asynchronously
for i = 1:nModels
    futures(i) = parfeval(@returnProb, 1, ...
        modelFiles{i}, Rmatrix, Rdotmatrix, ...
        R_stats, Rdot_stats, tmatrix, betaVals);
end

% Track progress and collect results
results = cell(1, nModels);
for i = 1:nModels
    [completedIdx, value] = fetchNext(futures);
    results{completedIdx} = value;
end

% Assign results to named variables in base workspace
for i = 1:nModels
    assignin('base', modelNames{i}, results{i});
end

% Save results
save(material, modelNames{:});
toc

% ================== return probability ==================
function modelProb = returnProb(model, Rexps, Rdotexps, R_stats, Rdot_stats, tmatrix, betaVals)
    % modelProb is a 1×nBeta struct with fields:
    %   .beta  - the sigma scaling used
    %   .prob  - cell ND array matching Sims grid with scalar likelihoods

    % --- Load sims (ND cell grid) ---
    m = load(model);
    Sims = m.Sims;
    dims = size(Sims);
    nd   = numel(dims);

    % --- All ND index combinations ---
    indices = cell(1, nd);
    for d = 1:nd, indices{d} = 1:dims(d); end
    [grids{1:nd}] = ndgrid(indices{:}); 
    allSubs = cell2mat(cellfun(@(x) x(:), grids, 'UniformOutput', false));

    % --- Baseline sigmas from stats (2nd column = stdev) ---
    if size(R_stats,2) < 2 || size(Rdot_stats,2) < 2
        error('R_stats and Rdot_stats must have at least 2 columns (mean, std).');
    end
    sigmaR_base    = R_stats(:,2);
    sigmaRdot_base = Rdot_stats(:,2);

    % --- Output struct over beta ---
    nBeta = numel(betaVals);
    modelProb = struct('beta', cell(1, nBeta), 'prob', []);

    % --- Loop over beta values ---
    for b = 1:nBeta
        beta = betaVals(b);

        % Scale per-timestep sigmas
        sigmaR    = max(beta .* sigmaR_base,    1e-12);
        sigmaRdot = max(beta .* sigmaRdot_base, 1e-12);

        % Likelihoods on the Sims grid
        probArray = cell(dims);
        for idx = 1:size(allSubs, 1)
            subs = num2cell(allSubs(idx, :));
            currentSims = Sims{subs{:}};   % [Nt x 3] -> [t, R, Rdot]

            probArray{subs{:}} = modelProbComp(currentSims, Rexps, Rdotexps, ...
                                               sigmaR, sigmaRdot, tmatrix, 3.0);
        end

        modelProb(b).beta = beta;
        modelProb(b).prob = probArray;
    end
end