shared_helper.m

function varargout = shared_helper(helper_name, varargin)
    %   Shared numerical diagnostics for the seven FBCIR scripts. Select 
    % 'regularise_posterior_covariance_history' or 
    % 'report_heatmap_display_diagnostic' through HELPER_NAME.

    switch helper_name
        case 'regularise_posterior_covariance_history'
            varargout{1} = regularise_posterior_covariance_history(varargin{:});
        case 'report_heatmap_display_diagnostic'
            report_heatmap_display_diagnostic(varargin{:});
        otherwise
            error('FBCIR:UnknownSharedHelper', 'Unknown shared helper ''%s''.', helper_name)
    end
end

function covariance_history = regularise_posterior_covariance_history(covariance_history, dt, covariance_label)
    % Regularise finite non-real, non-symmetric, or non-positive-semi-definite 
    % posterior filter/smoother covariances. Scalar (co)variances are simply 
    % made real and non-positive values are lifted to a machine-scale positive 
    % floor. Matrices are made real, symmetrised when required, and then 
    % projected to a positive definite matrix by positive-clipping only the 
    % eigenvalues below that same scale floor. Valid entries remain exactly 
    % unchanged, and no filtering or smoothing covariance recursion is altered.

    is_scalar_history = ndims(covariance_history) <= 2 && size(covariance_history, 1) == 1;
    if is_scalar_history
        history_length = numel(covariance_history);
    else
        history_length = size(covariance_history, 3);
    end

    corrected_count = 0;
    nonreal_count = 0;
    nonsymmetric_count = 0;
    nonspd_count = 0;
    first_corrected_index = NaN;
    last_corrected_index = NaN;
    worst_pre_repair_value = Inf;
    largest_spectral_floor = 0;

    for covariance_index = 1:history_length
        if is_scalar_history
            covariance_value = covariance_history(covariance_index);
        else
            covariance_value = covariance_history(:, :, covariance_index);
        end
        if any(~isfinite(covariance_value), 'all')
            error('FBCIR:NonFinitePosteriorCovariance', ...
                '%s posterior covariance is non-finite at time index %d (t = %0.6g); it cannot be regularised before ACI and CIR calculations.', ...
                covariance_label, covariance_index, (covariance_index-1)*dt)
        end

        if is_scalar_history
            is_nonreal = imag(covariance_value) ~= 0;
            repaired_value = real(covariance_value); % Make real if non-real via a simple real projection.
            is_nonsymmetric = false;
            is_nonspd = repaired_value <= 0;
            if ~(is_nonreal || is_nonspd)
                continue
            end
            spectral_floor = eps(max(1, abs(repaired_value))); % Machine-scale positive floor for non-positive scalar variance.
            worst_pre_repair_value = min(worst_pre_repair_value, repaired_value);
            if is_nonspd
                repaired_value = spectral_floor; % Lift non-positive scalar variance to the machine-scale positive floor.
            end
        else
            is_nonreal = any(imag(covariance_value) ~= 0, 'all');
            repaired_value = real(covariance_value); % Make real if non-real via a simple real projection.
            covariance_scale = max(1, norm(repaired_value, 'fro')); % Machine-scale positive floor for checking symmetry of matrix covariance.
            is_nonsymmetric = norm(repaired_value - repaired_value', 'fro') > eps(covariance_scale); 
            if is_nonreal || is_nonsymmetric
                repaired_value = 0.5*(repaired_value + repaired_value'); % Make symmetric if non-symmetric via a simple symmetrisation transformation.
            end
            % Check if the matrix is positive semi-definite via Cholesky 
            % decomposition. A non-zero flag indicates non-PSD, see flag output
            % of chol() in the MATLAB documentation:
            %   ➤ https://www.mathworks.com/help/matlab/ref/chol.html
            [~, positive_definite_flag] = chol(repaired_value);
            is_nonspd = positive_definite_flag ~= 0;
            if ~(is_nonreal || is_nonsymmetric || is_nonspd)
                continue
            end
            [eigenvectors, eigenvalues] = eig(repaired_value, 'vector'); % Compute eigenvalues and eigenvectors for the symmetric matrix.
            worst_pre_repair_value = min(worst_pre_repair_value, min(eigenvalues));
            spectral_floor = 0;
            if is_nonspd
                spectral_floor = eps(max(1, max(abs(eigenvalues)))); % Machine-scale positive floor for non-positive-semi-definite matrix covariance.
                eigenvalues(eigenvalues < spectral_floor) = spectral_floor; % Lift non-positive eigenvalues to the machine-scale positive floor.
                repaired_value = real(eigenvectors*diag(eigenvalues)*eigenvectors'); % Reconstruct the positive definite matrix from the lifted eigenvalues and eigenvectors using the spectral decomposition.
                repaired_value = 0.5*(repaired_value + repaired_value'); % Ensure the reconstructed matrix is symmetric via a simple symmetrisation transformation.
            end
        end

        corrected_count = corrected_count + 1;
        nonreal_count = nonreal_count + is_nonreal;
        nonsymmetric_count = nonsymmetric_count + is_nonsymmetric;
        nonspd_count = nonspd_count + is_nonspd;
        first_corrected_index = min(first_corrected_index, covariance_index);
        last_corrected_index = covariance_index;
        largest_spectral_floor = max(largest_spectral_floor, spectral_floor);
        if is_scalar_history
            covariance_history(covariance_index) = repaired_value;
        else
            covariance_history(:, :, covariance_index) = repaired_value;
        end
    end

    if corrected_count > 0
        fprintf(['Posterior covariance regularisation (%s): corrected %d/%d covariance(s) ' ...
                '(non-real=%d, non-symmetric=%d, non-SPD=%d); first and last t = %0.6g & %0.6g, ' ...
                'worst pre-repair eigenvalue/variance=%0.6g, largest spectral floor=%0.6g.\n'], ...
                covariance_label, corrected_count, history_length, nonreal_count, nonsymmetric_count, nonspd_count, ...
                (first_corrected_index-1)*dt, (last_corrected_index-1)*dt, ...
                worst_pre_repair_value, largest_spectral_floor);
    end
end

function report_heatmap_display_diagnostic(heatmap_metric, plot_time, direction)
    % Report display-only numerical adjustments before the optional logarithmic 
    % scale CIR metric heatmap plotting.
    nonfinite_mask = ~isfinite(real(heatmap_metric)) | ~isfinite(imag(heatmap_metric));
    nonreal_mask = ~nonfinite_mask & imag(heatmap_metric) ~= 0;
    nonpositive_mask = ~nonfinite_mask & ~nonreal_mask & real(heatmap_metric) <= 0;
    transformed_mask = nonreal_mask | nonfinite_mask | nonpositive_mask;

    raw_imag_values = abs(imag(heatmap_metric));
    raw_imag_values = raw_imag_values(isfinite(raw_imag_values));
    if isempty(raw_imag_values)
        max_imaginary_magnitude = NaN;
    else
        max_imaginary_magnitude = max(raw_imag_values);
    end
    if strcmp(direction, 'Forward')
        time_axis_label = 'natural time t';
    else
        time_axis_label = 'observational time T';
    end

    affected_time_indices = find(any(transformed_mask, 1));
    if isempty(affected_time_indices)
        first_affected_time = NaN;
        last_affected_time = NaN;
    else
        affected_physical_times = plot_time(affected_time_indices);
        first_affected_time = min(affected_physical_times);
        last_affected_time = max(affected_physical_times);
    end
    fprintf(['%s CIR metric heatmap display diagnostic: transformed %d/%d values ' ...
            '(non-real=%d, non-finite=%d, non-positive=%d); affected %s: first=%0.6g, last=%0.6g, ' ...
            'max |imag|=%0.6g.\n'], ...
            direction, nnz(transformed_mask), numel(heatmap_metric), nnz(nonreal_mask), ...
            nnz(nonfinite_mask), nnz(nonpositive_mask), time_axis_label, ...
            first_affected_time, last_affected_time, max_imaginary_magnitude);
end