Main_Calcium_Code.m

%
delete(gcp('nocreate'))
close all
clearvars
clc
%% Add folder containing this code and all functions in the folder to the path

% --- Robust Path Detection (Works with "Run" and "Run Section") ---

% Get the full path of the currently active script in the MATLAB editor
try
    active_editor_doc = matlab.desktop.editor.getActive;
    full_script_path = active_editor_doc.Filename;
catch
    error('Could not get active editor document. Please ensure the script is saved and active in the editor.');
end

% Check if the editor returned a valid path
if isempty(full_script_path)
    error('The active script has no filename. Please save the file first.');
end

% Get the directory containing the script
scriptPath = fileparts(full_script_path);

% Add the script's folder and all its subfolders to the MATLAB path
addpath(genpath(scriptPath)); 

fprintf('Successfully added path and subfolders for: %s\n', scriptPath);

%% Pre-Processing Settings App
% Uncomment line 15 when wanting to change the signal processing parameters
% used throughout the code
%
% PreProcessingSettingsApp

%% DETERMINE WHAT KIND OF DATA MATLAB IS ANALYZING

prompt = 'How many time-points does this samples dataset contain? \n 1 - 1 timepoint \n 2 - 2 timepoints \n 3 - 3 timepoints \n 4 - 4 timepoints \n etc. \n';
num_timepoints_expected = input(prompt); % Renamed for clarity

if isempty(num_timepoints_expected) || ~isnumeric(num_timepoints_expected) || num_timepoints_expected < 1
    error('Invalid number of timepoints entered.');
end

%% INTEGRATE SETTINGS SAVED FROM THE PRE-PROCESSING SETTINGS APP

settingsFilename = 'preprocessing_settings.mat';
settingsFile = fullfile(scriptPath, settingsFilename); % Construct the full path to the settings file

% Load smoothing factor from settings or fallback
if isfile(settingsFile) % Use isfile for clarity
    settings = load(settingsFile);
    fprintf('Loaded settings from %s\n', settingsFile);
else
    settings = struct();
    warning('Preprocessing settings file (%s) not found. Using default values.', settingsFile);
end

% Define defaults
defaultSettings.n = 60;
defaultSettings.ProminenceValue = 1.5;
defaultSettings.use_fixed_threshold = false; % Default to auto-calculating threshold in main script if not set
defaultSettings.signalStrengthThresholdValue = 5; % Default fixed threshold value if fixed mode is on but no value set
defaultSettings.useBounds = false;
defaultSettings.minBound = -50;
defaultSettings.maxBound = 50;
defaultSettings.analysisDurationSeconds = 120; % Default analysis duration
defaultSettings.timepointOrderInApp = {}; % Default to empty if not in file

% Helper function to get setting or default
getSetting = @(fieldName, defaultValue) getfield(settings, fieldName, 'ifnofieldfound', defaultValue);
% This custom helper handles missing fields more gracefully
function value = getfield_with_default(s, field, default)
    if isfield(s, field)
        value = s.(field);
    else
        value = default;
    end
end

% Apply settings, using defaults if fields are missing
settings.n = getfield_with_default(settings, 'n', defaultSettings.n);
settings.ProminenceValue = getfield_with_default(settings, 'ProminenceValue', defaultSettings.ProminenceValue);
settings.use_fixed_threshold = getfield_with_default(settings, 'use_fixed_threshold', defaultSettings.use_fixed_threshold);
settings.signalStrengthThresholdValue = getfield_with_default(settings, 'signalStrengthThresholdValue', defaultSettings.signalStrengthThresholdValue);
settings.useBounds = getfield_with_default(settings, 'useBounds', defaultSettings.useBounds);
settings.minBound = getfield_with_default(settings, 'minBound', defaultSettings.minBound);
settings.maxBound = getfield_with_default(settings, 'maxBound', defaultSettings.maxBound);
settings.analysisDurationSeconds = getfield_with_default(settings, 'analysisDurationSeconds', defaultSettings.analysisDurationSeconds);
settings.timepointOrderInApp = getfield_with_default(settings, 'timepointOrderInApp', defaultSettings.timepointOrderInApp);

