%% Sort Data 12WP
% Script to organise data from axio statistics compiler. This script is
% specific to the 12 Well plate format, and is specific to Gemma's
% micrfluidc channel experiment.

%% Headers

% if ~exist('PLATE_FORMAT', 'var')
%     temppath = 'Header*.m';
%     if exist('my_path', 'var'); if my_path; temppath = [my_path temppath]; end; end
%     [header_file, my_path, indx] = uigetfile(temppath, 'Select Plate Format', 'MultiSelect', 'off');
%     if indx; run([my_path header_file]); end
%     clearvars temppath indx
% end

if ~exist('PLATE_LAYOUT', 'var')
    temppath = 'Plate*.m';
    if exist('my_path', 'var'); if my_path; temppath = [my_path temppath]; end; end
    [header_file, my_path, indx] = uigetfile(temppath, 'Select Plate Layout', 'MultiSelect', 'off');
    if indx; run([my_path header_file]); disp(['Plate layout: ' header_file]); end
    clearvars temppath indx
end

% Check Compatibility
% if ~strcmp(PLATE_FORMAT, EXPECTED_FORMAT)
%     errordlg('Filetype mismatch. Check plate layouts and try again.')
% end

% Generate some lists from the input formats.
% Groups = unique(Well_Group);
% no_groups = length(Groups);

% Col_Labels = 1:no_cols;
% Row_Labels = 1:no_rows;


%% Defaults
% Default inputs and constants.

% Define a minimum value of activity for a set to be included. To be included in the set averages. Only the Set_Start index is checked.
Cutoff_Activity = 100;

% Common Recording duration in seconds.
Recording_Duration = 300; % Time in seconds. ####################################################<<>>


%% Get The Spike Data
% Things that need to be declared by the user.

if ~exist('SpikeTables', 'var')
    [SpikeTables, InfoTables, my_path, filenames] = Spike_List_Import;
end
no_tables = size(SpikeTables, 2);
% no_tables = 1; % Limit no_tables for debugging.


%% Clean the Spike Data
% minimum activity for an electrode to be active etc...

SpikeTablesClean = SpikeTables;                                             % Table to store the cleaned electrodes.
disp('Cleaning tables...')
for i_table = 1:1:no_tables
    fprintf('Table %i: ', i_table)
    current_table = SpikeTablesClean{i_table};
    q_nonactive = current_table.No_Spikes < Cutoff_Activity;                % Find Electrodes that do not reach the miniumum activity.
    q_enabled = ismember(current_table.Well, Well_Name(Well_Enabled));      % Find Wells that are not enabled.
    current_table(q_nonactive | ~q_enabled, :) = [];                        % Remove selected Electrodes.
    fprintf('%i Electrodes removed.\n', sum(q_nonactive | ~q_enabled))
    SpikeTablesClean{i_table} = current_table;                              % Update the cleaned table list.
end
disp('  All Cleaned!')
clearvars i_table current_table q_nonactive q_enabled


%% Define Channels

disp('Defining Channels')
Table_Channel_Electrodes = table();
% Table_Channel_Electrodes
%  Ch_ID, Well, Ch#, Electrodes
Ch_ID = 0;
for i_w = 1:no_wells
    
    Well = Well_Name(i_w);
    Temp_Ch_Start = Ch_Start(i_w);
    Temp_Ch_End   = Ch_End(i_w);
    Temp_Ch_Range = Temp_Ch_Start:Temp_Ch_End;
    
    for i_ch = 1:no_channels             % Scan the channels in order.
        Ch_ID = Ch_ID + 1;
        
        % Collect data based on organoid placement in the well.
        switch upper(Well_Placement(i_w))   % Build the i_querry index.
            case "L"    % Organoid placed at the Left.
                Elec_List = Elec_Labels(i_ch, Temp_Ch_Range);
                
            case "R"    % Organoid placed at the Right.
                Elec_List = Elec_Labels(i_ch, 9 - Temp_Ch_Range);
                
            case "B"    % Organoid placed at the Bottom.
                Elec_List = Elec_Labels(9 - Temp_Ch_Range, i_ch)';
                
            case "T"    % Organoid placed at the Top.
                Elec_List = Elec_Labels(Temp_Ch_Range, i_ch)';
        end
        
        Channel = i_ch;
        Elec_List_Names = {append([char(Well) '_'], Elec_List)};
