% Function collects the receptive fields (X,Y) (On, Off) from all neurons 
% into a cell array 4 x nChannels with rows: (xOff, yOff, xOn, yOn)
% This reads only first layer recorded in a fly.
% For this only data with 5 stimuli can be used:
%       OnOff flashes plus 4 bar stimuli
function receptiveFieldsAll = collectRFs(dataEdges, tuningStat, varargin)

if nargin >= 3
    deleteLastROI = varargin{1};
else
    deleteLastROI = 1;
end

% Define the array to compare to randomized recordings.
xOffStimFileName = 'StandingStripe_1s_XAxis_5degWide_2degSep_m1.0Con_rand_USEFRUSTUM.txt';
yOffStimFileName = 'StandingStripe_1s_YAxis_5degWide_2degSep_m1.0Con_rand_USEFRUSTUM.txt';
xOnStimFileName = 'StandingStripe_1s_XAxis_5degWide_2degSep_p1.0Con_rand_USEFRUSTUM.txt';
yOnStimFileName = 'StandingStripe_1s_YAxis_5degWide_2degSep_p1.0Con_rand_USEFRUSTUM.txt';
barStimFileNames = {xOffStimFileName, yOffStimFileName, xOnStimFileName, yOnStimFileName};

nChannels = numel(dataEdges{1}.roiTCsInterp);
receptiveFieldsAll = cell(4, nChannels);
for iFly = 1: size(dataEdges, 1)
%     flyPathInd = iFly;
    for jStim = 2:5
        % Find index of matching stimulus.
        rfInd = find(strcmp(dataEdges{iFly, 1, jStim}.stimulusName, barStimFileNames));
        if ~isempty(rfInd)
            for kChannel = 1: nChannels
                nEpochs = numel(dataEdges{iFly, 1, jStim}.stimFrameInds); 
                nCells = size(dataEdges{iFly, 1, jStim}.roiTCs{kChannel}, 2);
                roiTCs = dataEdges{iFly, 1, jStim}.roiTCs{kChannel};
                receptiveFields = nan * ones(nCells, nEpochs);
                switch tuningStat
                    case 'maxAbs'                        
                        for iEpoch = 1: nEpochs
                            epochInds = dataEdges{iFly, 1, jStim}.stimFrameInds{iEpoch}(3: end - 1);
                            [~, maxEpochInd] = max(abs(squeeze(roiTCs(1, :, epochInds))), [], 2);
                            receptiveFields(:, iEpoch) = vec(arrayfun(@(x) squeeze(roiTCs(1, x, epochInds(maxEpochInd(x)))), 1: size(roiTCs, 2)));
                        end  
                    case 'median'
                        receptiveFields = cell2mat(cellfun(@(x) median(squeeze(roiTCs(1, :, x)), 2), ...
                                                                dataEdges{iFly, 1, jStim}.stimFrameInds, 'UniformOutput', false));
                    case 'mean'
                        receptiveFields = cell2mat(cellfun(@(x) mean(squeeze(roiTCs(1, :, x)), 2), ...
                                                                dataEdges{iFly, 1, jStim}.stimFrameInds, 'UniformOutput', false));
                end
                receptiveFieldsAll{rfInd, kChannel}  = [receptiveFieldsAll{rfInd, kChannel}; receptiveFields(1: end - 1 * deleteLastROI, :)];                
            end     
        end
    end
    clear receptiveFields
end

% Hard-coded error correction for some extra ROI.
if contains(dataEdges{1}.flyPath, 'Tm9GCaMP6f-Tm4jRGECO1a')
    receptiveFieldsAll{1,1}(10, :) =[];
end

end