% For convenience in processTimepoint, ensure signal_strength_threshold exists even if not fixed
% This will be overridden if use_fixed_threshold is false later
if ~isfield(settings, 'signal_strength_threshold') % Legacy or direct use
    settings.signal_strength_threshold = settings.signalStrengthThresholdValue;
end
%% LOAD FDATA GENERATED BY NEUROCA FOR ALL TIMEPOINTS OF ONE SAMPLE

fprintf('---------------------------------------------------------\n');
fprintf('LOADING DATA FOR CURRENT SAMPLE\n');
fprintf('---------------------------------------------------------\n');

% Prompt user to select all timepoint folders
rootFolder = uigetdir(pwd, 'Select the ROOT FOLDER for the CURRENT SAMPLE containing timepoint subfolders');
if rootFolder == 0
    error('Sample folder selection cancelled. Exiting.');
end
fprintf('Selected sample root folder: %s\n', rootFolder);

% Initialize containers
fdata_original_all_loaded = {};
timepointTags_loaded = {}; % Store the tags of successfully loaded timepoints

% Check timepoint order (imperative to identify baseline/timepoint 1)
if ~isempty(settings.timepointOrderInApp)
    % --- Preferred method: Use order from settings ---
    fprintf('Attempting to load timepoints based on order from settings: %s\n', strjoin(settings.timepointOrderInApp, ', '));
    
    if num_timepoints_expected > numel(settings.timepointOrderInApp)
        warning(['Number of timepoints to process for this sample (%d) is greater than ',...
                 'the number of timepoint names defined in settings (%d). Will only load up to %d.'], ...
                 num_timepoints_expected, numel(settings.timepointOrderInApp), numel(settings.timepointOrderInApp));
    end
    
    num_to_try_load = min(num_timepoints_expected, numel(settings.timepointOrderInApp));

    for i = 1:num_to_try_load
        tag = settings.timepointOrderInApp{i};
        expected_subfolder_path = fullfile(rootFolder, tag);
        fdataPath = fullfile(expected_subfolder_path, 'fdata.mat');
        
        fprintf('  Looking for Timepoint "%s" in %s\n', tag, expected_subfolder_path);
        if ~isfolder(expected_subfolder_path)
            warning('  Subfolder "%s" not found. Skipping.', tag);
            continue;
        end
        if isfile(fdataPath)
            try
                fdata_struct = load(fdataPath); % Returns a struct
                fnames = fieldnames(fdata_struct);
                if isempty(fnames) error('fdata.mat for %s is empty', tag); end
                fdata_original_all_loaded{end+1} = fdata_struct.(fnames{1});
                timepointTags_loaded{end+1} = tag;
                fprintf('    Successfully loaded "%s".\n', tag);
            catch ME
                warning('    Error loading fdata.mat from %s: %s. Skipping.', tag, ME.message);
            end
        else
            warning('    No fdata.mat found in %s. Skipping.', expected_subfolder_path);
        end
    end
else
    % --- Fallback: Your original method using dir() if settings.timepointOrderInApp is empty ---
    warning('No timepoint order in settings. Using subfolders found in root (order might be OS dependent).');
    subfolders = dir(rootFolder);
    subfolders = subfolders([subfolders.isdir] & ~startsWith({subfolders.name}, '.'));

    if numel(subfolders) < num_timepoints_expected
        warning(['Found %d subfolders, but expected %d timepoints. ', ...
                 'Will attempt to load from available subfolders. Ensure correct naming and presence.'], ...
                 numel(subfolders), num_timepoints_expected);
    end
    
    num_to_try_load = min(num_timepoints_expected, numel(subfolders));

    for i = 1:num_to_try_load % Loop through found subfolders up to num_timepoints_expected
        tag = subfolders(i).name;
        fdataPath = fullfile(rootFolder, tag, 'fdata.mat');
        fprintf('  Looking for Timepoint "%s" in %s (fallback method)\n', tag, fullfile(rootFolder,tag));
        if isfile(fdataPath)
            try
                fdata_struct = load(fdataPath); % Changed from importdata for consistency
                fnames = fieldnames(fdata_struct);
                if isempty(fnames) error('fdata.mat for %s is empty', tag); end
                fdata_original_all_loaded{end+1} = fdata_struct.(fnames{1});
                timepointTags_loaded{end+1} = tag;
                fprintf('    Successfully loaded "%s".\n', tag);
            catch ME
                 warning('    Error loading fdata.mat from %s: %s. Skipping.', tag, ME.message);
            end
        else
            warning('    No fdata.mat found in %s. Skipping.', tag);
        end
    end