%         Elec_List = {Elec_List};
        Table_Channel_Electrodes(Ch_ID, :) = table(Well, Channel, Elec_List_Names);
        
    end
end
Table_Channel_Electrodes(~ismember(Table_Channel_Electrodes.Well, Well_Name(Well_Enabled)), :) = [];
clearvars i_ch i_w Ch_ID Well Channel Set_Label Elec_List Elec_List_Names Temp_Ch_Start Temp_Ch_End Temp_Ch_Range


%% Define Sets

% disp('Defining Sets')
% Table_Set_Electrodes = table();
% % Table_Set_Electrodes
% %  Ch_ID, Well, Ch#, Electrodes
% Set_ID = 0;
% for i_w = 1:no_wells
%     
%     Well = Well_Name(i_w);
%     Temp_Set_Start = Ch_Start(i_w);
%     Temp_Set_End   = Ch_End(i_w);
%     Temp_Ch_Range = 1:no_channels;
% 
%     for i_set = Temp_Set_Start:Temp_Set_End             % Scan the Sets in order.
%         Set_ID = Set_ID + 1;
%         
%         % Collect data based on organoid placement in the well.
%         switch upper(Well_Placement(i_w))   % Build the i_querry index.
%             case "L"    % Organoid placed at the Left.
%                 Elec_List = Elec_Labels(Temp_Ch_Range, i_set)';
%                 
%             case "R"    % Organoid placed at the Right.
%                 Elec_List = Elec_Labels(Temp_Ch_Range, no_sets + 1 - i_set)';
%                 
%             case "B"    % Organoid placed at the Bottom.
%                 Elec_List = Elec_Labels(no_sets + 1 - i_set, Temp_Ch_Range);
%                 
%             case "T"    % Organoid placed at the Top.
%                 Elec_List = Elec_Labels(i_set, Temp_Ch_Range);
%         end
%         
%         Set = i_set - Temp_Set_Start + 1;
%         Elec_List_Names = {append([char(Well) '_'], Elec_List)};
% %         Elec_List = {Elec_List};
%         Table_Set_Electrodes(Set_ID, :) = table(Well, Set, Elec_List_Names);
%         
%     end
% end
% Table_Set_Electrodes(~ismember(Table_Set_Electrodes.Well, Well_Name(Well_Enabled)), :) = [];
% clearvars i_set i_w Ch_ID Well Channel Set_Label Elec_List Elec_List_Names Temp_Set_Start Temp_Set_End Temp_Ch_Range


%% Find All Spikes from a single channel/set
% Store this in a table.

SpikeTables_Channels    = cell(1, no_tables);
% SpikeTables_Sets        = cell(1, no_tables);

for i_table = 1:1:no_tables
    fprintf('Grouping Table %i', i_table)
    SpikeTables_Channels(i_table)   = {GroupSpikeTimes(Table_Channel_Electrodes, SpikeTablesClean{i_table})};
%     SpikeTables_Sets(i_table)       = {GroupSpikeTimes(Table_Set_Electrodes,     SpikeTablesClean{i_table})};
    fprintf(' complete.\n')
end
clearvars i_table


