B02_calAllParams.m

%% 
% HT-LMD - High-Throughput Learning, Memory, and Drug-Testing Platform for Small Aquatic Models
% Copyright (C) 2025 Gokul Rajan
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.


    %% Initialization
    clear; clc; close all;
    set(0,'DefaultFigureVisible','off');
    
    % Set your base directory here
    rootDir = uigetdir;
    if rootDir == 0
        error('No directory selected. Please select a valid directory.');
    end

    % Recursively find all subfolders
    subFolders = getSubfolders(rootDir);

    % Initialize the animal data structure and index counter
    animal = struct;
    lastAnimalIndex = 0;  % Initialize animal index counter

    %% Process each subfolder containing .mat files
    for k = 1:length(subFolders)
        currentFolder = subFolders{k};
        
        files = dir(fullfile(currentFolder, '*.mat'));
        
        % Process each file
        for j = 1:length(files)
            filename = files(j).name;

                % Extract keywords from the base filename
                keywords = regexp(filename, '[_-]', 'split');
            
                % Initialize variables to store group, phase, and day
                group = '';
                phase = '';
                day = '';
            
                % Dynamically detect group, phase, and day from the keywords
                for idx = 1:length(keywords)
                    keyword = keywords{idx};
            
                    % Check for group
                    if any(strcmp(keyword, {'c1', 'c2','c3', 't1', 't2','t3','mix1','mix2','mix3'}))
                    group = keyword;
                end
            
                    % Check for phase
                    if strcmp(keyword, 'pre') || strcmp(keyword, 'post')
                        phase = keyword;
                    end
            
                    % Check for day (d1, d2, d3, etc.)
                    if startsWith(keyword, 'd') && isstrprop(keyword(2), 'digit')
                        day = keyword;
                    end
                  end
    
        % Validate that group, phase, and day were found
        if isempty(group) || isempty(phase) || isempty(day)
            error(['Incomplete designation in filename %s. ', ...
                   'Expected format includes group (c1/t1), phase (pre/post), and day (d1, d2, d3, etc.).'], filename);
        end
                
                % Load the file
            data = load(fullfile(currentFolder, filename));
            camscale_px_per_mm = 8.4;
            data_rate_hz = 1;
            window_speed = 10 * data_rate_hz; % 1 minute window for speed smoothing
            bin_size_pos = 10 * data_rate_hz; % Binning size for position data (30 seconds)
            bin_size_other = 10 * data_rate_hz; % Binning size for other data (10 seconds)
            cs1_start_idx=data.cs1_start_idx;
        
            % Process data for each arena
            if isfield(data, 'arena')  % Check if 'arena' field is present
                for ii = 1:length(data.arena)
                    animalIndex = lastAnimalIndex + 1;  % Increment animal index
                    arenaData = processArenaData(data.arena(ii), camscale_px_per_mm, data_rate_hz, window_speed, bin_size_pos, bin_size_other,cs1_start_idx);
                    animal(animalIndex).data = arenaData;
                    animal(animalIndex).day = day;
                    animal(animalIndex).group = group;
                    animal(animalIndex).phase = phase;
                    animal(animalIndex).roi = ii;
                    lastAnimalIndex=animalIndex;
                end
            end
          end
      end
        
        % Save the animal structure
        save(fullfile(rootDir, 'processedAnimalData.mat'), 'animal');


%% Helper function to process data for each arena (stub function)
function arenaData = processArenaData(arena, camscale_px_per_mm, data_rate_hz, window_speed, bin_size_pos, bin_size_other,cs1_start_idx)
    
    bin_data = @(data, bin_size) arrayfun(@(k) mean(data(max(1, k-bin_size+1):min(end, k))), bin_size:bin_size:length(data));
    
    required_length = cs1_start_idx + 600;

    % Check and pad X if necessary
    if length(arena.X) < required_length
        padded_length = required_length - length(arena.X);
        % Check if X is a row or column vector and pad with NaNs accordingly
        if isrow(arena.X)
            arena.X = [arena.X, NaN(1, padded_length)];  % Pad with NaNs as a row
        else
            arena.X = [arena.X; NaN(padded_length, 1)];  % Pad with NaNs as a column
        end
    end

    % Check and pad Y if necessary
    if length(arena.Y) < required_length
        padded_length = required_length - length(arena.Y);
        % Check if Y is a row or column vector and pad with NaNs accordingly
        if isrow(arena.Y)
            arena.Y = [arena.Y, NaN(1, padded_length)];  % Pad with NaNs as a row
        else
            arena.Y = [arena.Y; NaN(padded_length, 1)];  % Pad with NaNs as a column
        end
    end
    
    % Convert and normalize coordinates
    X_mm = arena.X(cs1_start_idx:(cs1_start_idx+600))/ camscale_px_per_mm;
    Y_mm = arena.Y(cs1_start_idx:(cs1_start_idx+600))/ camscale_px_per_mm;

    X_normalized = (X_mm - min(X_mm)) / (max(X_mm) - min(X_mm));
    Y_normalized = (Y_mm - min(Y_mm)) / (max(Y_mm) - min(Y_mm));

    % Calculate and smooth speed
    deltaX_mm = diff(X_mm);
    deltaY_mm = diff(Y_mm);
    deltaDistance_mm = sqrt(deltaX_mm.^2 + deltaY_mm.^2);
    speed_mmps = [0; deltaDistance_mm] * data_rate_hz;  % Append zero to align lengths

    speed_mmps = fillmissing(speed_mmps, 'linear');
    speed_mmps_avg = movmean(speed_mmps, window_speed);

    % Midpoint calculations
    midpoint_mm = arena.center(1) / camscale_px_per_mm;
    midpoint_crossings = [0; diff(X_mm > midpoint_mm) ~= 0];
    midpoint_crossings_freq = movmean(midpoint_crossings, window_speed) * data_rate_hz * 60;

    % Performance index calculation
    performance_index = 2 * (X_mm > midpoint_mm) - 1;

    % Binning data using anonymous function defined outside
    X_normalized_binned = bin_data(X_normalized, bin_size_pos);
    Y_normalized_binned = bin_data(Y_normalized, bin_size_pos);
    speed_mmps_avg_binned = bin_data(speed_mmps_avg, bin_size_other);
    midpoint_crossings_freq_binned = bin_data(midpoint_crossings_freq, bin_size_other);
    performance_index_binned = bin_data(performance_index, bin_size_pos);

    % Pack data into a structure for output
    arenaData = struct;
    arenaData.X_mm=X_mm;
    arenaData.Y_mm=Y_mm;
    arenaData.X_normalized_binned = X_normalized_binned;
    arenaData.Y_normalized_binned = Y_normalized_binned;
    arenaData.X_normalized = X_normalized;
    arenaData.Y_normalized = Y_normalized;
    arenaData.speed_mmps_avg_binned = speed_mmps_avg_binned;
    arenaData.midpoint_crossings_freq_binned = midpoint_crossings_freq_binned;
    arenaData.performance_index_binned = performance_index_binned;
end

%% Recursive function to get all subfolders
function folders = getSubfolders(baseFolder)
    folders = {baseFolder}; % Include the base folder itself
    files = dir(baseFolder);
    for k = 1:length(files)
        if files(k).isdir && ~strcmp(files(k).name, '.') && ~strcmp(files(k).name, '..')
            folderPath = fullfile(baseFolder, files(k).name);
            folders = [folders, getSubfolders(folderPath)];  % Recursive call
        end
    end
end