end

% Check after attempting to load
if isempty(fdata_original_all_loaded)
    error('No valid fdata.mat files found for this sample. Exiting.');
end
% Update num_timepoints_expected to reflect what was actually loaded for processing
if numel(fdata_original_all_loaded) < num_timepoints_expected
    warning('Successfully loaded %d timepoints, but expected %d. Processing available data.', numel(fdata_original_all_loaded), num_timepoints_expected);
end
num_timepoints_to_process = numel(fdata_original_all_loaded);
if num_timepoints_to_process == 0, error('No timepoints available to process.'); end

%% PRE-PROCESS IMPORTED FDATA MATRICES AND SET NOISE FLOOR THRESHOLD (IF NOT SET IN PREPROCESSING_SETTINGS_APP)

fprintf('---------------------------------------------------------\n');
fprintf('CONFIGURING NOISE FLOOR THRESHOLD\n');
fprintf('---------------------------------------------------------\n');

if settings.use_fixed_threshold
    % settings.signal_strength_threshold is already set from settings.signalStrengthThresholdValue
    fprintf('Using fixed signal strength threshold from settings: %.2f\n', settings.signal_strength_threshold);
else
    fprintf('Calculating signal strength threshold automatically across loaded timepoints for this sample...\n');
    noise_values_calculated = NaN(1, num_timepoints_to_process); % YOUR ORIGINAL: noise = zeros(...)
    for i = 1:num_timepoints_to_process % Iterate over successfully loaded data
        fprintf('  Calculating for timepoint: %s...\n', timepointTags_loaded{i});
        current_fdata_original = fdata_original_all_loaded{i}; % Use loaded data
        
        % Call ProcessNeuroCa and smoothCalciumSignals for threshold calculation
        temp_signal_info = ProcessNeuroCa(current_fdata_original, settings.analysisDurationSeconds); % Pass duration
        if isempty(temp_signal_info) || ~isfield(temp_signal_info, 'fdata')
            warning('  Failed to process %s for threshold calculation. Skipping.', timepointTags_loaded{i});
            continue; % noise_values_calculated(i) remains NaN
        end
        temp_smooth_signal = smoothCalciumSignals(temp_signal_info, settings.n);
         if isempty(temp_smooth_signal) || size(temp_smooth_signal,2) == 0
            warning('  Failed to smooth %s (or got empty result) for threshold calculation. Skipping.', timepointTags_loaded{i});
            continue; % noise_values_calculated(i) remains NaN
        end
        
        try
            noise_values_calculated(i) = SignalStrengthThreshold(temp_smooth_signal, temp_signal_info.fdata);
        catch ME_thresh
            warning('  Error during SignalStrengthThreshold for %s: %s. Skipping.', timepointTags_loaded{i}, ME_thresh.message);
            % noise_values_calculated(i) remains NaN
        end
    end
    
    valid_noise_values = noise_values_calculated(~isnan(noise_values_calculated));
    if isempty(valid_noise_values)
        warning(['Could not calculate any valid noise thresholds for this sample. ', ...
                 'Using default fixed threshold value: %.2f'], defaultSettings.signalStrengthThresholdValue);
        settings.signal_strength_threshold = defaultSettings.signalStrengthThresholdValue;
    else
        settings.signal_strength_threshold = ceil(mean(valid_noise_values));
    end
    fprintf('Auto-calculated signal strength threshold for this sample: %.2f\n', settings.signal_strength_threshold);
end

%% PROCESS DATA FROM ALL TIMEPOINTS
parpool

fprintf('---------------------------------------------------------\n');
fprintf('PROCESSING TIMEPOINTS FOR THIS SAMPLE\n');
fprintf('---------------------------------------------------------\n');

% Initialize cell arrays to store results
all_calcium_spike_calculations = cell(1, num_timepoints_to_process);
all_decay_rates = cell(1, num_timepoints_to_process);
all_CWT_results = cell(1, num_timepoints_to_process);