%%  Filter the spikes for directionality 
% Method, A vs A+1, and B vs B-1. (No longer used since A and B definitions were sometimes not accurate enough.
% Method, find two highest firing electrodes in the channel and directionally filter each based on the other.

% Maximum sample propogation window between electrodes.
SPWindow = 3; % <<This should be moved to defaults or plate definitions?>>
SPExclude = 1; % <<This should be moved to defaults or plate definitions?>>

for i_table = 1:1:no_tables
    fprintf('Directional Filtering Table %i...', i_table)
    current_table = SpikeTables_Channels{i_table};
    
    total_chs = size(current_table, 1);
    Dir_AB = cell(total_chs, 1);
    Dir_BA = cell(total_chs, 1);

    eA = zeros(total_chs, 1);
    eB = zeros(total_chs, 1);
    eDiff = zeros(total_chs, 1);
    eActivity = cell(total_chs, 1);
    for i_ch = 1:total_chs
        SpikeTimes = current_table.G_Times{i_ch};
        eActivity{i_ch} = cellfun('size', SpikeTimes, 2);

        [~, eActive] = maxk(eActivity{i_ch}, 2);
        temp_eA = min(eActive);
        temp_eB = max(eActive);
        xDiff = temp_eB - temp_eA;

        if ~isempty(xDiff)
            eA(i_ch, 1) = temp_eA;
            eB(i_ch, 1) = temp_eB;
            eDiff(i_ch, 1) = xDiff;
            SampleWindow = eDiff(i_ch) * SPWindow;
            ExcludeWindow = max([1 (eDiff(i_ch) * SPExclude)]); % max([0 ((eDiff(i_ch) * SPExclude) - 1)]); %Changed this as the filtering function uses the exclusion window inclusive.

            ST_A = SpikeTimes{eA(i_ch)};
            ST_B = SpikeTimes{eB(i_ch)};
    
            [ST_A, ST_B] = BiDirFilterSpikeTimes(ST_A, ST_B, SampleWindow, ExcludeWindow);
            
            Dir_AB(i_ch) = {ST_A};
            Dir_BA(i_ch) = {ST_B};
        end
    end
    
    current_table = [current_table table(eA, eB, eActivity, Dir_AB, Dir_BA)];
    
    SpikeTables_Channels{i_table} = current_table;
    fprintf('\b\b\b completed.\n')
end
clearvars i_table current_table total_chs Dir_AB Dir_BA eA eB eDiff xDiff i_ch SpikeTimes eActivity eActive temp_eA temp_eB SampleWindow ST_A ST_B


%% Convert channels to sides
% Restructure the spike table into the filtered table format.

FilteredTables = cell(1, no_tables);
for i_table = 1:1:no_tables
    DirFilteredTable = table();
    current_table = SpikeTables_Channels{i_table};
    total_chs = size(current_table, 1);
    for i_ch = 1:total_chs
        outTableLine = table();
        inTableLine = current_table(i_ch, :);
        Side = {'AB'};
        Measurement = {char(join([inTableLine.Well, num2str(inTableLine.Channel), Side{:}], '_'))};
        FilteredSpikes = inTableLine.Dir_AB;
        outTableLine(1, :) = [table(Measurement, Side) inTableLine(:, {'Well', 'Channel'}) table(FilteredSpikes)];

        Side = {'BA'};
        Measurement = {char(join([inTableLine.Well, num2str(inTableLine.Channel), Side{:}], '_'))};
        FilteredSpikes = inTableLine.Dir_BA;
        outTableLine(2, :) = [table(Measurement, Side) inTableLine(:, {'Well', 'Channel'}) table(FilteredSpikes)];


        if isempty(DirFilteredTable)
            DirFilteredTable = outTableLine;
        else
            DirFilteredTable = [DirFilteredTable; outTableLine];
        end
    end
    FilteredTables(i_table) = {DirFilteredTable};
end
clearvars i_table current_table total_chs i_ch outTableLine inTableLine Side Measurement FilteredSpikes


%% Produce Stats Summary
% Go through each filtered data set and find the standard stats available.

disp('Producing Summary Information:')
StatsTables = cell(1, no_tables);
BurstTables = cell(1, no_tables);
for i_table = 1:1:no_tables
    fprintf('    Table %i\n', i_table);
    myStatsTable = table();                 % Generate/clear a table to hold the output stats.
    myBurstTable = table();
    input_table = FilteredTables{i_table};  % Select the input table.
    total_chs = size(input_table, 1);
    for i_ch = 1:total_chs                  % Go through all associated channels in the input table to produce the associated stats.
        SpikeTimes = input_table.FilteredSpikes{i_ch};
        
        [outputTable, outputBursts] = Times2Stats(SpikeTimes, Recording_Duration);  % Find the standard stats. % <<Possibly could add the cross corellation stats here.
        if isempty(myStatsTable)                                    % Store the stats if returned successfuly.
            myStatsTable = outputTable;
            myBurstTable = [outputTable {{outputBursts}}];
        else
            myBurstTable = [myBurstTable; [outputTable {{outputBursts}}]];
            myStatsTable = [myStatsTable; outputTable];
        end
    end
    myStatsTable = [input_table myStatsTable];  
    myBurstTable = [input_table myBurstTable];                      % Add the output stats to the input table. Thus labelling is preserved.
    myStatsTable = renamevars(myStatsTable, 'Well', 'WellName');    % To avoid conflict later, change Well to Wellname.
    StatsTables(i_table) = {myStatsTable};                          % Store the overall stats table.
    BurstTables(i_table) = {myBurstTable};  
    fprintf('\b completed.\n');
end
clearvars i_table myStatsTable input_table total_chs i_ch SpikeTimes outputTable


%% Create output stats sheet
% Generate a cell array to represent a spreadsheet to store the output data.
StatsSheets = cell(1, no_tables);
for i_sheet = 1:1:no_tables
    temp_sheet = [StatsTables{i_sheet}.Properties.VariableNames' table2cell(StatsTables{i_sheet})'];
    temp_sheet(ismember(temp_sheet(:, 1), {'FilteredSpikes', 'Spike_Times', 'Spike_Mags'}), :) = []; % Remove the spike time lists, they would lead to a very awkward sheet.
    
    temp_sheet = [InfoTables{i_sheet} cell(size(InfoTables{i_sheet}, 1), (size(temp_sheet, 2) - size(InfoTables{i_sheet}, 2))); temp_sheet];
    ts_str = cellfun(@(temp_sheet) any(isstring(temp_sheet)), temp_sheet);              % Find all string Values.
    temp_sheet(ts_str) = cellfun(@char, temp_sheet(ts_str), 'UniformOutput', false);    % Convert them to char Values.
    temp_sheet(cellfun(@(temp_sheet) any(isnan(temp_sheet)), temp_sheet)) = {[]};       % Replace all NaN values with blanks.

    StatsSheets(i_sheet) = {temp_sheet};
end
clearvars i_sheet temp_sheet ts_str


%% Save output stats sheets
% Save the output stats and all of scripts used to generate them.
Timestamp = char(datetime("now", "Format", '_yyyyMMdd-HHmm'));

Save_Sheets_Folder = [my_path 'DirSheets' Timestamp filesep];
if ~exist(Save_Sheets_Folder, 'dir');   mkdir(Save_Sheets_Folder);  end

disp('Saving Converted Sheets...')
for i_sheet = 1:1:no_tables
    my_name = [Save_Sheets_Folder 'DirF_' filenames{i_sheet}];
    writecell(StatsSheets{i_sheet}, my_name);
    disp(['    ' my_name])
end
disp('  Sheets Saved!')
disp(['<a href="matlab: winopen(''' Save_Sheets_Folder ''')">Link to save location...</a>']);
clearvars i_sheet my_name

% Find and save dependancies ###################################################################################################
% This will save a copy of this script, and other dependencies, in a folder that has the same name as the main sheet just saved.

Save_Script_Folder = [my_path 'Scripts' Timestamp filesep];
if ~exist(Save_Script_Folder, 'dir');   mkdir(Save_Script_Folder);  end

disp('Saving Scripts...')
% Locate this script.
ScriptPath = matlab.desktop.editor.getActiveFilename;

% Locate the dependancies.
req_files = matlab.codetools.requiredFilesAndProducts(ScriptPath);

% Save copies of the dependancies.
for i = 1:1:length(req_files)
    [~, name, exten] = fileparts(req_files{i});
    copyfile(req_files{i}, [Save_Script_Folder name exten]);
end
clearvars i name exten

disp('  Scripts Saved!')
clearvars ScriptPath req_files Timestamp