% Process Baseline Timepoint (index 1 of loaded data)
fprintf('Processing Baseline Timepoint: %s...\n', timepointTags_loaded{1});
[calc_tp1, decay_tp1, cwt_tp1, calcium_spikes_tp1] = ... % Add calcium_spikes_tp1 as an output
        processTimepoint(rootFolder, timepointTags_loaded{1}, fdata_original_all_loaded{1}, settings, 0, 0, 0, 1);

all_calcium_spike_calculations{1} = calc_tp1;
all_calcium_spikes_structs{1} = calcium_spikes_tp1;
all_decay_rates{1} = decay_tp1;
all_CWT_results{1} = cwt_tp1;

% If more timepoints exist, process subsequent timepoints
for i = 2:num_timepoints_to_process
    fprintf('Processing Timepoint %d: %s...\n', i, timepointTags_loaded{i});
    % Retrieve baseline results for passing to processTimepoint
    baseline_calc_results = all_calcium_spike_calculations{1}; 
    baseline_decay_results = all_decay_rates{1};             
    baseline_cwt_results = all_CWT_results{1};  

    [calc_tpi, decay_tpi, cwt_tpi, calcium_spikes_tpi] = ...
        processTimepoint(rootFolder, timepointTags_loaded{i}, fdata_original_all_loaded{i}, settings, ...
                         baseline_calc_results.network_avg_spike_rate, ...
                         baseline_calc_results.network_avg_spike_intensity, ...
                         baseline_decay_results.network_avg_rate_cst, ... 
                         baseline_cwt_results.max_0minus); 
    all_calcium_spike_calculations{i} = calc_tpi;
    all_calcium_spikes_structs{i} = calcium_spikes_tpi;
    all_decay_rates{i} = decay_tpi;
    all_CWT_results{i} = cwt_tpi;
end

%% SAVE SETTINGS AND CALCULATIONS FOR EASY ACCESS BY POST-PROCESSING APP

fprintf('---------------------------------------------------------\n');
fprintf('SAVING RESULTS FOR POST-PROCESSING\n');
fprintf('---------------------------------------------------------\n');

% Save the main settings file used for this analysis run into the rootFolder
% This ensures the PostProcessingApp can access the same prominence, etc.
mainScriptPath = fileparts(mfilename('fullpath')); % Path to Main_Calcium_Code.m
originalSettingsFile = fullfile(mainScriptPath, 'preprocessing_settings.mat'); % The file next to the script
destinationSettingsFile = fullfile(rootFolder, 'preprocessing_settings.mat'); % Target: inside the sample's folder
if isfile(originalSettingsFile)
    copyfile(originalSettingsFile, destinationSettingsFile);
end


for i = 1:num_timepoints_to_process
    tp_tag = timepointTags_loaded{i};
    tp_output_folder = fullfile(rootFolder, tp_tag); % Results will be saved within each timepoint's folder

    if ~isfolder(tp_output_folder)
        mkdir(tp_output_folder);
    end

    fprintf('Saving data for timepoint: %s into %s\n', tp_tag, tp_output_folder);

    % 1. Calcium Spike Calculations
    calcium_spikes_output_for_saving = all_calcium_spikes_structs{i}; % This is the direct output of findCalciumSpikes
    save(fullfile(tp_output_folder, 'calcium_spikes_raw_output.mat'), 'calcium_spikes_output_for_saving'); % Save it with a distinct name

    % 2. Decay Rates (originally all_decay_rates{i})
    decay_rates_data = all_decay_rates{i};
    save(fullfile(tp_output_folder, 'decay_rates_data.mat'), 'decay_rates_data');

    % 3. CWT Results (if you want to load them later, originally all_CWT_results{i})
    % cwt_data = all_CWT_results{i};
    % save(fullfile(tp_output_folder, 'cwt_data.mat'), 'cwt_data');

    % 4. Signal Info (fdata, time, etc.) and Smoothed Signal
    % These were generated per timepoint. We need to ensure they are saved correctly.
    % The easiest is to re-generate them here or ensure 'processTimepoint' returns them
    % or that they are stored during the loop. Assuming fdata_original_all_loaded{i} and settings are available:

    current_fdata_original = fdata_original_all_loaded{i};
    signal_info_struct = ProcessNeuroCa(current_fdata_original, settings.analysisDurationSeconds);
    smooth_signal_matrix = smoothCalciumSignals(signal_info_struct, settings.n);

    save(fullfile(tp_output_folder, 'signal_info_struct.mat'), 'signal_info_struct');
    save(fullfile(tp_output_folder, 'smooth_signal_matrix.mat'), 'smooth_signal_matrix');

    % 5. If manualModelCheck was run for this timepoint, save its .mat file also
    %    Your current manualModelCheck saves files like 'ManualModelCheck_0minus.mat'
    %    It would be better if manualModelCheck saved into the tp_output_folder
    %    with a consistent name, e.g., 'manual_validation_results.mat'.
    %    For now, the app will assume these might be in the rootFolder or not exist yet.
end

fprintf('All necessary data for post-processing saved.\n');
%% PLOT HISTOGRAMS OF AVERAGE CELL SPIKE RATE AND AVERAGE CELL DECAY RATE CONSTANT

fprintf('---------------------------------------------------------\n');
fprintf('GENERATING PLOTS\n');
fprintf('---------------------------------------------------------\n');

num_for_hist_plots = numel(all_calcium_spike_calculations);

% Check if timepointTags_loaded has the expected number of elements
if numel(timepointTags_loaded) ~= num_timepoints_to_process
    warning(['Mismatch between number of loaded timepoint tags (%d) and number of processed results (%d). ', ...
             'Histogram labels might be affected. Using generic labels if needed.'], ...
             numel(timepointTags_loaded), num_timepoints_to_process);
    % Create generic tags as a fallback if there's a severe mismatch
    % (though this shouldn't happen if loading and processing are consistent)
    temp_tags = cell(1, num_timepoints_to_process);
    for i_tag = 1:num_timepoints_to_process
        temp_tags{i_tag} = sprintf('Timepoint %d', i_tag);
    end
    current_timepoint_tags = temp_tags;
else
    current_timepoint_tags = timepointTags_loaded; % Use the actual loaded tags
end

% --- Spike Rate Histogram ---
plot_inputs_spike_rate = cell(1, num_for_hist_plots);
labels_for_spike_hist = cell(1, num_timepoints_to_process); % To store corresponding labels

for k = 1:num_timepoints_to_process
    if ~isempty(all_calcium_spike_calculations{k}) && isfield(all_calcium_spike_calculations{k}, 'cell_spike_rate')
        current_data = all_calcium_spike_calculations{k}.cell_spike_rate;
        if ~isempty(current_data) && sum(~isnan(current_data(:))) > 1 % Check for actual data points
            plot_inputs_spike_rate{k} = current_data;
            labels_for_spike_hist{k} = current_timepoint_tags{k}; % Store the corresponding tag
        else
            plot_inputs_spike_rate{k} = []; % Mark as invalid/empty if not enough data
            labels_for_spike_hist{k} = []; % Also mark label as invalid
            warning('Spike rate data for processed timepoint index %d ("%s") is empty or has insufficient points for histogram.', k, current_timepoint_tags{k});
        end
    else
        plot_inputs_spike_rate{k} = []; % Mark as invalid/empty
        labels_for_spike_hist{k} = [];
        warning('Spike rate data missing or field missing for processed timepoint index %d ("%s") for histogram.', k, current_timepoint_tags{k});
    end
end

% Filter out empty datasets AND their corresponding labels
valid_indices_spike = ~cellfun('isempty', plot_inputs_spike_rate);
valid_plot_inputs_spike_rate = plot_inputs_spike_rate(valid_indices_spike);
valid_labels_for_spike_hist = labels_for_spike_hist(valid_indices_spike);

if ~isempty(valid_plot_inputs_spike_rate)
    [SpikeRateHist] = CalciumDataHistogram(valid_plot_inputs_spike_rate, ...
                                           valid_labels_for_spike_hist, ... % Pass the filtered labels
                                           'Spike Rate', 'spikes/sec', 0.006, 'Lognormal');
    if ~isempty(SpikeRateHist) && isgraphics(SpikeRateHist,'figure')
        savefig(SpikeRateHist, fullfile(rootFolder, 'SpikeRateHistogram.fig')); % Save in sample folder
        saveas(SpikeRateHist, fullfile(rootFolder, 'SpikeRateHistogram.png'));  % Save in sample folder
        fprintf('Spike Rate Histogram saved to sample folder.\n');
    else
        warning('Spike Rate Histogram was not generated or is not a valid figure handle.');
    end
else
    warning('No valid spike rate data to generate histogram.');
end

% --- Decay Rate Constant Histogram ---
% num_for_decay_hist_plots is num_timepoints_to_process
plot_inputs_decay_rate = cell(1, num_timepoints_to_process);
labels_for_decay_hist = cell(1, num_timepoints_to_process); % To store corresponding labels

for k = 1:num_timepoints_to_process
    if ~isempty(all_decay_rates{k}) && isfield(all_decay_rates{k}, 'cell_avg_rate_cst')
        current_data = all_decay_rates{k}.cell_avg_rate_cst;
        if ~isempty(current_data) && sum(~isnan(current_data(:))) > 1 % Check for actual data points
            plot_inputs_decay_rate{k} = current_data;
            labels_for_decay_hist{k} = current_timepoint_tags{k}; % Store the corresponding tag
        else
            plot_inputs_decay_rate{k} = []; % Mark as invalid/empty if not enough data
            labels_for_decay_hist{k} = [];
            warning('Decay rate data for processed timepoint index %d ("%s") is empty or has insufficient points for histogram.', k, current_timepoint_tags{k});
        end
    else
        plot_inputs_decay_rate{k} = []; % Mark as invalid/empty
        labels_for_decay_hist{k} = [];
        warning('Decay rate data missing or field missing for processed timepoint index %d ("%s") for histogram.', k, current_timepoint_tags{k});
    end
end

% Filter out empty datasets AND their corresponding labels
valid_indices_decay = ~cellfun('isempty', plot_inputs_decay_rate);
valid_plot_inputs_decay_rate = plot_inputs_decay_rate(valid_indices_decay);
valid_labels_for_decay_hist = labels_for_decay_hist(valid_indices_decay);

if ~isempty(valid_plot_inputs_decay_rate)
    [DecayRateConstantHist] = CalciumDataHistogram(valid_plot_inputs_decay_rate, ...
                                                   valid_labels_for_decay_hist, ... % Pass the filtered labels
                                                   'Decay Rate Constant', 's^{-1}', 0.02, 'Normal');
    if ~isempty(DecayRateConstantHist) && isgraphics(DecayRateConstantHist,'figure')
        savefig(DecayRateConstantHist, fullfile(rootFolder, 'DecayRateConstantHistogram.fig')); % Save in sample folder
        saveas(DecayRateConstantHist, fullfile(rootFolder, 'DecayRateConstantHistogram.png'));  % Save in sample folder
        fprintf('Decay Rate Constant Histogram saved to sample folder.\n');
    else
        warning('Decay Rate Constant Histogram was not generated or is not a valid figure handle.');
    end
else
    warning('No valid decay rate data to generate histogram.');
end

fprintf('---------------------------------------------------------\n');
fprintf('ANALYSIS COMPLETE FOR SAMPLE: %s\n', rootFolder);
fprintf('---------------------------------------------------------\n');

%% OPTIONAL: GENERATE FILTERED HISTOGRAMS WITH OUTLIER REMOVAL
% Uncomment the line below to run an interactive tool that allows you to
% re-generate the Spike Rate and/or Decay Rate histograms after removing
% statistical outliers. This can help improve the visual fit of the
% distribution curves (e.g., 'Normal' or 'Lognormal').
%
generateFilteredHistograms(all_calcium_spike_calculations, all_decay_rates, timepointTags_loaded, rootFolder);
%% Post-Processing App
% Uncomment line 443 when wanting to visualize the analysis just run or to
% calculate code accuracy metrics.
%
% PostProcessingApp

%% Relative Baseline Intensity Change App
% Uncomment line 449 when wanting to calculate the relative baseline
% fluorescent intensity change between two time points.
%
% RelativeBaselineChangeApp